text:
"""
Compare an in-memory table to a table schema.
"""
module inmemory_table

using CategoricalArrays
using DataFrames
using Tables

using ..handle_validvalues
using ..types
using ..common

function compare(tableschema::TableSchema, indata, sorted_by_primarykey::Bool)
    # Init
    tablename     = tableschema.name
    outdata       = init_outdata(tableschema, size(indata, 1))
    issues_in     = init_issues(tableschema)  # Issues for indata
    issues_out    = init_issues(tableschema)  # Issues for outdata
    pk_colnames   = tableschema.primarykey
    pk_ncols      = length(pk_colnames)
    primarykey    = missings(String, pk_ncols)  # Stringified primary key
    pkvalues_in   = Set{String}()  # Values of the primary key
    pkvalues_out  = Set{String}()
    pk_issues_in  = pk_ncols == 1 ? issues_in[:columnissues][pk_colnames[1]]  : Dict{Symbol, Int}()
    pk_issues_out = pk_ncols == 1 ? issues_out[:columnissues][pk_colnames[1]] : Dict{Symbol, Int}()
    i_outdata     = 0
    nconstraints  = length(tableschema.intrarow_constraints)
    colname2colschema = tableschema.colname2colschema
    uniquevalues_in   = Dict(colname => Set{nonmissingtype(eltype(getproperty(indata, colname)))}() for (colname, colschema) in colname2colschema if colschema.isunique==true)
    uniquevalues_out  = Dict(colname => Set{colschema.datatype}() for (colname, colschema) in colname2colschema if colschema.isunique==true)

    # Row-level checks
    for inputrow in Tables.rows(indata)
        # Parse inputrow into outputrow according to ColumnSchema
        i_outdata += 1
        outputrow  = outdata[i_outdata, :]
        parserow!(outputrow, inputrow, colname2colschema)

        # Assess input row
        if pk_ncols == 1
            pk_n_missing   = pk_issues_in[:n_missing]    # Number of earlier rows (excluding the curent row) with missing primary key
            pk_n_notunique = pk_issues_in[:n_notunique]  # Number of earlier rows (excluding the curent row) with duplicated primary key
        end
        assess_row!(issues_in[:columnissues], inputrow, colname2colschema, uniquevalues_in)
        nconstraints > 0 && test_intrarow_constraints!(issues_in[:intrarow_constraints], tableschema, inputrow)
        if pk_ncols == 1
            pk_incomplete, pk_duplicated = assess_singlecolumn_primarykey!(issues_in, pk_issues_in, pk_n_missing, pk_n_notunique)
        else
            pk_incomplete, pk_duplicated = assess_multicolumn_primarykey!(issues_in, primarykey, pk_colnames, pkvalues_in, sorted_by_primarykey, inputrow)
        end

        # Assess output row
        if pk_ncols == 1
            pk_n_missing   = pk_issues_out[:n_missing]    # Number of earlier rows (excluding the curent row) with missing primary key
            pk_n_notunique = pk_issues_out[:n_notunique]  # Number of earlier rows (excluding the curent row) with duplicated primary key
        end
        assess_row_set_invalid_to_missing!(issues_out[:columnissues], outputrow, colname2colschema, uniquevalues_out)
        nconstraints > 0 && test_intrarow_constraints!(issues_out[:intrarow_constraints], tableschema, outputrow)
        if pk_ncols == 1
            pk_incomplete, pk_duplicated = assess_singlecolumn_primarykey!(issues_out, pk_issues_out, pk_n_missing, pk_n_notunique)
        else
            pk_incomplete, pk_duplicated = assess_multicolumn_primarykey!(issues_out, primarykey, pk_colnames, pkvalues_out, sorted_by_primarykey, outputrow)
        end
    end

    # Column-level checks
    for (colname, colschema) in colname2colschema
        !colschema.iscategorical && continue
        transform!(outdata, colname => categorical, renamecols=false)
    end
    datacols_match_schemacols!(issues_in, tableschema, Set(propertynames(indata)))  # By construction this issue doesn't exist for outdata
    compare_datatypes!(issues_in,  indata,  colname2colschema)
    compare_datatypes!(issues_out, outdata, colname2colschema)

    # Format result
    issues_in  = construct_issues_table(issues_in,  tableschema, i_outdata)
    issues_out = construct_issues_table(issues_out, tableschema, i_outdata)
    outdata, issues_in, issues_out
end

"""
Modified: issues.

Checks whether each column and its schema are both categorical or both not categorical, and whether they have the same data types.
"""
function compare_datatypes!(issues, table, colname2colschema)
    for (colname, colschema) in colname2colschema
        coldata    = getproperty(table, colname)
        data_eltyp = colschema.iscategorical ? eltype(levels(coldata)) : nonmissingtype(eltype(coldata))
        if data_eltyp != colschema.datatype  # Check data type matches that specified in the ColumnSchema
            issues[:columnissues][colname][:different_datatypes] = 1
        end
        if colschema.iscategorical && !(coldata isa CategoricalVector)  # Ensure categorical values
            issues[:columnissues][colname][:data_not_categorical] = 1
        end
        if !colschema.iscategorical && coldata isa CategoricalVector    # Ensure non-categorical values
            issues[:columnissues][colname][:data_is_categorical] = 1
        end
    end
end

"Sets invalid values to missing."
function assess_row_set_invalid_to_missing!(columnissues, outputrow, colname2colschema, uniquevalues)
    for (colname, colschema) in colname2colschema
        val = outputrow[colname]
        if ismissing(val)
            !colschema.isrequired && continue  # Checks are done before the function call to avoid unnecessary dict lookups
            assess_missing_value!(columnissues[colname])
        elseif value_is_valid(val, colschema.validvalues)
            !colschema.isunique && continue
            assess_nonmissing_value!(columnissues[colname], val, uniquevalues[colname])
        else
            outputrow[colname] = missing
            !colschema.isrequired && continue
            assess_missing_value!(columnissues[colname])
        end
    end
end

end


text:
using DWARF
@osx_only using MachO
@linux_only using ELF
import ObjFileBase: DebugSections

tag_filter(which_tag, arr) = filter(x->DWARF.tag(x) == which_tag, arr)
tag_filter(which_tag::Vector, arr) = filter(x->(DWARF.tag(x) in which_tag), arr)

immutable Subprogram
    name::AbstractString
    linkage_name::AbstractString
    low_pc::UInt
    high_pc::UInt
    fbreg::UInt
    variables::Dict
end

immutable CompilationUnit
    SPs::Vector{Subprogram}
    low_pc::UInt
    high_pc::UInt
end

function process_cus(debugoh)
    dbgs = debugsections(debugoh);
    strtab = load_strtab(dbgs.debug_str)
    process_cus(dbgs, strtab)
end

function process_SP(SP, strtab)
    sp_low = (-1 % UInt); sp_high = 0; fbreg = (-1 % UInt)
    linkage_name = ""; name = ""
    for at in DWARF.attributes(SP)
        tag = DWARF.tag(at)
        if tag == DWARF.DW_AT_low_pc
            sp_low = convert(UInt,at)
        elseif tag == DWARF.DW_AT_high_pc
            sp_high = convert(UInt,at)
        elseif tag == DWARF.DW_AT_linkage_name ||
               tag == DWARF.DW_AT_MIPS_linkage_name
            linkage_name = bytestring(at,strtab)
        elseif tag == DWARF.DW_AT_frame_base
            # Simple for now
            fbreg = at.content[1]-DWARF.DW_OP_reg0
        elseif tag == DWARF.DW_AT_name
            name = bytestring(at,strtab)
        end
    end
    variables = Dict{Symbol,Any}()
    for v in SP.children
        name = nothing
        value = 0
        for at in DWARF.attributes(v)
            tag = DWARF.tag(at)
            if tag == DWARF.DW_AT_location
                value = at
            elseif tag == DWARF.DW_AT_const_value
                value = (-1 % UInt)
            elseif tag == DWARF.DW_AT_name
                name = bytestring(at, strtab)
            end
        end
        if name !== nothing
            variables[symbol(name)] = value
        end
    end
    Subprogram(name, linkage_name, sp_low, sp_high, fbreg, variables)
