text:
const number_of_errors_happen_so_far = Threads.Atomic{Int}(0)
@inline function parse_data!(res, buffer, types, cc, en, current_line, char_buff, char_cnt, df, dt_cnt, int_cnt, j, informat, int_bases, string_trim)
        flag = 0
        if !isempty(informat)
            _infmt! = get(informat, j, identity)
            if _infmt! !== identity
                _infmt!(buffer, cc, en)
            end
        end

        @inbounds if types[j] <: Int64
            flag = buff_parser(res[j]::Vector{Union{Missing, Int64}}, buffer, cc, en, current_line, Int64; base = int_bases === nothing ? 10 : int_bases[int_cnt])
        elseif types[j] <: Float64
            flag = buff_parser(res[j]::Vector{Union{Missing, Float64}}, buffer.data, cc, en, current_line, Float64)
        elseif types[j] <: Bool
            flag = buff_parser(res[j]::Vector{Union{Missing, Bool}}, buffer, cc, en, current_line, Bool)
        elseif types[j] <: Date
            flag = buff_parser(res[j]::Vector{Union{Missing, Date}}, buffer, cc, en, current_line, df[dt_cnt], Date)
        elseif types[j] <: DateTime
            flag = buff_parser(res[j]::Vector{Union{Missing, DateTime}}, buffer, cc, en, current_line, df[dt_cnt], DateTime)
        elseif types[j] <: String
            if string_trim
                flag = buff_parser(res[j]::Vector{Union{Missing, String}}, buffer.data, cc, en, current_line, string_trim, String)
            else
                flag = buff_parser(res[j]::Vector{Union{Missing, String}}, buffer.data, cc, en, current_line, String)
            end
        elseif types[j] <: Int32
            flag = buff_parser(res[j]::Vector{Union{Missing, Int32}}, buffer, cc, en, current_line, Int32; base = int_bases === nothing ? 10 : int_bases[int_cnt])
        elseif types[j] <: Float32
            flag = buff_parser(res[j]::Vector{Union{Missing, Float32}}, buffer.data, cc, en, current_line, Float32)
        elseif types[j] <: Int8
            flag = buff_parser(res[j]::Vector{Union{Missing, Int8}}, buffer, cc, en, current_line, Int8; base = int_bases === nothing ? 10 : int_bases[int_cnt])
        elseif types[j] <: Int16
            flag = buff_parser(res[j]::Vector{Union{Missing, Int16}}, buffer, cc, en, current_line, Int16; base = int_bases === nothing ? 10 : int_bases[int_cnt])
        elseif types[j] <: String1
            flag = buff_parser(res[j]::Vector{Union{Missing,String1}}, buffer.data, cc, en, current_line,String1)
        elseif types[j] <: String3
            flag = buff_parser(res[j]::Vector{Union{Missing,String3}}, buffer.data, cc, en, current_line,String3)
        elseif types[j] <: String7
            flag = buff_parser(res[j]::Vector{Union{Missing,String7}}, buffer.data, cc, en, current_line,String7)
        elseif types[j] <: String15
            flag = buff_parser(res[j]::Vector{Union{Missing,String15}}, buffer.data, cc, en, current_line,String15)
        elseif types[j] <: Characters{3, UInt8}
            flag = buff_parser(res[j]::Vector{Union{Missing,Characters{3, UInt8}}}, buffer.data, cc, en, current_line, char_buff[char_cnt]::Vector{UInt8}, Characters{3, UInt8})
        elseif types[j] <: Characters{5, UInt8}
            flag = buff_parser(res[j]::Vector{Union{Missing,Characters{5, UInt8}}}, buffer.data, cc, en, current_line, char_buff[char_cnt]::Vector{UInt8}, Characters{5, UInt8})
        elseif types[j] <: Characters{8, UInt8}
            flag = buff_parser(res[j]::Vector{Union{Missing,Characters{8, UInt8}}}, buffer.data, cc, en, current_line, char_buff[char_cnt]::Vector{UInt8}, Characters{8, UInt8})
        elseif types[j] <: Characters{10, UInt8}
            flag = buff_parser(res[j]::Vector{Union{Missing,Characters{10, UInt8}}}, buffer.data, cc, en, current_line, char_buff[char_cnt]::Vector{UInt8}, Characters{10, UInt8})
        elseif types[j] <: Characters{11, UInt8}
            flag = buff_parser(res[j]::Vector{Union{Missing,Characters{11, UInt8}}}, buffer.data, cc, en, current_line, char_buff[char_cnt]::Vector{UInt8}, Characters{11, UInt8})
        elseif types[j] <: Characters{12, UInt8}
            flag = buff_parser(res[j]::Vector{Union{Missing,Characters{12, UInt8}}}, buffer.data, cc, en, current_line, char_buff[char_cnt]::Vector{UInt8}, Characters{12, UInt8})
        elseif types[j] <: Characters{13, UInt8}
            flag = buff_parser(res[j]::Vector{Union{Missing,Characters{13, UInt8}}}, buffer.data, cc, en, current_line, char_buff[char_cnt]::Vector{UInt8}, Characters{13, UInt8})
        elseif types[j] <: Characters{14, UInt8}
            flag = buff_parser(res[j]::Vector{Union{Missing,Characters{14, UInt8}}}, buffer.data, cc, en, current_line, char_buff[char_cnt]::Vector{UInt8}, Characters{14, UInt8})
        elseif types[j] <: Characters{15, UInt8}
            flag = buff_parser(res[j]::Vector{Union{Missing,Characters{15, UInt8}}}, buffer.data, cc, en, current_line, char_buff[char_cnt]::Vector{UInt8}, Characters{15, UInt8})
        elseif types[j] <: Characters{1, UInt8}
            flag = buff_parser(res[j]::Vector{Union{Missing,Characters{1, UInt8}}}, buffer.data, cc, en, current_line, char_buff[char_cnt]::Vector{UInt8}, Characters{1, UInt8})
        elseif types[j] <: Characters{2, UInt8}
            flag = buff_parser(res[j]::Vector{Union{Missing,Characters{2, UInt8}}}, buffer.data, cc, en, current_line, char_buff[char_cnt]::Vector{UInt8}, Characters{2, UInt8})
        elseif types[j] <: Characters{4, UInt8}
            flag = buff_parser(res[j]::Vector{Union{Missing,Characters{4, UInt8}}}, buffer.data, cc, en, current_line, char_buff[char_cnt]::Vector{UInt8}, Characters{4, UInt8})
        elseif types[j] <: Characters{6, UInt8}
            flag = buff_parser(res[j]::Vector{Union{Missing,Characters{6, UInt8}}}, buffer.data, cc, en, current_line, char_buff[char_cnt]::Vector{UInt8}, Characters{6, UInt8})
        elseif types[j] <: Characters{7, UInt8}
            flag = buff_parser(res[j]::Vector{Union{Missing,Characters{7, UInt8}}}, buffer.data, cc, en, current_line, char_buff[char_cnt]::Vector{UInt8}, Characters{7, UInt8})
        elseif types[j] <: Characters{9, UInt8}
            flag = buff_parser(res[j]::Vector{Union{Missing,Characters{9, UInt8}}}, buffer.data, cc, en, current_line, char_buff[char_cnt]::Vector{UInt8}, Characters{9, UInt8})
        elseif types[j] <: Characters{16, UInt8}
            flag = buff_parser(res[j]::Vector{Union{Missing,Characters{16, UInt8}}}, buffer.data, cc, en, current_line, char_buff[char_cnt]::Vector{UInt8}, Characters{16, UInt8})
        elseif types[j] <: Time
            flag = buff_parser(res[j]::Vector{Union{Missing, Time}}, buffer, cc, en, current_line, df[dt_cnt], Time)
        elseif types[j] <: String31
            flag = buff_parser(res[j]::Vector{Union{Missing,String31}}, buffer.data, cc, en, current_line,String31)
        elseif types[j] <: String63
            flag = buff_parser(res[j]::Vector{Union{Missing,String63}}, buffer.data, cc, en, current_line,String63)
        elseif types[j] <: String127
            flag = buff_parser(res[j]::Vector{Union{Missing,String127}}, buffer.data, cc, en, current_line,String127)
        elseif types[j] <: String255
            flag = buff_parser(res[j]::Vector{Union{Missing,String255}}, buffer.data, cc, en, current_line,String255)
        elseif types[j] <: String255
            flag = buff_parser(res[j]::Vector{Union{Missing,String255}}, buffer.data, cc, en, current_line,String255)
        elseif types[j] <: UInt8
            flag = buff_parser(res[j]::Vector{Union{Missing, UInt8}}, buffer, cc, en, current_line, UInt8; base = int_bases === nothing ? 10 : int_bases[int_cnt])
        elseif types[j] <: UInt16
            flag = buff_parser(res[j]::Vector{Union{Missing, UInt16}}, buffer, cc, en, current_line, UInt16; base = int_bases === nothing ? 10 : int_bases[int_cnt])
        elseif types[j] <: UInt32
            flag = buff_parser(res[j]::Vector{Union{Missing, UInt32}}, buffer, cc, en, current_line, UInt32; base = int_bases === nothing ? 10 : int_bases[int_cnt])
        elseif types[j] <: UInt64
            flag = buff_parser(res[j]::Vector{Union{Missing, UInt64}}, buffer, cc, en, current_line, UInt64; base = int_bases === nothing ? 10 : int_bases[int_cnt])
        elseif types[j] <: Int128
            flag = buff_parser(res[j]::Vector{Union{Missing, Int128}}, buffer, cc, en, current_line, Int128; base = int_bases === nothing ? 10 : int_bases[int_cnt])
        elseif types[j] <: BigFloat
            flag = buff_parser(res[j]::Vector{Union{Missing, BigFloat}}, buffer, cc, en, current_line, BigFloat)
        else # anything else
            flag = buff_parser(res[j]::Vector{Union{Missing, String}}, buffer.data, cc, en, current_line, String)
        end
        flag
end

@inline function _process_iobuff!(res, buffer, types, dlm, eol, cnt_read_bytes, buffsize, current_line, last_line, last_valid_buff, charbuff, df, fixed, dlmstr, informat, quotechar, escapechar, warn, colnames, int_bases, string_trim, ignorerepeated, limit)
    n_cols = length(types)
    line_start = 1
    current_cursor_position = 1

    #keeping track of number of warnings
    warn_pass_end_of_line = 0

    dlmstr === nothing ? dlm_length = 1  : dlm_length = length(dlm)
    anything_is_wrong = 0
    any_problem_with_parsing = 0
    while true

        # keep track of Characters and DateTime columns
        char_cnt = 0
        dt_cnt = 0
        int_cnt = 0

        line_end = find_end_of_line(buffer.data, line_start, last_valid_buff, eol)
        field_start = line_start
        any_problem_with_parsing = 0
        for j in 1:n_cols
            if types[j] <: Characters
                char_cnt += 1
            elseif types[j] <: TimeType
                dt_cnt += 1
            elseif types[j] <: Integer
                int_cnt += 1
            end
            # if there is no fixed width information for the current column
            if fixed === 0:0 || fixed[j].start === 0
                if quotechar !== nothing
                    dlm_pos, new_lo, new_hi = find_next_delim(buffer.data, field_start, line_end, dlm, dlmstr, quotechar, escapechar, ignorerepeated)
                    if new_lo != 0 && new_hi != 0 && new_lo <= new_hi
                        new_lo, new_hi = clean_escapechar!(buffer, new_lo, new_hi, quotechar, escapechar)
                    end
                else
                    dlm_pos, new_lo, new_hi = find_next_delim(buffer.data, field_start, line_end, dlm, dlmstr, ignorerepeated)
                end
                # we should have a strategy for this kind of problem, for now just let the end of line as endpoint
                dlm_pos == 0 ? dlm_pos = line_end + dlm_length : nothing
                anything_is_wrong = parse_data!(res, buffer, types, new_lo == 0 ? field_start : new_lo, new_hi == 0 ? dlm_pos - dlm_length : new_hi, current_line, charbuff, char_cnt, df, dt_cnt, int_cnt, j, informat, int_bases, string_trim)
                field_start = dlm_pos + 1
            else # we have a fixed width information for the current column
                offset = line_start - 1
                dlm_pos = fixed[j].stop + offset
                if dlm_pos > line_end
                    dlm_pos = line_end
                    warn_pass_end_of_line += 1
                    # TODO use better wording
                    if warn_pass_end_of_line < 5
                        @warn "for column $(j) in line $(current_line[]) the cursor goes beyond the line, and the value is truncated"
                    end
                end
                # dlm_pos doesn't contain dlm so it shouldn't be dlm_pos-1 like the non-fixed case
                anything_is_wrong = parse_data!(res, buffer, types, fixed[j].start + offset, dlm_pos, current_line, charbuff, char_cnt, df, dt_cnt, int_cnt, j, informat, int_bases, string_trim)
                field_start = dlm_pos + 1
            end
            any_problem_with_parsing += anything_is_wrong
            dlm_pos > line_end && break
        end
        if any_problem_with_parsing>0
            Threads.atomic_add!(number_of_errors_happen_so_far, 1)
            if number_of_errors_happen_so_far[] <= warn
                if colnames !== nothing
                    @warn "There are problems with parsing data in line $(current_line[]), the values are set as missing:$(_write_warn_detail(buffer.data, line_start, line_end, res, current_line[], colnames))"
                end
            end
        end

        line_start = line_end + length(eol) + 1
        current_line[] += 1
        current_line[] > limit && return
        line_start > last_valid_buff && break
    end

