text:
"""
load xplor or charmm (namd) format dcd file
"""
function readdcd(filename::String; index=nothing, stride=1)
    #TODO: endian
    header_ischarmm_4dims = false
    header_ischarmm_extrablock = false
    x = y = z = Matrix{Float64}(undef, 0, 0)
    boxsize = Matrix{Float64}(undef, 0, 0)

    open(filename, "r") do io
        seekend(io)
        file_size = position(io) # get file size in byte
        seekstart(io)
        blocksize = read(io, Int32) # should be 84
        # header (4 chars) either "CORD" or "VELD"; position(io) == 4 bytes
        header_hdr = Array{Char, 1}(undef, 4)
        for i in 1:4
            header_hdr[i] = read(io, Char)
        end
        # the total # of frames (snapshots); position(io) == 8 bytes
        header_nset = read(io, Int32)
        # starting time-step; position(io) == 12 bytes
        header_istrt = read(io, Int32)
        # frequency to save trajectory; position(io) == 16 bytes
        header_nsavc = read(io, Int32)
        # the total # of simulation step; position(io) == 20 bytes
        header_nstep = read(io, Int32)
        # null4 (int*4); position(io) == 24 bytes
        header_null4 = Array{Int32, 1}(undef, 4)
        read!(io, header_null4)
        # of degrees of freedom; position(io) == 40 bytes
        header_nfreat = read(io, Int32)
        # step-size of simulation; position(io) == 44 bytes
        header_delta = read(io, Float32)
        # null9 (int*9); position(io) == 48 bytes
        header_null9 = Array{Int32, 1}(undef, 9)
        read!(io, header_null9)
        # version; position(io) == 84 bytes
        header_version = read(io, Int32)
        # charmm extension format
        if header_version > 0
            header_ischarmm = true
        end
        # delta is double precision in xplor format
        if header_ischarmm == false
            cof = position(io)
            # reread delta in double precision
            fseek(fid, 44)
            header_delta = read(io, Float64)
            seek(io, cof)
        end
        # check charmm extensions
        if header_ischarmm == true
            cof = position(io)
            # charmm extrablock extension
            seek(io, 48)
            n = read(io, Int32)
            if n == 1
                header_ischarmm_extrablock = true
            end
            # charmm 4dims extension
            seek(io, 52);
            n = read(io, Int32)
            if n == 1
                header_ischarmm_4dims = true
            end
            seek(io, cof)
        end
        # blocksize1; position(io) == 88 bytes
        blocksize1 = read(io, Int32)
        #### read block 2 (title)
        # blocksize2; position(io) == 92 bytes
        blocksize2 = read(io, Int32)
        # of title lines; position(io) == 96 bytes
        header_ntitle = read(io, Int32)
        # title
        header_title = []
        for i in 1:header_ntitle
            t = Array{Char, 1}(undef, 80)
            for i in 1:80
                t[i] = read(io, Char)
            end
            push!(header_title, String(t))
        end
        # blocksize2
        blocksize2 = read(io, Int32)
        #### read block 3 (natom)
        # blocksize3
        blocksize3 = read(io, Int32)
        # # of atoms
        header_natom = read(io, Int32)
        # blocksize3
        blocksize3 = read(io, Int32)
        ## read coordinates
        header_size = position(io)
        if header_ischarmm_extrablock == true
            extrablocksize = 4*2 + 8*6
        else
            extrablocksize = 0
        end
        coordblocksize = (4*2 + 4*header_natom)*3
        nframe = (file_size - header_size) / (extrablocksize + coordblocksize)
        nframe = Int64(floor(nframe))
        nframe_stride = length(1:stride:nframe)
        if typeof(index) <: AbstractVector{Bool}
            index2 = findall(index)
        elseif index == nothing
            index2 = collect(1:header_natom)
        else
            index2 = index
        end
        x = zeros(Float64, (nframe_stride, length(index2)))
        y = zeros(Float64, (nframe_stride, length(index2)))
        z = zeros(Float64, (nframe_stride, length(index2)))
        if header_ischarmm_extrablock == true
            boxsize = zeros(Float64, (nframe_stride, 3))
        end
        # read next frames
        dummy = Array{Float64, 1}(undef, 6)
        crd_x = Array{Float32, 1}(undef, header_natom)
        crd_y = Array{Float32, 1}(undef, header_natom)
        crd_z = Array{Float32, 1}(undef, header_natom)
        for iframe in 1:nframe_stride
            # read charmm extrablock (unitcell info)
            if header_ischarmm_extrablock == true
                blocksize_box = read(io, Int32)
                read!(io, dummy)
                blocksize_box = read(io, Int32)
            end
            # read x coordinates
            blocksize_x = read(io, Int32)
            read!(io, crd_x)
            blocksize_x = read(io, Int32)
            # read y coordinates
            blocksize_y = read(io, Int32)
            read!(io, crd_y)
            blocksize_y = read(io, Int32)
            # read z coordinates
            blocksize_z = read(io, Int32)
            read!(io, crd_z)
            blocksize_z = read(io, Int32)
            # skip charmm 4dims extension
            if header_ischarmm_4dims == true
                blocksize_4dims = read(io, Int32)
                skip(io, blocksize_4dims)
                blocksize_4dims = read(io, Int32)
            end
            # get data
            x[iframe, :] = crd_x[index2];
            y[iframe, :] = crd_y[index2];
            z[iframe, :] = crd_z[index2];
            if header_ischarmm_extrablock == true
                boxsize[iframe, :] = dummy[[1 3 6]];
            end

            # skip stride-1 frames
            if iframe < nframe_stride
                for ii = 1:(stride-1)
                    if header_ischarmm_extrablock == true
                        blocksize_box = read(io, Int32)
                        skip(io, blocksize_box)
                        blocksize_box = read(io, Int32)
                    end
                    # read x coordinates
                    blocksize_x = read(io, Int32)
                    skip(io, blocksize_x)
                    blocksize_x = read(io, Int32)
                    # read y coordinates
                    blocksize_y = read(io, Int32)
                    skip(io, blocksize_y)
                    blocksize_y = read(io, Int32)
                    # read z coordinates
                    blocksize_z = read(io, Int32)
                    skip(io, blocksize_z)
                    blocksize_z = read(io, Int32)
                    # skip charmm 4dims extension
                    if header_ischarmm_4dims == true
                        blocksize_4dims = read(io, Int32)
                        skip(io, blocksize_4dims)
                        blocksize_4dims = read(io, Int32)
                    end
                end
            end
        end
    end

    if isempty(boxsize)
        TrjArray{Float64, Int64}(x=x, y=y, z=z)
    else
        TrjArray{Float64, Int64}(x=x, y=y, z=z, boxsize=boxsize)
    end
end

############################################################################
"""
read netcdf file
"""
function readnetcdf(filename::String; index=nothing)
    #finfo = ncinfo(filename)
    #attributes = finfo.gatts
    #dimensions = finfo.dims
    nc = NetCDF.open(filename, mode=NetCDF.NC_64BIT_OFFSET);
    nframe = Int64(nc.dim["frame"].dimlen)
    natom = Int64(nc.dim["atom"].dimlen)

    is_trj = haskey(nc.vars, "coordinates") ? true : false
    is_box = haskey(nc.vars, "cell_lengths") ? true : false
    is_vel = haskey(nc.vars, "velocities") ? true : false
    is_temp = haskey(nc.vars, "temp0") ? true : false

    NetCDF.close(nc)

    if typeof(index) <: AbstractVector{Bool}
        index2 = findall(index)
    elseif index == nothing
        index2 = collect(1:natom)
    else
        index2 = index
    end

    start_atom = minimum(index2)
    count_atom = maximum(index2) - start_atom + 1
    index3 = index2 .- start_atom .+ 1
    start_time = 1
    count_time = nframe

    if is_trj
        d = ncread(filename, "coordinates", start=[1, start_atom, start_time], count=[3, count_atom, count_time])
        d = map(Float64, d)
        x = d[1, index3, :]'
        y = d[2, index3, :]'
        z = d[3, index3, :]'
    else
        x = y = z = Matrix{Float64}(undef, 0, 0)
    end

    if is_box
        d = ncread(filename, "cell_lengths", start=[1, start_time], count=[3, count_time]);
        d = map(Float64, d)
        boxsize = d'
    else
        boxsize = Matrix{Float64}(undef, 0, 0)
    end

    if is_vel
        d = ncread(filename, "velocities", start=[1, start_atom, start_time], count=[3, count_atom, count_time])
        d = map(Float64, d)
        vx = d[1, index3, :]'
        vy = d[2, index3, :]'
        vz = d[3, index3, :]'
    else
        vx = vy = vz = Matrix{Float64}(undef, 0, 0)
    end

    if is_temp
        d = ncread(filename, "temp0", start=[start_time], count=[count_time]);
        temp = map(Float64, d)
    else
        temp = Vector{Float64}(undef, 0)
    end

    TrjArray{Float64, Int64}(x=x, y=y, z=z, boxsize=boxsize)
end

############################################################################
"""
write netcdf file
"""
function writenetcdf(filename::String, ta::TrjArray; velocity = nothing, force = nothing)
    scale_factor = 20.455
    natom = ta.natom
    nframe = ta.nframe

    # global attributes
    gatts = Dict(
        "ConventionVersion" => "1.0",
        "Conventions" => "AMBER",
        "application" => "Julia",
        "program" => "MDToolbox.jl",
        "programVersion" => "0.1",
        "title" => "Created by MDToolbox.jl"
    )

    # dimensions
    ncdim_frame = NcDim("frame", nframe, unlimited=true)
    ncdim_atom = NcDim("atom", natom)
    ncdim_spatial = NcDim("spatial", 3)
    ncdim_cell_spatial = NcDim("cell_spatial", 3)
    ncdim_cell_angular = NcDim("cell_angular", 3)
    ncdim_label = NcDim("label", 5)

    # variables
    ncvar_spatial = NcVar("spatial", ncdim_spatial, t=NetCDF.NC_CHAR)
    ncvar_time = NcVar("time", ncdim_frame, t=NetCDF.NC_FLOAT, atts=Dict("units" => "picosecond"))
    ncvar_coordinates = NcVar("coordinates", [ncdim_spatial, ncdim_atom, ncdim_frame], t=NetCDF.NC_FLOAT, atts=Dict("units" => "angstrom"))
    if !isempty(ta.boxsize)
        ncvar_cell_spatial = NcVar("cell_spatial", ncdim_cell_spatial, t=NetCDF.NC_CHAR)
        ncvar_cell_angular = NcVar("cell_angular", [ncdim_label, ncdim_cell_angular], t=NetCDF.NC_CHAR)
        ncvar_cell_lengths = NcVar("cell_lengths", [ncdim_cell_spatial, ncdim_frame], t=NetCDF.NC_FLOAT, atts=Dict("units" => "angstrom"))
        ncvar_cell_angles = NcVar("cell_angles", [ncdim_cell_angular, ncdim_frame], t=NetCDF.NC_FLOAT, atts=Dict("units" => "degree"))
    end
    if velocity != nothing
        ncvar_velocities = NcVar("velocities", [ncdim_spatial, ncdim_atom, ncdim_frame], t=NetCDF.NC_FLOAT, atts=Dict("units" => "angstrom/picosecond"))
    end
    if force != nothing
        ncvar_forces = NcVar("forces", [ncdim_spatial, ncdim_atom, ncdim_frame], t=NetCDF.NC_FLOAT, atts=Dict("units" => "amu*angstrom/picosecond^2"))
    end
    #ncvar_temp0 = NcVar("temp0", ncdim_frame, t=NetCDF.NC_FLOAT, atts=Dict("units" => "kelvin"))

    # create the NetCDF file
    varlist = [ncvar_spatial, ncvar_time, ncvar_coordinates]
    if !isempty(ta.boxsize)
        push!(varlist, ncvar_cell_spatial)
        push!(varlist, ncvar_cell_angular)
        push!(varlist, ncvar_cell_lengths)
        push!(varlist, ncvar_cell_angles)
    end
    if velocity != nothing
        push!(varlist, ncvar_velocities)
    end
    if force != nothing
        push!(varlist, ncvar_forces)
    end
    nc = NetCDF.create(filename, varlist, gatts=gatts, mode=NetCDF.NC_64BIT_OFFSET)

    # write data
    # check by ncread("tmp.nc", "spatial")
    c = zeros(UInt8, 3); copyto!(c, 1, "xyz", 1); #@show c
    NetCDF.putvar(nc, "spatial", c)
    NetCDF.putvar(nc, "time", zeros(Float32, nframe))
    trj = zeros(Float32, 3, natom, nframe)
    for iframe in 1:nframe
        for iatom in 1:natom
            trj[1, iatom, iframe] = ta.x[iframe, iatom]
            trj[2, iatom, iframe] = ta.y[iframe, iatom]
            trj[3, iatom, iframe] = ta.z[iframe, iatom]
        end
    end
    NetCDF.putvar(nc, "coordinates", trj)

    if !isempty(ta.boxsize)
        c = zeros(UInt8, 3); copyto!(c, 1, "abc", 1); #@show c
        NetCDF.putvar(nc, "cell_spatial", c)
        c = zeros(UInt8, 5, 3); copyto!(c, 1, "alpha", 1); copyto!(c, 6, "beta ", 1); copyto!(c, 11, "gamma", 1); #@show c
        NetCDF.putvar(nc, "cell_angular", c)
        d = zeros(Float32, 3, nframe)
        d[:, :] .= ta.boxsize'
        NetCDF.putvar(nc, "cell_lengths", d)
        d[:, :] .= 90.0
        NetCDF.putvar(nc, "cell_angles", d)
    end

    if velocity != nothing
        for iframe in 1:nframe
            for iatom in 1:natom
                trj[1, iatom, iframe] = velocity.x[iframe, iatom]/scale_factor
                trj[2, iatom, iframe] = velocity.y[iframe, iatom]/scale_factor
                trj[3, iatom, iframe] = velocity.z[iframe, iatom]/scale_factor
            end
        end
        NetCDF.putvar(nc, "velocities", trj)
    end
    if force != nothing
        for iframe in 1:nframe
            for iatom in 1:natom
                trj[1, iatom, iframe] = force.x[iframe, iatom]; #TODO: scale_factor?
                trj[2, iatom, iframe] = force.y[iframe, iatom]
                trj[3, iatom, iframe] = force.z[iframe, iatom]
            end
        end
        NetCDF.putvar(nc, "forces", trj)
    end

    NetCDF.close(nc)