end

function process_tree(tree, strtab)
    low = (-1 % UInt); high = 0
    for at in DWARF.attributes(tree)
        tag = DWARF.tag(at)
        if tag == DWARF.DW_AT_low_pc
            low = convert(UInt,at)
        elseif tag == DWARF.DW_AT_high_pc
            high = convert(UInt,at)
        end
    end
    subprograms = map(tag_filter(DWARF.DW_TAG_subprogram,tree.children)) do SP
        process_SP(SP, strtab)
    end
    CompilationUnit(subprograms, low, high)
end
process_tree(tree::DWARF.DIETreeRef) = process_tree(tree.tree, tree.strtab)

process_cus(dbgs::DebugSections, strtab) = map(x->process_tree(x,strtab), DIETrees(dbgs))

function sections_for_oh(oh)
    dorelocate = false
    if isa(oh,MachO.MachOHandle)
        dorelocate = readheader(oh).filetype == MachO.MH_OBJECT
        sects = MachO.Sections(first(filter(LoadCmds(oh)) do lc
            isa(lc.cmd, MachO.segment_commands) && segname(lc) == "__TEXT"
        end))
        sects
    elseif isa(oh,ELF.ELFHandle)
        sects = ELF.Sections(oh)
    else
        error("Object file type unsupported")
    end
    (dorelocate, sects)
end

ext: jl


text:
mutable struct IntegratedFactorsGibbsOperator <: OperatorXMLElement
    el::XMLOrNothing
    sampled_tree_likelihood::TraitLikelihoodXMLElement
    integrated_tree_likelihood::TraitLikelihoodXMLElement
    factor_likelihood::IntegratedFactorsXMLElement
    weight::Float64

    function IntegratedFactorsGibbsOperator(
                    sampled_tree_likelihood::TraitLikelihoodXMLElement,
                    tree_likelihood::TraitLikelihoodXMLElement,
                    factor_likelihood::IntegratedFactorsXMLElement,
                    weight::Float64)
        return new(nothing, sampled_tree_likelihood,
                    tree_likelihood, factor_likelihood, weight)
    end
end

function make_xml(ifg::IntegratedFactorsGibbsOperator)
    el = new_element(bn.INTEGRATED_FACTOR_GIBBS)
    set_attribute(el, bn.WEIGHT, ifg.weight)

    sampled_el = make_xml(ifg.sampled_tree_likelihood)
    tp_el = find_element(sampled_el, bn.TRAIT_PARAMETER)
    fac_mat_el = find_element(tp_el, bn.PARAMETER)
    add_ref_el(el, fac_mat_el)

    integrated_el = make_xml(ifg.integrated_tree_likelihood)
    add_ref_el(el, integrated_el)

    fac_el = make_xml(ifg.factor_likelihood)
    add_ref_el(el, fac_el)

    ifg.el = el
    return el
end


text:
@testset "make_chain" begin

    @inferred make_chain(2)

    chain = make_chain(2)

    genesis_block = chain[0]
    @test double_sha256(genesis_block) ==
        UInt256("000000000019d6689c085ae165831e934ff763ae46a2a6c172b3f1b60a8ce26f")

    @test length(chain) == 2
    @test length(chain[1:1]) == 1
    @test chain[1] isa Link

    @test chain[0].hash == double_sha256(chain[0])
    @test chain[1].hash == double_sha256(chain[1])

    @test BTCParser.get_file_pos(chain[0]) == 0x0
    @test BTCParser.get_file_num(chain[0]) == 0
    @test BTCParser.get_file_pos(chain[1]) == 0x125
    @test BTCParser.get_file_num(chain[1]) == 0

    @test double_sha256(Header(chain[0])) == double_sha256(chain[0])
    @test double_sha256(Header(chain[1])) == double_sha256(chain[1])
    @test double_sha256(Block(chain[0]))  == double_sha256(chain[0])
    @test double_sha256(Block(chain[1]))  == double_sha256(chain[1])

    @test Header(chain[1]).previous_hash == double_sha256(chain[0])

    # TODO: should these be == ? I think the reason that these are longer is
    # that make chain finishes the out of order blocks.
    chain100 = make_chain(chain, 100)
    @test length(chain100) >= 100
    @test chain100[0] == genesis_block

    chain4000 = make_chain(deepcopy(chain100), 4000)
    @test length(chain4000) >= 4000
    @test chain4000[0] == genesis_block
    @test chain100[end] == chain4000[99]

    chain8000 = make_chain(chain4000, 8000)
    @test length(chain8000) >= 8000
    @test chain8000[0] == genesis_block
    @test chain100[end] == chain8000[99]
    @test chain4000[end] == chain8000[3999]
end


text:
using PyCall
@pyimport qutip as qt
@pyimport qutip.qip.algorithms.qft as qft
@pyimport qutip.control.pulseoptim as cpo

function grape_qutip(Ut,Hd,Hc,N,t,n)
    # start point for the gate evolution
    U0 = qt.identity(N)
    res = cpo.optimize_pulse_unitary(Hd,Hc,U0,Ut,n,t,init_pulse_type="ZERO",fid_err_targ=1E-16,min_grad=1E-16,max_iter=1000)
    uf = res[:final_amps]
    Uf = res[:evo_full_final][:full]()
    return uf, Uf, res
end

function opt_hadamard_qutip(n)
    Hd = qt.sigmaz()
    Hc = [qt.sigmax()]
    t = 10.0
    Ut = qt.hadamard_transform(1)
    grape_qutip(Ut,Hd,Hc,2,t,n)
end

function opt_2Q_QFT_qutip(n)
    Sx = qt.sigmax()
    Sy = qt.sigmay()
    Sz = qt.sigmaz()
    Si = qt.identity(2)
    # Drift Hamiltonian and the (four) control Hamiltonians
    H_d = (qt.tensor(Sx, Sx) + qt.tensor(Sy, Sy) + qt.tensor(Sz, Sz))
    H_c = [qt.tensor(Sx, Si), qt.tensor(Sy, Si), qt.tensor(Si, Sx), qt.tensor(Si, Sy)]
    # Target for the gate evolution - Quantum Fourier Transform gate
    U_targ = qft.qft(2)
    # Evolution time
    t = 6.0
    grape_qutip(U_targ,H_d,H_c,4,t,n)
end

#"import qutip as qt; import qutip.qip.algorithms.qft as qft; import qutip.control.pulseoptim as cpo; Sx = qt.sigmax();Sy = qt.sigmay();Sz = qt.sigmaz();Si = qt.identity(2);H_d = (qt.tensor(Sx, Sx) + qt.tensor(Sy, Sy) + qt.tensor(Sz, Sz));H_c = [qt.tensor(Sx, Si), qt.tensor(Sy, Si), qt.tensor(Si, Sx), qt.tensor(Si, Sy)];U_0 = qt.identity(4);U_targ = qft.qft(2);" "cpo.optimize_pulse_unitary(H_d, H_c, U_0, U_targ, 100, 6, init_pulse_type='ZERO')"


text:
struct Prefix
    by::Union{Nothing, Vector{Symbol}}
    sort::Union{Nothing, Vector{Symbol}}
end

struct Command
    by::Union{Vector{Symbol}, Nothing}
    sort::Union{Vector{Symbol}, Nothing}
    command::String
    arguments::Union{Vector{Symbol},Union{Expr, Nothing}}
    filter::Union{Expr, Nothing}
    use::Union{String,Nothing}
    options::Union{Dict{String,Any}, Nothing}
end


text:
struct Interval{T<:Real}
    a::T
    b::T
    # Interval(a,b) = a > b ? error("out of order") : new(a,b)
end

function ∈(x::T, interval::Interval) where T<:Real
    if interval.a ≤ x ≤ interval.b
        return true
    else
        return false
    end
end

function ∩(I::Interval, J::Interval)
    min = maximum([I.a, J.a])
    max = minimum([I.b, J.b])

    return Interval(min, max)
end

mutable struct MapperInformations
    data_number
    original_coordinate
    original_subset
    clustered_subset
    filtered_coordinate
    interval_subset
end