end

# different function for reading multiple observations per line - the main different is here we push one observation at a time
@inline function _process_iobuff_multiobs!(res, buffer, types, dlm, eol, cnt_read_bytes, buffsize, current_line, last_line, last_valid_buff, charbuff, df, fixed, dlmstr, informat, quotechar, escapechar, warn, colnames, int_bases, string_trim, ignorerepeated, limit)
    n_cols = length(types)
    line_start = 1
    current_cursor_position = 1

    #keeping track of number of warnings
    warn_pass_end_of_line = 0

    dlmstr === nothing ? dlm_length = 1  : dlm_length = length(dlm)
    anything_is_wrong = 0
    any_problem_with_parsing = 0
    j = 1
    field_start = 1
    read_one_obs = false
    while true

        # keep track of Characters and DateTime columns
        char_cnt = 0
        dt_cnt = 0
        int_cnt = 0

        line_end = find_end_of_line(buffer.data, line_start, last_valid_buff, eol)
        any_problem_with_parsing = 0
        if types[j] <: Characters
            char_cnt += 1
        elseif types[j] <: TimeType
            dt_cnt += 1
        elseif types[j] <: Integer
            int_cnt += 1
        end

        if quotechar !== nothing
            dlm_pos, new_lo, new_hi = find_next_delim(buffer.data, field_start, line_end, dlm, dlmstr, quotechar, escapechar, ignorerepeated)
            if new_lo != 0 && new_hi != 0 && new_lo <= new_hi
                new_lo, new_hi = clean_escapechar!(buffer, new_lo, new_hi, quotechar, escapechar)
            end
        else
            dlm_pos, new_lo, new_hi = find_next_delim(buffer.data, field_start, line_end, dlm, dlmstr, ignorerepeated)
        end
        # we should have a strategy for this kind of problem, for now just let the end of line as endpoint
        dlm_pos == 0 ? dlm_pos = line_end + dlm_length : nothing
        anything_is_wrong = parse_data!(res, buffer, types, new_lo == 0 ? field_start : new_lo, new_hi == 0 ? dlm_pos - dlm_length : new_hi, current_line, charbuff, char_cnt, df, dt_cnt, int_cnt, j, informat, int_bases, string_trim)
        field_start = dlm_pos + 1
        read_one_obs = true
        any_problem_with_parsing += anything_is_wrong
        j = j + 1
        if j > n_cols
            map(x->push!(x, missing), res)
            read_one_obs = false
            current_line[] += 1
            current_line[]> limit && return read_one_obs
            j = 1
        end

        if j == n_cols && any_problem_with_parsing>0
            Threads.atomic_add!(number_of_errors_happen_so_far, 1)
            if number_of_errors_happen_so_far[] <= warn
                if colnames !== nothing
                    @warn "There are problems with parsing data in line $(current_line[]), the values are set as missing:$(_write_warn_detail(buffer.data, line_start, line_end, res, current_line[], colnames))"
                end
            end
        end
        if dlm_pos > line_end
            line_start = line_end + length(eol) + 1
            field_start = line_start
        end
        line_start > last_valid_buff && break
    end
    read_one_obs
end

# lo is the begining of the read and hi is the end of read. hi should be end of file or a linebreak
function readfile_chunk!(res, llo, lhi, charbuff, path, types, n, lo, hi, colnames; delimiter = ',', linebreak = '\n', lsize = 2^15, buffsize = 2^16, fixed = 0:0, df = dateformat"yyyy-mm-dd", dlmstr = nothing, informat = Dict{Int, Function}(), escapechar = nothing, quotechar = nothing, warn = 20, eolwarn = true, int_bases = nothing, string_trim = false, ignorerepeated = false, multiple_obs = false, limit = typemax(Int))
    read_one_obs = true
    f = OUR_OPEN(path, read = true)
    try
        if dlmstr === nothing
            dlm = UInt8.(delimiter)
        elseif dlmstr isa AbstractString
            dlm = collect(codeunits(dlmstr))
        else
            throw(ArgumentError("`dlmstr` must be a string of ASCII characters"))
        end
        eol = UInt8.(linebreak)
        eol_first = first(eol)
        eol_last = last(eol)
        eol_len = length(eol)

        buffsize > hi-lo+1 ? buffsize = hi-lo+1 : nothing

        # later we will use LineBuffer for this
        buffer = LineBuffer(Vector{UInt8}(undef, buffsize))
        current_line = Ref{Int}(llo)
        last_line = false
        last_valid_buff = buffsize
        # position which reading should be started
        seek(f, max(0, lo - 1))

        # TODO we should have a strategy when line is to long
        while true
            cnt_read_bytes = readbytes!(f, buffer.data)
            if cnt_read_bytes == 0
                return res
            end
            cur_position = position(f)

            if !eof(f) && cur_position < hi
                if buffer.data[end] !== eol_last || buffer.data[end - eol_len + 1] !== eol_first
                    #this means the buffer is not ended with a eol char, so we move back into buffer to have complete line
                    back_cnt = 0
                    for i in buffsize:-1:1
                        last_valid_buff = i
                        buffer.data[i] == eol_last && buffer.data[i - eol_len + 1] == eol_first && break
                        back_cnt += 1
                    end
                    cur_position = position(f)
                    seek(f, cur_position - back_cnt)
                else
                    last_valid_buff = buffsize
                end
            elseif cur_position == hi
                last_line = true
                if buffer.data[cnt_read_bytes] != eol_last || buffer.data[cnt_read_bytes - eol_len + 1] != eol_first
                    if eolwarn
                        @warn "the last line is not ended with `end of line` character"
                    end
                    # we read it anyway
                    # if !eof(f)
                    #     for i in cnt_read_bytes:-1:1
                    #         last_valid_buff = i
                    #         buffer.data[i] == eol && break
                    #     end
                    # end
                    last_valid_buff = cnt_read_bytes
                else
                    last_valid_buff = cnt_read_bytes
                end
            else cur_position > hi
                last_line = true
                last_valid_buff = buffsize - (cur_position - hi + 1)
            end
            if multiple_obs
                read_one_obs = _process_iobuff_multiobs!(res, buffer, types, dlm, eol, cnt_read_bytes, buffsize, current_line, last_line, last_valid_buff, charbuff, df, fixed, dlmstr, informat, quotechar, escapechar, warn, colnames, int_bases, string_trim, ignorerepeated, limit)
            else
                _process_iobuff!(res, buffer, types, dlm, eol, cnt_read_bytes, buffsize, current_line, last_line, last_valid_buff, charbuff, df, fixed, dlmstr, informat, quotechar, escapechar, warn, colnames, int_bases, string_trim, ignorerepeated, limit)
            end
                # we need to break at some point
            current_line[] > limit && break
            last_line && break
        end
    catch e
        CLOSE(f)
        rethrow(e)
    end
    CLOSE(f)
    if !read_one_obs
        map(x->pop!(x), res)
    end
    res
end

