text:
"""
    HarmonicPowerModels.replicate(data::Dict{String,<:Any};
                                  harmonics::Array{Int}=Int[],
                                  xfmr_exc::Bool=false)

Creates a multinetwork data-file for all considered harmonics based on the data
specified for the fundamental harmonic.

The considered harmonics are provided by the user through the named
argument harmonics, and collected based on the harmonic bus data.

Additionally, a Bool `xfmr_exc` may be passed to include the excitation current
of the transformers.

NOTES:
i) branch
-   The susceptance b of any branch is assumed to be purely capacitive, i.e.,
    b = ω * C.
-   The resistance r of any branch is scaled using the square of the harmonic.
ii) transformer
-   The transformer configuration defaults to Yy0.
-   The resistance r of any transformer is scaled using the square of the harmonic.
"""
function _HPM.replicate(data::Dict{String, Any};
                        harmonics::Array{Int}=Int[],
                        xfmr_exc::Dict{String, Any}=Dict{String, Any}())
    # extend the user-provided harmonics based on the data
    collect_harmonics!(data, harmonics, xfmr_exc)

    # for (t, xfmr) in data["xfmr"]
    #     xfmr["voltage_harmonics"] = []
    #     xfmr["current_harmonics"] = []
    #     xfmr["voltage_harmonics_ntws"] = []
    #     xfmr["current_harmonics_ntws"] = []
    # end
    @show harmonics
    # create a multinetwork data structure
    Nh = length(harmonics)
    data = _PMs.replicate(data, Nh)

    # add the harmonics to the overall data
    data["harmonics"] = Dict("$nw" => harmonics[nw] for nw in 1:Nh)

    # add the xfmr excitation
    sample_xfmr_excitation(data, xfmr_exc)

    # re-evaluate the data for each harmonic
    for nw in keys(data["nw"])
        nh = data["harmonics"][nw]
        mh = 2 * pi * nh

        # re-evaluate load power
        for load in values(data["nw"][nw]["load"])
            nb = load["source_id"][2]
            bus = data["nw"][nw]["bus"]["$nb"]
            multiplier = haskey(bus, "nh_$nh") ? bus["nh_$nh"] : 0.0 ;
            haskey(load, "pd") ? load["pd"] *= multiplier : ~ ;
            haskey(load, "qd") ? load["qd"] *= multiplier : ~ ;
            load["multiplier"] = multiplier
        end

        # re-evaluate gen
        for gen in values(data["nw"][nw]["gen"])
            if haskey(gen, "isfilter") && gen["isfilter"] == 1
                # do nothing, is handled by constraints constraint_active_filter
                gen["pmin"] = -abs(gen["pmax"])
                gen["qmin"] = -abs(gen["qmax"])
            else #is true generator
                if nw !="1" #cost of harmonics set to 0
                    gen["cost"] *= 0
                    #harmonics can be injected/absorbed to match load
                    gen["pmin"] = -abs(gen["pmax"])
                    gen["qmin"] = -abs(gen["qmax"])
                end
            end

        end

        # re-evaluate the bus data
        for bus in values(data["nw"][nw]["bus"])
            #use fundamental as limit for rms
            bus["vminrms"] = bus["vmin"]
            bus["vmaxrms"] = bus["vmax"]

            # true harmonics don't have minimum voltage
            if nw !="1"
                bus["vmin"] = 0
            end
            rm = Dict(3=> 0.05, 5=>0.06, 7=> 0.05, 9=> 0.015, 11=>0.035)
            # inject relative magnitude limits for harmonics
            bus["rm"] = rm
        end

        # re-evaluate the branch data
        for branch in values(data["nw"][nw]["branch"])
            haskey(branch, "br_r") ? branch["br_r"] *= sqrt(nh) : ~ ;
            haskey(branch, "br_x") ? branch["br_x"] *= nh : ~ ;
            haskey(branch, "b_fr") ? branch["b_fr"] *= nh : ~ ;
            haskey(branch, "b_to") ? branch["b_to"] *= nh : ~ ;
        end

        # re-evaluate the transformer data
        if haskey(data["nw"][nw], "xfmr")
            for xfmr in values(data["nw"][nw]["xfmr"])
                haskey(xfmr, "xsc") ? xfmr["xsc"] *= nh : ~ ;
                haskey(xfmr, "r1")  ? xfmr["r1"] *= sqrt(nh) : ~ ;
                haskey(xfmr, "r2")  ? xfmr["r2"] *= sqrt(nh) : ~ ;

                haskey(xfmr, "xe1") ? xfmr["xe1"] *= nh : ~ ;
                haskey(xfmr, "xe2") ? xfmr["xe2"] *= nh : ~ ;
                haskey(xfmr, "re1") ? xfmr["re1"] *= sqrt(nh) : ~ ;
                haskey(xfmr, "re2") ? xfmr["re2"] *= sqrt(nh) : ~ ;

                xfmr["cnf1"] = haskey(xfmr, "vg") ? uppercase(xfmr["vg"][1]) : "Y" ;
                xfmr["cnf2"] = haskey(xfmr, "vg") ? uppercase(xfmr["vg"][2]) : "Y" ;

                shift = haskey(xfmr, "vg") ? parse(Int, xfmr["vg"][3]) : 0 ;
                if is_pos_sequence(nh)
                    xfmr["tr"] = cosd(-30.0 * shift)
                    xfmr["ti"] = sind(-30.0 * shift)
                elseif is_neg_sequence(nh)
                    xfmr["tr"] = cosd(30.0 * shift)
                    xfmr["ti"] = sind(30.0 * shift)
                elseif is_zero_sequence(nh)
                    xfmr["tr"] = 1.0
                    xfmr["ti"] = 0.0
                end
            end
        end
    end

    return data
end

function collect_harmonics!(data::Dict{String, Any}, harmonics::Array{Int}, xfmr_exc::Dict{String, Any})
    bus_harmonics = [parse(Int,ni[4:end])   for ni in keys(data["bus"]["1"])
                                            if  ni[1:2] == "nh"]
    push!(harmonics, bus_harmonics...)
    if haskey(xfmr_exc, "voltage_harmonics")
        push!(harmonics, xfmr_exc["voltage_harmonics"]...)
    end
    if haskey(xfmr_exc, "current_harmonics")
        push!(harmonics, xfmr_exc["current_harmonics"]...)
    end
    unique!(sort!(harmonics))
end

is_pos_sequence(nh::Int) = nh % 3 == 1
is_neg_sequence(nh::Int) = nh % 3 == 2
is_zero_sequence(nh::Int) = nh % 3 == 0

### TRANSFORMER EXCITATION ###
""
excitation_flux_polar(V, θ, w, t) =
    sum(V[n] ./ w[n] .* sin.(w[n] .* t .+ θ[n]) for n in 1:length(V))
excitation_flux_rectangular(Vre, Vim, w, t) =
    excitation_flux_polar(hypot.(Vre,Vim), atan.(Vim,Vre), w, t)
excitation_current_sigmoid(inom, ψmax, ψ) =
    -inom .* log.(2.0 ./ (ψ ./ ψmax .+ 1.0) .- 1.0)

""
function sample_voltage_polar(voltage_harmonics, dv, vmin, vmax, dθ, θmin, θmax)
    S = reduce(vcat,[[vmin[ni]:dv[ni]:vmax[ni],θmin[ni]:dθ[ni]:θmax[ni]]
                      for (ni,nh) in enumerate(voltage_harmonics)])
    R = [1:length(s) for s in S]
    return S, R
end
function sample_voltage_rectangular(voltage_harmonics, dv, vmin, vmax)
    @show voltage_harmonics, dv, vmin, vmax
    S = reduce(vcat,[[vmin[ni]:dv[ni]:vmax[ni],vmin[ni]:dv[ni]:vmax[ni]]
                      for (ni,nh) in enumerate(voltage_harmonics)])
    R = [1:length(s) for s in S]
    return S, R
end

""
function sample_xfmr_excitation(data::Dict{String, <:Any}, xfmr_exc::Dict{String, <:Any})
    if isempty(xfmr_exc)
        ## HARMONICS
        voltage_harmonics = []
        current_harmonics = []
    else
        ## HARMONICS
        voltage_harmonics = xfmr_exc["voltage_harmonics"]
        current_harmonics = xfmr_exc["current_harmonics"]

        ## TIME
        N = xfmr_exc["N"]

        ## EXCITATION CURRENT
        current_type = xfmr_exc["current_type"]

        ## EXCITATION FLUX
        excitation_type = xfmr_exc["excitation_type"]
        inom = xfmr_exc["inom"]
        ψmax = xfmr_exc["ψmax"]

        ## EXCITATION VOLTAGE
        voltage_type = xfmr_exc["voltage_type"]
        dv   = xfmr_exc["dv"]
        vmin = xfmr_exc["vmin"]
        vmax = xfmr_exc["vmax"]
        dθ   = xfmr_exc["dθ"]
        θmin = xfmr_exc["θmin"]
        θmax = xfmr_exc["θmax"]

        ## TIME
        w = (2.0 * pi * freq) .* voltage_harmonics
        t = 0.0:(1 / (freq * N * maximum(current_harmonics))):(5.0 / freq)

        ## DECOMPOSITION
        fq = _SDC.Sinusoidal(freq .* current_harmonics)

        ## SAMPLE VOLTAGE
        if voltage_type == :polar
            S, R = sample_voltage_polar(voltage_harmonics, dv, vmin, vmax, dθ, θmin, θmax)
        elseif voltage_type == :rectangular
            S, R = sample_voltage_rectangular(voltage_harmonics, dv, vmin, vmax)
        end

        ## INITIALIZE EXCITATION CURRENT DICTIONARIES
        Ia = Dict(nh => zeros(R...) for nh in current_harmonics)
        Ib = Dict(nh => zeros(R...) for nh in current_harmonics)

        ## SAMPLING LOOP
        @showprogress for nr in Iterators.product(R...)
            sample = [S[ni][ns] for (ni,ns) in enumerate(nr)]

            if voltage_type == :polar
                V, θ = sample[1:2:end], sample[2:2:end]
                ψexc = excitation_flux_polar(V, θ, w, t)
            else voltage_type == :rectangular
                Vre, Vim = sample[1:2:end], sample[2:2:end]
                ψexc = excitation_flux_rectangular(Vre, Vim, w, t)
            end

            if excitation_type == :sigmoid
                I_exc = excitation_current_sigmoid(inom, ψmax, ψexc)
            end

            _SDC.decompose(t, I_exc, fq)
            # angle convention is reversed -> therefore - sign in the I, φ expressions below
            if current_type == :polar
                I, φ = fq.A[2:end], -fq.φ[2:end]
                for (ni,nh) in enumerate(current_harmonics)
                    Ia[nh][nr...], Ib[nh][nr...] = I[ni], φ[ni]
                end
            elseif current_type == :rectangular
                Ire, Iim = fq.A[2:end] .* sin.(-fq.φ[2:end]), fq.A[2:end] .* cos.(-fq.φ[2:end])
                for (ni,nh) in enumerate(current_harmonics)
                    Ia[nh][nr...], Ib[nh][nr...] = Ire[ni], Iim[ni]
                end
            end
        end
    end
    reverse_harmonics = Dict(value => key for (key, value) in data["harmonics"])
    current_harmonics_ntws =
        [parse(Int,reverse_harmonics[nc]) for nc in current_harmonics]
    voltage_harmonics_ntws =
        [parse(Int,reverse_harmonics[nv]) for nv in voltage_harmonics]

    method = _INT.BSpline(_INT.Cubic(_INT.Line(_INT.OnGrid())))
    for nw in keys(data["nw"])
        nih = parse(Int, nw)
        nh = data["harmonics"][nw]
        ni = findfirst(x->x==nh, voltage_harmonics)
        if haskey(data["nw"][nw], "xfmr")
            for xfmr in values(data["nw"][nw]["xfmr"])
                xfmr["current_harmonics_ntws"] = current_harmonics_ntws
                xfmr["voltage_harmonics_ntws"] = voltage_harmonics_ntws
                if nh in current_harmonics
                    xfmr["EXC_A"]  = _INT.extrapolate(_INT.scale(_INT.interpolate(Ia[nh], method), S...), _INT.Line())
                    xfmr["EXC_B"]  = _INT.extrapolate(_INT.scale(_INT.interpolate(Ib[nh], method), S...), _INT.Line())
                    xfmr["INT_A"]  = (x...) -> xfmr["EXC_A"](x...)
                    xfmr["INT_B"]  = (x...) -> xfmr["EXC_B"](x...)
                    xfmr["GRAD_A"] = (x...) -> _INT.gradient(xfmr["EXC_A"], x...)
                    xfmr["GRAD_B"] = (x...) -> _INT.gradient(xfmr["EXC_B"], x...)
                end
                if nh in voltage_harmonics
                    @show voltage_harmonics,  vmin
                    xfmr["ert_min"], xfmr["ert_max"] = vmin[ni], vmax[ni]
                    xfmr["eit_min"], xfmr["eit_max"] = vmin[ni], vmax[ni]
                else
                    xfmr["ert_min"], xfmr["ert_max"] = -1.1, 1.1
                    xfmr["eit_min"], xfmr["eit_max"] = -1.1, 1.1
                end
            end
        end
    end
