text:
###################################################################################################
# CONSTANTS
###################################################################################################
const BG_BUFFER = 50000                     # Bytes of bedGraph records until write
const GFF_BUFFER = 500000                   # Bytes of GFF records until write
const THRESH_MAPQ = 39                      # MAPQ threshold (-10*log(p)) only true uni-reads
const FLAGS_ALLOWED = [0,16,83,99,147,163]  # Flags allowed in BAM recs
###################################################################################################
# STRUCTS
###################################################################################################
struct AlignTemp
    strand::String      # Methylation call strand
    R1::BAM.Record      # First record from left to right
    R2::BAM.Record      # Second record from left to right
end
mutable struct AllAlignTemps
    paired::Bool
    templates::Array{AlignTemp,1}
end
###################################################################################################
# FUNCTIONS
###################################################################################################
"""
    `get_align_strand(PAIRED_END,FLAG1,FLAG2)`

Function that returns the strand of methylation call based on Bismark's logic. In single-end mode,
OT and CTOT reads will both receive a FLAG of 0 while OB and CTOB both receive a FLAG of 16. In
paired end mode:

             Read 1       Read 2
  OT:         99          147
  OB:         83          163
  CTOT:      147           99
  CTOB:      163           83

This table was extracted from

    https://github.com/FelixKrueger/Bismark/issues/151

# Examples
```julia-repl
julia> CpelAsm.get_align_strand(true,UInt16(99),UInt16(147))
"OT"
```
"""
function get_align_strand(pe::Bool,flag1::UInt16,flag2::UInt16)::String

    # Treat SE and PE separately
    if !pe
        if flag1==0
            s="OT"
        elseif flag1==16
            s="OB"
        else
            println(stderr,"[$(now())]: SE: was expecting FLAG 0 or and")
            println(stderr,"[$(now())]: encountered $(flag1) instead.")
            println(stderr,"[$(now())]: Exiting julia ...")
            exit(1)
        end
    else
        if flag1==99 && flag2==147
            s="OT"
        elseif flag1==83 && flag2==163
            s="OB"
        elseif flag1==147 && flag2==99
            s="CTOT"
        elseif flag1==163 && flag2==83
            s="CTOB"
        else
            println(stderr, "[$(now())]: PE: unexpected flag combination.")
            println(stderr, "[$(now())]: Expected 99/147, 147/99, 83/163, 163/83.")
            println(stderr, "[$(now())]: Exiting julia ...")
            exit(1)
        end
    end

    # Return strand
    return s

end # end get_align_strand
"""
    `order_bams(PAIRED_END,RECORDS)`

Function that returns an AlignTemp object with R1 as the first record in the forward strand and
with the methylation call strand taken from `get_align_strand`.

# Examples
```julia-repl
julia> CpelAsm.order_bams(true,RECORDS)
```
"""
function order_bams(pe::Bool,records::Vector{BAM.Record})::AlignTemp

    # Check which record appears first in forward strand
    if BAM.position(records[1])<=BAM.position(records[2])
        s = get_align_strand(pe, BAM.flag(records[1]), BAM.flag(records[2]))
        return AlignTemp(s,records[1],records[2])
    else
        s = get_align_strand(pe, BAM.flag(records[2]), BAM.flag(records[1]))
        return AlignTemp(s,records[2],records[1])
    end

end
"""
    `clean_records(PAIRED_END,RECORDS)`

Function that takes a set of records and returns an AllAlignTemps object that contains all the
properly aligned templates as an array of AlignTemp, which contains information about the
methylation call strand as well as the relevant BAM records. In the PE case, R1 corresponds to
the BAM record that appears before in the forward strand.

# Examples
```julia-repl
julia> CpelAsm.clean_records(true,RECORDS)
```
"""
function clean_records(pe::Bool,records::Vector{BAM.Record})::AllAlignTemps

    # Initialize struct
    out = AllAlignTemps(pe,[])

    # Consider PE vs. SE
    if pe
        # Handle PE case, first retrieve unique template names
        temp_names = unique([BAM.tempname(x) for x in records])
        for temp_name in temp_names
            # Get records with same template name
            temp_recs = filter(x->BAM.tempname(x)==temp_name,records)

            # There should only be two records with the same template name
            length(temp_recs)==2 && push!(out.templates,order_bams(pe,temp_recs))
        end
    else
        # Handle SE case
        out.templates = [AlignTemp(get_align_strand(pe,BAM.flag(x),UInt16(0)),x,BAM.Record()) for
            x in records]
    end

    # Return struct
    return out

end # end clean_records
"""
    `try_olaps(READER,CHR,WINDOW)`

Function that tries to find BAM records overlaping with `CHR` at positions `WINDOW`.

# Examples
```julia-repl
julia> CpelAsm.try_olaps(reader,chr,win)
```
"""
function try_olaps(reader::BAM.Reader,chr::String,win::Vector{Int64})::Vector{BAM.Record}

    # NOTE: known bug that has to do w/ BioAlignments when close to end?
    records_olap =
    try
        collect(eachoverlap(reader, chr, win[1]:win[2]))
    catch
        Array{BAM.Record,1}()
    end

    return records_olap

end # try_olaps
"""
    `read_bam(BAM_PATH,CHR,FEAT_ST,FEAT_END,CPG_POS,CHR_SIZE,PE,TRIM)`

Function that reads in BAM file in `BAM_PATH` and returns methylation vectors for those records
that overlap with (1-based) genomic coordinates `chr:chrSt-chrEnd` at `cpg_pos`. A trimming given
by TRIM=(Rf_5p,Rf_3p,Rr_5p,Rr_3p) is applied to the reads. The information was taken from:

  XM: meth calls (https://github.com/FelixKrueger/Bismark/tree/master/Docs)
  XS: uniqueness (http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml)

For info on OT, CTOT, OB, CTOB nomenclature see

    http://software.broadinstitute.org/software/igv/book/export/html/37.

# Examples
```julia-repl
julia> CpelAsm.read_bam(BAM_PATH,"chr1",30,80,[40,60],1000,false,(0,0,0,0))
```
"""
function read_bam(bam::String,chr::String,hap_st::Int64,hap_end::Int64,cpg_pos::Vector{Int64},
                  chr_size::Int64,pe::Bool,trim::NTuple{4,Int64})::Array{Vector{Int64},1}

    # Number of CpG sites is determined by that in the region
    N = length(cpg_pos)

    # Expand window in PE case to include pair even if it's outside the window
    exp_win = pe ? [max(1,hap_st-75),min(chr_size,hap_end+75)] : [hap_st,hap_end]

    # Get records overlapping window.
    reader = open(BAM.Reader,bam,index=bam*".bai")
    records_olap = try_olaps(reader,chr,exp_win)
    length(records_olap)>0 || return Vector{Int64}[]
    close(reader)

    # Relevant flags in BAM file (both Bismark and Arioc)
    filter!(x-> (BAM.ismapped(x)) && (haskey(x,"XM")) && (!haskey(x,"XS")) && (BAM.flag(x) in
           FLAGS_ALLOWED) && (BAM.mappingquality(x)>THRESH_MAPQ),records_olap)

    # Clean records & organize them
    records_org = clean_records(pe,records_olap)

    # Loop over organized records
    xobs = Vector{Int64}[]
    for record in records_org.templates
        # Get methylation call and offset (depends on strand where call is made)
        if !(records_org.paired)
            # Obtain methylation call from single end ()
            meth_call = record.R1[:"XM"]
            meth_call = SubString(meth_call,(1+trim[1]):(length(meth_call)-trim[2]))
            OFFSET = record.strand in ["OT","CTOT"] ? 1 : 2
            OFFSET -= BAM.position(record.R1)+trim[1]
        else
            # Obtain methylation call
            R1_call = record.R1[:"XM"]
            R2_call = record.R2[:"XM"]
            R1_call = SubString(R1_call,(1+trim[1]):(length(R1_call)-trim[2]))
            R2_call = SubString(R2_call,(1+trim[4]):(length(R2_call)-trim[3]))
            dist_st = abs(BAM.position(record.R2)+trim[4] - BAM.position(record.R1)-trim[1])
            meth_call = R1_call * "."^max(0,dist_st-length(R1_call))
            meth_call *= R2_call[max(1,(length(R1_call)+1-dist_st)):end]
            OFFSET = record.strand in ["OT","CTOT"] ? 1 : 2
            OFFSET -= BAM.position(record.R1)+trim[1]
        end

        # Cross positions of CpG sites if template contains CpG sites
        obs_cpgs = [x.offset for x in eachmatch(r"[zZ]",meth_call)] .- OFFSET
        length(obs_cpgs)>0 || continue

        # If overlapping CpG sites then store and add to xobs
        olap_cpgs = findall(x-> x in obs_cpgs, cpg_pos)
        length(olap_cpgs)>0 || continue
        x = zeros(Int64,N)
        obs_cpgs = obs_cpgs[findall(x-> x in cpg_pos,obs_cpgs)] .+ OFFSET
        x[olap_cpgs] = reduce(replace,["Z"=>1,"z"=>-1],init=split(meth_call[obs_cpgs],""))

        # Add to set of observations
        push!(xobs,x)

    end # end loop over templates sequenced

    # Return
    return xobs

end # end read_bam
"""
    `write_gff!(OUT_GFF,GFF_RECORDS)`

Function that appends `GFF_RECORDS` into `OUT_GFF`.

# Examples
```julia-repl
julia> CpelAsm.write_gff!(gff, gff_records)
```
"""
function write_gff!(gff::String,gff_records::Vector{GFF3.Record})

    # Create folder if it doesn't exist
    isdir(dirname(gff)) || mkdir(dirname(gff))

    # Append gff records
    output = open(gff, "a")
    writer = GFF3.Writer(output)
    while length(gff_records)>0
        write(writer,popfirst!(gff_records))
    end
    close(writer)

    # Return
    return nothing

end # end write_gff
"""
    `next_record(READER,CHR_NAMES,RECORD)`

Recursive function that returns next record in stream if in chr_names, and false if the stream
reached its end. This function is only used in `gen_gffs()`.

# Examples
```julia-repl
julia> CpelAsm.next_record(reader, chr_names, record)
```
"""
function next_record(reader::VCF.Reader,chr_names::Vector{String},record::VCF.Record)::VCF.Record

    # Check if record in chr_names and end of reader
    if VCF.chrom(record) in chr_names
        # If we're good, return record
        return record
    elseif !eof(reader)
        # If not end of file, call function again
        read!(reader, record)
        next_record(reader, chr_names, record)
    else
        # return empty record if finished stream
        return VCF.Record()
    end

end # end next_record
"""
    `is_het_cpg!(VAR,SEQ,H1,H2)`

Function that adds position of heterozygous CpG site, if existant, to the correct haplotype.
Assumes that REF is the one pertaining to genome 1, and ALT is the one pertaining to genome 2.

# Examples
```julia-repl
julia> CpelAsm.is_het_cpg!(var,seq,h1,h2)
```
"""
function is_het_cpg!(var::VCF.Record,seq::FASTA.Record,h1::Vector{Int64},h2::Vector{Int64})

    # Initialize
    ref_var = join(VCF.ref(var))
    alt_var = join(VCF.alt(var))

    # Check if heterozygous CpG site in XG context
    if ((ref_var=="C") && (alt_var in ["A","G","T"]))
        # If followed by G, then heterozygous CpG in haplotype 1
        if FASTA.sequence(seq,VCF.pos(var):(VCF.pos(var)+1))[2]==DNA_G
            push!(h1,VCF.pos(var))
        end
    elseif ((alt_var=="C") && (ref_var in ["A","G","T"]))
        # If followed by G, then heterozygous CpG in haplotype 2
        if FASTA.sequence(seq,VCF.pos(var):(VCF.pos(var)+1))[2]==DNA_G
            push!(h2,VCF.pos(var))
        end
    end

    # Check if het CpG site in CX context
    if ((ref_var=="G") && (alt_var in ["A","C","T"]))
        # If preceeded by C, then heterozygous CpG in haplotype 1
        if FASTA.sequence(seq,(VCF.pos(var)-1):VCF.pos(var))[1]==DNA_C
            push!(h1,VCF.pos(var)-1)
        end
    elseif ((alt_var=="G") && (ref_var in ["A","C","T"]))
        # If preceeded by C, then heterozygous CpG in haplotype 2
        if FASTA.sequence(seq,(VCF.pos(var)-1):VCF.pos(var))[1]==DNA_C
            push!(h2,VCF.pos(var)-1)
        end
    end

    # Return
    return nothing

end # end is_het_cpg!
"""
    `get_records_ps!(VCF_READER,SEQ,RECORD,CURR_PS,WIN_END,H1,H2)`

Recursive function that returns the end of current PS phased SNPs along the next SNP. This
function is only used in `gen_gffs()`.

# Examples
```julia-repl
julia> CpelAsm.get_records_ps!(reader,seq,record,curr_ps,wEnd,h1,h2)
```
"""
function get_records_ps!(reader::VCF.Reader,seq::FASTA.Record,record::VCF.Record,curr_ps::String,
                         wEnd::Int64,h1::Vector{Int64},h2::Vector{Int64})::Tuple{VCF.Record,Int64}

    # Check if record has same PS
    ind_PS = findfirst(x->x=="PS",VCF.format(record))
    VCF.genotype(record)[1][ind_PS]==curr_ps || return record,wEnd

    # Check GT string and update het. CpG sites
    ind_GT = findfirst(x->x=="GT",VCF.format(record))
    curr_gt = VCF.genotype(record)[1][ind_GT]
    curr_gt in ["0/1","0|1"] ? is_het_cpg!(record,seq,h1,h2) : is_het_cpg!(record,seq,h2,h1)

    # Continue reading
    wEnd = VCF.pos(record)
    if !eof(reader)
        read!(reader, record)
        if "PS" in VCF.format(record)
            record,wEnd = get_records_ps!(reader,seq,record,curr_ps,wEnd,h1,h2)
        else
            return record,wEnd
        end
    else
        # return empty record if finished stream
        return VCF.Record(),wEnd
    end