# main distributer
function distribute_file(path, types; delimiter = ',', linebreak = '\n', header = true, threads = true, guessingrows = 20, fixed = 0:0, buffsize = 2^16, quotation = nothing, dtformat = dateformat"yyyy-mm-dd", lsize = 2^15, dlmstr = nothing, informat = Dict{Int, Function}(), escapechar = nothing, quotechar = nothing, warn = 20, eolwarn = true, emptycolname = false, int_bases = nothing, string_trim = false, makeunique = false, ignorerepeated = false, multiple_obs = false, skipto = 1, limit = typemax(Int))
    eol = UInt8.(linebreak)
    eol_first = first(eol)
    eol_last = last(eol)
    eol_len = length(eol)
    colwidth = 0:0

    f = OUR_OPEN(path, read = true)
    # generating varnames
    f_pos = 0
    if skipto != 1
        _fsize_ = FILESIZE(f)
        l_length, f_pos = read_multiple_lines(path, f_pos, _fsize_, eol, skipto-1)
    else
        f_pos = 0
    end

    if !multiple_obs

        lsize_estimate = estimate_linesize(path, eol, lsize, f_pos, guessingrows = guessingrows)
        lsize_estimate > lsize && throw(ArgumentError("the lines are larger than $lsize, increase buffers by setting `lsize` and `buffsize` arguments, they are currently set as $lsize and $buffsize respectively"))
        if fixed != 0:0
            colwidth = Vector{UnitRange{Int}}(undef, length(types))
            for i in 1:length(colwidth)
                colwidth[i] = get(fixed, i, 0:0)
            end
        else
            colwidth = 0:0
        end
    end

    if (header isa AbstractVector) && (eltype(header) <: Union{AbstractString, Symbol})
        # colnames = header
        colnames = InMemoryDatasets.make_unique(header; makeunique = makeunique)
    elseif header === true
        f_pos, colnames = _generate_colname_based(path, eol, f_pos+1, lsize, lsize, types, delimiter, linebreak, buffsize, colwidth, dlmstr, quotechar, escapechar, emptycolname, ignorerepeated, multiple_obs)
        isempty(colnames) && throw(ArgumentError("Detecting column names return empty array, this can happen if your file is empty or the column names are not detectable."))
        colnames = InMemoryDatasets.make_unique(colnames; makeunique = makeunique)
    elseif header === false
        colnames = ["x"*string(k) for k in 1:length(types)]
    else
        throw(ArgumentError("`header` can be true or false, or a list of variable names"))
    end
    if !multiple_obs
        # how many bytes we should skip - we should use this information to have a better distribution of file into nt chunks.
        skip_bytes = f_pos
        fs = FILESIZE(f)
        initial_guess_for_row_num = (fs-skip_bytes)/lsize_estimate
        _check_nt_possible_size = div(fs-skip_bytes, buffsize)
        if !threads || fs-skip_bytes < 10^6 || div(fs-skip_bytes, Threads.nthreads()) < lsize_estimate || path isa IOBuffer || big(limit)*lsize_estimate < 10^7
            nt = 1
        else
            if _check_nt_possible_size > 0
                nt = min(Threads.nthreads(), _check_nt_possible_size)
            else
                nt = 1
            end
        end
    else
        nt = 1
    end

    charbuff = Vector{Vector{UInt8}}()

    # this is temp approach to TimeType parsing, later we should use DT{N} structure to parse TimeType
    if dtformat isa Dict
        dtfmt = Vector{valtype(dtformat)}()
    else
        dtfmt = Vector{typeof(_todate(dtformat))}()
    end

    for j in 1:length(types)
        if types[j] <: Characters
            if types[j] == Characters
                types[j] = Characters{8, UInt8}
                push!(charbuff, Vector{UInt8}(undef, 8))
            else
                types[j] = Characters{sizeof(types[j]), UInt8}
                push!(charbuff, Vector{UInt8}(undef, sizeof(types[j])))
            end
        elseif types[j] <: TimeType
            if dtformat isa Dict
                push!(dtfmt, _todate(get(dtformat, j, ISODateFormat)))
            else
                push!(dtfmt, _todate(dtformat))
            end
        end
    end
    charbuff = [deepcopy(charbuff) for _ in 1:nt]

    if !multiple_obs
        cz = div(fs-skip_bytes, nt)
        lo = [(i-1)*cz+skip_bytes+1 for i in 1:nt]
        hi = [i*cz+skip_bytes for i in 1:nt]
        hi[end] = fs
        # lo[1] += skip_bytes
        cur_value = Vector{UInt8}(undef, eol_len)
        for i in 1:length(hi)-1
            seek(f, hi[i] - eol_len)
            for k in 1:eol_len
                cur_value[k] = read(f, UInt8)
            end
            if last(cur_value) != eol_last || first(cur_value) != eol_first
                while true
                    if eol_len == 2
                        seek(f, position(f) - 1)
                    end
                    read!(f, cur_value)
                    last(cur_value) == eol_last && first(cur_value) == eol_first && break
                end
                hi[i] = position(f)
                lo[i+1] = hi[i] + 1
            end
        end

        ns = fill(1, nt)
        if nt > 1
            Threads.@threads for i in 1:nt
                ns[i] = count_lines_of_file(path, lo[i], hi[i], eol)
            end
        else
            for i in 1:nt
                ns[i] = count_lines_of_file(path, lo[i], hi[i], eol; limit = limit)
            end
        end

        line_hi = cumsum(ns)
        if issorted(line_hi)
            last_chunk_to_read = searchsortedfirst(line_hi, limit)
        else
            throw(ArgumentError("The current buffer size and line size are too small increase them by setting `lsize` and `buffsize`"))
        end
        res = Any[allocatecol_for_res(types[i], min(sum(ns), limit)) for i in 1:length(types)]

        line_lo = [1; line_hi[1:end] .+ 1]
        CLOSE(f)
        if nt > 1
            Threads.@threads for i in 1:min(nt, last_chunk_to_read)
                line_lo[i]>line_hi[i] && continue
                readfile_chunk!(res, line_lo[i], line_hi[i], charbuff[i], path, types, ns[i], lo[i], hi[i], colnames; delimiter = delimiter, linebreak = linebreak, lsize = lsize, buffsize = buffsize, fixed = colwidth, df = dtfmt, dlmstr = dlmstr, informat = informat, escapechar = escapechar, quotechar = quotechar, warn = warn, eolwarn = eolwarn, int_bases = int_bases, string_trim = string_trim, ignorerepeated = ignorerepeated, limit = limit)
            end
        else
            for i in 1:1
                readfile_chunk!(res, line_lo[i], line_hi[i], charbuff[i], path, types, ns[i], lo[i], hi[i], colnames; delimiter = delimiter, linebreak = linebreak, lsize = lsize, buffsize = buffsize, fixed = colwidth, df = dtfmt, dlmstr = dlmstr, informat = informat, escapechar = escapechar, quotechar = quotechar, warn = warn, eolwarn = eolwarn, int_bases = int_bases, string_trim = string_trim, ignorerepeated = ignorerepeated, limit = limit)
            end
        end
    else
        res = Any[allocatecol_for_res(types[i], 1) for i in 1:length(types)]
        lo = 0
        readfile_chunk!(res, 1, 1, charbuff[1], path, types, 1, 1, FILESIZE(path), colnames; delimiter = delimiter, linebreak = linebreak, lsize = lsize, buffsize = buffsize, fixed = colwidth, df = dtfmt, dlmstr = dlmstr, informat = informat, escapechar = escapechar, quotechar = quotechar, warn = warn, eolwarn = eolwarn, int_bases = int_bases, string_trim = string_trim, ignorerepeated = ignorerepeated, multiple_obs = true, limit = limit)
    end

    Dataset(res, colnames, copycols = false, makeunique = makeunique)
end

function guess_structure_of_delimited_file(path, delimiter; linebreak = nothing , header = true, guessingrows = 20, fixed = 0:0, dtformat = nothing, dlmstr = nothing, lsize = 2^15, buffsize = 2^16, informat = Dict{Int, Function}(), escapechar = nothing, quotechar = nothing, eolwarn = false, ignorerepeated = false, skipto = 1)

    if linebreak === nothing
        linebreak = (guess_eol_char(path))
    end
    eol = UInt8.(linebreak)
    eol_last = last(eol)
    eol_first = first(eol)
    eol_len = length(eol)
    if dlmstr === nothing
        dlm = UInt8.(delimiter)
    elseif dlmstr isa AbstractString
        dlm = collect(codeunits(dlmstr))
    else
        throw(ArgumentError("`dlmstr` must be a string of ASCII characters"))
    end
    read_byte_val = Vector{UInt8}(undef, eol_len)

    f = OUR_OPEN(path, read = true)

    f_pos = 0
    if skipto != 1
        _fsize_ = FILESIZE(f)
        l_length, f_pos = read_multiple_lines(path, f_pos, _fsize_, eol, skipto-1)
    else
        f_pos = 0
    end

    cnt = 1
    n_cols = 0
    col_width = 0:0

    if fixed isa Dict
        colwidth = UnitRange{Int}[]
    else
        colwidth = 0:0
    end
    l_length, _tmp_f_pos = read_one_line(path, f_pos+1, FILESIZE(f), eol)
    if l_length > lsize
        throw(ArgumentError("very wide delimited file! you need to set `lsize` and `buffsize` argument with larger values,  they are currently set as $lsize and $buffsize respectively. It is also recommended to use lower number of `guessingrows`"))
    end

    seek(f, f_pos)

    a_line_buff = Vector{UInt8}(undef, l_length)
    nb = readbytes!(f, a_line_buff)
    nb_pos = 0
    seen_dlm = true
    while true
        if fixed isa Dict
            col_width = get(fixed, n_cols + 1, 0:0)
        end
        if col_width != 0:0
            nb_pos = col_width.stop
            n_cols += 1
            seen_dlm = false
            push!(colwidth, col_width)
        else
            if quotechar !== nothing
                new_dlm_pos, new_lo, new_hi = find_next_delim(a_line_buff, nb_pos+1, l_length, dlm, dlmstr, quotechar, escapechar, ignorerepeated)
            else
                new_dlm_pos,new_lo, new_hi = find_next_delim(a_line_buff, nb_pos+1, l_length, dlm, dlmstr, ignorerepeated)
            end
            if new_dlm_pos > 0
                n_cols += 1
                nb_pos = new_dlm_pos
                seen_dlm = true
                fixed isa Dict && push!(colwidth, col_width)
            else
                if seen_dlm
                    n_cols += 1
                end
                fixed isa Dict && push!(colwidth, col_width)
                break
            end

        end
    end
    if header === true
        lo = _tmp_f_pos+1
    else
        lo = 0
    end
    seek(f, lo)
    rows_in = 0
    l_length = 1
    file_pos = max(0, lo - 1)
    while l_length>0
        l_length, file_pos = read_one_line(path, file_pos+1, FILESIZE(f), eol)
        if l_length > 0
            rows_in += 1
        end
        rows_in > guessingrows && break
    end
    hi = file_pos
    types = repeat([String], n_cols)
    res = [allocatecol_for_res(String, rows_in) for _ in 1:n_cols]
    readfile_chunk!(res, 1, rows_in, [], path, types, rows_in, lo, hi, nothing; delimiter = delimiter, linebreak = linebreak, buffsize = buffsize, fixed = colwidth, dlmstr = dlmstr, lsize = lsize, informat = informat, quotechar = quotechar, escapechar = escapechar, eolwarn = eolwarn, ignorerepeated = ignorerepeated)
    outtypes = Vector{DataType}(undef, n_cols)
    if !(dtformat isa Dict)
        for j in 1:n_cols
            old_cnt_vals = 0
            if all(isequal.(missing, res[j]))
                outtypes[j] = String
                continue
            end
            #TODO should we apply format here or in readfile_chunk! in the above?
            outtypes[j] = r_type_guess(res[j], identity)
        end
    else
        for j in 1:n_cols
            old_cnt_vals = 0
            if haskey(dtformat, j)
                hasdpart = occursin(r"\([ymd]", string(dtformat[j].tokens))
                hastpart = occursin(r"\([HMSs]", string(dtformat[j].tokens))
                if hasdpart && hastpart
                    outtypes[j] = DateTime
                elseif hasdpart
                    outtypes[j] = Date
                elseif hastpart
                    outtypes[j] = Time
                else
                    throw(ArgumentError("the date time format for column $(j) is not valid"))
                end
            else
                if all(isequal.(missing, res[j]))
                    outtypes[j] = String
                    continue
                end
                outtypes[j] = r_type_guess(res[j], identity)
            end
        end
    end
    CLOSE(f)
    linebreak, outtypes
end

function filereader(path; types = nothing, delimiter = ',', linebreak = nothing, header = true, threads::Bool = true, guessingrows::Int = 20, fixed = 0:0, buffsize::Int = 2^16, quotechar = nothing, escapechar = nothing, dtformat = dateformat"yyyy-mm-dd", dlmstr = nothing, lsize::Int = 2^15, informat = Dict{Int, Function}(), warn::Int = 20, eolwarn::Bool = true, emptycolname::Bool = false, int_base = Dict{Int, Tuple{DataType, Int}}(), string_trim::Bool = false, makeunique::Bool = false, ignorerepeated::Bool = false, multiple_obs::Bool = false, skipto::Int = 1, limit::Int = typemax(Int))
    supported_types = [Bool, Int8, Int16, Int32, Int64, Int8, UInt16, UInt32, UInt64, Float16, Float32, Float64, Int128, UInt128, BigFloat, String1, String3, String7, String15, String31, String63, String127, InlineString1, InlineString3, InlineString7, InlineString15, InlineString31, InlineString63, InlineString127,  Characters{1, UInt8},  Characters{2, UInt8}, Characters{3, UInt8}, Characters{4, UInt8}, Characters{5, UInt8}, Characters{6, UInt8}, Characters{7, UInt8}, Characters{8, UInt8}, Characters{9, UInt8}, Characters{10, UInt8}, Characters{11, UInt8}, Characters{12, UInt8}, Characters{13, UInt8}, Characters{14, UInt8}, Characters{15, UInt8}, Characters{16, UInt8},  Characters{1},  Characters{2}, Characters{3}, Characters{4}, Characters{5}, Characters{6}, Characters{7}, Characters{8}, Characters{9}, Characters{10}, Characters{11}, Characters{12}, Characters{13}, Characters{14}, Characters{15}, Characters{16}, TimeType, DateTime, Date, Time, String]

    lsize > buffsize && throw(ArgumentError("`lsize` must not be larger than `buffsize`"))
    ignorerepeated && dlmstr !== nothing && throw(ArgumentError("`ignorerepeated` option cannot be used when `dlmstr` is set"))

    guessingrows = min(guessingrows, limit)

    number_of_errors_happen_so_far[] = 0

    if quotechar !== nothing
        quotechar = UInt8(quotechar)
        if escapechar !== nothing
            escapechar = UInt8(escapechar)
        else
            escapechar = UInt8('\\')
        end
    end

    if multiple_obs
        !(types isa AbstractVector) && throw(ArgumentError("For reading multiple observations per line the types of each column must be specified"))
        fixed != 0:0 && throw(ArgumentError("When multiple observations per line is `true` the `fixed` keyword argument cannot be set"))
    end

    if types === nothing || types isa Dict
        linebreak, intypes = guess_structure_of_delimited_file(path, delimiter; linebreak = linebreak, header = header, guessingrows = guessingrows, fixed = fixed, buffsize = buffsize, dtformat = dtformat, dlmstr = dlmstr, lsize = lsize, informat = informat, escapechar = escapechar, quotechar = quotechar, eolwarn = false, ignorerepeated = ignorerepeated, skipto = skipto)
        if types isa Dict
            for (k, v) in types
                intypes[k] = v
            end
        end
    elseif types isa AbstractVector && eltype(types) <: Union{DataType,Type}
        intypes = types
        if linebreak === nothing
            linebreak = guess_eol_char(path)
        end
    else
        throw(ArgumentError("types should be a vector of types"))
    end
    # update intypes with information about integer base - by default we assume integer are base 10
    if !isempty(int_base)
        for (k, v) in int_base
            intypes[k] = v[1]
        end
        int_bases = fill(10, count(x -> x <: Integer, intypes))
        cnt = 1
        for (k, v) in int_base
            int_bases[cnt] = v[2]
            cnt += 1
        end
    else
        int_bases = nothing
    end
    !all(intypes .∈ Ref(Set(supported_types))) && throw(ArgumentError("DLMReaser only supports the following types(and their Subtypes): Bool, Integers, Floats, BigFloat, Characters, InlineStrings, TimeType, String"))
    !all(isascii.(delimiter)) && throw(ArgumentError("delimiter must be ASCII"))

    if multiple_obs
        if header == true
            throw(ArgumentError("`filereader` doesn't support reading header from file when `multiple_obs = true`"))
        end
        if skipto != 1
            throw(ArgumentError("`filereader` doesn't support `skipto` option when `multiple_obs = true`"))
        end
    end

    distribute_file(path, intypes; delimiter = delimiter, linebreak = linebreak, header = header, threads = threads, guessingrows = guessingrows, fixed = fixed, buffsize = buffsize, dtformat = dtformat, dlmstr = dlmstr, lsize = lsize, informat = informat, escapechar = escapechar, quotechar = quotechar, warn = warn, eolwarn = eolwarn, emptycolname = emptycolname, int_bases = int_bases, string_trim = string_trim, makeunique = makeunique, ignorerepeated = ignorerepeated, multiple_obs = multiple_obs, skipto = skipto, limit = limit)