end


text:
function write(f::Function, path::AbstractString, ds::LazIO.LazDataset)
    reader = ds.filehandle
    # create writer
    writer = Ref{Ptr{Cvoid}}(C_NULL)
    LazIO.@check writer[] LazIO.laszip_create(writer)
    # copy header from reader to writer (to be updated later)
    LazIO.@check writer[] LazIO.laszip_set_header(writer[], Ref(ds.header))
    # open writer
    LazIO.@check writer[] LazIO.laszip_open_writer(writer[], path, Cint(1))
    try
        f(writer[])
    finally
        # close files and destroy pointers
        LazIO.@check writer[] LazIO.laszip_close_writer(writer[])
        LazIO.@check writer[] LazIO.laszip_destroy(writer[])
    end
end

function write(f::Function, path::AbstractString, header::LazIO.LazHeader)
    # create writer
    writer = Ref{Ptr{Cvoid}}(C_NULL)
    LazIO.@check writer[] LazIO.laszip_create(writer)
    # copy header to writer (to be updated later)
    LazIO.@check writer[] LazIO.laszip_set_header(writer[], Ref(header))
    # open writer
    LazIO.@check writer[] LazIO.laszip_open_writer(writer[], path, Cint(1))
    try
        f(writer[])
    finally
        # close files and destroy pointers
        LazIO.@check writer[] LazIO.laszip_close_writer(writer[])
        LazIO.@check writer[] LazIO.laszip_destroy(writer[])
    end
end

function writepoint(writer::Ptr{Cvoid}, p::LazIO.LazPoint)
    LazIO.@check writer LazIO.laszip_set_point(writer, Ref(p))
    LazIO.@check writer LazIO.laszip_write_point(writer)
    LazIO.@check writer LazIO.laszip_update_inventory(writer)
end

ext: jl

max_stars_repo_name: evetion/LazIO

max_issues_repo_name: evetion/LazIO

max_forks_repo_name: evetion/LazIO


text:
# update_alloc.jl

export llik_numerator, llik_numerator_Σx_diag;

function lNymat(model::Model_BNP_WMReg_Joint)

    out = zeros(Float64, model.n, model.H)

    if model.γ_type in (:fixed, :global)
        γ_indx = findall(model.state.γ)
    end

    for h = 1:model.H

        if model.γ_type == :local
            γ_indx = findall(model.state.γ[h,:])
        end

        σ2 = deepcopy(model.state.δ_y[h])
        aa = -0.5 * log( 2.0π * σ2 )

        for i = 1:model.n

            μ = deepcopy(model.state.μ_y[h])

            for k in γ_indx
                μ -= model.state.β_y[h, k] * (model.X[i,k] - model.state.μ_x[h, k])
            end

            out[i,h] = aa - 0.5*( model.y[i] - μ )^2 / σ2

        end
    end

    return out
end

function llik_numerator(model::Model_BNP_WMReg_Joint)

    lNy = lNymat(model)
    size(lNy) == size(model.state.lNX) == (model.n, model.H) || throw("lNy and lNX dimension mismatch. lNy: $(size(lNy)) and lNX: $(size(model.state.lNX))")

    ldens_num = lNy .+ model.state.lNX
    lW = broadcast(+, permutedims(model.state.lω), ldens_num)

    return lW # n by H matrix
end

## The alternates below are deprecated, but may be useful sometime?
# function llik_numerator(y::Array{T,1}, X::Array{T,2}, K::Int, H::Int,
#     μ_y::Array{T, 1}, β_y::Array{T, 2}, δ_y::Array{T, 1},
#     μ_x::Array{T, 2},
#     β_x::Union{Array{Array{Float64, 2}, 1}, Nothing},
#     δ_x::Array{Float64,2},
#     γ::Union{BitArray{1}, BitArray{2}}, γδc::Union{Float64, Array{T, 1}, Nothing},
#     lω::Array{T, 1}) where T <: Real

#     if isnothing(β_x) && K > 1 # indicates that Σx_type == :diag
#         return llik_numerator_Σx_diag(y, X, K, H, μ_y, β_y, δ_y, μ_x, δ_x, γ, lω)
#     else
#         yX = hcat(y, X)
#         return llik_numerator(yX, K, H, μ_y, β_y, δ_y, μ_x, β_x, δ_x, γ, γδc, lω)
#     end
# end
# function llik_numerator(yX::Array{T,2}, K::Int, H::Int,
#     μ_y::Array{T, 1}, β_y::Array{T, 2}, δ_y::Array{T, 1},
#     μ_x::Array{T, 2},
#     β_x::Union{Array{Array{Float64, 2}, 1}, Nothing},
#     δ_x::Array{Float64,2},
#     γ::BitArray{1}, γδc::Array{T, 1}, # inflated variance method
#     lω::Array{T, 1}) where T <: Real

#     βγ_y = β_y_modify_γ(β_y, γ)

#     if K > 1
#         βγ_x, δγ_x = βδ_x_modify_γ(β_x, δ_x, γ, γδc)
#     else
#         δγ_x = δ_x_modify_γ(δ_x, γ, γδc)
#     end

#     μ = hcat(μ_y, μ_x)
#     βγ = ( K > 1 ? [βγ_y, βγ_x...] : [βγ_y] )
#     δγ = hcat(δ_y, δγ_x)

#     # the rest could be done in parallel
#     lW = hcat([ lω[h] .+
#                 lNX_sqfChol( Matrix(yX'), μ[h,:], [ βγ[k][h,:] for k = 1:K ], δγ[h,:] )
#                 for h = 1:H
#               ]...) # lW is a n by H matrix

#     return lW
# end
# function llik_numerator(yX::Array{T,2}, K::Int, H::Int,
#     μ_y::Array{T, 1}, β_y::Array{T, 2}, δ_y::Array{T, 1},
#     μ_x::Array{T, 2},
#     β_x::Union{Array{Array{Float64, 2}, 1}, Nothing},
#     δ_x::Array{Float64,2},
#     γ::BitArray{1}, γδc::Float64, # γδc == Inf (subset method)
#     lω::Array{T, 1}) where T <: Real

#     γδc == Inf || throw("Subset method for variable selection requires γδc = Inf.")
#     size(yX,2) == (K+1) || throw("llik_numerator assumes a full X matrix.")

#     γindx = findall(γ)
#     nγ = length(γindx)

#     if nγ == 0

#         lW = hcat( [ lω[h] .+ logpdf.( Normal(μ_y[h], sqrt(δ_y[h])), yX[:,1] )  for h = 1:H ]... )

#     elseif nγ > 0

#         βγ_y = deepcopy(β_y[:, γindx]) # H by nγ matrix

#         if nγ > 1
#             βγ_x, δγ_x = βδ_x_modify_γ(β_x, δ_x, γ, γδc)
#         else
#             δγ_x = δ_x_modify_γ(δ_x, γ, γδc)
#         end

#         μγ = hcat(μ_y, μ_x[:,γindx])
#         βγ = ( nγ > 1 ? [βγ_y, βγ_x...] : [βγ_y] )
#         δγ = hcat(δ_y, δγ_x)

#         yXγ = yX[:,vcat(1, (γindx .+ 1))]

#         # the rest could be done in parallel
#         lW = hcat([ lω[h] .+
#             lNX_sqfChol( Matrix(yXγ'), μγ[h,:], [ βγ[k][h,:] for k = 1:nγ ], δγ[h,:] )
#                 for h = 1:H ]...) # lW is a n by H matrix
#         end

#     return lW # lW is a n by H matrix
# end
# function llik_numerator(yX::Array{T,2}, K::Int, H::Int,
#     μ_y::Array{T, 1}, β_y::Array{T, 2}, δ_y::Array{T, 1},
#     μ_x::Array{T, 2},
#     β_x::Union{Array{Array{Float64, 2}, 1}, Nothing},
#     δ_x::Array{Float64,2},
#     γ::BitArray{1}, γδc::Nothing, # γδc == nothing (integration method)
#     lω::Array{T, 1}) where T <: Real

#     size(yX,2) == (K+1) || throw("llik_numerator assumes a full X matrix.")

#     γindx = findall(γ)
#     γindx_withy = vcat(1, γindx .+ 1)
#     nγ = length(γindx)

#     if nγ == 0

#         lW = hcat( [ lω[h] .+ logpdf.( Normal(μ_y[h], sqrt(δ_y[h])), yX[:,1] )  for h = 1:H ]... )

#     elseif nγ > 0

#         βγ_y = β_y_modify_γ(β_y, γ) # H by K matrix

#         μ = hcat(μ_y, μ_x)
#         β = ( K > 1 ? [βγ_y, β_x...] : [βγ_y] )
#         δ = hcat(δ_y, δ_x)

#         # the rest could be done in parallel
#         Σxs = [ PDMat( sqfChol_to_Σ( [ β[k][h,:] for k = 1:K ], δ[h,:] ).mat[γindx_withy, γindx_withy] ) for h = 1:H ]
#         lW = hcat( [ lω[h] .+ logpdf(MultivariateNormal(μ[h,γindx_withy], Σxs[h]), Matrix(yX[:,γindx_withy]')) for h = 1:H ]... ) # n by H matrix

#         end

#     return lW # lW is a n by H matrix
# end