end # end get_records_ps
"""
    `expand_win(FIRST_SNP,LAST_SNP,WIN_EXP,CHR_SIZE)`

Function that returns an expanded window given the position of the first and last SNPs, the
expansion to be done to the left and right, of the first and last SNP respectively, to include
nearby CpG sites, and the size of the subregions into which the whole region will be divided into.

# Examples
```julia-repl
julia> CpelAsm.expand_win(50,80,10,2000)
2-element Array{Int64,1}:
 40
 90
```
"""
function expand_win(l_snp::Int64,r_snp::Int64,win_exp::Int64,chr_size::Int64)::Vector{Int64}

    # Return extended window
    return [max(1,l_snp-win_exp),min(chr_size,r_snp+win_exp)]

end # end expand_win
"""
    `find_num_models(CPG_POS,NMAX,NMOD)`

Function that returns the number of models required so as to avoid joint model with N>NMAX.

# Examples
```julia-repl
julia> CpelAsm.find_num_models([50,55,60,65],2,1)
2
```
"""
function find_num_models(cpg_pos::Vector{Int64},n_max::Int64,n_mod::Int64)::Int64

    # Check if resulting models are acceptable
    return length(cpg_pos)/n_mod>n_max ? find_num_models(cpg_pos,n_max,n_mod+1) : n_mod

end # end find_num_models
"""
    `gen_gffs([HET_GFF_PATH,HOM_GFF_PATH],FA_PATH,VCF_PATH,WIN_EXP,NMAX)`

Function that creates a GFF file containing the genomic regions to be analyzed by `CpelAsm`. If a
set of SNPs is phased, this function will create a single window for all. In case, the set of SNPs
is not phased, then an individual window will be created for each SNP. The phasing of the VCF
records should be specified by means of the standard `PS` tag. In addition, a GFF file containing
the homozygous portion of the DNA is also generated together with the position of the CpG sites
in it.

# Examples
```julia-repl
julia> CpelAsm.gen_gffs([HET_GFF_PATH,HOM_GFF_PATH],FA_PATH,VCF_PATH,WIN_EXP,NMAX)
```
"""
function gen_gffs(gff::Vector{String},fa::String,vcf::String,win_exp::Int64,n_max::Int64)

    # Initialize VCF variables
    het_records = Vector{GFF3.Record}()
    hom_records = Vector{GFF3.Record}()
    reader_vcf = open(VCF.Reader,vcf)
    record = VCF.Record()
    read!(reader_vcf,record)

    # Initialize FASTA variables
    prev_end = 0
    reader_fa = open(FASTA.Reader, fa, index=fa*".fai")
    chr_names = reader_fa.index.names
    curr_chr = chr_names[1]
    fa_record = reader_fa[curr_chr]
    close(reader_fa)

    # Chrom sizes
    chr_sizes = reader_fa.index.lengths
    chr_size = chr_sizes[findfirst(x->x==curr_chr,chr_names)]

    # Loop over variants
    while true

        # Check if chromosome is in FASTA file. If not, obtain the next valid 1
        record = next_record(reader_vcf,chr_names,record)

        # Check we have loaded the right chromosome from FASTA
        if !(curr_chr == VCF.chrom(record))
            curr_chr = VCF.chrom(record)
            reader_fa = open(FASTA.Reader,fa,index=fa*".fai")
            fa_record = reader_fa[curr_chr]
            chr_size = chr_sizes[findfirst(x->x==curr_chr,chr_names)]
            close(reader_fa)
            prev_end=0
        end

        # Get entire (contiguous) haplotype using the PS field
        wSt = VCF.pos(record)
        wEnd = VCF.pos(record)
        het1 = Vector{Int64}()
        het2 = Vector{Int64}()
        if "PS" in VCF.format(record)
            # If PS tag, then phased
            ind_PS = findfirst(x->x=="PS",VCF.format(record))
            curr_ps = VCF.genotype(record)[1][ind_PS]
            record,wEnd = get_records_ps!(reader_vcf,fa_record,record,curr_ps,wEnd,het1,het2)
        else
            # If no PS tag, then single SNP (marked as NOPS)
            curr_ps = "NOPS"
            is_het_cpg!(record,fa_record,het1,het2)
            eof(reader_vcf) ? record=VCF.Record() : read!(reader_vcf, record)
        end

        # Expand window win_exp on each side
        win = expand_win(wSt,wEnd,win_exp,chr_size)

        # Obtain DNA sequence from reference and homozygous CpG loci from it
        wSeq = convert(String,FASTA.sequence(fa_record,win[1]:win[2]))
        cpg_pos = map(x->getfield(x,:offset),eachmatch(r"CG",wSeq)).+win[1].-1

        # Store heterozygous haplotype & corresponding CpG sites if homozygous CpG site/s
        if length(cpg_pos)>0

            # Obtain required number of models based on homozygous CpG sites
            n_mod = find_num_models(cpg_pos,n_max,1)
            mod_size = div(length(cpg_pos),n_mod)
            mod_rem = length(cpg_pos)-n_mod*mod_size

            # Print each set separately
            for i=1:n_mod
                # Index of subset of CpG sites in cpg_pos
                hom_rng = (mod_size*(i-1)+1):min((mod_size*i),length(cpg_pos))
                # Add remainders to last block
                i==n_mod && mod_rem>0 && (hom_rng=extrema(hom_rng)[1]:length(cpg_pos))
                # Get subset of CpG sites for each
                hom_mod = cpg_pos[hom_rng]
                # Get submodel boundaries
                inf_bound = i==1 ? win[1] : cpg_pos[max(1,extrema(hom_rng)[1]-1)]
                sup_bound = i==n_mod ? win[2] : hom_mod[end]
                # Get heterozygous CpG sites in current submodel
                het1_mod = het1[findall(x-> inf_bound <= x <= sup_bound,het1)]
                het2_mod = het2[findall(x-> inf_bound <= x <= sup_bound,het2)]
                # Check there is no overlap between homozygous and heterozygous CpG sites
                if !isempty(intersect(union(het1_mod,het2_mod),hom_mod))
                    print_log("Found overlap between homozygous and heterozygous CpG sites.")
                    print_log("Check FASTA file is masked. Exiting julia ...")
                    exit(1)
                end
                # Make submodel window
                win_mod = expand_win(inf_bound,sup_bound,0,chr_size)
                # Push output string
                out_str = "$(curr_chr)\t.\t$(curr_ps)-$(i)/$(n_mod)\t$(win_mod[1])\t$(win_mod[2])"
                out_str *= "\t.\t.\t.\tN=$(length(hom_mod));CpGs=$(hom_mod)"
                length(het1_mod)>0 && (out_str*=";hetCpGg1=$(het1_mod)")
                length(het2_mod)>0 && (out_str*=";hetCpGg2=$(het2_mod)")
                # Push GFF record
                push!(het_records,GFF3.Record(out_str))
            end

        end

        # Check if het_records object is too big and empty it if so
        sizeof(het_records)>GFF_BUFFER && write_gff!(gff[1],het_records)

        # Obtain homozygous portion
        hom_win = [prev_end+1,win[1]-1]
        hom_win = [max(1,hom_win[1]),min(chr_size,hom_win[2])]

        # Obtain DNA sequence from reference and homozygous CpG loci from it
        wSeq = convert(String,FASTA.sequence(fa_record,hom_win[1]:hom_win[2]))
        cpg_pos = map(x->getfield(x,:offset),eachmatch(r"CG",wSeq)).+hom_win[1].-1

        # Store homozygous stretch of CpG sites found
        if length(cpg_pos)>0
            out_str = "$(curr_chr)\t.\t.\t$(hom_win[1])\t$(hom_win[2])\t$(length(cpg_pos))\t.\t.\t"
            out_str *= "CpGs=$(cpg_pos)"
            push!(hom_records,GFF3.Record(out_str))
        end

        # Update previous end
        prev_end = win[2]

        # Check if hom_records object is too big and empty it if so
        sizeof(hom_records)>GFF_BUFFER && write_gff!(gff[2],hom_records)

        # If new record empty break while loop
        VCF.haspos(record) || break

    end

    # Dump remaining records
    write_gff!(gff[1],het_records)
    write_gff!(gff[2],hom_records)

    # Return
    return nothing

end # end gen_gffs
"""
    `read_gff_chr(GFF_PATH,CHR)`

Function that reads in a GFF3 file in `GFF_PATH` and returns an array of GenomicFeatures Record
objects contained in chromosome `CHR`. The format is the standard defined by ENSEMBL in
https://useast.ensembl.org/info/website/upload/gff3.html.

# Examples
```julia-repl
julia> CpelAsm.read_gff_chr(GFF_PATH,"chr1")
```
"""
function read_gff_chr(gff::String,chr::String)::Vector{GFF3.Record}

    # Load genomic features from a GFF3 file
    features = open(collect,GFF3.Reader,gff)

    # Keep features in chr
    filter!(x -> GFF3.seqid(x)==chr, features)

    # Return
    return features

end # end read_gff_chr
"""
    `read_model_chr(MODEL_PATH,CHR)`

Function that reads in model file and returns dictionary with key the region ID and entry the model.

# Examples
```julia-repl
julia> CpelAsm.read_model_chr(model_path,"chr1")
```
"""
function read_model_chr(mod_path::String,chr::String)::Dict{String,Tuple{Vector{Float64},Vector{Int64}}}

    # Init output dictionary
    out_dic = Dict{String,Tuple{Vector{Float64},Vector{Int64}}}()

    # Go over file and add entries in chromosome
    open(mod_path) do f
        while !eof(f)

            # Read line
            x = readline(f)

            # Get data
            x = split(x,"\t")
            chr_reg = convert(String,x[1])

            # If not right chromosome continue
            chr_reg == chr || continue

            # Get rest of data
            chrst = convert(String,x[2])
            chrend = convert(String,x[3])
            θ = parse.(Float64,split(x[4],","))
            nvec = parse.(Int64,split(x[5],","))

            # Add entry to dictionary
            out_dic[join([chr_reg,chrst,chrend],"-")] = (θ,nvec)

        end
    end

    # Return dictionary
    return out_dic

end # end read_model_chr
"""
    `write_tobs(RECORDS,CHR,PATH)`

Function that writes records in `RECORDS` into `PATH`.

# Examples
```julia-repl
julia> CpelAsm.write_tobs(RECORDS,CHR,PATH)
```
"""
function write_tobs(recs::Vector{Tuple{Int64,Int64,Float64,Int64,Int64}},chr::String,path::String)

    # Open output bedGraph file in append mode (no need to close it)
    open(path, "a") do f
        for i=1:length(recs)
            recs[i][1]!=0 || continue
            write(f,"$(chr)\t"*join(string.(collect(recs[i])),"\t"),"\n")
        end
    end

    # Return
    return nothing

end # end write_tobs
"""
    `write_model(RECORDS,CHR,PATH)`

Function that writes θ parameter vectors in `RECORDS` into `PATH`.

# Examples
```julia-repl
julia> CpelAsm.write_model(RECORDS,CHR,PATH)
```
"""
function write_model(recs::Vector{Tuple{Int64,Int64,Vector{Float64},Vector{Int64}}},chr::String,path::String)::Nothing

    # Open output bedGraph file in append mode (no need to close it)
    open(path, "a") do f
        for i=1:length(recs)

            # Skip if no data
            recs[i][1]!=0 || continue

            # Get model
            θ = join(recs[i][3],",")
            nvec = join(recs[i][4],",")

            # Write line
            write(f,"$(chr)\t$(recs[i][1])\t$(recs[i][2])\t$(θ)\t$(nvec)\n")

        end
    end

    # Return
    return nothing

end # end write_model
"""
    `get_cpg_pos(FEAT_ATTS)`

Function that takes the attributes in the GFF file and returns the homozygous and heterozygous CpG
sites in it.

# Examples
```julia-repl
julia> CpelAsm.get_cpg_pos(feat_atts)
```
"""
function get_cpg_pos(atts::Dict{String,Vector{String}})::Array{Vector{Int64},1}

    # Retrieve homozygous CpG sites
    hom = [parse(Int64,replace(s,r"[\[\] ]"=>"")) for s in atts[:"CpGs"]]

    # Retrieve heterozygous CpG sites for each haplotype (if existant)
    het1=Vector{Int64}()
    het2=Vector{Int64}()
    if haskey(atts,"hetCpGg1")
        het1=[parse(Int64,replace(s,r"[\[\] ]"=>"")) for s in atts[:"hetCpGg1"]]
    end
    if haskey(atts,"hetCpGg2")
        het2=[parse(Int64,replace(s,r"[\[\] ]"=>"")) for s in atts[:"hetCpGg2"]]
    end

    return [hom,sort!(append!(het1,hom)),sort!(append!(het2,hom))]

end # end get_cpg_pos
"""
    `get_ns(CPG_POS,G_MAX,HAP_ST,HAP_END)`

Functions that returns [N1,...,NK] given the position of the CpG sites & maximum granularity G_MAX.

# Examples
```julia-repl
julia> CpelAsm.get_ns([100,200,300,350],200,90,400)
2-element Array{Int64,1}:
 2
 2
```
"""
function get_ns(cpg_pos::Vector{Int64},g_max::Int64,hap_st::Int64,hap_end::Int64)::Vector{Int64}

    # Check if need to partition
    (hap_end-hap_st)<=g_max && return [length(cpg_pos)]

    # Get K of model
    k = ceil((hap_end-hap_st+1)/g_max)

    # Get delimiters
    delimiters = vcat(collect(hap_st:(hap_end-hap_st+1)/k:hap_end),hap_end)

    # Partition CpG sites
    ids = [findfirst(y->y>=x,delimiters)-1 for x in cpg_pos]

    # Return [N1,...,NK] vector
    return [sum(ids.==i) for i in unique(ids)]

end # end get_ns
"""
    `mean_cov(XOBS)`

Function returns the average coverage per CpG given some observations `XOBS`.

# Examples
```julia-repl
julia> xobs=[[1,-1] for i=1:10]; append!(xobs,[[1,0] for i=1:10]);
julia> CpelAsm.mean_cov(xobs)
15.0
```
"""
function mean_cov(xobs::Array{Vector{Int64},1})::Float64

    # Return 0 if no observations
    return length(xobs)>0 ? norm(hcat(xobs...),1)/length(xobs[1]) : 0.0