end


text:
module ArrayMeta

using MacroTools

include("util.jl")
include("lowering.jl")
include("impl/abstract.jl")
# include("impl/dagger.jl")

end # module

ext: jl


text:
module WebAppRKHS

# Write your package code here.
import Genie

# include("./simple_rest.jl")
# include("./simple_query.jl")

end

max_stars_repo_name: RoyCCWang/WebAppRKHS

max_issues_repo_name: RoyCCWang/WebAppRKHS

max_forks_repo_name: RoyCCWang/WebAppRKHS


text:
module Statistics
  using ..DataParser
  files = ["extrastats"]
  for f in files
    include("$f.jl")
  end

end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

module VirtualNetworkGatewayClient

using Random
using Dates
using HTTP
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_LocalNetworkGatewaysApi.jl")
include("api_VirtualNetworkGatewayConnectionsApi.jl")
include("api_VirtualNetworkGatewaysApi.jl")

# export models
export convert, AddressSpace
export convert, BgpPeerStatus
export convert, BgpPeerStatusListResult
export convert, BgpSettings
export convert, ConnectionResetSharedKey
export convert, ConnectionSharedKey
export convert, GatewayRoute
export convert, GatewayRouteListResult
export convert, IpsecPolicy
export convert, LocalNetworkGatewayListResult
export convert, LocalNetworkGatewayPropertiesFormat
export convert, Resource
export convert, SubResource
export convert, TunnelConnectionHealth
export convert, VirtualNetworkGatewayConnectionListResult
export convert, VirtualNetworkGatewayConnectionPropertiesFormat
export convert, VirtualNetworkGatewayIPConfigurationPropertiesFormat
export convert, VirtualNetworkGatewayListResult
export convert, VirtualNetworkGatewayPropertiesFormat
export convert, VirtualNetworkGatewaySku
export convert, VpnClientConfiguration
export convert, VpnClientParameters
export convert, VpnClientRevokedCertificatePropertiesFormat
export convert, VpnClientRootCertificatePropertiesFormat
export convert, LocalNetworkGateway
export convert, VirtualNetworkGateway
export convert, VirtualNetworkGatewayConnection
export convert, VirtualNetworkGatewayIPConfiguration
export convert, VpnClientRevokedCertificate
export convert, VpnClientRootCertificate

# export operations
export convert, LocalNetworkGatewaysApi, VirtualNetworkGatewayConnectionsApi, VirtualNetworkGatewaysApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
# Plots target motion and sampled rang observations as circles

using SubmodularMaximization
using PyPlot

close("all")

steps = 100
grid_size = 20
num_robots = 3

grid = Grid(grid_size, grid_size)

robots = map(1:num_robots) do _
  random_state(grid)
end

sensor = RangingSensor()

initial = random_state(grid)
states = Array{typeof(initial)}(undef, steps)
states[1] = initial

plot_state_space(grid)
xlim([0, grid_size+1])
ylim([0, grid_size+1])

# precomputation
histogram_filter = Filter(grid)

for ii = 2:steps
  println("Step ", ii)
  state = target_dynamics(grid, states[ii-1])
  states[ii] = state

  range_observations = map(robots) do robot
    generate_observation(grid, sensor, robot, state)
  end

  # compute filter updates in place for this script
  process_update!(histogram_filter, transition_matrix(grid))
  for (robot, observation) in zip(robots, range_observations)
    measurement_update!(histogram_filter, robot, get_states(grid), sensor, grid,
                        observation)
  end

  plots=[]
  for (robot, observation) in zip(robots, range_observations)
    append!(plots, plot_robot(robot))
    append!(plots, plot_observation(robot, observation))
  end
  append!(plots, plot_trajectory(states[1:ii]))
  append!(plots, visualize_filter(histogram_filter))

  readline()

  if ii < steps
    foreach(x->x.remove(), plots)
  end
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
# startup.jl
#
# Author: Tom Price
# Date: Dec 2018
#
# This script starts up the chaipcr development environment
# in the Julia REPL from the bioinformatics/test directory

# NB.
# The following are not implemented yet:
# ensure_ci in shared.jl
# get_mysql_data_well in shared.jl

# ;cd ~/chaipcr/bioinformatics/QpcrAnalysis
# ;sudo mount -t vboxsf shared /mnt/share
#
push!(LOAD_PATH,pwd()*"/../QpcrAnalysis")
LOAD_FROM_DIR=pwd()*"/../QpcrAnalysis"

using QpcrAnalysis

import DataStructures.OrderedDict
import JuMP: Model, @variable, @objective, @NLobjective, @constraint, @NLconstraint,
    solve, getvalue, getobjectivevalue
import Ipopt.IpoptSolver
import DataArrays.DataArray
import Clustering: ClusteringResult, kmeans!, kmedoids!, silhouettes
import Combinatorics.combinations
import JSON, BSON, JLD.load
import Dierckx: Spline1D, derivative
import DataFrames: DataFrame, by

# development & testing
import Base.Test
import FactCheck: facts, context, @fact, clear_results, exitstatus, less_than_or_equal
FactCheck.clear_results()

include("../QpcrAnalysis/shared.jl")

# dispatch
include("../QpcrAnalysis/action_types.jl")
include("../QpcrAnalysis/dispatch.jl")

# data format verification
include("../QpcrAnalysis/verify_request.jl")
include("../QpcrAnalysis/verify_response.jl")

# calibration

include("../QpcrAnalysis/deconv.jl")
include("../QpcrAnalysis/adj_w2wvaf.jl")
include("../QpcrAnalysis/calib.jl")

# amplification
include("../QpcrAnalysis/amp_models/types_for_amp_models.jl")
include("../QpcrAnalysis/amp_models/sfc_models.jl")
include("../QpcrAnalysis/amp_models/MAKx.jl")
include("../QpcrAnalysis/amp_models/MAKERGAUL.jl")
include("../QpcrAnalysis/types_for_allelic_discrimination.jl")
include("../QpcrAnalysis/amp.jl")
include("../QpcrAnalysis/allelic_discrimination.jl")

#include("../QpcrAnalysis/standard_curve.jl")

# melt standard_curve
include("../QpcrAnalysis/multi_channel.jl")
include("../QpcrAnalysis/supsmu.jl")
include("../QpcrAnalysis/meltcrv.jl")

# analyze_customized
include("../QpcrAnalysis/analyze_customized/thermal_consistency.jl")
include("../QpcrAnalysis/analyze_customized/thermal_performance_diagnostic.jl")
include("../QpcrAnalysis/analyze_customized/optical_test_single_channel.jl")
include("../QpcrAnalysis/analyze_customized/optical_test_dual_channel.jl")
include("../QpcrAnalysis/analyze_customized/optical_cal.jl")

k = JLD.load("$LOAD_FROM_DIR/k4dcv_ip84_calib79n80n81_vec.jld")["k4dcv"]
kk = K4Deconv(k.k_s,k.k_inv_vec,k.inv_note)
const K4DCV = kk

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
using Cfdsof
using Test

@testset "Cfdsof.jl" begin
    # Write your tests here.
end


text:
"""
    GXBeamAssembly

Model which describes the behavior of an assembly of beam elements, as modeled by the
`GXBeam` package.  State variables are as defined by GXBeam.  Inputs correspond to the
external forces ``F_{x,i}, F_{y,i}, F_{z,i}, M_{x,i}, M_{y,i}, M_{z,i}`` or
displacements ``u_{x,i}, u_{y,i}, u_{z,i}, \\theta_{x,i}, \\theta_{y,i},
\\theta_{z,i}`` applied to each node, followed by the distributed loads
``f_{x,i}, f_{y,i}, f_{z,i}, m_{x,i}, m_{y,i}, m_{z,i}`` applied to each beam
element, followed by the point mass properties ``m, p, I_{11}, I_{22}, I_{33}, I_{12},
I_{13}, I_{23}`` corresponding to each beam element, followed by the gravity vector, linear
velocity, angular velocity, linear acceleration, and angular acceleration of the system.
Parameters correspond to the location ``p_{x}, p_{y}, p_{z}`` of each node followed by each
beam element's properties. Each beam element's properties are defined by a triad which defines
the orientation of the beam element ``e_{1,x}, e_{1,y}, e_{1,z}, e_{2,x}, e_{2,y}, e_{2,z},
e_{3,x}, e_{3,y}, e_{3,z}``, followed by the 21 independent entries of the compliance matrix
``C_{11}, C_{12}, C_{13}, C_{14}, C_{15}, C_{16}, C_{22}, C_{23}, C_{24}, C_{25}, C_{26},
C_{33}, C_{34}, C_{35}, C_{36}, C_{44}, C_{45}, C_{46}, C_{55}, C_{56}, C_{66}``, followed
by the beam element's inertial properties ``\\mu, x_{m,2}, x_{m,3}, i_{22}, i_{33}, i_{23}``
"""
struct GXBeamAssembly{TF,TI,TC,TB}
    # beam element connectivity
    start::TC
    stop::TC
    # location and degree of freedom of displacement constraints
    displacement::TB
    # scaling parameters for the system
    force_scaling::TF
    # indices for accessing governing equations for each point and beam element
    irow_point::TI
    irow_elem::TI
    irow_elem1::TI
    irow_elem2::TI
    # indices for accessing state variables for each point and beam element
    icol_point::TI
    icol_elem::TI
end

# --- Constructors --- #

"""
    GXBeamAssembly(start, stop, displacement; kwargs...)

Construct a geometrically exact beam theory structural model with beam elements
which extend from the point indices in `start` to the point indices in `stop` and
points with prescribed displacements as specified in `displacement`.

# Arguments
 - `start`: Vector containing point index where each beam element starts
 - `stop`: Vector containing point index where each beam element stops
 - `displacement`: Boolean matrix indicating the point index and degree of
    freedom of the system's displacement constraints.  Rows correspond to the
    degrees of freedom ``x, y, z, \\theta_x, \\theta_y, \\theta_z`` and
    columns correspond to the point indices ``1, 2, \\dots, N_p`` where ``N_p``
    is the total number of points.

# Keyword Arguments
 - `force_scaling = 1.0`: Factor used to scale forces and moments internally
"""
function GXBeamAssembly(start, stop, displacement; force_scaling=1.0)

    static = false

    N, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem =
        GXBeam.system_indices(start, stop, static)

    return GXBeamAssembly(start, stop, displacement, force_scaling,
        irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)