function interval_test(Info, Intervals)
    M = length(Info)
    N = length(Intervals)

    X_subset = [[] for i ∈ 1:N]
    for i ∈ 1:M
        for j ∈ 1:N
            if Info[i].filtered_coordinate ∈ Intervals[j]
                Info[i].interval_subset = Intervals[j]
                Info[i].original_subset = push!(Info[i].original_subset, j)
                push!(X_subset[j], Info[i].original_coordinate)
            end
        end
    end

    return X_subset
end

function domain_clustering(Info, X_subset)

    X_subset_clustering = []

    for i ∈ 1:length(X_subset)
        # i=3
        temp = zeros(length(X_subset[i][1]), length(X_subset[i]))
        for j ∈ 1:length(X_subset[i])
            temp[:, j] = X_subset[i][j]
        end
        scatter(temp[1,:], temp[2,:])

        D = pairwise(Euclidean(), temp, temp)

        SL = hclust(D, linkage=:single)
        plot(SL)

        dbsacn_temp = dbscan(reshape(SL.heights, (1,length(SL.heights))), 0.2maximum(SL.heights))
        cutree_temp = cutree(SL, h=SL.heights[dbsacn_temp[1].core_indices[end]])

        X_subset_clustering_temp = [[] for j ∈ 1:maximum(cutree_temp)]

        for n ∈ 1:length(X_subset[i])
            clustering_index = cutree_temp[n]
            data_index = findfirst(==(X_subset[i][n]), original_X)
            push!(X_subset_clustering_temp[clustering_index], data_index)
            Info[data_index].clustered_subset = push!(Info[data_index].clustered_subset, (i, clustering_index))
        end

        push!(X_subset_clustering, X_subset_clustering_temp)

    end
    return X_subset_clustering
end

function coor_mean(X_subset_clustering, Info)
    center_ = [ [] for i ∈ 1:length(X_subset_clustering)]
    for i ∈ 1:length(X_subset_clustering)
        for j ∈ 1:length(X_subset_clustering[i])
            temp = [Info[v].original_coordinate for v ∈ X_subset_clustering[i][j]]
            center_vector = sum(temp)./length(temp)
            push!(center_[i], center_vector)
        end
    end

    return center_
end

function mapper(X, filter, interval, overlab)

    original_X = []
    for i in 1:size(X,2)
        push!(original_X, X[:,i])
    end
    # findfirst(==(Info[500].original_coordinate), original_X)

    Info = [MapperInformations(i, X[:, i], Int64[], Tuple{Int64, Int64}[], :None, :None) for i ∈ 1:length(X[1,:])]

    filtered_X = filter(X)
    filtered_coordinate.(Info, filtered_X)
    Info

    min_X = minimum(filtered_X)
    max_X = maximum(filtered_X)
    filter_range = Interval(min_X, max_X)

    length_ = (max_X - min_X)/interval
    ϵ = overlab*length_

    Intervals = [∩(Interval(min_X +(i-1)*length_ - ϵ, min_X +i*length_ + ϵ), filter_range) for i ∈ 1:interval]

    X_subset = interval_test(Info, Intervals)

    Info

    X_subset_clustering = domain_clustering(Info, X_subset)

    # function myshowall(io, x, limit = false)
    #     println(io, summary(x), ":")
    #     Base.print_matrix(IOContext(io, :limit => limit), x)
    # end
    # myshowall(Base.stdout, Info, false)

    # X_subset_clustering[1]
    # X_subset_clustering[2]
    # X_subset_clustering[3]
    # X_subset_clustering[4]
    # X_subset_clustering[5]

    link = []

    for i ∈ 1:length(X_subset)-1
        for j ∈ 1:length(X_subset_clustering[i])
            for k ∈ 1:length(X_subset_clustering[i+1])
                if length(∩(X_subset_clustering[i][j], X_subset_clustering[i+1][k])) > 0
                    push!(link, [[i,j], [i+1,k]])
                end
            end
        end
    end

    center_ = coor_mean(X_subset_clustering, Info)

    center_number = 0
    for i ∈ 1:length(X_subset)
        for j ∈ 1:length(X_subset_clustering[i])
            center_number +=1
        end
    end

    center_matrix = zeros(2, center_number)
    k_temp = 1
    for i ∈ 1:length(X_subset)
        for j ∈ 1:length(X_subset_clustering[i])
            center_matrix[:,k_temp] = center_[i][j]
            k_temp += 1
        end
    end

    zcolor_temp = filter(cat(Data, center_matrix, dims=2))

    zcolor = [[] for i ∈ 1:length(X_subset)]
    k_temp = 1
    for i ∈ 1:length(X_subset)
        for j ∈ 1:length(X_subset_clustering[i])
            push!(zcolor[i], zcolor_temp[size(Data)[2]+k_temp])
            k_temp +=1
        end
    end

    return link, X_subset_clustering, Info, center_, zcolor
end

function mapper_plot(link, X_subset_clustering, center_, zcolor, xlim, ylim)

    linked_center = zeros(length(link), 4)
    for i ∈ 1:length(link)
        linked_center[i, 1:2] = center_[link[i][1][1]][link[i][1][2]]
        linked_center[i, 3:4] = center_[link[i][2][1]][link[i][2][2]]
    end
    linked_center

    p = plot([linked_center[1,1], linked_center[1,3]], [linked_center[1,2], linked_center[1,4]], color=:black, xlim=xlim, ylim=ylim, label="")
    for i ∈ 2:length(link)
        plot!([linked_center[i,1], linked_center[i,3]], [linked_center[i,2], linked_center[i,4]], color=:black, label="")
    end

    for i ∈ 1:length(X_subset_clustering)
        for j ∈ 1:length(X_subset_clustering[i])
            scatter!([center_[i][j][1]], [center_[i][j][2]], markersize=25*(452.106/181)*sqrt(length(X_subset_clustering[i][j])/500), zcolor=zcolor[i][j], label="", color=:jet, clim=(0,Inf))
            length_temp = length(X_subset_clustering[i][j])
            annotate!([center_[i][j][1]], [center_[i][j][2]], "$length_temp")
        end
    end

    return p
end


text:
""" We store the info about the cloud composition in these types.
The general procedure to fetch the cloud composition at agiven location has two
steps: first we `probe` the cloud at a given location at find the local
properties relevant to the cloud coposition model (for example this may be the
droplet radius).  Then we call functions `density`, `radius`, `mie_*` to
obtain the particular parameter required for some computation.  The reason to
split this procedure into several stages is to avoid repetitive computations
to obtain e.g. the droplet radius.

A cloud-composition type must implement these methods:

 - probe: Return an object that stores the local composition state.  For example, it may return a Float64 with the
    radius of droplets at a given point.  It may also return nothing if we do not need to use any local info.
 - density: droplet density.
 - mie_qext: extinction parameter.
 - mie_g: asymmetry parameter.
 - mie_ω0: Single-scattering albedo.
"""

abstract type AbstractComposition end

const DEF_REFINDEX_FILE = "Hale.dat"

# Linear regression
linreg(x, y) = hcat(fill!(similar(x), 1), x) \ y
linreg1(x, y) = x \ y

# Fast computation of x^(-3/4)
power34(x) = 1. / sqrt(x) / sqrt(sqrt(x))

""" A type for a fixed-composition cloud.  It may be initiated with the standard constructor if you know
g, ω0, Qext or with fixedradius if you want to solve the Mie problem for a fixed radius.
"""
struct Homogeneous <: AbstractComposition
    n::Float64
    radius::Float64
    g::Float64
    ω0::Float64
    Qext::Float64
    νmiemax::Float64

    function Homogeneous(n, radius, g, ω0, Qext)
        new(n, radius, g, ω0, Qext, Qext * π * radius^2 * n)
    end
end

probe(comp::Homogeneous, r::Point) = nothing
density(comp::Homogeneous, r::Point, ::Nothing) = comp.n
radius(comp::Homogeneous, r::Point, ::Nothing) = comp.radius
mie_qext(comp::Homogeneous, r::Point, ::Nothing) = comp.Qext
mie_g(comp::Homogeneous, r::Point, ::Nothing) = comp.g
mie_ω0(comp::Homogeneous, r::Point, ::Nothing) = comp.ω0
νmiemax(comp::Homogeneous) = comp.νmiemax

