text:
#function for ArgParse

"""
    test_parse_main(ARGS::Vector{String})
Defines argument parsing rules for viva script.
"""
function test_parse_main(ARGS::Vector{String})

    # initialize the settings
    s = ArgParseSettings(

    # text in the description is printed to screen when help is called with help flag (--help, -h)
    description = "VIVA VCF Visualization Tool is a tool for creating publication quality plots of data contained within VCF files. For a complete description of features with examples read the docs here https://github.com/compbiocore/VariantVisualization.jl",
    suppress_warnings = true,
    epilog = "Thank you for using VIVA. Please submit any bugs to https://github.com/compbiocore/VariantVisualization.jl/issues "
    )

    @add_arg_table s begin
        "--vcf_file", "-f"
        help = "vcf filename in format: file.vcf"
        arg_type = String
        required = true

        "--output_directory", "-o"
        help =" function checks if directory exists and saves there, if not creates and saves here"
        arg_type = String
        default = "output"

        "--save_format", "-s"
        help = "file format you wish to save graphics as (eg. pdf, html, png). Defaults to html"
        arg_type = String
        default = "html"

        "--genomic_range", "-r"
        help = "select rows within a given chromosome range. Provide chromosome range after this flag in format chr4:20000000-30000000."
        arg_type = String

        "--pass_filter", "-p"
        help = "select rows with PASS in the FILTER field."
        action = :store_true

        "--positions_list", "-l"
        help = "select variants matching list of chromosomal positions. Provide filename of text file formatted with two columns in csv format: 1,2000345."
        arg_type = String

        "--group_samples", "-g"
        help = "group samples by common trait using user generated matrix key of traits and sample names following format guidelines in documentation. Provide file name of .csv file"
        nargs = 2
        arg_type = String

        "--select_samples"
        help = "select samples to include in visualization by providing tab delimited list of sample names (eg. samplenames.txt). Works for heatmap visualizations and numeric array generation only (not average dp plots)"
        arg_type = String

        "--heatmap", "-m"
        help = "genotype field to visualize (eg. genotype, read_depth, or 'genotype,read_depth' to visualize each separately)"
        arg_type = String
        default = "genotype,read_depth"

        "--y_axis_labels", "-y"
        help = "specify whether to label y-axis with all chromosome positions (options = positions / chromosome) separators. Defaults to chromosome separators."
        default = "chromosomes"
        arg_type = String

        "--x_axis_labels", "-x"
        help = "flag to specify whether to label x-axis with sample ids from vcf file. Defaults to FALSE."
        action = :store_true

        "--num_array", "-n"
        help = "flag to save numeric array of categorical genotype values or read depth values before heatmap plotting. Must be used with --heatmap set."
        action = :store_true

        "--heatmap_title", "-t"
        help = "Specify filename for heatmap with underscores for spaces."
        arg_type = String

        "--avg_dp"
        help = "visualize average read depths as line chart. Options: average sample read depth, average variant read depth, or both. eg. =sample, =variant, =sample,variant"

        "--save_remotely"
        help = "Save html support files online rather than locally so files can be shared between systems. Files saved in this way require internet access to open."
        action = :store_true

    end

    parsed_args = parse_args(s)

    #Developers can display parsed args during new feature development:
    #println("Parsed args:")

#activate block to show all argument keys

#=
    for (key,val) in parsed_args
        println("  $key  =>  $(repr(val))")
    end
    =#

    return parsed_args
end

#functions for loading VCF file, vcf stats, clean vcf, load siglist

"""
    clean_column1!(matrix_with_chr_column)
Replace String "X","Y","M" from chromosome column with 23,24,25 respectively so variants can be sorted by descending chr position for plotting
"""
function clean_column1!(matrix_with_chr_column)

    for i = 1:(size(matrix_with_chr_column, 1))
        matrix_with_chr_column[i, 1] = matrix_with_chr_column[i, 1] == "X" ? "23" : matrix_with_chr_column[i, 1]
        matrix_with_chr_column[i, 1] = matrix_with_chr_column[i, 1] == "Y" ? "24" : matrix_with_chr_column[i, 1]
        matrix_with_chr_column[i, 1] = matrix_with_chr_column[i, 1] == "M" ? "25" : matrix_with_chr_column[i, 1]
        matrix_with_chr_column[i, 1] = matrix_with_chr_column[i, 1] == "MT" ? "25" : matrix_with_chr_column[i, 1]
        matrix_with_chr_column[i, 1] = matrix_with_chr_column[i, 1] == "Mt" ? "25" : matrix_with_chr_column[i, 1]
    end
end

"""
    clean_column1_siglist!(siglist)
Replaces strings "X","Y","M" with 23,24,25 {Int} in array generated in load_siglist()
use in load_siglist() because X and Y need to be replaced with Int
"""
function clean_column1_siglist!(siglist)

    for i = 1:(size(siglist, 1))
        siglist[i, 1] = siglist[i, 1] == "X" ? 23 : siglist[i, 1]
        siglist[i, 1] = siglist[i, 1] == "Y" ? 24 : siglist[i, 1]
        siglist[i, 1] = siglist[i, 1] == "M" ? 25 : siglist[i, 1]
        siglist[i, 1] = siglist[i, 1] == "MT" ? 25 : siglist[i, 1]
        siglist[i, 1] = siglist[i, 1] == "Mt" ? 25 : siglist[i, 1]
    end
end

"""
    returnXY_column1!(chr_label_vector)
Replace String "23","24","25" with "X","Y","M" in chromosome label vector used for plot labels
"""
function returnXY_column1!(chr_label_vector)

    for i = 1:(size(chr_label_vector, 1))
        chr_label_vector[i, 1] = chr_label_vector[i, 1] == 23 ? "X" : chr_label_vector[i, 1]
        chr_label_vector[i, 1] = chr_label_vector[i, 1] == 24 ? "Y" : chr_label_vector[i, 1]
        chr_label_vector[i, 1] = chr_label_vector[i, 1] == 25 ? "M" : chr_label_vector[i, 1]
        chr_label_vector[i, 1] = chr_label_vector[i, 1] == 25 ? "MT" : chr_label_vector[i, 1]
        chr_label_vector[i, 1] = chr_label_vector[i, 1] == 25 ? "Mt" : chr_label_vector[i, 1]
    end
end

"""
    returnXY_column1_siglist!(siglist_sorted)
Replace String "23","24","25" with "X","Y","M" in siglist for filtering
"""
function returnXY_column1_siglist!(siglist_sorted)

    for i = 1:(size(siglist_sorted, 1))
        siglist_sorted[i, 1] = siglist_sorted[i, 1] == 23 ? "X" : siglist_sorted[i, 1]
        siglist_sorted[i, 1] = siglist_sorted[i, 1] == 24 ? "Y" : siglist_sorted[i, 1]
        siglist_sorted[i, 1] = siglist_sorted[i, 1] == 25 ? "M" : siglist_sorted[i, 1]
        siglist_sorted[i, 1] = siglist_sorted[i, 1] == 25 ? "MT" : siglist_sorted[i, 1]
        siglist_sorted[i, 1] = siglist_sorted[i, 1] == 25 ? "Mt" : siglist_sorted[i, 1]
    end
end

"""
    load_siglist(filename::AbstractString)
where x = filename of significant SNP variant location list in comma delimited format (saved as .csv)
"""
function load_siglist(filename::AbstractString)

siglist_unsorted = readdlm(filename, ',', skipstart=1)

VariantVisualization.clean_column1_siglist!(siglist_unsorted)

siglist = sortslices(siglist_unsorted, dims=1, by=x->(x[1],x[2]))

returnXY_column1_siglist!(siglist)

if typeof(siglist) == Array{Float64,2}
    siglist = trunc.(Int,siglist)
    return siglist
end

return siglist

end

"""
    build_set_from_list(sig_list)
build set of tuples of chrom and pos of each record in vcf for use in sig_list_filters.
Method 1: build list from input variant locations in Int64 and String chromosomes in variant list (chr1-XYM)
"""
function build_set_from_list(sig_list) #did not define type because throws errors. Tried defining sig_list::Array(Any,2), ::Array(Int,2), and ::Array(String,2) and only stored one
    position_set = Set{Tuple{Any,Int}}()

    for row in 1:size(sig_list,1)
        chr = sig_list[row,1]
        pos = sig_list[row,2]
        push!(position_set, (string(chr), Int(pos)))
    end

    return position_set
end

#functions for variant filters

"""
io_genomic_range_vcf_filter(chr_range::String, vcf_filename::AbstractString)
create subarray of vcf variant records matching user specified chromosome range in format: (e.g. chr1:0-30000000)
"""
function io_genomic_range_vcf_filter(chr_range::String,vcf_filename::AbstractString)

        a=split(chr_range,":")
        chrwhole=a[1]
        chrnumber=split(chrwhole,"r")
        string_chr=chrnumber[2]
        chr=String(string_chr)
        range=a[2]
        splitrange=split(range, "-")
        lower_limit=splitrange[1]
        chr_range_low=Meta.parse(lower_limit)
        upper_limit=splitrange[2]
        chr_range_high=Meta.parse(upper_limit)

        reader = VCF.Reader(open(vcf_filename, "r"))
        record1=VCF.read(reader)

        open(VCF.Reader, vcf_filename) do reader

        vcf_record = VCF.Record()
        vcf_subarray = Array{Any}(undef,0)

        if occursin("chr",VCF.chrom(record1))

            while !eof(reader)
                read!(reader, vcf_record)

                   if (VCF.chrom(vcf_record) == chrwhole) && (chr_range_high >= VCF.pos(vcf_record) >= chr_range_low)
                          push!(vcf_subarray,copy(vcf_record))
                   end
            end

        else

            while !eof(reader)
                read!(reader, vcf_record)

                   if (VCF.chrom(vcf_record) == chr) && (chr_range_high >= VCF.pos(vcf_record) >= chr_range_low)
                          push!(vcf_subarray,copy(vcf_record))
                   end
            end
        end

        return vcf_subarray

    end

end

"""
    io_sig_list_vcf_filter(sig_list,vcf_filename)
returns subarray of variant records matching a list of variant positions returned from load_siglist()
"""
function io_sig_list_vcf_filter(sig_list, vcf_filename)
        open(VCF.Reader, vcf_filename) do reader

        vcf_subarray = Array{Any}(undef,0)

        posset = build_set_from_list(sig_list)
        vcf_record = VCF.Record()

        while !eof(reader)
            read!(reader, vcf_record)
            chr = VCF.chrom(vcf_record)
            pos = VCF.pos(vcf_record)

            if (chr, pos) in posset
                push!(vcf_subarray, copy(vcf_record))
            end
        end

        return vcf_subarray
    end
end

"""
    io_pass_filter(vcf_filename)
returns subarray of vcf records including only records with FILTER status = PASS
"""
function io_pass_filter(vcf_filename)

    open(VCF.Reader, vcf_filename) do reader

    vcf_subarray = Array{Any}(undef,0)

    vcf_record = VCF.Record()

    while !eof(reader)
        read!(reader, vcf_record)

       if VCF.hasfilter(vcf_record) && VCF.filter(vcf_record) == String["PASS"]
              push!(vcf_subarray,copy(vcf_record))
       end
end

return vcf_subarray
end
end

"""
    pass_genomic_range_siglist_filter(vcf_filename,sig_list,chr_range::AbstractString)
returns subarray of vcf records with io_pass_filter, io_sig_list_vcf_filter, and io_genomic_range_vcf_filter applied.
"""
function pass_genomic_range_siglist_filter(vcf_filename,sig_list,chr_range::AbstractString)

    a=split(chr_range,":")
    chrwhole=a[1]
    chrnumber=split(chrwhole,"r")
    string_chr=chrnumber[2]
    chr=String(string_chr)
    range=a[2]
    splitrange=split(range, "-")
    lower_limit=splitrange[1]
    chr_range_low=Meta.parse(lower_limit)
    upper_limit=splitrange[2]
    chr_range_high=Meta.parse(upper_limit)

    posset = build_set_from_list(sig_list)

    vcf_subarray = Array{Any}(undef,0)
    vcf_record = VCF.Record()

    reader = VCF.Reader(open(vcf_filename, "r"))
    record1=VCF.read(reader)

    open(VCF.Reader, vcf_filename) do reader

    while !eof(reader)
        read!(reader, vcf_record)

        if size(vcf_subarray,1) == size(sig_list,1)
                       break
        end

    for row = 1:size(sig_list,1)
           dimension = size(sig_list,1)

           chr_sig=(sig_list[row,1])
           pos_sig=(sig_list[row,2])

                  if typeof(VCF.chrom(vcf_record)) == String

                      chr_sig = string(chr_sig)

                         if occursin("chr",VCF.chrom(record1))

                             if (VCF.chrom(vcf_record) == chr_sig) && (VCF.pos(vcf_record) == pos_sig) && (VCF.hasfilter(vcf_record)) && (VCF.filter(vcf_record) == String["PASS"]) && ((VCF.chrom(vcf_record) == chrwhole)) && ((chr_range_high >= VCF.pos(vcf_record) >= chr_range_low))
                                 push!(vcf_subarray,copy(vcf_record))
                                 break
                             end

                        else

                            if (VCF.chrom(vcf_record) == chr_sig) && (VCF.pos(vcf_record) == pos_sig) && (VCF.hasfilter(vcf_record)) && (VCF.filter(vcf_record) == String["PASS"]) && ((VCF.chrom(vcf_record) == chr)) && ((chr_range_high >= VCF.pos(vcf_record) >= chr_range_low))
                                push!(vcf_subarray,copy(vcf_record))
                                break
                            end
                        end

                 else

                     if occursin("chr",VCF.chrom(record1))

                         if (VCF.chrom(vcf_record) == chr_sig) && (VCF.pos(vcf_record) == pos_sig) && (VCF.hasfilter(vcf_record)) && (VCF.filter(vcf_record) == String["PASS"]) && ((VCF.chrom(vcf_record) == chrwhole)) && ((chr_range_high >= VCF.pos(vcf_record) >= chr_range_low))
                             push!(vcf_subarray,copy(vcf_record))
                             break
                         end

                     else
                        if (VCF.chrom(vcf_record) == chr_sig) && (VCF.pos(vcf_record) == pos_sig) && (VCF.hasfilter(vcf_record)) && (VCF.filter(vcf_record) == String["PASS"]) && ((VCF.chrom(vcf_record) == chr)) && ((chr_range_high >= VCF.pos(vcf_record) >= chr_range_low))
                            push!(vcf_subarray,copy(vcf_record))
                            break
                        end
                     end

                 end
             end
           end
    end

    return vcf_subarray
    end