end

"""
    GXBeamAssembly(assembly, prescribed_conditions; kwargs...)

Construct a geometrically exact beam theory structural model with connectivity
as specified in `assembly` and displacement constraints as specified in
`prescribed_conditions`.

# Keyword Arguments
 - `force_scaling = 1.0`: Factor used to scale forces and moments internally.  If
    not specified, a suitable default will be chosen based on the entries of the
    compliance matrix.
"""
function GXBeamAssembly(assembly, prescribed_conditions;
    force_scaling = GXBeam.default_force_scaling(assembly))

    npoint = length(assembly.points)
    nelem = length(assembly.elements)

    start = assembly.start
    stop = assembly.stop

    static = false

    N, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem =
        GXBeam.system_indices(start, stop, static)

    displacement = zeros(Bool, 6, npoint)
    for key in keys(prescribed_conditions)
        displacement[:,key] .= prescribed_conditions[key].isforce .== false
    end

    return GXBeamAssembly(start, stop, displacement, force_scaling,
        irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)
end

# --- Submodel Creation --- #

function Submodel(model::GXBeamAssembly)

    # unpack constants
    start = model.start
    stop = model.stop
    displacement = model.displacement
    force_scaling = model.force_scaling
    irow_point = model.irow_point
    irow_elem = model.irow_elem
    irow_elem1 = model.irow_elem1
    irow_elem2 = model.irow_elem2
    icol_point = model.icol_point
    icol_elem = model.icol_elem

    # number of points and elements
    npoint = length(icol_point)
    nelem = length(icol_elem)

    # residual function
    fresid = (resid, dx, x, y, p, t) -> gxbeam_residual!(resid, dx, x, y, p, t;
        start, stop, displacement, force_scaling, irow_point, irow_elem,
        irow_elem1, irow_elem2, icol_point, icol_elem)

    # number of states, inputs, and parameters
    nx = 6*npoint + 18*nelem
    ny = 6*npoint + 6*nelem + 10*nelem + 15
    np = 3*npoint + 36*nelem

    # jacobian definitions
    ratejac = Linear((J, dx, x, y, p, t) -> gxbeam_rate_jacobian!(J, dx, x, y, p, t;
        start, stop, displacement, force_scaling, irow_point, irow_elem, irow_elem1,
        irow_elem2, icol_point, icol_elem))
    statejac = Nonlinear((J, dx, x, y, p, t) -> gxbeam_state_jacobian!(J, dx, x, y, p, t;
        start, stop, displacement, force_scaling, irow_point, irow_elem, irow_elem1,
        irow_elem2, icol_point, icol_elem))
    inputjac = Nonlinear() # TODO: define inputjac function
    paramjac = Nonlinear() # TODO: define paramjac function
    tgrad = Zeros()

    # convenience functions for setting states, inputs, and parameters
    setstate = (x; kwargs...) -> gxbeam_set_states!(x, displacement, icol_point, icol_elem,
        force_scaling; kwargs...)

    setinput = (y; kwargs...) -> gxbeam_set_inputs!(y, displacement, icol_point, icol_elem;
        kwargs...)

    setparam = (p; kwargs...) -> gxbeam_set_parameters!(p, icol_point, icol_elem; kwargs...)

    # convenience functions for separating states, inputs, and parameters
    sepstate = (x) -> gxbeam_separate_states(x, displacement, icol_point, icol_elem, force_scaling)
    sepinput = (y) -> gxbeam_separate_inputs(y, displacement, icol_point, icol_elem)
    sepparam = (p) -> gxbeam_separate_parameters(p, start, stop, icol_point, icol_elem)

    # model definition
    return Submodel{true}(fresid, nx, ny, np;
        ratejac = ratejac,
        statejac = statejac,
        inputjac = inputjac,
        paramjac = paramjac,
        tgrad = tgrad,
        setstate = setstate,
        setinput = setinput,
        setparam = setparam,
        sepstate = sepstate,
        sepinput = sepinput,
        sepparam = sepparam,
    )
end

# --- Internal Methods for this Model --- #

# residual function
function gxbeam_residual!(resid, dx, x, y, p, t; start, stop, displacement, force_scaling,
    irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)

    # number of points and elements
    np = length(icol_point)
    ne = length(icol_elem)

    # construct assembly from parameters
    assembly = gxbeam_assembly(p, np, ne, start, stop)

    # construct point loads, distributed loads, and point masses
    prescribed, distributed, point_masses = gxbeam_loads(y, np, ne, displacement, assembly.elements)

    # gravity vector
    gvec = SVector(y[end-14], y[end-13], y[end-12])

    # linear/angular velocity and acceleration
    x0 = @SVector zeros(3)
    v0 = SVector(y[end-11], y[end-10], y[end-9])
    ω0 = SVector(y[end-8], y[end-7], y[end-6])
    a0 = SVector(y[end-5], y[end-4], y[end-3])
    α0 = SVector(y[end-2], y[end-1], y[end])

    # return mass matrix multiplied state rates
    return GXBeam.dynamic_system_residual!(resid, dx, x, assembly, prescribed, distributed,
        point_masses, gvec, force_scaling, irow_point, irow_elem, irow_elem1, irow_elem2,
        icol_point, icol_elem, x0, v0, ω0, a0, α0)
end

# rate jacobian function
function gxbeam_rate_jacobian!(J, dx, x, y, p, t; start, stop, displacement, force_scaling,
    irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)

    # zero out jacobian entries
    J .= 0

    # number of points and elements
    np = length(icol_point)
    ne = length(icol_elem)

    # construct assembly from parameters
    assembly = gxbeam_assembly(p, np, ne, start, stop)

    # construct point loads, distributed loads, and point masses
    prescribed, distributed, point_masses = gxbeam_loads(y, np, ne, displacement, assembly.elements)

    # populate mass matrix entries
    return GXBeam.system_mass_matrix!(J, x, assembly, point_masses, force_scaling,
        irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)
end

# state jacobian function
function gxbeam_state_jacobian!(J, dx, x, y, p, t; start, stop, displacement, force_scaling,
    irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)

    # zero out jacobian entries
    J .= 0

    # number of points and elements
    np = length(icol_point)
    ne = length(icol_elem)

    # construct assembly from parameters
    assembly = gxbeam_assembly(p, np, ne, start, stop)

    # construct point and element loads from inputs
    prescribed, distributed, point_masses = gxbeam_loads(y, np, ne, displacement, assembly.elements)

    # gravity vector
    gvec = SVector(y[end-14], y[end-13], y[end-12])

    # linear/angular velocity and acceleration
    x0 = @SVector zeros(3)
    v0 = SVector(y[end-11], y[end-10], y[end-9])
    ω0 = SVector(y[end-8], y[end-7], y[end-6])
    a0 = SVector(y[end-5], y[end-4], y[end-3])
    α0 = SVector(y[end-2], y[end-1], y[end])

    # return jacobian of right hand side with respect to state variables
    return GXBeam.dynamic_system_jacobian!(J, dx, x, assembly, prescribed, distributed,
        point_masses, gvec, force_scaling, irow_point, irow_elem, irow_elem1, irow_elem2,
        icol_point, icol_elem, x0, v0, ω0, a0, α0)
end

# function for creating an assembly from the provided parameters
function gxbeam_assembly(p, np, ne, start, stop)
    points = [gxbeam_point(view(p, 3*(ip-1) + 1 : 3*(ip-1) + 3)) for ip = 1:np]
    elements = [gxbeam_element(view(p, 3*np + 36*(ie-1) + 1 : 3*np + 36*(ie-1) + 36),
        points, start[ie], stop[ie]) for ie = 1:ne]
    return GXBeam.Assembly(points, start, stop, elements)
end

# function for creating a point
gxbeam_point(p) = SVector(p[1], p[2], p[3])

# function for creating a beam element
function gxbeam_element(p, points, start, stop)
    # separate element parameters
    e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z, C11, C12, C13, C14, C15, C16,
        C22, C23, C24, C25, C26, C33, C34, C35, C36, C44, C45, C46, C55, C56,
        C66, μ, xm2, xm3, i22, i33, i23 = p
    # element length
    ΔL = norm(points[stop] - points[start])
    # element location
    x = (points[start] + points[stop])/2
    # element compliance matrix
    C = @SMatrix [
        C11 C12 C13 C14 C15 C16;
        C12 C22 C23 C24 C25 C26;
        C13 C23 C33 C34 C35 C36;
        C14 C24 C34 C44 C45 C46;
        C15 C25 C35 C45 C55 C56;
        C16 C26 C36 C46 C56 C66
        ]
    # element mass matrix
    mass = @SMatrix [
            μ       0     0       0 μ*xm3 -μ*xm2;
            0       μ     0  -μ*xm3     0      0;
            0       0     μ   μ*xm2     0      0;
            0  -μ*xm3 μ*xm2 i22+i33     0      0;
         μ*xm3      0     0       0   i22   -i23;
        -μ*xm2      0     0       0  -i23    i33;
    ]
    # element triad
    Cab = @SMatrix [
        e1x e2x e3x;
        e1y e2y e3y;
        e1z e2z e3z;
    ]

    return GXBeam.Element(ΔL, x, C, mass, Cab)
end

# define prescribed conditions, distributed loads, and point masses
function gxbeam_loads(y, np, ne, d, elements)
    prescribed = Dict(ip => gxbeam_prescribed_condition(view(y,
        6*(ip-1) + 1 : 6*(ip-1) + 6), view(d, :, ip)) for ip = 1:np)
    distributed = Dict(ie => gxbeam_distributed_load(view(y,
        6*np + 6*(ie-1) + 1 : 6*np + 6*(ie-1) + 6), elements[ie].L) for ie = 1:ne)
    point_masses = Dict(ie => gxbeam_point_mass(view(y,
        6*np + 6*ne + 10*(ie-1) + 1 : 6*np + 6*ne + 10*(ie-1) + 10)) for ie = 1:ne)
    return prescribed, distributed, point_masses
end

# define prescribed conditions at a point
function gxbeam_prescribed_condition(y, d)
    displacement = SVector{6}(d)
    force = displacement .== false
    value = SVector{6}(y)
    follower = @SVector zeros(6)
    return PrescribedConditions(force, value, follower)
end

# define distributed loads on an element
function gxbeam_distributed_load(y, ΔL)
    # prescribed distributed loads
    f = SVector(y[1], y[2], y[3])
    m = SVector(y[4], y[5], y[6])
    # non-follower distributed loads
    f1 = f2 = ΔL*f/2
    m1 = m2 = ΔL*m/2
    # follower distributed loads
    f1_follower = f2_follower = @SVector zeros(3)
    m1_follower = m2_follower = @SVector zeros(3)
    # distributed loads
    return DistributedLoads(f1, f2, m1, m2, f1_follower, f2_follower, m1_follower,
        m2_follower)
end

# define point mass attached to an element
function gxbeam_point_mass(y)
    # point mass
    m = y[1]
    # point location
    p = SVector(y[2], y[3], y[4])
    # point inertia
    I11 = y[5]
    I22 = y[6]
    I33 = y[7]
    I12 = y[8]
    I13 = y[9]
    I23 = y[10]
    I = @SMatrix [I11 I12 I13; I12 I22 I23; I13 I23 I33]
    # point mass
    return PointMass(m, p, I)
end