"""
  Create a Homogeneous composition instance for a given wavelength `λ`and
  fixed droplet density `n` and droplet radius `r`.
"""
function Homogeneous(λ, n, radius; refindex=DEF_REFINDEX_FILE, verbose=true)
    interp = m_interpolator(refindex)
    m = interp(λ)

    g, ω0, Qext = MieParams.compute(radius, λ, m)
    verbose && @info "Mie parameters" g ω0 Qext

    Homogeneous(n, radius, g, ω0, Qext)
end

struct VariableNR{T} <: AbstractComposition
    nrfetch::T

    # The assymetry parameter is fit as g = mr / (r + r0)
    g0::Float64
    r0::Float64

    # The single-scattering albedo is fit as 1 - ω0 = ar
    a::Float64

    # The Extinction coefficient is fit as Qext = 2 + cr^-3/4
    c::Float64

    νmiemax::Float64
end

"""
    m = VariableNR(λ)

Compute the fit parameters to compute g, ω0 and Qext for an arbitrary radius
in the range 1 μm < r < 100 μm.
"""
function VariableNR(λ, nrfetch, νmiemax; refindex=DEF_REFINDEX_FILE)
    g0, r0, a, c = miefitparams(λ, refindex)
    VariableNR(nrfetch, g0, r0, a, c, νmiemax)
end

function VariableNR(λ, nrfetch, rmax, nmax; refindex=DEF_REFINDEX_FILE)
    g0, r0, a, c = miefitparams(λ, refindex)
    Qext_max = 2. + c * power34(rmax)
    νmiemax = Qext_max * π * rmax^2 * nmax
    VariableNR(nrfetch, g0, r0, a, c, νmiemax)
end

function miefitparams(λ, refindex)
    interp = m_interpolator(refindex)
    m = interp(λ)

    r = 1 * co.micro:1 * co.micro:100 * co.micro
    g, ω0, Qext = compute(r, λ, m; order=10000)

    # Linear regression always more robust so we fit r/g ~ r/m + r0 / m
    β = linreg(r, r ./ g)
    g0, r0 = 1 / β[2], β[1] / β[2]

    a = linreg1(r, 1 .- ω0)[1]
    c = linreg1(power34.(r), Qext .- 2)[1]

    (g0, r0, a, c)
end

probe(comp::VariableNR, r::Point) = radius(comp.nrfetch, r)
density(comp::VariableNR, r::Point, _) = density(comp.nrfetch, r)
radius(comp::VariableNR, r::Point, rad) = rad
mie_qext(comp::VariableNR, r::Point, rad) = 2. + comp.c * power34(rad)
mie_g(comp::VariableNR, r::Point, rad) = comp.g0 * rad / (rad + comp.r0)
mie_ω0(comp::VariableNR, r::Point, rad) = 1. - comp.a * rad
νmiemax(comp::VariableNR) = comp.νmiemax

# Compute a linear interpolator for the refraction index from a datafile
function m_interpolator(fname)
    # If the provided filename does not exist we look for it in our
    # data directory.
    isfile(fname) || (fname = joinpath(@__DIR__, "..", "data", fname))

    λ1, n1, k1 = eachcol(readdlm(fname, comments=true))

    # Lambda is in micrometers in the input file.
    λ1 *= co.micro

    interp_n1 = LinearInterpolation(λ1, n1)
    interp_k1 = LinearInterpolation(λ1, k1)

    function interp(λ)
        # Here we ensure type stability of the return value
        a::Float64 = interp_n1(λ)
        b::Float64 = interp_k1(λ)
        complex(a, b)
    end
end


text:
"""Checks output folders to determine a 4-digit run id number."""
function get_run_id_path(P::Params)

    @unpack output,outpath = P

    if output
        runlist = filter(x->startswith(x,"run"),readdir(outpath))
        existing_runs = [parse(Int,id[4:end]) for id in runlist]
        run_id = maximum(existing_runs)+1
        runpath = joinpath(outpath,@sprintf("run%04d",run_id))
        mkdir(runpath)

        return run_id,runpath
    else
        return 0,"no runpath"
    end
end

function output_initialise( Diag::DiagnosticVariables,   # output grid-space variables only
                            G::GeoSpectral,
                            P::Params)

    @unpack nlon,nlat,nlev = G.geometry
    @unpack lon,lat = G.geometry
    @unpack output_startdate, compression_level = P

    # INITIALISE FOLDER
    run_id,runpath = get_run_id_path(P)
    file_name = @sprintf("run%04d",run_id)

    # TIME
    time_string = "hours since $(Dates.format(output_startdate, "yyyy-mm-dd HH:MM:0.0"))"
    timedim = NcDim("time",0,unlimited=true)
    timevar = NcVar("time",timedim,t=Int32,atts=Dict("units"=>time_string))

    # SPACE
    londim  = NcDim("lon",nlon,values=lon)
    latdim  = NcDim("lat",nlat,values=lat)
    levdim  = NcDim("lev",nlev,values=collect(1:nlev))
    lonvar  = NcVar("lon",londim,t=Float32,atts=Dict("long_name"=>"longitude"))
    latvar  = NcVar("lat",latdim,t=Float32,atts=Dict("long_name"=>"latitude"))
    levvar  = NcVar("lev",levdim,t=Int32,atts=Dict("long_name"=>"model level"))

    # VARIABLES
    uvar        = NcVar("u",[londim,latdim,levdim,timedim],t=Float32,compress=compression_level,
                    atts=Dict("long_name"=>"zonal wind","units"=>"m/s","missing_value"=>-999999f0))
    vvar        = NcVar("v",[londim,latdim,levdim,timedim],t=Float32,compress=compression_level,
                    atts=Dict("long_name"=>"meridional wind","units"=>"m/s","missing_value"=>-999999f0))
    Tvar        = NcVar("T",[londim,latdim,levdim,timedim],t=Float32,compress=compression_level,
                    atts=Dict("long_name"=>"temperature","units"=>"K","missing_value"=>-999999f0))
    humidvar    = NcVar("humid",[londim,latdim,levdim,timedim],t=Float32,compress=compression_level,
                    atts=Dict("long_name"=>"specific humidity","units"=>"1","missing_value"=>-999999f0))
    geopotvar   = NcVar("geopot",[londim,latdim,levdim,timedim],t=Float32,compress=compression_level,
                    atts=Dict("long_name"=>"geopotential height","units"=>"m","missing_value"=>-999999f0))
    logp0var    = NcVar("logp0",[londim,latdim,timedim],t=Float32,compress=compression_level,
                    atts=Dict("long_name"=>"surface pressure","units"=>"Pa","missing_value"=>-999999f0))

    # CREATE NETCDF FILE
    netcdffile = NetCDF.create("$file_name",
                    [timevar,lonvar,latvar,levvar,uvar,vvar,Tvar,humidvar,geopotvar,logp0var],mode=NetCDF.NC_NETCDF4)

    # WRITE INITIAL CONDITIONS TO FILE
    write_output!(0,0,netcdffile,Diag,P)

    return netcdffile
end

function write_output!( iout::Int,
                        t::Real,
                        ncfile::NcFile,
                        Diag::DiagnosticVariables,
                        P::Params)

    @unpack u,v,Tabs,humid,logp0 = Diag.gridvars
    @unpack keepbits = P

    # always convert to float32 for output
    u_output = Float32.(u)
    v_output = Float32.(v)
    T_output = Float32.(Tabs)
    humid_output = Float32.(humid)
    logp0 = Float32.(logp0)

    # ROUNDING FOR ROUND+LOSSLESS COMPRESSION
    for var in (u_output,v_output,T_output,humid_output)
        round!(var,keepbits)
    end

    # WRITE TO FILE
    NetCDF.putvar(ncfile,"u",u_output,start=[1,1,1,iout],count=[-1,-1,-1,1])
    NetCDF.putvar(ncfile,"v",v_output,start=[1,1,1,iout],count=[-1,-1,-1,1])
    NetCDF.putvar(ncfile,"T",T_output,start=[1,1,1,iout],count=[-1,-1,-1,1])
    NetCDF.putvar(ncfile,"humid",humid_output,start=[1,1,1,iout],count=[-1,-1,-1,1])
    NetCDF.putvar(ncfile,"geopot",geopot_output,start=[1,1,1,iout],count=[-1,-1,-1,1])
    NetCDF.putvar(ncfile,"logp0",logp0_output,start=[1,1,iout],count=[-1,-1,1])

    # WRITE TIME
    NetCDF.putvar(ncfile,"time",Int32[t],start=[iout])