"""
    pass_genomic_range_filter(reader::GeneticVariation.VCF.Reader,chr_range::AbstractString,vcf_filename)
returns subarray of vcf records with io_pass_filter and io_genomic_range_vcf_filter applied.
"""
    function pass_genomic_range_filter(reader::GeneticVariation.VCF.Reader,chr_range::AbstractString,vcf_filename)

        a=split(chr_range,":")
        chrwhole=a[1]
        chrnumber=split(chrwhole,"r")
        string_chr=chrnumber[2]
        chr=String(string_chr)
        range=a[2]
        splitrange=split(range, "-")
        lower_limit=splitrange[1]
        chr_range_low=Meta.parse(lower_limit)
        upper_limit=splitrange[2]
        chr_range_high=Meta.parse(upper_limit)

        vcf_subarray = Array{Any}(undef,0)

        reader = VCF.Reader(open(vcf_filename, "r"))
        record1=VCF.read(reader)

        vcf_record = VCF.Record()

        open(VCF.Reader, vcf_filename) do reader

            while !eof(reader)
                read!(reader, vcf_record)

                 if occursin("chr",VCF.chrom(record1))

                      if typeof(VCF.chrom(vcf_record)) == String

                             chr = string(chr)

                             if (VCF.hasfilter(vcf_record)) && (VCF.filter(vcf_record) == String["PASS"]) && ((VCF.chrom(vcf_record) == chrwhole)) && ((chr_range_high >= VCF.pos(vcf_record) >= chr_range_low))
                                 push!(vcf_subarray,copy(vcf_record))

                             end

                     else

                             if (VCF.hasfilter(vcf_record)) && (VCF.filter(vcf_record) == String["PASS"]) && ((VCF.chrom(vcf_record) == chrwhole)) && ((chr_range_high >= VCF.pos(vcf_record) >= chr_range_low))
                                 push!(vcf_subarray,copy(vcf_record))
                             end
                     end

                 else

                     if typeof(VCF.chrom(vcf_record)) == String
                            chr = string(chr)

                            if (VCF.hasfilter(vcf_record)) && (VCF.filter(vcf_record) == String["PASS"]) && ((VCF.chrom(vcf_record) == chr)) && ((chr_range_high >= VCF.pos(vcf_record) >= chr_range_low))
                                push!(vcf_subarray,copy(vcf_record))
                            end

                    else

                            if (VCF.hasfilter(vcf_record)) && (VCF.filter(vcf_record) == String["PASS"]) && ((VCF.chrom(vcf_record) == chr)) && ((chr_range_high >= VCF.pos(vcf_record) >= chr_range_low))
                                push!(vcf_subarray,copy(vcf_record))

                            end
                    end

               end
           end
       end

        return vcf_subarray
    end

"""
    pass_siglist_filter(vcf_filename,sig_list,chr_range::AbstractString)
returns subarray of vcf records with io_pass_filter, io_sig_list_vcf_filter, and io_genomic_range_vcf_filter applied.
"""
function pass_siglist_filter(vcf_filename,sig_list)

    vcf_subarray = Array{Any}(undef,0)
    open(VCF.Reader, vcf_filename) do reader

    posset = build_set_from_list(sig_list)
    vcf_record = VCF.Record()

    while !eof(reader)
        read!(reader, vcf_record)

        if size(vcf_subarray,1) == size(sig_list,1)
                       break
        end

        for row= 1:size(sig_list,1)
               dimension = size(sig_list,1)

               chr=(sig_list[row,1])
               pos=(sig_list[row,2])

                      if typeof(VCF.chrom(vcf_record)) == String
                             chr = string(chr)

                             if (VCF.chrom(vcf_record) == chr) && (VCF.pos(vcf_record) == pos) && (VCF.hasfilter(vcf_record)) && (VCF.filter(vcf_record) == String["PASS"])
                                push!(vcf_subarray, copy(vcf_record))
                                 break
                             end

                     else

                             if (VCF.chrom(vcf_record) == chr) && (VCF.pos(vcf_record) == pos) && (VCF.hasfilter(vcf_record)) && (VCF.filter(vcf_record) == String["PASS"])
                                push!(vcf_subarray, copy(vcf_record))
                                 break
                             end
                     end
        end
    end
end

    return vcf_subarray
end

"""
    genomic_range_siglist_filter(vcf_filename,sig_list,chr_range::AbstractString)
returns subarray of vcf records with io_pass_filter, io_sig_list_vcf_filter, and io_genomic_range_vcf_filter applied.
"""
function genomic_range_siglist_filter(vcf_filename,sig_list,chr_range::AbstractString)

    a=split(chr_range,":")
    chrwhole=a[1]
    chrnumber=split(chrwhole,"r")
    string_chr=chrnumber[2]
    chr=String(string_chr)
    range=a[2]
    splitrange=split(range, "-")
    lower_limit=splitrange[1]
    chr_range_low=Meta.parse(lower_limit)
    upper_limit=splitrange[2]
    chr_range_high=Meta.parse(upper_limit)

    vcf_subarray = Array{Any}(undef,0)

    reader = VCF.Reader(open(vcf_filename, "r"))
    record1=VCF.read(reader)

    open(VCF.Reader, vcf_filename) do reader

    vcf_record = VCF.Record()

    while !eof(reader)
        read!(reader, vcf_record)

        if size(vcf_subarray,1) == size(sig_list,1)
                       break
        end

        for row= 1:size(sig_list,1)
            dimension = size(sig_list,1)

            chr_sig=(sig_list[row,1])
            pos_sig=(sig_list[row,2])

            if occursin("chr",VCF.chrom(record1))

                    if typeof(VCF.chrom(vcf_record)) == String
                             chr_sig = string(chr_sig)

                             if (VCF.chrom(vcf_record) == chr_sig) && (VCF.pos(vcf_record) == pos_sig) && ((VCF.chrom(vcf_record) == chrwhole)) && ((chr_range_high >= VCF.pos(vcf_record) >= chr_range_low))
                                 push!(vcf_subarray,copy(vcf_record))
                                 break
                             end

                     else

                             if (VCF.chrom(vcf_record) == chr_sig) && (VCF.pos(vcf_record) == pos_sig) && ((VCF.chrom(vcf_record) == chrwhole)) && ((chr_range_high >= VCF.pos(vcf_record) >= chr_range_low))
                                 push!(vcf_subarray,copy(vcf_record))
                                 break
                             end
                     end

            else

                     if typeof(VCF.chrom(vcf_record)) == String
                              chr_sig = string(chr_sig)

                              if (VCF.chrom(vcf_record) == chr_sig) && (VCF.pos(vcf_record) == pos_sig) && ((VCF.chrom(vcf_record) == chr)) && ((chr_range_high >= VCF.pos(vcf_record) >= chr_range_low))
                                  push!(vcf_subarray,copy(vcf_record))
                                  break
                              end

                      else

                              if (VCF.chrom(vcf_record) == chr_sig) && (VCF.pos(vcf_record) == pos_sig) && ((VCF.chrom(vcf_record) == chr)) && ((chr_range_high >= VCF.pos(vcf_record) >= chr_range_low))
                                  push!(vcf_subarray,copy(vcf_record))
                                  break
                              end
                      end
            end
        end
    end
    end

    return vcf_subarray
end

#functions for converting vcf record array to numerical array, manipulating numerical array for plotting, and generating data structures for plotting

"""
    create_chr_dict()
creates dict for use in combined_all_genotype_array_functions() for removing 'chr' from chromosome labels to allow sorting variant records by chromosome position.
"""
function create_chr_dict()

    chr_dict = Dict()

    chr_dict = Dict("chr1"=> 1,"chr2"=> 2,"chr3"=> 3,"chr4"=> 4,"chr5"=> 5,"chr6"=> 6,"chr7"=> 7,"chr8"=> 8,"chr9"=> 9,"chr10"=> 10,"chr11"=> 11,"chr12"=> 12,"chr13"=> 13,"chr14"=> 14,"chr15"=> 15,"chr16"=> 16,"chr17"=> 17,"chr18"=> 18,"chr19"=> 19,"chr20"=> 20,"chr21"=> 21,"chr22"=> 22,"chrX"=> X,"chrY"=> Y,"chrM"=> M)

    return chr_dict
end

"""
    combined_all_genotype_array_functions(sub)
convert sub from variant filters to gt_num_array and gt_chromosome_labels for plot functions.
"""
function combined_all_genotype_array_functions(sub)
    genotype_array = generate_genotype_array(sub,"GT")
    map!(s->replace(s, "chr" => ""), genotype_array, genotype_array)
    clean_column1!(genotype_array)
    genotype_array=VariantVisualization.sort_genotype_array(genotype_array)
    geno_dict = define_geno_dict()
    gt_num_array,gt_chromosome_labels = translate_genotype_to_num_array(genotype_array, geno_dict)

    return gt_num_array,gt_chromosome_labels
end

"""
    combined_all_read_depth_array_functions(sub)
convert sub from variant filters to dp_num_array and dp_chromosome_labels for plot functions.
"""
function combined_all_read_depth_array_functions(sub)

    read_depth_array = VariantVisualization.generate_genotype_array(sub,"DP")
    map!(s->replace(s, "chr"=>""), read_depth_array, read_depth_array)
    clean_column1!(read_depth_array)
    read_depth_array=VariantVisualization.sort_genotype_array(read_depth_array)
    dp_num_array,dp_chromosome_labels = translate_readdepth_strings_to_num_array(read_depth_array)
    return dp_num_array,dp_chromosome_labels
end

"""
    combined_all_read_depth_array_functions_for_avg_dp(sub)
convert sub from variant filters to dp_num_array and dp_chromosome_labels for plot functions.
"""
function combined_all_read_depth_array_functions_for_avg_dp(sub)

    read_depth_array = VariantVisualization.generate_genotype_array(sub,"DP")
    map!(s->replace(s, "chr"=>""), read_depth_array, read_depth_array)
    clean_column1!(read_depth_array)
    read_depth_array=VariantVisualization.sort_genotype_array(read_depth_array)
    dp_num_array,dp_chromosome_labels = translate_readdepth_strings_to_num_array_for_avg_dp(read_depth_array)
    return dp_num_array,dp_chromosome_labels
end

"""
    generate_genotype_array(record_sub::Array{Any,1},genotype_field::String)
Returns numerical array of genotype values (either genotype or read_depth values) which are translated by another function into num_array
Where genotype_field is either GT or DP to visualize genotype or read_depth
"""
function generate_genotype_array(record_sub::Array{Any,1},y)

       num_samples = length(VCF.genotype(record_sub[1]))

       df = DataFrame(String, 0, num_samples+2)

       for record in record_sub

              genotype_data_per_variant = VCF.genotype(record, 1:num_samples, y)
              vcf_chrom = VCF.chrom(record)
              vcf_pos = string(VCF.pos(record))
              genotype_data_per_variant = vcat(vcf_pos,genotype_data_per_variant)
              genotype_data_per_variant = vcat(vcf_chrom,genotype_data_per_variant)
              push!(df, genotype_data_per_variant)
       end

       pre_num_array = Matrix(df)

       return pre_num_array
end

"""
    sort_genotype_array(genotype_array)
sorts genotype data array (which can contain genotype or read depth values) for GT or DP by chromosomal location
"""
function sort_genotype_array(genotype_array)

    data=genotype_array[:,3:size(genotype_array,2)]
    chrom_positions = [parse(Int, i) for i in genotype_array[:,1:2]]
    genotype_array = hcat(chrom_positions,data)
    genotype_array = sortslices(genotype_array, dims=1, by=x->(x[1],x[2]))

return genotype_array
end

"""
    define_geno_dict()
returns dictionary of values for use in replace_genotype_with_vals()
"""
function define_geno_dict()
    geno_dict = Dict()

    homo_variant = ["1/1" "1/2" "2/2" "1/3" "2/3" "3/3" "1/4" "2/4" "3/4" "4/4" "1/5" "2/5" "3/5" "4/5" "5/5" "1/6" "2/6" "3/6" "4/6" "5/6" "6/6" "1|1" "1|2" "2|2" "1|3" "2|3" "3|3" "1|4" "2|4" "3|4" "4|4" "1|5" "2|5" "3|5" "4|5" "5|5" "1|6" "2|6" "3|6" "4|6" "5|6" "6|6" "2/1" "3/1" "4/1" "5/1" "6/1" "3/2" "4/2" "5/2" "6/2" "4/3" "5/3" "6/3" "5/4" "6/4" "6/5" "2|1" "3|1" "4|1" "5|1" "6|1" "3|2" "4|2" "5|2" "6|2" "4|3" "5|3" "6|3" "5|4" "6|4" "6|5"]

    hetero_variant = ["0/1" "0/2" "0/3" "0/4" "0/5" "0/6" "1/0" "2/0" "3/0" "4/0" "5/0" "6/0" "0|1" "0|2" "0|3" "0|4" "0|5" "0|6" "1|0" "2|0" "3|0" "4|0" "5|0" "6|0" "1/0" "2/0" "3/0" "4/0" "5/0" "6/0" "1|0" "2|0" "3|0" "4|0" "5|0" "6|0"]

    no_data = ["./." ".|." "."]

    ref = ["0/0" "0|0"]

    for item in homo_variant
           geno_dict[item] = 3
    end
    for item in hetero_variant
           geno_dict[item] = 2
    end
    for item in no_data
           geno_dict[item] = 0
    end
    for item in ref
           geno_dict[item] = 1
    end

    return geno_dict