end

############################################################################
function parse_line(line::String, index, mytype::DataType, default_value)
    try
        if mytype == String
            return strip(string(line[index]))
        else
            return parse(mytype, line[index])
        end
    catch
        return default_value
    end
end

############################################################################
"""
read charmm or xplor type psf file
"""
function readpsf(filename::String)
    isPSF = false
    isEXT = false
    isCMAP = false
    isCHEQ = false

    lines = open(filename, "r" ) do fp
        readlines(fp)
    end

    # first line
    line = lines[1]
    line_size = length(line)
    if line_size >= 3
        if line[1:3] == "PSF"
            isPSF = true
        end
    end
    for i = 1:4:line_size
        i_max = minimum([i+3, line_size])
        isPSF = strip(line[i:i_max]) == "PSF" ? true : isPSF
        isEXT = strip(line[i:i_max]) == "EXT" ? true : isEXT
        isCMAP = strip(line[i:i_max]) == "CMAP" ? true : isCMAP
        isCHEQ = strip(line[i:i_max]) == "CHEQ" ? true : isCHEQ
    end
    if !isPSF
        print("Sorry, this seems not be a PSF file")
        return false
    end
    if isEXT
        #fmt_atom = "%10d| %8s| %8d| %8s| %8s| %4s| %14f|%14f|%8d %*f %*f" %octave
        fmt_atom = [1:10, 12:19, 21:28, 30:37, 39:46, 48:51, 53:66, 67:80, 81:88]
        #fmt_list = "%10d"
    else
        #fmt_atom = "%8d %4s %4d %4s| %4s| %4s| %14f%14f%8d %*f %*f" #octave
        fmt_atom = [1:8, 10:13, 15:18, 20:23, 25:28, 30:33, 34:47, 48:61, 62:69]
        #fmt_list = "%8d"
        fmt_list = "%8d"
    end

    psf_title = Vector{String}(undef, 0)
    psf_segment_name = Vector{String}(undef, 0)
    psf_residue_name = Vector{String}(undef, 0)
    psf_residue_id = Vector{Int64}(undef, 0)
    psf_atom_name = Vector{String}(undef, 0)
    psf_atom_type_letter = Vector{String}(undef, 0)
    psf_atom_type_number = Vector{Int64}(undef, 0)
    psf_atom_id = Vector{Int64}(undef, 0)
    psf_charge = Vector{Float64}(undef, 0)
    psf_mass = Vector{Float64}(undef, 0)
    list_bond = Array{Int64}(undef, 0, 0)
    list_angle = Array{Int64}(undef, 0, 0)
    list_dihedral = Array{Int64}(undef, 0, 0)
    list_improper = Array{Int64}(undef, 0, 0)
    list_cmap = Array{Int64}(undef, 0, 0)

    iline = 2
    while iline <= length(lines)
        line = lines[iline]; iline += 1

        if occursin(r".*\d.*!\w.*", line) && !occursin("NGRP", line)
            line_splitted = split(line, "!")
            #num = parse(Int64, line_splitted[1])
            num = parse_line(String(line_splitted[1]), 1:length(line_splitted[1]), Int64, 0)
            key = strip(line_splitted[2])
            if occursin(r"^NATOM", key)
                for ii = iline:(iline+num-1)
                    line = lines[ii]
                    push!(psf_atom_id, parse_line(line, fmt_atom[1], Int64, 0))
                    push!(psf_segment_name, parse_line(line, fmt_atom[2], String, "None"))
                    push!(psf_residue_id, parse_line(line, fmt_atom[3], Int64, 0))
                    push!(psf_residue_name, parse_line(line, fmt_atom[4], String, "None"))
                    push!(psf_atom_name, parse_line(line, fmt_atom[5], String, "None"))
                    if isletter(strip(line[fmt_atom[6]])[1])
                        push!(psf_atom_type_letter, parse_line(line, fmt_atom[6], String, "None")) # CHARMM format
                    else
                        push!(psf_atom_type_number, parse_line(line, fmt_atom[6], Int64, 0)) # XPLOR format
                    end
                    push!(psf_charge, parse_line(line, fmt_atom[7], Float64, 0.0))
                    push!(psf_mass, parse_line(line, fmt_atom[8], Float64, 0.0))
                end
                iline += num

            elseif occursin(r"^NBOND", key)
                a = Vector{Int64}(undef, 0)
                while length(a) < (num*2)
                    line = lines[iline]; iline += 1
                    append!(a, [parse(Int64, y) for y in split(line)])
                end
                list_bond = collect(reshape(a, 2, :)')

            elseif occursin(r"^NTHETA", key)
                a = Vector{Int64}(undef, 0)
                while length(a) < (num*3)
                    line = lines[iline]; iline += 1
                    append!(a, [parse(Int64, y) for y in split(line)])
                end
                list_angle = collect(reshape(a, 3, :)')

            elseif occursin(r"^NPHI", key)
                a = Vector{Int64}(undef, 0)
                while length(a) < (num*4)
                    line = lines[iline]; iline += 1
                    append!(a, [parse(Int64, y) for y in split(line)])
                end
                list_dihedral = collect(reshape(a, 4, :)')

            elseif occursin(r"^NIMPHI", key)
                a = Vector{Int64}(undef, 0)
                while length(a) < (num*4)
                    line = lines[iline]; iline += 1
                    append!(a, [parse(Int64, y) for y in split(line)])
                end
                list_improper = collect(reshape(a, 4, :)')

            elseif occursin(r"^NCRTERM", key)
                a = Vector{Int64}(undef, 0)
                while length(a) < (num*8)
                    line = lines[iline]; iline += 1
                    append!(a, [parse(Int64, y) for y in split(line)])
                end
                list_cmap = collect(reshape(a, 8, :)')

            end
        end
    end

    TrjArray{Float64, Int64}(chainname=psf_segment_name,
             resname=psf_residue_name, resid=psf_residue_id,
             atomname=psf_atom_name, atomid=psf_atom_id,
             mass=psf_mass, charge=psf_charge,
             list_bond=list_bond, list_angle=list_angle,
             list_dihedral=list_dihedral, list_improper=list_improper,
             list_cmap=list_cmap)
end

############################################################################
"""
write psf file
"""
function writepsf(filename::String, ta::TrjArray)
    natom = ta.natom

    open(filename, "w") do io
        Printf.@printf(io, "%4s", "PSF ")
        if !isempty(ta.list_cmap)
            Printf.@printf(io, "%4s", "CMAP")
        end
        Printf.@printf(io, "\n\n")

        Printf.@printf(io, "%8d !NTITLE\n", 1)
        Printf.@printf(io, "CREATED by MDToolbox.jl\n")
        Printf.@printf(io, "\n")

        Printf.@printf(io, "%8d !NATOM\n", natom)
        for iatom = 1:natom
            Printf.@printf(io, "%8d", iatom)
            Printf.@printf(io, " %4s", isempty(ta.chainname) ? rpad("NONE", 4) : rpad(ta.chainname[iatom], 4))
            Printf.@printf(io, " %-4d", isempty(ta.resid) ? 0 : ta.resid[iatom])
            Printf.@printf(io, " %4s", isempty(ta.resname) ? rpad("NONE", 4) : rpad(ta.resname[iatom], 4))
            Printf.@printf(io, " %4s", isempty(ta.atomname) ? rpad("NONE", 4) : rpad(ta.atomname[iatom], 4))
            Printf.@printf(io, " %4s", isempty(ta.atomname) ? rpad("NONE", 4) : rpad(ta.atomname[iatom], 4))
            Printf.@printf(io, "%14.9f", isempty(ta.charge) ? 0.0 : ta.charge[iatom])
            Printf.@printf(io, "%14.7f", isempty(ta.mass) ? 0.0 : ta.mass[iatom])
            Printf.@printf(io, "%9d", 0)
            Printf.@printf(io, "\n")
        end

        if !isempty(ta.list_bond)
            Printf.@printf(io, "%8d !NBOND\n", size(ta.list_bond, 1))
            for i = 1:size(ta.list_bond, 1)
                for j = 1:size(ta.list_bond, 2)
                    Printf.@printf(io, "%8d", ta.list_bond[i, j])
                end
                if (i % 4) == 0
                    Printf.@printf(io, "\n")
                end
            end
            Printf.@printf(io, "\n")
        end

        if !isempty(ta.list_angle)
            Printf.@printf(io, "%8d !NTHETA\n", size(ta.list_angle, 1))
            for i = 1:size(ta.list_angle, 1)
                for j = 1:size(ta.list_angle, 2)
                    Printf.@printf(io, "%8d", ta.list_angle[i, j])
                end
                if (i % 3) == 0
                    Printf.@printf(io, "\n")
                end
            end
            Printf.@printf(io, "\n")
        end

        if !isempty(ta.list_dihedral)
            Printf.@printf(io, "%8d !NPHI\n", size(ta.list_dihedral, 1))
            for i = 1:size(ta.list_dihedral, 1)
                for j = 1:size(ta.list_dihedral, 2)
                    Printf.@printf(io, "%8d", ta.list_dihedral[i, j])
                end
                if (i % 2) == 0
                    Printf.@printf(io, "\n")
                end
            end
            Printf.@printf(io, "\n")
        end

        if !isempty(ta.list_improper)
            Printf.@printf(io, "%8d !NIMPHI\n", size(ta.list_improper, 1))
            for i = 1:size(ta.list_improper, 1)
                for j = 1:size(ta.list_improper, 2)
                    Printf.@printf(io, "%8d", ta.list_improper[i, j])
                end
                if (i % 2) == 0
                    Printf.@printf(io, "\n")
                end
            end
            Printf.@printf(io, "\n")
        end

        if !isempty(ta.list_cmap)
            Printf.@printf(io, "%8d !NCRTERM\n", size(ta.list_cmap, 1))
            for i = 1:size(ta.list_cmap, 1)
                for j = 1:size(ta.list_cmap, 2)
                    Printf.@printf(io, "%8d", ta.list_cmap[i, j])
                end
                Printf.@printf(io, "\n")
            end
            Printf.@printf(io, "\n")
        end
    end
end

############################################################################
"""
read protein data bank (PDB) file
"""
function readpdb(filename::String)
    lines = open(filename, "r" ) do fp
        readlines(fp)
    end

    model = []
    lines_clean = []
    for i = 1:length(lines)
        line = lines[i]
        if match(r"^ATOM", line) != nothing
            push!(lines_clean, line)
        end
        if match(r"^HETATM", line) != nothing
            push!(lines_clean, line)
        end
        if match(r"^ENDMDL", line) != nothing || (i == length(lines) && !isempty(lines_clean))
            push!(model, lines_clean)
            lines_clean = []
        end
    end

    natom = length(model[1])
    pdb_record = Vector{String}(undef, natom)
    pdb_serial = Vector{Int64}(undef, natom)
    pdb_name = Vector{String}(undef, natom)
    pdb_altloc = Vector{String}(undef, natom)
    pdb_resname = Vector{String}(undef, natom)
    pdb_chainid = Vector{String}(undef, natom)
    pdb_resseq = Vector{Int64}(undef, natom)
    pdb_x = zeros(Float64, length(model), natom)
    pdb_y = zeros(Float64, length(model), natom)
    pdb_z = zeros(Float64, length(model), natom)
    pdb_occupancy = Vector{Float64}(undef, natom)
    pdb_tempfactor = Vector{Float64}(undef, natom)
    pdb_element = Vector{String}(undef, natom)
    pdb_charge = Vector{Float64}(undef, natom)

    lines = model[1]
    for iatom = 1:natom
        line = lines[iatom]
        pdb_record[iatom] = parse_line(line, 1:6, String, "None")
        num = parse_line(line, 7:12, String, "0")
        if iatom == 1 || pdb_serial[iatom-1] < 99999
            pdb_serial[iatom] = parse(Int64, num)
        elseif match(r"\*", num) != nothing
            pdb_serial[iatom] = iatom == 1 ? 1 : pdb_serial[iatom-1]
        else
            pdb_serial[iatom] = parse(Int64, num)
        end
        pdb_name[iatom] = parse_line(line, 13:16, String, "None")
        pdb_altloc[iatom] = parse_line(line, 17:17, String, "None")
        pdb_resname[iatom] = parse_line(line, 18:21, String, "None")
        pdb_chainid[iatom] = parse_line(line, 22:22, String, "None")
        pdb_resseq[iatom] = parse_line(line, 23:28, Int64, 0)
        pdb_x[1, iatom] = parse_line(line, 31:38, Float64, 0.0)
        pdb_y[1, iatom] = parse_line(line, 39:46, Float64, 0.0)
        pdb_z[1, iatom] = parse_line(line, 47:54, Float64, 0.0)
        pdb_occupancy[iatom] = parse_line(line, 55:60, Float64, 0.0)
        pdb_tempfactor[iatom] = parse_line(line, 61:66, Float64, 0.0)
        pdb_element[iatom] = parse_line(line, 77:78, String, "None")
        pdb_charge[iatom] = parse_line(line, 79:80, Float64, 0.0)
    end

    for imodel = 2:length(model)
        lines = model[imodel]
        for iatom = 1:natom
            line = lines[iatom]
            pdb_x[imodel, iatom] = parse_line(line, 31:38, Float64, 0.0)
            pdb_y[imodel, iatom] = parse_line(line, 39:46, Float64, 0.0)
            pdb_z[imodel, iatom] = parse_line(line, 47:54, Float64, 0.0)
        end
    end

    TrjArray{Float64, Int64}(x=pdb_x, y=pdb_y, z=pdb_z,
             chainname=pdb_chainid,
             resid=pdb_resseq, resname=pdb_resname,
             atomid=pdb_serial, atomname=pdb_name)
end

############################################################################
"""
write pdb file
"""
function writepdb(io::IO, ta::TrjArray; format_type="vmd")
    natom = ta.natom
    nframe = ta.nframe

    if format_type == "namd"
        Printf.@printf io "CRYST1    0.000    0.000    0.000  90.00  90.00  90.00 P 1           1\n"
    end

    for iframe = 1:nframe
        for iatom = 1:natom
            Printf.@printf(io, "%6s", "ATOM  ")
            if isempty(ta.atomid)
                Printf.@printf(io, "%5d", mod(iatom, 100000))
            else
                Printf.@printf(io, "%5d", mod(ta.atomid[iatom], 100000))
            end
            Printf.@printf(io, "%1s", " ")
            if isempty(ta.atomname)
                Printf.@printf(io, "%4s", " CA ")
            else
                Printf.@printf(io, "%4s", rpad(ta.atomname[iatom], 3))
            end
            Printf.@printf(io, "%1s", " ")
            #Printf.@printf(io, "%3s", ta.resname(iatom, :))
            #Printf.@printf(io, "%1s", " ")
            if isempty(ta.resname)
                Printf.@printf(io, "%4s", "ALA ")
            else
                Printf.@printf(io, "%4s", rpad(ta.resname[iatom], 4))
            end
            if isempty(ta.chainname)
                Printf.@printf(io, "%1s", " ")
            else
                Printf.@printf(io, "%1s", ta.chainname[iatom])
            end
            if isempty(ta.chainname)
                Printf.@printf(io, "%4d", mod(iatom, 10000))
            else
                Printf.@printf(io, "%4d", mod(ta.resid[iatom], 10000))
            end
            #Printf.@printf(io, "%1s", ta.icode(iatom, :))
            Printf.@printf(io, "%1s", " ")
            Printf.@printf(io, "%1s", " ")
            Printf.@printf(io, "%1s", " ")
            Printf.@printf(io, "%1s", " ")
            #Printf.@printf(io, "%8.3f", ta.xyz(iatom, 1))
            #Printf.@printf(io, "%8.3f", ta.xyz(iatom, 2))
            #Printf.@printf(io, "%8.3f", ta.xyz(iatom, 3))
            Printf.@printf(io, "%8.3f", ta.x[iframe, iatom])
            Printf.@printf(io, "%8.3f", ta.y[iframe, iatom])
            Printf.@printf(io, "%8.3f", ta.z[iframe, iatom])
            Printf.@printf(io, "%6.2f", 0.0)
            Printf.@printf(io, "%6.2f", 0.0)
            Printf.@printf(io, "%1s", " ")
            Printf.@printf(io, "%1s", " ")
            Printf.@printf(io, "%1s", " ")
            Printf.@printf(io, "%1s", " ")
            Printf.@printf(io, "%1s", " ")
            Printf.@printf(io, "%1s", " ")
            Printf.@printf(io, "%1s", " ")
            Printf.@printf(io, "%1s", " ")
            Printf.@printf(io, "%1s", " ")
            Printf.@printf(io, "%1s", " ")
            Printf.@printf(io, "%2s", "  "); #element
            Printf.@printf(io, "%2s", "  "); #charge
            Printf.@printf(io, "\n")
        end
        if nframe > 1
            Printf.@printf(io, "TER\n");
            Printf.@printf(io, "ENDMDL\n");
        end
        #Printf.@printf " %8.2f" X[i,j]
    end
end

"""
write pdb file
"""
function writepdb(filename::String, ta::TrjArray; format_type="vmd")
    open(filename, "w") do io
        writepdb(io, ta, format_type=format_type)
    end
end

############################################################################
"""
read amber coordinates file
"""
function readcrd(filename::String)
    lines = open(filename, "r" ) do fp
        readlines(fp)
    end

    title = lines[1]
    natom = parse_line(lines[2], 1:length(lines[2]), Int64, 0)
    #pdb_name[iatom] = parse_line(line, 13:16, String, "None")

    xyz = Vector{Float64}(undef, natom*3)

    iatom3 = 1
    icount = 3
    while true
        n = length(lines[icount])
        for i = 1:12:n
            xyz[iatom3] = parse_line(lines[icount], i:(i+11), Float64, 0.0)
            iatom3 += 1
        end
        icount += 1
        if iatom3 > (natom*3)
            break
        end
    end

    x = zeros(Float64, 1, natom)
    y = zeros(Float64, 1, natom)
    z = zeros(Float64, 1, natom)
    x[:] .= xyz[1:3:end]
    y[:] .= xyz[2:3:end]
    z[:] .= xyz[3:3:end]

    TrjArray{Float64, Int64}(x=x, y=y, z=z)
end


text:
using Test
using CthulhuVision.Random
using CthulhuVision.Math

using CUDAnative, CuArrays

@inline function randomunitvectors_gpu(result, n)
    index = threadIdx().x + (blockIdx().x - 1) * blockDim().x

    seedgen = SplitMix64(index)
    s0 = next(seedgen)
    s1 = next(seedgen)
    s2 = next(seedgen)
    s3 = next(seedgen)
    rng = UniformRNG(Xoshiro256pp(s0, s1, s2, s3))

    if index <= n
        @inbounds result[index] = randominunitsphere(rng)
    end

    nothing
end

@testset "Math GPU                " begin
    @testset "Random vectors in unit sphere have length < 1.0f0" begin
        threads = 32
        blocks = 8
        n = threads * blocks
        result = CuArray{Vec3}(undef, n)

        CUDAnative.@sync begin
            @cuda threads=threads blocks=blocks randomunitvectors_gpu(result, n)
        end

        result_host = Vector{Vec3}(result)

        for i = 1:n
            v = result_host[i]
            l = sqrt(v.x * v.x + v.y * v.y + v.z * v.z)
            @test l < 1.0f0
        end
    end
end

struct TTV
    v::Vec3
    t::Transform
    expected::Vec3
end

@testset "Math Transform          " begin
    vec3_transforms = [
        TTV(
            Vec3(0f0, 0f0, 0f0),
            translation(1f0, 0f0, 0f0),
            Vec3(1f0, 0f0, 0f0),
        ),

        TTV(
            Vec3(0f0, 0f0, 0f0),
            translation(0f0, 1f0, 0f0),
            Vec3(0f0, 1f0, 0f0),
        ),

        TTV(
            Vec3(0f0, 0f0, 0f0),
            translation(0f0, 0f0, 1f0),
            Vec3(0f0, 0f0, 1f0),
        ),

        TTV(
            Vec3(0f0, 0f0, 0f0),
            translation(1f0, 2f0, 3f0),
            Vec3(1f0, 2f0, 3f0),
        ),

        TTV(
            Vec3(4f0, 5f0, 6f0),
            translation(1f0, 2f0, 3f0),
            Vec3(5f0, 7f0, 9f0),
        ),

        TTV(
            Vec3(2f0, 0f0, 0f0),
            translation(1f0, 0f0, 0f0),
            Vec3(3f0, 0f0, 0f0),
        ),
    ]

    @testset "Transform * Vector" begin
        for ttv in vec3_transforms
            actual = ttv.t * ttv.v

            @test actual == ttv.expected
        end
    end
end

ext: jl


text:
const ALLOCATED = []

struct CUDAMalloc
end

function (::CUDAMalloc)(::Type{T}, dims::Dims{N}) where {T,N}
    bytesize = sizeof(T) * prod(dims)
    @assert bytesize > 0
    ptr = Ptr{T}(memalloc(bytesize))
    arr = CuArray(ptr, dims, getdevice())
    push!(ALLOCATED, arr)
    finalizer(x -> memfree(x.ptr), arr)
    arr
end

function sss()
    ref = Ref{Ptr{Cvoid}}()
    @apicall :cuMemAlloc (Ptr{Ptr{Cvoid}},Csize_t) ref n*sizeof(T)
    dev = getdevice()
    mb = MemBlock(Ptr{T}(ref[]), n, dev)
    finalizer(mb, memfree)
end

function rrr()
    ref = Ref{Ptr{Cvoid}}()
    status = @unsafe_apicall :cuMemAlloc (Ptr{Ptr{Cvoid}},Csize_t) ref bytesize
    if status != CUDA_SUCCESS
        gc(false)
        status = @unsafe_apicall :cuMemAlloc (Ptr{Ptr{Cvoid}},Csize_t) ref bytesize
        if status != CUDA_SUCCESS
            gc()
            memalloc(bytesize)
        end
    end
    ref[]
end


text:
module TotalLeastSquares
export tls, tls!, wtls, wls, rtls, irls, sls, rowcovariance, hankel, ishankel, unhankel, flts
export rpca, lowrankfilter, rpca_ga, entrywise_median, entrywise_trimmed_mean, μ!
using FillArrays, Printf, LinearAlgebra, SparseArrays, Statistics, StatsBase

include("flts.jl")

"""
    wls(A,y,Σ)

Solves the weigted standard least-squares problem Ax = y + e, e ~ N(0,Σ)
# Arguments
- `A ∈ R(n,u)` Design matrix
- `y ∈ R(n)` RHS
- `Σ ∈ R(n,n)` Covariance matrix of the residuals (can be sparse or already factorized).
"""
function wls(A,y,Σ)
    Symmetric(A'*(Σ\A))\(A'*(Σ\y))
end

"""
    wls!(zA, A, y, w::Vector)

Overwrites `zA` *and* `w`.
- `zA`: storage of same size as `A`
"""
function wls!(zA,A,y,w)
    zA .= A ./ w
    w  .= y ./ w
    Symmetric(A'*zA)\(A'w)
end

wls(A,y,Σ::Union{Matrix, SparseMatrixCSC}) = wls(A,y,factorize(Hermitian(Σ)))

"""
    tls(A,y)

Solves the total least-squares problem Ax=y using the SVD method
# Arguments
- `A` Design matrix
- `y` RHS
"""
function tls(A::AbstractArray,y::AbstractArray)
    AA  = [A y]
    s   = svd(AA) # Not in-place for Zygote compatability
    n   = size(A,2)
    V21 = s.V[1:n,n+1:end]
    V22 = s.V[n+1:end,n+1:end]
    x   = -V21/V22
end

"""
    tls!(Ay::AbstractArray, n::Integer)
    tls!(s::SVD, n::Integer)

Inplace version of `tls`. `Ay` is `[A y]` and `n` is the number of columns in `A`. Also accepts a precomputed SVD of `Ay`.
"""
tls!(Ay::AbstractArray, n::Integer)  = tls!(svd!(Ay), n)

function tls!(s::SVD, n::Integer)
    V21 = s.V[1:n,n+1:end]
    V22 = s.V[n+1:end,n+1:end]
    x   = -V21/V22
end

a ⊗ b = kron(a,b)

"""
    x = wtls(A,y,Qaa,Qay,Qyy; iters = 10)

Solves min nᵀQ⁻¹n s.t. (A+E)x = y + v
where Q = [Qaa Qay; Qay' Qyy], n = [vec(E); v]

Uses algorithm 1 from
Weighted total least squares: necessary and sufficient conditions,
fixed and random parameters, Fang 2013
https://link.springer.com/article/10.1007/s00190-013-0643-2

# Arguments
- `A` Design matrix
- `y` RHS
- `Qaa` Covariance matrix of `e = vec(E)`
- `Qay` Covariance between A and y

# Keyword Arguments
- ` iters = 10` Number of iterations
"""
function wtls(A,y,Qaa,Qay,Qyy; iters = 10)
    n,u    = size(A)
    eyef   = issparse(Qaa) ? n->SparseMatrixCSC(1I,n,n) : Eye
    Iₙ,Iᵤ  = eyef(n),eyef(u)
    x      = wls(A,y,Qyy) # Initialize with weigted LS slution
    QₐₐQₐy = [Qaa Qay]
    QΠ     = [Qaa Qay; [Qay' Qyy]]
    for i = 1:iters
        B    = [(x' ⊗ Iₙ) -Iₙ]
        BQBᵀ = factorize(Symmetric(B*QΠ*B'))
        λ    = BQBᵀ\(y-A*x)
        v    = QₐₐQₐy*B'λ
        x    = (A'*(BQBᵀ\A))\((Iᵤ ⊗ λ')*v + A'*(BQBᵀ\y))
    end
    x
end

"""
    Qaa,Qay,Qyy = rowcovariance(rowQ::AbstractVector{<:AbstractMatrix})

Takes row-wise covariance matrices `QAy[i]` and returns the full (sparse) covariance matrices. `rowQ = [cov([A[i,:] y[i]]) for i = 1:length(y)]`
"""
function rowcovariance(rowQ::AbstractVector{<:AbstractMatrix})
    n = length(rowQ)
    u = size(rowQ[1],1)-1
    Qaa,Qay,Qyy = spzeros(n*u,n*u), spzeros(n*u,n), Diagonal(spzeros(n,n))
    for (i,Q) = enumerate(rowQ)
        for col = 1:u
            for row = 1:u
                aind1 = i + (row-1)*n
                aind2 = i + (col-1)*n
                Qaa[aind1,aind2] = Q[row,col]
            end
        end
        yind = i
        for row = 1:u
            aind = i + (row-1)*n
            Qay[aind,yind] = Q[row,end]
        end
        Qyy[yind,yind] = Q[end,end]
    end
    Qaa,Qay,Qyy
end

include("robustPCA.jl")

"""
    rtls(A,y; nukeA=false, kwargs...)

Solves a robust total least-squares problem Ax=y using the robust PCA method of "The Augmented Lagrange Multiplier Method for Exact Recovery of Corrupted Low-Rank Matrices", https://people.eecs.berkeley.edu/~yima/psfile/Lin09-MP.pdf

# Arguments
- `A` Design matrix
- `y` RHS
- `kwargs` are passed to [`rpca`](@ref), see the docstring for [`rpca`](@ref) for more help.
"""
function rtls(A::AbstractArray,y::AbstractArray; kwargs...)
    AA  = [A y]
    Ah,Eh,s = rpca(AA; nukeA=false, kwargs...)
    tls!(s,size(A,2))
end

# function rtls_ga(A::AbstractArray,y::AbstractArray; μ = μ!, kwargs...)
#     AA  = [A y]
#     m = mean(AA, dims=1)
#     AA .-= m
#     V = rpca_ga(AA'; μ = μ, kwargs...)
#     n   = size(A,2)
#     V21 = V[1:n,n+1:end]
#     V22 = V[n+1:end,n+1:end]
#     x   = -V21/V22
# end

"""
    irls(A, y; tolx=0.001, tol=1.0e-6, verbose=false, iters=100)

Iteratively reweighted least squares. Solves
minimizeₓ ||Ax-y||₁

#Arguments:
- `A`: Design matrix
- `y`: Measurements
- `tolx`: Minimum change in `x` before quitting
- `tol`: Minimum change in error before quitting
- `iters`: Maximum number of iterations
"""
function irls(A,y; tolx=1e-4, tol=1e-6, verbose=false, iters=100)
    x = A\y
    zA = similar(A)
    w = similar(y)
    mul!(w,A,x) # w used for storage
    xold = copy(x)
    e = abs.(y - w)
    me = median(e)
    meold = me
    verbose && println("0: Error: ", round(me, sigdigits=3))
    w .= max.(e, 1e-5)
    for i = 1:iters
        x = wls!(zA,A,y,w)
        mul!(w,A,x) # w used for storage
        e .= abs.(y .- w)
        w .= max.(e, 1e-5)
        d = norm(x-xold)/norm(x)
        me = median(e)
        verbose && println("$i: Error: ", round(me, sigdigits=3), " relative step: ", round(d, sigdigits=3))
        if d < tolx || meold-me < tol
            verbose && println("Success")
            return x
        end
        meold = me
        xold .= x
    end
    verbose && println("Maximum number of iterations reached")
    x
end

"""
    sls(A, y; r = 1, iters = 100, verbose = false, tol = 1.0e-8, α0 = 0.1)

Simplex least-squares: minimizeₓ ||Ax-y||₂ s.t. sum(x) = r

# Arguments:
- `A`: Design matrix
- `y`: RHS
- `r`: Simplex radius. the default (1) is the probability simplex.
- `iters`: Maximum number of iterations of projected gradient
- `verbose`: Print stuff
- `tol`: Tolerance (change in x between iterations).
- `α0`: initial step size. This parameter influences speed of convergence and also to which point the algorithm converges.
"""
function sls(A, y; r=1, iters=100, verbose=false, tol=1e-8, α0 = 0.1)

    x = A\y
    proj_simplex!(x; r=r)
    xo = copy(x)
    g = similar(x)
    e = similar(y)
    verbose && @info "Iter 0 cost: $(norm(A*x-y))"
    for iter = 1:iters
        mul!(e,A,x)
        e .= y .- e
        mul!(g, A', e)
        ng = norm(g)
        α = α0 / sqrt(iter)
        x .+= α .* g
        proj_simplex!(x; r=r)
        step = sqrt(sum(abs2(x-xo) for (x,xo) in zip(x,xo)))
        verbose && @info "Iter $iter norm(g): $ng norm(x-xo): $step, cost: $(norm(e))"
        step < tol && break
        xo .= x

    end
    verbose &&  @info "Converged - cost function: $(norm(e))"
    x
end

"""
    proj_simplex!(x; iters = 1000, r = 1, tol = 1.0e-8)

Project x onto the simplex with radius `r` such that `sum(x) = r` and all `x >= 0`
"""
function proj_simplex!(x; iters=1000, r=1, tol=1e-8)
    μ = minimum(x) - r
    for iter = 1:iters
        cost = sum(max(x - μ, 0)  for x in x) - r
        df   = sum(-((x - μ) > 0) for x in x)
        μ   -= cost / df
        abs(cost) < tol && break
    end
    @. x = max(x - μ, 0)
end

end # module


text:
using Deepgroebner
using Test

@testset "Deepgroebner.jl" begin
    f = [term(1.0, (1, 2, 3)), term(3.0, (3, 2, 1))]
    g = [term(4.0, (1, 2, 3)), term(1.0, (1, 2, 0))]
    h = [term(2.0, (1, 2, 3)), term(3.0, (3, 2, 1))]
    t = term(2.0, (1, 2, 2))
    s = term(4.0, (2, 2, 1))
    @test LT(f) == term(3.0, (3, 2, 1))
    @test ! Deepgroebner.div(t, LT(f))
    @test   Deepgroebner.div(s, LT(f))
    @test minus(f, g) == [term(-3.0, (1, 2, 3)), term(-1.0, (1, 2, 0)), term(3.0, (3, 2, 1))]
    @test mdiv([term(1.0, (1, 2, 3)), term(1.0, (0, 1, 1))], [[term(1.0, (1, 0, 0))]]) == [term(1.0, (0, 1, 1))]
end


text:
function QSLLT_lautat(surf :: ThreeDSurfSimple, field :: ThreeDFieldSimple, nsteps :: Int64 = 500, dtstar :: Float64 = 0.015, startflag = 0, writeflag = 0, writeInterval = 1000., delvort = delNone(); maxwrite = 50, nround = 6)

    # If a restart directory is provided, read in the simulation data
    if startflag == 0
        mat = zeros(0, 4+surf.nspan*8)
        t = 0.
    elseif startflag == 1
        dirvec = readdir()
        dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec)
        latestTime = maximum(dirresults)
        mat = DelimitedFiles.readdlm("resultsSummary")
        t = mat[end,1]
    else
        throw("invalid start flag, should be 0 or 1")
    end
    mat = mat'

    dt = dtstar*surf.cref/surf.uref

    # if writeflag is on, determine the timesteps to write at
    if writeflag == 1
        writeArray = Int64[]
        tTot = nsteps*dt
        for i = 1:maxwrite
            tcur = writeInterval*real(i)
            if t > tTot
                break
            else
                push!(writeArray, Int(round(tcur/dt)))
            end
        end
    end

    cl = zeros(surf.nspan)
    cd = zeros(surf.nspan)
    cm = zeros(surf.nspan)

    #Initialise flowfield
    for is = 1:surf.nspan
        push!(field.f2d, TwoDFlowField())
    end

    for istep = 1:nsteps
        #Udpate current time
        t = t + dt

        for is = 1:surf.nspan
            #Update external flowfield
            update_externalvel(field.f2d[is], t)

            #Update kinematic parameters
            update_kinem(surf.s2d[is], t)

            #Update bound vortex positions
            update_boundpos(surf.s2d[is], dt)

            #Add a TEV with dummy strength
            place_tev(surf.s2d[is], field.f2d[is], dt)

            update_indbound(surf.s2d[is], field.f2d[is])
        end

        kelv = KelvinConditionQSLLT(surf, field)

        #Solve for TEV strength to satisfy Kelvin condition

        soln = nlsolve(not_in_place(kelv), 0.0001*ones(surf.nspan))

        for is = 1:surf.nspan
            field.f2d[is].tev[end].s = soln.zero[is]
        end

        calc_a2toan3d(surf)

        for is = 1:surf.nspan
            #Update 3D effect on A0 and A0dot
            surf.s2d[is].a0[1] += surf.a03d[is]

            #Update rest of Fourier terms
            update_a2toan(surf.s2d[is])

            #Update 3D effect on An
            for ia = 1:surf.naterm
                surf.s2d[is].aterm[ia] += surf.aterm3d[ia,is]
            end

            #Update derivatives
            update_adot(surf.s2d[is],dt)

            #Set previous values of aterm to be used for derivatives in next time step
            surf.s2d[is].a0prev[1] = surf.s2d[is].a0[1]
            for ia = 1:3
                surf.s2d[is].aprev[ia] = surf.s2d[is].aterm[ia]
            end

            #Calculate bound vortex strengths
            update_bv(surf.s2d[is])
        end

        #Wake rollup
        wakeroll(surf, field, dt)

        cl3d, cd3d, cm3d, cl, cd, cm = calc_forces(surf, field, dt)

        # write flow details if required
        if writeflag == 1
            if istep in writeArray
                dirname = "$(round(t, digits=nround))"
                writeStamp(dirname, t, surf, field)
            end
        end

        matvect = [t, cl3d, cd3d, cm3d]
        for is = 1:surf.nspan
            matvect = [matvect; [surf.s2d[is].kinem.alpha, surf.s2d[is].kinem.h,
                                 surf.s2d[is].kinem.u, surf.s2d[is].a0[1], cl[is]*surf.s2d[is].c/surf.cref, cd[is]*surf.s2d[is].c/surf.cref, cm[is]*surf.s2d[is].c*surf.s2d[is].c/(surf.cref*surf.cref), surf.a03d[is]]]
        end
        mat = hcat(mat,matvect)
    end
    mat = mat'

    f = open("resultsSummary", "w")
    Serialization.serialize(f, ["#time \t", "CL3D \t", "CD3D \t", "CM3D \t", "alpha -1 (deg) \t", "h/c -1 \t", "u/uref -1 \t", "A0 -1 \t", "Cl -1 \t", "Cd -2\t", "Cm -2 \t", "A03D-1 \t", "alpha -2 ...\n"])
    DelimitedFiles.writedlm(f, mat)
    close(f)

    mat, surf, field, cl, cd, cm

end

function QSLLT_ldvm(surf :: ThreeDSurfSimple, field :: ThreeDFieldSimple, nsteps :: Int64 = 500, dtstar :: Float64 = 0.015, startflag = 0, writeflag = 0, writeInterval = 1000., delvort = delNone(); maxwrite = 50, nround = 6)

    # If a restart directory is provided, read in the simulation data
    if startflag == 0
        mat = zeros(0, 4+surf.nspan*8)
        t = 0.
    elseif startflag == 1
        dirvec = readdir()
        dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec)
        latestTime = maximum(dirresults)
        mat = DelimitedFiles.readdlm("resultsSummary")
        t = mat[end,1]
    else
        throw("invalid start flag, should be 0 or 1")
    end
    mat = mat'

    dt = dtstar*surf.cref/surf.uref

    # if writeflag is on, determine the timesteps to write at
    if writeflag == 1
        writeArray = Int64[]
        tTot = nsteps*dt
        for i = 1:maxwrite
            tcur = writeInterval*real(i)
            if t > tTot
                break
            else
                push!(writeArray, Int(round(tcur/dt)))
            end
        end
    end

    cl = zeros(surf.nspan)
    cd = zeros(surf.nspan)
    cm = zeros(surf.nspan)
    lesp = zeros(surf.nspan)

    #Initialise flowfield
    for is = 1:surf.nspan
        push!(field.f2d, TwoDFlowField())
    end

    for istep = 1:nsteps
        #Udpate current time
        t = t + dt

        for is = 1:surf.nspan
            #Update external flowfield
            update_externalvel(field.f2d[is], t)

            #Update kinematic parameters
            update_kinem(surf.s2d[is], t)

            #Update bound vortex positions
            update_boundpos(surf.s2d[is], dt)

            #Add a TEV with dummy strength
            place_tev(surf.s2d[is], field.f2d[is], dt)

            update_indbound(surf.s2d[is], field.f2d[is])
        end

        kelv = KelvinConditionQSLLT(surf, field)

        #Solve for TEV strength to satisfy Kelvin condition

        soln = nlsolve(not_in_place(kelv), 0.0001*ones(surf.nspan))

        for is = 1:surf.nspan
            field.f2d[is].tev[end].s = soln.zero[is]
        end

        for is = 1:surf.nspan
            lesp[is] = surf.s2d[is].a0[1]
        end

        shed_ind = Int[]

        #Check for LESP condition
        for is = 1:surf.nspan
            if abs(lesp[is]) > surf.s2d[is].lespcrit[1]
                push!(shed_ind, is)
            end
        end
        nshed = length(shed_ind)

        #println(map(q->q.a0[1], surf.s2d))
        #println(shed_ind)
        #error("h")

        if nshed > 0
            #2D iteration at section if LESP_crit is exceeded
            #Add LEVs with dummy strength
            for is = 1:surf.nspan
                if is in shed_ind
                    place_lev(surf.s2d[is], field.f2d[is], dt)
                end
            end

            kelvkutta = KelvinKuttaQSLLT(surf, field, shed_ind)
            soln = nlsolve(not_in_place(kelvkutta), [-0.01*ones(surf.nspan); 0.01*ones(nshed)])

            levcount = 0
            for is = 1:surf.nspan
                if is in shed_ind
                    field.f2d[is].tev[end].s = soln.zero[is]
                    levcount +=1
                    field.f2d[is].lev[end].s = soln.zero[levcount+surf.nspan]
                    surf.s2d[is].levflag[1] = 1
                else
                    field.f2d[is].tev[end].s = soln.zero[is]
                    surf.s2d[is].levflag[1] = 0
                end
            end
        else
            for is = 1:surf.nspan
                surf.s2d[is].levflag[1] = 0
            end
        end

        calc_a2toan3d(surf)

        for is = 1:surf.nspan
            #Update 3D effect on A0 and A0dot
            surf.s2d[is].a0[1] += surf.a03d[is]

            #Update rest of Fourier terms
            update_a2toan(surf.s2d[is])

            #Update 3D effect on An
            for ia = 1:surf.naterm
                surf.s2d[is].aterm[ia] += surf.aterm3d[ia,is]
            end

            #Update derivatives
            update_adot(surf.s2d[is],dt)

            #Set previous values of aterm to be used for derivatives in next time step

            surf.s2d[is].a0prev[1] = surf.s2d[is].a0[1]
            for ia = 1:3
                surf.s2d[is].aprev[ia] = surf.s2d[is].aterm[ia]
            end

            #Calculate bound vortex strengths
            update_bv(surf.s2d[is])

            # Delete or merge vortices if required
            controlVortCount(delvort, surf.s2d[is].bnd_x[Int(round(surf.ndiv/2))], surf.s2d[is].bnd_z[Int(round(surf.s2d[is].ndiv/2))], field.f2d[is])

            #Wake rollup
            wakeroll(surf.s2d[is], field.f2d[is], dt)
        end

        cl3d, cd3d, cm3d, cl, cd, cm = calc_forces(surf, field, dt)

        # write flow details if required
        if writeflag == 1
            if istep in writeArray
                dirname = "$(round(t, digits=nround))"
                writeStamp(dirname, t, surf, field)
            end
        end

        matvect = [t, cl3d, cd3d, cm3d]
        for is = 1:surf.nspan
            matvect = [matvect; [surf.s2d[is].kinem.alpha, surf.s2d[is].kinem.h,
                                 surf.s2d[is].kinem.u, surf.s2d[is].a0[1], cl[is]*surf.s2d[is].c/surf.cref, cd[is]*surf.s2d[is].c/surf.cref, cm[is]*surf.s2d[is].c*surf.s2d[is].c/(surf.cref*surf.cref), surf.a03d[is]]]
        end
        mat = hcat(mat,matvect)
    end
    mat = mat'

    f = open("resultsSummary", "w")
    Serialization.serialize(f, ["#time \t", "CL3D \t", "CD3D \t", "CM3D \t", "alpha -1 (deg) \t", "h/c -1 \t", "u/uref -1 \t", "A0 -1 \t", "Cl -1 \t", "Cd -2\t", "Cm -2 \t", "A03D-1 \t", "alpha -2 ...\n"])
    DelimitedFiles.writedlm(f, mat)
    close(f)

    mat, surf, field, cl, cd, cm

end

# function QSLLT_ldvm(surf :: ThreeDSurfSimple, field :: ThreeDFieldSimple, nsteps :: Int64 = 500, dtstar :: Float64 = 0.015, startflag = 0, writeflag = 0, writeInterval = 1000., delvort = delNone(); maxwrite = 50, nround = 6)

#     # If a restart directory is provided, read in the simulation data
#     if startflag == 0
#         mat = zeros(0, 4+surf.nspan*8)
#         t = 0.
#     elseif startflag == 1
#         dirvec = readdir()
#         dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec)
#         latestTime = maximum(dirresults)
#         mat = DelimitedFiles.readdlm("resultsSummary")
#         t = mat[end,1]
#     else
#         throw("invalid start flag, should be 0 or 1")
#     end
#     mat = mat'

#     # if writeflag is on, determine the timesteps to write at
#     if writeflag == 1
#         writeArray = Int64[]
#         tTot = nsteps*dtstar
#         for i = 1:maxwrite
#             tcur = writeInterval*real(i)
#             if t > tTot
#                 break
#             else
#                 push!(writeArray, Int(round(tcur/dtstar)))
#             end
#         end
#     end

#     dt = dtstar*surf.cref/surf.uref

#     cl = zeros(surf.nspan)
#     cd = zeros(surf.nspan)
#     cm = zeros(surf.nspan)
#     lesp = zeros(surf.nspan)
#     levstr = zeros(surf.nspan)

#     #Initialise flowfield
#     for is = 1:surf.nspan
#         push!(field.f2d, TwoDFlowField())
#     end

#     for istep = 1:nsteps
#         #Udpate current time
#         t = t + dt

#         for is = 1:surf.nspan
#             #Update external flowfield
#             update_externalvel(field.f2d[is], t)

#             #Update kinematic parameters
#             update_kinem(surf.s2d[is], t)

#             #Update bound vortex positions
#             update_boundpos(surf.s2d[is], dt)

#             #Add a TEV with dummy strength
#             place_tev(surf.s2d[is], field.f2d[is], dt)

#             update_indbound(surf.s2d[is], field.f2d[is])
#         end

#         kelv = KelvinConditionQSLLT(surf, field)

#         #Solve for TEV strength to satisfy Kelvin condition

#         soln = nlsolve(not_in_place(kelv), 0.0001*ones(surf.nspan))

#         for is = 1:surf.nspan
#             field.f2d[is].tev[end].s = soln.zero[is]
#         end

#         calc_a2toan3d(surf)

#         for is = 1:surf.nspan
#             #Update 3D effect on A0 and A0dot
#             surf.s2d[is].a0[1] += surf.a03d[is]

#             #Update rest of Fourier terms
#             update_a2toan(surf.s2d[is])

#             #Update 3D effect on An
#             for ia = 1:surf.naterm
#                 surf.s2d[is].aterm[ia] += surf.aterm3d[ia,is]
#             end

#             #Update derivatives
#             update_adot(surf.s2d[is],dt)

#             lesp[is] = surf.s2d[is].a0[1]
#         end

#         shed_ind = Int[]

#         #Check for LESP condition
#         for is = 1:surf.nspan
#             if abs(lesp[is]) > surf.s2d[is].lespcrit[1]
#                 push!(shed_ind, is)
#             end
#         end
#         nshed = length(shed_ind)

#         #println(map(q->q.a0[1], surf.s2d))
#         #println(shed_ind)
#         #error("h")

#         if nshed > 0
#             #2D iteration at section if LESP_crit is exceeded
#             #Add LEVs with dummy strength
#             for is = 1:surf.nspan
#                 if is in shed_ind
#                     place_lev(surf.s2d[is], field.f2d[is], dt)
#                 else
#                     push!(field.f2d[is].lev, TwoDVort(0., 0., 0., 0., 0., 0.))
#                 end
#             end

#             #Start global iteration
#             iter = 0
#             tev_previter = zeros(surf.nspan)
#             while true
#                 iter += 1

#                 for is = 1:surf.nspan
#                     if is in shed_ind
#                         levstr[is] -= field.f2d[is].lev[end].s
#                         kelvkutta = KelvinKuttaQSLLT_sep(surf, field, is)
#                         soln = nlsolve(not_in_place(kelvkutta), [-0.01; 0.01])
#                         field.f2d[is].tev[end].s = soln.zero[1]
#                         field.f2d[is].lev[end].s = soln.zero[2]
#                         levstr[is] += field.f2d[is].lev[end].s
#                     else
#                         kelv = KelvinConditionQSLLT_sep(surf, field, is)
#                         soln = nlsolve(not_in_place(kelv), [-0.01;])
#                         field.f2d[is].tev[end].s = soln.zero[1]
#                         levstr[is] = 0.
#                     end
#                 end

#                 #Update cumulative LEV strength
#                 calc_a0a13d_wlev2(surf, levstr)

#                 if sum(abs.(map(q->q.tev[end].s, field.f2d) - tev_previter)) < 1e-4
#                     break
#                 end
#                 tev_previter = map(q->q.tev[end].s, field.f2d)
#             end
#             println(iter)

#             calc_a2toan3d(surf)

#             for is = 1:surf.nspan
#                 #Update 3D effect on A0 and A0dot
#                 surf.s2d[is].a0[1] += surf.a03d[is]

#                 #Update rest of Fourier terms
#                 update_a2toan(surf.s2d[is])

#                 #Update 3D effect on An
#                     for ia = 1:surf.naterm
#                     surf.s2d[is].aterm[ia] += surf.aterm3d[ia,is]
#                 end
#             end

#         else
#             for is = 1:surf.nspan
#                 surf.s2d[is].levflag[1] = 0
#             end
#         end

#         #Set previous values of aterm to be used for derivatives in next time step
#         for is = 1:surf.nspan
#             surf.s2d[is].a0prev[1] = surf.s2d[is].a0[1]
#             for ia = 1:3
#                 surf.s2d[is].aprev[ia] = surf.s2d[is].aterm[ia]
#             end

#             #Calculate bound vortex strengths
#             update_bv(surf.s2d[is])

#             # Delete or merge vortices if required
#             controlVortCount(delvort, surf.s2d[is].bnd_x[Int(round(surf.ndiv/2))], surf.s2d[is].bnd_z[Int(round(surf.s2d[is].ndiv/2))], field.f2d[is])

#             #Wake rollup
#             wakeroll(surf.s2d[is], field.f2d[is], dt)
#         end

#         cl3d, cd3d, cm3d, cl, cd, cm = calc_forces(surf)

#         # write flow details if required
#         if writeflag == 1
#             if istep in writeArray
#                 dirname = "$(round(t, digits=nround))"
#                 writeStamp(dirname, t, surf, field)
#             end
#         end

#         matvect = [t, cl3d, cd3d, cm3d]
#         for is = 1:surf.nspan
#             matvect = [matvect; [surf.s2d[is].kinem.alpha, surf.s2d[is].kinem.h,
#                                  surf.s2d[is].kinem.u, surf.s2d[is].a0[1], cl[is]*surf.s2d[is].c/surf.cref, cd[is]*surf.s2d[is].c/surf.cref, cm[is]*surf.s2d[is].c*surf.s2d[is].c/(surf.cref*surf.cref), surf.a03d[is]]]
#         end
#         mat = hcat(mat,matvect)
#     end
#     mat = mat'

#     f = open("resultsSummary", "w")
#     Serialization.serialize(f, ["#time \t", "CL3D \t", "CD3D \t", "CM3D \t", "alpha -1 (deg) \t", "h/c -1 \t", "u/uref -1 \t", "A0 -1 \t", "Cl -1 \t", "Cd -2\t", "Cm -2 \t", "A03D-1 \t", "alpha -2 ...\n"])
#     DelimitedFiles.writedlm(f, mat)
#     close(f)

#     mat, surf, field, cl, cd, cm

# end

function stripLDVM(surf :: ThreeDSurfSimple, field :: ThreeDFieldSimple, nsteps :: Int64 = 500, dtstar :: Float64 = 0.015, startflag = 0, writeflag = 0, writeInterval = 1000., delvort = delNone(); maxwrite = 50, nround = 6)

    # If a restart directory is provided, read in the simulation data
    if startflag == 0
        mat = zeros(0, 4+surf.nspan*8)
        t = 0.
    elseif startflag == 1
        dirvec = readdir()
        dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec)
        latestTime = maximum(dirresults)
        mat = DelimitedFiles.readdlm("resultsSummary")
        t = mat[end,1]
    else
        throw("invalid start flag, should be 0 or 1")
    end
    mat = mat'

    dt = dtstar*surf.cref/surf.uref

    # if writeflag is on, determine the timesteps to write at
    if writeflag == 1
        writeArray = Int64[]
        tTot = nsteps*dt
        for i = 1:maxwrite
            tcur = writeInterval*real(i)
            if t > tTot
                break
            else
                push!(writeArray, Int(round(tcur/dt)))
            end
        end
    end

    cl = zeros(surf.nspan)
    cd = zeros(surf.nspan)
    cm = zeros(surf.nspan)
    lesp = zeros(surf.nspan)

    #Initialise flowfield
    for is = 1:surf.nspan
        push!(field.f2d, TwoDFlowField())
    end

    for istep = 1:nsteps
        #Udpate current time
        t = t + dt

        for is = 1:surf.nspan
            #Update external flowfield
            update_externalvel(field.f2d[is], t)

            #Update kinematic parameters
            update_kinem(surf.s2d[is], t)

            #Update bound vortex positions
            update_boundpos(surf.s2d[is], dt)

            #Add a TEV with dummy strength
            place_tev(surf.s2d[is], field.f2d[is], dt)

            update_indbound(surf.s2d[is], field.f2d[is])

            kelv = KelvinCondition(surf.s2d[is], field.f2d[is])

            #Solve for TEV strength to satisfy Kelvin condition

            soln = nlsolve(not_in_place(kelv), [-0.01;])

            field.f2d[is].tev[end].s = soln.zero[1]

            #Update adot
            update_a2a3adot(surf.s2d[is],dt)

            lesp[is] = surf.s2d[is].a0[1]

            #Check for LESP condition
            if abs(lesp[is]) > surf.s2d[is].lespcrit[1]

                place_lev(surf.s2d[is], field.f2d[is], dt)

                kelvkutta = KelvinKutta(surf.s2d[is], field.f2d[is])
                soln = nlsolve(not_in_place(kelvkutta), [-0.01; 0.01])

                field.f2d[is].tev[end].s = soln.zero[1]

                field.f2d[is].lev[end].s = soln.zero[2]

                surf.s2d[is].levflag[1] = 1
            else
                surf.s2d[is].levflag[1] = 0
            end

            #Update rest of Fourier terms
            update_a2toan(surf.s2d[is])

            #Set previous values of aterm to be used for derivatives in next time step

            surf.s2d[is].a0prev[1] = surf.s2d[is].a0[1]
            for ia = 1:3
                surf.s2d[is].aprev[ia] = surf.s2d[is].aterm[ia]
            end

            #Calculate bound vortex strengths
            update_bv(surf.s2d[is])

            # Delete or merge vortices if required
            #controlVortCount(delvort, surf.s2d[is].bnd_x[Int(round(surf.ndiv/2))], surf.s2d[is].bnd_z[Int(round(surf.s2d[is].ndiv/2))], field.f2d[is])

            #Wake rollup
            wakeroll(surf.s2d[is], field.f2d[is], dt)

            cl[is], cd[is], cm[is] = calc_forces(surf.s2d[is])
        end

        cl3d = 0; cd3d = 0; cm3d = 0

        for is = 1:surf.nspan
            if is == 1
                yl = -0.5*surf.AR*surf.cref
                yu = 0.5*(surf.yle[is] + surf.yle[is+1])
            elseif is == surf.nspan
                yl = 0.5*(surf.yle[is] + surf.yle[is-1])
                yu = 0.
            else
                yl = 0.5*(surf.yle[is] + surf.yle[is-1])
                yu = 0.5*(surf.yle[is] + surf.yle[is+1])
            end
            cl3d += cl[is]*(yu-yl)
            cd3d += cd[is]*(yu-yl)
            cm3d += cm[is]*(yu-yl)
        end
        cl3d = cl3d/(0.5*surf.AR*surf.cref)
        cd3d = cd3d/(0.5*surf.AR*surf.cref)
        cm3d = cm3d/(0.5*surf.AR*surf.cref)

        # write flow details if required
        if writeflag == 1
            if istep in writeArray
                dirname = "$(round(t, digits=nround))"
                writeStamp(dirname, t, surf, field)
            end
        end

        matvect = [t, cl3d, cd3d, cm3d]
        for is = 1:surf.nspan
            matvect = [matvect; [surf.s2d[is].kinem.alpha, surf.s2d[is].kinem.h,
                                 surf.s2d[is].kinem.u, surf.s2d[is].a0[1], cl[is]*surf.s2d[is].c/surf.cref, cd[is]*surf.s2d[is].c/surf.cref, cm[is]*surf.s2d[is].c*surf.s2d[is].c/(surf.cref*surf.cref), surf.a03d[is]]]
        end
        mat = hcat(mat,matvect)
    end
    mat = mat'

    f = open("resultsSummary", "w")
    Serialization.serialize(f, ["#time \t", "CL3D \t", "CD3D \t", "CM3D \t", "alpha -1 (deg) \t", "h/c -1 \t", "u/uref -1 \t", "A0 -1 \t", "Cl -1 \t", "Cd -2\t", "Cm -2 \t", "A03D-1 \t", "alpha -2 ...\n"])
    DelimitedFiles.writedlm(f, mat)
    close(f)

    mat, surf, field, cl, cd, cm

end

function stripCantilever(surf::ThreeDSurfSimple, field::ThreeDFieldSimple, data::Dict{String,Any}, kinem::KinemCantilever, nsteps::Int64 = 500, dtstar::Float64 = 0.015, startflag = 0, writeflag = 0, writeInterval = 1000., delvort = delNone(); maxwrite = 50, nround = 6)

    # If a restart directory is provided, read in the simulation data
    if startflag == 0
        mat = zeros(0, 4+surf.nspan*8)
        t = 0.
    elseif startflag == 1
        dirvec = readdir()
        dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec)
        latestTime = maximum(dirresults)
        mat,_ = DelimitedFiles.readdlm("resultsSummary", header=true)
        t = mat[end,1]
    else
        throw("invalid start flag, should be 0 or 1")
    end
    mat = mat'

    dt = dtstar*surf.cref/surf.uref

    # if writeflag is on, determine the timesteps to write at
    if writeflag == 1
        writeArray = Int64[]
        tTot = nsteps*dt
        for i = 1:maxwrite
            tcur = writeInterval*real(i)
            if t > tTot
                break
            else
                push!(writeArray, Int(round(tcur/dt)))
            end
        end
    end

    cl = zeros(surf.nspan)
    cd = zeros(surf.nspan)
    cm = zeros(surf.nspan)
    lesp = zeros(surf.nspan)

    #Extract required data
    ulist = get(data, "U", "fail")
    ulist = ulist[1,:]
    unum = findfirst(isequal(surf.uref), ulist)
    mmod = get(data, "mmod", "fail")

    m_eta_s_etadd = get(mmod, "m_eta_s_etadd", "fail")
    m_eta_s_etadd = m_eta_s_etadd[unum]
    m_eta_s_etad = get(mmod, "m_eta_s_etad", "fail")
    m_eta_s_etad = m_eta_s_etad[unum]
    m_eta_s_eta = get(mmod, "m_eta_s_eta", "fail")
    m_eta_s_eta = m_eta_s_eta[unum]
    Phi = get(data, "Phi", "fail")

    #Initialise flowfield
    for is = 1:surf.nspan
        push!(field.f2d, TwoDFlowField())
    end

    for istep = 1:nsteps
        #Udpate current time
        t = t + dt

        for is = 1:surf.nspan
            #Update external flowfield
            update_externalvel(field.f2d[is], t)
        end

        #Update kinematic parameters
        update_kinemCantilever(surf, dt, kinem, Phi, m_eta_s_etadd, m_eta_s_etad, m_eta_s_eta, cl, cm)

        for is = 1:surf.nspan
            #Update bound vortex positions
            update_boundpos(surf.s2d[is], dt)

            #Add a TEV with dummy strength
            place_tev(surf.s2d[is], field.f2d[is], dt)

            update_indbound(surf.s2d[is], field.f2d[is])

            kelv = KelvinCondition(surf.s2d[is], field.f2d[is])

            #Solve for TEV strength to satisfy Kelvin condition

            soln = nlsolve(not_in_place(kelv), [-0.01;])

            field.f2d[is].tev[end].s = soln.zero[1]

            #Update adot
            update_a2a3adot(surf.s2d[is],dt)

            lesp[is] = surf.s2d[is].a0[1]

            #Check for LESP condition
            if abs(lesp[is]) > surf.s2d[is].lespcrit[1]

                place_lev(surf.s2d[is], field.f2d[is], dt)

                kelvkutta = KelvinKutta(surf.s2d[is], field.f2d[is])
                soln = nlsolve(not_in_place(kelvkutta), [-0.01; 0.01])

                field.f2d[is].tev[end].s = soln.zero[1]

                field.f2d[is].lev[end].s = soln.zero[2]

                surf.s2d[is].levflag[1] = 1
            else
                surf.s2d[is].levflag[1] = 0
            end

            #Update rest of Fourier terms
            update_a2toan(surf.s2d[is])

            #Set previous values of aterm to be used for derivatives in next time step

            surf.s2d[is].a0prev[1] = surf.s2d[is].a0[1]
            for ia = 1:3
                surf.s2d[is].aprev[ia] = surf.s2d[is].aterm[ia]
            end

            #Calculate bound vortex strengths
            update_bv(surf.s2d[is])

            # Delete or merge vortices if required
            controlVortCount(delvort, surf.s2d[is].bnd_x[Int(round(surf.ndiv/2))], surf.s2d[is].bnd_z[Int(round(surf.s2d[is].ndiv/2))], field.f2d[is])

            #Wake rollup
            wakeroll(surf.s2d[is], field.f2d[is], dt)

            cl[is], cd[is], cm[is] = calc_forces(surf.s2d[is])
        end

        cl3d = 0; cd3d = 0; cm3d = 0

        for is = 1:surf.nspan
            if is == 1
                yl = -0.5*surf.AR*surf.cref
                yu = 0.5*(surf.yle[is] + surf.yle[is+1])
            elseif is == surf.nspan
                yl = 0.5*(surf.yle[is] + surf.yle[is-1])
                yu = 0.
            else
                yl = 0.5*(surf.yle[is] + surf.yle[is-1])
                yu = 0.5*(surf.yle[is] + surf.yle[is+1])
            end
            cl3d += cl[is]*(yu-yl)
            cd3d += cd[is]*(yu-yl)
            cm3d += cm[is]*(yu-yl)
        end
        cl3d = cl3d/(0.5*surf.AR*surf.cref)
        cd3d = cd3d/(0.5*surf.AR*surf.cref)
        cm3d = cm3d/(0.5*surf.AR*surf.cref)

        # write flow details if required
        if writeflag == 1
            if istep in writeArray
                dirname = "$(round(t, digits=nround))"
                writeStamp(dirname, t, surf, field)
            end
        end

        matvect = [t, cl3d, cd3d, cm3d]
        for is = 1:surf.nspan
            matvect = [matvect; [surf.s2d[is].kinem.alpha, surf.s2d[is].kinem.h,
                                 surf.s2d[is].kinem.u, surf.s2d[is].a0[1], cl[is]*surf.s2d[is].c/surf.cref, cd[is]*surf.s2d[is].c/surf.cref, cm[is]*surf.s2d[is].c*surf.s2d[is].c/(surf.cref*surf.cref), surf.a03d[is]]]
        end
        mat = hcat(mat,matvect)
    end
    mat = mat'

    f = open("resultsSummary", "w")
    Serialization.serialize(f, ["#time \t", "CL3D \t", "CD3D \t", "CM3D \t", "alpha -1 (deg) \t", "h/c -1 \t", "u/uref -1 \t", "A0 -1 \t", "Cl -1 \t", "Cd -2\t", "Cm -2 \t", "A03D-1 \t", "alpha -2 ...\n"])
    DelimitedFiles.writedlm(f, mat)
    close(f)

    mat, surf, field, cl, cd, cm

end


text:
"""
Samples from the posterior distribution
* `model`: a model containing likelihood function with data and priors
* `de`: differential evolution object
* `n_iter`: number of iterations or samples

Function signature
```@example
    sample(model::DEModel, de::DE, n_iter::Int; kwargs...)
```
"""
sample(model::DEModel, de::DE, n_iter::Int; kwargs...) = _sample(model::DEModel, de::DE, n_iter::Int; stepfun=step!, kwargs...)

function _sample(model::DEModel, de::DE, n_iter::Int; progress=false, stepfun=step!, kwargs...)
    meter = Progress(n_iter)
    # initialize particles based on prior distribution
    groups = sample_init(model, de, n_iter)
    for iter in 1:n_iter
        de.iter = iter
        # explicitly pass groups so parallel works
        groups = stepfun(model, de, groups)
        progress ? next!(meter) : nothing
    end
    # convert to chain object
    chain = bundle_samples(model, de, groups, n_iter)
    return chain
end

"""
Samples from the posterior distribution with each group of particles on a seperarate thread for
the mutation and crossover steps.
* `model`: a model containing likelihood function with data and priors
* `de`: differential evolution object
* `n_iter`: number of iterations or samples

Function signature
```@example
    psample(model::DEModel, de::DE, n_iter::Int; kwargs...)
```
"""
function sample(model::DEModel, de::DE, ::MCMCThreads, n_iter::Int; progress=false, kwargs...)
    _sample(model::DEModel, de::DE, n_iter::Int; progress=progress, stepfun=pstep!, kwargs...)
end

"""
Perform a single step for DE-MCMC.
* `model`: model containing a likelihood function with data and priors
* `de`: DE-MCMC sampler object
* `groups`: Array of vectors of particles
"""
function step!(model::DEModel, de::DE, groups)
    rand() <= de.α ? migration!(de, groups) : nothing
    groups = mutate_crossover!(model, de, groups)
    store_samples!(de, groups)
    return groups
end

function pstep!(model::DEModel, de::DE, groups)
    rand() <= de.α ? migration!(de, groups) : nothing
    groups = pmutate_crossover!(model, de, groups)
    store_samples!(de, groups)
    return groups
end

_f(x...; y...) = nothing

"""
On a single core, selects between mutation and crossover step with probability β
* `model`: model containing a likelihood function with data and priors
* `de`: differential evolution object
* `group`: a vector of interacting particles (e.g. chains)
"""
function mutate_crossover!(model, de, groups)
    seeds = rand(UInt, length(groups))
    groups = map((group,s) -> mutate_or_crossover!(model, de, group, s), groups, seeds)
    return groups
end

"""
On multiple cores, selects between mutation and crossover step with probability β
* `model`: model containing a likelihood function with data and priors
* `de`: differential evolution object
* `groups`: a vector of interacting particles (e.g. chains)
"""
function pmutate_crossover!(model, de, groups)
    seeds = rand(UInt, de.n_groups)
    Threads.@threads for g in 1:de.n_groups
        group = mutate_or_crossover!(model, de, groups[g], seeds[g])
    end
    return groups
end

"""
Selects between mutation and crossover step with probability β
* `model`: model containing a likelihood function with data and priors
* `de`: differential evolution object
* `group`: a vector of interacting particles (e.g. chains)
"""
function mutate_or_crossover!(model, de, group, seed)
    Random.seed!(seed)
    rand() <= de.β ? mutation!(model, de, group) : crossover!(model, de, group)
    return group
end

"""
Converts group particles to a chain object capable of generating convergence diagnostics
and posterior summaries
* `model`: model containing a likelihood function with data and priors
* `de`: differential evolution object
* `groups`: a vector of groups of particles
* `n_iter`: number of iterations
"""
function bundle_samples(model::DEModel, de::DE, groups, n_iter)
    particles = vcat(groups...)
    Np = length(particles)
    Ns = n_iter- de.burnin
    all_names = get_names(model, particles[1])
    n_parms = length(all_names)
    Nnames = length(model.names)
    v = fill(0.0, Ns, n_parms, Np)
    for (c,p) in enumerate(particles)
        for s in 1:Ns
            temp = Float64[]
            for ni in 1:Nnames
                push!(temp, p.samples[s,ni]...)
            end
            push!(temp, p.accept[s], p.lp[s])
            v[s,:,c] = temp'
        end
    end
    return Chains(v, all_names, (internals=["acceptance", "lp"],); sorted=true)
end

"""
Creates vectors of particles and samples initial parameter values from priors
* `model`: model containing a likelihood function with data and priors
* `de`: differential evolution object
* `n_iter`: number of iterations
"""
function sample_init(model::DEModel, de::DE, n_iter)
    groups = [[Particle(Θ=sample_prior(model.priors)) for p in 1:de.Np]
        for c in 1:de.n_groups]
    for group in groups
        for p in group
            init_particle!(model, de, p, n_iter)
        end
    end
    return groups
end


text:
"""
    rand(::Type{DQMC}, m::Model, nslices::Int)

Draw random configuration.
"""
Base.rand(::Type{DQMC}, m::Model, nslices::Int) = error("Model $(typeof(m)) doesn't implement `rand(::Type{DQMC}, m::Model, nslices::Int)`!")

"""
	hopping_matrix(mc::DQMC, m::Model)

Calculates the hopping matrix \$T_{i\\sigma, j\\sigma '}\$ where \$i, j\$ are site indices and \$\\sigma , \\sigma '\$
are flavor indices (e.g. spin indices). The hopping matrix should also contain
potential chemical potential terms on the diagonal.

A matrix element is the hopping amplitude for a hopping process: \$j,\\sigma ' \\rightarrow i,\\sigma\$.

Regarding the order of indices, if `T[i, σ, j, σ']` is your desired 4D hopping array, then `reshape(T, (n_sites * n_flavors, :))`
is the hopping matrix.
"""
hopping_matrix(mc::DQMC, m::Model) = error("Model has no implementation of `hopping_matrix(mc::DQMC, m::Model)`.")

"""
    interaction_matrix_exp!(mc::DQMC, m::Model, result::Matrix, conf, slice::Int, power::Float64=1.) -> nothing

Calculate the interaction matrix exponential `expV = exp(- power * delta_tau * V(slice))`
and store it in `result::Matrix`. Potential chemical potential terms should be
part of the `hopping_matrix` and not the interaction.

This is a performance critical method and one might consider efficient in-place (in `result`) construction.
"""
interaction_matrix_exp!(mc::DQMC, m::Model, result::Matrix, conf, slice::Int, power::Float64=1.) = error("Model has no implementation of `interaction_matrix_exp!(mc::DQMC, m::Model, result::Matrix, conf, slice::Int, power::Float64=1.)`.")

"""
    propose_local(mc::DQMC, m::Model, i::Int, conf) -> detratio, ΔE_boson, Δ

Propose a local move for lattice site `i` of current configuration `conf` . Returns the
Green's function determinant ratio, the boson energy difference `ΔE_boson = E_boson_new - E_boson`,
and potentially additional local move information `Δ` (will be forwarded to `accept_local!`).

See also [`accept_local!`](@ref).
"""
propose_local(mc::DQMC, m::Model, i::Int, conf) =
    error("Model has no implementation of `propose_local(mc::DQMC, m::Model, i::Int, conf)`!")

"""
    accept_local(mc::DQMC, m::Model, i::Int, slice::Int, conf, Δ, detratio, ΔE_boson)

Accept a local move for site `i` at imaginary time slice `slice` of current configuration `conf`.
Arguments `Δ`, `detratio` and `ΔE_boson` correspond to output of `propose_local` for that local move.

See also [`propose_local`](@ref).
"""
accept_local(mc::DQMC, m::Model, i::Int, slice::Int, conf, Δ, detratio, ΔE_boson) =
    error("Model has no implementation of `accept_local(mc::DQMC, m::Model, i::Int, slice::Int, conf, Δ, detratio, ΔE_boson)`!")

max_stars_repo_name: yosefcap/MC

max_issues_repo_name: yosefcap/MC

max_forks_repo_name: yosefcap/MC


text:
# ref: https://github.com/apache/parquet-format/blob/master/Encodings.md

macro bitwidth(i)
    quote
        ceil(Int, log(2, $(esc(i))+1))
    end
end
macro bit2bytewidth(i)
    quote
        ceil(Int, $(esc(i))/8)
    end
end
macro byt2itype(i)
    quote
        ($(esc(i)) <= 4) ? Int32 : ($(esc(i)) <= 8) ? Int64 : Int128
    end
end
macro byt2uitype(i)
    quote
        ($(esc(i)) <= 4) ? UInt32 : ($(esc(i)) <= 8) ? UInt64 : UInt128
    end
end
macro byt2uitype_small(i)
    quote
        ($(esc(i)) <= 1) ? UInt8 : ($(esc(i)) <= 2) ? UInt16 : ($(esc(i)) <= 4) ? UInt32 : ($(esc(i)) <= 8) ? UInt64 : UInt128
    end
end

const MSB = 0x80
const MASK7 = 0x7f
const MASK8 = 0xff
const MASK3 = 0x07
function MASKN(nbits::UInt8)
    byte_width = @bit2bytewidth(nbits)
    type_small = @byt2uitype_small(byte_width)
    MASKN(nbits, type_small)
end
function MASKN(nbits::UInt8, ::Type{T}) where {T}
    O = convert(T, 0x1)
    (O << nbits) - O
end

read_fixed(io::IO, typ::Type{UInt32}) = _read_fixed(io, convert(UInt32,0), 4)
read_fixed(io::IO, typ::Type{UInt64}) = _read_fixed(io, convert(UInt64,0), 8)
read_fixed(io::IO, typ::Type{Int32}) = reinterpret(Int32, _read_fixed(io, convert(UInt32,0), 4))
read_fixed(io::IO, typ::Type{Int64}) = reinterpret(Int64, _read_fixed(io, convert(UInt64,0), 8))
read_fixed(io::IO, typ::Type{Int128}) = reinterpret(Int128, _read_fixed(io, convert(UInt128, 0), 12))   # INT96: 12 bytes little endian
read_fixed(io::IO, typ::Type{Float32}) = reinterpret(Float32, _read_fixed(io, convert(UInt32,0), 4))
read_fixed(io::IO, typ::Type{Float64}) = reinterpret(Float64, _read_fixed(io, convert(UInt64,0), 8))
function _read_fixed(io::IO, ret::T, N::Int) where {T <: Unsigned}
    for n in 0:(N-1)
        byte = convert(T, read(io, UInt8))
        ret |= (byte << *(8,n))
    end
    ret
end

function _read_varint(io::IO, ::Type{T}) where {T <: Integer}
    res = zero(T)
    n = 0
    byte = UInt8(MSB)
    while (byte & MSB) != 0
        byte = read(io, UInt8)
        res |= (convert(T, byte & MASK7) << (7*n))
        n += 1
    end
    # in case of overflow, consider it as missing field and return default value
    if (n-1) > sizeof(T)
        #@debug("overflow reading $T. returning 0")
        return zero(T)
    end
    res
end

# parquet types:   BOOLEAN, INT32, INT64,  INT96,   FLOAT,  DOUBLE, BYTE_ARRAY, FIXED_LEN_BYTE_ARRAY
# enum values:           0,     1,     2,      3,       4,       5,          6,                    7
const PLAIN_JTYPES = (Bool, Int32, Int64, Int128, Float32, Float64,      UInt8,               UInt8)

# read plain encoding (PLAIN = 0)
function read_plain_byte_array(io)
    count = read_fixed(io, Int32)
    read!(io, Array{UInt8}(undef, count))
end

function plain_values_eltype(typ::Int32)
    if typ === _Type.BOOLEAN
        Bool
    elseif typ === _Type.BYTE_ARRAY
        Vector{UInt8}
    elseif typ === _Type.FIXED_LEN_BYTE_ARRAY
        error("not implemented") # this is likely same as BYTE_ARRAY for decoding purpose
    else
        PLAIN_JTYPES[typ+1]
    end
end

# read plain values or dictionary (PLAIN_DICTIONARY = 2)
function read_plain_values(io, count::Int32, typ::Int32)
    T = plain_values_eltype(typ)
    arr = Array{T}(undef, count)
    read_plain_values(io, count, typ, arr)
end
function read_plain_values(io, count::Int32, typ::Int32, arr::Vector{T}, offset::Int=0) where {T}
    #@debug("reading plain values", type=typ, count=count)
    if typ === _Type.BOOLEAN
        arr = read_bitpacked_booleans(io, count, arr)
    elseif typ === _Type.BYTE_ARRAY
        @inbounds for i in 1:count
            arr[i+offset] = read_plain_byte_array(io)
        end
    elseif typ === _Type.FIXED_LEN_BYTE_ARRAY
        error("not implemented") # this is likely same as BYTE_ARRAY for decoding purpose
    else
        @inbounds for i in 1:count
            arr[i+offset] = read_fixed(io, T)
        end
    end
    #@debug("read $(length(arr)) plain values")
    arr
end

function read_bitpacked_booleans(io::IO, count::Int32, arr::Vector{Bool}, offset::Int=0)::Vector{Bool}
    #@debug("reading bitpacked booleans", count)
    arrpos = 1
    bits = UInt8(0)
    bitpos = 9
    @inbounds while arrpos <= count
        if bitpos > 8
            bits = read(io, UInt8)
            #@debug("bits", bits, bitstring(bits))
            bitpos = 1
        end
        arr[arrpos+offset] = Bool(bits & 0x1)
        arrpos += 1
        bits >>= 1
        bitpos += 1
    end
    arr
end

# read rle dictionary (RLE_DICTIONARY = 8, or PLAIN_DICTIONARY = 2 in a data page)

function read_rle_dict(io, count::Int32)
    bits = read(io, UInt8)
    read_hybrid(io, count, Val{bits}(); read_len=false)
end

# read RLE or bit backed format (RLE = 3)
function read_hybrid(io, count::Int32, ::Val{W}; read_len::Bool=true) where {W}
    byte_width = @bit2bytewidth(W)
    typ = @byt2itype(byte_width)
    arr = Array{typ}(undef, count)

    read_hybrid(io, count, W, byte_width, arr, 0; read_len=read_len)
end
function read_hybrid(io, count::Int32, bits::UInt8, arr::Vector{T}, offset::Int=0; read_len::Bool=true) where {T}
    byte_width = @bit2bytewidth(bits)
    read_hybrid(io, count, bits, byte_width, arr, offset; read_len=read_len)
end
function read_hybrid(io, count::Int32, bits::UInt8, byt::Int, arr::Vector{T}, offset::Int=0; read_len::Bool=true) where {T}
    len = read_len ? read_fixed(io, Int32) : Int32(0)
    #@debug("reading hybrid data", len, count, bits)
    mask = MASKN(bits)
    arrpos = 1
    offset -= 1    # to counter arrpos starting at 1
    while arrpos <= count
        runhdr = _read_varint(io, Int)
        isbitpack = ((runhdr & 0x1) == 0x1)
        runhdr >>= 1
        nitems = min(isbitpack ? runhdr*8 : runhdr, count - arrpos + 1)

        if isbitpack
            read_bitpacked_run(io, runhdr, bits, byt, arr, offset+arrpos, mask)
        else # rle
            read_type = @byt2uitype(byt)
            read_rle_run(io, nitems, bits, byt, arr, read_type, offset+arrpos)
        end
        arrpos += nitems
    end
    arr
end

function read_rle_run(io, count::Int, bits::UInt8, byt::Int, arr::Vector{T}, read_type::Type{V}, offset::Int=0) where {T,V}
    #@debug("read_rle_run", count, T, bits, byt)
    val = reinterpret(T, _read_fixed(io, zero(V), byt))
    @assert length(arr) >= (count+offset)
    @inbounds for idx in 1:count
        arr[idx+offset] = val
    end
    arr
end

function read_bitpacked_run(io, grp_count::Int, ::Val{W}) where {W}
    byte_width = @bit2bytewidth(W)
    typ = @byt2itype(byte_width)
    arr = Array{typ}(undef, grp_count)
    read_bitpacked_run(io, grp_count, W, byte_width, arr)
end
function read_bitpacked_run(io, grp_count::Int, bits::UInt8, byt::Int, arr::Vector{T}, offset::Int=0, mask::V=MASKN(bits)) where {T,V}
    count = min(grp_count * 8, length(arr)-offset)
    # multiple of 8 values at a time are bit packed together
    nbytes = bits * grp_count # same as: round(Int, (bits * grp_count * 8) / 8)
    #@debug("read_bitpacked_run. grp_count:$grp_count, count:$count, nbytes:$nbytes, nbits:$bits, available:$(bytesavailable(io))")
    data = Array{UInt8}(undef, min(nbytes, bytesavailable(io)))
    read!(io, data)

    bitbuff = zero(V)
    nbitsbuff = UInt8(0)
    shift = UInt8(0)

    arridx = 1
    dataidx = 1
    while arridx <= count
        #@debug("arridx:$arridx nbitsbuff:$nbitsbuff shift:$shift bits:$bits")
        if nbitsbuff > 0
            # we have leftover bits, which must be appended
            if nbitsbuff < bits
                # but only append if we need to read more in this cycle
                @inbounds arr[arridx+offset] = bitbuff & MASKN(nbitsbuff, V)
                shift = nbitsbuff
                nbitsbuff = UInt8(0)
                bitbuff = zero(V)
            end
        end

        # fill buffer
        while (nbitsbuff + shift) < bits
             # shift 8 bits and read directly into bitbuff
            bitbuff |= (V(data[dataidx]) << nbitsbuff)
            dataidx += 1
            nbitsbuff += UInt8(8)
        end

        # set values
        while ((nbitsbuff + shift) >= bits) && (arridx <= count)
            if shift > 0
                remshift = bits - shift
                #@debug("setting part from bitbuff nbitsbuff:$nbitsbuff, shift:$shift, remshift:$remshift")
                arr[arridx+offset] |= convert(T, (bitbuff << shift) & mask)
                bitbuff >>= remshift
                nbitsbuff -= remshift
                shift = UInt8(0)
            else
                #@debug("setting all from bitbuff nbitsbuff:$nbitsbuff")
                arr[arridx+offset] = convert(T, bitbuff & mask)
                bitbuff >>= bits
                nbitsbuff -= bits
            end
            arridx += 1
        end
    end
    arr
end

# read bit packed in deprecated format (BIT_PACKED = 4)
function read_bitpacked_run_old(io, count::Int32, ::Val{W}) where {W}
    byte_width = @bit2bytewidth(W)
    typ = @byt2itype(byte_width)
    arr = Array{typ}(undef, count)
    read_bitpacked_run_old(io, count, W, arr)
end
function read_bitpacked_run_old(io, count::Int32, bits::UInt8, arr::Vector{T}, offset::Int32=Int32(0), mask::V=MASKN(bits)) where {T <: Integer, V <: Integer}
    # multiple of 8 values at a time are bit packed together
    nbytes = round(Int, (bits * count) / 8)
    #@debug("read_bitpacked_run. count:$count, nbytes:$nbytes, nbits:$bits")
    data = Array{UInt8}(undef, nbytes)
    read!(io, data)

    # the mask is of the smallest bounding type for bits
    # T is one of the types that map on to the appropriate Julia type in Parquet (which may be larger than the mask type)
    bitbuff = zero(V)
    nbitsbuff = 0

    arridx = Int32(1)
    dataidx = Int32(1)
    while arridx <= count
        diffnbits = bits - nbitsbuff
        while diffnbits > 8
            # shift 8 bits and read directly into bitbuff
            bitbuff <<= 8
            bitbuff |= data[dataidx]
            dataidx += Int32(1)
            nbitsbuff += 8
            diffnbits -= 8
        end

        if diffnbits > 0
            # read next byte from input
            nxtdata = data[dataidx]
            dataidx += Int32(1)
            # shift bitbuff by diffnbits, add diffnbits and set result
            nbitsbuff = 8 - diffnbits
            arr[arridx+offset] = convert(T, ((bitbuff << diffnbits) | (nxtdata >> nbitsbuff)) & mask)
            arridx += Int32(1)
            # keep remaining bits in bitbuff
            bitbuff <<= 8
            bitbuff |= nxtdata
        else
            # set result
            arr[arridx+offset] = convert(T, (bitbuff >> abs(diffnbits)) & mask)
            arridx += Int32(1)
            nbitsbuff -= bits
        end
    end
    arr
end

function logical_timestamp(barr; offset::Dates.Period=Dates.Second(0))
    nanos = read(IOBuffer(barr[1:8]), Int64)
    julian_days = read(IOBuffer(barr[9:12]), Int32)
    Dates.julian2datetime(julian_days) + Dates.Nanosecond(nanos) + offset
end

function logical_timestamp(i128::Int128; offset::Dates.Period=Dates.Second(0))
    iob = IOBuffer()
    write(iob, i128)
    logical_timestamp(take!(iob); offset=offset)
end

logical_string(bytes::Vector{UInt8}) = String(copy(bytes))


text:
#==========================================#
# Associative array construction and display.
#==========================================#
using D4M

# Create comma separated lists of row keys, column keys and values.
# Note: last entry in each list is the delimiter.
const r = "a,a,a,a,a,a,a,aa,aaa,b,bb,bbb,a,aa,aaa,b,bb,bbb,"
const c = "a,aa,aaa,b,bb,bbb,a,a,a,a,a,a,a,aa,aaa,b,bb,bbb,"
const v = "a-a,a-aa,a-aaa,a-b,a-bb,a-bbb,a-a,aa-a,aaa-a,b-a,bb-a,bbb-a,a-a,aa-aa,aaa-aaa,b-b,bb-bb,bbb-bbb,"

# Create assoc array and list triples.
A = Assoc(r,c,v)

# Show data structure.
dump(A);

# Display in tabular form.
printtable(A);

# Write to a CSV file.
D4M.writecsv(open("A.csv","w"), A) # or

# which is a wrapper around
open("A.csv","w") do io
    D4M.write(io, MIME"text/csv", A)
end

open("A.csv","w") do io
    D4M.write(io, MIME"text/tab-separated-values", A)
end

#==================================================#
# Show different ways to index an associative array.
#==================================================#

A = D4M.readcsv(open("A.csv", "r"))

A1r = A["a,b,",:]               # Get rows a and b.
A2r = A['a *,',1:3]             # Get rows containing a and columns 1 thru 3.
A3r = A['a,:,b,',:]             # Get rows a to b.
A4r = A[startswith"a,c,",:]     # Get rows starting with a or c.
A1c = A[:,"a,b,"]               # Get cols a and b.
A2c = A[1:3,"a *,"]             # Get rows 1 thru 3 and cols containing a.
A3c = A[:,"a,:,b,"]             # Get cols a to b.
A4c = A[:,startswith"a,c,"]     # Get cols starting with a or b.
A1v = A .< "b,"                 # Get all values less than b.

#=======================#
# Simple Math Example
#=======================#
A = D4M.readcsv("A.csv")    # Read CSV file into associative array.
A = float64(bool(A))        # Convert to numeric. Short for double(logical()).

# Show in tabular form
printtable(A);

show(sum(A,1))         # Sum into row vector.
printtable(sum(A,2))   # Sum into column vector.

# Compute a simple join.
Aa = A[:,"a,"];                  # Get column a.
Ab = A[:,"b,"];                  # Get column b.
Aab = nocol(Aa) && nocol(Ab)     # Find both a and b.

# Dataframes.jl supports complex joins, it would be nice to leverage this
# Aab = D4M.join(Aa, Ab)

#=============================================#
# Compute histogram (facets) of other columns
# that are in rows with both a and b.
#=============================================#
F = transpose( nocol(A[:,"a,"]) && nocol(A[:,"b,"]) ) * A
printtable(F.')     # Show in tabular form.

Fn = F ./ sum(A,1)  # Normalized facets.
printtable(Fn.')    # Show in tabular form.

#=============================================#
# Compute correlation.
#=============================================#
AtA = A'*A                          # this form gets special cased in julia to be fast
d = diag(adj(AtA))                  # Get diagonal of adjacency matrix.
AtA = adj(AtA) - diag(d)
setadj!(AtA, adj(AtA) - diag(d));   # Subtract diagonal, works on square matrix.
printtable(AtA)                     #  Show in tabular form.

#==================================================#
# Mixed string and numeric associative arrays.
#==================================================#

# String vectors must be row vectors
istr = "01,02,03,04,21,22,23,24,41,51,61,62,63,64,"

# Numeric vectors must be Column vectors
inum = [1,1,1,1,4,3,2,1,4,5,6,6,6,6]

A00 = Assoc("", "","")            # All empty
A01 = Assoc("", istr, inum)       # Empty row
A02 = Assoc(inum, istr, "")       # Empty value
A03 = Assoc(inum, [], "")         # Empty column and value
A04 = Assoc("", inum, [])         # Empty column and value
A05 = Assoc('','',[]);            # All empty.
A06 = Assoc('',iNum,[]);          # Row and value empty.

# Mixed type and non-empty strings
A = Assoc("a,", istr, inum)
A = Assoc(inum, istr, "a,")
A = Assoc(inum, 1, "a,")
A = Assoc("a,", 1, 1)
A = Assoc("a,", "a,", 1)
A = Assoc("a,", "a,", inum)
A = Assoc("a,", inum, 1)


text:
println("Simple example of transforming person's face from one shape to another")

using Images
using ImageView
using MAT
using VoronoiDelaunay
using PiecewiseAffineTransforms
using Colors
using FixedPointNumbers

const IMG_HEIGHT = 480
package_dir = joinpath(Pkg.dir(), "PiecewiseAffineTransforms")
data_dir = joinpath(package_dir, "examples", "cootes")

# some auxilary functions for reading data
rawdata(img) = convert(Array{Float64, 3}, data(separate(img)))
xy2ij(shape, height) = [height .- shape[:, 2] shape[:, 1]]

read_image(path) = rawdata(convert(Image{RGB{U8}}, imread(path)))
read_shape(path) = xy2ij(matread(path)["annotations"], IMG_HEIGHT)

src_img_path = joinpath(data_dir, "107_0764.bmp")
dst_img_path = joinpath(data_dir, "107_0779.bmp")
src_shape_path = src_img_path * ".mat"
dst_shape_path = dst_img_path * ".mat"

# read 2 images: src_img will be transformed to resemble dst_img
src_img = read_image(src_img_path)
# dst_img = convert(Array, separate(imread(joinpath(data_dir, dst_img_name))))
dst_img = read_image(dst_img_path)
# read corresponding shape data
src_shape = read_shape(src_shape_path)
dst_shape = read_shape(dst_shape_path)

# get triangulation
trigs = delaunayindexes(src_shape)

println("Original images and shapes (press Enter to continue)")
triplot(src_img, src_shape, trigs)
triplot(dst_img, dst_shape, trigs)
readline(STDIN)

println("Simple warp (press Enter to continue)")
@time warped_simple = pa_warp(src_img, src_shape, dst_shape, trigs)
view(src_img)
view(warped_simple)
readline(STDIN)

println("Preparing warp parameters")
@time pa_params = pa_warp_params(dst_shape, trigs, (480, 640))

println("Prepared warp (press Enter to finish)")
@time warped_prepared = pa_warp(pa_params, src_img, src_shape)
view(src_img)
view(warped_prepared)
readline(STDIN)


text:
function RDConvolution(d::Dataset;out=1024,σ=NNlib.relu)
    in_channels = prod(length ∘ values, values(d.similaritys))
    encoder = Flux.Chain(
        Conv((2, 2), in_channels => 16, σ, stride=(1, 1), pad=(1, 1)),
        MaxPool((2,2),pad=(1,1),stride=(2,2)),
        # Conv((2, 2), 16 => 16, σ, stride=(2, 2), pad=(1, 1)),
        Conv((2, 2), 16 => 32, σ, stride=(1, 1), pad=(1, 1)),
        # Conv((2, 2), 32 => 32, σ, stride=(2, 2), pad=(1, 1))
        MaxPool((2,2),pad=(1,1),stride=(2,2))
    )
    outsize = Flux.outdims(encoder, (2, d.ndiseases + d.ndrugs, in_channels,1))
    classifier = Flux.Chain(
        Flux.flatten,
        Dense(prod(outsize),out,σ),
    )
    Flux.Chain(
        encoder,
        classifier
    )
end

"注意力机制"
struct Attention
    W
    b
    H
end

""
Attention(in::Integer,out::Integer=in;σ=tanh,hid::Integer=in÷2,init=Flux.glorot_uniform)=Attention(
    init(hid,in),
    init(hid,),
    init(out,hid)
)

function (a::Attention)(x)
    softmax(a.H*tanh.(a.W*x.+a.b))
end

Flux.@functor Attention

"属性注意力机制"
struct AttrAttention
    encoders
    attention
end

""
function AttrAttention(in::Integer,out::Integer,as::Integer,init=Flux.glorot_uniform)
    AttrAttention(
        init(out,in,as),
        Attention(out,out)
    )
end

Flux.@functor AttrAttention

function (m::AttrAttention)(x)
    t,d,a,s,b=size(x)
    x1=reduce(
        vcat,
        map(eachslice(reshape(x,t*d,1,a,s*b),dims=4))do scale
            c=batched_mul(m.encoders,scale)
            e=Flux.flatten(c)
            sum(m.attention(e).*e,dims=2)
        end
    )
    reshape(x1,:,s,b)
end

struct ScaleAttention
    W
    b
    h
end

ScaleAttention(in::Integer;init=Flux.glorot_uniform)=ScaleAttention(init(in÷2,in),init(in÷2),init(1,in÷2))

function (m::ScaleAttention)(x)
    d,s,b=size(x)
    α=reduce(
        vcat,
        map(eachslice(x,dims=3))do t
            e=softmax(m.h*tanh.(m.W*t.+m.b),dims=2)
            vec(e)
        end
    )
    reshape(α,1,s,b).*x
end

Flux.@functor ScaleAttention

struct MSRD
    attratt
    scaleatt
    forward
    backward
    encoder
end

function MSRD(data::Dataset;out_features=512,hidden_features=out_features÷2)
    in_features=2*(data.ndiseases+data.ndrugs)
    as=prod(length ∘ values, values(data.similaritys))
    MSRD(
        AttrAttention(in_features,hidden_features,as),
        ScaleAttention(hidden_features),
        LSTM(hidden_features,out_features),
        LSTM(hidden_features,out_features),
        Dense(2out_features,out_features,sigmoid)
    )
end

Flux.@functor MSRD

function (m::MSRD)(x)
    x1=m.attratt(x)
    x2=m.scaleatt(x1)
    x3=reduce(
        hcat,
            map(eachslice(x2,dims=3))do t
            Flux.reset!(m.forward)
            for e in eachslice(t,dims=2)
                m.forward(e)
            end
            h,_=m.forward.state
            h
        end
    )
    x4=reduce(
        hcat,
        map(eachslice(x2,dims=3))do t
            Flux.reset!(m.backward)
            for e in eachslice(reverse(t,dims=2),dims=2)
                m.backward(e)
            end
            h,_=m.backward.state
            h
        end
    )
    x5=m.encoder(cat(x3,x4,dims=1))
end

struct CombineModel
    left
    right
    combine
end

function CombineModel(data::Dataset;feature_left=512,feature_right=512)
    CombineModel(
        MSRD(data;out_features=feature_left),
        RDConvolution(data;out=feature_right),
        Dense(feature_left+feature_right,2)
    )
end

function (m::CombineModel)(xₗ,xᵣ)
    yₗ=m.left(xₗ)
    yᵣ=m.right(xᵣ)
    y=m.combine(cat(yₗ,yᵣ,dims=1))
end

Flux.@functor CombineModel

export RDConvolution,Attention,AttrAttention,ScaleAttention,MSRD,CombineModel

max_stars_repo_name: yeruoforever/DrugRepositioning

max_issues_repo_name: yeruoforever/DrugRepositioning

max_forks_repo_name: yeruoforever/DrugRepositioning