# convenience function for defining this model's state vector
function gxbeam_set_states!(x, displacement, icol_point, icol_elem, force_scaling;
    u_e = nothing, theta_e = nothing, F_e = nothing, M_e = nothing,
    V_e = nothing, Omega_e = nothing, u_p = nothing, theta_p = nothing,
    F_p = nothing, M_p = nothing)

    nelem = length(icol_elem)
    npoint = length(icol_point)

    for ielem = 1:nelem
        icol = icol_elem[ielem]
        if !isnothing(u_e)
            GXBeam.set_element_deflection!(x, icol, u_e[ielem])
        end
        if !isnothing(theta_e)
            GXBeam.set_element_rotation!(x, icol, theta_e[ielem])
        end
        if !isnothing(F_e)
            GXBeam.set_element_forces!(x, icol, F_e[ielem], force_scaling)
        end
        if !isnothing(M_e)
            GXBeam.set_element_moments!(x, icol, M_e[ielem], force_scaling)
        end
        if !isnothing(V_e)
            GXBeam.set_element_linear_velocity!(x, icol, V_e[ielem])
        end
        if !isnothing(Omega_e)
            GXBeam.set_element_angular_velocity!(x, icol, Omega_e[ielem])
        end
    end

    for ipoint = 1:npoint
        icol = icol_point[ipoint]
        prescribed_forces = SVector{6}(view(displacement, :, ipoint)) .== false
        if !isnothing(u_p)
            GXBeam.set_point_deflections!(x, icol, u_p[ipoint], prescribed_forces)
        end
        if !isnothing(theta_p)
            GXBeam.set_point_rotations!(x, icol, theta_p[ipoint], prescribed_forces)
        end
        if !isnothing(F_p)
            GXBeam.set_point_forces!(x, icol, F_p[ipoint], prescribed_forces, force_scaling)
        end
        if !isnothing(M_p)
            GXBeam.set_point_moments!(x, icol, M_p[ipoint], prescribed_forces, force_scaling)
        end
    end

    return x
end

# convenience function for defining this model's input vector
function gxbeam_set_inputs!(y, displacement, icol_point, icol_elem;
    u_p = nothing, theta_p = nothing, F_p = nothing, M_p = nothing,
    f_d = nothing, m_d = nothing,
    m_m = nothing, p_m = nothing, I_m = nothing,
    gvec = nothing, v = nothing, omega = nothing, a = nothing, alpha = nothing)

    np = length(icol_point)
    ne = length(icol_elem)

    for ip = 1:np
        prescribed_forces = SVector{6}(view(displacement, :, ip)) .== false
        for i = 1:3
            if prescribed_forces[i]
                if !isnothing(F_p)
                    y[6*(ip-1) + i] = F_p[ip][i]
                end
            else
                if !isnothing(u_p)
                    y[6*(ip-1) + i] = u_p[ip][i]
                end
            end
        end
        for i = 1:3
            if prescribed_forces[i]
                if !isnothing(M_p)
                    y[6*(ip-1) + 3 + i] = M_p[ip][i]
                end
            else
                if !isnothing(theta_p)
                    y[6*(ip-1) + 3 + i] = theta_p[ip][i]
                end
            end
        end
    end

    for ie = 1:ne
        if !isnothing(f_d)
            for i = 1:3
                y[6*np + 6*(ie-1) + i] = f_d[ie][i]
            end
        end
        if !isnothing(m_d)
            for i = 1:3
                y[6*np + 6*(ie-1) + 3 + i] = m_d[ie][i]
            end
        end
    end

    for ie = 1:ne
        if !isnothing(m_m)
            y[6*np + 6*ne + 10*(ie-1) + 1] = m_m[ie]
        end
        if !isnothing(p_m)
            y[6*np + 6*ne + 10*(ie-1) + 2] = p_m[ie][1]
            y[6*np + 6*ne + 10*(ie-1) + 3] = p_m[ie][2]
            y[6*np + 6*ne + 10*(ie-1) + 4] = p_m[ie][3]
        end
        if !isnothing(I_m)
            y[6*np + 6*ne + 10*(ie-1) + 5] = I_m[ie][1,1]
            y[6*np + 6*ne + 10*(ie-1) + 6] = I_m[ie][2,2]
            y[6*np + 6*ne + 10*(ie-1) + 7] = I_m[ie][3,3]
            y[6*np + 6*ne + 10*(ie-1) + 8] = I_m[ie][1,2]
            y[6*np + 6*ne + 10*(ie-1) + 9] = I_m[ie][1,3]
            y[6*np + 6*ne + 10*(ie-1) + 10] = I_m[ie][2,3]
        end
    end

    if !isnothing(gvec)
        y[6*np + 6*ne + 10*ne + 1] = gvec[1]
        y[6*np + 6*ne + 10*ne + 2] = gvec[2]
        y[6*np + 6*ne + 10*ne + 3] = gvec[3]
    end

    if !isnothing(v)
        y[6*np + 6*ne + 10*ne + 4] = v[1]
        y[6*np + 6*ne + 10*ne + 5] = v[2]
        y[6*np + 6*ne + 10*ne + 6] = v[3]
    end

    if !isnothing(omega)
        y[6*np + 6*ne + 10*ne + 7] = omega[1]
        y[6*np + 6*ne + 10*ne + 8] = omega[2]
        y[6*np + 6*ne + 10*ne + 9] = omega[3]
    end

    if !isnothing(a)
        y[6*np + 6*ne + 10*ne + 10] = a[1]
        y[6*np + 6*ne + 10*ne + 11] = a[2]
        y[6*np + 6*ne + 10*ne + 12] = a[3]
    end

    if !isnothing(alpha)
        y[6*np + 6*ne + 10*ne + 13] = alpha[1]
        y[6*np + 6*ne + 10*ne + 14] = alpha[2]
        y[6*np + 6*ne + 10*ne + 15] = alpha[3]
    end

    return y
end

# convenience function for defining this model's parameter vector
function gxbeam_set_parameters!(p, icol_point, icol_elem; assembly)

    np = length(icol_point)
    ne = length(icol_elem)

    for ip = 1:np
        p[3*(ip-1) + 1] = assembly.points[ip][1]
        p[3*(ip-1) + 2] = assembly.points[ip][2]
        p[3*(ip-1) + 3] = assembly.points[ip][3]
    end

    for ie = 1:ne
        element = assembly.elements[ie]
        p[3*np + 36*(ie-1) + 1] = element.Cab[1,1]
        p[3*np + 36*(ie-1) + 2] = element.Cab[2,1]
        p[3*np + 36*(ie-1) + 3] = element.Cab[3,1]
        p[3*np + 36*(ie-1) + 4] = element.Cab[1,2]
        p[3*np + 36*(ie-1) + 5] = element.Cab[2,2]
        p[3*np + 36*(ie-1) + 6] = element.Cab[3,2]
        p[3*np + 36*(ie-1) + 7] = element.Cab[1,3]
        p[3*np + 36*(ie-1) + 8] = element.Cab[2,3]
        p[3*np + 36*(ie-1) + 9] = element.Cab[3,3]
        p[3*np + 36*(ie-1) + 10] = element.compliance[1,1]
        p[3*np + 36*(ie-1) + 11] = element.compliance[1,2]
        p[3*np + 36*(ie-1) + 12] = element.compliance[1,3]
        p[3*np + 36*(ie-1) + 13] = element.compliance[1,4]
        p[3*np + 36*(ie-1) + 14] = element.compliance[1,5]
        p[3*np + 36*(ie-1) + 15] = element.compliance[1,6]
        p[3*np + 36*(ie-1) + 16] = element.compliance[2,2]
        p[3*np + 36*(ie-1) + 17] = element.compliance[2,3]
        p[3*np + 36*(ie-1) + 18] = element.compliance[2,4]
        p[3*np + 36*(ie-1) + 19] = element.compliance[2,5]
        p[3*np + 36*(ie-1) + 20] = element.compliance[2,6]
        p[3*np + 36*(ie-1) + 21] = element.compliance[3,3]
        p[3*np + 36*(ie-1) + 22] = element.compliance[3,4]
        p[3*np + 36*(ie-1) + 23] = element.compliance[3,5]
        p[3*np + 36*(ie-1) + 24] = element.compliance[3,6]
        p[3*np + 36*(ie-1) + 25] = element.compliance[4,4]
        p[3*np + 36*(ie-1) + 26] = element.compliance[4,5]
        p[3*np + 36*(ie-1) + 27] = element.compliance[4,6]
        p[3*np + 36*(ie-1) + 28] = element.compliance[5,5]
        p[3*np + 36*(ie-1) + 29] = element.compliance[5,6]
        p[3*np + 36*(ie-1) + 30] = element.compliance[6,6]
        p[3*np + 36*(ie-1) + 31] =  element.mass[1,1] # μ
        p[3*np + 36*(ie-1) + 32] =  element.mass[3,4]/element.mass[1,1] # xm2
        p[3*np + 36*(ie-1) + 33] =  element.mass[1,5]/element.mass[1,1] # xm3
        p[3*np + 36*(ie-1) + 34] =  element.mass[5,5] # i22
        p[3*np + 36*(ie-1) + 35] =  element.mass[6,6] # i33
        p[3*np + 36*(ie-1) + 36] = -element.mass[5,6] # i23
    end

    return p
end

# convenience function for separating this model's state vector
function gxbeam_separate_states(x, displacement, icol_point, icol_elem, force_scaling)

    # floating point type
    TF = eltype(x)

    # number of points and elements
    np = length(icol_point)
    ne = length(icol_elem)

    # separate point state variables
    u_p = Vector{SVector{3,TF}}(undef, np)
    theta_p = Vector{SVector{3,TF}}(undef, np)
    F_p = Vector{SVector{3,TF}}(undef, np)
    M_p = Vector{SVector{3,TF}}(undef, np)

    for ip = 1:np

        icol = icol_point[ip]

        prescribed_forces = SVector{6}(view(displacement, :, ip)) .== false

        # get the displacement and rotations of the point
        u_p[ip] = SVector(ifelse(prescribed_forces[1], x[icol  ], NaN),
                    ifelse(prescribed_forces[2], x[icol+1], NaN),
                    ifelse(prescribed_forces[3], x[icol+2], NaN))
        theta_p[ip] = SVector(ifelse(prescribed_forces[4], x[icol+3], NaN),
                    ifelse(prescribed_forces[5], x[icol+4], NaN),
                    ifelse(prescribed_forces[6], x[icol+5], NaN))

        # overwrite external forces/moments with solved for forces/moments
        F_p[ip] = SVector(ifelse(prescribed_forces[1], NaN, x[icol  ] * force_scaling),
                    ifelse(prescribed_forces[2], NaN, x[icol+1] * force_scaling),
                    ifelse(prescribed_forces[3], NaN, x[icol+2] * force_scaling))
        M_p[ip] = SVector(ifelse(prescribed_forces[4], NaN, x[icol+3] * force_scaling),
                    ifelse(prescribed_forces[5], NaN, x[icol+4] * force_scaling),
                    ifelse(prescribed_forces[6], NaN, x[icol+5] * force_scaling))

        # convert rotation parameter to Wiener-Milenkovic parameters
        scaling = GXBeam.rotation_parameter_scaling(theta_p[ip])
        theta_p[ip] *= scaling

    end

    # element state variables
    u_e = Vector{SVector{3,TF}}(undef, ne)
    theta_e = Vector{SVector{3,TF}}(undef, ne)
    F_e = Vector{SVector{3,TF}}(undef, ne)
    M_e = Vector{SVector{3,TF}}(undef, ne)
    V_e = Vector{SVector{3,TF}}(undef, ne)
    Omega_e = Vector{SVector{3,TF}}(undef, ne)

    for ie = 1:ne

        icol = icol_elem[ie]

        u_e[ie] = SVector(x[icol], x[icol+1], x[icol+2])
        theta_e[ie] = SVector(x[icol+3], x[icol+4], x[icol+5])
        F_e[ie] = SVector(x[icol+6], x[icol+7], x[icol+8]) .* force_scaling
        M_e[ie] = SVector(x[icol+9], x[icol+10], x[icol+11]) .* force_scaling
        V_e[ie] = SVector(x[icol+12], x[icol+13], x[icol+14])
        Omega_e[ie] = SVector(x[icol+15], x[icol+16], x[icol+17])

        # convert rotation parameter to Wiener-Milenkovic parameters
        scaling = GXBeam.rotation_parameter_scaling(theta_e[ie])
        theta_e[ie] *= scaling

    end

    return (u_p = u_p, theta_p = theta_p, F_p = F_p, M_p = M_p,
        u_e = u_e, theta_e = theta_e, F_e = F_e, M_e = M_e, V_e = V_e, Omega_e = Omega_e)
end