end

"""
    translate_genotype_to_num_array(genotype_array,geno_dict)
returns a tuple of num_array for plotting, and chromosome labels for plotting as label bar.
Translates array of genotype values to numerical array of categorical values.
Genotype values are converted to categorical values. No_call=0, 0/0=1, heterozygous_variant=2, homozygous_variant=3
"""
function translate_genotype_to_num_array(genotype_array,geno_dict)

    array_for_plotly = genotype_array[:,3:size(genotype_array,2)]
    chromosome_labels = genotype_array[:,1:2]

        function translate(c)
               geno_dict[c]
        end

    num_array = map(translate, array_for_plotly)

    return num_array,chromosome_labels
end

"""
    translate_readdepth_strings_to_num_array(read_depth_array::Array{Any,2})
Returns array of read_depth as int for plotting and average calculation.
Ceiling of dp=100 is set to prevent high dp value from hiding (or "blowing out") low dp values. (see read_depth_threshhold() ).
Where read_depth_array is output of generate_genotype_array() for DP option
returns a tuple of num_array type Int for average calculation and plotting, and chromosome labels for plotting as label bar
"""
function translate_readdepth_strings_to_num_array(read_depth_array::Array{Any,2})

       chromosome_labels = read_depth_array[:,1:2]

       dp_array_for_plotly = read_depth_array[:,3:size(read_depth_array,2)]

       map!(s->replace(s, "."=>"-20"), dp_array_for_plotly, dp_array_for_plotly)

       dp_array_for_plotly = [parse(Int, i) for i in dp_array_for_plotly]

       return dp_array_for_plotly, chromosome_labels
end

"""
    translate_readdepth_strings_to_num_array_for_avg_dp(read_depth_array::Array{Any,2})
Returns array of read_depth as int for plotting and average calculation.
'read_depth_array' is output of generate_genotype_array() for DP option
returns a tuple of num_array type Int for average calculation and plotting, and chromosome labels for plotting as label bar
No call is replaced with 0 for avg_calculation. Ceiling of dp=100 is set to prevent high dp value from hiding (or "blowing out") low dp values.
"""
function translate_readdepth_strings_to_num_array_for_avg_dp(read_depth_array::Array{Any,2})

       chromosome_labels = read_depth_array[:,1:2]

       dp_array_for_plotly = read_depth_array[:,3:size(read_depth_array,2)]

       map!(s->replace(s, "."=>"0"), dp_array_for_plotly, dp_array_for_plotly)

       dp_array_for_plotly = [parse(Int, i) for i in dp_array_for_plotly]

       return dp_array_for_plotly, chromosome_labels
end

#functions for sample filters

"""
    get_sample_names(reader)
returns sample ids of vcf file as a vector of symbols for naming columns of num_array dataframe object for column filter functions
"""
function get_sample_names(reader)

    hdr=VCF.header(reader)
    sample_names=hdr.sampleID

    sample_names=reshape(sample_names,1,length(sample_names))

    sample_names = [Symbol(i) for i in sample_names]

    return sample_names
end

"""
    find_group_label_indices(pheno,trait_to_group_by,row_to_sort_by)
find indices and determines names for group 1 and group 2 labels on plots. finds index of center of each sample group to place tick mark and label.
"""
function find_group_label_indices(pheno,trait_to_group_by,row_to_sort_by)

    group1_label=split(trait_to_group_by,",")[1]
    group2_label=split(trait_to_group_by,",")[2]

    group1_index = (something(findlast(isequal(1),pheno[row_to_sort_by,:]),0) -
    something(findfirst(isequal(1),pheno[row_to_sort_by,:]),0))/2

    group2_index = (something(findlast(isequal(2),pheno[row_to_sort_by,:]),0) -
    something(findfirst(isequal(2),pheno[row_to_sort_by,:]),0))/2 +
    something(findlast(isequal(1),pheno[row_to_sort_by,:]),0) + 0.5

    group_dividing_line = something(findfirst(isequal(2),pheno[row_to_sort_by,:]),0) - 0.5

    group_label_pack = [group1_index, group2_index, group_dividing_line, group1_label, group2_label]

    return group_label_pack
end

"""
    sortcols_by_phenotype_matrix(pheno_matrix_filename::String,trait_to_group_by::String,num_array::Array{Int64,2}, sample_names::Array{Symbol,2})
group samples by a common trait using a user generated key matrix ("phenotype matrix") returns num_array,group_label_pack,
"""
function sortcols_by_phenotype_matrix(pheno_matrix_filename::String,trait_to_group_by::String, num_array::Array{Int64,2}, sample_names)

    pheno = readdlm(pheno_matrix_filename, ',')

    #get row numer of user-chosen phenotype characteristic to sort columns by
    row_to_sort_by = findall(x -> x == trait_to_group_by, pheno)
    row_to_sort_by = row_to_sort_by[1]

    #remove phenotype_row_labels used to identify row to sort by, so row can be sorted without strings causing issues
    pheno_no_trait_labels = pheno[:,2:size(pheno,2)]
    trait_labels=pheno[2:size(pheno,1),1]

    pheno_no_trait_labels = sortslices(pheno_no_trait_labels, dims=2, by = x -> x[row_to_sort_by], rev = false)

    id_list = pheno_no_trait_labels[1,:]

    sample_ids=vec(sample_names)

    col_new_order=vec(id_list)

    col_new_order = [Symbol(i) for i in col_new_order]

    df1_vcf = DataFrame(num_array)

    names!(df1_vcf, sample_ids)

    ordered_num_array = df1_vcf[:, col_new_order]

    ordered_num_array = Matrix(ordered_num_array)

    group_label_pack = find_group_label_indices(pheno_no_trait_labels,trait_to_group_by,row_to_sort_by)

    id_list=reshape(id_list,1,length(id_list))

    return ordered_num_array,group_label_pack,pheno_no_trait_labels,id_list,trait_labels

end

"""
    select_columns(filename_sample_list::AbstractString, num_array::Array{Int64,2}, sample_names)
returns num_array with columns matching user generated list of sample ids to select for analysis. num_array now has sample ids in first row.
"""
function select_columns(filename_sample_list::AbstractString, num_array::Array{Int64,2}, sample_names)

    selectedcolumns=readdlm(filename_sample_list)

    header_as_strings = selectedcolumns

    for item = 1:size(selectedcolumns,2)
        selectedcolumns[item] = Symbol(selectedcolumns[item])
    end

    col_selectedcolumns=vec(selectedcolumns)

    selectedcolumns=Symbol(col_selectedcolumns)

    df_num_array = DataFrame(num_array)

    rename!(df_num_array, f => t for (f, t) = zip(names(df_num_array), sample_names))

    df_selected_samples_num_array = df_num_array[:, col_selectedcolumns]

    selected_samples_num_array = Matrix(df_selected_samples_num_array)

    sample_names=permutedims(col_selectedcolumns)

    return selected_samples_num_array,sample_names
end

#functions for mathematic analysis
"""
    avg_dp_samples(dp_num_array::Array{Int64,2})
create sample_avg_list vector that lists averages of read depth for each sample for input into avg_sample_dp_line_chart(sample_avg_list)
dp_num_array must contain dp values as Int64 and be without chromosome position columns
"""
function avg_dp_samples(dp_num_array::Array{Int64,2})

    avg_dps_all = Array{Float64}(undef,0)

    for column = 1:size(dp_num_array,2)

        all_dps_patient = dp_num_array[1:size(dp_num_array,1),column]
        sum_dp = sum(all_dps_patient)
        avg_dp = sum_dp/(size(dp_num_array,1))
        push!(avg_dps_all,avg_dp)
    end

    return avg_dps_all
end

"""
    avg_dp_variant(dp_num_array::Array{Int64,2})
create variant_avg_list vector that lists averages of read depth for each variant for input into avg_variant_dp_line_chart(variant_avg_list)
"""
function avg_dp_variant(dp_num_array::Array{Int64,2})

    avg_dps_all = Array{Float64}(undef,0)

    for row = 1:size(dp_num_array,1)

        all_dps_variant = dp_num_array[row,1:size(dp_num_array,2)]
        sum_dp = sum(all_dps_variant)
        avg_dp = sum_dp/(size(dp_num_array,2))
        push!(avg_dps_all,avg_dp)
    end

    return avg_dps_all
end

"""
    list_sample_names_low_dp(sample_avg_list::Array{Float64,2},sample_names)
returns list of sample ids that have an average read depth less than 15 across all variant positions. Developers can implement this in the VIVA script - search script for function name and read notes.
"""
function list_sample_names_low_dp(sample_avg_list::Array{Float64,1},sample_names)

    low_dp_index_list = Array{Int64}(undef,0)

        for item = 1:length(sample_avg_list)
            if sample_avg_list[item] < 15 #developers can adjust this read depth threshold value to list positions with dp below threshold
                push!(low_dp_index_list,item)
            end
        end

    low_dp_sample_names = Array{String}(undef,0)

        for i in low_dp_index_list
            sample_name = sample_names[i]
            string_sample_name=string(sample_name)
            push!(low_dp_sample_names,string_sample_name)
        end

        return low_dp_sample_names
end

"""
    list_variant_positions_low_dp(variant_avg_list::Array{Float64,2},chrom_labels)
finds variant positions that have an average read depth less than 15 across all patients. Developers can implement this in the VIVA script - search script for function name and read notes.
"""
function list_variant_positions_low_dp(variant_avg_list::Array{Float64,1},chrom_labels)

    low_dp_index_list = Array{Int64}(undef,0)

        for item = 1:length(variant_avg_list)
            if variant_avg_list[item] < 15 #developers can adjust this read depth threshold value to list positions with dp below threshold
                push!(low_dp_index_list,item)
            end
        end

    low_dp_positions = Array{Tuple{Any,Int64}}(undef,0)

        for i in low_dp_index_list
            chrom_position = chrom_labels[i,1],chrom_labels[i,2]
            push!(low_dp_positions,chrom_position)
        end

        return low_dp_positions
end

#functions for producing objects for plot functions

"""
    read_depth_threshhold(dp_array::Array{Int64,2})
Caps read depth outlier values at user defined threshhold.
threshhold defaults to dp = 100. All dp over 100 are set to 100 to visualize
read depth values between 0 < dp > 100 in better definition.
"""
function read_depth_threshhold(dp_array::Array{Int64,2})

    dp_array[dp_array[:,:].>100].=100

    return dp_array

end

"""
    save_numerical_array(num_array::Matrix{Any},sample_names,chr_labels,title,output_directory)
save numerical array with chr labels and sample ids to working directory
"""
function save_numerical_array(num_array,sample_names,chr_labels,title,output_directory)

      headings = hcat("chr,position")
      sample_names = hcat(headings,sample_names)

      chr_labeled_array_for_plotly=hcat(chr_labels, num_array)

      labeled_value_matrix_withsamplenames= vcat(sample_names,
                                            chr_labeled_array_for_plotly)

      writedlm(joinpath(output_directory ,"$(title)_matrix.csv"), labeled_value_matrix_withsamplenames, ",")
end

"""
    chromosome_label_generator(chromosome_labels::Array{Any,1})
Returns vector of chr labels and indices to mark chromosomes in plotly heatmap
Specifically, saves indexes and chrom labels in vectors to pass into heatmap
function to ticvals and tictext respectively.
Input is either gt_chromosome_labels or dp_chromosome_labels from
translate_gt/dp_to_num_array()
"""
function chromosome_label_generator(chromosome_labels::Array{Any,1})
    chrom_label_indices = findfirst.(map(a -> (y -> isequal(a, y)),
    unique(chromosome_labels)), [chromosome_labels])
    chrom_labels = unique(chromosome_labels)
    chrom_labels = [string(i) for i in chrom_labels]

    if length(chrom_labels) > 1
        for item=2:(length(chrom_labels))

            ratio=((chrom_label_indices[item])-(chrom_label_indices[item-1]))/(length(chromosome_labels))

            if ratio < 0.2
                font_size = "8"
                return chrom_labels,chrom_label_indices,font_size

            else
                font_size = "10"
                return chrom_labels,chrom_label_indices,font_size

            end
        end

    else

        font_size = "10"

        return chrom_labels,chrom_label_indices,font_size
    end
end

"""
    checkfor_outputdirectory(path::String)
Checks to see if output directory exists already. If it doesn't, it creates the
new directory to write output files to.
"""
function checkfor_outputdirectory(path::String)
    if isdir(path) == true
    else
            mkdir(path)
    end
end

"""
    generate_chromosome_positions_for_hover_labels(chr_labels::Array{Any,2})
creates tuple of genomic locations to set as tick labels. This is automatically
store chromosome positions in hover labels. However tick labels are set to hidden with showticklabels=false so they will not crowd the y axis.
"""
function generate_chromosome_positions_for_hover_labels(chr_labels::Array{Any,2})

returnXY_column1!(chr_labels)

chr_pos_tuple_list=Array{String}(undef,0)

    for row = 1:size(chr_labels,1)

        chr=chr_labels[row,1]
        pos=chr_labels[row,2]
        chr_pos_tuple="chr$chr,$pos"
        push!(chr_pos_tuple_list,chr_pos_tuple)

    end

    return chr_pos_tuple_list