end


text:
@testset "Testing Bipartite Matching" begin
    using LightGraphs

    ############################################################################
    # Auxiliary functions to help with testing.
    ############################################################################
    function build_testgraph(d)
        # create empty array to count the number of elements on RHS
        value_array = []
        for value in values(d)
            for entry in value
                !(in(entry, value_array)) && push!(value_array, entry)
            end
        end
        # number of nodes equals to entries on RHS + entries on LHS + source node
        # + sink node
        g = DiGraph(length(d) + length(value_array) + 2)
        # add edges
        for (key, values) in d
            add_edge!(g, 1 => key)
            for value in values
                add_edge!(g, (key) => (value))
                add_edge!(g, (value) => 2)
            end
        end
        return g
    end

    # this function converts the matched graph into a dictionary format used in
    # verification
    function graph2dict(g)
        dict = Dict{Int64,Int64}()
        for a in outneighbors(g, 1)
            for b in inneighbors(g, a)
                if (b != 1)
                    dict[a] = b
                end
            end
        end
        return dict
    end

    ############################################################################
    # Perform Test
    ############################################################################

    dict = Dict(
        3 => [9], 4 => [9, 10, 11], 5 => [11, 12], 6 => [12], 7 => [12, 13], 8 => [9, 13]
    )
    # Convert dict to LightGraph
    graph = build_testgraph(dict)
    # Do matching on the graph
    matched_graph = Mapper2.SA.bipartite_match!(graph)
    # Convert matched graph back to dictonary
    d = graph2dict(matched_graph)
    # Generate expected result
    d_ref = Dict(3 => 9, 4 => 10, 5 => 11, 6 => 12, 7 => 13)

    @test length(keys(d_ref)) == length(keys(d))
end


text:
"""
Bayesian Gaussian Mixture Model
"""
module GaussianMixtureModel
using StatsFuns.logsumexp, SpecialFunctions.digamma
using Distributions
using PDMats

export GW, BGMM, Gauss, GMM
export sample_GMM, sample_data, winner_takes_all
export learn_GS, learn_CGS, learn_VI

####################
## Types
struct GW
    # Parameters of Gauss Wisahrt distribution
    beta::Float64
    m::Vector{Float64}
    nu::Float64
    W::Matrix{Float64}
end

struct BGMM
    # Parameters of Bayesian Gaussian Mixture Model
    D::Int
    K::Int
    alpha::Vector{Float64}
    cmp::Vector{GW}
end

struct Gauss
    # Parameters of Gauss Distribution
    mu::Vector{Float64}
    Lambda::Matrix{Float64}
end

struct GMM
    # Parameters of Gauss Mixture Model
    D::Int
    K::Int
    phi::Vector{Float64}
    cmp::Vector{Gauss}
end

####################
## Common functions
"""
Sample a GMM given hyperparameters.
"""
function sample_GMM(bgmm::BGMM)
    cmp = Vector{Gauss}()
    for c in bgmm.cmp
        Lambda = rand(Wishart(c.nu, PDMats.PDMat(Symmetric(c.W))))
        mu = rand(MvNormal(c.m, PDMats.PDMat(Symmetric(inv(c.beta*Lambda)))))
        push!(cmp, Gauss(mu, Lambda))
    end
    phi = rand(Dirichlet(bgmm.alpha))
    return GMM(bgmm.D, bgmm.K, phi, cmp)
end

"""
Sample data from a specific GMM model.
"""
function sample_data(gmm::GMM, N::Int)
    X = zeros(gmm.D, N)
    S = categorical_sample(gmm.phi, N)
    for n in 1 : N
        k = indmax(S[:, n])
        X[:,n] = rand(MvNormal(gmm.cmp[k].mu, PDMats.PDMat(Symmetric(inv(gmm.cmp[k].Lambda)))))
    end
    return X, S
end

categorical_sample(p::Vector{Float64}) = categorical_sample(p, 1)[:,1]
function categorical_sample(p::Vector{Float64}, N::Int)
    K = length(p)
    S = zeros(K, N)
    S_tmp = rand(Categorical(p), N)
    for k in 1 : K
        S[k,find(S_tmp.==k)] = 1
    end
    return S
end

function sumdigamma(nu, D)
    ret = 0.0
    for d in 1 : D
        ret += digamma.(0.5*(nu + 1 - d))
    end
    return ret
end

function init_S(X::Matrix{Float64}, bgmm::BGMM)
    N = size(X, 2)
    K = bgmm.K
    S = categorical_sample(ones(K)/K, N)
    return S
end