end # end mean_cov
"""
    `mean_cov_sr(XOBS)`

Function returns the average coverage per CpG given some observations `XOBS` per subregion.

# Examples
```julia-repl
julia> n=[1,1]; xobs=[[1,-1] for i=1:10]; append!(xobs,[[1,0] for i=1:10]);
julia> CpelAsm.mean_cov_sr(xobs,n)
2-element Array{Float64,1}:
 20.0
 10.0
```
"""
function mean_cov_sr(xobs::Array{Vector{Int64},1},n::Vector{Int64})::Vector{Float64}

    # Obtain columns
    rows = [(sum(n[1:(i-1)])+1):sum(n[1:i]) for i=2:length(n)]
    pushfirst!(rows,1:n[1])

    # Return 0 if no observations
    return length(xobs)>0 ? [norm(hcat(xobs...)[rows[i],:],1)/n[i] for i=1:length(rows)] :
        zeros(Float64,length(n))

end # end mean_cov_sr
"""
    `get_outpaths(OUTDIR,PREFIX)`

Function that given outdir and prefix returns paths.

# Examples
```julia-repl
julia> CpelAsm.get_outpaths(outdir,prefix)
```
"""
function get_outpaths(outdir::String,prefix::String)

    # Paths
    tobs_path = "$(outdir)/$(prefix)" .* ["_mml1","_mml2","_nme1","_nme2","_pdm"] .* ".bedGraph"
    tnull_path = "$(outdir)/$(prefix)" .* ["_tmml_null","_tnme_null","_tpdm_null"] .* ".bedGraph"
    p_path = "$(outdir)/$(prefix)" .* ["_tmml_pvals","_tnme_pvals","_tpdm_pvals"] .* ".bedGraph"

    # Return path for Tobs and Tnull
    return tobs_path,tnull_path,p_path

end # end get_outpaths
"""
    `sort_bedgraphs(BEDGRAPH_FILES)`

Function that sorts bedGraph files in vector BEDGRAPH_FILES.

# Examples
```julia-repl
julia> CpelAsm.sort_bedgraphs(bg_files)
```
"""
function sort_bedgraphs(bg_files::Vector{String})

    # Loop over files
    for f in bg_files
        Sys.iswindows() || run(`sort -V -k 1,1 -k 2,2n -o $f $f`)
    end

    # Return
    return nothing

end # end sort_bedgraphs
"""
    `print_log(MESSAGE)`

Function that prints MESSAGE to stderr.

# Examples
```julia-repl
julia> CpelAsm.print_log("Hello")
[2020-03-30 16:24:18]: Hello
```
"""
function print_log(mess::String)

    # Right format for date
    date = Dates.format(now(), "yyyy-mm-dd HH:MM:SS")
    println(stderr,"[$(date)]: "  * mess)
    flush(stderr)

    # Return
    return nothing

end # end print_log
"""
    `proc_obs_hap(HAP,CHR,CHR_SIZE,BAM1,BAM2,GFF,FA,OUT_PATHS,PE,GMAX,COV_THS,TRIM)`

Function that computes observed MML, NME, and PDM for given a haplotype.

```julia-repl
julia> proc_obs_hap(HAP,CHR,CHR_SIZE,BAM1,BAM2,GFF,FA,OUT_PATHS,PE,GMAX,COV_THS,TRIM)
```
"""
function proc_obs_hap(hap::GFF3.Record,chr::String,chr_size::Int64,bam1::String,bam2::String,
                      gff::String,fa::String,out_paths::Vector{String},pe::Bool,g_max::Int64,
                      cov_ths::Int64,trim::NTuple{4,Int64},bound_check::Bool)::Vector{Tuple{Int64,
                      Int64,Float64,Int64,Int64}}

    # Empty output
    nan_out = [(0,0,0.0,0,0),(0,0,0.0,0,0),(0,0,0.0,0,0),(0,0,0.0,0,0),(0,0,0.0,0,0)]

    # Get window
    hap_st = GFF3.seqstart(hap)
    hap_end = GFF3.seqend(hap)

    # Get homozygous & heterozygous CpG sites (1:hom, 2:hap1, 3: hap2)
    cpg_pos = get_cpg_pos(Dict(GFF3.attributes(hap)))
    n = get_ns(cpg_pos[1],g_max,hap_st,hap_end)
    n1 = get_ns(cpg_pos[2],g_max,hap_st,hap_end)
    n2 = get_ns(cpg_pos[3],g_max,hap_st,hap_end)
    length(n)>0 || return nan_out

    # Get binary vector with homozygous CpG sites
    z1 = BitArray([p in cpg_pos[1] ? true : false for p in cpg_pos[2]])
    z2 = BitArray([p in cpg_pos[1] ? true : false for p in cpg_pos[3]])

    # Get vectors from BAM1 overlapping haplotype, assign md to het CpGs, & compute average coverage
    xobs1 = read_bam(bam1,chr,hap_st,hap_end,cpg_pos[2],chr_size,pe,trim)
    xobs1 = map(x->x.*z1,xobs1)
    xobs1_hom = map(x->x[z1],xobs1)
    obs_per_cpg1 = get_obs_per_cpg(xobs1_hom)
    (cov_ths<=mean_cov(xobs1)<=400) && (sum(obs_per_cpg1.==0)<=1.0/5.0*sum(n)) || return nan_out

    # Get vectors from BAM2 overlapping haplotype, assign md to het CpGs, & compute average coverage
    xobs2 = read_bam(bam2,chr,hap_st,hap_end,cpg_pos[3],chr_size,pe,trim)
    xobs2 = map(x->x.*z2,xobs2)
    xobs2_hom = map(x->x[z2],xobs2)
    obs_per_cpg2 = get_obs_per_cpg(xobs2_hom)
    (cov_ths<=mean_cov(xobs2)<=400) && (sum(obs_per_cpg2.==0)<=1.0/5.0*sum(n)) || return nan_out

    # Estimate each single-allele model and check if on boundary of parameter space
    θ1 = est_theta_sa(n1,xobs1)
    bound_check && check_boundary(θ1) && return nan_out
    θ2 = est_theta_sa(n2,xobs2)
    bound_check && check_boundary(θ2) && return nan_out

    # Estimate intermediate quantities
    if all(z1)
        # In case all are homozygous CpG sites
        ∇1 = get_grad_logZ(n1,θ1)
    else
        # In case NOT all are homozygous CpG sites
        ex1 = comp_ex(n1,θ1[1:(end-1)],θ1[end])
        exx1 = comp_exx(n1,θ1[1:(end-1)],θ1[end])
    end
    if all(z2)
        # In case all are homozygous CpG sites
        ∇2 = get_grad_logZ(n2,θ2)
    else
        # In case NOT all are homozygous CpG sites
        ex2 = comp_ex(n2,θ2[1:(end-1)],θ2[end])
        exx2 = comp_exx(n2,θ2[1:(end-1)],θ2[end])
    end

    # Compute output
    mml1 = all(z1) ? comp_mml_∇(n1,∇1) : comp_mml(z1,ex1)
    mml2 = all(z2) ? comp_mml_∇(n2,∇2) : comp_mml(z2,ex2)
    nme1 = all(z1) ? comp_nme_∇(n1,θ1,∇1) : comp_nme(z1,n1,θ1[1:(end-1)],θ1[end],ex1,exx1)
    nme2 = all(z2) ? comp_nme_∇(n2,θ2,∇2) : comp_nme(z2,n2,θ2[1:(end-1)],θ2[end],ex2,exx2)
    pdm = comp_uc(z1,z2,n1,n2,θ1,θ2,nme1,nme2)

    # Return output
    return [(hap_st,hap_end,mml1,sum(n),length(n)),
            (hap_st,hap_end,mml2,sum(n),length(n)),
            (hap_st,hap_end,nme1,sum(n),length(n)),
            (hap_st,hap_end,nme2,sum(n),length(n)),
            (hap_st,hap_end,pdm,sum(n),length(n))]

end
"""
    `comp_tobs(BAM1_PATH,BAM2_PATH,GFF_PATH,FA_PATH,OUT_PATHS)`

Function that computes MML1/2, NME1/2, and PDM on a pair of BAM files (allele 1, allele 2) given a
(phased) VCF file that contains the heterozygous SNPs and a FASTA file that contains the reference
genome.

# Examples
```julia-repl
julia> comp_tobs(BAM1_PATH,BAM2_PATH,GFF_PATH,FA_PATH,OUT_PATHS)
```
"""
function comp_tobs(bam1::String,bam2::String,gff::String,fa::String,out_paths::Vector{String};
                   pe::Bool=true,g_max::Int64=300,cov_ths::Int64=5,trim::NTuple{4,Int64}=
                   (0,0,0,0),bound_check::Bool=true)::Nothing

    # BigWig output files
    mml1_path,mml2_path,nme1_path,nme2_path,pdm_path = out_paths

    # Find chromosomes
    reader_fa = open(FASTA.Reader,fa,index=fa*".fai")
    chr_sizes = reader_fa.index.lengths
    chr_names = reader_fa.index.names
    close(reader_fa)

    # Loop over chromosomes
    for chr in sort(chr_names)

        # Get windows pertaining to current chromosome
        print_log("Processing chromosome $(chr) ...")
        haps_chr = read_gff_chr(gff,chr)
        chr_size = chr_sizes[findfirst(x->x==chr, chr_names)]

        # Process them
        out_pmap = pmap(hap -> proc_obs_hap(hap,chr,chr_size,bam1,bam2,gff,fa,out_paths,pe,g_max,cov_ths,trim,bound_check),haps_chr)

        length(out_pmap)>0 || continue

        # Add last to respective bedGraph file
        write_tobs([x[1] for x in out_pmap],chr,mml1_path)
        write_tobs([x[2] for x in out_pmap],chr,mml2_path)
        write_tobs([x[3] for x in out_pmap],chr,nme1_path)
        write_tobs([x[4] for x in out_pmap],chr,nme2_path)
        write_tobs([x[5] for x in out_pmap],chr,pdm_path)

    end

    # Sort
    sort_bedgraphs([mml1_path,mml2_path,nme1_path,nme2_path,pdm_path])

    # Return
    return nothing

end # end comp_tobs
"""
    `write_tnull(RECORDS,PATH)`

Function that writes null records in `RECORDS` into `PATH`.

# Examples
```julia-repl
julia> CpelAsm.write_tnull(RECORDS,PATH)
```
"""
function write_tnull(recs::Vector{Tuple{Int64,Int64,Float64}},path::String)

    # Open output bedGraph file in append mode (no need to close it)
    open(path, "a") do f
        for i=1:length(recs)
            recs[i][1]!=0 || continue
            write(f,join(string.(collect(recs[i])),"\t"),"\n")
        end
    end

    # Return
    return nothing

end # end write_tnull
"""
    `get_kstar_table_gff(GFF_PATH,CHR_NAMES,BLK_SIZE)`

Function that returns a table with maximum K for each N in GFF_PATH.

# Examples
```julia-repl
julia> CpelAsm.get_kstar_table_gff(GFF_PATH,CHR_NAMES,BLK_SIZE)
```
"""
function get_kstar_table_gff(gff::String,chr_names::Vector{String},g_max::Int64)::Dict{Int64,Int64}

    # Loop over chromosomes
    table = Dict{Int64,Int64}()
    for chr in chr_names
        hom_win = read_gff_chr(gff,chr)
        # Loop over windows in chromosome
        for win in hom_win
            # Get window info
            win_st = GFF3.seqstart(win)
            win_end = GFF3.seqend(win)
            cpg_pos = get_cpg_pos(Dict(GFF3.attributes(win)))
            n = get_ns(cpg_pos[1],g_max,win_st,win_end)
            # Update table if necessary
            haskey(table,sum(n)) || (table[sum(n)]=length(n))
            length(n)>table[sum(n)] && (table[sum(n)]=length(n))
        end
    end

    # Return
    return table

end # end get_kstar_table_gff
"""
    `get_kstar_table_tobs(TOBS_PATH)`

Function that returns a table with maximum K for each N in TOBS_PATH.

# Examples
```julia-repl
julia> CpelAsm.get_kstar_table_tobs(TOBS_PATH)
```
"""
function get_kstar_table_tobs(tobs_path::String)::Dict{Int64,Int64}

    # Get (N,K) pairs from Tobs
    nk_pairs = readdlm(tobs_path,'\t',Any)[:,5:6]

    # Loop over pairs
    out = Dict{Int64,Int64}()
    @inbounds for i=1:size(nk_pairs)[1]
        if haskey(out,nk_pairs[i,1])
            # If key exists update if necessary
            out[nk_pairs[i,1]] = nk_pairs[i,2]>out[nk_pairs[i,1]] ? nk_pairs[i,2] : out[nk_pairs[i,1]]
        else
            # If key doesn't exist create it
            out[nk_pairs[i,1]] = nk_pairs[i,2]
        end
    end

    # Return
    return out

end # end get_kstar_table_tobs
"""
    `get_obs_per_cpg(XOBS)`

Function that returns the number of observations per CpG site as a vector.

# Examples
```julia-repl
julia> CpelAsm.get_obs_per_cpg(xobs)
```
"""
function get_obs_per_cpg(xobs::Array{Vector{Int64},1})::Vector{Int64}

    # Get observations per CpG site
    out = length(xobs)>0 ? vcat(sum(abs.(hcat(xobs...)),dims=2)...) : [0]

    # Return
    return out