end

"""
    index_vcf(vcf_filename)
Creates and saves index file with three column array of vcf chrom, position, and row number to be used by significant list filter functions.
"""
function index_vcf(vcf_filename)

 reader = VCF.Reader(open(vcf_filename, "r"))

 chr_pos_col = Array{Any}(undef,0)
 row_col = Array{Any}(undef,0)

 row=1

  for record in reader

   chr = VCF.chrom(record)
   pos = VCF.pos(record)
   chr_pos_tuple = chr,pos
   row = row+1

   push!(chr_pos_col,chr_pos_tuple)
   push!(row_col,row)

  end

  index_matrix=hcat(chr_pos_col,row_col)
  return index_matrix

end

"""
    match_siglist_to_index(sig_list,vcf_index)
Returns vcf row indices of each variant position in sig_list for reader function to allow fast filtering in significant list filter funcitons.
"""
function match_siglist_to_index(sig_list,index_matrix)

 sig_list_tuples = Array{Any}(undef,0)

  for row = 1:size(sig_list,1)

   chrom = sig_list[row,1]
   pos = sig_list[row,2]
   sig_list_tuple = chrom,pos
   push!(sig_list_tuples,sig_list_tuple)

  end

  sig_list_index = findin(index_matrix,sig_list_tuples)

  return sig_list_index

end

"""
    make_chromosome_labels(chrom_label_info)
Returns vector of values to use as tick vals to show first chromosome label per chromosome with blank spaces between each first chromosome position for use with --y_axis_labels=chromosomes. duplicate_last_label tells if last chrom label is single or mutiple which affects number_to_fill value.
"""
function make_chromosome_labels(chrom_label_info)

       labels=chrom_label_info[1]
       index=chrom_label_info[2]
       duplicate_last_label=chrom_label_info[4]

       x=Array{Any}(undef,0)

       for n = 1:size(labels,1)

              push!(x, labels[n])
              counter=0
              number_to_fill=index[n+1]-index[n]

              while counter < number_to_fill - 1
                  counter = counter+1
                  push!(x, labels[n])
              end
        end

    if duplicate_last_label != "true"
      push!(x,labels[size(labels,1)])
    end

    return x

end

"""
    add_pheno_matrix_to_gt_data_for_plotting(pheno_matrix,gt_num_array,trait_labels,chrom_label_info,number_rows)
add the pheno matrix used to group samples to the data array for input into plotting functions. Resizes the pheno matrix to maintain correct dimensions for heatmap viz by finding value=0.05*number_rows_data to multiply each pheno row by before vcat.
"""
function add_pheno_matrix_to_gt_data_for_plotting(gt_num_array,pheno_matrix,trait_labels,chrom_label_info,number_rows)

    pheno_matrix=pheno_matrix[2:size(pheno_matrix,1),:]

    #find multiplyer value
    pheno_row_multiplyer=0.05*(size(gt_num_array,1))

    pheno_row_multiplyer=(pheno_row_multiplyer)

    pheno_row_multiplyer=round(pheno_row_multiplyer/(size(pheno_matrix,1)))

    if pheno_row_multiplyer < 1
        pheno_row_multiplyer = 1
    end

    for i = 1:(size(pheno_matrix, 1))
        for n=1:(size(pheno_matrix, 2))
            pheno_matrix[i, n] = pheno_matrix[i, n] == 1 ? -1 : pheno_matrix[i, n]
            pheno_matrix[i, n] = pheno_matrix[i, n] == 2 ? -2 : pheno_matrix[i, n]
        end
    end

    #resize pheno matrix and create trait_label_array for labeling

    resized_pheno_matrix=pheno_matrix[1:1, :]

    trait_label_array=trait_labels[1]

    for row = 1:size(pheno_matrix,1)

        trait_row=pheno_matrix[row:row, :]
        trait = trait_labels[row]
        i=0

        while i < (pheno_row_multiplyer)

            i=i+1

            resized_pheno_matrix=vcat(resized_pheno_matrix,trait_row)
            trait_label_array=vcat(trait_label_array,trait)
        end

    end

    resized_pheno_matrix = resized_pheno_matrix[2:end,:]

    #if too few variants are selected for visualization, there is only 1 row of trait labels, so trait label must be stored in array for unique() below

    if typeof(trait_label_array) == SubString{String}

    trait_label_array=[trait_label_array]

    end

    trait_label_array = trait_label_array[2:end]

    trait_label_indices = findfirst.(map(a -> (y -> isequal(a, y)),
    unique(trait_label_array)), [trait_label_array])

    pheno_trait_labels = unique(trait_label_array)

    pheno_trait_labels = [String(i) for i in pheno_trait_labels]
    trait_label_indices = [i+number_rows for i in trait_label_indices]

    chrom_labels=chrom_label_info[1]
    chrom_label_indices=chrom_label_info[2]

    chrom_labels = vcat(chrom_labels,pheno_trait_labels)
    chrom_label_indices = vcat(chrom_label_indices,trait_label_indices)

    chrom_label_info = chrom_labels,chrom_label_indices,chrom_label_info[3]

    data_and_pheno_matrix=vcat(gt_num_array,resized_pheno_matrix)

    data_and_pheno_matrix=convert(Array{Int64,2},data_and_pheno_matrix)

    return data_and_pheno_matrix,trait_label_array,chrom_label_info

end

"""
    add_pheno_matrix_to_dp_data_for_plotting(pheno_matrix,dp_num_array,trait_labels,chrom_label_info,number_rows)
add the pheno matrix used to group samples to the data array for input into plotting functions. Resizes the pheno matrix to maintain correct dimensions for heatmap viz by finding value=0.05*number_rows_data to multiply each pheno row by before vcat.
"""
function add_pheno_matrix_to_dp_data_for_plotting(dp_num_array,pheno_matrix,trait_labels,chrom_label_info,number_rows)

    pheno_matrix=pheno_matrix[2:size(pheno_matrix,1),:]

    #find multiplyer value
    pheno_row_multiplyer=0.05*(size(dp_num_array,1))
    pheno_row_multiplyer=round(pheno_row_multiplyer)
    pheno_row_multiplyer=pheno_row_multiplyer/(size(pheno_matrix,1))

    if pheno_row_multiplyer < 1
        pheno_row_multiplyer = 1
    end

    for i = 1:(size(pheno_matrix, 1))
        for n=1:(size(pheno_matrix, 2))
            pheno_matrix[i, n] = pheno_matrix[i, n] == 1 ? -40 : pheno_matrix[i, n]
            pheno_matrix[i, n] = pheno_matrix[i, n] == 2 ? -60 : pheno_matrix[i, n]
        end
    end

    #resize pheno matrix and create trait_label_array for labeling
    resized_pheno_matrix=pheno_matrix[1:1, :]

    trait_label_array=trait_labels[1]

    for row = 1:size(pheno_matrix,1)

        trait_row=pheno_matrix[row:row, :]
        trait = trait_labels[row]

        i=0

        while i < (pheno_row_multiplyer)

            i=i+1

            resized_pheno_matrix=vcat(resized_pheno_matrix,trait_row)
            trait_label_array=vcat(trait_label_array,trait)
        end

    end

        #if too few variants are selected for visualization, there is only 1 row of trait labels, so trait label must be stored in array for unique() below
        resized_pheno_matrix = resized_pheno_matrix[2:end,:]

        if typeof(trait_label_array) == SubString{String}

        trait_label_array=[trait_label_array]

        end

        trait_label_array = trait_label_array[2:end]

    trait_label_indices = findfirst.(map(a -> (y -> isequal(a, y)),
    unique(trait_label_array)), [trait_label_array])

    pheno_trait_labels = unique(trait_label_array)
    pheno_trait_labels = [String(i) for i in pheno_trait_labels]
    trait_label_indices = [i+number_rows for i in trait_label_indices]

    chrom_labels=chrom_label_info[1]
    chrom_label_indices=chrom_label_info[2]

    chrom_labels = vcat(chrom_labels,pheno_trait_labels)
    chrom_label_indices = vcat(chrom_label_indices,trait_label_indices)

    chrom_label_info = chrom_labels,chrom_label_indices,chrom_label_info[3]

    #vcat pheno matrix to gt_num_array
    data_and_pheno_matrix=vcat(dp_num_array,resized_pheno_matrix)

    data_and_pheno_matrix=convert(Array{Int64,2},data_and_pheno_matrix)

    return data_and_pheno_matrix,trait_label_array,chrom_label_info
end

"""
    generate_hover_text_array(chr_pos_tuple_list,sample_names,input,mode)
Generate array of data for hovertext to use as custom hover text for ungrouped heatmaps. Where mode is GT or DP.
"""
function generate_hover_text_array(chr_pos_tuple_list,sample_names,input,mode)

hover_text_array=Array{Any,1}(undef,0)

    for n=1:length(chr_pos_tuple_list)
        pos=chr_pos_tuple_list[n]

        row_data_vector=Array{Any,1}(undef,0)

        for m=1:length(sample_names)
            id=sample_names[m]

            z=input[n,m]

            if mode == "GT"
            z_dict = Dict(0 => "No Call", 1 => "Homozygous Reference", 2 => "Heterozygous Variant", 3 => "Homozygous Variant")
            z=z_dict[z]

            elseif mode == "DP"
                if z==-20
                    z="No Call"
                end
            end

            text_cell="ID: $id \n $pos \n $mode: $z"

            push!(row_data_vector,text_cell)
        end

        push!(hover_text_array,row_data_vector)
    end

return(hover_text_array)

end

"""
    generate_hover_text_array_grouped(chr_pos_tuple_list,sample_names,input,mode)
Generate array of data for hovertext to use as custom hover text for grouped heatmaps. Where mode is GT or DP.
"""
function generate_hover_text_array_grouped(chr_pos_tuple_list,sample_names,input,mode,number_rows)

    hover_text_array=Array{Any,1}(undef,0)

    for n=1:number_rows
        pos=chr_pos_tuple_list[n]

        row_data_vector=Array{Any,1}(undef,0)

        for m=1:length(sample_names)
            id=sample_names[m]

            z=input[n,m]

            if mode == "GT"
            z_dict = Dict(0 => "No Call", 1 => "Homozygous Reference", 2 => "Heterozygous Variant", 3 => "Homozygous Variant")
            z=z_dict[z]

            elseif mode == "DP"
                if z==-1
                    z="No Call"
                elseif  z==100
                    z="100+"
                end

            end

            text_cell="ID: $id \n $pos \n $mode: $z"

            push!(row_data_vector,text_cell)
        end

        push!(hover_text_array,row_data_vector)
    end

    for n=number_rows+1:length(chr_pos_tuple_list)

        trait=chr_pos_tuple_list[n]

        split_trait=split(trait,',')

        row_data_vector=Array{Any,1}(undef,0)

        for m=1:length(sample_names)
            id=sample_names[m]

            z=input[n,m]

            if mode == "GT"
                if z==-1
                    z="Trait 1: $(split_trait[1])"
                elseif z==-2
                    z="Trait 2: $(split_trait[2])"
                end

            elseif mode == "DP"
                if z==-40
                    z="Trait 1: $(split_trait[1])"
                elseif z==-60
                    z="Trait 2: $(split_trait[2])"
                end
            end

            text_cell="ID: $id \n $z"

            push!(row_data_vector,text_cell)
        end

        push!(hover_text_array,row_data_vector)
    end

return(hover_text_array)

end

"""
    save_graphic(graphic,output_directory,save_ext,title,remote_option)
Save plot in either html or static image formats incuding eps, png, svg, and pdf
"""
function save_graphic(graphic,output_directory,save_ext,title,remote_option)

    if save_ext=="html" && remote_option == true
        PlotlyJS.savehtml(graphic, joinpath("$(output_directory)" ,"$title.$(save_ext)"), :remote)
    elseif save_ext=="html" && remote_option != true
        PlotlyJS.savehtml(graphic, joinpath("$(output_directory)" ,"$title.$(save_ext)"))
    elseif save_ext != "html"
        PlotlyJS.savefig(graphic, joinpath("$(output_directory)" ,"$title.$(save_ext)"))
    end

    display(graphic)
end


text:
module Types

function supertypes(type::DataType)
    types = Vector{DataType}()
    a_type = type
    while a_type.super ≠ Any
        push!(types, a_type.super)
        a_type = a_type.super
    end
    return types
end

function isasupertype(type::DataType, stype::DataType)
    return any(type -> type == stype, supertypes(type))
end

export supertypes
export isasupertype

end  # module Types

ext: jl


text:
struct BitFields{N,U} <: AbstractBitFields  # N BitField{U} fields
   bitfields::NTuple{N,BitField{U}}

   function BitFields(bitfields::Vararg{BitField{U}, N}) where {N,U}
       return new{N,U}(bitfields)
   end
end

Base.eltype(::Type{BitFields{N,U}}) where {N,U} = U
Base.eltype(x::BitFields{N,U}) where {N,U} = U

Base.length(x::BitFields{N,U}) where {N,U} = N

Base.lastindex(x::BitFields{N,U}) where {N,U} =
    lastindex(x.bitfields)

Base.getindex(x::BitFields{N,U}, idx::Int) where {N,U} =
    getindex(x.bitfields, idx)

Base.getindex(x::BitFields{N,U}, idxs::AbstractRange{Int}) where {N,U} =
    getindex(x.bitfields, idxs)

function Base.iterate(x::BitFields{N,U}) where {N,U<:UBits}
    if length(x) > 0
       x[1], 2
    else
       nothing
    end
end

function Base.iterate(x::BitFields{N,U}, state::Int) where {N,U<:UBits}
    if state > length(x)
       nothing
    else
       x[state], state+1
    end
end