function calc_ELBO(X::Array{Float64, 2}, pri::BGMM, pos::BGMM)
    function logCw(nu, W)
        D = size(W, 1)
        return -0.5*nu*logdet(W) - 0.5*nu*D*log.(2) - 0.25*D*(D-1)*log.(pi) - sum([lgamma.(0.5*(nu+1-d)) for d in 1 : D])
    end

    ln_expt_S = update_S(pos, X)
    expt_S = exp.(ln_expt_S)
    K, N = size(expt_S)
    D = size(X, 1)

    expt_ln_lkh = 0
    for k in 1 : K
        expt_Lambda = pos.cmp[k].nu * pos.cmp[k].W
        expt_Lambda_mu = pos.cmp[k].nu * pos.cmp[k].W * pos.cmp[k].m
        expt_mu_Lambda_mu = (pos.cmp[k].nu * pos.cmp[k].m' * pos.cmp[k].W * pos.cmp[k].m)[1] + D/pos.cmp[k].beta
        expt_ln_Lambda = sumdigamma(pos.cmp[k].nu, D) + D*log.(2) + logdet(pos.cmp[k].W)
        expt_ln_pi = digamma.(pos.alpha) - digamma.(sum(pos.alpha))
        for n in 1 : N
            # <ln p(X|S, mu, Lambda)>
            expt_ln_lkh += -0.5 * expt_S[k,n]*(trace(X[:,n]*X[:,n]'*expt_Lambda)
                                               - 2*(X[:,n]'*expt_Lambda_mu)[1]
                                               + expt_mu_Lambda_mu
                                               - expt_ln_Lambda
                                               + D * log.(2*pi)
                                               )
            # <ln p(S|pi)>
            expt_ln_lkh += expt_S[k,n]*expt_ln_pi[k]
        end
    end
    # -<ln q(S)>
    expt_ln_lkh -= sum(expt_S.*ln_expt_S)

    KL_mu_Lambda = [(0.5*D*(log.(pos.cmp[k].beta) - log.(pri.cmp[k].beta) + pri.cmp[k].beta/pos.cmp[k].beta - pos.cmp[k].nu - 1)
                     + 0.5*(pos.cmp[k].nu-pri.cmp[k].nu)*(sumdigamma(pos.cmp[k].nu, D) + D*log.(2) + logdet(pos.cmp[k].W))
                     + logCw(pos.cmp[k].nu, pos.cmp[k].W) - logCw(pri.cmp[k].nu, pri.cmp[k].W)
                     + 0.5*pos.cmp[k].nu*trace((pri.cmp[k].beta*(pos.cmp[k].m-pri.cmp[k].m)*(pos.cmp[k].m-pri.cmp[k].m)'
                                                +inv(pri.cmp[k].W))*pos.cmp[k].W)) for k in 1 : K]

    KL_pi = (lgamma.(sum(pos.alpha)) - lgamma.(sum(pri.alpha))
             - sum(lgamma.(pos.alpha)) + sum(lgamma.(pri.alpha))
             + (pos.alpha - pri.alpha)' * (digamma.(pos.alpha) - digamma.(sum(pos.alpha)))
             )[1]

    VB = expt_ln_lkh - (sum(KL_mu_Lambda) + KL_pi)
    return VB
end

function add_stats(bgmm::BGMM, X::Matrix{Float64}, S::Matrix{Float64})
    D = bgmm.D
    K = bgmm.K
    sum_S = sum(S, 2)
    alpha = [bgmm.alpha[k] + sum_S[k] for k in 1 : K]
    cmp = Vector{GW}()

    XS = X*S';
    for k in 1 : K
        beta = bgmm.cmp[k].beta + sum_S[k]
        m = (1.0/beta)*(vec(X*S[[k],:]') + bgmm.cmp[k].beta*bgmm.cmp[k].m)
        nu = bgmm.cmp[k].nu + sum_S[k]
        W = inv(X*diagm(S[k,:])*X'
                       - beta*m*m'
                       + bgmm.cmp[k].beta*bgmm.cmp[k].m*bgmm.cmp[k].m'
                       + inv(bgmm.cmp[k].W))

        push!(cmp, GW(beta, m, nu, W))
    end
    return BGMM(D, K, alpha, cmp)
end

remove_stats(bgmm::BGMM, X::Matrix{Float64}, S::Matrix{Float64}) = add_stats(bgmm, X, -S)

####################
## used for Variational Inference
function update_S(bgmm::BGMM, X::Matrix{Float64})
    D, N = size(X)
    K = bgmm.K
    ln_S = zeros(K, N)
    tmp = zeros(K)

    tmp = NaN * zeros(K)
    sum_digamma_tmp = digamma.(sum(bgmm.alpha))
    for k in 1 : K
        tmp[k] = -0.5*(bgmm.cmp[k].nu*trace(bgmm.cmp[k].m*bgmm.cmp[k].m'*bgmm.cmp[k].W)
                       + D*(1.0/bgmm.cmp[k].beta)
                       - (sumdigamma(bgmm.cmp[k].nu, D) + logdet(bgmm.cmp[k].W)))
        tmp[k] += digamma.(bgmm.alpha[k]) - sum_digamma_tmp
    end
    for n in 1 : N
        tmp_ln_pi = NaN * zeros(K)
        for k in 1 : K
            tmp_ln_pi[k] = tmp[k] -0.5*bgmm.cmp[k].nu*trace((X[:,n]*X[:,n]' - 2*bgmm.cmp[k].m*X[:,n]')*bgmm.cmp[k].W)
        end
        ln_S[:,n] = tmp_ln_pi - logsumexp(tmp_ln_pi)
    end
    return ln_S
end

"""
Pick single states having a max probability.
"""
function winner_takes_all(S::Matrix{Float64})
    S_ret = zeros(size(S))
    for n in 1 : size(S_ret, 2)
        idx = indmax(S[:,n])
        S_ret[idx,n] = 1
    end
    return S_ret
end

####################
## used for Gibbs Sampling
function sample_S_GS(gmm::GMM, X::Matrix{Float64})
    D, N = size(X)
    K = gmm.K
    S = zeros(K, N)

    tmp = [0.5*logdet(gmm.cmp[k].Lambda) + log.(gmm.phi[k]) for k in 1 : K]

    for n in 1 : N
        tmp_ln_phi = [-0.5*trace(gmm.cmp[k].Lambda*(X[:,n] - gmm.cmp[k].mu)*(X[:,n] - gmm.cmp[k].mu)') + tmp[k] for k in 1 : K]
        tmp_ln_phi = tmp_ln_phi - logsumexp(tmp_ln_phi)
        S[:,n] = categorical_sample(exp.(tmp_ln_phi))
    end

    return S
end

####################
## used for Collapsed Gibbs Sampling
function calc_ln_ST(Xn::Vector{Float64}, gw::GW)
    # TODO; need to check value?
    D = size(Xn, 1)
    W = ((1 - D + gw.nu)*gw.beta / (1 + gw.beta)) * gw.W
    #ln_lkh = logpdf(MvTDist(1 - D + gw.nu, gw.m, (gw.nu/(gw.nu - 2))*inv(W)), Xn)
    ln_lkh = logpdf(MvTDist(1 - D + gw.nu, gw.m, PDMats.PDMat(Symmetric(inv(W)))), Xn)
    return sum(ln_lkh)
end

function sample_Sn(Xn::Vector{Float64}, bgmm::BGMM)
    ln_tmp = [(calc_ln_ST(Xn, bgmm.cmp[k]) + log.(bgmm.alpha[k])) for k in 1 : bgmm.K]
    ln_tmp = ln_tmp -  logsumexp(ln_tmp)
    Sn = categorical_sample(exp.(ln_tmp))
    return Sn
end

function sample_S_CGS(S::Matrix{Float64}, X::Matrix{Float64}, bgmm::BGMM)
    D, N = size(X)
    K = size(S, 1)
    for n in randperm(N)
        # remove
        bgmm = remove_stats(bgmm, X[:,[n]], S[:,[n]])
        # sample
        S[:,n] = sample_Sn(X[:,n], bgmm)
        # insert
        bgmm = add_stats(bgmm, X[:,[n]], S[:,[n]])
    end
    return S, bgmm
end

####################
## Algorithm main
"""
Compute posterior distributions via variational inference.
"""
function learn_VI(X::Matrix{Float64}, prior_bgmm::BGMM, max_iter::Int)
    # initialisation
    expt_S = init_S(X, prior_bgmm)
    bgmm = add_stats(prior_bgmm, X, expt_S)
    VB = NaN * zeros(max_iter)

    # inference
    for i in 1 : max_iter
        # E-step
        expt_S = exp.(update_S(bgmm, X))
        # M-step
        bgmm = add_stats(prior_bgmm, X, expt_S)
        # calc VB
        VB[i] = calc_ELBO(X, prior_bgmm, bgmm)
    end

    # assign binary values
    S = winner_takes_all(expt_S)
    return S, bgmm, VB
end

"""
Compute posterior distributions via Gibbs sampling.
"""
function learn_GS(X::Matrix{Float64}, prior_bgmm::BGMM, max_iter::Int)
    # initialisation
    S = init_S(X, prior_bgmm)
    bgmm = add_stats(prior_bgmm, X, S)
    VB = NaN * zeros(max_iter)

    # inference
    for i in 1 : max_iter
        # sample parameters
        gmm = sample_GMM(bgmm)
        # sample latent variables
        S = sample_S_GS(gmm, X)
        # update current model
        bgmm = add_stats(prior_bgmm, X, S)
        # calc VB
        VB[i] = calc_ELBO(X, prior_bgmm, bgmm)
    end

    return S, bgmm, VB
end

"""
Compute posterior distributions via collapsed Gibbs sampling.
"""
function learn_CGS(X::Matrix{Float64}, prior_bgmm::BGMM, max_iter::Int)
    # initialisation
    S = init_S(X, prior_bgmm)
    bgmm = add_stats(prior_bgmm, X, S)
    VB = NaN * zeros(max_iter)

    # inference
    for i in 1 : max_iter
        # directly sample S
        S, bgmm = sample_S_CGS(S, X, bgmm)
        # calc VB
        VB[i] = calc_ELBO(X, prior_bgmm, bgmm)
    end

    return S, bgmm, VB
end

end


text:
module PanX

using Printf
using JSON, GZip
using PanGraph

struct Gene
    id     :: Int
    tag    :: String
    low    :: Int
    high   :: Int
    strand :: Bool
end

function find(database, tag)
    for (i,cluster) in enumerate(database)
        tag ∈ cluster && return i
    end
    return nothing
end

askey(name,tag) = "$(name)|$(tag)"

function gene(line, name, database)
    elt = split(line)

    strand = elt[1] == "+"
    low, high = parse.(Int, split(elt[2],".."))
    tag = elt[3]
    id  = find(database,askey(name,tag))
    id !== nothing || return nothing

    return Gene(id, tag, low, high, strand)
end

struct Genome
    name     :: String
    sequence :: Array{UInt8,1}
    gene     :: Array{Gene,1}
end

Genome(name::String, sequence::Array{UInt8,1}) = Genome(name, sequence, Gene[])

function database(path)
    data = Set{String}[]
    meta = String[]
    re   = r"(N[A-Z]_[A-Z0-9]+)_"
    sub  = s"\1|"
    open(path) do io
        genes = JSON.parse(io)
        for (i,gene) in enumerate(genes)
            loci = map(split(gene["locus"])) do g
                replace(g, re => sub)
            end
            push!(data, Set(loci))
            push!(meta, gene["msa"])
        end
    end
    return (
        loci=data,
        meta=meta,
        path=dirname(path),
    )
end

function alignment(io::IO)
    record = [ (name=rec.name, seq=rec.seq) for rec in PanGraph.Utility.read_fasta(io)]
    return hcat(map(r->r.seq,record)...), map(record) do r
        entry = split(r.name, "-")
        return "$(entry[1])|$(entry[2])"
    end
end

function pangraph(genomes, genes)
    vertex = [
        let
            path = "$(genes.path)/geneCluster/$(base)_na_aln.fa.gz"
            aln, isos = GZip.open(alignment, path)

            ref = PanGraph.Utility.make_consensus(aln)
            blk = PanGraph.Graphs.Block(ref)
            iso = [ PanGraph.Graphs.Node(blk, true) for _ in isos ]

            blk.uuid = @sprintf("Gene_%05d", i)
            blk.gaps,
            blk.mutate,
            blk.delete,
            blk.insert,
            blk.sequence = PanGraph.Utility.alignment_alleles(ref, aln, iso)

            (block=blk, node=Dict( name => node for (name,node) in zip(isos,iso)))
        end for (i,(gene,base)) in enumerate(zip(genes.loci,genes.meta))
    ]

    path = [
        let
            nodes = [
             let
                 node = pop!(vertex[gene.id].node, askey(genome.name, gene.tag))
                 node.strand = gene.strand
                 node
             end for gene in genome.gene
            ]
            breaks = [ x for gene in genome.gene for x in [gene.low, gene.high] ]
            PanGraph.Graphs.Path(genome.name, nodes, 0, true, breaks)
        end for genome in values(genomes)
    ]

    # remove any unpaired nodes we initially added
    # delete empty blocks
    index = Int[]
    for (i,v) in enumerate(vertex)
        for n in values(v.node)
            pop!(v.block, n)
        end

        if length(v.block.mutate) == 0
            push!(index, i)
        end
    end
    deleteat!(vertex, index)

    block = map(v->v.block, vertex)
    graph = PanGraph.Graph(Dict(b.uuid => b for b in block), Dict(p.name => p for p in path) )

    return graph
end

function main(database)
    genomes = Dict{String,Genome}()
    while !eof(stdin)
        line = readline(stdin)
        line[1] == '>' || error("invalid stream detected")

        name = line[2:end]
        sbuf = IOBuffer()

        # grab sequence
        @label moreseq
        line = readline(stdin)
        if line != "ENDSEQ"
            print(sbuf, line)
            @goto moreseq
        end
        genome = Genome(name, take!(sbuf))

        # grab loci
        @label moregene
        line = readline(stdin)
        if line != "ENDGENE"
            g = gene(line, name, database.loci)
            g !== nothing && push!(genome.gene, g)
            @goto moregene
        end

        genomes[name] = genome
    end

    return pangraph(genomes, database)
end

# NOTE: just for debugging
function echo()
    while !eof(stdin)
        println(readline(stdin))
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    database(ARGS[1]) |> main |> (G) -> PanGraph.Graphs.marshal(stdout, G; fmt=:json)
end

end

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
# binary search in interval (low,high)
# assumes that if pred(x)=true, then pred(y)=true ∀y>x
# assumes pred(low)=false and pred(high)=true
function binarysearchinterval(pred, low::T, high::T, precision::T) where {T<:AbstractFloat}
    while high-low > precision
        mid = (low+high)/2

        if pred(mid)
            high = mid
        else
            low = mid
        end
    end
    (low+high)/2
end
binarysearchinterval(pred, low, high, precision) = binarysearchinterval(pred, promote(low, high, precision)...)

# binary search predicate function ||log₂(f+α)-log₂(r+α)||_(RMS) < threshold
function lodpred(α, f, r, threshold)
    N = length(f)
    norm(log2.(f.+α)-log2.(r.+α)) <= threshold*sqrt(N)
end
function _limitofdetection(freqsF::Vector{T}, freqsR::Vector{T}, threshold, lgLow, lgHigh) where {T}
    @assert size(freqsF)==size(freqsR)
    lodpred(10.0^lgLow,freqsF,freqsR,log2(threshold)) && return 10.0^lgLow # early out
    10.0^binarysearchinterval( α->lodpred(10.0^α,freqsF,freqsR,log2(threshold)), lgLow, lgHigh, 1e-3) # logspace search for α
end

# Two 4x4x4xN arrays, one per strand. N=Number of samples.
function _limitofdetection(freqs::Array{T,4}, freqsF::Array{T,4}, freqsR::Array{T,4}, threshold, lgLow, lgHigh) where {T}
    N = size(freqs,4)
    @assert size(freqs)==size(freqsF)==size(freqsR)==(4,4,4,N)

    # group samples by consensus codon
    consensusInd = map(i->argmax(freqs[:,:,:,i][:]), 1:N)

    lod = (10.0^lgLow)*ones(T,4,4,4)

    # compute per group and per variant
    for c in unique(consensusInd)
        mask = consensusInd.==c
        fF = freqsF[:,:,:,mask]
        fR = freqsR[:,:,:,mask]
        for k=1:4,j=1:4,i=1:4
            x = _limitofdetection(fF[i,j,k,:][:], fR[i,j,k,:][:], threshold, lgLow, lgHigh)
            lod[i,j,k] = max(x, lod[i,j,k])
        end
    end
    lod
end

# All swarms should be from the same run.
function _limitofdetection(swarms::AnnotatedArray, swarmsF::AnnotatedArray, swarmsR::AnnotatedArray, threshold, lgLow, lgHigh)
    @assert size(swarms)==size(swarmsF)==size(swarmsR)
    @assert swarms[:id]==swarmsF[:id]==swarmsR[:id]
    P = size(swarms,4)
    N = size(swarms,5)

    lod = zeros(4,4,4,P)
    # TODO: parallelize here???
    for p=1:P
        lod[:,:,:,p] = _limitofdetection(convert(Array,swarms)[:,:,:,p,:],
                                         convert(Array,swarmsF)[:,:,:,p,:],
                                         convert(Array,swarmsR)[:,:,:,p,:],
                                         threshold, lgLow, lgHigh)
    end
    lod
end

function limitofdetection(swarms::AnnotatedArray, swarmsF::AnnotatedArray, swarmsR::AnnotatedArray; groupID=ones(Int,size(swarms,5)), threshold=1.5, lgLow=-6, lgHigh=1)
    @assert length(groupID)==size(swarms,5)

    P = size(swarms,4)
    uniqueGroups = unique(groupID)

    lods = zeros(4,4,4,P,length(uniqueGroups))
    # TODO: parallelize here???
    for (i,g) in enumerate(uniqueGroups)
        mask = groupID.==g
        lods[:,:,:,:,i] = _limitofdetection(swarms[:,:,:,:,mask], swarmsF[:,:,:,:,mask], swarmsR[:,:,:,:,mask], threshold, lgLow, lgHigh)
    end
    lod = dropdims(maximum(lods; dims=5); dims=5)
    lod = AnnotatedArray(lod)
    annotate!(lod, :segment,   dropdims(swarms[:segment]; dims=5))
    annotate!(lod, :position,  dropdims(swarms[:position]; dims=5))
    annotate!(lod, :n1,        dropdims(swarms[:n1]; dims=5))
    annotate!(lod, :n2,        dropdims(swarms[:n2]; dims=5))
    annotate!(lod, :n3,        dropdims(swarms[:n3]; dims=5))
    annotate!(lod, :aminoacid, dropdims(swarms[:aminoacid]; dims=5))
    annotate!(lod, :codon,     dropdims(swarms[:codon]; dims=5))
    lod
end


text:
"""
Compute the principal eigenvector and eigenvalue of a linear operator 𝕋, where 𝕋 is a Metzler matrix (i.e. off-diagonal components are nonnegative)

Denote a = -minimum(Diagonal(V)), which implies 𝕋 + a * I has all positive entries. Applying Perron Frobenus, there a unique largest eigenvalue for aI + 𝕋, which is real, and the correspondind eigenctor is strictly positive.
Note that, in particular, it is the eigenvalue with largest real part, and so this also correspoinds to the eigenvalue with largest real part of 𝕋, which happens to be real.
Denote η(𝕋) the eigenvalue with largest real part of a matrix and ρ(𝕋) the eigenvalue with largest modulus. We have

    η(𝕋) = ρ(𝕋 + a * I) - a

Moreover, the associated eigenvector is real and strictly positive.

Note that, when 𝕋 is generator, its rows sum to zero. This means that eigenvalue with largest real part is
    η(𝕋) = 0
In other words, all eigenvalues of 𝕋 have real part <= 0. This means that 𝕋 is a singular M matrix.
(another proof is to say that, for any s, sI - 𝕋 is a non-singular M-Matrix for any s> 0, since there exists x = e such that (sI - 𝕋) * x > 0).

"""
function principal_eigenvalue(𝕋; r0 = ones(size(𝕋, 1)))
    η, r = 0.0, r0
    a = - minimum(diag(𝕋))
    try
        vals, vecs = Arpack.eigs(𝕋 + a * I; v0 = collect(r0), nev = 1, which = :LM)
        η = vals[1]
        r = vecs[:, 1]
    catch
        vals, vecs = KrylovKit.eigsolve(𝕋 + a * I, collect(r0), 1, :LM; maxiter = size(𝕋, 1))
        η = vals[1]
        r = vecs[1]
    end
    abs(imag(η)) <= eps() || @warn "Principal Eigenvalue has an imaginary part"
    maximum(abs.(imag.(r))) <= eps() || @warn "Principal Eigenvector has an imaginary part"
    real(η) - a, abs.(r)
end


text:
using GPFlux
using LinearAlgebra
using GaussianProcesses
using Statistics
using Calculus
using Zygote
using Test, Random

function K_check(mykernel, gpkernel, X)
  myK = mykernel(X)
	gpK = cov(gpkernel, X) + Diagonal(1e-6*ones(size(X, 2)))
	@show isposdef(myK)
	@show isposdef(gpK)
	isposdef(myK) == isposdef(gpK) &&  isapprox(myK, gpK, rtol=1e-3, atol=1e-8)
end

function dKdθ_check(mykernel, gpkernel, X, v)
  N = size(X, 2)

  myps = GPFlux.params(mykernel)
  f = () -> v'*mykernel(X)*v
  mygs = gradient(f, myps)
  mygv = [mygs.grads[p][1] for p in myps]

  n_gpps = GaussianProcesses.num_params(gpkernel)
  gpps = GaussianProcesses.get_params(gpkernel)
  gpdKdθ = zeros(N, N, n_gpps)
  gpmetric = GaussianProcesses.KernelData(gpkernel, X, X)
  GaussianProcesses.grad_stack!(gpdKdθ, gpkernel, X, X, gpmetric)
  gpgv = [tr(v*v'*gpdKdθ[:, :, i]) for i in 1:n_gpps]

	# gpgv = [tr(v*v'*gpdKdθ[:, :, i]).*exp(-gpps[i]) for i in 1:n_gpps]

#  calps = [p[1] for p in myps]
#  calgv = Calculus.gradient(calps) do ps
#      ps_ = [[p] for p in ps]
#      kernel = typeof(mykernel)(ps_...)
#      v'*kernel(X)*v
#  end

	@show sort(mygv)
# 	@show sort(calgv)
	@show sort(gpgv)
#  isapprox(sort(mygv),sort(calgv),rtol=1e-3,atol=1e-8) &&
	isapprox(sort(mygv),sort(gpgv),rtol=1e-3,atol=1e-8)
end

@testset "kernel check" begin
	rng = MersenneTwister(4)
	d = 4
	n = 10
	X = rand(rng, d, n)
	v = rand(rng, n)
	kernel_pairs = [(IsoGaussKernel([4.0], [1.0]), SE(4.0, 1.0)),
									(IsoPeriodKernel([0.0], [0.0], [1.0]), Periodic(0.0, 1.0, 0.0)),
									(IsoRQKernel([0.0], [-1.0], [0.0]), RQ(0.0, 0.0, -1.0)),
									(IsoLinearKernel([0.0]), Lin(0.0)),
								  (SE_add_PeriodKernel(IsoGaussKernel([4.0],[0.0]),IsoPeriodKernel([0.0],[0.0],[1.0])), SE(4.0,0.0)+Periodic(0.0,1.0,0.0)),
									(SE_mul_PeriodKernel(IsoGaussKernel([4.0],[0.0]),IsoPeriodKernel([0.0],[0.0],[1.0])), SE(4.0,0.0)*Periodic(0.0,1.0,0.0)),
									]

	for (mykernel, gpkernel) in kernel_pairs
		@show (mykernel, gpkernel)
		@test K_check(mykernel, gpkernel, X)
		@test dKdθ_check(mykernel, gpkernel, X, v)
	end

end

max_issues_repo_name: HamletWantToCode/GPFlux

max_forks_repo_name: HamletWantToCode/GPFlux


text: const enable_backprop = Ref(true)


text:
################################################################################
#
#   module LatPhysLuttingerTisza
#   -> LatticePhysics_Base
#   -> LatticePhysics_ReciprocalSpace
#   -> LinearAlgebra
#
#   --> TYPE DEFINITIONS
#           - LTBandstructure
#           - LTGroundstateManifold
#
#   --> IO OF TYPES
#
#   --> SPIN INTERACTION MATRICES
#
#   --> CALCULATION OF LT CONSTRAINTS
#
#   --> CALCULATION OF LT BANDSTRUCTURES AND EIGENVALUES
#           - in real space (pure eigenvalues)
#           - in momentum space (along cardinal directions)
#           - in momentum space (along reciprocal lattice directions)
#           - in momentum space (along path)
#
#   --> CALCULATION OF LT GROUNDSTATE MANIFOLDS
#           - calculation for 2D unitcells
#           - calculation for 3D unitcells
#
################################################################################

# start of module
module LatPhysLuttingerTisza

# used libraries
using LatPhysBase
using LatPhysReciprocal
using LinearAlgebra

# end of module
end


text:
# meanconfid
function meanconfid_test()
	x = [2.,3.,5.,6.,9.]
	x_mean,x_confid = meanconfid(x,sigma=2.5,alpha=.05)
	@test x_mean ≈ 5.0
	@test x_confid[1] ≈ 2.81 atol=0.02
	@test x_confid[2] ≈ 7.19 atol=0.02

	# Unknown sigma
	x_mean,x_confid = meanconfid(x)
	@test x_mean ≈ 5.0
	@test x_confid[1] ≈ 1.60 atol=0.02
	@test x_confid[2] ≈ 8.4 atol=0.02

	# with NaN
	push!(x,NaN);
	x_mean,x_confid = meanconfid(x)
	@test x_mean ≈ 5.0
	@test x_confid[1] ≈ 1.60 atol=0.02
	@test x_confid[2] ≈ 8.4 atol=0.02
end

# meanconfid
function stdconfid_test()
	x = [1.95,1.96,1.98,1.99,2.00,2.00,2.01,2.01,2.01,2.01,
		2.02,2.02,2.03,2.03,2.04,2.04,2.05,2.06,2.08,2.10];
	# unknown mean value
	s,s_interval = stdconfid(x,alpha=0.05)
	@test s^2 ≈ 0.001352 atol = 0.00001
	@test s_interval[1] ≈ 0.027964 atol = 0.00001
	@test s_interval[2] ≈ 0.053710 atol = 0.00001

	# known mean
	x = [19.01,19.02,18.99,19.00,19.05,19.05,19.0,18.98,18.99,19.0];
	s,s_interval = stdconfid(x,alpha=0.05,x_mean=19.)
	@test s^2 ≈ 0.00061 atol = 0.000001
end

# Run
meanconfid_test()
stdconfid_test()

max_stars_repo_name: emenems/HypoTest

max_issues_repo_name: emenems/HypoTest

max_forks_repo_name: emenems/HypoTest