end # end get_obs_per_cpg
"""
    `screen_haps(GFF,BAM,FA,PE,COV_THS,NTOT,TRIM,CHR,CHR_SIZE)`

Function that returns a set of haplotypes with N=Ntot that satisfies the coverage requirements.

# Examples
```julia-repl
julia> CpelAsm.screen_haps(gff,bam,fa,pe,cov_ths,ntot,trim,chr,chr_size)
```
"""
function screen_haps(gff::String,bam::String,fa::String,pe::Bool,cov_ths::Int64,n::Vector{Int64},
                     trim::NTuple{4,Int64},chr::String,chr_size::Int64;fold::Int64=5)::Vector{GFF3.Record}

    # Ntot
    ntot = sum(n)

    # Load genomic features from a GFF3 file
    out_haps = Vector{GFF3.Record}()
    hom_haps = open(collect,GFF3.Reader,gff)

    # Filter based on N and chromosomes
    filter!(x -> (GFF3.score(x)>=ntot) && (GFF3.seqid(x) == chr),hom_haps)
    length(hom_haps)>0 || return out_haps

    # Loop over homozygous blocks
    for hap in hom_haps

        # All CpG sites in haplotype
        cpg_pos = get_cpg_pos(Dict(GFF3.attributes(hap)))[1]

        # Get ntot contiguous CpG sites at random
        ind = sample(1:(length(cpg_pos)-ntot+1))
        cpgs = cpg_pos[ind:(ind+ntot-1)]

        # Obtain observations overlapping with CpG sites
        xobs = read_bam(bam,chr,cpgs[1],cpgs[end],cpgs,chr_size,pe,trim)
        length(xobs)>0 || continue

        # Obtain observations per CpG site
        obs_per_cpg = get_obs_per_cpg(xobs)

        # NOTE: we leave xFOLD as much per allele to ensure variety. This determines the rate at
        # which null statistics are produced. However, we might have trouble finding regions
        # that satisfy both conditions. Consider lowering FOLD in that case.

        # Accept candidate haplotype if coverage is OK
        if (2*fold*cov_ths<=mean_cov(xobs)<=400) && (sum(obs_per_cpg.==0)<=1.0/5.0*ntot)
            new_hap = "$(chr)\t.\t.\t$(cpgs[1])\t$(cpgs[end])\t$(ntot)\t.\t.\tCpGs=$(cpgs)"
            push!(out_haps,GFF3.Record(new_hap))
        end

        # Break if enough haps for given chromosome
        length(out_haps)>500 && break

    end

    # Return
    return out_haps

end # end screen_haps
"""
    `get_nvec_kstar(NTOT,KSTAR)`

Function that returns vector n with Ntot CpG sites divided into Kstar subregions.

# Examples
```julia-repl
julia> Random.seed!(1234)
julia> CpelAsm.get_nvec_kstar(10,4)
4-element Array{Int64,1}:
 3
 2
 2
 3
```
"""
function get_nvec_kstar(ntot::Int64,kstar::Int64)::Vector{Int64}

    # Partition CpG sites into Kstar subregions
    n = div(ntot,kstar) * ones(Int64,kstar)
    if sum(n)<ntot
        cpg_ids = rand(Categorical(1.0/kstar*ones(kstar)),ntot-sum(n))
        n += [length(findall(x->x==i,cpg_ids)) for i=1:kstar]
    end

    # Return n
    return n

end # end get_nvec_kstar
"""
    `cov_obs_part(XOBS,N,COV_THS,COV_A,COV_B)`

Function that partitions observations into two sets of observations that satisfy a minimum
coverage within [COV_THS-COV_A,COV_THS+COV_B] for each subregion. If not partition is found that
satisfies it, then two empty vectors are returned.

# Examples
```julia-repl
julia> CpelAsm.cov_obs_part(xobs,n,cov_ths,cov_a,cov_b)
```
"""
function cov_obs_part(xobs::Vector{Vector{Int64}},n::Vector{Int64},cov_ths::Int64,cov_a::Float64,
                      cov_b::Float64)::Tuple{Vector{Vector{Int64}},Vector{Vector{Int64}}}

    ## NOTE: Part that takes care of minimum coverage. After this step an allele can have
    ## more than cov_ths+b. Here we just guarantee mean_cov(xobsi)≧cov_ths.

    # Divide into two sets of equal size while checking bot have at least cov_ths
    ct = 1
    fail = false
    ntot = sum(n)
    xobs1 = Vector{Int64}()
    xobs2 = Vector{Int64}()
    while true
        # Partition
        ind = sample(1:length(xobs),div(length(xobs),2),replace=false)
        xobs1 = xobs[ind]
        xobs2 = xobs[setdiff(1:length(xobs),ind)]
        obs_per_cpg_1 = get_obs_per_cpg(xobs1)
        obs_per_cpg_2 = get_obs_per_cpg(xobs2)
        # If coverage condition are OK break loop
        (sum(obs_per_cpg_1.==0)<=1.0/5.0*ntot) && (sum(obs_per_cpg_2.==0)<=1.0/5.0*ntot) &&
            (mean_cov(xobs1)>=cov_ths) && (mean_cov(xobs2)>=cov_ths) && break
        # Check we haven't exceeded the maximum number of permutations
        if ct>20
            # We failed and break
            fail=true
            break
        else
            # Increase counter
            ct += 1
        end

    end

    # Return empty arrays if fail=true
    fail==false || return [],[]

    ## NOTE: Part extra to limit number of observations used. The goal is to reduce the size of
    ## xobs1 & xobs2 until both coverages are within [cov_ths-cov_a,cov_ths+cov_b].

    # Check for [cov_ths-cov_a,cov_ths+cov_b] for xobs1
    less_ok = false
    while true
        # Check if we can keep decreasing number
        obs_per_cpg_1 = get_obs_per_cpg(xobs1[2:end])
        less_ok = mean_cov(xobs1[2:end])>=(cov_ths-cov_a) && mean_cov(xobs1)>=(cov_ths+cov_b) &&
            (sum(obs_per_cpg_1.==0)<=1.0/5.0*ntot)
        # If less is okay, decrease in one read, otherwise that's it
        xobs1 = less_ok ? xobs1=xobs1[2:end] : xobs1
        less_ok || break
    end

    # Check if resulting observations meet b-bound as well
    fail = mean_cov(xobs1)<=(cov_ths+cov_b) ? false : true
    fail==false || return [],[]

    # Check for [cov_ths-cov_a,cov_ths+cov_b] for xobs2
    less_ok = false
    while true
        # Check if we can keep decreasing number
        obs_per_cpg_2 = get_obs_per_cpg(xobs2[2:end])
        less_ok = mean_cov(xobs2[2:end])>=(cov_ths-cov_a) && mean_cov(xobs2)>=(cov_ths+cov_b) &&
            (sum(obs_per_cpg_2.==0)<=1.0/5.0*ntot)
        # If less is okay, decrease in one read, otherwise that's it
        xobs2 = less_ok ? xobs2=xobs2[2:end] : xobs2
        less_ok || break
    end

    # Check if resulting observations meet b-bound as well
    fail = mean_cov(xobs2)<=(cov_ths+cov_b) ? false : true
    fail==false || return [],[]
        # println("Third OK: $(mean_cov(xobs2))")

    # Return partition
    return xobs1,xobs2

end # end cov_obs_part
"""
    `swap_pair(XOBS1,XOBS2)`

Function that swaps a pair of reads between XOBS1 and XOBS2.

# Examples
```julia-repl
julia> CpelAsm.swap_pair(xobs1,xobs2)
```
"""
function swap_pair(xobs1::Vector{Vector{Int64}},xobs2::Vector{Vector{Int64}})::Tuple{Vector{Vector{Int64}},Vector{Vector{Int64}}}

    # Sample indices to swap
    ind1 = sample(1:length(xobs1))
    ind2 = sample(1:length(xobs2))

    # Swap elements
    push!(xobs1,xobs2[ind2])
    push!(xobs2,xobs1[ind1])

    # Delete them in their origin
    deleteat!(xobs1,ind1)
    deleteat!(xobs2,ind2)

    # Return partition
    return xobs1,xobs2

end # end swap_pair
"""
    `proc_null_hap(HAP,NTOT,BAM,HET_GFF,HOM_GFF,FA,KSTAR,OUT_PATHS,PE,GMAX,COV_THS,COV_A,COV_B,TRIM,NNULL,NMAX,NSUBSET,CHR_DICT)`

Function that computes a null Tmml, Tnme, and Tpdm given a homozygous haplotype.

# Examples
```julia-repl
julia> proc_null_hap(HAP,NTOT,BAM,HET_GFF,HOM_GFF,FA,KSTAR,OUT_PATHS,PE,GMAX,COV_THS,COV_A,COV_B,TRIM,NNULL,NMAX,NSUBSET,CHR_DICT)
```
"""
function proc_null_hap(hap::GFF3.Record,ntot::Int64,bam::String,het_gff::String,hom_gff::String,
                       fa::String,kstar::Int64,out_paths::Vector{String},pe::Bool,g_max::Int64,
                       cov_ths::Int64,cov_a::Float64,cov_b::Float64,trim::NTuple{4,Int64},
                       n_null::Int64,n_max::Int64,n_subset::Vector{Int64},chr_dic::Dict{String,
                       Int64},bound_check::Bool)::Vector{Tuple{Float64,Float64,Float64}}

    # Get homozygous window
    chr = GFF3.seqid(hap)
    chr_size = chr_dic[chr]

    # Sample ntot contiguous CpG sites and partition into Kstar
    cpg_pos = get_cpg_pos(Dict(GFF3.attributes(hap)))[1]
    n = get_nvec_kstar(ntot,kstar)

    # Check average coverage is within normal limits (watch for repetitive regions)
    xobs = read_bam(bam,chr,cpg_pos[1],cpg_pos[end],cpg_pos,chr_size,pe,trim)

    # Obtain a number of null statistics with different permutations
    success = false
    xobs1 = Vector{Vector{Int64}}()
    xobs2 = Vector{Vector{Int64}}()
    stats = Vector{Tuple{Float64,Float64,Float64}}()
    for i=1:10

        # If previous partition successful, just swap two reads. Otherwise, re-partition
        xobs1,xobs2 = success ? swap_pair(xobs1,xobs2) : cov_obs_part(xobs,n,cov_ths,cov_a,cov_b)
        (length(xobs1)>0) && (length(xobs2)>0) || continue

        # Reset success
        success = false

        # Estimate each single-allele model and check if on boundary of parameter space
        θ1 = est_theta_sa(n,xobs1)
        bound_check && check_boundary(θ1) && continue
        θ2 = est_theta_sa(n,xobs2)
        bound_check && check_boundary(θ2) && continue

        # Estimate moments
        ∇1 = get_grad_logZ(n,θ1)
        ∇2 = get_grad_logZ(n,θ2)

        # Estimate output quantities
        nme1 = round(comp_nme_∇(n,θ1,∇1);digits=8)
        nme2 = round(comp_nme_∇(n,θ2,∇2);digits=8)
        pdm = round(comp_uc(trues(ntot),trues(ntot),n,n,θ1,θ2,nme1,nme2);digits=8)

        # Append stats
        push!(stats,(round(abs(comp_mml_∇(n,∇1)-comp_mml_∇(n,∇2));digits=8),abs(nme1-nme2),pdm))

        # Update success for next iteration
        success = true

    end

    # Return output
    return stats

end
"""
    `comp_tnull(BAM_PATH,GFF_PATH,FA_PATH,OUT_PATH)`

Function that computes null Tmml, Tnme, and Tpdm from a BAM files at locations given by GFF that
contains the windows with not genetic variants and a FASTA file that contains the reference genome.

# Examples
```julia-repl
julia> comp_tnull(BAM_PATH,GFF_PATH,FA_PATH,OUT_PATH)
```
"""
function comp_tnull(bam::String,het_gff::String,hom_gff::String,fa::String,tobs_path::Vector{String},
                    out_paths::Vector{String};pe::Bool=true,g_max::Int64=500,cov_ths::Int64=5,
                    cov_a::Float64=0.0,cov_b::Float64=2.0,trim::NTuple{4,Int64}=(5,0,5,0),
                    n_null::Int64=1000,n_max::Int64=20,n_subset::Vector{Int64}=collect(1:n_max),
                    bound_check::Bool=true)::Nothing

    # BigWig output files
    tmml_path,tnme_path,pdm_path = out_paths

    # Find relevant chromosomes and sizes
    reader_fa = open(FASTA.Reader,fa,index=fa*".fai")
    chr_dic = Dict(zip(reader_fa.index.names,reader_fa.index.lengths))
    chr_names = collect(keys(chr_dic))
    close(reader_fa)

    # Obtain Kstar for every N & maximum N to analyze
    print_log("Generating Kstar table ...")
    kstar_table = get_kstar_table_tobs(tobs_path[1])

    # Cap Nmax
    n_max = min(maximum(keys(kstar_table)),n_max)

    # Loop over Ns of interest
    for ntot in sort(collect(intersect(keys(kstar_table),n_subset)))

        # Skip if greater than n_max
        ntot>n_max && continue

        # Get kstar
        kstar = kstar_table[ntot]

        # Get n vector
        n = get_nvec_kstar(ntot,kstar)

        # Print current N
        out_tmml = Vector{Tuple{Int64,Int64,Float64}}()
        out_tnme = Vector{Tuple{Int64,Int64,Float64}}()
        out_pdm = Vector{Tuple{Int64,Int64,Float64}}()

        # Screen all haplotypes for those with enough coverage
        print_log("Screening haplotypes ...")
        haps = vcat(pmap(chr->screen_haps(hom_gff,bam,fa,pe,cov_ths,n,trim,chr,chr_dic[chr]),chr_names)...)

        # Keep only n_null for memory efficiency
        haps = length(haps)>n_null ? haps[sample(1:length(haps),n_null)] : haps

        # Sample null stats until desired number is reached or maximum number of attempts
        print_log("Sampling a total of $(n_null) null statistics for N=$(ntot) & K*=$(kstar) ...")
        print_log("Candidate haplotypes with N=$(ntot) are $(length(haps)) ...")
        i = 1
        while length(out_tmml)<n_null

            # Process them in parallel
            out_pmap = vcat(pmap(hap->proc_null_hap(hap,ntot,bam,het_gff,hom_gff,fa,kstar,out_paths,pe,g_max,
                            cov_ths,cov_a,cov_b,trim,n_null,n_max,n_subset,chr_dic,bound_check),haps)...)

            # Save succesful runs
            if length(out_pmap)>0
                print_log("Succesfully generated $(length(out_pmap)) null stats ...")
                append!(out_tmml,map(stat->(ntot,kstar,stat[1]),out_pmap))
                append!(out_tnme,map(stat->(ntot,kstar,stat[2]),out_pmap))
                append!(out_pdm,map(stat->(ntot,kstar,stat[3]),out_pmap))
            else
                print_log("No null statistics computed in $(i)-th attempt ...")
            end

            # Break if run for too long
            if i>1000
                print_log("Exceeded $(i-1) attempts in comp_tnull ...")
                break
            end

            # Increase counter
            i += 1

        end

        # Add last to respective bedGraph file
        write_tnull(out_tmml,tmml_path)
        write_tnull(out_tnme,tnme_path)
        write_tnull(out_pdm,pdm_path)

    end # end N for loop

    # Sort files
    sort_bedgraphs([tmml_path,tnme_path,pdm_path])

    # Print message
    print_log("Done with null statistics ...")

    # Return nothing
    return nothing