# convenience function for separating this model's input vector
function gxbeam_separate_inputs(y, displacement, icol_point, icol_elem)

    TF = eltype(y)

    np = length(icol_point)
    ne = length(icol_elem)

    # separate prescribed conditions
    u_p = Vector{SVector{3,TF}}(undef, np)
    theta_p = Vector{SVector{3,TF}}(undef, np)
    F_p = Vector{SVector{3,TF}}(undef, np)
    M_p = Vector{SVector{3,TF}}(undef, np)
    for ip = 1:np
        prescribed_forces = SVector{6}(view(displacement, :, ip)) .== false
        F_p[ip] = SVector(
            ifelse(prescribed_forces[1], y[6*(ip-1) + 1], NaN),
            ifelse(prescribed_forces[2], y[6*(ip-1) + 2], NaN),
            ifelse(prescribed_forces[3], y[6*(ip-1) + 3], NaN))
        M_p[ip] = SVector(
            ifelse(prescribed_forces[4], y[6*(ip-1) + 4], NaN),
            ifelse(prescribed_forces[5], y[6*(ip-1) + 5], NaN),
            ifelse(prescribed_forces[6], y[6*(ip-1) + 6], NaN))
        u_p[ip] = SVector(
            ifelse(prescribed_forces[1], NaN, y[6*(ip-1) + 1]),
            ifelse(prescribed_forces[2], NaN, y[6*(ip-1) + 2]),
            ifelse(prescribed_forces[3], NaN, y[6*(ip-1) + 3]))
        theta_p[ip] = SVector(
            ifelse(prescribed_forces[4], NaN, y[6*(ip-1) + 4]),
            ifelse(prescribed_forces[5], NaN, y[6*(ip-1) + 5]),
            ifelse(prescribed_forces[6], NaN, y[6*(ip-1) + 6]))
    end

    # separate distributed loads
    f_d = Vector{SVector{3,TF}}(undef, ne)
    m_d = Vector{SVector{3,TF}}(undef, ne)
    for ie = 1:ne
        f_d[ie] = SVector(
            y[6*np + 6*(ie-1) + 1],
            y[6*np + 6*(ie-1) + 2],
            y[6*np + 6*(ie-1) + 3])
        m_d[ie] = SVector(
            y[6*np + 6*(ie-1) + 4],
            y[6*np + 6*(ie-1) + 5],
            y[6*np + 6*(ie-1) + 6])
    end

    # separate point masses
    m_m = Vector{TF}(undef, ne)
    p_m = Vector{SVector{3,TF}}(undef, ne)
    I_m = Vector{SMatrix{3,3,TF}}(undef, ne)
    for ie = 1:ne
        m_m[ie] = y[6*np + 6*ne + 10*(ie-1) + 1]
        p_m[ie] = SVector(
            y[6*np + 6*ne + 10*(ie-1) + 2],
            y[6*np + 6*ne + 10*(ie-1) + 3],
            y[6*np + 6*ne + 10*(ie-1) + 4])
        I11 = y[6*np + 6*ne + 10*(ie-1) + 5]
        I22 = y[6*np + 6*ne + 10*(ie-1) + 6]
        I33 = y[6*np + 6*ne + 10*(ie-1) + 7]
        I12 = y[6*np + 6*ne + 10*(ie-1) + 8]
        I13 = y[6*np + 6*ne + 10*(ie-1) + 9]
        I23 = y[6*np + 6*ne + 10*(ie-1) + 10]
        I_m[ie] = @SMatrix [I11 I12 I13; I12 I22 I23; I13 I23 I33]
    end

    # gravity
    gvec =  view(y, 6*np + 6*ne + 10*ne +  1 : 6*np + 6*ne + 10*ne +  3)

    # linear/angular velocity
    v =     view(y, 6*np + 6*ne + 10*ne +  4 : 6*np + 6*ne + 10*ne +  6)
    omega = view(y, 6*np + 6*ne + 10*ne +  7 : 6*np + 6*ne + 10*ne +  9)

    # linear/angular acceleration
    a =     view(y, 6*np + 6*ne + 10*ne + 10 : 6*np + 6*ne + 10*ne + 12)
    alpha = view(y, 6*np + 6*ne + 10*ne + 13 : 6*np + 6*ne + 10*ne + 15)

    return (u_p=u_p, theta_p=theta_p, F_p=F_p, M_p=M_p,
        f_d=f_d, m_d=m_d, m_m=m_m, p_m=p_m, I_m=I_m,
        gvec=gvec, v=v, omega=omega, a=a, alpha=alpha)
end

# convenience function for separating this model's parameter vector
function gxbeam_separate_parameters(p, start, stop, icol_point, icol_elem)
    np = length(icol_point)
    ne = length(icol_elem)
    return (assembly = gxbeam_assembly(p, np, ne, start, stop),)
end

# --- Internal Methods for Couplings with this Model --- #

# transformation from body to local frame
gxbeam_body_to_element(elem, θ) = elem.Cab'*GXBeam.get_C(θ)

# local element velocities (in local frame)
function gxbeam_element_velocities(V, Ω)
    vi = V
    ωi = Ω
    return vi, ωi
end

# local element accelerations (in local frame)
function gxbeam_element_accelerations(elem, u, θ, Vdot, Ωdot, a0, α0)
    x0 = @SVector zeros(3)
    C = gxbeam_body_to_element(elem, θ)
    ai = C*(a0 + cross(α0, elem.x - x0) + cross(α0, u)) + Vdot
    αi = C*α0 + Ωdot
    return ai, αi
end


text:
AXIS = ["y","z"]
NEV = 4

gx,gy,gz = 3,2,1
αx,αy,αz = 1,1,1
Cx,Cy,Cz = (1,1,1).*100

hw = (sqrt(gx/Cx),sqrt(gy/Cy),sqrt(gz/Cz))
I = (αx/sqrt(2)/(gx*Cx)^(1/4),αy/sqrt(2)/(gy*Cy)^(1/4),αz/sqrt(2)/(gz*Cz)^(1/4))

ϕ = linspace(0,1*5.4,100)

@everywhere N = (27,25,51)
@everywhere include("junction-hamiltonian.jl")