## local variable selection
function llik_numerator(yX::Array{T,2}, K::Int, H::Int,
    μ_y::Array{T, 1}, β_y::Array{T, 2}, δ_y::Array{T, 1},
    μ_x::Array{T, 2},
    β_x::Union{Array{Array{Float64, 2}, 1}, Nothing},
    δ_x::Array{Float64,2},
    γ::BitArray{2}, # H by K bit matrix
    γδc::Float64, # γδc == Inf (subset method)
    lω::Array{T, 1}) where T <: Real

    γδc == Inf || throw("Subset method for variable selection requires γδc = Inf.")
    size(yX,2) == (K+1) || throw("llik_numerator assumes a full X matrix.")
    size(γ) == (H,K) || throw("Incorrect dimensions for γ.")

    lW = zeros(Float64, size(yX,1), H)
    lenβx = length(β_x)

    # this could be done in parallel
    for h = 1:H

        γindx = findall(γ[h,:])
        nγ = length(γindx)

        if nγ == 0

            lW[:,h] = lω[h] .+ logpdf.( Normal(μ_y[h], sqrt(δ_y[h])), yX[:,1] )

        elseif nγ > 0

            βγ_y = deepcopy(β_y[h, γindx]) # nγ length vector

            if nγ > 1
                βγ_x, δγ_x = βδ_x_h_modify_γ( [ β_x[j][h,:] for j = 1:lenβx ], δ_x[h,:], γ[h,:], γδc )
            else
                δγ_x = δ_x_h_modify_γ(δ_x[h,:], γ[h,:], γδc)
            end

            μγ = vcat(μ_y[h], μ_x[h,γindx])
            βγ = ( nγ > 1 ? [βγ_y, βγ_x...] : [βγ_y] )
            δγ = vcat(δ_y[h], δγ_x)

            yXγ = yX[:,vcat(1, (γindx .+ 1))]

            lW[:,h] =  lω[h] .+ lNX_sqfChol( Matrix(yXγ'), μγ, βγ, δγ )
        end

    end

    return lW # lW is a n by H matrix
end
function llik_numerator(yX::Array{T,2}, K::Int, H::Int,
    μ_y::Array{T, 1}, β_y::Array{T, 2}, δ_y::Array{T, 1},
    μ_x::Array{T, 2},
    β_x::Union{Array{Array{Float64, 2}, 1}, Nothing},
    δ_x::Array{Float64,2},
    γ::BitArray{2}, # H by K bit matrix
    γδc::Nothing, # γδc == nothing (integration method)
    lω::Array{T, 1}) where T <: Real

    size(yX,2) == (K+1) || throw("llik_numerator assumes a full X matrix.")
    size(γ) == (H,K) || throw("Incorrect dimensions for γ.")

    lW = zeros(Float64, size(yX,1), H)
    lenβx = length(β_x)

    # this could be done in parallel
    for h = 1:H

        γindx = findall(γ[h,:])
        γindx_withy = vcat(1, γindx .+ 1)
        nγ = length(γindx)

        if nγ == 0

            lW[:,h] = lω[h] .+ logpdf.( Normal(μ_y[h], sqrt(δ_y[h])), yX[:,1] )

        elseif nγ > 0

            βγ_y = β_y_modify_γ(β_y[h,:], γ[h,:]) # length K vector

            μ = vcat(μ_y[h], μ_x[h,:])
            β = ( K > 1 ? [βγ_y, [ β_x[j][h,:] for j = 1:lenβx ]...] : [βγ_y] )
            δ = vcat(δ_y[h], δ_x[h,:])

            Σx = PDMat( sqfChol_to_Σ( β, δ ).mat[γindx_withy, γindx_withy] )
            lW[:,h] = lω[h] .+ logpdf(MultivariateNormal(μ, Σx), Matrix(yX[:,γindx_withy]'))
        end
    end

    return lW # lW is a n by H matrix
end

function llik_numerator_Σx_diag(y::Array{T,1}, X::Array{T,2}, K::Int, H::Int,
    μ_y::Array{T, 1}, β_y::Array{T, 2}, δ_y::Array{T, 1},
    μ_x::Array{T, 2},
    δ_x::Array{Float64,2},
    γ::BitArray{1}, # compatible with both subset and integration methods
    lω::Array{T, 1}) where T <: Real

    n = length(y)
    γindx = findall(γ)

    lW = zeros(Float64, n, H)

    for i = 1:n
        for h = 1:H
            mean_y = deepcopy(μ_y[h])
            for k in γindx
                mean_y -= β_y[h, k] * (X[i,k] - μ_x[h, k])
                lW[i,h] += -0.5*log( 2.0π * δ_x[h,k] ) - 0.5*( X[i,k] - μ_x[h,k] )^2 / δ_x[h,k]
            end
            lW[i,h] += -0.5*log( 2.0π * δ_y[h] ) - 0.5*( y[i] - mean_y )^2 / δ_y[h]
            lW[i,h] += lω[h]
        end
    end

    return lW # lW is a n by H matrix
end
function llik_numerator_Σx_diag(y::Array{T,1}, X::Array{T,2}, K::Int, H::Int,
    μ_y::Array{T, 1}, β_y::Array{T, 2}, δ_y::Array{T, 1},
    μ_x::Array{T, 2},
    δ_x::Array{Float64,2},
    γ::BitArray{2}, # H by K bit matrix, compatible with both subset and integration methods
    lω::Array{T, 1}) where T <: Real

    size(γ) == (H,K) || throw("Incorrect dimensions for γ.")

    n = length(y)
    γindxes = [ findall(γ[h,:]) for h = 1:H ]

    lW = zeros(Float64, n, H)

    for i = 1:n
        for h = 1:H
            mean_y = deepcopy(μ_y[h])
            for k in γindxes[h]
                mean_y -= β_y[h, k] * (X[i,k] - μ_x[h, k])
                lW[i,h] += -0.5*log( 2.0π * δ_x[h,k] ) - 0.5*( X[i,k] - μ_x[h,k] )^2 / δ_x[h,k]
            end
            lW[i,h] += -0.5*log( 2.0π * δ_y[h] ) - 0.5*( y[i] - mean_y )^2 / δ_y[h]
            lW[i,h] += lω[h]
        end
    end

    return lW # lW is a n by H matrix
end

## Standard Gibbs
# function update_alloc!(model::Model_BNP_WMReg_Joint) where T <: Real

#     lW = llik_numerator(model)

#     ms = maximum(lW, dims=2) # maximum across columns
#     bc_lWmimusms = broadcast(-, lW, ms)
#     W = exp.(bc_lWmimusms)

#     alloc_new = [ sample(StatsBase.Weights(W[i,:])) for i = 1:model.n ]
#     model.state.S = deepcopy(alloc_new)

#     model.state.n_occup = length(unique(alloc_new))

#     return lW # for llik calculation
# end

## Metropolized full conditional update
function update_alloc!(model::Model_BNP_WMReg_Joint) where T <: Real

    lW = llik_numerator(model)

    ms = maximum(lW, dims=2) # maximum across columns
    bc_lWmimusms = broadcast(-, lW, ms)
    W = exp.(bc_lWmimusms)

    ## Metropolized discrete Gibbs
    for i in 1:model.n

        Si = deepcopy(model.state.S[i])
        w_cand = deepcopy(W[i,:])
        w_cand[Si] = 0.0
        cand = sample(StatsBase.Weights(w_cand))
        lar = logsumexp( bc_lWmimusms[i, 1:end .!= Si ] ) - logsumexp( bc_lWmimusms[i, 1:end .!= cand ] )

        if log(rand()) < lar
            model.state.S[i] = deepcopy(cand)
        end

    end

    model.state.n_occup = length(unique(model.state.S))

    return lW # for llik calculation
end

## These alternates are deprecated, but may be useful sometime?
# function update_alloc!(model::Model_BNP_WMReg_Joint, yX::Array{T,2}) where T <: Real

#     lW = llik_numerator(yX, model.K, model.H,
#             model.state.μ_y, model.state.β_y, model.state.δ_y,
#             model.state.μ_x, model.state.β_x, model.state.δ_x,
#             model.state.γ, model.state.γδc, model.state.lω)

#     ms = maximum(lW, dims=2) # maximum across columns
#     bc_lWmimusms = broadcast(-, lW, ms)
#     W = exp.(bc_lWmimusms)

#     alloc_new = [ sample(StatsBase.Weights(W[i,:])) for i = 1:model.n ]
#     model.state.S = deepcopy(alloc_new)

#     model.state.n_occup = length(unique(alloc_new))

#     return lW # for llik calculation
# end
# function update_alloc!(model::Model_BNP_WMReg_Joint, y::Array{T,1}, X::Array{T,2}) where T <: Real

#     lW = llik_numerator(y, X, model.K, model.H,
#             model.state.μ_y, model.state.β_y, model.state.δ_y,
#             model.state.μ_x, model.state.β_x, model.state.δ_x,
#             model.state.γ, model.state.γδc, model.state.lω)

#     ms = maximum(lW, dims=2) # maximum across columns
#     bc_lWmimusms = broadcast(-, lW, ms)
#     W = exp.(bc_lWmimusms)

#     alloc_new = [ sample(StatsBase.Weights(W[i,:])) for i = 1:model.n ]
#     model.state.S = deepcopy(alloc_new)

#     model.state.n_occup = length(unique(alloc_new))

#     return lW # for llik calculation
# end

function lmargy(Λ::PDMat, a::T, b::T) where T <: Real
    return -0.5 * logdet(Λ) + logabsgamma(a)[1] - a * log(b)
end

## Integrate out eta_y and Metropolize
# function update_alloc!(model::Model_BNP_WMReg_Joint) where T <: Real

#     ## These remain fixed throughout
#     Λβ0star_ηy = model.state.Λ0star_ηy * model.state.β0star_ηy
#     βΛβ0star_ηy = PDMats.quad(model.state.Λ0star_ηy, model.state.β0star_ηy)

#     if model.γ_type == :local
#         D_vec = [ construct_Dh(h, model.X, model.state.μ_x[h,:], model.state.γ[h,:] ) for h = 1:model.H ] # this is evaluated for ALL i and h, so it never changes
#     elseif model.γ_type in (:fixed, :global)
#         D_vec = [ construct_Dh(h, model.X, model.state.μ_x[h,:], model.state.γ ) for h = 1:model.H ] # this is evaluated for ALL i and h, so it never changes
#     end

#     ## These will evolve
#     indx_h_vec = [ findall(model.state.S .== h) for h = 1:model.H ]
#     n_h_vec = [ length(indx_h_vec[h]) for h = 1:model.H ]

#     ## Calculate Lam1, a1, b1 for all h; these will evolve
#     a1_vec = 0.5 .* ( model.state.ν_δy .+ n_h_vec )
#     Λ1star_ηy_vec = [ n_h_vec[h] > 0 ? get_Λ1star_ηy_h(D_vec[h][indx_h_vec[h],:], model.state.Λ0star_ηy) : deepcopy(model.state.Λ0star_ηy) for h = 1:model.H ]
#     β1star00_vec = [ n_h_vec[h] > 0 ? (Λβ0star_ηy + D_vec[h][indx_h_vec[h],:]'model.y[indx_h_vec[h]]) : deepcopy(Λβ0star_ηy) for h = 1:model.H ]
#     β1star_ηy_vec = [ n_h_vec[h] > 0 ? Λ1star_ηy_vec[h] \ β1star00_vec[h] : deepcopy(model.state.β0star_ηy) for h = 1:model.H ]
#     b1_vec = [ n_h_vec[h] > 0 ? get_b1_δy_h(model.state.ν_δy, model.state.s0_δy, model.y[indx_h_vec[h]], βΛβ0star_ηy, Λ1star_ηy_vec[h], β1star_ηy_vec[h]) : 0.5*(model.state.ν_δy * model.state.s0_δy)  for h = 1:model.H ]

#     lmargy_n0 = lmargy(model.state.Λ0star_ηy, 0.5*model.state.ν_δy, 0.5*(model.state.ν_δy * model.state.s0_δy)) # this remains constant throughout
#     lmargy_out = [ n_h_vec[h] > 0 ? lmargy(Λ1star_ηy_vec[h], a1_vec[h], b1_vec[h]) : deepcopy(lmargy_n0) for h = 1:model.H ] # this currently has y_i in its component

#     for i = 1:model.n

#         Si = model.state.S[i]

#         ### collect Lam1, a1, b1 for component Si without y_i
#         lmargy_in = deepcopy(lmargy_out)

#         ## create tmps (these will reflect the 'alternate' state for each h=1:H)
#         Λ_tmp_vec = deepcopy(Λ1star_ηy_vec)
#         β1star00_tmp_vec = deepcopy(β1star00_vec)
#         β1star_tmp_vec = deepcopy(β1star_ηy_vec)
#         b1_tmp_vec = deepcopy(b1_vec)

#         ## remove y_i from its component
#         if n_h_vec[Si] == 1 # if it was the only one in the component
#             lmargy_out[ Si ] = deepcopy(lmargy_n0)
#             Λ_tmp_vec[ Si ] = deepcopy(model.state.Λ0star_ηy)
#             β1star00_tmp_vec[ Si ] = deepcopy(Λβ0star_ηy)
#             β1star_tmp_vec[ Si ] = deepcopy(model.state.β0star_ηy)
#             b1_tmp_vec[ Si ] = 0.5*(model.state.ν_δy * model.state.s0_δy)
#         else
#             d_tmp = deepcopy(D_vec[Si][i,:])
#             Λ_tmp_vec[ Si ] = PDMat( Λ_tmp_vec[ Si ] + (-1.0) .* d_tmp * d_tmp' ) # PDMats doesn't have a subtract function
#             β1star00_tmp_vec[ Si ] -= (d_tmp .* model.y[i])
#             β1star_tmp_vec[ Si ] = Λ_tmp_vec[ Si ] \ β1star00_tmp_vec[ Si ]
#             b1_tmp_vec[ Si ] = get_b1_δy_h(model.state.ν_δy, model.state.s0_δy, model.y[setdiff(indx_h_vec[Si], i)], βΛβ0star_ηy,
#                 Λ_tmp_vec[ Si ], β1star_tmp_vec[ Si ] )
#             lmargy_out[ Si ] = lmargy(Λ_tmp_vec[ Si ], a1_vec[ Si ] - 0.5,  b1_tmp_vec[ Si ])
#         end

#         # now lmargy_out and all tmps are as though y_i didn't exist

#         ### evaluate all other Lam1, a1, b1, as if y_i assigned to its component
#         for h in setdiff(1:model.H, Si)
#             d_tmp = deepcopy(D_vec[h][i,:])
#             Λ_tmp_vec[ h ] = PDMat( Λ_tmp_vec[ h ] + d_tmp * d_tmp' )
#             β1star00_tmp_vec[ h ] += (d_tmp .* model.y[i])
#             β1star_tmp_vec[ h ] = Λ_tmp_vec[h] \ β1star00_tmp_vec[h]
#             b1_tmp_vec[ h ] = get_b1_δy_h(model.state.ν_δy, model.state.s0_δy, vcat(model.y[indx_h_vec[h]], model.y[i]), βΛβ0star_ηy,
#                 Λ_tmp_vec[h], β1star_tmp_vec[h]) # may concatenate y since that argument only calulates y'y
#             lmargy_in[ h ] = lmargy(Λ_tmp_vec[h], a1_vec[h] + 0.5, b1_tmp_vec[h])
#         end

#         # now lmargy_in and all tmp[1:H \ Si] are are as though y_i were assigned to its h slot
#         # lmargy_out is still as though y_i didn't exist

#         lw = [ sum(lmargy_out[ 1:end .!= h ]) + lmargy_in[h] for h = 1:model.H ]

#         lw += model.state.lω
#         lw += model.state.lNX[i,:]

#         lw .-= maximum(lw)
#         w = exp.(lw)

#         ## Metropolized discrete Gibbs (Liu, 1996)
#         w_cand = deepcopy(w)
#         w_cand[Si] = 0.0
#         cand = sample(StatsBase.Weights(w_cand))
#         lar = logsumexp( lw[ 1:end .!= Si ] ) - logsumexp( lw[ 1:end .!= cand ] )
#         if log(rand()) < lar
#             # update the running stats
#             indx_h_vec[ Si ] = setdiff( indx_h_vec[ Si ], i ) # take i out of its currently assigned set
#             push!( indx_h_vec[ cand ], i ) # and add it to its newly assigned set (order doesn't matter)
#             n_h_vec[ Si ] -= 1
#             n_h_vec[ cand ] += 1
#             a1_vec[ Si ] -= 0.5
#             a1_vec[ cand ] += 0.5
#             Λ1star_ηy_vec[ Si ] = deepcopy( Λ_tmp_vec[ Si ] )
#             Λ1star_ηy_vec[ cand ] = deepcopy( Λ_tmp_vec[ cand ] )
#             β1star00_vec[ Si ] = deepcopy( β1star00_tmp_vec[ Si ] )
#             β1star00_vec[ cand ] = deepcopy( β1star00_tmp_vec[ cand ] )
#             β1star_ηy_vec[ Si ] = deepcopy( β1star_tmp_vec[ Si ] )
#             β1star_ηy_vec[ cand ] = deepcopy( β1star_tmp_vec[ cand ] )
#             b1_vec[ Si ] = deepcopy( b1_tmp_vec[ Si ] )
#             b1_vec[ cand ] = deepcopy( b1_tmp_vec[ cand ] )

#             lmargy_out[ cand ] = deepcopy( lmargy_in[ cand ] ) # at the start of each iteration of the loop, lmargy_out reflects the current lmargy vector

#             # update S[i]
#             model.state.S[i] = cand
#         else
#             lmargy_out[ Si ] = deepcopy( lmargy_in[ cand ] ) # put y_i back in so lmargy_out reflects the lmargy vector before this iteration of the loop
#             # nothing else?
#         end

#     end

#     model.state.n_occup = length(unique(model.state.S))

#     return nothing
# end


text:
# load the include -
include("Include.jl")

# load the training data -
path_to_training_data = joinpath(_PATH_TO_DATA, "Training-Synthetic-Thrombin-TF-1K.csv")
training_df = CSV.read(path_to_training_data, DataFrame)

# size of training set -
(R,C) = size(training_df)

# build the model structure -
path_to_model_file = joinpath(_PATH_TO_MODEL, "Coagulation.net")
model_buffer = read_model_file(path_to_model_file)

# build the default model structure -
model = build_default_model_dictionary(model_buffer)

# main simulation loop -
SF = 1e9
for i ∈ 1:10

    # build new model -
    dd = deepcopy(model)

    # setup static -
    sfa = dd["static_factors_array"]
    sfa[1] = training_df[i,:TFPI]       # 1 TFPI
    sfa[2] = training_df[i,:AT]         # 2 AT
    sfa[3] = (5e-12) * SF               # 3 TF
    sfa[6] = 0.005                      # 6 TRAUMA

    # grab the multiplier from the data -
    ℳ = dd["number_of_dynamic_states"]
    xₒ = zeros(ℳ)
    xₒ[1] = training_df[i, :II]         # 1 FII
    xₒ[2] = training_df[i, :VII]        # 2 FVII
    xₒ[3] = training_df[i, :V]          # 3 FV
    xₒ[4] = training_df[i, :X]          # 4 FX
    xₒ[5] = training_df[i, :VIII]       # 5 FVIII
    xₒ[6] = training_df[i, :IX]         # 6 FIX
    xₒ[7] = training_df[i, :XI]         # 7 FXI
    xₒ[8] = training_df[i, :XII]        # 8 FXII
    xₒ[9] = (1e-14)*SF                  # 9 FIIa
    dd["initial_condition_vector"] = xₒ

    # generate a κ vector -
    # κ = [
    #     0.061   ; # 1
    #     1.0     ; # 2
    #     1.0     ; # 3
    #     1.0     ; # 4
    #     1.0     ; # 5
    #     1.0     ; # 6
    #     1.0     ; # 7
    #     1.0     ; # 8
    #     0.70    ; # 9
    #     0.1     ; # 10
    #     0.045   ; # 11
    #     -0.065  ; # 12
    # ];

    κ = [
        0.7478700663653719
        1.7676181108566378
        1.6647348616139388
        1.777505405930172
        1.7730138521648193
        1.7272541526011016
        1.691479136564128
        1.7180148264997843
        1.4655502130381282
        0.3795056851485291
        0.6631825691926526
        0.7002036750754579
    ];

    # what is the output array?
    Y = Array{Float64,1}(undef,5)
    Y[1] = training_df[i,:Lagtime]
    Y[2] = training_df[i,:Peak]
    Y[3] = training_df[i,Symbol("T.Peak")]
    Y[4] = training_df[i,:Max]
    Y[5] = training_df[i,:AUC]

    # compute the loss -
    l = loss(κ, Y, dd)

    @show l
end

max_stars_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation

max_issues_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation

max_forks_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation


text:
function ueberpruefeSatzvervollstaendigung(c,d)

if cmp(c,"abgeleitet") == 0 #Null bedeutet richtig, stimmt überein
    if cmp(d,"0") == 0
        println("Eure Antwort für Wort3 und Wort4 sind beide korrekt!");
    elseif cmp(d,"Null") == 0
        println("Eure Antwort für Wort3 und Wort4 sind beide korrekt!");
    elseif cmp(d,"null") == 0
        println("Eure Antwort für Wort3 und Wort4 sind beide korrekt!");
    elseif cmp(d,"zero") == 0
        println("Eure Antwort für Wort3 und Wort4 sind beide korrekt!");
    elseif cmp(d,"Zero") == 0
        println("Eure Antwort für Wort3 und Wort4 sind beide korrekt!");
    else
        println("Eure Antwort für Wort3 ist korrekt, nicht aber die für Wort4. Lest euch noch einmal den Satz durch und überprüft eure Rechtschreibung.")
    end
elseif cmp(c,"abgelitten") == 0 #Null bedeutet richtig, stimmt überein
    if cmp(d,"0") == 0
        println("Eure Antwort für Wort3 und Wort4 sind beide korrekt!");
    elseif cmp(d,"Null") == 0
        println("Eure Antwort für Wort3 und Wort4 sind beide korrekt!");
    elseif cmp(d,"null") == 0
        println("Eure Antwort für Wort3 und Wort4 sind beide korrekt!");
    elseif cmp(d,"zero") == 0
        println("Eure Antwort für Wort3 und Wort4 sind beide korrekt!");
    elseif cmp(d,"Zero") == 0
        println("Eure Antwort für Wort3 und Wort4 sind beide korrekt!");
    else
        println("Eure Antwort für Wort3 ist korrekt, nicht aber die für Wort4. Lest euch noch einmal den Satz durch und überprüft eure Rechtschreibung.")
    end
elseif cmp(c,"differenziert") == 0 #Null bedeutet richtig, stimmt überein
    if cmp(d,"0") == 0
        println("Eure Antwort für Wort3 und Wort4 sind beide korrekt!");
    elseif cmp(d,"Null") == 0
        println("Eure Antwort für Wort3 und Wort4 sind beide korrekt!");
    elseif cmp(d,"null") == 0
        println("Eure Antwort für Wort3 und Wort4 sind beide korrekt!");
    elseif cmp(d,"zero") == 0
        println("Eure Antwort für Wort3 und Wort4 sind beide korrekt!");
    elseif cmp(d,"Zero") == 0
        println("Eure Antwort für Wort3 und Wort4 sind beide korrekt!");
    else
        println("Eure Antwort für Wort3 ist korrekt, nicht aber die für Wort4. Lest euch noch einmal den Satz durch und überprüft eure Rechtschreibung.")
    end
else
    if cmp(d,"0") == 0
        println("Eure Antwort für Wort4 ist korrekt, nicht aber die für Wort3. Lest euch noch einmal den Satz durch und überprüft eure Rechtschreibung.");
    elseif cmp(d,"Null") == 0
        println("Eure Antwort für Wort4 ist korrekt, nicht aber die für Wort3. Lest euch noch einmal den Satz durch und überprüft eure Rechtschreibung.");
    elseif cmp(d,"null") == 0
        println("Eure Antwort für Wort4 ist korrekt, nicht aber die für Wort3. Lest euch noch einmal den Satz durch und überprüft eure Rechtschreibung.");
    elseif cmp(d,"zero") == 0
        println("Eure Antwort für Wort4 ist korrekt, nicht aber die für Wort3. Lest euch noch einmal den Satz durch und überprüft eure Rechtschreibung.");
    elseif cmp(d,"Zero") == 0
        println("Eure Antwort für Wort4 ist korrekt, nicht aber die für Wort3. Lest euch noch einmal den Satz durch und überprüft eure Rechtschreibung.");
    else
        println("Eure Antworten für Wort3 und Wort4 sind noch nicht korrekt. Lest euch noch einmal den Satz durch und überprüft eure Rechtschreibung.")
    end

end

end


text:
using Multisets
using Combinatorics
using DataStructures

abstract type Tree end

mutable struct Branch <: Tree
    left :: Union{Tree,Nothing}
    right :: Union{Tree, Nothing}
    _clades :: Multiset{Int}
    _parent :: Union{Tree, Nothing}
end

mutable struct Leaf <: Tree
    value :: Int
    _clades :: Multiset{Int}
    _parent :: Union{Tree, Nothing}
end

isroot(t) = isnothing(t._parent)
isleaf(x) = isa(x, Leaf)

function Base.show(io::IO, t::Tree)
    print(io, "Rooted[")
    notatetree(io, t)
    println(io, ";]")
end

function Base.:(==)(x :: Tree, y :: Tree)
    if (isleaf(x))
        isleaf(y) && y.value == x.value
    elseif isleaf(y)
        isleaf(x) && x.value == y.value
    else
        (x.left == y.left && x.right == y.right) || (x.left == y.right && x.right == y.left)
    end
end

function Base.hash(x :: Tree, h :: UInt)
    if isleaf(x)
        return hash(x.value, h)
    else
        return hash(hash((x.left, x.right), h) + hash((x.right, x.left), h), h)
    end
end

# reroot at node t by a cut between t and barring
function reroot_exceptchild_(t::Branch, barring::Tree)
    otherchild(t, barring) = t.left == barring ? t.right : t.left
    if isroot(t)
        return otherchild(t, barring)
    else
        oc = otherchild(t, barring)
        rerooted = reroot_exceptchild_(t._parent, t)
        return Branch(oc, rerooted, Multiset{Int}(), nothing)
    end
end

function fixroot_(t::Tree)
    if isa(t, Leaf) return t end
    t.left._parent = t
    t.right._parent = t
    fixroot_(t.left)
    fixroot_(t.right)
    t
end

function destroy_root(t :: Tree)
    t._parent = nothing
    t
end

fixroot = destroy_root ∘ fixroot_

function fixclades(t::Tree)
    if isa(t, Leaf)
        t._clades = Multiset([t.value])
    else
        fixclades(t.left)
        fixclades(t.right)
        t._clades = t.left._clades + t.right._clades
    end
    t
end

reroot_exceptchild = fixclades ∘ fixroot ∘ reroot_exceptchild_

function reroot_root_clades(t :: Branch, barring :: Tree, root_clades :: Multiset{Int})
    tre_wo_clades = reroot_exceptchild_(t, barring) |> fixroot
    tre_wo_clades._clades = root_clades - barring._clades
    tre_wo_clades
end

function cuts(f, t :: Tree)
    let root_clade = nothing
        function cuts_(f, t :: Tree)
            if isleaf(t)

            else
                if (isroot(t))
                    root_clade = t._clades
                    f(t)
                end
                f(t.left)
                f(t.right)
                if !isroot(t)
                    f(reroot_root_clades(t, t.left, root_clade))
                    f(reroot_root_clades(t, t.right, root_clade))
                end
                cuts_(f, t.left)
                cuts_(f, t.right)
            end
        end
        cuts_(f, t)
    end
end

function cladeat(t :: Tree)
    t._clades
end

function cuts(t :: Tree)
    res = Tree[]
    cuts(t) do s
        push!(res, s)
    end
    res
end

function clade_induced_trees(tre, f)
    f(tre._clades)
end

function isset(M :: Multiset{T}) where {T}
    all(M.data[x] <= 1 for x in keys(M.data))
end

function issingly(t :: Tree)
    isset(t._clades)
end

function issingly(t :: Multiset{Int})
    isset(t)
end

unrooted_subtrees = cuts
subtrees = unrooted_subtrees
unrooted_clades(t) = cladeat.(cuts(t))
clades = unrooted_clades

function unrooted_clades_trees(t)
    c = cuts(t)
    zip(cladeat.(c), c)
end

function isbelow(l :: Tree, r :: Tree)
    l._clades < r._clades
end

function issinglysubtree(l :: Tree)
    if isleaf(l) return false end
    if isnothing(l._parent) return false end
    if issingly(l._parent)
        return true
    end
    return false
end

function maximal_elements(coll, rel=isbelow)
    n = coll |> length
    rg = 1:n
    outdegs = zeros(n)
    for (l, r) in combinations(rg, 2)
        lhs = coll[l]
        rhs = coll[r]
        if outdegs[l] == 0 && rel(lhs, rhs)
            outdegs[l] += 1
        end
        if outdegs[r] == 0 && rel(rhs, lhs)
            outdegs[r] += 1
        end
    end
    [e for (i,e) = enumerate(coll) if outdegs[i] == 0 ]
end

max_stars_repo_name: RuneBlaze/clade-phylo

max_issues_repo_name: RuneBlaze/clade-phylo

max_forks_repo_name: RuneBlaze/clade-phylo


text:
export rand_rot, set_rot, get_rot, check_rot, faces, euler_char
export embed_rot, NF, dual

"""
`rand_rot(G::SimpleGraph)` creates a random rotation system for
the graph.
"""
function rand_rot(G::SimpleGraph{T}) where {T}
    d = Dict{T,RingList{T}}()
    for v in G.V
        a = shuffle(RingList(G[v]))
        d[v] = a
    end
    set_rot(G, d)
end

"""
`_default_rot(G::SimpleGraph)` creates a default
rotation system for the graph.
"""
function _default_rot(G::SimpleGraph{T}) where {T}
    d = Dict{T,RingList{T}}()
    for v in G.V
        d[v] = RingList(G[v])
    end
    return d
end

"""
`set_rot(G::SimpleGraph,d)` makes `d` the rotation system
for this graph (held in the graph's cache).

If `d` is omitted, the a default rotation is used.
"""
function set_rot(G::SimpleGraph, d)
    if !check_rot(G, d)
        error("Not a valid rotation system for this graph")
    end
    cache_save(G, :RotationSystem, d)
    nothing
end
set_rot(G::SimpleGraph) = set_rot(G, _default_rot(G))

"""
`get_rot(G::SimpleGraph,v)` returns a `RingList` of the neighbors of `v`.
This assumes that `G` has an associate rotation system. If not, one is
generated and an alert is generated.

`get_rot(G::SimpleGraph)` returns  a copy of the rotation system associated with `G`
(or an error if there is no rotation system).
"""
function get_rot(G::SimpleGraph, v)
    if !has(G, v)
        error("No such vertex $v in this graph")
    end
    d = get_rot(G)
    return d[v]
end
function get_rot(G::SimpleGraph)
    if cache_check(G, :RotationSystem)
        return cache_recall(G, :RotationSystem)
    end

    if cache_check(G, :GraphEmbedding)
        embed_rot(G)
        @info "Giving this graph, $G, a rotation system based on its xy-embedding."
    else
        set_rot(G)
        @info "Giving this graph, $G, a default rotation system."
    end
    return get_rot(G)
end

"""
`embed_rot(G::SimpleGraph)` assigns a rotation system to `G`
corresponding to its current `xy` embedding.
"""
function embed_rot(G::SimpleGraph{T}) where {T}
    xy = getxy(G)
    d = Dict{T,RingList{T}}()
    for v in G.V
        o = Complex(xy[v]...)    # location of v as a complex number
        Nv = G[v]  # neighbors of v as a list

        Z = [Complex(xy[w]...) for w in Nv]  # neighbor locations as complex nums
        th = [mod(angle(o - z), 2pi) for z in Z]  # angles of edges from v
        p = sortperm(th)
        d[v] = RingList(Nv[p])
    end
    set_rot(G, d)
end

"""
`check_rot(G::SimpleGraph,d::Dict)`  checks if `d` is a valid
rotation system for `G`.
"""
function check_rot(G::SimpleGraph{T}, d::Dict{T,RingList{T}}) where {T}
    for v in G.V
        if !haskey(d, v)
            return false
        end
        if Set(d[v]) != Set(G[v])
            return false
        end
    end
    return true
end

function _next_edge(G::SimpleGraph{T}, uv::Tuple{T,T}) where {T}
    u, v = uv
    r = get_rot(G, v)
    w = r(u)
    return (v, w)
end

"""
`_trace_face(G::SimpleGraph, v,w)` uses the graph's rotation system to find a
face starting with the edge `(u,v)` (in that order). Also may be called
by `_trace_face(G,(u,v))`.
"""
function _trace_face(G::SimpleGraph{T}, uv::Tuple{T,T}) where {T}
    data = T[]
    u, v = uv
    push!(data, u)
    push!(data, v)
    e = uv
    while true
        vw = _next_edge(G, e)
        if vw == uv
            break
        end
        push!(data, vw[2])
        e = vw
    end

    boundary = [(data[i], data[i+1]) for i = 1:length(data)-1]

    return RingList(boundary)
end
_trace_face(G::SimpleGraph{T}, u::T, v::T) where {T} = _trace_face(G, (u, v))

"""
`faces(G::SimpleGraph)` returns the set of faces of this graph
(given its rotation system). The rotation system is a
planar embedding iff this returns `2`.

Each face is a `RingList` of the (directed) edges bordering the face.

*Requires that the graph is connected and has at least one edge.*
"""
function faces(G::SimpleGraph{T}) where {T}
    FT = RingList{Tuple{T,T}}   # type of an oriented face
    result = Set{FT}()

    for e in G.E
        u, v = e
        F = _trace_face(G, u, v)
        push!(result, F)
        F = _trace_face(G, v, u)
        push!(result, F)
    end

    return result
end

"""
`NF(G)` returns the number of faces in the graph `G`
given its current rotation system.

*Requires that the graph is connected and has at least one edge.*
"""
NF(G::SimpleGraph) = length(faces(G))

"""
`euler_char(G::SimpleGraph)` computes the Euler characteristic
of the graph `G` with its associated rotation system.
Specifically, `euler_char(G)` returns `NV(G) - NE(G) + NF(G)`.

*Requires that the graph is connected and has at least one edge.*
"""
euler_char(G::SimpleGraph) = NV(G) - NE(G) + NF(G)

_shorten(F::RingList) = first.(F)

"""
`dual(G::SimpleGraph)` returns the dual graph of `G`.
The vertices of the dual are the faces of `G` and they are
adjacent if and only if they share a common edge.

*Requires that the graph is connected and has at least one edge.*
"""
function dual(G::SimpleGraph{T}) where {T}
    Flist = collect(faces(G))    # list of faces of the graph
    FT = RingList{Tuple{T,T}}    # data type of faces
    VT = Vector{T}               # data type of shortened faces
    GG = SimpleGraph{VT}()       # graph to return

    # the faces of G are the vertices of GG
    for f in Flist
        add!(GG, _shorten(f))
    end

    d = Dict{Tuple{T,T},FT}()    # mapping from (directed) edges to faces
    for f in Flist
        for e in f
            d[e] = f
        end
    end

    # add edges
    for e in G.E
        F1 = d[e]
        F2 = d[reverse(e)]
        add!(GG, _shorten(F1), _shorten(F2))
    end

    try
        # create a rotation system
        R = Dict{VT,RingList{VT}}()   # mapping from a face to a cycle of its neighbors
        for f in Flist
            bdy = [_shorten(d[reverse(e)]) for e in f]
            R[_shorten(f)] = RingList(bdy)
        end

        set_rot(GG, R)
    catch
        @info "Unable to transfer rotation system to the dual; assigning a default instead."
        set_rot(GG)
    end

    try
        embed(GG, :tutte)
    catch
    end

    name(GG, "Dual of $(name(G))")

    return GG
end


text:
using DynamicMixtures
using Test

@testset "DynamicMixtures.jl" begin
    include("dynamic.jl")
end


text:
using LaTeXStrings
import AlgebraicMultigrid
import DelimitedFiles
import DPFEHM
import FEHM
import HDF5
import ILUZero
import IterativeSolvers
import JLD
import Preconditioners
import PyPlot

include("utilities.jl")

network_dirs = ["thomas_tpl_p2_x01", "thomas_tpl_p3_x01", "thomas_tpl_p5_x01", "thomas_tpl_p10_x01", "thomas_tpl_p20_x01"]
network_dirs = [network_dirs; ["var_b_0.0", "var_b_0.5", "var_b_1.0"]]
network_dirs = [network_dirs; ["TSA250_50"]]
need_to_download_data = false
for dir in network_dirs
	if !isdir(dir)
		need_to_download_data = true
	end
end
if need_to_download_data
	download("https://zenodo.org/record/5213727/files/results.tar.gz?download=1", "./results.tar.gz")
	run(`tar xzf results.tar.gz`)
	download("https://zenodo.org/record/5213727/files/meshes.tar.gz?download=1", "./meshes.tar.gz")
	run(`tar xzf meshes.tar.gz`)
end
times = Dict()
chs = Dict()
for (i, network_dir) in enumerate(network_dirs)
	resultsfilename = "results/" * replace(network_dir, "/"=>"_") * "_results.jld"
	if isfile(resultsfilename)
		thesetimes, thesechs = JLD.load(resultsfilename, "times", "chs")
		times[network_dir] = thesetimes
		chs[network_dir] = thesechs
	else
		times[network_dir] = Dict()
		chs[network_dir] = Dict()
	end
	@time A, b, others = loaddir(network_dir)
	for (j, (solver_name, solver)) in enumerate(preconditioned_solvers)
		for (k, (preconditioner_name, preconditioner)) in enumerate(preconditioners)
			combo_name = string(solver_name, "+", preconditioner_name)
			if !haskey(times[network_dir], combo_name)
				t = @elapsed begin
					pl = preconditioner(A)
					x, ch = solver(A, b, pl)
				end
				ch = ch / ch[end] * sqrt(sum((A * x .- b) .^ 2) / sum(b .^ 2))
				@show combo_name, t
				times[network_dir][combo_name] = t
				chs[network_dir][combo_name] = ch
				JLD.save(resultsfilename, "times", times[network_dir], "chs", chs[network_dir])
			end
		end
	end
	for (j, (solver_name, solver)) in enumerate(unpreconditioned_solvers)
		if !haskey(times[network_dir], solver_name)
			t = @elapsed begin
				x, ch = solver(A, b)
			end
			@show solver_name, t
			ch = ch / ch[end] * sqrt(sum((A * x .- b) .^ 2) / sum(b .^ 2))
			times[network_dir][solver_name] = t
			chs[network_dir][solver_name] = ch
			JLD.save(resultsfilename, "times", times[network_dir], "chs", chs[network_dir])
		end
	end
end


text:
using PkgBenchmark, Reachability
results = benchmarkpkg("ARCH2018_RE")
export_markdown("results.md", results)


text:
###############################################################################
#
# Recursive FT test
#
###############################################################################
"""
Using consistent floating number type may accelerate calculations, particularly for GPUs. Thus, we provided a function to test the floating number type
    recursively:

$(METHODLIST)

"""
function FT_test end

"""
If the testing variable is an array, the function will test if element type is float number:
- If true, the function tests if the element type is the same as given `FT`
- If false, the function tests each element recursively

    FT_test(para::Array, FT)

Test the the floating point type of a variable, given
- `para` An array
- `FT` given FT

---
Example
```julia
ft_1 = FT_test([1, 2, 3], Float64);
ft_2 = FT_test(Any[1, 1.0f0, 1.0e0], Float64);
ft_3 = FT_test([1, 2.0, "a"], Float64);
```
"""
FT_test(para::Array, FT) =
(
    passed = true;

    # fail if para is float but not FT
    if eltype(para) <: AbstractFloat
        if eltype(para) != FT
            passed = false;
        end;
    else
        if !all(FT_test.(para, FT))
            passed = false;
        end;
    end;

    return passed
)

"""
When the given variable is a number, the following criteria are used based on the type of the number
- If the variable is a float, then compare it to `FT`
- Otherwise, return true (because integers can be mix-used with floats)

    FT_test(para::Number, FT)

Test the the floating point type of a variable, given
- `para` A number
- `FT` given FT

Example
```julia
ft_1 = FT_test(1, Float64);
ft_2 = FT_test(1.0, Float64);
ft_3 = FT_test(1.0f0, Float64);
```
"""
FT_test(para::Number, FT) =
(
    passed = true;

    # fail if para is float but not FT
    if typeof(para) <: AbstractFloat
        if typeof(para) != FT
            passed = false;
        end;
    end;

    return passed
)

"""
If the variable is a function, module, or symbol, then nothing will be done, and a true will be returned. This method is to avoid endless nested loop
    within the test of Function or Module:

    FT_test(para::Union{Function,Module,Symbol})

Test the the floating point type of a variable, given
- `para` An array
- `FT` given FT

---
Example
```julia
f(x) = x^2;
m = CLIMAParameters;

ft_1 = FT_test(f, Float64);
ft_2 = FT_test(m, Float64);
ft_3 = FT_test("haha", Float64);
```
"""
FT_test(para::Union{Function,Module,Symbol}) = true

"""
The variable to test maybe a struct, but `FT_test` does not know the struct type name a priori. Thus, we try to read out the fields of the variable:

- If succeeds, the function test the fields recursively
- If fails, then do nothing

    FT_test(para::Any, FT)

Test the the floating point type of a variable, given
- `para` A function, module, or symbol
- `FT` given FT

---
Example
```julia
struct SA
    a
    b
end
sa = SA(1, 2.0);

ft_1 = FT_test(sa);
ft_2 = FT_test();
```
"""
FT_test(para::Any, FT) =
(
    passed = true;

    # try to detech struct
    if !(typeof(para) <: DataType)
        try
            arr = [];
            for fn in fieldnames( typeof(para) )
                push!(arr, FT_test( getfield(para, fn), FT ));
            end;
            if !all(arr)
                passed = false;
            end;
        catch e
            nothing
        end;
    end;

    return passed
)

###############################################################################
#
# Recursive NaN test
#
###############################################################################
"""
Like [`FT_test`](@ref), same logic is used to test if all the elements within the tested variable are not NaN:

$(METHODLIST)

"""
function NaN_test end

"""
When an array is passed to the function, the function first determines if the variable is an array of number:

- If true, the function tests each element within the array directly
- If false, the function tests each element recursively

    NaN_test(para::Array)

Test the the floating point type of para is not NaN, given
- `para` An array

---
Example
```julia
nan_1 = NaN_test([1,2,NaN]);
nan_2 = NaN_test([1,3,4]);
nan_3 = NaN_test([1,2,"a"]);
```
"""
NaN_test(para::Array) =
(
    passed = true;

    # fail if para is number
    if eltype(para) <: Number
        if !all(.!isnan.(para))
            passed = false;
        end;
    else
        if !all(NaN_test.(para))
            passed = false;
        end;
    end;

    return passed
)

"""
If the variable is a number, the function tests if the number is NaN directly:

    NaN_test(para::Number)

Test if the variable is NaN, given
- `para` A number

---
Example
```julia
nan_1 = NaN_test(1);
nan_2 = NaN_test(NaN);
```
"""
NaN_test(para::Number) = !isnan(para)

"""
If the variable is a function, module, or symbol, then nothing will be done, and a true will be returned. This method is to avoid endless nested loop
    within the test of Function or Module:

    NaN_test(para::Union{Function,Module,Symbol})

Test the the floating point type of a variable, given
- `para` A function, module, or symbol

---
Example
```julia
f(x) = x^2;
m = CLIMAParameters;

nan_1 = NaN_test(f);
nan_2 = NaN_test(m);
nan_3 = NaN_test("haha");
```
"""
NaN_test(para::Union{Function,Module,Symbol}) = true

"""
If the given variable is a struct, the function will test the fields recursively:

    NaN_test(para::Any)

Test if the variable contains any NaN, given
- `para` Tested variable

---
Example
```julia
struct SA
    a
    b
end

nan_1 = NaN_test(SA(1,2));
nan_2 = NaN_test(SA(1,NaN));
```
"""
NaN_test(para::Any) =
(
    passed = true;

    # try to detech struct
    if !(typeof(para) <: DataType)
        try
            arr = [];
            for fn in fieldnames( typeof(para) )
                push!(arr, NaN_test( getfield(para, fn) ));
            end;
            if !all(arr)
                passed = false;
            end;
        catch e
            nothing
        end;
    end;

    return passed
)


text:
#!/usr/bin/env julia
println("Initialising libraries...")
using NeRCA
using Plots
using PlotThemes
using Dates
using Measures
GR.inline("png")

if length(ARGS) < 3
    println("Usage: ./live_royfit.jl DETX LIGIER_PORT TIME_RES")
    exit(1)
end

const calib = Calibration("latest.detx")
const LIGIER_PORT = parse(Int, ARGS[2])
const TIME_RES = ARGS[3]
const DOWNSAMPLE = 0.1  # fraction to keep

function main()
    println("Starting live ROyFit")

    sparams = NeRCA.SingleDURecoParams()

    for message in CHClient(ip"127.0.0.1", LIGIER_PORT, ["IO_EVT"])
        if rand() > DOWNSAMPLE
            print("x")
            continue
        else
            print(".")
        end
        event = read(IOBuffer(message.data), NeRCA.DAQEvent)

        hits = calibrate(calib, event.hits)
        triggered_hits = triggered(hits)
        dus = sort(unique(map(h->h.du, hits)))
        triggered_dus = sort(unique(map(h->h.du, triggered_hits)))
        n_dus = length(dus)
        n_triggered_dus = length(triggered_dus)
        n_doms = length(unique(h->h.dom_id, triggered_hits))

        if n_doms < 4
            continue
        end

        colours = palette(:default)
        plot()
        Q = []
        for (idx, du) in enumerate(dus)
            du_hits = filter(h->h.du == du, hits)
            if length(triggered(du_hits))== 0
                continue
            end
            fit = NeRCA.single_du_fit(du_hits, sparams)
            push!(Q, fit.Q)
            plot!(du_hits, fit, markercolor=colours[idx], label="DU $(du)", max_z=calib.max_z)
            write_time_residuals(TIME_RES, event, fit.selected_hits, fit)
        end
        if sum(Q) < 200 && n_doms > 12 && n_dus > 1
            println("\nPlotting...")
            fit_params = "ROy live reconstruction (combined single line): Q=$([round(_Q,digits=2) for _Q in Q])"
            event_params = "Det ID $(event.det_id), Run $(event.run_id), FrameIndex $(event.timeslice_id), TriggerCounter $(event.trigger_counter), Overlays $(event.overlays)"
            time_params = "$(unix2datetime(event.timestamp)) UTC"
            trigger_params = "Trigger: $(is_mxshower(event) ? "MX " : "")$(is_3dmuon(event) ? "3DM " : "")$(is_3dshower(event) ? "3DS " : "")"
            time_params = "$(unix2datetime(event.timestamp)) UTC"

            plot!(title="$(fit_params)\n$(event_params), $(trigger_params)\n$(time_params)", titlefontsize=8, margin=5mm)

            savefig("plots/ztplot_roy.png")
        end

        #= if n_triggered_dus > 2 && n_doms > 14 =#
        #=     selected_hits = unique(h->h.dom_id, triggered_hits) =#
        #=     println("\nStarting multiline fit with $(n_dus) DUs and $(length(selected_hits)) selected hits") =#
        #=     prefit_track = NeRCA.prefit(selected_hits) =#
        #=     println(prefit_track) =#
        #=     plot(selected_hits, prefit_track) =#
        #=     savefig("plots/ztplot_roy_prefit.png") =#
        #= end =#
    end
end

function write_time_residuals(filename, event, hits, fit)
    if !isfile(filename)
        fobj = open(filename, "w")
        write(fobj, "run,timestamp,du,floor,dom_id,t_res,Q\n")
    else
        fobj = open(filename, "a")
    end
    dγ, ccalc = NeRCA.make_cherenkov_calculator(fit.sdp)
    for hit in hits
        Δt = hit.t - ccalc(hit.pos.z)
        write(fobj, "$(event.run_id),$(event.timestamp),$(hit.du),$(hit.floor),$(hit.dom_id),$(Δt),$(fit.Q)\n")
    end
    close(fobj)
end

main()


text:
using ExoplanetsSysSim
using StatsFuns
using JLD
using DataFrames
#import ExoplanetsSysSim.StellarTable.df
#import ExoplanetsSysSim.StellarTable.usable
import Compat: UTF8String, ASCIIString

## simulation_parameters
function setup_sim_param_christiansen(args::Vector{ASCIIString} = Array{ASCIIString}(0) )   # allow this to take a list of parameter (e.g., from command line)
  sim_param = setup_sim_param_demo()

  set_inactive(sim_param,["log_eta_pl","power_law_P","power_law_r"])
  add_param_fixed(sim_param,"star_table_setup",setup_star_table_christiansen)
  add_param_fixed(sim_param,"stellar_catalog","q1_q17_christiansen.jld")
  add_param_fixed(sim_param,"generate_num_planets",generate_num_planets_christiansen)
  add_param_fixed(sim_param,"generate_period_and_sizes", generate_period_and_sizes_christiansen)
  add_param_fixed(sim_param,"num_targets_sim_pass_one",150000)  # For faster simulated catalogs

  p_lim_arr_num = [0.5, 1.25, 2.5, 5., 10., 20., 40., 80., 160., 320.]
  r_lim_arr_num = [1., 1.25, 1.5, 1.75, 2.]
  rate_tab_init = [0.027 0.149 0.354 1.37 2.56 1.96 2.85 4.77 0.0001 ;
                   0.02 0.063 0.432 0.9 1.72 1.46 2.4 5.42 3.63 ;
                   0.0074 0.032 0.389 0.545 1.32 1.65 1.96 1.69 1.17 ;
                   0.0001 0.017 0.156 0.542 0.777 1.35 0.965 1.22 2.49]*0.01

  add_param_fixed(sim_param, "p_lim_arr", p_lim_arr_num)
  add_param_fixed(sim_param, "r_lim_arr", r_lim_arr_num*ExoplanetsSysSim.earth_radius)
  add_param_active(sim_param,"obs_par", rate_tab_init)

  return sim_param
end

## planetary_system
function generate_num_planets_christiansen(s::Star, sim_param::SimParam)
  const max_tranets_in_sys::Int64 = get_int(sim_param,"max_tranets_in_sys")
  rate_tab::Array{Float64,2} = get_any(sim_param, "obs_par", Array{Float64,2})
  lambda = sum_kbn(rate_tab)
  #println("# lambda= ", lambda)
  ExoplanetsSysSim.generate_num_planets_poisson(lambda,max_tranets_in_sys)
end

function generate_period_and_sizes_christiansen(s::Star, sim_param::SimParam; num_pl::Integer = 1)
  rate_tab::Array{Float64,2} = get_any(sim_param, "obs_par", Array{Float64,2})

  limitP::Array{Float64,1} = get_any(sim_param, "p_lim_arr", Array{Float64,1})
  limitRp::Array{Float64,1} = get_any(sim_param, "r_lim_arr", Array{Float64,1})

  @assert ((length(limitP)-1) == size(rate_tab, 2))
  @assert ((length(limitRp)-1) == size(rate_tab, 1))
  Plist = []
  Rplist = []
  rate_tab_1d = reshape(rate_tab,length(rate_tab))
  #logmaxcuml = logsumexp(rate_tab_1d)
  #cuml = cumsum_kbn(exp(rate_tab_1d-logmaxcuml))
  maxcuml = sum(rate_tab_1d)
  cuml = cumsum_kbn(rate_tab_1d/maxcuml)

  for n in 1:num_pl
    rollp = rand()
    idx = findfirst(x -> x > rollp, cuml)
    i_idx = (idx-1)%size(rate_tab,1)+1
    j_idx = floor(Int64,(idx-1)//size(rate_tab,1))+1
    ### TODO: Keep uniform log sampling here?
    Rp = exp(rand()*(log(limitRp[i_idx+1])-log(limitRp[i_idx]))+log(limitRp[i_idx]))
    P = exp(rand()*(log(limitP[j_idx+1])-log(limitP[j_idx]))+log(limitP[j_idx]))
    push!(Plist, P)
    push!(Rplist, Rp)
  end
  return Plist, Rplist
end

## stellar_table
function setup_christiansen(sim_param::SimParam; force_reread::Bool = false)
  #global df
  df = ExoplanetsSysSim.StellarTable.df
  if haskey(sim_param,"read_stellar_catalog") && !force_reread
     return df
     #return data
  end
  stellar_catalog_filename = convert(ASCIIString,joinpath(Pkg.dir("ExoplanetsSysSim"), "data", convert(ASCIIString,get(sim_param,"stellar_catalog","q1_q17_dr25_stellar.csv")) ) )
  df = setup_christiansen(stellar_catalog_filename)
  add_param_fixed(sim_param,"read_stellar_catalog",true)
  set_star_table(df)
  return df
end

function setup_christiansen(filename::ASCIIString; force_reread::Bool = false)
  #global df, usable
  df = ExoplanetsSysSim.StellarTable.df
  usable = ExoplanetsSysSim.StellarTable.usable
  if ismatch(r".jld$",filename)
  try
    data = load(filename)
    df::DataFrame = data["stellar_catalog"]
    usable::Array{Int64,1} = data["stellar_catalog_usable"]
    set_star_table(df, usable)
  catch
    error(string("# Failed to read stellar catalog >",filename,"< in jld format."))
  end
  else
  try
    df = readtable(filename)
  catch
    error(string("# Failed to read stellar catalog >",filename,"< in ascii format."))
  end

  has_mass = ! (isna(df[:mass]) | isna(df[:mass_err1]) | isna(df[:mass_err2]))
  has_radius = ! (isna(df[:radius]) | isna(df[:radius_err1]) | isna(df[:radius_err2]))
  has_dens = ! (isna(df[:dens]) | isna(df[:dens_err1]) | isna(df[:dens_err2]))
  has_rest = ! (isna(df[:rrmscdpp04p5]) | isna(df[:dataspan]) | isna(df[:dutycycle]))
  in_Q1Q12 = []
  for x in df[:st_quarters]
    subx = string(x)
    subx = ("0"^(17-length(subx)))*subx
    indQ = search(subx, '1')
    if ((indQ < 1) | (indQ > 12))
      push!(in_Q1Q12, false)
    else
      push!(in_Q1Q12, true)
    end
  end
  is_FGK = []
  for x in 1:length(df[:teff])
    if ((df[x,:teff] > 4000.0) & (df[x,:teff] < 7000.0) & (df[x,:logg] > 4.0))
      push!(is_FGK, true)
    else
      push!(is_FGK, false)
    end
  end
  is_usable = has_radius & is_FGK & has_mass & has_rest #& in_Q1Q12 # & has_dens
  if contains(filename,"q1_q12_christiansen.jld")
    is_usable = is_usable & in_Q1Q12
  end
  # See options at: http://exoplanetarchive.ipac.caltech.edu/docs/API_keplerstellar_columns.html
  # TODO SCI DETAIL or IMPORTANT?: Read in all CDPP's, so can interpolate?
  symbols_to_keep = [ :kepid, :mass, :mass_err1, :mass_err2, :radius, :radius_err1, :radius_err2, :dens, :dens_err1, :dens_err2, :rrmscdpp04p5, :dataspan, :dutycycle ]
  delete!(df, [~(x in symbols_to_keep) for x in names(df)])    # delete columns that we won't be using anyway
  usable = find(is_usable)
  df = df[usable, symbols_to_keep]
  set_star_table(df, usable)
  end
  return df
  #global data = convert(Array{Float64,2}, df) # df[usable, symbols_to_keep] )
  #global colid = Dict(zip(names(df),[1:length(names(df))]))
  #return data
end

setup_star_table_christiansen(sim_param::SimParam; force_reread::Bool = false) = setup_christiansen(sim_param, force_reread=force_reread)
setup_star_table_christiansen(filename::ASCIIString; force_reread::Bool = false) = setup_christiansen(filename, force_reread=force_reread)

## summary_statistics
function calc_summary_stats_sim_pass_one_binned_rates(cat_obs::KeplerObsCatalog, cat_phys::KeplerPhysicalCatalog, param::SimParam )      # Version for simulated data, since includes cat_phys
  ssd = Dict{ASCIIString,Any}()
  cache = Dict{ASCIIString,Any}()

  max_tranets_in_sys = get_int(param,"max_tranets_in_sys")    # Demo that simulation parameters can specify how to evalute models, too
  @assert max_tranets_in_sys >= 1
  idx_tranets = find(x::KeplerTargetObs-> length(x.obs) > 0, cat_obs.target)::Array{Int64,1}             # Find indices of systems with at least 1 tranet = potentially detectable transiting planet

  # Count total number of tranets and compile indices for N-tranet systems
  num_tranets = 0
  idx_n_tranets = Vector{Int64}[ Int64[] for m = 1:max_tranets_in_sys]
  for n in 1:max_tranets_in_sys-1
    idx_n_tranets[n] = find(x::KeplerTargetObs-> length(x.obs) == n, cat_obs.target[idx_tranets] )
    num_tranets += n*length(idx_n_tranets[n])
  end
  idx_n_tranets[max_tranets_in_sys] = find(x::KeplerTargetObs-> length(x.obs) >= max_tranets_in_sys, cat_obs.target[idx_tranets] )

  num_tranets += max_tranets_in_sys*length(idx_n_tranets[max_tranets_in_sys])  # WARNING: this means we need to ignore planets w/ indices > max_tranets_in_sys
  num_tranets  = convert(Int64,num_tranets)            # TODO OPT: Figure out why isn't this already an Int.  I may be doing something that prevents some optimizations

  cache["num_tranets"] = num_tranets
  cache["idx_tranets"] = idx_tranets                                   # We can save lists of indices to summary stats for pass 2, even though we won't use these for computing a distance or probability
  #cache["idx_n_tranets"] = idx_n_tranets

  expected_num_detect = 0.0
  expected_num_sys_n_tranets = zeros(max_tranets_in_sys)
  period_list = zeros(num_tranets)
  weight_list = zeros(num_tranets)
  radius_list = zeros(num_tranets)

  n = 1    # tranet id
  for i in idx_tranets
    for j in 1:num_planets(cat_obs.target[i])
      p_tr_and_det = ExoplanetsSysSim.prob_detect(cat_obs.target[i].prob_detect,j)
      expected_num_detect += p_tr_and_det
      (s,p) = cat_obs.target[i].phys_id[j]

      period_list[n] = cat_phys.target[i].sys[s].orbit[p].P
      weight_list[n] = p_tr_and_det
      radius_list[n] = cat_phys.target[i].sys[s].planet[p].radius
      n = n+1
    end
    for k in 1:max_tranets_in_sys
      expected_num_sys_n_tranets[k] += ExoplanetsSysSim.prob_detect_n_planets(cat_obs.target[i].prob_detect,k)
    end
  end
  ssd["expected planets detected"] = expected_num_detect
  ssd["num_sys_tranets"] = expected_num_sys_n_tranets
  ssd["num targets"] = get_int(param,"num_targets_sim_pass_one")
  #println("expected planets = ",expected_num_detect,", num_sys_tranets = ",expected_num_sys_n_tranets,", num targets = ",ssd["num targets"])

  limitP::Array{Float64,1} = get_any(param, "p_lim_arr", Array{Float64,1})
  limitRp::Array{Float64,1} = get_any(param, "r_lim_arr", Array{Float64,1})

  np_bin = zeros((length(limitP)-1) * (length(limitRp)-1))
  np_bin_idx = 1
  for i in 1:(length(limitP)-1)
    P_match = find(x -> ((x > limitP[i]) && (x < limitP[i+1])), period_list)
    for j in 1:(length(limitRp)-1)
      R_match = find(x -> ((x > limitRp[j]) && (x < limitRp[j+1])), radius_list)

      bin_match = intersect(P_match, R_match)

      np_bin[np_bin_idx] = sum(weight_list[bin_match])
      np_bin_idx += 1
    end
  end

  #ssd["expected planets detected"] = sum(np_bin)
  ssd["expected planets table"] = np_bin

  return CatalogSummaryStatistics(ssd, cache)
end

function calc_summary_stats_obs_binned_rates(cat_obs::KeplerObsCatalog, param::SimParam; trueobs_cat::Bool = false)
  ssd = Dict{ASCIIString,Any}()
  cache = Dict{ASCIIString,Any}()

  if !trueobs_cat
    ssd["num targets"] = get_int(param,"num_targets_sim_pass_one")
  else
    ssd["num targets"] = get_int(param,"num_kepler_targets")
  end

  max_tranets_in_sys = get_int(param,"max_tranets_in_sys")    # Demo that simulation parameters can specify how to evalute models, too
  @assert max_tranets_in_sys >= 1
  idx_tranets = find(x::KeplerTargetObs-> length(x.obs) > 0, cat_obs.target)::Array{Int64,1}             # Find indices of systems with at least 1 tranet = potentially detectable transiting planet

  # Count total number of tranets and compile indices for N-tranet systems
  num_tranets = 0
  idx_n_tranets = Vector{Int64}[ Int64[] for m = 1:max_tranets_in_sys]
  for n in 1:max_tranets_in_sys-1
    idx_n_tranets[n] = find(x::KeplerTargetObs-> length(x.obs) == n, cat_obs.target[idx_tranets] )
    num_tranets += n*length(idx_n_tranets[n])
  end
  idx_n_tranets[max_tranets_in_sys] = find(x::KeplerTargetObs-> length(x.obs) >= max_tranets_in_sys, cat_obs.target[idx_tranets] )

  num_tranets += max_tranets_in_sys*length(idx_n_tranets[max_tranets_in_sys])  # WARNING: this means we need to ignore planets w/ indices > max_tranets_in_sys
  if ( length( find(x::KeplerTargetObs-> length(x.obs) > max_tranets_in_sys, cat_obs.target[idx_tranets] ) ) > 0)   # Make sure max_tranets_in_sys is at least big enough for observed systems
    warn("Observational data has more transiting planets in one systems than max_tranets_in_sys allows.")
  end
  num_tranets  = convert(Int64,num_tranets)            # TODO OPT: Figure out why isn't this already an Int.  I may be doing something that prevents some optimizations

  num_sys_tranets = zeros(max_tranets_in_sys)                           # Since observed data, don't need to calculate probabilities.
  for n in 1:max_tranets_in_sys                                         # Make histogram of N-tranet systems
    num_sys_tranets[n] = length(idx_n_tranets[n])
  end
  ssd["num_sys_tranets"] = num_sys_tranets
  ssd["planets detected"] = num_tranets

  period_list = zeros(num_tranets)
  weight_list = zeros(num_tranets)
  radius_list = zeros(num_tranets)

  n = 1    # tranet id
  for i in idx_tranets
    for j in 1:num_planets(cat_obs.target[i])
      period_list[n] = cat_obs.target[i].obs[j].period
      weight_list[n] = 1.0
      radius_list[n] = sqrt(cat_obs.target[i].obs[j].depth)*cat_obs.target[i].star.radius
      n = n+1
    end
  end

  limitP::Array{Float64,1} = get_any(param, "p_lim_arr", Array{Float64,1})
  limitRp::Array{Float64,1} = get_any(param, "r_lim_arr", Array{Float64,1})

  np_bin = zeros((length(limitP)-1) * (length(limitRp)-1))
  np_bin_idx = 1
  for i in 1:(length(limitP)-1)
    P_match = find(x -> ((x > limitP[i]) && (x < limitP[i+1])), period_list)
    for j in 1:(length(limitRp)-1)
      R_match = find(x -> ((x > limitRp[j]) && (x < limitRp[j+1])), radius_list)

      bin_match = intersect(P_match, R_match)

      np_bin[np_bin_idx] = sum(weight_list[bin_match])
      np_bin_idx += 1
    end
  end

  #ssd["planets detected"] = sum(np_bin)
  ssd["planets table"] = np_bin

  return CatalogSummaryStatistics(ssd, cache)
end

## abc_distance
function calc_distance_vector_binned(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, pass::Int64, sim_param::SimParam ; verbose::Bool = false)
  d = Array{Float64}(0)
  if pass == 1
    if verbose
      println("# Summary 1, pass 1: ",summary1)
      println("# Summary 2, pass 1: ",summary2)
    end
    d = zeros(3)
    # Since observed and simulated catalogs can have different summary statistics for the number of planets, prefer detections if avaliable (e.g., after pass2), otherwise use expected (e.g., from pass 1)
    #np1 = haskey(summary1.stat,"planets detected") ? summary1.stat["planets detected"] : summary1.stat["expected planets detected"]
    #np2 = haskey(summary2.stat,"planets detected") ? summary2.stat["planets detected"] : summary2.stat["expected planets detected"]
    #d[1] = dist_L1_abs(np1/summary1.stat["num targets"],np2/summary2.stat["num targets"])    #  Normalize so different statistics weighted appropriately and not dominated by this one
    #println("np1 = ",np1,", np2 = ",np2)
    #println("np1 (normalized) = ",np1/summary1.stat["num targets"],", np2 (normalized) = ",np2/summary2.stat["num targets"],", d[1] = ",d[1])

    np1 = haskey(summary1.stat,"planets table") ? summary1.stat["planets table"] : summary1.stat["expected planets table"]
    np2 = haskey(summary2.stat,"planets table") ? summary2.stat["planets table"] : summary2.stat["expected planets table"]

    np_bin = zeros(length(np1))
    for n in 1:length(np1)
      np_bin[n] = dist_L1_abs(np1[n]/summary1.stat["num targets"], np2[n]/summary2.stat["num targets"])
      #println("True # [Bin ", n,"] = ",np1[n],", Expected # [Bin ", n,"] = ",np2[n])
    end
    d[1] = maximum(np_bin)
    else
    println("# calc_distance_vector_demo doesn't know what to do for pass= ", pass)
  end
  return d
end

## eval_model
function test_christiansen()
  global sim_param_closure = setup_sim_param_christiansen()
  cat_phys = generate_kepler_physical_catalog(sim_param_closure)
  cat_obs = observe_kepler_targets_single_obs(cat_phys,sim_param_closure)
  global summary_stat_ref_closure = calc_summary_stats_obs_demo(cat_obs,sim_param_closure)
  global cat_phys_try_closure  = generate_christiansen_catalog(sim_param_closure)
  global cat_obs_try_closure  = observe_kepler_targets_sky_avg(cat_phys_try_closure,sim_param_closure)
  global summary_stat_try_closure  = calc_summary_stats_sim_pass_one_demo(cat_obs_try_closure,cat_phys_try_closure,sim_param_closure)
  summary_stat_try_closure   = calc_summary_stats_sim_pass_two_demo(cat_obs_try_closure,cat_phys_try_closure,summary_stat_try_closure,sim_param_closure)
  param_guess = make_vector_of_sim_param(sim_xparam_closure)
  evaluate_model_scalar_ret( param_guess)
end

## inverse_detection
function inv_det_prob(param::SimParam)
  cat_obs = setup_actual_planet_candidate_catalog(setup_star_table_christiansen(param), param)
  inv_det_prob(cat_obs, param)
end

function inv_det_prob(cat_obs::KeplerObsCatalog, param::SimParam)
  num_targ = num_usable_in_star_table()
  idx_tranets = find(x::KeplerTargetObs-> length(x.obs) > 0, cat_obs.target)::Array{Int64,1}

  limitP::Array{Float64,1} = get_any(param, "p_lim_arr", Array{Float64,1})
  limitRp::Array{Float64,1} = get_any(param, "r_lim_arr", Array{Float64,1})

  np_bin = zeros((length(limitP)-1) * (length(limitRp)-1))
  w2_bin = zeros((length(limitP)-1) * (length(limitRp)-1))
  cnt_bin = zeros((length(limitP)-1) * (length(limitRp)-1))
  pl_idx = 1

  for i in idx_tranets
    for j in 1:num_planets(cat_obs.target[i])
      pper = cat_obs.target[i].obs[j].period
      prad = sqrt(cat_obs.target[i].obs[j].depth)*cat_obs.target[i].star.radius

      pbin = findfirst(x -> ((pper > limitP[x]) && (pper < limitP[x+1])), collect(1:(length(limitP)-1)))
      rbin = findfirst(x -> ((prad > limitRp[x]) && (prad < limitRp[x+1])), collect(1:(length(limitRp)-1)))
      if (pbin > 0 && rbin > 0)
        cnt_bin[(pbin-1)*(length(limitRp)-1) + rbin] += 1
        pgeo = ExoplanetsSysSim.calc_transit_prob_single_planet_approx(pper, cat_obs.target[i].star.radius, cat_obs.target[i].star.mass)
	pdet = 0.0
	for star_id in 1:num_targ
	  star = SingleStar(star_table(star_id,:radius),star_table(star_id,:mass),1.0, star_id)
	  cdpp = 1.0e-6 * star_table(star_id, :rrmscdpp04p5) * sqrt(4.5/24.0 / ExoplanetsSysSim.LC_duration )
	  contam = 0.0
	  data_span = star_table(star_id, :dataspan)
	  duty_cycle = star_table(star_id, :dutycycle)
	  pl_arr = Array{Planet}( 1)
	  orbit_arr = Array{Orbit}( 1)
          incl = acos(rand()*star.radius*ExoplanetsSysSim.rsol_in_au/ExoplanetsSysSim.semimajor_axis(pper, star.mass))
	  orbit_arr[1] = Orbit(pper, 0., incl, 0., 0., rand()*2.*pi)
	  pl_arr[1] = Planet(prad, 1.0e-6)
	  kep_targ = KeplerTarget([PlanetarySystem(star, pl_arr, orbit_arr)], fill(cdpp,ExoplanetsSysSim.num_cdpp_timescales,ExoplanetsSysSim.num_quarters),contam,data_span,duty_cycle)

	  duration_central = ExoplanetsSysSim.calc_transit_duration(kep_targ,1,1)
	  if duration_central <= 0.
	    continue
	  end
	  ntr = ExoplanetsSysSim.calc_expected_num_transits(kep_targ, 1, 1, param)
	  depth = ExoplanetsSysSim.calc_transit_depth(kep_targ,1,1)
          snr_central = ExoplanetsSysSim.calc_snr_if_transit(kep_targ, depth, duration_central, param, num_transit=ntr)
	  pdet += ExoplanetsSysSim.calc_prob_detect_if_transit(kep_targ, snr_central, param, num_transit=ntr)
	end
        np_bin[(pbin-1)*(length(limitRp)-1) + rbin] += 1.0/pgeo/(pdet/num_targ)
	w2_bin[(pbin-1)*(length(limitRp)-1) + rbin] += (1.0/pgeo/(pdet/num_targ))^2
	println("Planet ",pl_idx," => Bin ", (pbin-1)*(length(limitRp)-1) + rbin, ", C = ", 1.0/pgeo/(pdet/num_targ))
	pl_idx += 1
      end
    end
  end

  println()
  ess_bin = stellar_ess(param)
  println()

  println("Christiansen Rates:")
  for i in 1:(length(limitP)-1)
    for j in 1:(length(limitRp)-1)
      rate_f = np_bin[(i-1)*(length(limitRp)-1) + j]/num_targ*100.
      if cnt_bin[(i-1)*(length(limitRp)-1) + j] > 0.
        println(rate_f,
                " +/- ", rate_f/sqrt(cnt_bin[(i-1)*(length(limitRp)-1) + j]), " %")
      else
        println(rate_f,
                " +/- N/A %")
      end
    end
  end

  println()
  println("New Rates:")
  for i in 1:(length(limitP)-1)
    for j in 1:(length(limitRp)-1)
      rate_f = (1.0+cnt_bin[(i-1)*(length(limitRp)-1) + j])/(1.0+ess_bin[(i-1)*(length(limitRp)-1) + j])*100.
      up_quant = quantile(Gamma(1.0+cnt_bin[(i-1)*(length(limitRp)-1) + j], 1.0/(1.0+ess_bin[(i-1)*(length(limitRp)-1) + j])), 0.8413)*100.
      low_quant = quantile(Gamma(1.0+cnt_bin[(i-1)*(length(limitRp)-1) + j], 1.0/(1.0+ess_bin[(i-1)*(length(limitRp)-1) + j])), 0.1587)*100.
#      println(rate_f,
#                " +/- ", rate_f/sqrt(1.+cnt_bin[(i-1)*(length(limitRp)-1) + j]), " %")
       println(rate_f,
               " + ", up_quant - rate_f,
               " - ", rate_f - low_quant, " %")
    end
  end

end

## stellar catalog ess
function stellar_ess(param::SimParam)
  num_realiz = 100
  setup_star_table_christiansen(param)
  num_targ = num_usable_in_star_table()

  limitP::Array{Float64,1} = get_any(param, "p_lim_arr", Array{Float64,1})
  limitRp::Array{Float64,1} = get_any(param, "r_lim_arr", Array{Float64,1})

  np_bin = zeros((length(limitP)-1) * (length(limitRp)-1))
  w2_bin = zeros((length(limitP)-1) * (length(limitRp)-1))

  for star_id in 1:num_targ
    star = SingleStar(star_table(star_id,:radius),star_table(star_id,:mass),1.0, star_id)
    cdpp = 1.0e-6 * star_table(star_id, :rrmscdpp04p5) * sqrt(4.5/24.0 / ExoplanetsSysSim.LC_duration )
    contam = 0.0
    data_span = star_table(star_id, :dataspan)
    duty_cycle = star_table(star_id, :dutycycle)

    for i_idx in 1:(length(limitP)-1)
      for j_idx in 1:(length(limitRp)-1)
        temp_bin = 0.0
        for n_test in 1:num_realiz
          pper = exp(rand()*(log(limitP[i_idx+1])-log(limitP[i_idx]))+log(limitP[i_idx]))
	  prad = exp(rand()*(log(limitRp[j_idx+1])-log(limitRp[j_idx]))+log(limitRp[j_idx]))

          pgeo = ExoplanetsSysSim.calc_transit_prob_single_planet_approx(pper, star.radius, star.mass)
	  pdet = 0.0

	  pl_arr = Array{Planet}(1)
	  orbit_arr = Array{Orbit}(1)
          incl = acos(rand()*star.radius*ExoplanetsSysSim.rsol_in_au/ExoplanetsSysSim.semimajor_axis(pper, star.mass))
	  orbit_arr[1] = Orbit(pper, 0., incl, 0., 0., rand()*2.*pi)
	  pl_arr[1] = Planet(prad, 1.0e-6)
	  kep_targ = KeplerTarget([PlanetarySystem(star, pl_arr, orbit_arr)], fill(cdpp,ExoplanetsSysSim.num_cdpp_timescales,ExoplanetsSysSim.num_quarters),contam,data_span,duty_cycle)

	  duration_central = ExoplanetsSysSim.calc_transit_duration(kep_targ,1,1)
	  if duration_central <= 0.
	    continue
	  end
	  ntr = ExoplanetsSysSim.calc_expected_num_transits(kep_targ, 1, 1, param)
	  depth = ExoplanetsSysSim.calc_transit_depth(kep_targ,1,1)
          snr_central = ExoplanetsSysSim.calc_snr_if_transit(kep_targ, depth, duration_central, param, num_transit=ntr)
	  #println(semimajor_axis(kep_targ.sys[1], 1))
	  #println(pper, " ", duration, " ", depth, " ", snr, " ", ntr)
	  pdet = ExoplanetsSysSim.calc_prob_detect_if_transit(kep_targ, snr_central, param, num_transit=ntr)

	  temp_bin += (pgeo*pdet)
          #np_bin[(i_idx-1)*(length(limitRp)-1) + j_idx] += (pgeo*pdet)
	  #w2_bin[(i_idx-1)*(length(limitRp)-1) + j_idx] += (pgeo*pdet)^2
        end
        np_bin[(i_idx-1)*(length(limitRp)-1) + j_idx] += temp_bin/num_realiz
      end
    end
    if (rem(star_id, 1000) == 0.)
      println(string("Star #", star_id, " finished"))
    end
  end

  #ess_bin = zeros((length(limitP)-1) * (length(limitRp)-1))

  for i in 1:(length(limitP)-1)
    for j in 1:(length(limitRp)-1)
  #    ess_bin[(i-1)*(length(limitRp)-1) + j] = (np_bin[(i-1)*(length(limitRp)-1) + j]^2)/w2_bin[(i-1)*(length(limitRp)-1) + j]
      println("Period limits: ", limitP[i:i+1], " / Radius limits: ", limitRp[j:j+1]/ExoplanetsSysSim.earth_radius, " / ESS = ", np_bin[(i-1)*(length(limitRp)-1) + j])
    end
  end
  return np_bin
end