end # end comp_tnull
"""
    `comp_pvals_stat(TOBS_PATH,TNULL_PATH,N_MAX)`

Function that computes adjusted p-values for each haplotype in TOBS_PATH from the empirical
distribution of the pertaining statistic in TNULL_PATH. The output file has the genomic
coordinates, the value of the statistic, and the adjusted p-value.

# Examples
```julia-repl
julia> CpelAsm.comp_pvals_stat(tobs_path,tnull_path)
```
"""
function comp_pvals_stat(tobs_path::Vector{String},tnull_path::String,p_path::String,n_max::Int64,
                         n_null::Int64)

    # Get null stats
    filesize(tnull_path)>0 || return nothing
    tnull = readdlm(tnull_path,'\t',Any)

    # Consider case Tmml/Tnme VS Tpdm
    if length(tobs_path)>1
        # Tmml/Tnme. NOTE: both files need to have the same order of haplotypes!
        tobs = readdlm(tobs_path[1],'\t',Any)
        tobs[:,4] = abs.(tobs[:,4]-readdlm(tobs_path[2],'\t',Any)[:,4])
    else
        # Tpdm
        tobs = readdlm(tobs_path[1],'\t',Any)
    end

    # Initialize p-value matrix
    pvals = tobs[:,1:5]
    pvals[:,5] .= "NO NULL DATA"

    # Compute p-value for each N
    for n in unique(tobs[:,5])
        # Check we can compute p-value
        (n<=n_max) && (sum(tnull[:,1].==n)>=n_null) || continue
        # Get indexes of Tobs with N=n
        ind = tobs[:,5].==n
        # Compute p-values
        pvals[ind,5] = pmap(t -> (1.0+sum(tnull[tnull[:,1].==n,3].>=t))/(1.0+sum(tnull[:,1].==n)),tobs[ind,4])
    end

    # Multiple hypothesis testing correction. NOTE: should we apply BH on each N independently?
    ind = pvals[:,5].!="NO NULL DATA"
    pvals[ind,5] = MultipleTesting.adjust(convert(Vector{Float64},pvals[ind,5]),BenjaminiHochberg())

    # Write to output
    open(p_path,"w") do io
        writedlm(io,pvals,'\t')
    end

    # Return
    return nothing

end # end comp_pvals_stat
"""
    `comp_pvals(TOBS_PATH,TNULL_PATH,PVALS_PATH,N_MAX)`

Function that computes p-values from empirical null distribution.

# Examples
```julia-repl
julia> comp_pvals(tobs_path,tnull_path,p_path,n_max)
```
"""
function comp_pvals(tobs_path::Vector{String},tnull_path::Vector{String},p_path::Vector{String},
                    n_max::Int64,n_null::Int64)

    # Compute three sets of pvalues
    print_log("Computing p-values Tmml ...")
    comp_pvals_stat(tobs_path[1:2],tnull_path[1],p_path[1],n_max,n_null)
    print_log("Computing p-values Tnme ...")
    comp_pvals_stat(tobs_path[3:4],tnull_path[2],p_path[2],n_max,n_null)
    print_log("Computing p-values Tpdm ...")
    comp_pvals_stat([tobs_path[5]],tnull_path[3],p_path[3],n_max,n_null)
    print_log("Done with p-values ...")

    # Return
    return nothing

end # end comp_pvals
"""
    `run_analysis(BAM1_PATH,BAM2_PATH,BAMU_PATH,VCF_PATH,FA_PATH,OUT_PATH)`

Function that estimates MML1/2, NME1/2, and UC on a pair of BAM files (allele 1, allele 2) given a
(phased) VCF file that contains the heterozygous SNPs and a FASTA file that contains the reference
genome. An empirical null distribution for each quantity is estimated from a set of representative
homozygous regions from unassigned BAM records in BAMU_PATH. A p-value is computed for each
haplotype in the VCF file.

# Examples
```julia-repl
julia> run_analysis(BAM1_PATH,BAM2_PATH,BAMU_PATH,VCF_PATH,FA_PATH,OUT_PATH)
```
"""
function run_analysis(bam1::String,bam2::String,bamu::String,vcf::String,fa::String,outdir::String;
                      pe::Bool=true,g_max::Int64=500,win_exp::Int64=100,cov_ths::Int64=5,
                      cov_a::Float64=0.0,cov_b::Float64=2.0,trim::NTuple{4,Int64}=(5,0,5,0),
                      n_null::Int64=1000,n_max::Int64=25,n_subset::Vector{Int64}=collect(1:n_max),
                      bound_check::Bool=true)::Nothing

    # Print initialization of juliASM
    print_log("Starting CpelAsm analysis ...")

    # Check index files exist
    if !(isfile.(bam1*".bai",bam2*".bai",bamu*".bai",fa*".fai"))
        print_log("Index files for BAM or FASTA missing. Exiting julia ...")
        exit(1)
    end

    # Create output folder if it doesn't exist
    isdir(outdir) || mkdir(outdir)

    # BedGraph output files
    prefix_sample = String(split(basename(bam1),".")[1])
    tobs_path,tnull_path,p_path = get_outpaths(outdir,prefix_sample)

    # Check for existance of at least an output files
    if all(isfile.(vcat(tobs_path,tnull_path)))
        print_log("At least an output file already exists. Exiting julia ...")
        exit(1)
    end

    # Create gff file with heterozygous loci
    print_log("Reading in VCF & FASTA files ...")
    het_gff = outdir * split(basename(vcf),".")[1] * "_het.cpelasm.gff"
    hom_gff = outdir * split(basename(vcf),".")[1] * "_hom.cpelasm.gff"
    if !(isfile(het_gff))
        print_log("Generating CpelAsm GFF files ...")
        gen_gffs([het_gff,hom_gff],fa,vcf,win_exp,n_max)
    end

    # Compute observed statistics from heterozygous loci
    print_log("Computing observed statistics in heterozygous loci...")
    comp_tobs(bam1,bam2,het_gff,fa,tobs_path;pe=pe,g_max=g_max,cov_ths=cov_ths,trim=trim,bound_check=bound_check)

    # Compute null statistics from homozygous loci
    print_log("Generating null statistics in homozygous loci ...")
    comp_tnull(bamu,het_gff,hom_gff,fa,tobs_path,tnull_path;pe=pe,g_max=g_max,cov_ths=cov_ths,cov_a=cov_a,
               cov_b=cov_b,trim=trim,n_null=n_null,n_max=n_max,n_subset=n_subset,bound_check=bound_check)

    # Compute null statistics from heterozygous loci
    print_log("Computing p-values in heterozygous loci ...")
    comp_pvals(tobs_path,tnull_path,p_path,n_max,n_null)

    # Print number of features interrogated and finished message
    print_log("Done.")

    # Return
    return nothing

end # end run_analysis


text:
# Advent of Code 2020 - day 18

list = readlines("day18_input.txt")

function reduce_left_first(line)
    while occursin("(", line)
        part = match(r".*(\([^\)]*\)).*", line)[1]
        line = replace(line, part => reduce_left_first(part[2:end-1]))
    end
    while occursin("+", line) || occursin("*", line)
        parts = match(r"(\d+\s[\+|\*]\s\d+)(.*)", line)
        line = string(eval(Meta.parse(parts[1]))) * parts[2]
    end
    return parse(Int, line)
end

function reduce_plus_first(line)
    while occursin("(", line)
        part = match(r".*(\([^\)]*\)).*", line)[1]
        line = replace(line, part => reduce_plus_first(part[2:end-1]))
    end
    while occursin("+", line)
        parts = match(r"([^\+]*\s|)(\d+\s\+\s\d+)(.*)", line)
        line = parts[1] * string(eval(Meta.parse(parts[2]))) * parts[3]
    end
    return eval(Meta.parse(line))
end

function evaluate_list_left_first(list)
    total = 0
    for i in 1:length(list)
        total += reduce_left_first(list[i])
    end
    return total
end

function evaluate_list_plus_first(list)
    total = 0
    for i in 1:length(list)
        total += reduce_plus_first(list[i])
    end
    return total
end

@show reduce_left_first("1 + 2 * 3 + 4 * 5 + 6") == 71
@show reduce_left_first("1 + (2 * 3) + (4 * (5 + 6))") == 51
@show reduce_left_first("2 * 3 + (4 * 5)") == 26
@show reduce_left_first("5 + (8 * 3 + 9 + 3 * 4 * 3)") == 437
@show reduce_left_first("5 * 9 * (7 * 3 * 3 + 9 * 3 + (8 + 6 * 4))") == 12240
@show reduce_left_first("((2 + 4 * 9) * (6 + 9 * 8 + 6) + 6) + 2 + 4 * 2") == 13632

@show evaluate_list_left_first(list) == 1402255785165

@show reduce_plus_first("1 + 2 * 3 + 4 * 5 + 6") == 231
@show reduce_plus_first("1 + (2 * 3) + (4 * (5 + 6))") == 51
@show reduce_plus_first("2 * 3 + (4 * 5)") == 46
@show reduce_plus_first("5 + (8 * 3 + 9 + 3 * 4 * 3)") == 1445
@show reduce_plus_first("5 * 9 * (7 * 3 * 3 + 9 * 3 + (8 + 6 * 4))") == 669060
@show reduce_plus_first("((2 + 4 * 9) * (6 + 9 * 8 + 6) + 6) + 2 + 4 * 2") == 23340

@show evaluate_list_plus_first(list) == 119224703255966

nothing

ext: jl


text:
using BenchmarkTools, Suppressor

import BenchmarkTools: prunekwargs, hasevals

# Adapted from BenchmarkTools.btime() https://github.com/JuliaCI/BenchmarkTools.jl/blob/5f8b7f2880b68e019d8cd65317f275b1b1d1ca03/src/execution.jl#L474
# Returns the *median* (instead of the *minimum* in btime) elapsed time (in
# seconds instead in pretty time in btime).
# No output about memory and no tuning phase.
macro mybtimes(args...)
    _, params = prunekwargs(args...)
    bench, trial, result = gensym(), gensym(), gensym()
    trialmin = gensym()
    trialmed = gensym()
    trialmean = gensym()
    trialmax = gensym()
    return esc(quote
        local $bench = $BenchmarkTools.@benchmarkable $(args...)
        $BenchmarkTools.warmup($bench)
        local $trial, $result = $BenchmarkTools.run_result($bench)
        local $trialmin = $BenchmarkTools.minimum($trial)
        local $trialmed = $BenchmarkTools.median($trial)
        local $trialmean = $BenchmarkTools.mean($trial)
        local $trialmax = $BenchmarkTools.maximum($trial)
        println(" ", $BenchmarkTools.time($trialmin)/1e9,
        		" ", $BenchmarkTools.time($trialmed)/1e9,
        		" ", $BenchmarkTools.time($trialmean)/1e9,
        		" ", $BenchmarkTools.time($trialmax)/1e9)
        $result
    end)
end

# Returns benchmark results as array of floats instead of a string
macro mybtimes2(args...)
    _, params = prunekwargs(args...)
    bench, trial, result = gensym(), gensym(), gensym()
    trialmin = gensym()
    trialmed = gensym()
    trialmean = gensym()
    trialmax = gensym()
    btimes = gensym()
    return esc(quote
        local $bench = $BenchmarkTools.@benchmarkable $(args...)
        $BenchmarkTools.warmup($bench)
        local $trial, $result = $BenchmarkTools.run_result($bench)
        local $trialmin = $BenchmarkTools.minimum($trial)
        local $trialmed = $BenchmarkTools.median($trial)
        local $trialmean = $BenchmarkTools.mean($trial)
        local $trialmax = $BenchmarkTools.maximum($trial)
        local $btimes = [
        	$BenchmarkTools.time($trialmin)/1e9,
        	$BenchmarkTools.time($trialmed)/1e9,
        	$BenchmarkTools.time($trialmean)/1e9,
        	$BenchmarkTools.time($trialmax)/1e9
        ]
        $result, $btimes
    end)
end

# Create an array of vector lengths of size nrun
# starting with N0 and increasing in geometrica progression.
# ppomag denotes the number of
# elements per order of magnitude.
# Try to keep full powers of 10 at once.
# From: https://github.com/j-fu/julia-tests/blob/3a3d3646e41e13b1ec23f2dea7df7eee5967c551/parallel/jlvtriad.jl#L60
function get_sizes_gp(N0::Int, ppomag::Int, nrun::Int)
    vsz = [N0]
    N = N0
    N0 *= 10
    for irun=1:nrun-1
        N = N * 10^(1.0 / ppomag)
        if (irun%ppomag == 0)
            N = N0
            N0 *= 10
        end
        push!(vsz, Int(ceil(N)))
    end
    return vsz
end

# Create an array of vector lengths of size nrun
# starting with N0 and increasing in arithmetic progression.
# ppomag denotes the number of
# elements per order of magnitude.
# Try to keep full powers of 10 at once.
function get_sizes_ap(N0::Int, ppomag::Int, nrun::Int)
	NN = Int[]
	for i in 1:ceil(nrun/ppomag)+1
		push!(NN, Int(N0 * 10^i))
	end
    sizes = []
    s = N0
    push!(sizes, s)
    irun = 0
	for N in NN
		commdiff = N / ppomag
		for j in 1:ppomag
			if irun == nrun
				break
			end
		 	s = j*commdiff
		 	if sizes[length(sizes)] != s
				push!(sizes, Int(round(s)))
				irun += 1
			end
		end
	end
    return sizes
end

# Return sizes in geometric or arithmetic progression
function get_sizes(prog::String, N0::Int, ppomag::Int, nrun::Int)
	if prog == "gp"
		sizes = get_sizes_gp(N0::Int, ppomag::Int, nrun::Int)
	elseif prog == "ap"
		sizes = get_sizes_ap(N0::Int, ppomag::Int, nrun::Int)
	end
	return sizes
end