Base.get(bitfields::BitFields{N,U}, bits::Ref{U}) where {N,U<:UBits} = [get(bitfield, bits) for bitfield in bitfields]
Base.get(bitfields::BitFields{N,U}, bits::U) where {N,U<:UBits} = [get(bitfield, bits) for bitfield in bitfields]

names(bitfields::BitFields{N,U}) where {N,U<:UBits} = ((name(bitfield) for bitfield in bitfields)...,)

function Base.NamedTuple(bitfields::BitFields{N,U}) where {N,U<:UBits}
   bitfieldnames  = names(bitfields)
   values = ((bitfields)...,)
   nt = NamedTuple{bitfieldnames,NTuple{N,BitField}}(values)
   return nt
end

"""
    valuesoffields

obtain values assigned to fields of a struct type (in field order)
<from NamedTupleTools>
"""
valuesoffields(x::T) where {T<:NamedTuple} =
    ((getfield(x, name) for name in fieldnames(T))...,)

Base.NamedTuple(x::NamedTuple{S,NTuple{N,BitField}}, z::Ref{T}) where {S,N,T} =
    NamedTuple{fieldnames(typeof(x))}(get.(valuesoffields(x), z[]))

Base.NamedTuple(x::NamedTuple{S,NTuple{N,BitField}}, z::ByRef{T}) where {S,N,T} =
    NamedTuple{fieldnames(typeof(x))}(get.(valuesoffields(x), z.ref[]))

function BitFields(bitfields::NamedTuple, nt::NamedTuple)
    z = ByRef(eltype(bitfields[1]))
    for (fld,val) in zip(bitfields, nt)
        set!(fld, val, z)
    end
    return z
end


text:
function test_point_sampling()
    radius = 1.0f0
    origin = Point3f0(0, 0, 0)
    mesh = icosphere(4, radius)
    verts = coordinates(mesh)
    triangles = faces(mesh)
    areas_array = area.(Ref(verts), triangles)
    cumsum_area = similar(areas_array)
    cumsum_area = cumsum!(cumsum_area, areas_array)
    total_area = cumsum_area[end]
    bins = cumsum_area / total_area
    point1 = point_sampling(verts, triangles, bins)
    d_point1 = norm(point1 - origin)
    @test isapprox(d_point1, radius, atol = 0.001)
    points1 = points_sampling(mesh, 10000)
    points_dist = norm.(points1)
    mean_dist = mean(points_dist)
    sd_dist = std(points_dist)
    @test isapprox(mean_dist, 1.0f0, atol = 0.001)
    @test isapprox(sd_dist, 0.0f0, atol = 0.001)
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _calz_unred{T<:AbstractFloat}(wave::T, flux::T, ebv::T, r_v::T)
    x  = 10000.0/wave # Wavelength in inverse microns
    if 6300.0 <= wave <= 22000.0
        klam = 2.659*(-1.857 + 1.040*x) + r_v
    elseif 912.0 <= wave < 6300.0
        klam = 2.659*(@evalpoly(x, -2.156, 1.509, -0.198, 0.011)) + r_v
    else
        return flux
    end
    return flux*10.0^(0.4*klam*ebv)
end

"""
    calz_unred(wave, flux, ebv[, r_v]) -> deredden_wave

### Purpose ###

Deredden a galaxy spectrum using the Calzetti et al. (2000) recipe.

### Explanation ###

Calzetti et al.  (2000, ApJ 533, 682;
http://adsabs.harvard.edu/abs/2000ApJ...533..682C) developed a recipe for
dereddening the spectra of galaxies where massive stars dominate the radiation
output, valid between \$0.12\$ to \$2.2\$ microns.  (`calz_unred` extrapolates
between \$0.12\$ and \$0.0912\$ microns.)

### Arguments ###

* `wave`: wavelength vector (Angstroms)
* `flux`: calibrated flux vector, same number of elements as `wave`.
* `ebv`: color excess E(B-V), scalar.  If a negative `ebv` is supplied, then
  fluxes will be reddened rather than deredenned.  Note that the supplied color
  excess should be that derived for the stellar continuum, EBV(stars), which is
  related to the reddening derived from the gas, EBV(gas), via the Balmer
  decrement by EBV(stars) = 0.44*EBV(gas).
* `r_v` (optional): scalar ratio of total to selective extinction, default is
 4.05.  Calzetti et al. (2000) estimate \$r_v = 4.05 \pm 0.80\$ from optical-IR
 observations of 4 starbursts.

### Output ###

Unreddened flux vector, same units and number of elements as `flux`.  Flux
values will be left unchanged outside valid domain (\$0.0912\$ - \$2.2\$
microns).

### Example ###

Estimate how a flat galaxy spectrum (in wavelength) between \$1200 Å\$ and
\$3200 Å\$ is altered by a reddening of E(B-V) = 0.1.

``` julia
wave = reshape(1200:50:3150,40);
flux = ones(wave);
flux_new = calz_unred(wave, flux, -0.1);
```

Using a plotting tool you can visualize the unreddend flux.  For example, with
[PyPlot.jl](https://github.com/stevengj/PyPlot.jl)

``` julia
using PyPlot
plot(wave, flux_new)
```

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
calz_unred(wave::Real, flux::Real, ebv::Real, r_v::Real=4.05) =
    _calz_unred(promote(float(wave), float(flux), float(ebv), float(r_v))...)

function calz_unred{W<:Real, F<:Real}(wave::AbstractArray{W},
                                      flux::AbstractArray{F},
                                      ebv::Real, r_v::Real=4.05)
    @assert length(wave) == length(flux)
    flux_unred  = similar(flux, typeof(float(one(W))))
    for i in eachindex(flux)
        flux_unred[i] = calz_unred(wave[i], flux[i], ebv, r_v)
    end
    return flux_unred
end


text:
module FailBuild

end # module


text:
# Example of neural network
# using Revise
using PiecewiseDeterministicMarkovProcesses, LinearAlgebra, Random, SparseArrays
const PDMP = PiecewiseDeterministicMarkovProcesses

const N = 100

function f(x)
	return	x^8
end

function Phi(out::Array{Float64,2}, xc, xd, parms, t::Array{Float64})
	# vector field used for the continuous variable
	# for this particular model, the empirical mean is constant between jumps
	λ = 0.24
	xbar::Float64 = sum(xc) / N
	out[1,:] .= xc
	out[2,:] .= xbar .+ exp(-λ*(t[2]-t[1])) .* (xc .- xbar)
	nothing
end

function R_mf_rejet(rate, xc, xd, parms, t::Float64, issum::Bool)
	bound = N * f(1.201)#1.5 works well
	# rate function
	if issum == false
		for i=1:N
			rate[i] = f(xc[i])
		end
		return -1., bound
	else
		res = 0.
		for i=1:N
			res += f(xc[i])
		end
		return res, bound
	end
end

function Delta_xc_mf(xc, xd, parms, t::Float64, ind_reaction::Int64)
	# this function return the jump in the continuous component
	J = 0.98
	for i=1:N
		xc[i] += J/N
	end
	xc[ind_reaction] = 0.0
	xd[ind_reaction] += 1
	return true
end

Random.seed!(1234)
xc0 = rand(N)*0.2 .+ 0.5
xd0 = zeros(Int64, N)

nu_neur = spzeros(Int64,N,N)
parms = [0.1]
tf = 10_050.

problem = PDMP.PDMPProblem(Phi,R_mf_rejet,Delta_xc_mf,nu_neur, xc0, xd0, parms, (0.0, tf))
	Random.seed!(8)
	result = PDMP.solve(problem, PDMP.RejectionExact(); n_jumps = 10_000, ind_save_d = 1:2, ind_save_c = 1:2)

result = PDMP.solve(problem, PDMP.RejectionExact(); n_jumps = 10_000, ind_save_d = 1:2, ind_save_c = 1:2)


text:
using Images, ImageShow
using Test

@testset "show (MIME)" begin
    # Test that we remembered to turn off Colors.jl's colorswatch display
    @test !showable(MIME("image/svg+xml"), rand(Gray{N0f8}, 5, 5))
    @test !showable(MIME("image/svg+xml"), rand(RGB{N0f8},  5, 5))
    @test showable(MIME("image/png"), rand(Gray{N0f8}, 5, 5))
    @test showable(MIME("image/png"), rand(RGB{N0f8},  5, 5))
    workdir = joinpath(tempdir(), "Images")
    if !isdir(workdir)
        mkdir(workdir)
    end
    @testset "no compression or expansion" begin
        A = N0f8[0.01 0.99; 0.25 0.75]
        fn = joinpath(workdir, "writemime.png")
        open(fn, "w") do file
            show(file, MIME("image/png"), Gray.(A), minpixels=0, maxpixels=typemax(Int))
        end
        b = load(fn)
        @test b == A
        @test Images.metadata(fn) == (reverse(size(A)), Gray{N0f8})

        img = fill(RGB{N0f16}(1,0,0), 1, 1)
        open(fn, "w") do file
            show(IOContext(file, :full_fidelity=>true), MIME("image/png"), img, minpixels=0, maxpixels=typemax(Int))
        end
        b = load(fn)
        @test b == img
        if VERSION >= v"1.3"
            # ImageMagick backend is broken here.
            # For Julia >=v1.3 FileIO will use ImageIO backend
            @test Images.metadata(fn) == (reverse(size(img)), RGB{N0f16})
        end

        A = N0f8[0.01 0.99; 0.25 0.75]
        fn = joinpath(workdir, "writemime.png")
        open(fn, "w") do file
            show(IOContext(file, :full_fidelity=>true), MIME("image/png"), Gray.(A), minpixels=5, maxpixels=typemax(Int))
        end
        @test load(fn) == A
        @test Images.metadata(fn) == (reverse(size(A)), Gray{N0f8})

        A = N0f8[0.01 0.4 0.99; 0.25 0.8 0.75; 0.6 0.2 0.0]
        fn = joinpath(workdir, "writemime.png")
        open(fn, "w") do file
            show(IOContext(file, :full_fidelity=>true), MIME("image/png"), Gray.(A), minpixels=0, maxpixels=5)
        end
        @test load(fn) == A
        @test Images.metadata(fn) == (reverse(size(A)), Gray{N0f8})

        # a genuinely big image (tests the defaults)
        abig = colorview(Gray, normedview(rand(UInt8, 1024, 1023)))
        fn = joinpath(workdir, "big.png")
        open(fn, "w") do file
            show(IOContext(file, :full_fidelity=>true), MIME("image/png"), abig, maxpixels=10^6)
        end
        b = load(fn)
        @test b == abig
        @test Images.metadata(fn) == (reverse(size(abig)), Gray{N0f8})
    end
    @testset "colorspace normalization" begin
        img = fill(HSV{Float64}(0.5, 0.5, 0.5), 1, 1)
        fn = joinpath(workdir, "writemime.png")
        open(fn, "w") do file
            show(file, MIME("image/png"), img, minpixels=0, maxpixels=typemax(Int))
        end
        b = load(fn)
        @test b == convert(Array{RGB{N0f8}}, img)
        @test Images.metadata(fn) == (reverse(size(img)), RGB{N0f8})
        img = fill(RGB{N0f16}(1,0,0), 1, 1)
        open(fn, "w") do file
            show(file, MIME("image/png"), img, minpixels=0, maxpixels=typemax(Int))
        end
        b = load(fn)
        @test eltype(b) <: AbstractRGB && eltype(eltype(b)) == N0f8 && b[1] == RGB(1,0,0)
        @test Images.metadata(fn) == (reverse(size(img)), RGB{N0f8})
        img = fill(RGBA{Float32}(1,0,0,0.5), 1, 1)
        open(fn, "w") do file
            show(file, MIME("image/png"), img, minpixels=0, maxpixels=typemax(Int))
        end
        b = load(fn)
        @test isa(b, Matrix{RGBA{N0f8}}) && b[1] == RGBA{N0f8}(1,0,0,0.5)
        @test Images.metadata(fn) == (reverse(size(img)), RGBA{N0f8})
        img = Gray.([0.1 0.2; -0.5 0.8])
        open(fn, "w") do file
            show(file, MIME("image/png"), img, minpixels=0, maxpixels=typemax(Int))
        end
        b = load(fn)
        @test isa(b, Matrix{Gray{N0f8}}) && b == Gray{N0f8}[0.1 0.2; 0 0.8]
        @test Images.metadata(fn) == (reverse(size(img)), Gray{N0f8})
    end
    @testset "small images (expansion)" begin
        A = N0f8[0.01 0.99; 0.25 0.75]
        fn = joinpath(workdir, "writemime.png")
        open(fn, "w") do file
            show(file, MIME("image/png"), Gray.(A), minpixels=5, maxpixels=typemax(Int))
        end
        @test load(fn) == A[[1,1,2,2],[1,1,2,2]]
        @test Images.metadata(fn) == ((4,4), Gray{N0f8})
    end
    @testset "big images (use of restrict)" begin
        A = N0f8[0.01 0.4 0.99; 0.25 0.8 0.75; 0.6 0.2 0.0]
        Ar = restrict(A)
        fn = joinpath(workdir, "writemime.png")
        open(fn, "w") do file
            show(file, MIME("image/png"), Gray.(A), minpixels=0, maxpixels=5)
        end
        @test load(fn) == N0f8.(Ar)
        @test Images.metadata(fn) == (reverse(size(Ar)), Gray{N0f8})
        # a genuinely big image (tests the defaults)
        abig = colorview(Gray, normedview(rand(UInt8, 1024, 1023)))
        fn = joinpath(workdir, "big.png")
        open(fn, "w") do file
            show(file, MIME("image/png"), abig, maxpixels=10^6)
        end
        b = load(fn)
        @test b == N0f8.(restrict(abig, (1,2)))
        @test Images.metadata(fn) == (reverse(size(restrict(abig,(1,2)))), Gray{N0f8})
    end
    @testset "display matrix of images" begin
        img() = colorview(Gray, rand([0.250 0.5; 0.75 1.0], rand(2:10), rand(2:10)))
        io = IOBuffer()
        # test that these methods don't fail
        show(io, MIME"text/html"(), [img() for i=1:2])
        show(io, MIME"text/html"(), [img() for i=1:2, j=1:2])
        show(io, MIME"text/html"(), [img() for i=1:2, j=1:2, k=1:2])
    end
    @testset "display matrix of 1-D images" begin
        flat_imgs = [zeros(Gray{Float32}, 1)]
        io = IOBuffer()
        # These methods should not invoke the Images.jl display code, but they
        # used to throw errors: https://github.com/JuliaImages/Images.jl/issues/623
        @test !applicable(ImageShow._show_odd, io, MIME"text/html"(), flat_imgs)
        @test !applicable(ImageShow._show_even, io, MIME"text/html"(), flat_imgs)
    end
    rm(workdir, recursive=true)
end

nothing


text:
using OptimalEstimationProject
using AstroTime

# GPS Simulation Span
startWeek   = 2170
startDay    = 0
endWeek     = 2174
endDay      = 6

# Create GPS simulation object
gpssim = GPSSim(startWeek, startDay, endWeek, endDay)

# Generate measurement
measEpoch = TAIEpoch("2021-08-25T12:00:00")
σρ        = 0.1/1000;
r         = [6700.0, 0.0, 0.0]
mSecTIA   = Float64(measEpoch.second)
meas = OptimalEstimationProject.gpsMeasurement(gpssim, r, mSecTIA, σρ)

# Spacecraft simulator
t0        = TAIEpoch("2021-08-08T12:00:00")
scsim     = OptimalEstimationProject.SpacecraftSim(t0.second)

# Get state and control
(xs,us)   = OptimalEstimationProject.GetStateAndControl(scsim, t0.second .+ [t for t in 0:10:1000])

# Plotting
OptimalEstimationProject.PlotTrajectory(scsim; frame = :inertial)
#OptimalEstimationProject.PlotTrajectory(scsim, range(t0.second + 400; stop = t0.second + 800, length = 100))

max_stars_repo_name: GrantHecht/OptimalEstimationProject


text:
using Documenter, FrankWolfe
using SparseArrays
using LinearAlgebra

makedocs(
    modules=[FrankWolfe],
    sitename="FrankWolfe.jl",
    format=Documenter.HTML(prettyurls=get(ENV, "CI", nothing) == "true"),
    pages=[
        "Home" => "index.md",
        "Examples" => "examples.md",
        "References" => "reference.md",
        "Index" => "indexlist.md",
    ],
)

deploydocs(repo="github.com/ZIB-IOL/FrankWolfe.jl.git", push_preview=true)


text:
module QuantumESPRESSOCommands

using AbInitioSoftwareBase: parentdir
using AbInitioSoftwareBase.Commands: mpiexec
@static if VERSION >= v"1.6"
    using Preferences: @load_preference, @set_preferences!
end

@static if VERSION >= v"1.6"
    function get_path(exe)
        if exe == "pw"
            return @load_preference("pw.x path", "pw.x")
        elseif exe == "ph"
            return @load_preference("ph.x path", "ph.x")
        elseif exe == "q2r"
            return @load_preference("q2r.x path", "q2r.x")
        elseif exe == "matdyn"
            return @load_preference("matdyn.x path", "matdyn.x")
        elseif exe == "dynmat"
            return @load_preference("dynmat.x path", "dynmat.x")
        else
            throw(ArgumentError("invalid option $exe."))
        end
    end
    function set_path(exe, path::String)
        @assert ispath(path)
        if exe == "pw"
            @set_preferences!("pw.x path" => path)
        elseif exe == "ph"
            @set_preferences!("ph.x path" => path)
        elseif exe == "q2r"
            @set_preferences!("q2r.x path" => path)
        elseif exe == "matdyn"
            @set_preferences!("matdyn.x path" => path)
        elseif exe == "dynmat"
            @set_preferences!("dynmat.x path" => path)
        else
            throw(ArgumentError("invalid option $exe."))
        end
    end
else
    function get_path(exe)
        if exe == "pw"
            return "pw.x"
        elseif exe == "ph"
            return "ph.x"
        elseif exe == "q2r"
            return "q2r.x"
        elseif exe == "matdyn"
            return "matdyn.x"
        elseif exe == "dynmat"
            return "dynmat.x"
        else
            throw(ArgumentError("invalid option $exe."))
        end
    end
end

export pw, ph, q2r, matdyn, dynmat

"""
    pw(input, output; chdir, nimage, npool, ntg, nyfft, nband, ndiag, np, env, kwargs...)