@everywhere function energies(hw,I,phi,N,NEV)
    println("phi = $phi")
    run(`hostname`)
    BLAS.set_num_threads(1)
    p = Parameters(hw,I,phi)
    @time Hm = FillHamiltonian(N...,p)
    @time solution = eigs(Hm,nev=NEV,tol=1e-18,maxiter=1000,which=:SR)

    spinx = []
    spiny = []
    spinz = []
    for i in 1:length(solution[1])
        up=solution[2][:,i][1:2:end]
        down=solution[2][:,i][2:2:end]

        spinxi = real(up'*down + down'*up)
        spinyi = real(-1im*up'*down + 1im*down'*up)
        spinzi = up'*up - down'*down

        push!(spinx,spinxi)
        push!(spiny,spinyi)
        push!(spinz,spinzi)
    end

    phix = []
    phiy = []
    phiz = []
    for i in 1:length(solution[1])
        phixi = real(solution[2][:,i]'*Phix*solution[2][:,i])
        phiyi = real(solution[2][:,i]'*Phiy*solution[2][:,i])
        phizi = real(solution[2][:,i]'*Phiz*solution[2][:,i])

        push!(phix,phixi)
        push!(phiy,phiyi)
        push!(phiz,phizi)
    end

    return Dict(:EVALS=>real(solution[1]),:spinx=>spinx,:spiny=>spiny,:spinz=>spinz,:phix=>phix,:phiy=>phiy,:phiz=>phiz) #real(solution[1]),spin
end

knotsx = Iterators.product((ϕ,0:0,0:0)...)
knotsy = Iterators.product((0:0,ϕ,0:0)...)
knotsz = Iterators.product((0:0,0:0,ϕ)...)
knotsxy = zip(((ϕ, ϕ, zeros(ϕ)) ./ sqrt(2))...)
knotsyz = zip(((zeros(ϕ), ϕ, ϕ) ./ sqrt(2))...)
knotszx = zip(((ϕ, zeros(ϕ), ϕ) ./ sqrt(2))...)

knots = []
for i in AXIS
    vi = eval(parse("knots$i"))
    push!(knots,collect(vi))
end

knots = vcat([i[:] for i in knots]...)
#knots = vcat(knots...)

#knots = [knotsx;knotsy;knotsz;knotsxy;knotsyz;knotszx]
results = pmap(x->energies(hw,I,x,N,NEV),knots)

parlist = []

Lright=0
for i in AXIS
    push!(parlist,"knots$i")
    push!(parlist,eval(parse("knots$i")))
    push!(parlist,"results$i")

    Li = length(eval(parse("knots$i")))
    Lright += Li
    Lleft = 1+Lright-Li
    push!(parlist,results[Lleft:Lright])
end

using JLD
save("outdir/energies-$N $hw $I.jld","knots",knots,"results",results,parlist...)

max_stars_repo_name: akels/WeylDisc

max_issues_repo_name: akels/WeylDisc

max_forks_repo_name: akels/WeylDisc


text:
"""
    function turtle()

A mysterious turtle...

# Example
```jldoctest
julia> cowsay("Where is that pesky rabbit?", cow=Cowsay.turtle)
 _____________________________
< Where is that pesky rabbit? >
 -----------------------------
    \\                                  ___-------___
     \\                             _-~~             ~~-_
      \\                         _-~                    /~-_
             /^\\__/^\\         /~  \\                   /    \\
           /|  O|| O|        /      \\_______________/        \\
          | |___||__|      /       /                \\          \\
          |          \\    /      /                    \\          \\
          |   (_______) /______/                        \\_________ \\
          |         / /         \\                      /            \\
           \\         \\^\\\\         \\                  /               \\     /
             \\         ||           \\______________/      _-_       //\\__//
               \\       ||------_-~~-_ ------------- \\ --/~   ~\\    || __/
                 ~-----||====/~     |==================|       |/~~~~~
                  (_(__/  ./     /                    \\_\\      \\.
                         (_(___/                         \\_____)_)
```
"""
function turtle(;eyes="oo", tongue="  ", thoughts="\\")
the_cow =
"""
    $thoughts                                  ___-------___
     $thoughts                             _-~~             ~~-_
      $thoughts                         _-~                    /~-_
             /^\\__/^\\         /~  \\                   /    \\
           /|  O|| O|        /      \\_______________/        \\
          | |___||__|      /       /                \\          \\
          |          \\    /      /                    \\          \\
          |   (_______) /______/                        \\_________ \\
          |         / /         \\                      /            \\
           \\         \\^\\\\         \\                  /               \\     /
             \\         ||           \\______________/      _-_       //\\__//
               \\       ||------_-~~-_ ------------- \\ --/~   ~\\    || __/
                 ~-----||====/~     |==================|       |/~~~~~
                  (_(__/  ./     /                    \\_\\      \\.
                         (_(___/                         \\_____)_)
"""
return the_cow
end


text:
using AutoRisk

function analyze_risk_dataset(output_filepath)
    dataset = h5open(output_filepath)
    features = read(dataset["risk/features"])
    targets = read(dataset["risk/targets"])
    println("avg features: $(mean(features, (2,3)))")
    avg_targets = mean(clamp(sum(targets, 2), 0, 1), 3)
    println("avg targets: $(avg_targets)")
    println("size of dataset features: $(size(features))")
    println("size of dataset targets: $(size(targets))")
    if exists(dataset, "risk/weights")
        weights = read(dataset["risk/weights"])
        inds = find(weights .!= 1.)

        if length(inds) > 0
                avg_prop_weight = mean(weights[1, inds])
                println("avg proposal weight: $(avg_prop_weight)")
                med_prop_weight = median(weights[1, inds])
                println("median proposal weight: $(med_prop_weight)")
        end
    end
end

function build_dataset_collector(
        output_filepath::String,
        col_id::Int = 0,
        num_scenarios = 1;
        num_lanes = 1,
        min_num_veh = 30,
        max_num_veh = 30,
        min_base_speed = 12.,
        max_base_speed = 12.,
        min_vehicle_length = 3.5,
        max_vehicle_length = 3.5,
        min_vehicle_width = 1.5,
        max_vehicle_width = 1.5,
        min_init_dist = 10.,
        lon_accel_std_dev = 0.5,
        lat_accel_std_dev = 0.1,
        overall_response_time = 0.,
        err_p_a_to_i = 0.,
        err_p_i_to_a = 1.,
        sampling_time = .1,
        sampling_period = .1,
        num_runs = 1,
        veh_idx_can_change = false,
        feature_step_size = 5,
        feature_timesteps = 50,
        chunk_dim = 1,
        roadway_length = 120.,
        roadway_radius = 30.,
        burn_in_time = 60.
    )

    # roadway generator
    roadway = gen_stadium_roadway(num_lanes, length = roadway_length, radius = roadway_radius)
    roadway_gen = StaticRoadwayGenerator(roadway)

    # scene gen
    scene = Scene(max_num_veh)
    scene_gen = HeuristicSceneGenerator(
        min_num_veh,
        max_num_veh,
        min_base_speed,
        max_base_speed,
        min_vehicle_length,
        max_vehicle_length,
        min_vehicle_width,
        max_vehicle_width,
        min_init_dist
    )

    # beh generator
    passive = get_passive_behavior_params(
                lon_σ = lon_accel_std_dev,
                lat_σ = lat_accel_std_dev,
                overall_response_time = overall_response_time,
                err_p_a_to_i = err_p_a_to_i,
                err_p_i_to_a = err_p_i_to_a)
    aggressive = get_aggressive_behavior_params(
                lon_σ = lon_accel_std_dev,
                lat_σ = lat_accel_std_dev,
                overall_response_time = overall_response_time,
                err_p_a_to_i = err_p_a_to_i,
                err_p_i_to_a = err_p_i_to_a)

    params = [aggressive, passive]
    weights = StatsBase.Weights([.5,.5])
    behavior_gen = PredefinedBehaviorGenerator(params, weights)
    gen = FactoredGenerator(roadway_gen, scene_gen, behavior_gen)

    subexts = [
        CoreFeatureExtractor(),
        TemporalFeatureExtractor(),
        CarLidarFeatureExtractor(2, extract_carlidar_rangerate = true),
        BehavioralFeatureExtractor()
    ]
    ext = MultiFeatureExtractor(subexts)

    feature_dim = length(ext)
    target_ext = TargetExtractor()
    target_dim = length(target_ext)
    target_timesteps = Int(ceil(sampling_time / sampling_period))

    prime_time = (feature_timesteps * feature_step_size) * .1 + burn_in_time
    max_num_scenes = Int(ceil((prime_time + sampling_time) / sampling_period))
    rec = SceneRecord(max_num_scenes, sampling_period, max_num_veh)
    features = Array{Float64}(feature_dim, feature_timesteps, max_num_veh)
    targets = Array{Float64}(target_dim, target_timesteps, max_num_veh)
    agg_targets = Array{Float64}(target_dim, target_timesteps, max_num_veh)

    eval = MonteCarloEvaluator(
        ext,
        target_ext,
        num_runs,
        prime_time,
        sampling_time,
        veh_idx_can_change,
        rec,
        features,
        targets,
        agg_targets,
        feature_step_size = feature_step_size
    )

    max_num_samples = num_scenarios * max_num_veh
    feature_dim = length(ext)
    use_weights = typeof(weights) == Void ? false : true
    attrs = Dict()
    attrs["feature_names"] = feature_names(ext)
    attrs["target_names"] = feature_names(target_ext)
    dataset = Dataset(
        output_filepath,
        feature_dim,
        feature_timesteps,
        target_dim,
        target_timesteps,
        max_num_samples,
        chunk_dim = chunk_dim,
        init_file = false,
        attrs = attrs
    )

    seeds = Vector{Int}() # seeds are replaced by parallel collector
    scene = Scene(max_num_veh)
    models = Dict{Int, DriverModel}()
    col = DatasetCollector(seeds, gen, eval, dataset, scene, models, roadway, id = col_id)
    return col
end

function get_filepaths(filepath, n)
    dir = dirname(filepath)
    filename = basename(filepath)
    return [string(dir, "/proc_$(i)_$(filename)") for i in 1:n]
end

function build_parallel_dataset_collector(
        num_scenarios,
        num_col,
        output_filepath
    )
    filepaths = get_filepaths(output_filepath, num_col)
    num_scenarios_each = Int(ceil(num_scenarios / num_col))
    cols = [build_dataset_collector(filepaths[i], i, num_scenarios_each) for i in 1:num_col]
    seeds = collect(1:(num_scenarios))
    pcol = ParallelDatasetCollector(cols, seeds, output_filepath)
    return pcol
end


text:
using PyCall
using DifferentialEquations
const cn = PyNULL()
const signal = PyNULL()
run(`$(PyCall.python) -m pip install colorednoise`)
copy!(cn, pyimport("colorednoise"))
copy!(signal, pyimport_conda("scipy.signal", "scipy"))

function seed(theSeed=nothing) # Seed the rng, but return the seed. If no, nothing or NaN argument, randomly seed rng
    if isnothing(theSeed)
        theSeed = abs(Random.rand(Int64))
    end
    Random.seed!(theSeed)
    return theSeed
end

"""Define a function that, if it gets a Discontinuity, fills in tstops"""
function dsolve(prob, alg; kwargs...)
    if prob.p isa Discontinuous
        DifferentialEquations.solve(prob, alg; kwargs..., tstops=sort(collect(prob.p.d))) # May need to check tstops isn't in args in the future
    else
        DifferentialEquations.solve(prob, alg; kwargs...)
    end
end

"""
Whip up an ODEproblem from a Process
"""
function process2problem(P::Process, jac=nothing)
    args = [P.process, P.X0, (P.transient_t0, P.tmax), tuplef2ftuple(P.parameter_profile, P.parameter_profile_parameters)]
    isnothing(jac) || append!(ags, jac)
    if getalg(P) isa FunctionMap
        DiscreteProblem(args...)
    else
        ODEProblem(args...)
    end
end
export process2problem

process2solution(P::Process, jac=nothing) = dsolve(process2problem(P, jac), P.alg; dt = P.dt, saveat=P.savedt, P.solver_opts...)
export process2solution

function tuplef2ftuple(f, params)
    # turn a tuple of functions into a function of tuples
    if all(isempty.(params)) # The f's are just functions on their own, no need to add parameters
        # Be warned that you can't mix these; either use all parameter functions, or all standard functions. Don't be greedy.
        if f isa Tuple
            ds = [fi isa Discontinuous ? fi.d : [] for fi in f]
            ds = reduce(∪, ds) |> Set
            pp = Discontinuous(t->[x(t) for x in f], ds)
        else
            pp = f
        end
        return pp
    elseif f isa Tuple
        ps = Vector{Function}(undef, length(f))
        for i = 1:length(f)
            ps[i] = f[i](params[i]...)
        end
        ds = [fi isa Discontinuous ? fi.d : [] for fi in ps]
        ds = reduce(∪, ds) |> Set
        p = Discontinuous(t->map((x, g) -> g(x), fill(t, length(ps)), ps), ds) # Something like that
    else
        p = f(params...)
    end
    return p
end
export tuplef2ftuple

include("Processes/ARMA.jl")
include("Processes/ChaoticFlows.jl")
include("Processes/ChaoticMaps.jl")
include("Processes/DeterministicFlows.jl")
include("Processes/Noise.jl")
include("Processes/Signals.jl")
include("Processes/Stochastic.jl")
include("Processes/Transforms.jl")


text:
@testset "K-Fold cross-validation" begin
    target_R2 = [0.46405795117804827, 0.4314854720368352, 0.4471605188073827, 0.4650281295531955, 0.5353418610442862]
    target_ADJR2 = [0.45098619388970795, 0.41761926403773364, 0.433839085525633, 0.45213724110869413, 0.5241452793827027]
    target_TRAIN_MSE = [1.5167244045363427e6, 1.728038497400901e6, 1.766104797959807e6, 1.824818743281271e6, 1.6143508602587546e6]
    target_TRAIN_RMSE = [1231.5536547533536, 1314.548780913398, 1328.9487567095305, 1350.8585208234322, 1270.5710764293176]
    target_TEST_MSE = [2.430733983431065e6, 1.619864225581008e6, 1.447032419711302e6, 1.1981745824179968e6, 2.5641127272976018e6]
    target_TEST_RMSE = [1559.0811343323558, 1272.7388677890717, 1202.926606119967, 1094.61161259051, 1601.28471150436]

    year = [2011, 2011, 2011, 2011, 2012, 2010, 2011, 2010, 2011, 2010, 2010, 2011, 2011, 2010, 2011, 2011, 2010, 2010, 2011, 2011, 2010, 2010, 2011, 2010, 2009, 2010, 2010, 2010, 2010, 2009, 2010, 2009, 2011, 2011, 2009, 2010, 2011, 2010, 2010, 2010, 2010, 2010, 2010, 2009, 2010, 2010, 2010, 2010, 2009, 2009, 2010, 2010, 2010, 2009, 2010, 2010, 2010, 2009, 2009, 2009, 2010, 2010, 2009, 2009, 2009, 2010, 2010, 2010, 2009, 2009, 2009, 2009, 2010, 2009, 2010, 2009, 2009, 2010, 2010, 2010, 2009, 2009, 2007, 2010, 2010, 2010, 2009, 2009, 2009, 2008, 2009, 2009, 2009, 2010, 2010, 2008, 2008, 2009, 2009, 2008, 2009, 2010, 2008, 2009, 2010, 2008, 2008]
    mileage = [7413, 10926, 7351, 11613, 8367, 25125, 27393, 21026, 32655, 36116, 40539, 9199, 9388, 32058, 15367, 16368, 19926, 36049, 11662, 32069, 16035, 39943, 36685, 24920, 20019, 29338, 7784, 35636, 22029, 33107, 36306, 34419, 4867, 18948, 24030, 33036, 23967, 37905, 28955, 11165, 44813, 36469, 22143, 34046, 32703, 35894, 38275, 24855, 29501, 35394, 36447, 35318, 24929, 23785, 15167, 13541, 20278, 46126, 53733, 21108, 21721, 26716, 26887, 36252, 9450, 31414, 37185, 48174, 50533, 36713, 34888, 38380, 35574, 27528, 33302, 43369, 64055, 41342, 34503, 16573, 32403, 34846, 39665, 21325, 32743, 40058, 42325, 44518, 53902, 127327, 27136, 45813, 31538, 29517, 35871, 49787, 36323, 39211, 44789, 45996, 54988, 49720, 36322, 39222, 44089, 45993, 50988]
    price = [21992, 20995, 19995, 17809, 17500, 17495, 17000, 16995, 16995, 16995, 16995, 16992, 16950, 16950, 16000, 15999, 15999, 15995, 15992, 15992, 15988, 15980, 15899, 15889, 15688, 15500, 15499, 15499, 15298, 14999, 14999, 14995, 14992, 14992, 14992, 14990, 14989, 14906, 14900, 14893, 14761, 14699, 14677, 14549, 14499, 14495, 14495, 14480, 14477, 14355, 14299, 14275, 14000, 13999, 13997, 13995, 13995, 13995, 13995, 13992, 13992, 13992, 13992, 13991, 13950, 13950, 13950, 13895, 13888, 13845, 13799, 13742, 13687, 13663, 13599, 13584, 13425, 13384, 13383, 13350, 12999, 12998, 12997, 12995, 12995, 12995, 12995, 12995, 12995, 12995, 12992, 12990, 12988, 12849, 12780, 12777, 12704, 12595, 12507, 12500, 12500, 12773, 12704, 12595, 12519, 12333, 12399]
    df = DataFrame(Year = year, Mileage = mileage , Price = price)

    res = kfold(@formula(Price ~ 1 + Year + Mileage ), df, 5, 1, false)
    @test isapprox(target_R2, res.R2)
    @test isapprox(target_ADJR2, res.ADJR2)
    @test isapprox(target_TRAIN_MSE, res.TRAIN_MSE)
    @test isapprox(target_TRAIN_RMSE, res.TRAIN_RMSE)
    @test isapprox(target_TEST_MSE, res.TEST_MSE)
    @test isapprox(target_TEST_RMSE, res.TEST_RMSE)

end