# Initialize the matrix, transpose the matrix and compute rhs
function get_linsys(matrix::String, n::Int)
	if n <= 0
		# SuiteSparse/MatrixMarket without rhs
		if n == 0
		    Random.seed!(1)
		    if matrix == "sprand_sdd!/1e6"
		        A = fem2d(1000000)
		    elseif matrix == "sprand_sdd!/1e7"
		        A = fem2d(10000000)
		    elseif matrix == "sprand/1e6"
		        n = 1000000
		        A = sprand(n, n, 5/n)
		    elseif matrix == "sprand/1e7"
		        n = 10000000
		        A = sprand(n, n, 5/n)
		    else
			    A = @suppress matrixdepot(matrix)
			end
			b = A * ones(A.n)
		# SuiteSparse/MatrixMarket including rhs
		else
			md = @suppress mdopen(matrix)
			A = md.A
			b = md.b[:, -n]
		end
	else
		Random.seed!(1)
		# myrand.jl
		if matrix == "fem2d"
			A = fem2d(n)
		# MatrixDepot matrix generator
		else
			A = @suppress matrixdepot(matrix, n)
		end
		b = A * ones(A.n)
	end
	transA = transpose(sparse(transpose(A)))
	return A, transA, b
end

# Get benchmark results from CSV file and return as arrays
function get_results(filename::String, pos::Int)
    # Read in CSV file as dataframe
    df = DataFrame!(CSV.File(filename, header=true))
    # Select nthreads/matsize and runtime columns
    select!(df, [1, pos+1])
    # Convert columns to arrays
    nthreads_or_matsize = convert(Vector{Int}, df[:, 1])
    seconds = convert(Vector{Float64}, df[:, 2])
    # Return nthreads/matsize and runtime
    nthreads_or_matsize, seconds
end

max_stars_repo_name: jkrch/JuliaParGMRES-BenchCLI

max_issues_repo_name: jkrch/JuliaParGMRES-BenchCLI

max_forks_repo_name: jkrch/JuliaParGMRES-BenchCLI


text:
import Base.method_argnames
# Long form, as in display(ws) or REPL: ws enter
function Base.show(io::IO, ::MIME"text/plain", ws::WebSocket{T}) where T
    ioc = IOContext(io, :wslog => true)
    print(ioc, "WebSocket{", nameof(T), "}(", ws.server ? "server, " : "client, ")
    _show(ioc, ws.state)
    print(ioc, "): ")
    _show(ioc, ws.socket)
end
# Short form with state, as in print(stdout, ws)
function Base.show(io::IO, ws::WebSocket{T}) where T
    ioc = IOContext(io, :compact=>true, :wslog => true)
    if T == TCPSocket
        print(ioc, "WebSocket(")
    else
        print(ioc, "WebSocket{", nameof(T), "}(")
    end
    print(ioc, ws.server ? "server, " : "client, ")
    _show(ioc, ws.state)
    print(ioc,  ")")
end

# The default Juno / Atom display works nicely with standard output
Base.show(io::IO, ::MIME"application/prs.juno.inline", ws::WebSocket) = Base.show(io, ws)

# A Base.show method is already defined by @enum
function _show(io::IO, state::ReadyState)
    kwargs, msg = _uv_status_tuple(state)
    printstyled(io, msg; kwargs...)
end

function _show(io::IO, stream::Base.LibuvStream)
    # To avoid accidental type piracy, a double check:
    if !get(IOContext(io), :wslog, false)
        show(io, stream)
    else
        kwargs, msg = _uv_status_tuple(stream)
        printstyled(io, msg; kwargs...)
        if !(stream isa HTTP.Servers.UDPSocket)
            nba = bytesavailable(stream.buffer)
            nba > 0 && print(io, ", ", nba, " bytes")
            nothing
        end
    end
end
function _show(io::IO, stream::IOStream)
    # To avoid accidental type piracy, a double check:
    if !get(IOContext(io), :wslog, false)
        show(io, stream)
    else
        kwargs, msg = _uv_status_tuple(stream)
        printstyled(io, msg; kwargs...)
    end
end
function _show(io::IO, buf::Base.GenericIOBuffer)
    # To avoid accidental type piracy, a double check:
    if !get(IOContext(io), :wslog, false)
        show(io, buf)
    else
        kwargs, msg = _uv_status_tuple(buf)
        printstyled(io, msg; kwargs...)
        nba = buf.size
        nba > 0 && print(io, ", ", nba, " bytes")
        nothing
    end
end

"For colorful printing"
function _uv_status_tuple(x)
    # adaption of base/stream.jl_uv_status_string
    s = x.status
    if x.handle == Base.C_NULL
        if s == Base.StatusClosed
            (color = :red,), "✘" #"closed"
        elseif s == Base.StatusUninit
            (color = :red,), "null"
        else
            (color = :red,), "invalid status"
        end
    elseif s == Base.StatusUninit
        (color = :blue,), "uninit"
    elseif s == Base.StatusInit
        (color = :yellow,), "init"
    elseif s == Base.StatusConnecting
        (color = :blue,), "connecting"
    elseif s == Base.StatusOpen
        (color = :green,), "✓"   # "open"
    elseif s == Base.StatusActive
        (color = :green,), "active"
    elseif s == Base.StatusPaused
        (color = :red,), "paused"
    elseif s == Base.StatusClosing
        (color = :blue,), "closing"
    elseif s == Base.StatusClosed
        (color = :red,), "✘" #"closed"
    elseif s == Base.StatusEOF
        (color = :yellow,), "eof"
    else
        (color = :red,), "invalid status"
    end
end
function _uv_status_tuple(bs::Union{Base.BufferStream, IOStream, Base.GenericIOBuffer})
    if isopen(bs)
        (color = :green,), "✓"   # "open"
    else
        (color = :red,), "✘" #"closed"
    end
end
function _uv_status_tuple(status::ReadyState)
    if status == CONNECTED
        (color = :green,), "CONNECTED"
    elseif status == CLOSING
        (color = :blue,), "CLOSING"
    elseif status == CLOSED
        (color = :red,), "CLOSED"
    end
end

### ServerWS
function Base.show(io::IO, sws::ServerWS)
    print(io, ServerWS, "(handler=")
    _show(io, sws.handler.func)
    print(io, ", wshandler=")
    _show(io, sws.wshandler.func)
    if sws.connection_count[] != 0
        print(io, ", connection_count=" * string(sws.connection_count[])  )
    end
    for dke in keys(default_options())
        if dke ∉ (:in, :out, :connection_count)
            dva = get(default_options(), dke, nothing)
            ava = getfield(sws, dke)
            if dva != ava
                # ServerWS field not default
                print(io, ", ")
                _showoptions(IOContext(io, :wslog=>true), sws)
                break
            end
        end
    end
    print(io, ")")
    if isready(sws.in)
        printstyled(io, ".in:", color= :yellow)
        print(io, sws.in, " ")
    end
    if isready(sws.out)
        printstyled(io, ".out:", color= :yellow)
        print(io, sws.out, " ")
    end
end

function _showoptions(io::IO, sws::ServerWS)
    fina = fieldnames(ServerWS)
    for field in fina
        if field ∉ (:handler, :wshandler, :in, :out, :connection_count)
            fiva = getfield(sws, field)
            print(io, field, "=")
            if fiva == nothing
                print(io, "nothing")
            else
                _show(io, fiva)
            end
            if field != last(fina)
                print(io, ", ")
            end
        end
    end
    nothing
end

_show(io::IO, x) = show(io, x)
function _show(io::IO, f::Function)
   m = methods(f)
   if length(m) > 1
       print(io, f, " has ", length(m), " methods: ")
       Base.show_method_table(io, m, 4, false)
   else
       method = first(m)
       argnames = join(method_argnames(method)[2:end], ", ")
       print(io, method.name, "(", argnames, ")")
   end
end


text:
const PATH_TO_MINIZINC = "C:/Program Files/MiniZinc/minizinc.exe"

function gendata(kx, ky, n, percentile_cutoff, rho=0)
    X = rand(1:kx, n)
    Y = rand(1:kx, n)
    Y = ifelse.(rand(Uniform(0,1), length(X)) .< rho, X, Y)
    A = Matrix(freqtable(X,Y))
    threshold = floor.(Int, percentile(vec(A),percentile_cutoff))
    independent_A = sum(A, dims=2) * sum(A, dims=1) / sum(A)
    observed_A = @. ifelse(A > threshold, A, missing)
    (obs=observed_A, gen=A, threshold=threshold, ind = independent_A)
end

function listoptions(data)
    locs = findall(ismissing, data.obs)
    open("listallsolutions.mzn", "w") do io
        for loc in locs
            println(io, "var 0..$(data.threshold): a$(loc[1])_$(loc[2]);")
        end
        for i in axes(data.obs, 1)
            print(io, "constraint ")
            l = []
            for j in axes(data.obs, 2)
                if ismissing(data.obs[i,j])
                    push!(l, "a$(i)_$j")
                else
                    push!(l, data.obs[i, j])
                end
            end
            println(io, join(l, "+"), "=", sum(data.gen[i, :]),";")
        end
        for j in axes(data.obs, 2)
            print(io, "constraint ")
            l = []
            for i in axes(data.obs, 1)
                if ismissing(data.obs[i,j])
                    push!(l, "a$(i)_$j")
                else
                    push!(l, data.obs[i, j])
                end
            end
            println(io, join(l, "+"), "=", sum(data.gen[:, j]),";")
        end
        println(io, "solve satisfy;")
        print(io, "output [")
        print(io, join(["show(a$(loc[1])_$(loc[2]))" for loc in locs], ",\" \","))
        println(io, "];")
    end
    df = DataFrame([Int[] for _ in 1:length(locs)],
         [Symbol("a", loc[1], "_", loc[2]) for loc in locs])
    res = read(`$PATH_TO_MINIZINC -a listallsolutions.mzn`, String)
    res_lines = split(res, "\n")
    for i in 1:2:length(res_lines)-2
        push!(df, parse.(Int, split(res_lines[i])))
    end
    return df
end

function compare(data, df)
    agg_df = by(stack(df, :), :variable, :value => mean)
    #closest_idx = argmin(sum(abs.(transpose(convert(Array, df)) .- agg_df.value_mean), dims = 2))[1]
    #agg_df.exact_case = vec(convert(Array,df[closest_idx,:]))
    #println("Closest to the mean is $closest_idx")
    agg_df.actual = [data.gen[parse.(Int, split(String(idx)[2:end], "_"))...] for idx in agg_df.variable]
    return agg_df
end

mse_compare(agg_df) = mean(@. (agg_df.value_mean - agg_df.actual) ^ 2)
mae_compare(agg_df) = mean(@. abs(agg_df.value_mean - agg_df.actual))

function method_exact(k, n, percentile_cutoff,rho; verbose=true)
    verbose && @info "generating data"
    data = gendata(k, k, n, percentile_cutoff,rho)
    verbose && @info "data has size $(size(data.gen)) with $(count(ismissing, data.obs)) missings"
    df = listoptions(data)
    verbose && @info "possible assignments: $(nrow(df))"
    agg_df = compare(data, df)
    agg_df.difference = agg_df.actual - agg_df.value_mean
    mae = mae_compare(agg_df)
    verbose && @info "mae: $mae"
    return (data=data, blanks=count(ismissing, data.obs), options=df, mae=mae, flatten_df = agg_df, vc = cramer_v(data.gen), fit = agg_df.value_mean)
end

max_stars_repo_name: Michal-Kot/FillMissingEntries

max_issues_repo_name: Michal-Kot/FillMissingEntries

max_forks_repo_name: Michal-Kot/FillMissingEntries


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE

module FEMBase

using LinearAlgebra

macro lintpragma(s)
end

import Base: getindex, setindex!, convert, length, size, isapprox,
             first, last, vec,
             ==, +, -, *, /, haskey, copy, push!, isempty, empty!,
             append!, read

using TimerOutputs
export @timeit, print_timer

import FEMBasis
import FEMQuad

include("fields.jl")
export interpolate, update!, DCTI, DCTV, DVTI, DVTV, CVTV, DVTId, DVTVd, field
include("types.jl")
include("sparse.jl")
include("elements.jl")
export element_info!
include("elements_lagrange.jl")
include("integrate.jl")

include("problems.jl")
export set_gdofs!, get_gdofs, get_formulation_type, add_element!, add_elements!

include("assembly.jl")
export assemble_prehook!, assemble_posthook!, assemble_elements!

include("solvers.jl")
export LinearSystem, AbstractLinearSystemSolver, solve!

include("analysis.jl")
export AbstractAnalysis, Analysis, add_problems!, get_problems, run!,
       AbstractResultsWriter, add_results_writer!, get_results_writers,
       write_results!, get_problem

export FieldProblem, BoundaryProblem, Problem, Element, Assembly
export add!, group_by_element_type,
       get_unknown_field_name, get_unknown_field_dimension,
       get_integration_points, initialize!, assemble!
export SparseMatrixCOO, SparseVectorCOO, Node, BasisInfo,
       IP, AbstractProblem, IntegrationPoint, AbstractField
export is_field_problem, is_boundary_problem, get_elements,
       get_connectivity,
       get_parent_field_name, get_reference_coordinates,
       get_reference_element_coordinates,
       get_assembly, get_nonzero_rows, get_nonzero_columns,
       eval_basis!, get_basis, get_dbasis, grad!,
       assemble_mass_matrix!, get_local_coordinates, inside,
       get_element_type, filter_by_element_type, get_element_id,
       resize_sparse, resize_sparsevec

include("test.jl")

function eliminate_boundary_conditions! end
export eliminate_boundary_conditions!

function add_slave_elements! end
export add_slave_elements!

function get_slave_elements end
export get_slave_elements

function add_master_elements! end
export add_master_elements!

function get_master_elements end
export get_master_elements

using FEMBasis
export Poi1, Seg2, Seg3, Tri3, Tri6, Tri7, Quad4, Quad8, Quad9,
       Tet4, Tet10, Pyr5, Wedge6, Wedge15, Hex8, Hex20, Hex27

include("deprecated.jl")

end


text:
_childsort(x::Union{Tuple, Vector, OrderedDict}) = x
_childsort(x::Dict) = sort!(OrderedDict(x))
_childsort(x::PairVec) = sort(x, by=first)
function _childsort(x::NamedTuple{T}) where T
    ks = tuple(sort(collect(T))...)
    NamedTuple{ks}(x[k] for k in ks)
end
_children_sorted(n) = _childsort(children(n))

_iter(x::PairVec) = last.(x)
_iter(x) = values(x)

_isnothing_iter(x) = isnothing.(_iter(x))

function _ith_child(m::Union{Tuple, Vector, NamedTuple}, i::Integer)
    return m[i]