Make commands for `pw.x`.
"""
pw(input, output = mktemp(parentdir(input))[1]; kwargs...) =
    cmdtemplate(get_path("pw"), input, output; kwargs...)
"""
    ph(input, output; chdir, nimage, npool, ntg, nyfft, nband, ndiag, np, env, kwargs...)

Make commands for `ph.x`.
"""
ph(input, output = mktemp(parentdir(input))[1]; kwargs...) =
    cmdtemplate(get_path("ph"), input, output; kwargs...)
"""
    q2r(input, output; chdir, nimage, npool, ntg, nyfft, nband, ndiag, np, env, kwargs...)

Make commands for `q2r.x`.
"""
q2r(input, output = mktemp(parentdir(input))[1]; kwargs...) =
    cmdtemplate(get_path("q2r"), input, output; kwargs...)
"""
    matdyn(input, output; chdir, nimage, npool, ntg, nyfft, nband, ndiag, np, env, kwargs...)

Make commands for `matdyn.x`.
"""
matdyn(input, output = mktemp(parentdir(input))[1]; kwargs...) =
    cmdtemplate(get_path("matdyn"), input, output; kwargs...)
"""
    dynmat(input, output; chdir, nimage, npool, ntg, nyfft, nband, ndiag, np, env, kwargs...)

Make commands for `dynmat.x`.
"""
dynmat(input, output = mktemp(parentdir(input))[1]; kwargs...) =
    cmdtemplate(get_path("dynmat"), input, output; kwargs...)

"""
    cmdtemplate(path, input, output; chdir, nimage, npool, ntg, nyfft, nband, ndiag, np, env, kwargs...)