end
_ith_child(m::OrderedDict, i::Integer) = _ith_child(collect(m), i)
function _ith_child(m::PairVec, i::Integer)
    return last(m[i])
end

_indexed(::Union{Vector, Tuple}) = true
_indexed(_) = false
_named(::Union{NamedTuple, Dict, OrderedDict, PairVec}) = true
_named(_) = false

function _children_pairs(ts, complete::Bool)
    chss1 = [_children_sorted(t) for t in ts if !(isnothing(t) || isleaf(t))]
    chss2 = [isnothing(t) || isleaf(t) ? nothing : _children_sorted(t) for t in ts]
    isempty(chss1) && return chss1
    # type stability here is impossible
    if all(isa.(chss2, PairVec))
        _children_pairs_aux(PairVec[], chss2, complete)
    else
        _children_pairs_aux(chss1, chss2, complete)
    end
end

# for performance reasons, this should always return the same type as input
function _children_pairs_aux(::Vector{<:Tuple}, chss, complete)
    chss = [isnothing(chs) ? tuple() : chs for chs in chss]
    ml = complete ? maximum(length.(chss)) : minimum(length.(chss))
    tuple(
          (tuple((i <= length(chs) ? chs[i] : nothing for chs in chss)...) for i in 1:ml)...
         )
end

function _children_pairs_aux(::Vector{<:Vector}, chss, complete)
    if any(isa.(chss, PairVec))
        throw_types(chss)
    end
    chss = [isnothing(chs) ? [] : collect(chs) for chs in chss]
    ml = complete ? maximum(length.(chss)) : minimum(length.(chss))
    [tuple((i <= length(chs) ? chs[i] : nothing for chs in chss)...) for i in 1:ml]
end

function _children_pairs_aux(::Vector{<:NamedTuple}, chss, complete)
    chss = [isnothing(chs) ? NamedTuple() : chs for chs in chss]
    ks = complete ? union(keys.(chss)...) : intersect(keys.(chss)...)
    (; (Symbol(k) => tuple(
                (k in keys(chss[i]) ? chss[i][k] : nothing for i in eachindex(chss))...
               ) for k in sort(ks))...)
end

function _children_pairs_aux(::Vector{<:OrderedDict}, chss, complete)
    chss = [isnothing(chs) ? OrderedDict() : chs for chs in chss]
    ks = complete ? union(collect.(keys.(chss))...) : intersect(collect.(keys.(chss))...)
    OrderedDict(k => tuple((k in keys(chss[i]) ? chss[i][k] : nothing for i in eachindex(chss))...)
        for k in sort(ks |> collect))
end

function _children_pairs_aux(::Vector{<:PairVec}, chss, complete)
    ks = [isnothing(chs) || length(Dict(chs)) == length(chs) for chs in chss]
    all(ks) || throw(ArgumentError("Two or more identical keys in children represented as a vector of pairs"))
    chss = [isnothing(chs) ? OrderedDict() : OrderedDict(chs) for chs in chss]
    _children_pairs_aux(chss, chss, complete) |> collect
end

# fallback definition, this should be avoided as much as possible
function _children_pairs_aux(chss1, chss2, complete)
    if all(_indexed.(chss1))
        chss2 = [isnothing(chs) ? chs : chs |> collect for chs in chss2]
        return _children_pairs_aux(Vector[], chss2, complete)
    elseif all(_named.(chss1))
        chss2 = [isnothing(chs) || chs isa PairVec ? chs : chs |> pairs |> collect for chs in chss2]
        return _children_pairs_aux(PairVec[], chss2, complete)
    end
    throw_types(chss1)
end

function throw_types(chss)
    throw(ArgumentError(
                        "Incompatible children types: " * join([typeof(chs) for chs in chss], ", ", " and ")
                       ))
end


text:
function visualize(dataset::AbstractStaticVTKData; color="", RGB=false, component=-1, opacity=1., window_size=(1000,600), ncolors=100, legend=false, legendtitle=color, scale_factor=1, representation=:simple)
    renderwindow, interactor = make_render_window_and_interactor(dataset, color, RGB, component, opacity, window_size, ncolors, legend, legendtitle, scale_factor, representation)
    renderwindow.Render()
    interactor.Start()
    renderwindow.Finalize()
end

function visualize_3ds(filename::AbstractString; window_size=(1000,600), background=(0.1, 0.2, 0.4), camera_position=(0, 1, 0), camera_focal=(0, 0, 0), camera_up=(0,0,1), camera_dolly=1.4)
    importer = vtk.vtk3DSImporter()
    importer.ComputeNormalsOn()
    importer.SetFileName(filename)
    importer.Read()
    renderwindow = importer.GetRenderWindow()
    interactor = vtk.vtkRenderWindowInteractor()
    interactor.SetRenderWindow(renderwindow)
    renderwindow.SetSize(window_size...);

    renderer = importer.GetRenderer()
    renderer.SetBackground(background);

    camera = renderer.GetActiveCamera()
    camera.SetFocalPoint(camera_focal) #Center
    camera.Azimuth(37.5)
    camera.Elevation(30)
    camera.Dolly(camera_dolly)

    renderer.ResetCamera();
    renderer.ResetCameraClippingRange();

    axes = vtk.vtkAxesActor()
    widget = vtk.vtkOrientationMarkerWidget()
    widget.SetOutlineColor(0.9300,0.5700,0.1300)
    widget.SetOrientationMarker(axes)
    widget.SetInteractor(interactor)
    widget.SetViewport(0., 0., 0.2, 0.2)
    widget.SetEnabled(1)
    widget.InteractiveOn()

    renderwindow.Render()
    interactor.Start();
end

function write_x3d(dataset::AbstractStaticVTKData, filepath_noext::String; color="", RGB=false, component=-1, opacity=1., window_size=(1000,600), ncolors=100, background=(.1, .2, .3), legend=false, legendtitle=color, scale_factor=1, representation=:simple)
    polydata = VTKPolyData(dataset)
    if color != "" && legend
        save_legend(polydata, string(filepath_noext, "_legend.png"), color=color, component=component, window_size=window_size, ncolors=ncolors, opacity=opacity, background=background, legendtitle=legendtitle)
    end

    renderwindow, interactor = make_render_window_and_interactor(polydata, color, RGB, component, opacity, window_size, ncolors, legend, legendtitle, scale_factor, representation)
    renderwindow.OffScreenRenderingOn()
    renderwindow.Render()

    writer = vtk.vtkX3DExporter()
    writer.SetFileName(string(filepath_noext, ".x3d"))
    writer.SetRenderWindow(renderwindow)
    writer.Update()

    renderwindow.Finalize()
end

function write_ply(dataset, filepath_noext; color="", component=-1, opacity=1., window_size=(800,600), ncolors=100, background=(.1, .2, .3), legend=false, legendtitle=color, scale_factor=1, representation=:simple)
    if representation == :wireframe
        throw("write_ply does not support wireframe representation.")
    end

    if color in keys(dataset.point_data)
        point_color = true
    elseif color in keys(dataset.cell_data)
        point_color = false
    end

    polydata = VTKPolyData(dataset)
    vtkobject = PyVTK(polydata)

    if color != ""
        if legend
            save_legend(polydata, string(filepath_noext, "_legend.png"), color=color, component=component, window_size=window_size, ncolors=ncolors, opacity=opacity, background=background, legendtitle=legendtitle)
        end
        jl_mapped_colors = get_jl_mapped_colors(vtkobject, polydata, color, component, opacity)[1]
    else
        jl_mapped_colors = fill(UInt8(255), (4, num_of_cells(polydata)))
    end

    if point_color
        mapped_colors = julia_to_vtk(jl_mapped_colors)
        mapped_colors.SetName("**_$(color)_**")
        mapped_colors.SetNumberOfComponents(size(jl_mapped_colors, 1))
        vtkobject.GetPointData().SetScalars(mapped_colors)
    else
        mapped_colors = julia_to_vtk(jl_mapped_colors)
        mapped_colors.SetName("**_$(color)_**")
        mapped_colors.SetNumberOfComponents(size(jl_mapped_colors, 1))
        vtkobject.GetCellData().SetScalars(mapped_colors)
    end

    if representation == :glyph
        vtkobject.GetPointData().SetActiveVectors(color)
        arrow_source = vtk.vtkArrowSource()
        glyph_filter = vtk.vtkGlyph3D()
        glyph_filter.SetSourceConnection(arrow_source.GetOutputPort())
        glyph_filter.SetScaleFactor(0.005*scale_factor)
        glyph_filter.SetColorModeToColorByScalar()
        glyph_filter.SetInputData(vtkobject)
        glyph_filter.OrientOn()
        glyph_filter.SetVectorModeToUseVector()
        glyph_filter.Update()
        writer = vtk.vtkPLYWriter()
        writer.SetFileName(string(filepath_noext,".ply"))
        writer.SetInputData(glyph_filter.GetOutput())
        writer.SetArrayName("**_$(color)_**")
        writer.Update()
    else
        writer = vtk.vtkPLYWriter()
        writer.SetFileName(string(filepath_noext,".ply"))
        writer.SetInputData(vtkobject)
        writer.SetArrayName("**_$(color)_**")
        writer.Update()
    end
end

function get_jl_mapped_colors(vtkobject, dataset::T, color_variable_name, component, opacity) where {T<:AbstractVTKSimpleData}
    if color_variable_name in keys(dataset.point_data)
        point_color = true
        _size = num_of_points(dataset)
        _color_variable = dataset.point_data[color_variable_name]
        _var_dim = var_dim(dataset, color_variable_name, "Point")
    elseif color_variable_name in keys(dataset.cell_data)
        point_color = false
        _size = num_of_cells(dataset)
        _color_variable = dataset.cell_data[color_variable_name]
        _var_dim = var_dim(dataset, color_variable_name, "Cell")
    end

    if T <: AbstractVTKStructuredData
        if _var_dim == 1
            color_variable = reshape(_color_variable, (_size,))
        else
            color_variable = reshape(_color_variable, (_var_dim, _size))
        end
    else
        color_variable = _color_variable
    end

    if component == -1 || component == 1 && _var_dim == 1
        if _var_dim == 1
            _values = color_variable
            cmin = min(_values...)
            cmax = max(_values...)
        else
            _values = [norm(color_variable[:,i]) for i in 1:_size]
            cmin = min(_values...)
            cmax = max(_values...)
        end
    else
        if 1 <= component <= _var_dim
            _values = color_variable[component,:]
            cmin = min(_values...)
            cmax = max(_values...)
        else
            throw("Cannot use component $component of a $color. $color only has $_var_dim components.")
        end
    end
    jl_mapped_colors = zeros(UInt8, 4, _size)
    if !(cmax ≈ cmin)
        scaled_values = (_values .- cmin) ./ (cmax - cmin)
        function fill_jl_mapped_colors(i)
            jl_mapped_colors[:,i] = UInt8[round.(color_map(scaled_values[i]) .* 255) ; round(255*opacity)]
        end
        map(fill_jl_mapped_colors, 1:_size)
    end
    return jl_mapped_colors, cmin, cmax
end

function save_legend(dataset::AbstractStaticVTKData, filepath_noext::String; color="", component=-1, window_size=(1000,600), ncolors=100, opacity=1., background=(.1, .2, .3), legendtitle=color)
    if color != ""
        renderwindow = make_legend_render_window(dataset, color, background, ncolors, legendtitle, component, opacity, window_size)
        imagefilter = vtk.vtkWindowToImageFilter()
        imagefilter.SetInput(renderwindow)
        imagefilter.SetViewport(0.8, 0, 1., 1.)

        imagewriter = vtk.vtkPNGWriter()
        imagewriter.SetFileName(string(filepath_noext,".png"))
        imagewriter.SetInputConnection(imagefilter.GetOutputPort())
        imagewriter.Write()

        renderwindow.Finalize()
    end
end

function make_legend_render_window(dataset, color, background, ncolors, legendtitle, component, opacity, window_size)
    if color in keys(dataset.point_data)
        point_color = true
        _size = num_of_points(dataset)
        _color_variable = dataset.point_data[color]
        _var_dim = var_dim(dataset, color, "Point")
    elseif color in keys(dataset.cell_data)
        point_color = false
        _size = num_of_cells(dataset)
        _color_variable = dataset.cell_data[color]
        _var_dim = var_dim(dataset, color, "Cell")
    end
    if component == -1 || component == 1 && _var_dim == 1
        if _var_dim == 1
            _values = reshape(_color_variable, (_size))
            cmin = min(_values...)
            cmax = max(_values...)
        else
            color_variable = reshape(_color_variable, (_var_dim, _size))
            _values = [norm(color_variable[:,i]) for i in 1:_size]
            cmin = min(_values...)
            cmax = max(_values...)
        end
    else
        if 1 <= component <= _var_dim
            color_variable = reshape(_color_variable, (_var_dim, _size))
            _values = color_variable[component,:]
            cmin = min(_values...)
            cmax = max(_values...)
        else
            throw("Cannot use component $component of a $color. $color only has $_var_dim components.")
        end
    end

    renderer = vtk.vtkRenderer()
    renderer.SetBackground(background...)
    renderwindow = vtk.vtkRenderWindow()
    renderwindow.SetSize(window_size...)
    renderwindow.AddRenderer(renderer)
    interactor = vtk.vtkRenderWindowInteractor()
    interactor.SetRenderWindow(renderwindow)

    add_scalar_bar!(interactor, cmin, cmax, ncolors, opacity, legendtitle)

    renderwindow.OffScreenRenderingOn()
    renderwindow.Render()

    return renderwindow
end

function add_scalar_bar!(interactor, cmin, cmax, ncolors, opacity, legendtitle)
    jl_legend_colors = zeros(UInt8, 4, ncolors)
    function fill_jl_legend_colors(i)
        jl_legend_colors[:,i] = UInt8[round.(color_map(0+(i-1)*1/ncolors).*255); round(255*opacity)]
    end
    map(fill_jl_legend_colors, 1:ncolors)
    legend_colors = julia_to_vtk(jl_legend_colors)

    lut = vtk.vtkLookupTable()
    lut.SetTableRange(cmin, cmax)
    lut.SetNumberOfColors(ncolors)
    lut.SetTable(legend_colors)

    scalar_bar = vtk.vtkScalarBarActor()
    scalar_bar.SetNumberOfLabels(8)
    scalar_bar.SetLabelFormat("%+#6.2e")
    scalar_bar.SetLookupTable(lut)
    scalar_bar.GetLabelTextProperty().SetFontFamilyToCourier()
    scalar_bar.GetLabelTextProperty().SetJustificationToRight()
    scalar_bar.GetLabelTextProperty().SetVerticalJustificationToBottom()
    scalar_bar.GetLabelTextProperty().BoldOff()
    scalar_bar.GetLabelTextProperty().ItalicOff()
    scalar_bar.GetLabelTextProperty().ShadowOff()
    scalar_bar.GetLabelTextProperty().SetColor(1, 1, 1)
    scalar_bar.SetTitle(legendtitle)
    scalar_bar.GetTitleTextProperty().SetFontFamilyToArial()
    scalar_bar.GetTitleTextProperty().ItalicOff()
    scalar_bar.GetTitleTextProperty().ShadowOff()

    # create the scalar bar widget
    scalar_bar_widget = vtk.vtkScalarBarWidget()
    scalar_bar_widget.SetInteractor(interactor)
    scalar_bar_widget.SetScalarBarActor(scalar_bar)
    scalar_bar_widget.On()
end

function color_map(v)
    if 0.0 <= v < 0.25
        b = 1.0
        g = v/0.25
        r = 0.0
    elseif 0.25 <= v < 0.5
        b = 1.0-(v-0.25)/0.25
        g = 1.0
        r = 0.0
    elseif 0.5 <= v < 0.75
        b = 0.0
        g = 1.0
        r = (v-0.5)/0.25
    else
        b = 0.0
        g = 1.0-(v-0.75)/0.25
        r = 1.0
    end
    return [r,g,b]
end

function setup_camera(camera, dataset)
    _bb = bb(dataset)
    calib = max([(_bb[i+1]-_bb[i]) for i in 1:2:length(_bb)]...)
    camera.SetFocalPoint(pseudo_center(dataset)...) #Center
    camera.Azimuth(37.5)
    camera.Elevation(30)
    camera.Dolly(0.2/calib)
end

function make_render_window_and_interactor(dataset::T, color,
    RGB, component, opacity, window_size, ncolors, legend, legendtitle, scale_factor,
    representation) where {T<:AbstractStaticVTKData}

    vtkobject = PyVTK(dataset)
    cmin = cmax = nothing
    if color in keys(dataset.point_data)
        if RGB
            jl_mapped_colors = UInt8.(round.(reshape(dataset.point_data[color],
                (size(dataset.point_data[color], 1), num_of_points(dataset)))))
            mapped_colors = julia_to_vtk(jl_mapped_colors)
            mapped_colors.SetNumberOfComponents(size(dataset.point_data[color], 1))
            vtkobject.GetPointData().SetScalars(mapped_colors)
        else
            jl_mapped_colors, cmin, cmax = get_jl_mapped_colors(vtkobject, dataset, color, component, opacity)
            if !(cmax ≈ cmin)
                mapped_colors = julia_to_vtk(jl_mapped_colors)
                vtkobject.GetPointData().SetScalars(mapped_colors)
            end
        end
    elseif color in keys(dataset.cell_data)
        if RGB
            jl_mapped_colors = UInt8.(round.(reshape(dataset.cell_data[color],
                (size(dataset.cell_data[color], 1), num_of_cells(dataset)))))
            mapped_colors = julia_to_vtk(jl_mapped_colors)
            mapped_colors.SetNumberOfComponents(size(dataset.cell_data[color], 1))
            vtkobject.GetCellData().SetScalars(mapped_colors)
        else
            jl_mapped_colors, cmin, cmax = get_jl_mapped_colors(vtkobject, dataset, color, component, opacity)
            if !(cmax ≈ cmin)
                mapped_colors = julia_to_vtk(jl_mapped_colors)
                vtkobject.GetCellData().SetScalars(mapped_colors)
            end
        end
    end

    if representation == :glyph
        vtkobject.GetPointData().SetActiveVectors(color)
        arrow_source = vtk.vtkArrowSource()
        glyph_filter = vtk.vtkGlyph3D()
        glyph_filter.SetSourceConnection(arrow_source.GetOutputPort())
        glyph_filter.SetScaleFactor(0.005*scale_factor)
        glyph_filter.SetColorModeToColorByScalar()
        glyph_filter.SetInputData(vtkobject)
        glyph_filter.OrientOn()
        glyph_filter.SetVectorModeToUseVector()
        glyph_filter.Update()
        mapper = vtk.vtkDataSetMapper()
        mapper.SetInputConnection(glyph_filter.GetOutputPort())
    else
        mapper = vtk.vtkDataSetMapper()
        mapper.SetInputData(vtkobject)
    end

    actor = vtk.vtkActor()
    actor.SetMapper(mapper)
    if color == ""
        actor.GetProperty().SetColor(1.0, 1.0, 1.0)
    end
    if cmin !== nothing && cmax !== nothing && cmin ≈ cmax
        actor.GetProperty().SetColor(0.0, 0.0, 1.0)
    end
    if representation == :wireframe
        property = vtk.vtkProperty()
        property.SetRepresentationToWireframe()
        actor.SetProperty(property)
    end

    renderer = vtk.vtkRenderer()
    renderer.SetBackground(0.1, 0.2, 0.3)
    renderer.AddActor(actor)

    camera = vtk.vtkCamera()
    setup_camera(camera, dataset)
    renderer.SetActiveCamera(camera)

    renderwindow = vtk.vtkRenderWindow()
    renderwindow.SetSize(window_size...)
    renderwindow.AddRenderer(renderer)
    interactor = vtk.vtkRenderWindowInteractor()
    interactor.SetRenderWindow(renderwindow)

    if legend && color != "" && !RGB && !(cmin ≈ cmax)
        add_scalar_bar!(interactor, cmin, cmax, ncolors, opacity, legendtitle)
    end

    axes = vtk.vtkAxesActor()
    widget = vtk.vtkOrientationMarkerWidget()
    widget.SetOutlineColor(0.9300,0.5700,0.1300)
    widget.SetOrientationMarker(axes)
    widget.SetInteractor(interactor)
    widget.SetViewport(0., 0., 0.2, 0.2)
    widget.SetEnabled(1)
    widget.InteractiveOn()

    return renderwindow, interactor
end


text:
# Load modules
using Optim

x = -10:0.01:10
len_x = length(x)

options = Optim.Options(g_abstol=1e-8, g_reltol=1e-8, iterations=10^5, show_trace=false, store_trace=false);

# Gaussian 1 peak
seed1 = [[0.0], [15.0, 0.0, sqrt(0.2)]]
# We generate some data using the Gaussian model
y_data = Utils.gaussian(x, seed1) .+ randn(len_x)
sol = fit_curve_model(:gaussian, x, y_data, seed1; options=options)
# Gaussian 2 peaks
seed2 = [[1.0], [15.0, -1.0, sqrt(0.2)], [10.0, 1.5, sqrt(0.1)]]
# We generate some data using the Gaussian model
y_data = Utils.gaussian(x, seed2) .+ randn(len_x)
sol = fit_curve_model(:gaussian, x, y_data, seed2; options=options)

# Lorentzian 1 peak
seed1 = [[0.0], [15.0, 0.0, 0.5]]
# We generate some data using the Lorentzian model
y_data = Utils.lorentzian(x, seed1) .+ randn(len_x)
sol = fit_curve_model(:lorentzian, x, y_data, seed1; options=options)
# Lorentzian 2 peaks
seed2 = [[0.0], [15.0, -1.0, 0.5], [30.0, 1.0, 0.5]]
# We generate some data using the Lorentzian model
y_data = Utils.lorentzian(x, seed2) .+ randn(len_x)
sol = fit_curve_model(:lorentzian, x, y_data, seed2; options=options)

# Voigtian 1 peak
seed1 = [[0.0], [18.0, 0.0, 0.0, 1.53]]
# We generate some data using the Voigtian model
y_data = Utils.voigtian(x, seed1) .+ randn(len_x)
sol = fit_curve_model(:voigtian, x, y_data, seed1; options=options)
# Voigtian 2 peaks
seed2 = [[0.0], [10.0, -1.6, 0.0, 1.53], [24.0, 5.0, 0.0, 1.53]]
# We generate some data using the Voigtian model
y_data = Utils.voigtian(x, seed2) .+ randn(len_x)
sol = fit_curve_model(:voigtian, x, y_data, seed2; options=options)


text:
## Copyright (c) 2015, Júlio Hoffimann Mendes <juliohm@stanford.edu>
##
## Permission to use, copy, modify, and/or distribute this software for any
## purpose with or without fee is hereby granted, provided that the above
## copyright notice and this permission notice appear in all copies.
##
## THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
## WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
## MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
## ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
## WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
## ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
## OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

function meanvoxreuse(training_image::AbstractArray,
                      tplsizex::Integer, tplsizey::Integer, tplsizez::Integer;
                      overlapx=1/6, overlapy=1/6, overlapz=1/6,
                      nreal=10, cut=:boykov, categorical=false, gpu=false)

  # calculate the overlap from given percentage
  ovx = ceil(Int, overlapx * tplsizex)
  ovy = ceil(Int, overlapy * tplsizey)
  ovz = ceil(Int, overlapz * tplsizez)

  # elementary raster path
  ntilex = ovx > 1 ? 2 : 1
  ntiley = ovy > 1 ? 2 : 1
  ntilez = ovz > 1 ? 2 : 1

  # simulation grid dimensions
  gridsizex = ntilex * (tplsizex - ovx) + ovx
  gridsizey = ntiley * (tplsizey - ovy) + ovy
  gridsizez = ntilez * (tplsizez - ovz) + ovz

  _, _, voxs = iqsim(training_image, tplsizex, tplsizey, tplsizez,
                     gridsizex, gridsizey, gridsizez,
                     overlapx=overlapx, overlapy=overlapy, overlapz=overlapz,
                     nreal=nreal, cut=cut, gpu=gpu, categorical=categorical, debug=true)

  mean(voxs)
end


text:
using SymPy: N, Sym, diff, expand, solve, subs

"""
Axisymmetric Solov'ev equilibra in (R/R₀,Z/R₀,phi) coordinates.
Based on Cerfon & Freidberg, Physics of Plasmas 17, 032502, 2010.

Parameters:
 * `R₀`: position of magnetic axis
 * `B₀`: B-field at magnetic axis
 * `ϵ`:  inverse aspect ratio
 * `κ`:  elongation
 * `δ`:  triangularity
 * `a`:  free constant, determined to match a given beta value
"""
struct Solovev{T <: Number} <: AbstractSolovevEquilibrium
    name::String
    R₀::T
    B₀::T
    ϵ::T
    κ::T
    δ::T
    a::T
    c::Vector{T}

    function Solovev{T}(R₀::T, B₀::T, ϵ::T, κ::T, δ::T, a::T, c::Vector{T}) where T <: Number
        new("SolovevEquilibrium", R₀, B₀, ϵ, κ, δ, a, c)
    end
end

function Solovev(R₀::T, B₀::T, ϵ::T, κ::T, δ::T, a::T) where T <: Number

    n = 7

    x = [Sym("x" * string(i)) for i in 1:3]
    c = [Sym("c" * string(i)) for i in 1:n]

    ψ = (ψ₀(x,a) + c[1] * ψ₁(x)
                 + c[2] * ψ₂(x)
                 + c[3] * ψ₃(x)
                 + c[4] * ψ₄(x)
                 + c[5] * ψ₅(x)
                 + c[6] * ψ₆(x)
                 + c[7] * ψ₇(x) )

    eqs = [
        expand(subs(subs(ψ, x[1], 1+ϵ), x[2], 0)),
        expand(subs(subs(ψ, x[1], 1-ϵ), x[2], 0)),
        expand(subs(subs(ψ, x[1], 1-δ*ϵ), x[2], κ*ϵ)),
        expand(subs(subs(diff(ψ, x[1]), x[1], 1-δ*ϵ), x[2], κ*ϵ)),
        expand(subs(subs(diff(ψ, x[2], 2), x[1], 1+ϵ), x[2], 0) - (1 + asin(δ))^2 / (ϵ * κ^2) * subs(subs(diff(ψ, x[1]), x[1], 1+ϵ), x[2], 0)),
        expand(subs(subs(diff(ψ, x[2], 2), x[1], 1-ϵ), x[2], 0) + (1 - asin(δ))^2 / (ϵ * κ^2) * subs(subs(diff(ψ, x[1]), x[1], 1-ϵ), x[2], 0)),
        expand(subs(subs(diff(ψ, x[1], 2), x[1], 1-δ*ϵ), x[2], κ*ϵ) - κ / (ϵ * (1 - δ^2)) * subs(subs(diff(ψ, x[2]), x[1], 1-δ*ϵ), x[2], κ*ϵ))
    ]

    csym = solve(eqs, c)
    cnum = [N(csym[c[i]]) for i in 1:n]

    Solovev{T}(R₀, B₀, ϵ, κ, δ, a, cnum)
end

SolovevITER() = Solovev(6.2, 5.3, 0.32, 1.7, 0.33, -0.155)
SolovevNSTX() = Solovev(0.85, 0.3, 0.78, 2.0, 0.35, 1.0)

function Base.show(io::IO, equ::Solovev)
    print(io, "Solovev Equilibrium with\n")
    print(io, "  R₀ = ", equ.R₀, "\n")
    print(io, "  B₀ = ", equ.B₀, "\n")
    print(io, "  ϵ  = ", equ.ϵ, "\n")
    print(io, "  κ  = ", equ.κ, "\n")
    print(io, "  δ  = ", equ.δ, "\n")
    print(io, "  a  = ", equ.a)
end

function A₃(x::AbstractArray{T,1}, equ::Solovev) where {T <: Number}
    (ψ₀(x, equ.a) + equ.c[1] * ψ₁(x)
                  + equ.c[2] * ψ₂(x)
                  + equ.c[3] * ψ₃(x)
                  + equ.c[4] * ψ₄(x)
                  + equ.c[5] * ψ₅(x)
                  + equ.c[6] * ψ₆(x)
                  + equ.c[7] * ψ₇(x) )
end

macro solovev_equilibrium(R₀, B₀, ϵ, κ, δ, a)
    generate_equilibrium_code(Solovev(R₀, B₀, ϵ, κ, δ, a); output=false)
end