Make commands for QuantumESPRESSO executables.
"""
function cmdtemplate(
    path,
    input,
    output = mktemp(parentdir(input))[1];
    chdir = true,
    nimage = 0,
    npool = 0,
    ntg = 0,
    nyfft = 0,
    nband = 0,
    ndiag = 0,
    np = 1,
    env = [],
    kwargs...,
)
    if ndiag^2 > np
        @error "`ndiag` square should be less than `np`!"
    end
    f = mpiexec(env; np = np, kwargs...)
    args = [path]
    for (key, value) in zip(
        (:nimage, :npool, :ntg, :nyfft, :nband, :ndiag),
        (nimage, npool, ntg, nyfft, nband, ndiag),
    )
        if !iszero(value)
            push!(args, "-$key", string(value))
        end
    end
    dir = abspath(chdir ? parentdir(input) : pwd())
    return pipeline(Cmd(f(args); dir = dir); stdin = input, stdout = output)
end

end


text:
#returns a mol volume from P and T
function v_zero(mt::SinglePT,model,args...;kwargs...)
    return v_zero(volume_solver_type(model),mt,model,args...;kwargs...)
end

function v_zero(mt::MultiPT,model,args...;kwargs...)
    return v_zero(volume_solver_type(model),mt,model,args...;kwargs...)
end

function v_zero(::SUVA,mt::SinglePT,model,p,t,v0=nothing;phase=:unspecified)
    _vt = QuickStates.vt()
    _p(_v) = pressure_impl(_vt,model,_v,t)
    if is_gas(phase)
        if isnothing(v0)
            v = gas_fzero(_p,p,t)
        else
            v = gas_fzero(_p,p,t,v0)
        end
    elseif is_liquid(phase)
        if isnothing(v0)
            v = v_zero(VolumeBisection(),mt,model,p,t;phase=phase)
            #v = liquid_fzero(_p,p,t)
            return v
        else
            try
                v= roots_fzero(_p,p,t,v0)
                return v
            catch
                v = v_zero(VolumeBisection(),mt,model,p,t;phase=phase)
                #v = liquid_fzero(_p,p,t)
                return v
            end
        end

    else #return the phase with least gibbs energy
        if isnothing(v0)
            return v_zero_general(mt,model,p,t)
        else
            try
                return roots_fzero(_p,p,t,v0)
            catch
                return v_zero_general(mt,model,p,t)
            end
        end
    end
end
#returns all mol volumes from P and T

function v_zero(::SUVA,mt::MultiPT,model,p,t,x,v0=nothing;phase=:unspecified)
    _vt = QuickStates.vtx()
    _p(_v) = pressure_impl(_vt,model,_v,t,x)
    if is_gas(phase)
        if isnothing(v0)
            v = gas_fzero(_p,p,t)
        else
            v = gas_fzero(_p,p,t,v0)
        end

    elseif is_liquid(phase)
        if isnothing(v0)
            v = v_zero(VolumeBisection(),mt,model,p,t,x;phase=phase)
            #v = liquid_fzero(_p,p,t)
            return v
        else
            try
                v= roots_fzero(_p,p,t,v0)
                return v
            catch
                v = v_zero(VolumeBisection(),mt,model,p,t,x;phase=phase)
                #v = liquid_fzero(_p,p,t)
                return v
            end
        end
    else #return the phase with least gibbs energy
        if isnothing(v0)
            return v_zero_general(mt,model,p,t,x)

        else
            try
                return roots_fzero(_p,p,t,v0)
            catch

                return v_zero_general(mt,model,p,t,x)
            end
        end
    end
end
#returns all mol volumes from P and T
function v_zeros(mt,model,args...;kwargs...)
    return v_zeros(volume_solver_type(model),mt,model,args...;kwargs...)
end
function v_zeros(::SUVA
    ,mt::SinglePT
    ,model::HelmholtzModel
    ,p
    ,t)

    _vt = QuickStates.vt()
    _p(_v) = pressure_impl(_vt,model,_v,t)
    v1 = liquid_fzero(_p,p,t)
    v2 = gas_fzero(_p,p,t)
    return (v1,v2)
end

function v_zeros(::SUVA
    ,mt::MultiPT
    ,model
    ,p
    ,t
    ,x)

    _vt = QuickStates.vtx()
    _p(_v) = pressure_impl(_vt,model,_v,t,x)
    v1 = liquid_fzero(_p,p,t)
    v2 = gas_fzero(_p,p,t)
    return (v1,v2)
end

function v_zeros(vmethod::VolumeBisection
    ,mt::SinglePT
    ,model
    ,p
    ,t
    ;phase=:unspecified)
    no_pts = vmethod.pts
    _p(z) = pressure_impl(QuickStates.vt(),model, z, t)
    fp(z) = _p(z) - p
        max_v = 40 * t / p #approx 5 times ideal gas
    #this is to be sure that (min_v,max_v) is a bracketing interval
    while (fp(max_v)) > 0 |  (_p(max_v) < 0)
        max_v *= 2
        (max_v > typemax(max_v)*0.125) && break
    end

    vv = Roots.find_zeros(fp, eps(typeof(max_v)), max_v, no_pts = no_pts)
    return (first(vv),last(vv))
end

function v_zeros(vmethod::VolumeBisection,
    mt::MultiPT,
    model,
    p,
    t,
    x
    ;phase=:unspecified)
    no_pts = vmethod.pts
    _p(z) = pressure_impl(QuickStates.vtx(),model, z, t,x)
    fp(z) = _p(z) - p
        max_v = 40 * t / p #approx 5 times ideal gas
    #this is to be sure that (min_v,max_v) is a bracketing interval
    while (fp(max_v)) > 0 |  (_p(max_v) < 0)
        max_v *= 2
        (max_v > typemax(max_v)*0.125) && break
    end

    vv = Roots.find_zeros(fp, eps(typeof(max_v)), max_v, no_pts = no_pts)
    v1 = first(vv)
    v2 = last(vv)
    #v1 =  Roots.find_zero(fp,first(vv))
    #v2 = Roots.find_zero(fp,last(vv))
    return (v1,v2)
end

function v_zero(vmethod::VolumeBisection
    ,mt::SinglePT
    ,model
    ,p
    ,t
    ,v0=nothing
    ;phase=:unspecified)

    _p(z) = pressure_impl(QuickStates.vt(),model, z, t)

    if !(v0 === nothing)
        try
            return roots_fzero(_p,p,t,v0)
        catch
            return v_zero_general(mt,model,p,t)
        end
    end
    gas = is_gas(phase)
    liquid = is_liquid(phase)
    if !gas & !liquid
        return v_zero_general(mt,model,p,t)
    end

    (vl,vg) = v_zeros(vmethod,mt,model,p,t)
    if gas
        return vg
    else
        return vl
    end
end

function v_zero(vmethod::VolumeBisection
    ,mt::MultiPT
    ,model
    ,p
    ,t
    ,x
    ,v0=nothing
    ;phase=:unspecified)

    _p(z) = pressure_impl(QuickStates.vtx(),model, z, t,x)

    if !(v0 === nothing)
        try
            return roots_fzero(_p,p,t,v0)
        catch
            return v_zero_general(mt,model,p,t,x)
        end
    end
    gas = is_gas(phase)
    liquid = is_liquid(phase)
    if !gas & !liquid
        return v_zero_general(mt,model,p,t,x)
    end

    (vl,vg) = v_zeros(vmethod,mt,model,p,t,x)
    if gas
        return vg
    else
        return vl
    end
end

function v_zero_general(mt::SinglePT,model,p,t)
    _vt = QuickStates.vt()
    _g(_v) = mol_gibbs_impl(_vt,model,_v,t)
    _p(_v) = pressure_impl(_vt,model,_v,t)

    v1,v2 = v_zeros(mt,model,p,t)
    #(v1 == v2) && return v1
    g1 = _g(v1)
    g2 = _g(v2)
    if g1 < g2
        return v1
    elseif g2 < g1
        return v2
    else #different volumes, but equal energies, equilibria conditions
        return v2
    end
end

#finds the root with lowest gibbs energy, corrected for gas
function v_zero_general(mt::MultiPT,model,p,t,x)
    _vt = QuickStates.vtx()
    _g(_v) = mol_gibbs_impl(_vt,model,_v,t,x)
    _p(_v) = pressure_impl(_vt,model,_v,t,x)

    v1,v2 = v_zeros(mt,model,p,t,x)
    #(v1 == v2) && return v1
    g1 = _g(v1)
    g2 = _g(v2)
    if g1 < g2
        return v1
    elseif g2 < g1
        return v2
    else #different volumes, but equal energies, equilibria conditions
        return v2
    end
end

function v_zero_general(mt::MultiPT,model::CubicModel,p,t,x)
    _vt = QuickStates.vtx()
    _g(_v) = αR_impl(_vt,model,_v,t,x) - p*_v/(RGAS*t)
    _p(_v) = pressure_impl(_vt,model,_v,t,x)

    v1,v2 = v_zeros(mt,model,p,t,x)
    #(v1 == v2) && return v1
    g1 = _g(v1)
    g2 = _g(v2)
    if g1 < g2
        return v1
    elseif g2 < g1
        return v2
    else #different volumes, but equal energies, equilibria conditions
        return v2
    end
end

function roots_fzero(p,pspec,t,v)
    f0(x) = p(x) - pspec
    return Roots.find_zero(f0,v)
end

#Sucessive van der waals approximations, patent pending (lol)
#aproximates locally the function to a van der waals
#gas_fzero aproximates supposing gas phase (a-aproximation)
#liquid_fzero aproximates supposing liquid phase (a-aproximation)

function gas_fzero(p,Pspec,t,v0 = 40 * t / Pspec;
    relax=0.9*one(float(Pspec)),
    maxevals=100,
    atol = zero(float(Pspec)),
    rtol = 8*eps(one(float(Pspec))))

    _Pspec,_v0,_t,_atol,_rtol,_relax = promote(Pspec,v0,t,atol,rtol,relax)
    return _gas_fzero(p,_Pspec,_v0,_t,_relax,maxevals,_atol,_rtol)
end

function _gas_fzero(P,Pspec::TT,v0::TT,T::TT,relax = TT(0.9),maxevals=100, atol=zero(TT), rtol=8*eps(TT)) where TT
    pspec = Pspec
    R = RGAS #constant defined outside, 8.314...
    v₀ = v0
    P₀= P(v0)

    λ = relax
    _1 = one(TT)
    _0 = zero(TT)
    P₁ = _0
    v₁ = _0

    nan = _0/_0
    Pold = nan
    vold = nan

    RT = R*T
    RT_P = RT/Pspec
    invRT = _1/RT
    preA = Pspec*invRT*invRT
    preB = Pspec*invRT
    count = 0
    uatol = atol / oneunit(atol) * oneunit(real(v₁))
    adjustunit = oneunit(real(P₀))/oneunit(real(v₀))

    while count < maxevals
        #at high volumes, P ≈ RT/v - a/v2
        a = v₀*(RT-P₀*v₀)
        Δ = RT_P^2 - 4*a/Pspec
        if Δ > 0
            v₁ = 0.5*(RT_P + sqrt(Δ))
            P₁ = P(v₁)
        else
            if !isnan(vold)
            v₁ = 0.5*(v₀+vold)
            P₁ = P(v₁)
            else
                v₁ = 20v₀
                P₁ = P(v₁)
            end
        end

        if (P₁ > Pspec)
            vmax = v₁
            Pmax = P₁
        end
        if (Pspec < P₀ < P₁) #step in the wrong direction, reduce volume directly
            v₁ = v₀*1.1
            P₁ = P(v₁)
        end
        ΔP₀ = Pspec - P₀
        ΔP = Pspec - P₁
        ΔPold = Pspec - Pold

        #@show ΔP , ΔP₀, ΔPold,v₁
        abs(ΔP) <= eps(Pspec) && return v₁
        abs(ΔP) <= adjustunit * max(uatol, abs(v₁) * rtol) && return v₁
        iszero(ΔP) && return v₁
        isnan(ΔP) && return nan
        if (ΔP == ΔP₀)
            println("weird, v₁ =",v₁," and v₀ =",v₀)
            sdp = sign(ΔP)
            v_sprev = sign(P(prevfloat(v₁))-Pspec)
            v_snext = sign(P(prevfloat(v₁))-Pspec)
            p_sprev = sign(prevfloat(ΔP))
            p_snext = sign(nextfloat(ΔP))
            sdp * v_sprev <= 0 && return v₁
            sdp * v_snext <= 0 && return v₁
            sdp * p_sprev <= 0 && return v₁
            sdp * p_snext <= 0 && return v₁
            return nan
        end
        Pold,vold,P₀,v₀ = P₀,v₀,P₁,v₁
        count +=1
    end
    return nan
end

#=
Z = pv/rt
ZRT/p

=#

function liquid_fzero(p,Pspec,t,v0 = 0.375*RGAS*t/Pspec
    ,relax=0.9*one(float(Pspec))
    ,maxevals=100
    ,atol = zero(float(Pspec))
    ,rtol = 8*eps(one(float(Pspec))))
    _Pspec,_v0,_t,_atol,_rtol,_relax = promote(Pspec,v0,t,atol,rtol,relax)
    return _liquid_fzero(p,_Pspec,_v0,_t,_relax,maxevals,_atol,_rtol)
end

#2.2V = b

function _liquid_fzero(P,Pspec::TT,v0::TT,T::TT,relax = TT(0.9),maxevals=100, atol=zero(TT), rtol=8*eps(TT)) where TT
    _1 = one(TT)
    v = v0
    count = 0
    maxcount = 10
    converged = false
    vmin = TT(Inf)
    vold=TT(NaN)
    pold = TT(NaN)
    function f0(v)
        res = _1/P(v) - _1/Pspec
        if res < 0
           return -10*res
        else
           return res
        end
    end

    f00(v) =  P(v) - Pspec
    while !converged
        v0 = Optim.optimize(f0,eps(TT),v0,iterations=16).minimizer
        p0  = P(v0)
        if abs(vold-v0)/v0 < 0.001
            break
        end
        b1 = abs(-RGAS*T/p0 + v0)
        b2 = abs(-RGAS*T/Pspec + b1)
        vold = v0
        vmin = min(vmin,b1,b2,v0)
        v0 = vmin
        #@show v0
        if count >= maxcount
            converged = true
        end

        count +=1
        pold = p0

    end
    p0 = P(v0)
    m = ForwardDiff.derivative(P,v0)
    #m = y2-y1/(x2-x1)
    #(yx-y0) = m(xx-x0)
    vx = (Pspec - p0)/m + v0
    res =  roots_fzero(P,Pspec,T,vx)
    #@show res
    return res
end

#=
p > RT/v-b
p(v-b) > RT
v - b > RT/P
v + b> v > RT/P + b
v> v - b> RT/P
b< v-RT/P

=#

function v_zeros(::CubicRoots,mt::MultiPT,model,p,t,x;phase=:unspecified)
    _poly = cubic_poly(mt,model,p,t,x)
    poly = Polynomials.ImmutablePolynomial(_poly)
    sols = Polynomials.roots(poly)
    real_sols = cardan_reals(sols)
    RTp = RGAS*t/p
    vl = minimum(real_sols)*RTp
    vv = maximum(real_sols)*RTp
    if isnan(vl) & !isnan(vv)
        return (vv,vv)
    elseif isnan(vv) & !isnan(vl)
        return (vl,vl)
    else
        return (vl,vv)
    end
end

function v_zeros(::CubicRoots,mt::SinglePT,model,p,t;phase=:unspecified)
    _poly = cubic_poly(mt,model,p,t)
    poly = Polynomials.ImmutablePolynomial(_poly)
    sols = Polynomials.roots(poly)
    real_sols = cardan_reals(sols)
    RTp = RGAS*t/p
    vl = minimum(real_sols)*RTp
    vv = maximum(real_sols)*RTp
    if isnan(vl) & !isnan(vv)
        return (vv,vv)
    elseif isnan(vv) & !isnan(vl)
        return (vl,vl)
    else
        return (vl,vv)
    end
end

function v_zero(method::CubicRoots,mt::MultiPT,model,p,t,x,v0=nothing;phase=:unspecified)
    vols = v_zeros(method,mt,model,p,t,x)
    if is_liquid(phase)
        return first(vols)
    elseif is_gas(phase)
        return last(vols)
    else
        return v_zero_general(mt,model,p,t,x)
    end
end

function v_zero(method::CubicRoots,mt::SinglePT,model,p,t,v0=nothing;phase=:unspecified)
    vols = v_zeros(method,mt,model,p,t)
    if is_liquid(phase)
        return first(vols)
    elseif is_gas(phase)
        return last(vols)
    else
        return v_zero_general(mt,model,p,t)
    end
end


text:
# Create a shared C library

const CC = "gcc"
const LIB_NAME = "PDFHighlightsWrapper"
const OBJ_FLAGS = ["-O3", "-fPIC"]

if Sys.islinux()
    const LIB_EXT = ".so"
    const LIB_FLAGS = ["-shared", "-Wl,--no-undefined"]
elseif Sys.isapple()
    const LIB_EXT = ".dylib"
    const LIB_FLAGS = ["-shared"]
end

const LIB = joinpath(@__DIR__, LIB_NAME * LIB_EXT)

const INCLUDE_FLAGS = split(chop(read(`pkg-config --cflags poppler-glib`, String)))
const LINK_FLAGS = split(
    chop(
        "-lgio-2.0 " * read(`pkg-config --libs poppler-glib`, String)
    )
)

const OBJECTS = joinpath.(@__DIR__, ["get_lines_comments_pages", "get_author_title"])
const OBJECTS_C = OBJECTS .* ".c"
const OBJECTS_O = OBJECTS .* ".o"

if !isfile(LIB)
    for i in eachindex(OBJECTS)
        run(`$(CC) $(OBJ_FLAGS) -c $(OBJECTS_C[i]) -o $(OBJECTS_O[i]) $(INCLUDE_FLAGS)`)
    end
    run(`$(CC) $(LIB_FLAGS) -o $(LIB) $(OBJECTS_O) $(LINK_FLAGS)`)
end


text:
immutable NaiveSolver{N, T} <: MathProgBase.AbstractMathProgSolver
    rate::T
    max_step::T
    precondition_divisors::SVector{N, T}
    iteration_limit::Int
    gradient_convergence_tolerance::T
end

NaiveSolver(num_vars;
    rate=1.0,
    max_step=1.0,
    precondition_divisors::AbstractVector=ones(SVector{num_vars, Float64}),
    iteration_limit=100,
    gradient_convergence_tolerance=1e-3) =
NaiveSolver{num_vars, Float64}(rate,
    max_step,
    convert(SVector{num_vars, Float64}, precondition_divisors),
    iteration_limit,
    gradient_convergence_tolerance)

type NaiveModel{N, T} <: MathProgBase.AbstractNonlinearModel
    solver::NaiveSolver{N, T}
    x::MVector{N, T}
    x_lb::SVector{N, T}
    x_ub::SVector{N, T}
    g::MVector{N, T}
    objval::T
    problem::MathProgBase.AbstractNLPEvaluator
    status::Symbol

    NaiveModel(solver::NaiveSolver{N, T}) = new(solver)
end

NaiveModel{N, T}(solver::NaiveSolver{N, T}) = NaiveModel{N, T}(solver)

MathProgBase.NonlinearModel{N, T}(solver::NaiveSolver{N, T}) = NaiveModel{N, T}(solver)

function MathProgBase.loadproblem!{N, T}(model::NaiveModel{N, T},
    num_vars::Integer, num_constraints::Integer,
    x_lb::AbstractVector, x_ub::AbstractVector,
    g_lb::AbstractVector, g_ub::AbstractVector,
    sense::Symbol,
    problem::MathProgBase.AbstractNLPEvaluator)

    @assert sense == :Min
    @assert num_vars == N
    @assert num_constraints == 0

    model.x_lb = convert(SVector{N, T}, x_lb)
    model.x_ub = convert(SVector{N, T}, x_ub)
    model.x = zeros(MVector{N, T})
    model.g = zeros(MVector{N, T})
    model.objval = zero(T)
    model.problem = problem
end

MathProgBase.setwarmstart!(model::NaiveModel, x) = copy!(model.x, x)

function MathProgBase.optimize!{N, T}(model::NaiveModel{N, T})
    model.status = :UserLimit
    conditioned_max_step = model.solver.max_step .* model.solver.precondition_divisors
    conditioned_rate = model.solver.rate .* model.solver.precondition_divisors
    for i in 1:model.solver.iteration_limit
        MathProgBase.eval_grad_f(model.problem, model.g, model.x)

        if maximum(abs.(model.g)) < model.solver.gradient_convergence_tolerance
            model.status = :Optimal
            break
        end

        step = max.(min.(conditioned_max_step, conditioned_rate .* model.g), -conditioned_max_step)
        model.x .= max.(min.(model.x .- step, model.x_ub), model.x_lb)

    end
end

MathProgBase.status(model::NaiveModel) = model.status
MathProgBase.getsolution(model::NaiveModel) = model.x


text:
# This is an implementation of rigid body dynamics
# RigidBody is a mutable structure which only contains data.
# We have functions that work on this data structure, updating the values of the member variables.
# Need to define springs, dampers, constraints.

using LinearAlgebra
# include("OrientationConversion.jl");

mutable struct extForces
    Forces::Matrix{T} where T<:Real# size n x 3 (n Forces)
    Positions::Matrix{T} where T<:Real # size n x 3 (same size as Forces)
    Torques::Matrix{T} where T<:Real # size m x 3 (m Torques)
end

mutable struct RigidBody
    m::T where T<:Real; # Mass of the rigid body.
    In::Matrix{T} where T<:Real; # Moment of inertia 3x3 matrix
    orientation::String; # "321", "quaternions". Only these two supported now.
    # This will determine how many states there are
    # 12 for euler angle, 13 for quaternions.
    bodyID::Integer # Body number in robot tree

    # Internal variables
    x::Vector{T} where T<:Real;    # State vector
    dcm::Matrix{T} where T<:Real;  # Inertial to body.
    invI::Matrix{T} where T<:Real; # inverse of moment of inertia, compute once during initialization.
    J::Matrix{T} where T<:Real # 4x4 Inertia with first elemtn a random positive number
    invJ::Matrix{T} where T<:Real # inverse of 4x4 inertia matrix
    ω::Vector{T} where T<:Real # Angular velocity in body frame
    function RigidBody(mass::T, Inertia::Matrix{T}, bodyID::Integer; orientation::String = "quaternions") where T<:Real
        # if orientation == "321"
        #     x0 = Vector{Float64}(undef,12); # 12 state vector
        # elseif lowercase(orientation) == "quaternions"
        x0 = Vector{T}(undef,14); # 14 state vector
        # else
        #     error("Unknown orientation specified")
        # end

        dcm = Matrix{T}(undef,3,3);

        if det(Inertia) != 0
            invI = inv(Inertia)
        else
            invI = Matrix{T}(I,3,3)
        end

        J = [rand(1) zeros(1,3)
             zeros(3,1) Inertia]
        invJ = [1/J[1,1] zeros(1,3)
                zeros(3,1) invI]
        xdot = function times2(in::T)::T
            return(2*in);
        end
        this = new();
        this.m = mass;
        this.In = Inertia;
        this.bodyID = bodyID
        this.orientation = lowercase(orientation);
        this.dcm = dcm
        this.x = x0;
        this.invI = invI;
        this.J = J
        this.invJ = invJ
        this.ω = Vector{T}(undef,3)
        return this;
    end
end

# mutable struct RigidBody
#     m::Real; # Mass of the rigid body.
#     In::Matrix{Real}# where T<:Real; # Moment of inertia 3x3 matrix
#     orientation::String; # "321", "quaternions". Only these two supported now.
#     # This will determine how many states there are
#     # 12 for euler angle, 13 for quaternions.
#     bodyID::Integer # Body number in robot tree
#
#     # Internal variables
#     x::Vector{Real}# where T<:Real;    # State vector
#     dcm::Matrix{Real}# where T<:Real;  # Inertial to body.
#     invI::Matrix{Real}# where T<:Real; # inverse of moment of inertia, compute once during initialization.
#     J::Matrix{Real}# where T<:Real # 4x4 Inertia with first elemtn a random positive number
#     invJ::Matrix{Real}# where T<:Real # inverse of 4x4 inertia matrix
#     ω::Vector{Real}# where T<:Real # Angular velocity in body frame
#     function RigidBody(mass::Real, Inertia::Matrix{Real}, bodyID::Integer; orientation::String = "quaternions") #where T<:Real
#         # if orientation == "321"
#         #     x0 = Vector{Float64}(undef,12); # 12 state vector
#         # elseif lowercase(orientation) == "quaternions"
#         x0 = Vector{Real}(undef,14); # 14 state vector
#         # else
#         #     error("Unknown orientation specified")
#         # end
#
#         dcm = Matrix{Real}(undef,3,3);
#
#         if det(Inertia) != 0
#             invI = inv(Inertia)
#         else
#             invI = Matrix{Real}(I,3,3)
#         end
#
#         J = [rand(1) zeros(Real,(1,3))
#              zeros(Real,(3,1)) Inertia]
#         invJ = [1/J[1,1] zeros(Real,(1,3))
#                 zeros(Real,(3,1)) invI]
#         xdot = function times2(in::Real)::Real
#             return(2*in);
#         end
#         this = new();
#         this.m = mass;
#         this.In = Inertia;
#         this.bodyID = bodyID
#         this.orientation = lowercase(orientation);
#         this.dcm = dcm
#         this.x = x0;
#         this.invI = invI;
#         this.J = J
#         this.invJ = invJ
#         this.ω = Vector{Real}(undef,3)
#         return this;
#     end
# end

function rbDynQuat(RB::RigidBody, unconstrainedF::Vector{Float64}, Fc::Vector{Float64})::Vector{Float64}

    # States are: 14 states
    # x = [x,y,z]:inertial, [u,v,w]:inertial, [β0,β1,β2,β3], [β̇0,β̇1,β̇2,β̇3]
    # unconstrainedF is a 7x1 vector, unique to each rigid body.
    # Fc is the  generalised force coming from constraints acting on the rigid body. Unique to each body. 7x1 vector.

    xdot = Vector{Float64}(undef,14)

    u = RB.x[8:10];   # Velocities in Inertial reference frame
    β = RB.x[4:7];  # Euler parameters
    βdot = RB.x[11:14] # Euler parameter derivatives

    TotalForce = unconstrainedF[1:3] + Fc[1:3]
    TotalMoment = unconstrainedF[4:7] + Fc[4:7]

    # # Udwadia's formulation using 14 States
    E = genE(β)
    xdot[1:3] = u
    xdot[8:10] = TotalForce/RB.m
    xdot[4:7] = βdot
    xdot[11:14] = 1/4*transpose(E)*RB.invJ*E*TotalMoment
    return xdot
end

function initialiseRigidBody!(RB::RigidBody,x0::Vector{T})::RigidBody where T<:Real
    # Initialise rigid body with x0
    RB.x = x0
    return RB
end

function updateRigidBody!(b::RigidBody,x::Vector{T})::RigidBody where T<:Real
    b.x = x
    b.dcm = quat2dcm(x[4:7])
    b.ω = angVel(x[4:7],x[11:14])
    return b
end

function InertialFrameAsRB()::RigidBody
    m = 0.0
    I = zeros(3,3)

    b = RigidBody(m,I,1)
    b.x = [zeros(3);1.0;zeros(3);zeros(3);zeros(4)]
    b.dcm = quat2dcm(b.x[4:7])
    b.ω = angVel(b.x[4:7],b.x[11:14])
    return b
end

function zeroExtForce()::extForces
    # External forces function to be changed later.
    Forces = zeros(1,3)
    Positions = zeros(1,3)
    Torques = zeros(1,3)
    return extF = extForces(Forces,Positions,Torques)
end

function zeroExtForceVec(l::Int64)::Vector{extForces}
    # l is the length of the vector
    v = Vector{extForces}(undef,l)
    for i=1:l
        v[i] = zeroExtForce()
    end
    return v
end

## Euler angle formulation (12 states)
# function rbDynEuler(xdot::Vector{Float64}, RB::RigidBody,
#     PositionList::Matrix{Float64}, ForceList::Matrix{Float64},
#     TorqueList::Matrix{Float64}, GravityInInertial::Vector{Float64})
#     # States are: 12 states
#     # x = [x,y,z]:inertial, [u,v,w]:body, [ϕ,θ,ψ], [ω1,ω2,ω3]:body
#     # Forces in ForceList are decribed in body reference, ForceList is nFx3.
#     # PositionList are positions where forces are acting relative to cm,
#     # decribed in body reference, PositionList is nFx3.
#     # Torques in TorqueList are in body reference. TorqueList is matrix nTx3
#
#     xpos = RB.x[1];  # x position in inertial frame
#     ypos = RB.x[2];  # y position in inertial frame
#     h    = RB.x[3];  # z position in inertial frame
#
#     u = RB.x[4:6];   # Velocities in body reference frame
#     ang = RB.x[7:9];     # Euler parameters
#     ω = RB.x[10:12]; # Angular velocity  in body reference frame
#
#     # DCM w.r.t Euler parameters. vb = C*vn, inertial to body.
#     ang2dcm(RB.dcm,ang,RB.orientation); # Have to code this up.
#
#     GravityInBody = RB.dcm*GravityInInertial;  # Convert to body reference
#
#     xdot[1:3] = RB.dcm'*u; # Velocities in inertial frame.
#     xdot[4:6] = sum(ForceList,dims=1)/RB.m + GravityInBody; # Euler's first law.
#     xdot[7:9] = # Code this
#
#     # ω dot equation -- Euler's second law
#     TotalMoment = sum(cross(PositionList[i,:],ForceList[i,:]) for i in 1:size(Fi)[1]) + sum(TorqueList,dims=1);
#     xdot[10:12] = RB.invI*(TotalMoment - cross(ω,RB.I*ω));
# end

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
# Warning: setupLLCgrid no longer exists
# ggebbie, 16-Apr-2021

using Revise
using MITgcmTools, MeshArrays, Statistics, ECCOtour, PyPlot

# get MITgcm / ECCOv4r4 LLC grid and depth information. Store in γ.
path_grid="../inputs/GRID_LLC90/"
γ = setupLLCgrid(path_grid)

# list of experiments on poseidon
exppath = "/batou/ECCOv4r4/MITgcm/exps/"
runpath,diagpath = listexperiments(exppath);

# print output here
path_out = "/home/gebbie/julia/outputs/"

# abbreviations for each experiment for labels, etc.
shortnames = expnames()

## SELECT EXPERIMENTS TO ANALYZE ##############
# manually choose from available experiments listed above.
exps = ("iter129_bulkformula","nointerannual")

# to do all experiments:
# exps = keys(shortnames)
###############################################

nexps = length(exps) # number of experiments

## DEFINE THE LIST OF SIGMA1 VALUES.
sig1grida = 24:0.05:31
sig1gridb = 31.02:0.02:33
sig1grid = vcat(sig1grida,sig1gridb)
sig1grid = sig1grid[1:3:end]
nσ = length(sig1grid)

# get index for sigma1 value of interest.
sig1 = 29.9 # a choice
latsig1 = -12.0 # lat of interest

junk,isig1grid = findmin(abs.(sig1grid.-sig1))

tecco= 1992+1/24:1/12:2018 # ecco years

ϕGreg,ϕCreg = latgridRegular(γ)
junk,ilat = findmin(abs.(ϕCreg.-latsig1))
λC = -179.5:179.5
# center the plot on indo-pac
cutoff = 181
λCcutoff = cat(λC[cutoff:end],λC[1:cutoff-1],dims = 1)

# make all longitudes positive
replace!(x -> x<0 ? x+=360 : x,λCcutoff)

nx = length(λC)
TSroot = "S_on_sigma1"
cmap_seismic =get_cmap("seismic")
figure(100)
clf()
plot(λCcutoff)

# cycle through all chosen experiments
#ex = exps[1]; for manual runs
for ex in exps

    # Get list of files for salinity on sigma1
    filelist = searchdir(runpath[ex]*"sigma1/",TSroot) # first filter for state_3d_set1
    datafilelist  = filter(x -> occursin("data",x),filelist) # second filter for "data"
    nt = length(datafilelist)
    Ssig1slice = Array{Float32,2}(undef,nx,nt)

    # Find a latitudinal slice. i.e., 12 South. Use convert2array? Or something else? Use regular poles method.

    # pre-allocate
    #TS = Array{Float32, 2}(undef, nt, nz*2)

    global tt = 0
    for Sname in datafilelist
        tt += 1
        println("year ",Int(floor(tecco[tt]))," month ",((tt-1)%12)+1)

        # get S on sigma1. Way to read a slice? (Didn't figure it out yet)
        @time S = γ.read(runpath[ex]*"sigma1/"*Sname,MeshArray(γ,Float32,nσ))

        # just take one sigma1 surface
        Ssig1 = S[:,isig1grid]

        # next slice at 12 South.
        Ssig1crop =  LLCcropC(Ssig1,γ) # regular grid by cropping
        Ssig1slice[:,tt] = Ssig1crop[:,ilat]
    end
    #Ssig1slice = Ssig1slice[cutoff:end 1:cutoff-1,:]
    Ssig1slice = cat(Ssig1slice[cutoff:end,:],Ssig1slice[1:cutoff-1,:],dims = 1)

    figure(10)
    clf()
    iplot = 140:290; # plot indo-pac only
    lims = 35.0:.05:36.0
    contourf(λCcutoff[iplot],tecco,Ssig1slice[iplot,:]',lims,cmap=cmap_seismic)
    colorbar(label="[PSS-78]",orientation="vertical",ticks=lims)
    contour(λCcutoff[iplot],tecco,Ssig1slice[iplot,:]',lims,colors="k")

    titlelbl = ex*", "*string(latsig1)*L"\degree"*"N, "*L"\sigma_1="*string(sig1)*" kg m"*L"^{-3}"
    xlbl = "longitude "*L"[\degree E]"
    ylbl = "years"

    title(titlelbl)
    xlabel(xlbl)
    ylabel(ylbl)

    # Put output (i.e., figures) locally.
    pathoutexp = path_out*ex
    !isdir(pathoutexp) ? mkdir(pathoutexp) : nothing;
    if latsig1 > 0
        outfname = pathoutexp*"/Sonsigma1_"*shortnames[ex]*"_"*string(latsig1)*"N_"*string(sig1)*".eps"
    else
        outfname = pathoutexp*"/Sonsigma1_"*shortnames[ex]*"_"*string(abs(latsig1))*"S_"*string(sig1)*"sig1.eps"
    end

    savefig(outfname)
end