text:
@testset "promotion" begin
    for (mr,br) in ((DynamicAxis(10),         OneTo(10)),)
        @testset "$(typeof(br)) lower range" begin
            @test last(promote(br, mr)) == mr
            @test last(promote(mr, br)) == mr
        end
        for (mr2,br2) in ((DynamicAxis(10),         OneTo(10)),
                             )
            if br != br2
                @testset "$(typeof(br)) & $(typeof(br2))" begin
                    br3 = last(promote(br, br2))
                    @test eltype(last(promote(mr, mr2))) == eltype(br3)
                    @test eltype(last(promote(mr2, mr))) == eltype(br3)

                end
            end
        end
    end
    @testset "promote_rule" begin
        for (t1,t2,p) in ((DynamicAxis, OneTo{Int}, OneTo{Int}),
                          (DynamicAxis, OneTo{Int}, OneTo{Int}))
            @testset "$t1 + $t2" begin
                @test (@inferred(promote_rule(t1, t2)) <: p) == true
                @test (@inferred(promote_rule(t2, t1)) <: p) == true
            end
        end
    end
end


text:
#
# simulate phenotypes from a LMM
#
using Revise
using MendelIHT
using SnpArrays
using Random
using GLM
using DelimitedFiles
using Distributions
using LinearAlgebra
using CSV
using DataFrames
using StatsBase
using TraitSimulation
using Knockoffs
BLAS.set_num_threads(1)

"""
    βi ~ Uniform([-0.5, -0.45, ..., -0.05, 0.05, ..., 0.5]) chosen uniformly
on the odd indices (i.e. real genotypes) across genome

k = Number of causal SNPs
p = Total number of SNPs
traits = Number of traits (phenotypes)
overlap = number of pleiotropic SNPs (affects each trait with effect size βi / r)
"""
function simulate_fixed_beta(k::Int, p::Int, traits::Int; overlap::Int=0)
    true_b = zeros(p, traits)
    effect_sizes = collect(0.05:0.05:0.5)
    k_indep = k - 2overlap # pleiotropic SNPs affect 2 phenotypes
    num_causal_snps = k_indep + overlap
    @assert num_causal_snps > 0 "number of causal SNPs should be positive but was $num_causal_snps"
    idx_causal_snps = sample(1:p, num_causal_snps, replace=false)
    @assert length(idx_causal_snps) == num_causal_snps "length(idx_causal_snps) = $(length(idx_causal_snps)) != num_causal_snps"
    shuffle!(idx_causal_snps)

    # pleiotropic SNPs affect 2 phenotypes
    for i in 1:overlap
        j = idx_causal_snps[i]
        rs = sample(1:traits, 2, replace=false)
        for r in rs
            true_b[j, r] = rand(-1:2:1) * effect_sizes[rand(1:10)]
        end
    end
    # non pleiotropic SNPs affect only 1 phenotype
    for i in (overlap+1):length(idx_causal_snps)
        idx = idx_causal_snps[i]
        true_b[idx, rand(1:traits)] = rand(-1:2:1) * effect_sizes[rand(1:10)]
    end

    @assert count(!iszero, true_b) == k "count(!iszero, true_b) = $(count(!iszero, true_b)) != k = $k"

    return true_b
end

"""
Trait covariance matrix is σg * Φ + σe * I where Φ is the GRM.
"""
function simulate_polygenic(
    plinkname::String, k::Int, r::Int;
    seed::Int=2021, σg=0.1, σe=0.9, βoverlap=2,
    )
    # set seed
    Random.seed!(seed)

    # simulate `.bed` file with no missing data
    x = SnpArray(plinkname * ".bed")
    xla = SnpLinAlg{Float64}(x, model=ADDITIVE_MODEL, impute=true, center=true, scale=true)
    n, p = size(x)

    # intercept is the only nongenetic covariate
    Z = ones(n, 1)
    intercepts = zeros(r)' # each trait have 0 intercept

    # simulate β
    B = simulate_fixed_beta(k, p, r, overlap=βoverlap)
    writedlm("sim$(seed)/trueb.txt", B)

    # between trait covariance matrix
    Σ = random_covariance_matrix(r)
    writedlm("sim$(seed)/true_cov.txt", Σ)

    # between sample covariance is identity + GRM
    Φ = readdlm(plinkname * ".grm")
    V = σg * Φ + σe * I

    # simulate y using TraitSimulations.jl (https://github.com/OpenMendel/TraitSimulation.jl/blob/master/src/modelframework.jl#L137)
    vc = @vc Σ ⊗ V
    μ = zeros(n, r)
    μ_null = zeros(n, r)
    LinearAlgebra.mul!(μ_null, Z, intercepts)
    mul!(μ, xla, B)
    BLAS.axpby!(1.0, μ_null, 1.0, μ)
    VCM_model = VCMTrait(Z, intercepts, xla, B, vc, μ)
    Y = Matrix(Transpose(simulate(VCM_model)))

    # simulate using Distributions.jl
    # μ = z * intercepts + xla * B
    # Y = rand(MatrixNormal(μ', Σ, V))

    return xla, Matrix(Z'), B, Σ, Y
end

"""
Computes power and false positive rates
- p: total number of SNPs
- pleiotropic_snps: Indices (or ID) of the true causal SNPs that affect >1 phenotype
- independent_snps: Indices (or ID) of the true causal SNPs that affect exactly 1 phenotype
- signif_snps: Indices (or ID) of SNPs that are significant after testing

returns: pleiotropic SNP's power, independent SNP's power, number of false positives, and false discovery rate
"""
function power_and_fdr(p::Int, pleiotropic_snps, independent_snps, signif_snps)
    pleiotropic_power = length(signif_snps ∩ pleiotropic_snps) / length(pleiotropic_snps)
    independent_power = length(signif_snps ∩ independent_snps) / length(independent_snps)
    correct_snps = pleiotropic_snps ∪ independent_snps
    FP = length(signif_snps) - length(signif_snps ∩ correct_snps) # number of false positives
    TN = p - length(signif_snps) # number of true negatives
    # FPR = FP / (FP + TN)
    FDR = FP / length(signif_snps)
    return pleiotropic_power, independent_power, FP, FDR
end

function make_grm(chr::Int)
    dir = "/scratch/users/bbchu/ukb/subset/"
    cd(dir)
    isfile(dir * "ukb.10k.chr$chr.bed") || error("PLINK file not present!")
    if !isfile(dir * "ukb.10k.chr$chr.grm")
        println("GRM file not present, generating robust GRM")
        Φ = SnpArrays.grm(SnpArray(dir * "ukb.10k.chr$chr.bed"), method=:Robust)
        writedlm("ukb.10k.chr$chr.grm", Φ)
    end
end

function one_simulation(
    k::Int, r::Int;
    seed::Int=2021, σg=0.1, σe=0.9, βoverlap=2,
    path=5:5:50, init_beta=false, model=:polygenic, debias=100, fdr=0.1
    )
    isdir("sim$seed") ? (return nothing) : mkdir("sim$seed")
    plinkname = "/scratch/users/bbchu/ukb/subset/ukb.10k.chr10"
    knockoff_file = "/scratch/users/bbchu/ukb/subset/ukb.10k.merged.chr10"
    grid = path[2] - path[1] - 1

    # simulate data
    Random.seed!(seed)
    if model == :polygenic
        xla, Z, B, Σ, Y = simulate_polygenic(plinkname, k, r,
            seed=seed, σg=σg, σe=σe, βoverlap=βoverlap)
    elseif model == :sparse
        xla, Z, B, Σ, Y = simulate_sparse(plinkname, k, r,
            seed=seed, σg=σg, σe=σe, βoverlap=βoverlap)
    else
        error("model misspecified!")
    end
    cd("sim$seed")
    writedlm("simulated_phenotypes.phen", Y', ',')

    correct_snps = unique([x[1] for x in findall(!iszero, B)])
    pleiotropic_snps, independent_snps = Int[], Int[]
    for snp in correct_snps
        count(x -> abs(x) > 0, @view(B[snp, :])) > 1 ?
            push!(pleiotropic_snps, snp) : push!(independent_snps, snp)
    end
    # snpdata = SnpData("../" * plinkname)
    # pleiotropic_snp_rsid = snpdata.snp_info[pleiotropic_snps, :snpid]
    # independent_snp_rsid = snpdata.snp_info[independent_snps, :snpid]

    # run GEMMA (GRM is precomputed already)
    # run(`cp ../../$(plinkname).bed .`)
    # run(`cp ../../$(plinkname).bim .`)
    # run(`cp ../../$(plinkname).cXX.txt .`)
    # make_GEMMA_fam_file(xla, Y, plinkname)
    # pheno_columns = [string(ri) for ri in 1:r]
    # gemma_time = @elapsed begin
    #     run(`../../gemma -bfile $plinkname -k $(plinkname).cXX.txt -notsnp -lmm 1 -n $pheno_columns -o gemma.sim$seed`)
    # end
    # gemma_pleiotropic_power, gemma_independent_power, gemma_FP, gemma_FPR, gemma_λ =
    #     process_gemma_result("output/gemma.sim$seed.assoc.txt", pleiotropic_snp_rsid, independent_snp_rsid)
    # println("GEMMA time = $gemma_time, pleiotropic power = $gemma_pleiotropic_power, independent power = $gemma_independent_power, FP = $gemma_FP, FDR = $gemma_FDR, gemma_λ=$gemma_λ")
    # mv("output/gemma.sim$seed.assoc.txt", "gemma.sim$seed.assoc.txt")
    # mv("output/gemma.sim$seed.log.txt", "gemma.sim$seed.log.txt")

    # run multivariate IHT
    mIHT_time = @elapsed begin
        mses = cross_validate(plinkname, MvNormal, path=path, phenotypes="simulated_phenotypes.phen";
            init_beta=init_beta, debias=debias)
        k_rough_guess = path[argmin(mses)]
        dense_path = (k_rough_guess - grid):(k_rough_guess + grid)
        mses_new = cross_validate(plinkname, MvNormal, path=dense_path, phenotypes="simulated_phenotypes.phen";
            init_beta=init_beta, debias=debias, cv_summaryfile="miht.cviht.summary.txt")
        iht_result = iht(plinkname, dense_path[argmin(mses_new)], MvNormal, phenotypes="simulated_phenotypes.phen";
            init_beta=init_beta, debias=debias, summaryfile="miht.summary.txt",
            betafile="miht.beta.txt")
    end
    detected_snps = Int[]
    for i in 1:r
        β = iht_result.beta[i, :]
        detected_snps = detected_snps ∪ findall(!iszero, β)
    end
    mIHT_pleiotropic_power, mIHT_independent_power, mIHT_FP, mIHT_FDR = power_and_fdr(size(B, 1), pleiotropic_snps, independent_snps, detected_snps)
    println("multivariate IHT time = $mIHT_time, pleiotropic power = $mIHT_pleiotropic_power, independent power = $mIHT_independent_power, FP = $mIHT_FP, FDR = $mIHT_FDR")

    # run knockoff + multivariate IHT
    mIHT_ko_time = @elapsed begin
        mses = cross_validate(knockoff_file, MvNormal, path=path, phenotypes="simulated_phenotypes.phen";
            init_beta=init_beta, debias=debias)
        k_rough_guess = path[argmin(mses)]
        dense_path = (k_rough_guess - grid):(k_rough_guess + grid)
        mses_new = cross_validate(knockoff_file, MvNormal, path=dense_path, phenotypes="simulated_phenotypes.phen";
            init_beta=init_beta, debias=debias, cv_summaryfile="miht.ko.cviht.summary.txt")
        iht_result = iht(knockoff_file, dense_path[argmin(mses_new)], MvNormal, phenotypes="simulated_phenotypes.phen";
            init_beta=init_beta, debias=debias, summaryfile="miht.ko.summary.txt",
            betafile="miht.knockoff.beta.txt")
    end
    W = coefficient_diff(iht_result.beta', :interleaved)
    τ = threshold(W, fdr)
    detected_snps = findall(W .> τ)
    mIHT_ko_pleiotropic_power, mIHT_ko_independent_power, mIHT_ko_FP, mIHT_ko_FDR = power_and_fdr(size(B, 1), pleiotropic_snps, independent_snps, detected_snps)
    println("multivariate IHT with knockoffs time = $mIHT_ko_time, pleiotropic power = $mIHT_ko_pleiotropic_power, independent power = $mIHT_ko_independent_power, FP = $mIHT_ko_FP, FDR = $mIHT_ko_FDR")

    # run multiple univariate IHT
    # detected_snps = Int[]
    # uIHT_time = @elapsed begin
    #     for trait in 1:r
    #         mses = cross_validate(plinkname, Normal, path=path, phenotypes=trait+5;
    #             init_beta=init_beta, debias=debias)
    #         k_rough_guess = path[argmin(mses)]
    #         dense_path = (k_rough_guess == 5) ? (0:5) : ((k_rough_guess - 4):(k_rough_guess + 4))
    #         mses_new = cross_validate(plinkname, Normal, path=dense_path, phenotypes=trait+5;
    #             init_beta=init_beta, debias=debias, cv_summaryfile="uiht.cviht.summary$trait.txt")
    #         best_k = dense_path[argmin(mses_new)]
    #         if best_k > 0
    #             iht_result = iht(plinkname, best_k, Normal, phenotypes=trait+5;
    #                 init_beta=init_beta, debias=debias, summaryfile="uiht.summary$trait.txt")
    #             β = iht_result.beta
    #         else
    #             β = zeros(size(B, 2))
    #         end

    #         # save results
    #         detected_snps = detected_snps ∪ findall(!iszero, β)
    #         writedlm("univariate_iht_beta$trait.txt", β)
    #     end
    # end
    # uIHT_pleiotropic_power, uIHT_independent_power, uIHT_FP, uIHT_FPR = power_and_fdr(size(B, 1), pleiotropic_snps, independent_snps, detected_snps)
    # println("univariate IHT time = $uIHT_time, pleiotropic power = $uIHT_pleiotropic_power, independent power = $uIHT_independent_power, FP = $uIHT_FP, FPR = $uIHT_FPR")

    # run MVPLINK
    # phenofile = plinkname * ".phen"
    # make_MVPLINK_fam_and_phen_file(xla, Y, plinkname)
    # mvplink_time = @elapsed run(`../../plink.multivariate --bfile $plinkname --noweb --mult-pheno $phenofile --mqfam`)
    # mvPLINK_pleitropic_power, mvPLINK_independent_power, mvPLINK_FP, mvPLINK_FPR, mvPLINK_λ =
    #     process_mvPLINK("plink.mqfam.total", pleiotropic_snps, independent_snps)
    # println("mvPLINK time = $mvplink_time, pleiotropic power = $mvPLINK_pleitropic_power, independent power = $mvPLINK_independent_power, FP = $mvPLINK_FP, FPR = $mvPLINK_FPR, mvPLINK_λ=$mvPLINK_λ")

    # clean up
    # rm("plink.hh", force=true)
    # rm("$(plinkname).fam", force=true)
    # rm("$(plinkname).bed", force=true)
    # rm("$(plinkname).bim", force=true)
    # rm("$(plinkname).cXX.txt", force=true)

    # save summary stats
    n, p = size(xla)
    open("summary.txt", "w") do io
        println(io, "Simulation $seed summary")
        println(io, "n = $n, p = $p, k = $k, r = $r, βoverlap=$βoverlap")
        println(io, "debias=$debias, init_beta=$init_beta")
        model == :polygenic ? println(io, "model = $model, σg=$σg, σe=$σe") : println(io, "model = $model")
        println(io, "")
        println(io, "mIHT time = $mIHT_time seconds, pleiotropic power = $mIHT_pleiotropic_power, independent power = $mIHT_independent_power, FP = $mIHT_FP, FDR = $mIHT_FDR, λ = NaN")
        println(io, "mIHT knockoff time = $mIHT_ko_time seconds, pleiotropic power = $mIHT_ko_pleiotropic_power, independent power = $mIHT_ko_independent_power, FP = $mIHT_ko_FP, FDR = $mIHT_ko_FDR, λ = NaN")
        # println(io, "uIHT time = $uIHT_time seconds, pleiotropic power = $uIHT_pleiotropic_power, independent power = $uIHT_independent_power, FP = $uIHT_FP, FDR = $uIHT_FDR, λ = NaN")
        # println(io, "mvPLINK time = $mvplink_time seconds, pleiotropic power = $mvPLINK_pleitropic_power, independent power = $mvPLINK_independent_power, FP = $mvPLINK_FP, FDR = $mvPLINK_FDR, λ = $mvPLINK_λ")
        # println(io, "GEMMA time = $gemma_time seconds, pleiotropic power = $gemma_pleiotropic_power, independent power = $gemma_independent_power, FP = $gemma_FP, FDR = $gemma_FDR, λ = $gemma_λ")
    end
    cd("../")

    return nothing
end

function run_simulation(set::Int, model::Symbol)
    σg = 0.1
    σe = 0.9
    fdr = 0.1
    init_beta = true
    debias = false
    βoverlap = [3, 5, 7]
    k = [10, 20, 100]
    r = [2, 3, 3]
    path = set ≥ 3 ? (10:10:200) : (5:5:50)

    println("Simulation model = $model, set $set has k = $(k[set]), r = $(r[set]), βoverlap = $(βoverlap[set])")

    cur_dir = pwd() * "/set$set"
    isdir(cur_dir) || mkdir(cur_dir)
    k_cur = k[set]
    r_cur = r[set]
    βoverlap_cur = βoverlap[set]

    for seed in 1:100
        try
            cd(cur_dir)
            one_simulation(k_cur, r_cur, seed = seed, path = path, βoverlap=βoverlap_cur,
                σg=σg, σe=σe, init_beta=init_beta, model=model, debias=debias)
        catch e
            bt = catch_backtrace()
            msg = sprint(showerror, e, bt)
            println("set $set sim $seed threw an error!")
            println(msg)
            continue
        end
    end
end

# set = parse(Int, ARGS[1])
# model = :polygenic
# run_simulation(set, model)

# set = 1
# model = :polygenic
# σg = 0.1
# σe = 0.9
# fdr=0.1
# init_beta = true
# debias = false
# path = 5:5:50
# βoverlap = 3
# k = 10
# r = 2

# cur_dir = pwd() * "/set$set"
# isdir(cur_dir) || mkdir(cur_dir)
# seed = 1
# cd(cur_dir)

ext: jl


text:
using DynamicalSystemsBase

function eom_ar1_unidir(x, p, n)
    a₁, b₁, c_xy, σ = (p...,)
    x, y = (x...,)
    ξ₁ = rand(Normal(0, σ))
    ξ₂ = rand(Normal(0, σ))

    dx = a₁*x + ξ₁
    dy = b₁*y + c_xy*x + ξ₂
    return SVector{2}(dx, dy)
end

function ar1_unidir(;uᵢ = rand(2), a₁ = 0.90693, b₁ = 0.40693, c_xy = 0.5, σ = 0.40662)
    p = [a₁, b₁, c_xy, σ]
    DiscreteDynamicalSystem(eom_ar1_unidir, uᵢ, p)
end

vars = (1, 2)
npts, tstep = 50, 50
d_xind = Uniform(2.5, 5.5)
d_yind = Uniform(2.5, 5.5)
d_xval = Uniform(0.01, 0.2)
d_yval = Uniform(0.01, 0.2)

X, Y = example_uncertain_indexvalue_datasets(ar1_unidir(c_xy = 0.5), npts, vars, tstep = tstep,
    d_xind = d_xind, d_yind = d_yind,
    d_xval = d_xval, d_yval = d_yval);

time_grid = -20:100:2540
n_draws = 10000 # draws per uncertain value
n_bins = length(time_grid) - 1

# Values in each bin represented as RawValues
b = BinnedResampling(RawValues, time_grid, n_draws)
bc, vs = bin(Y, b);

@test vs isa Vector{Vector{T}} where T
@test length(vs) == n_bins
@test sum(length.(vs)) == length(Y)*n_draws

# Values in each bin represented as UncertainScalarKDE
b_kde = BinnedResampling(UncertainScalarKDE, time_grid, n_draws)
Y_binned = bin(Y, b_kde);

@test Y_binned isa AbstractUncertainIndexValueDataset

# Values in each bin represented as UncertainScalarPopulation
b_pop = BinnedResampling(UncertainScalarPopulation, time_grid, n_draws)
Y_binned = bin(Y, b_pop);

@test Y_binned isa AbstractUncertainIndexValueDataset


text:
reshapeVec2Mat(vec::Vector, rows::Int) = reshape(vec, rows, round(Int,length(vec)/rows))

# """
#     $SIGNATURES
#
# Retrieve data structure stored in a node.
# """
# getData(v::DFGFactor)::GenericFunctionNodeData = v.data
# getData(v::DFGVariable; solveKey::Symbol=:default)::VariableNodeData = v.solverDataDict[solveKey]
# # For Bayes tree

# still used for Bayes Tree
import DistributedFactorGraphs: getData

getData(v::Graphs.ExVertex) = v.attributes["data"]

"""
    $SIGNATURES

Retrieve data structure stored in a variable.
"""
function getVariableData(dfg::AbstractDFG, lbl::Symbol; solveKey::Symbol=:default)::VariableNodeData
  return solverData(getVariable(dfg, lbl, solveKey=solveKey))
end

"""
    $SIGNATURES

Retrieve data structure stored in a factor.
"""
function getFactorData(dfg::T, lbl::Symbol)::GenericFunctionNodeData where {T <: AbstractDFG}
  return solverData(getFactor(dfg, lbl))
end
# TODO -- upgrade to dedicated memory location in Graphs.jl
# see JuliaArchive/Graphs.jl#233

# TODO: Intermediate for refactor. I'm sure we'll see this in 2024 though, it being 'temporary' and all :P
function setData!(v::DFGVariable, data::VariableNodeData; solveKey::Symbol=:default)::Nothing
  v.solverDataDict[solveKey] = data
  return nothing
end
function setData!(f::DFGFactor, data::GenericFunctionNodeData)::Nothing
  f.data = data
  return nothing
end
# For Bayes tree
function setData!(v::Graphs.ExVertex, data)
  # this is a memory gulp without replacement, old attr["data"] object is left to gc
  v.attributes["data"] = data
  nothing
end

## has been moved to DFG
# import DistributedFactorGraphs: getSofttype
# """
#    $(SIGNATURES)
#
# Variable nodes softtype information holding a variety of meta data associated with the type of variable stored in that node of the factor graph.
# """
# function getSofttype(vnd::VariableNodeData)
#   return vnd.softtype
# end
# function getSofttype(v::DFGVariable; solveKey::Symbol=:default)
#   return getSofttype(getData(v, solveKey=solveKey))
# end

"""
    $SIGNATURES
Return the manifolds on which variable `sym::Symbol` is defined.
"""
getManifolds(vd::VariableNodeData) = getSofttype(vd).manifolds
getManifolds(v::DFGVariable; solveKey::Symbol=:default) = getManifolds(solverData(v, solveKey))
function getManifolds(dfg::G, sym::Symbol; solveKey::Symbol=:default) where {G <: AbstractDFG}
  return getManifolds(getVariable(dfg, sym), solveKey=solveKey)
end
getManifolds(vartype::InferenceVariable) = vartype.manifolds
getManifolds(vartype::Type{<: InferenceVariable}) = getManifolds(vartype())

"""
    $(SIGNATURES)

Convenience function to get point values sampled i.i.d from marginal of `lbl` variable in the current factor graph.
"""
getVal(v::DFGVariable; solveKey::Symbol=:default) = v.solverDataDict[solveKey].val
getVal(v::DFGVariable, idx::Int; solveKey::Symbol=:default) = v.solverDataDict[solveKey].val[:,idx]
getVal(vnd::VariableNodeData) = vnd.val
getVal(vnd::VariableNodeData, idx::Int) = vnd.val[:, idx]
function getVal(dfg::T, lbl::Symbol; solveKey::Symbol=:default) where {T <: AbstractDFG}
  return getVariable(dfg, lbl).solverDataDict[solveKey].val
end

"""
    $(SIGNATURES)

Get the number of points used for the current marginal belief estimate represtation for a particular variable in the factor graph.
"""
function getNumPts(v::DFGVariable; solveKey::Symbol=:default)::Int
  return size(solverData(v, solveKey).val,2)
end

# import DistributedFactorGraphs: getfnctype
# # TODO: Refactor - was is das?
# function getfnctype(data::GenericFunctionNodeData)
#   if typeof(data).name.name == :VariableNodeData
#     return VariableNodeData
#   end
#   return data.fnc.usrfnc!
# end
#
# function getfnctype(fact::DFGFactor; solveKey::Symbol=:default)
#   data = getData(fact) # TODO , solveKey=solveKey)
#   return getfnctype(data)
# end
#
# function getfnctype(dfg::T, lbl::Symbol; solveKey::Symbol=:default) where T <: AbstractDFG
#   getfnctype(getFactor(dfg, exvertid))
# end

function getBW(vnd::VariableNodeData)
  return vnd.bw
end

# setVal! assumes you will update values to database separate, this used for local graph mods only
function getBWVal(v::DFGVariable; solveKey::Symbol=:default)
  return solverData(v, solveKey).bw
end
function setBW!(vd::VariableNodeData, bw::Array{Float64,2}; solveKey::Symbol=:default)::Nothing
  vd.bw = bw
  nothing
end
function setBW!(v::DFGVariable, bw::Array{Float64,2}; solveKey::Symbol=:default)::Nothing
  setBW!(solverData(v, solveKey), bw)
  nothing
end

function setVal!(vd::VariableNodeData, val::Array{Float64,2})::Nothing
    vd.val = val
    nothing
end
function setVal!(v::DFGVariable, val::Array{Float64,2}; solveKey::Symbol=:default)::Nothing
    setVal!(solverData(v, solveKey), val)
    nothing
end
function setVal!(vd::VariableNodeData, val::Array{Float64,2}, bw::Array{Float64,2})::Nothing
    setVal!(vd, val)
    setBW!(vd, bw)
    nothing
end
function setVal!(v::DFGVariable, val::Array{Float64,2}, bw::Array{Float64,2}; solveKey::Symbol=:default)::Nothing
  setVal!(v, val, solveKey=solveKey)
  setBW!(v, bw, solveKey=solveKey)
  nothing
end
function setVal!(vd::VariableNodeData, val::Array{Float64,2}, bw::Vector{Float64}; solveKey::Symbol=:default)
  setVal!(vd, val, reshape(bw,length(bw),1))
  nothing
end
function setVal!(v::DFGVariable, val::Array{Float64,2}, bw::Vector{Float64}; solveKey::Symbol=:default)
  setVal!(solverData(v, solveKey=solveKey), val, bw)
  nothing
end
function setVal!(dfg::AbstractDFG, sym::Symbol, val::Array{Float64,2}; solveKey::Symbol=:default)
  setVal!(getVariable(dfg, sym), val, solveKey=solveKey)
end

"""
    $SIGNATURES

Set the point centers and bandwidth parameters of a variable node, also set `isInitialized=true` if `setinit::Bool=true` (as per default).

Notes
- `initialized` is used for initial solve of factor graph where variables are not yet initialized.
- `inferdim` is used to identify if the initialized was only partial.
"""
function setValKDE!(vd::VariableNodeData,
                    pts::Array{Float64,2},
                    bws::Vector{Float64},
                    setinit::Bool=true,
                    inferdim::Float64=0 )::Nothing
  #
  setVal!(vd, pts, bws) # BUG ...al!(., val, . ) ## TODO -- this can be a little faster
  setinit ? (vd.initialized = true) : nothing
  vd.inferdim = inferdim
  nothing
end

function setValKDE!(vd::VariableNodeData,
                    p::BallTreeDensity,
                    setinit::Bool=true,
                    inferdim::Union{Float32, Float64, Int32, Int64}=0 )
  #
  pts = getPoints(p)
  bws = getBW(p)[:,1]
  setValKDE!(vd,pts,bws,setinit,inferdim )
  nothing
end

function setValKDE!(vd::VariableNodeData,
                    val::Array{Float64,2},
                    setinit::Bool=true,
                    inferdim::Float64=0 )::Nothing
  # recover softtype information
  sty = getSofttype(vd)
  p = AMP.manikde!(val, sty.manifolds)
  setValKDE!(vd, p, setinit, inferdim)
  nothing
end

function setValKDE!(v::DFGVariable,
                    val::Array{Float64,2},
                    bws::Vector{Float64},
                    setinit::Bool=true,
                    inferdim::Float64=0;
                    solveKey::Symbol=:default)::Nothing
  # recover softtype information
  setValKDE!(solverData(v, solveKey),val, bws, setinit, inferdim )

  nothing
end

function setValKDE!(v::DFGVariable,
                    val::Array{Float64,2},
                    setinit::Bool=true,
                    inferdim::Float64=0;
                    solveKey::Symbol=:default)::Nothing
  # recover softtype information
  setValKDE!(solverData(v, solveKey),val, setinit, inferdim )
  nothing
end
function setValKDE!(v::DFGVariable,
                    em::EasyMessage,
                    setinit::Bool=true,
                    # inferdim::Union{Float32, Float64, Int32, Int64}=0;
                    solveKey::Symbol=:default  )::Nothing
  #
  setValKDE!(v, em.pts, em.bws, setinit, em.inferdim, solveKey=solveKey)
  nothing
end
function setValKDE!(v::DFGVariable,
                    p::BallTreeDensity,
                    setinit::Bool=true,
                    inferdim::Union{Float32, Float64, Int32, Int64}=0;
                    solveKey::Symbol=:default  )
  #
  setValKDE!(solverData(v,solveKey),p,setinit,Float64(inferdim))
  nothing
end
function setValKDE!(dfg::G,
                    sym::Symbol,
                    p::BallTreeDensity,
                    setinit::Bool=true,
                    inferdim::Union{Float32, Float64, Int32, Int64}=0;
                    solveKey::Symbol=:default  ) where G <: AbstractDFG
    #
    setValKDE!(getVariable(dfg, sym), p, setinit, inferdim, solveKey=solveKey)
    nothing
end

# TODO: Confirm this is supposed to be a variable?
function setVal!(v::DFGVariable, em::EasyMessage; solveKey::Symbol=:default)
    @warn "setVal! deprecated, use setValKDE! instead"
    setValKDE!(v, em, solveKey=solveKey)
end
function setVal!(v::DFGVariable, p::BallTreeDensity; solveKey::Symbol=:default)
    @warn "setVal! deprecated, use setValKDE! instead"
    setValKDE!(v, p, solveKey=solveKey)
end

"""
    $(SIGNATURES)

Construct a BallTreeDensity KDE object from an IIF.EasyMessage object.

Related

manikde!, getKDE, getKDEMax, getKDEMean, EasyMessage
"""
function kde!(em::EasyMessage)
  return AMP.manikde!(em.pts, em.bws, em.manifolds)
end

"""
    $SIGNATURES

Set variable initialized status.
"""
function setVariableInitialized!(varid::VariableNodeData,
                                 status::Bool)
  #
  varid.initialized = status
end

setVariableInitialized!(vari::DFGVariable, status::Bool) = setVariableInitialized!(solverData(vari), status)

"""
    $SIGNATURES

Set method for the inferred dimension value in a variable.
"""
setVariableInferDim!(varid::VariableNodeData, val::Real) = varid.inferdim = convert(Float64,val)
setVariableInferDim!(vari::DFGVariable, val::Real) = setVariableInferDim!(solverData(vari), val)

"""
    $SIGNATURES

Reset the solve state of a variable to uninitialized/unsolved state.
"""
function resetVariable!(varid::VariableNodeData;
                        solveKey::Symbol=:default  )::Nothing
  #
  val = getKDE(varid)
  pts = getPoints(val)
  # TODO not all manifolds will initialize to zero
  fill!(pts, 0.0)
  pn = manikde!(pts, zeros(KDE.Ndim(val)), getManifolds(varid))
  setValKDE!(varid, pn, false, 0.0)
  # setVariableInferDim!(varid, 0)
  # setVariableInitialized!(vari, false)
  nothing
end

resetVariable!(vari::DFGVariable; solveKey::Symbol=:default  )::Nothing = resetVariable!(solverData(vari), solveKey=solveKey)

function resetVariable!(dfg::G,
                        sym::Symbol;
                        solveKey::Symbol=:default  )::Nothing where G <: AbstractDFG
  #
  resetVariable!(getVariable(dfg, sym), solveKey=solveKey)
end

# TODO -- there should be a better way, without retrieving full vertex
# TODO -- Deprecated for DFG -- must update
function getOutNeighbors(dfg::T, v::V; needdata::Bool=false, ready::Union{Nothing, Int}=nothing, backendset::Union{Nothing, Int}=nothing)::Vector{Symbol} where {T <: AbstractDFG, V <: DFGNode}
  @warn "TODO: needdata is currently ignored. Symbols are returned."
  nodes = getNeighbors(dfg, v, ready=ready, backendset=backendset)
  return nodes
end
function getOutNeighbors(dfg::T, vertSym::Symbol; needdata::Bool=false, ready::Int=1, backendset::Int=1 )::Vector{Symbol} where {T <: AbstractDFG, V <: DFGNode}
  @warn "TODO: needdata is currently ignored. Symbols are returned."
  nodes = getNeighbors(dfg, vertSym, ready=ready, backendset=backendset)
  return nodes
end

function setDefaultNodeData!(v::DFGVariable,
                             dodims::Int,
                             N::Int,
                             dims::Int;
                             gt=Dict(),
                             initialized::Bool=true,
                             dontmargin::Bool=false,
                             softtype=nothing)::Nothing
  # TODO review and refactor this function, exists as legacy from pre-v0.3.0
  # this should be the only function allocating memory for the node points (unless number of points are changed)
  data = nothing
  if initialized

      pN = AMP.manikde!(randn(dims, N), softtype.manifolds);

    sp = Int[0;] #round.(Int,range(dodims,stop=dodims+dims-1,length=dims))
    gbw = getBW(pN)[:,1]
    gbw2 = Array{Float64}(undef, length(gbw),1)
    gbw2[:,1] = gbw[:]
    pNpts = getPoints(pN)
    #initval, stdev
    setSolverData(v, VariableNodeData(pNpts,
                            gbw2, Symbol[], sp,
                            dims, false, :_null, Symbol[], softtype, true, 0.0, false, dontmargin))
  else
    sp = round.(Int,range(dodims,stop=dodims+dims-1,length=dims))
    setSolverData(v, VariableNodeData(zeros(dims, N),
                            zeros(dims,1), Symbol[], sp,
                            dims, false, :_null, Symbol[], softtype, false, 0.0, false, dontmargin))
  end
  return nothing
end
# if size(initval,2) < N && size(initval, 1) == dims
#   @warn "setDefaultNodeData! -- deprecated use of stdev."
#   p = AMP.manikde!(initval,diag(stdev), softtype.manifolds);
#   pN = resample(p,N)
# if size(initval,2) < N && size(initval, 1) != dims
  # @info "Node value memory allocated but not initialized"
# else
#   pN = AMP.manikde!(initval, softtype.manifolds)
# end
# dims = size(initval,1) # rows indicate dimensions

"""
    $SIGNATURES

Reference data can be stored in the factor graph as a super-solve.

Notes
- Intended as a mechanism to store reference data alongside the numerical computations.
"""
function setVariableRefence!(dfg::AbstractDFG,
                             sym::Symbol,
                             val::Array{Float64,2};
                             refKey::Symbol=:reference)
  #
  # which variable to update
  var = getVariable(dfg, sym)

  # Construct an empty VND object
  vnd = VariableNodeData(val,
                         zeros(getDimension(var),1),
                         Symbol[],
                         Int[0;],
                         getDimension(var),
                         false,
                         :_null,
                         Symbol[],
                         getSofttype(var),
                         true,
                         0.0,
                         false,
                         true  )
  #
  # set the value in the DFGVariable
  setSolverData(var, vnd, refKey)
end

"""
$(SIGNATURES)

Add a variable node `lbl::Symbol` to `fg::FactorGraph`, as `softtype<:InferenceVariable`.

Example
-------

```julia
fg = initfg()
addVariable!(fg, :x0, Pose2)
```
"""
function addVariable!(dfg::AbstractDFG,
                      lbl::Symbol,
                      softtype::InferenceVariable;
                      N::Int=100,
                      autoinit::Bool=true,  # does init need to be separate from ready? TODO
                      solvable::Int=1,
                      dontmargin::Bool=false,
                      labels::Vector{Symbol}=Symbol[],
                      smalldata=Dict{String, String}(),
                      checkduplicates::Bool=true  )::DFGVariable
  #
  v = DFGVariable(lbl, softtype)
  v.solvable = solvable
  # v.backendset = backendset
  v.tags = union(labels, Symbol.(softtype.labels), [:VARIABLE])
  v.smallData = smalldata
  setDefaultNodeData!(v, 0, N, softtype.dims, initialized=!autoinit, softtype=softtype, dontmargin=dontmargin) # dodims
  DFG.addVariable!(dfg, v)

  return v
end

function addVariable!(dfg::G,
                      lbl::Symbol,
                      softtype::Type{<:InferenceVariable};
                      N::Int=100,
                      autoinit::Bool=true,
                      solvable::Int=1,
                      dontmargin::Bool=false,
                      labels::Vector{Symbol}=Symbol[],
                      smalldata=Dict{String, String}())::DFGVariable where
                      {G <: AbstractDFG} #
  #
  sto = softtype()
  #TODO: Refactor
  if :ut in fieldnames(typeof(sto))
    sto.ut != -9999999999 ? nothing : error("please define a microsecond time (;ut::Int64=___) for $(softtype)")
  end
  return addVariable!(dfg,
               lbl,
               sto,
               N=N,
               autoinit=autoinit,
               solvable=solvable,
               dontmargin=dontmargin,
               labels=labels,
               smalldata=smalldata  )
end

"""
    $(SIGNATURES)

Fetch the variable marginal sample points without the KDE bandwidth parameter.  Use getVertKDE to retrieve the full KDE object.
"""
function getVal(vA::Vector{DFGVariable}, solveKey::Symbol=:default)::Array{Float64, 2}
  @warn "getVal(::Vector{DFGVariable}) is obsolete, use getVal.(DFGVariable) instead."
  len = length(vA)
  vals = Array{Array{Float64,2},1}()
  cols = Array{Int,1}()
  push!(cols,0)
  rows = Array{Int,1}()
  for v in vA
      push!(vals, getVal(v, solveKey=solveKey))
      c = size(vals[end],2)
      r = size(vals[end],1)
      push!(cols, floor(Int,c))
      push!(rows, floor(Int,r))
  end
  cols = cumsum(cols)
  sc = cols[end]
  rw = floor(Int,rows[1])
  val = Array{Float64,2}(undef,rw, sc)
  for i in 1:(len-1)
      val[:,(cols[i]+1):cols[i+1]] = vals[i]
  end
  val[:,(cols[len]+1):cols[len+1]] = vals[len] # and the last one
  return val
end

"""
    $(SIGNATURES)

Prepare the particle arrays `ARR` to be used for approximate convolution.
This function ensures that ARR has te same dimensions among all the parameters.
Function returns with ARR[sfidx] pointing at newly allocated deepcopy of the
existing values in getVal(Xi[.label==solvefor]).
Return values `sfidx` is the element in ARR where `Xi.label==solvefor` and
`maxlen` is length of all (possibly resampled) `ARR` contained particles.
Note `Xi` is order sensitive.
Note for initialization, solveFor = Nothing.
"""
function prepareparamsarray!(ARR::Array{Array{Float64,2},1},
            Xi::Vector{DFGVariable},
            N::Int,
            solvefor::Union{Nothing, Symbol}  ) # TODO: Confirm we can use symbols here
  #
  LEN = Int[]
  maxlen = N
  count = 0
  sfidx = 0

  for xi in Xi
    push!(ARR, getVal(xi))
    len = size(ARR[end], 2)
    push!(LEN, len)
    if len > maxlen
      maxlen = len
    end
    count += 1
    if xi.label == solvefor
      sfidx = count #xi.index
    end
  end
  SAMP=LEN.<maxlen
  for i in 1:count
    if SAMP[i]
      ARR[i] = KernelDensityEstimate.sample(getKDE(Xi[i]), maxlen)[1]
    end
  end

  # TODO --rather define reusable memory for the proposal
  # we are generating a proposal distribution, not direct replacement for existing memory and hence the deepcopy.
  if sfidx > 0 ARR[sfidx] = deepcopy(ARR[sfidx]) end
  return maxlen, sfidx
end

function parseusermultihypo(multihypo::Nothing)
  verts = Symbol[]
  mh = nothing
  return mh
end
function parseusermultihypo(multihypo::Union{Tuple,Vector{Float64}})
  mh = nothing
  if multihypo != nothing
    multihypo2 = Float64[multihypo...]
    # verts = Symbol.(multihypo[1,:])
    for i in 1:length(multihypo)
      if multihypo[i] > 0.999999
        multihypo2[i] = 0.0
      end
    end
    mh = Categorical(Float64[multihypo2...] )
  end
  return mh
end

# import IncrementalInference: prepgenericconvolution, convert

function prepgenericconvolution(
            Xi::Vector{DFGVariable},
            usrfnc::T;
            multihypo::Union{Nothing, Distributions.Categorical}=nothing,
            threadmodel=MultiThreaded  ) where {T <: FunctorInferenceType}
  #
  ARR = Array{Array{Float64,2},1}()
  maxlen, sfidx = prepareparamsarray!(ARR, Xi, 0, nothing) # Nothing for init.
  fldnms = fieldnames(T) # typeof(usrfnc)
  zdim = T != GenericMarginal ? size(getSample(usrfnc, 2)[1],1) : 0
  certainhypo = multihypo != nothing ? collect(1:length(multihypo.p))[multihypo.p .== 0.0] : collect(1:length(Xi))
  ccw = CommonConvWrapper(
          usrfnc,
          zeros(1,0),
          zdim,
          ARR,
          specialzDim = sum(fldnms .== :zDim) >= 1,
          partial = sum(fldnms .== :partial) >= 1,
          hypotheses=multihypo,
          certainhypo=certainhypo,
          threadmodel=threadmodel
        )
  #
  for i in 1:Threads.nthreads()
    ccw.cpt[i].factormetadata.variableuserdata = []
    ccw.cpt[i].factormetadata.solvefor = :null
    for xi in Xi
      push!(ccw.cpt[i].factormetadata.variableuserdata, solverData(xi).softtype)
    end
  end
  return ccw
end

function setDefaultFactorNode!(
      dfg::G,
      factor::DFGFactor,
      Xi::Vector{DFGVariable},
      usrfnc::T;
      multihypo::Union{Nothing,Tuple,Vector{Float64}}=nothing,
      threadmodel=SingleThreaded  )::GenericFunctionNodeData where
      {G <: AbstractDFG, T <: Union{FunctorInferenceType, InferenceType}}
  #
  ftyp = T
  mhcat = parseusermultihypo(multihypo)

  ccw = prepgenericconvolution(Xi, usrfnc, multihypo=mhcat, threadmodel=threadmodel)

  data_ccw = FunctionNodeData{CommonConvWrapper{T}}(Int[], false, false, Int[], Symbol(ftyp.name.module), ccw)
  factor.data = data_ccw

  return factor.data
end

"""
    $SIGNATURES

Returns state of vertex data `.initialized` flag.

Notes:
- used by Bayes tree clique logic.
- similar method in DFG
"""
function isInitialized(vert::TreeClique)::Bool
  return solverData(vert).initialized
end

"""
    $SIGNATURES

Return `::Bool` on whether at least one hypothesis is available for intended computations (assuming direction `sfidx`).
"""
function isLeastOneHypoAvailable(sfidx::Int,
                                 certainidx::Vector{Int},
                                 uncertnidx::Vector{Int},
                                 isinit::Vector{Bool})::Bool
  #
  # @show isinit
  # @show sfidx in certainidx, sum(isinit[uncertnidx])
  # @show sfidx in uncertnidx, sum(isinit[certainidx])
  return sfidx in certainidx && 0 < sum(isinit[uncertnidx]) ||
         sfidx in uncertnidx && sum(isinit[certainidx]) == length(certainidx)
end

"""
    $SIGNATURES

Return `(::Bool, ::OKVarlist, ::NotOkayVarList)` on whether all other variables (besides `loovar::Symbol`)
attached to factor `fct::Symbol` are all initialized -- i.e. `fct` is usable.

Notes:
- Special carve out for multihypo cases, see issue 427, where at least one hypothesis should be available, but not all required at first.

Development Notes
* TODO get faster version of isInitialized for database version

Related

doautoinit!, manualinit!, isInitialized, isMultihypo
"""
function factorCanInitFromOtherVars(dfg::AbstractDFG,
                                    fct::Symbol,
                                    loovar::Symbol)::Tuple{Bool, Vector{Symbol}, Vector{Symbol}}
  #
  # all variables attached to this factor
  varsyms = DFG.getNeighbors(dfg, fct)

  # which element is being solved for
  sfidx = (1:length(varsyms))[varsyms .== loovar][1]
  # list of factors to use in init operation
  fctlist = []
  # list fo variables that cannot be used
  faillist = Symbol[]
  isinit = Bool[]
  for vsym in varsyms
    # check each variable one by one
    xi = DFG.getVariable(dfg, vsym)
    isi = isInitialized(xi)
    push!(isinit, isi)
    if !isi
      push!(faillist, vsym)
    end
  end

  ## determine if this factor can be used
  # priors and general n-ary cases
  canuse = length(varsyms)==1 || (length(faillist)==1 && loovar in faillist)
  ## special multihypo case (at least one hypothesis is available or initializing first hypo)
  fctnode = getFactor(dfg, fct)
  # @show canuse, isMultihypo(fctnode), isinit
  if !canuse && isMultihypo(fctnode)
    # multihypo=[1;0.5;0.5] : sfidx=1, isinit=[0,1,0] -- true
    # multihypo=[1;0.5;0.5] : sfidx=1, isinit=[0,0,1] -- true
    # multihypo=[1;0.5;0.5] : sfidx=2|3, isinit=[1,0,0] -- true
    mhp = getMultihypoDistribution(fctnode).p
    allmhp,certainidx,uncertnidx = getHypothesesVectors(mhp)
    if isLeastOneHypoAvailable(sfidx, certainidx, uncertnidx, isinit)
       # special case works
       @info "allowing init from incomplete set of previously initialized hypotheses, fct=$fct"
       canuse = true
    end
  end

  # should add the factor for use?
  if canuse
    push!(fctlist, fct)
  end

  # return if can use, the factor in an array, and the non-initialized variables attached to the factor
  return (canuse, fctlist, faillist )
end

# wow, that was quite far off -- needs testing
# function factorCanInitFromOtherVars(dfg::T,
#                                     fct::Symbol,
#                                     loovar::Symbol)::Tuple{Bool, Vector{Symbol}, Vector{Symbol}} where T <: AbstractDFG
#   #
#   # all variables attached to this factor
#   varsyms = getNeighbors(dfg, fct)
#
#   # list of factors to use in init operation
#   useinitfct = Symbol[]
#   faillist = Symbol[]
#   for vsym in varsyms
#     xi = DFG.getVariable(dfg, vsym)
#     if (isInitialized(xi) && sum(useinitfct .== fct) == 0 ) || length(varsyms) == 1
#       push!(useinitfct, fct)
#     end
#   end
#
#   return (length(useinitfct)==length(varsyms)&&length(faillist)==0,
#           useinitfct,
#           faillist   )
# end

"""
    $(SIGNATURES)

EXPERIMENTAL: initialize target variable `xi` based on connected factors in the
factor graph `fgl`.  Possibly called from `addFactor!`, or `doCliqAutoInitUp!` (?).

Notes:
- Special carve out for multihypo cases, see issue 427.

Development Notes:
> Target factor is first (singletons) or second (dim 2 pairwise) variable vertex in `xi`.
* TODO use DFG properly with local operations and DB update at end.
* TODO get faster version of `isInitialized` for database version.
* TODO: Persist this back if we want to here.
"""
function doautoinit!(dfg::T,
                     xi::DFGVariable;
                     singles::Bool=true,
                     N::Int=100,
                     logger=ConsoleLogger() )::Bool where T <: AbstractDFG
  #
  didinit = false
  # don't initialize a variable more than once
  if !isInitialized(xi)
    with_logger(logger) do
      @info "try doautoinit! of $(xi.label)"
    end
    # get factors attached to this variable xi
    vsym = xi.label
    neinodes = DFG.getNeighbors(dfg, vsym)
    # proceed if has more than one neighbor OR even if single factor
    if (singles || length(neinodes) > 1)
      # Which of the factors can be used for initialization
      useinitfct = Symbol[]
      # Consider factors connected to $vsym...
      for xifct in neinodes
        canuse, usefct, notusevars = factorCanInitFromOtherVars(dfg, xifct, vsym)
        if canuse
          union!(useinitfct, usefct)
        end
      end
      with_logger(logger) do
        @info "init with useinitfct $useinitfct"
      end
      # println("Consider all singleton (unary) factors to $vsym...")
      # calculate the predicted belief over $vsym
      if length(useinitfct) > 0
        with_logger(logger) do
          @info "do init of $vsym"
        end
        pts,inferdim = predictbelief(dfg, vsym, useinitfct, logger=logger)
        setValKDE!(xi, pts, true, inferdim)
        #TODO test
        setVariablePosteriorEstimates!(xi)
        didinit = true
      end
    end
  end
  return didinit
end

function doautoinit!(dfg::T,
                     Xi::Vector{DFGVariable};
                     singles::Bool=true,
                     N::Int=100,
                     logger=ConsoleLogger() )::Bool where T <: AbstractDFG
  #
  #
  # Mighty inefficient function, since we only need very select fields nearby from a few neighboring nodes
  # do double depth search for variable nodes

  didinit = true

  # loop over all requested variables that must be initialized
  for xi in Xi
    didinit &= doautoinit!(dfg, xi, singles=singles, N=N, logger=logger)
  end
  return didinit
end

function doautoinit!(dfg::T,
                     xsyms::Vector{Symbol};
                     singles::Bool=true,
                     N::Int=100,
                     logger=SimpleLogger(logger)  )::Bool where T <: AbstractDFG
  #
  verts = getVariable.(dfg, xsyms)
  return doautoinit!(dfg, verts, singles=singles, N=N, logger=logger)
end
function doautoinit!(dfg::T,
                     xsym::Symbol;
                     singles::Bool=true,
                     N::Int=100,
                     logger=ConsoleLogger()  )::Bool where T <: AbstractDFG
  #
  return doautoinit!(dfg, [getVariable(dfg, xsym);], singles=singles, N=N, logger=logger)
end

"""
    $(SIGNATURES)

Workaround function when first-version (factor graph based) auto initialization fails.  Usually occurs when using factors that have high connectivity to multiple variables.
"""
function manualinit!(dfg::AbstractDFG, vert::DFGVariable, pX::BallTreeDensity)::Nothing
  setValKDE!(vert, pX, true)
  # getData(vert).initialized = true
  return nothing
end
function manualinit!(dfg::AbstractDFG, sym::Symbol, pX::BallTreeDensity)::Nothing
  vert = getVariable(dfg, sym)
  manualinit!(dfg, vert, pX)
  return nothing
end
function manualinit!(dfg::AbstractDFG, sym::Symbol, usefcts::Vector{Symbol})::Nothing
  @info "manualinit! $sym"
  pts = predictbelief(dfg, sym, usefcts)
  vert = getVariable(dfg, sym)
  Xpre = AMP.manikde!(pts, getSofttype(vert).manifolds )
  setValKDE!(vert, Xpre, true) # dfg, sym
  # getData(dfg, sym).initialized = true
  return nothing
end

function manualinit!(dfg::AbstractDFG, sym::Symbol, pts::Array{Float64,2})
  var = getVariable(dfg, sym)
  pp = manikde!(pts, getManifolds(var))
  manualinit!(dfg,sym,pp)
end

function ensureAllInitialized!(dfg::T; solvable::Int=1) where T <: AbstractDFG
  # allvarnodes = getVariables(dfg)
  syms = intersect(getAddHistory(dfg), ls(dfg, solvable=solvable) )
  # syms = ls(dfg, solvable=solvable) # |> sortDFG
  repeatCount = 0
  repeatFlag = true
  while repeatFlag
    repeatFlag = false
    repeatCount += 1
    if 10 < repeatCount
      @info "not able to initialize all variables via the factor graph, abort autoinit."
      break;
    end
    for sym in syms
      var = getVariable(dfg, sym)
      if !isInitialized(var)
        @info "$(var.label) is not initialized, and will do so now..."
        doautoinit!(dfg, [var;], singles=true)
        !isInitialized(var) ? (repeatFlag = true) : nothing
      end
    end
  end
  nothing
end

function assembleFactorName(dfg::T, Xi::Vector{DFGVariable}; maxparallel::Int=50)::Symbol where T <: AbstractDFG
  existingFactorLabels = getFactorIds(dfg)
  existingFactorLabelDict = Dict(existingFactorLabels .=> existingFactorLabels)
  namestring = ""
  for vert in Xi #f.Xi
    namestring = string(namestring,vert.label)
  end
  for i in 1:maxparallel
    tempnm = string(namestring, "f$i")
    if !haskey(existingFactorLabelDict, Symbol(tempnm))
      namestring = tempnm
      break
    end
    i != maxparallel ? nothing : error("Cannot currently add more than $(maxparallel) factors in parallel, please open an issue if this is too restrictive.")
  end
  return Symbol(namestring)
end

"""
    $(SIGNATURES)

Add factor with user defined type <: FunctorInferenceType to the factor graph
object. Define whether the automatic initialization of variables should be
performed.  Use order sensitive `multihypo` keyword argument to define if any
variables are related to data association uncertainty.
"""
function addFactor!(dfg::G,
                    Xi::Vector{DFGVariable},
                    usrfnc::R;
                    multihypo::Union{Nothing,Tuple,Vector{Float64}}=nothing,
                    solvable::Int=1,
                    labels::Vector{Symbol}=Symbol[],
                    autoinit::Bool=true,
                    threadmodel=SingleThreaded,
                    maxparallel::Int=50  ) where
                      {G <: AbstractDFG,
                       R <: Union{FunctorInferenceType, InferenceType}}
  #
  namestring = assembleFactorName(dfg, Xi, maxparallel=maxparallel)
  newFactor = DFGFactor{CommonConvWrapper{R}, Symbol}(Symbol(namestring))
  newFactor.tags = union(labels, [:FACTOR]) # TODO: And session info
  newFactor.solvable = solvable
  # addNewFncVertInGraph!(fgl, newvert, currid, namestring, solvable)
  newData = setDefaultFactorNode!(dfg, newFactor, Xi, deepcopy(usrfnc), multihypo=multihypo, threadmodel=threadmodel)

  # TODO: Need to remove this...
  for vert in Xi
    push!(newData.fncargvID, vert.label) # vert._internalId # YUCK :/ -- Yup, this is a problem
  end

  success = DFG.addFactor!(dfg, Xi, newFactor)

  # TODO: change this operation to update a conditioning variable
  autoinit && doautoinit!(dfg, Xi, singles=false)

  return newFactor
end

function addFactor!(dfg::AbstractDFG,
                    xisyms::Vector{Symbol},
                    usrfnc::Union{FunctorInferenceType, InferenceType};
                    multihypo::Union{Nothing,Tuple,Vector{Float64}}=nothing,
                    solvable::Int=1,
                    labels::Vector{Symbol}=Symbol[],
                    autoinit::Bool=true,
                    threadmodel=SingleThreaded,
                    maxparallel::Int=50  )
  #
  verts = map(vid -> DFG.getVariable(dfg, vid), xisyms)
  addFactor!(dfg, verts, usrfnc, multihypo=multihypo, solvable=solvable, labels=labels, autoinit=autoinit, threadmodel=threadmodel, maxparallel=maxparallel )
end

# """
#     $SIGNATURES
#
# Delete factor and its edges.
# """
# function deleteFactor!(fgl::FactorGraph, fsym::Symbol)
#   fid = fgl.fIDs[fsym]
#   eds = fgl.g.inclist[fid]
#   alledsids = Int[]
#   nedges = length(eds)
#   for eds in fgl.g.inclist[fid]
#     union!(alledsids, [eds.source.index; eds.target.index])
#   end
#   for edids in setdiff!(alledsids, fid)
#     count = 0
#     for eds in fgl.g.inclist[edids]
#       count += 1
#       if fid == eds.source.index || fid == eds.target.index
#         deleteat!(fgl.g.inclist[edids], count)
#         break
#       end
#     end
#   end
#   delete!(fgl.g.inclist, fid)
#   fgl.g.nedges -= nedges
#   delete!(fgl.g.vertices, fid)
#   delete!(fgl.fIDs, fsym)
#   deleteat!(fgl.factorIDs, findfirst(a -> a==fid, fgl.factorIDs))
#   nothing
# end

# """
#     $SIGNATURES
#
# Delete variables, and also the factors+edges if `andfactors=true` (default).
# """
# function deleteVariable!(fgl::FactorGraph, vsym::Symbol; andfactors::Bool=true)
#   vid = fgl.IDs[vsym]
#   vert = fgl.g.vertices[vid]
#   if andfactors
#     for ne in Graphs.out_neighbors(vert, fgl.g)
#       deleteFactor!(fgl, Symbol(ne.label))
#     end
#   end
#   delete!(fgl.g.vertices, vid)
#   delete!(fgl.IDs, vsym)
#   deleteat!(fgl.nodeIDs, findfirst(a -> a==vid, fgl.nodeIDs))
#   nothing
# end

function prtslperr(s)
  println(s)
  sleep(0.1)
  error(s)
end

"""
    $SIGNATURES

Determine the variable ordering used to construct both the Bayes Net and Bayes/Junction/Elimination tree.

Notes
- Heuristic method -- equivalent to QR or Cholesky.
- Are using Blas `QR` function to extract variable ordering.
- **NOT USING SUITE SPARSE** -- which would requires commercial license.
- For now `A::Array{<:Number,2}` as a dense matrix.
- Columns of `A` are system variables, rows are factors (without differentiating between partial or full factor).
- default is to use `solvable=1` and ignore factors and variables that might be used for dead reckoning or similar.

Future
- TODO: `A` should be sparse data structure (when we exceed 10'000 var dims)
- TODO: Incidence matrix is rectagular and adjacency is the square.
"""
function getEliminationOrder(dfg::G; ordering::Symbol=:qr, solvable::Int=1) where G <: AbstractDFG
  # Get the sparse adjacency matrix, variable, and factor labels
  adjMat, permuteds, permutedsf = DFG.getAdjacencyMatrixSparse(dfg, solvable=solvable)

  # Create dense adjacency matrix
  A = Array(adjMat)

  p = Int[]
  if ordering==:chol
      p = cholfact(A'A,:U,Val(true))[:p] #,pivot=true
  elseif ordering==:qr
    # this is the default
    q,r,p = qr(A, Val(true))
  else
    prtslperr("getEliminationOrder -- cannot do the requested ordering $(ordering)")
  end

  # Return the variable ordering that we should use for the Bayes map
  # reverse order checked in #475 and #499
  return permuteds[p] |> reverse
end

# lets create all the vertices first and then deal with the elimination variables thereafter
function addBayesNetVerts!(dfg::AbstractDFG, elimOrder::Array{Symbol,1})
  for pId in elimOrder
    vert = DFG.getVariable(dfg, pId)
    if solverData(vert).BayesNetVertID == nothing || solverData(vert).BayesNetVertID == :_null # Special serialization case of nothing
      @debug "[AddBayesNetVerts] Assigning $pId.data.BayesNetVertID = $pId"
      solverData(vert).BayesNetVertID = pId
    else
      @warn "addBayesNetVerts -- Something is wrong, variable '$pId' should not have an existing Bayes net reference to '$(solverData(vert).BayesNetVertID)'"
    end
  end
end

function addConditional!(dfg::AbstractDFG, vertId::Symbol, Si::Vector{Symbol})::Nothing
  bnv = DFG.getVariable(dfg, vertId)
  bnvd = solverData(bnv)
  bnvd.separator = Si
  for s in Si
    push!(bnvd.BayesNetOutVertIDs, s)
  end
  return nothing
end

function addChainRuleMarginal!(dfg::AbstractDFG, Si::Vector{Symbol}; maxparallel::Int=50)
  #
    @show Si
  lbls = String[]
  genmarg = GenericMarginal()
  Xi = map(v -> DFG.getVariable(dfg, v), Si)
  # @info "adding marginal to"
  # for x in Xi
  #   @info "x.index=",x.index
  # end
  addFactor!(dfg, Xi, genmarg, autoinit=false, maxparallel=maxparallel)
  nothing
end

function rmVarFromMarg(dfg::G,
                       fromvert::DFGVariable,
                       gm::Vector{DFGFactor};
                       maxparallel::Int=50 )::Nothing where G <: AbstractDFG
  @info " - Removing $(fromvert.label)"
  for m in gm
    @info "Looking at $(m.label)"
    for n in DFG.getNeighbors(dfg, m) #x1, x2
      if n == fromvert.label # n.label ==? x1
        @info "   - Breaking link $(m.label)->$(fromvert.label)..."
        @info "     - Original links: $(DFG.ls(dfg, m))"
        remvars = setdiff(DFG.ls(dfg, m), [fromvert.label])
        @info "     - New links: $remvars"

        DFG.deleteFactor!(dfg, m) # Remove it
        if length(remvars) > 0
          @info "$(m.label) still has links to other variables, readding it back..."
          addFactor!(dfg, remvars, solverData(m).fnc.usrfnc!, autoinit=false, maxparallel=maxparallel)
        else
          @info "$(m.label) doesn't have any other links, not adding it back..."
        end
      end
    end
    # Added back in chain rule.
    if DFG.exists(dfg, m) && length(DFG.getNeighbors(dfg, m)) <= 1
      @warn "removing vertex id=$(m.label)"
      DFG.deleteFactor!(dfg, m)
    end
  end
  return nothing
end

function buildBayesNet!(dfg::G,
                        elimorder::Vector{Symbol};
                        maxparallel::Int=50,
                        solvable::Int=1)::Nothing where G <: AbstractDFG
  #
  # addBayesNetVerts!(dfg, elimorder)
  for v in elimorder
    @info ""
    @info "Eliminating $(v)"
    @info "==============="
    @info ""
    # which variable are we eliminating

    # all factors adjacent to this variable
    fi = Symbol[]
    Si = Symbol[]
    gm = DFGFactor[]

    vert = DFG.getVariable(dfg, v)
    for fctId in DFG.getNeighbors(dfg, vert, solvable=solvable)
      fct = DFG.getFactor(dfg, fctId)
      if (solverData(fct).eliminated != true)
        push!(fi, fctId)
        for sepNode in DFG.getNeighbors(dfg, fct, solvable=solvable)
          # TODO -- validate !(sepNode.index in Si) vs. older !(sepNode in Si)
          if sepNode != v && !(sepNode in Si) # Symbol comparison!
            push!(Si,sepNode)
          end
        end
        solverData(fct).eliminated = true
      end

      if typeof(solverData(fct).fnc) == CommonConvWrapper{GenericMarginal}
        push!(gm, fct)
      end
    end

    if v != elimorder[end]
      addConditional!(dfg, v, Si)
      # not yet inserting the new prior p(Si) back into the factor graph
    end

    # mark variable
    solverData(vert).eliminated = true

    # TODO -- remove links from current vertex to any marginals
    rmVarFromMarg(dfg, vert, gm, maxparallel=maxparallel)

    #add marginal on remaining variables... ? f(xyz) = f(x | yz) f(yz)
    # new function between all Si (round the outside, right the outside)
    length(Si) > 0 && addChainRuleMarginal!(dfg, Si, maxparallel=maxparallel)

  end
  return nothing
end

function getKDE(vnd::VariableNodeData)
  AMP.manikde!(getVal(vnd), getBW(vnd)[:,1], getSofttype(vnd).manifolds)
end

"""
    $(SIGNATURES)

Get KernelDensityEstimate kde estimate stored in variable node.
"""
function getKDE(v::DFGVariable)
  return getKDE(solverData(v))
end

function getVert(dfg::G, sym::Symbol, nt::Symbol=:var) where G <: AbstractDFG
  @warn "IIF.getVert is deprecated, use DFG.getVariable or DFG.getFactor instead."
  if nt == :var
    return DFG.getVariable(dfg, sym)
  elseif nt == :fct
    return DFG.getFactor(dfg, sym)
  else
    error("unknown getVert request nt=$nt")
  end
end

"""
    $(SIGNATURES)

Get KernelDensityEstimate kde estimate stored in variable node.
"""
function getVertKDE(v::DFGVariable)
  return getKDE(v)
end
function getVertKDE(dfg::G, id::Int) where G <: AbstractDFG
  v = DFG.getVariable(dfg, id)
  return getKDE(v)
end
function getVertKDE(dfg::G, lbl::Symbol) where G <: AbstractDFG
  v = DFG.getVariable(dfg, lbl)
  return getKDE(v)
end
function getKDE(dfg::G, lbl::Symbol) where G <: AbstractDFG
  return getVertKDE(dfg, lbl)
end

function expandEdgeListNeigh!(fgl::FactorGraph,
                              vertdict::Dict{Int,TreeClique},
                              edgedict::Dict{Int,Graphs.Edge{TreeClique}})
  #asfd
  for vert in vertdict
    for newedge in out_edges(vert[2],fgl.g)
      if !haskey(edgedict, newedge.index)
        edgedict[newedge.index] = newedge
      end
    end
  end

  nothing
end

# dictionary of unique vertices from edgelist
function expandVertexList!(fgl::FactorGraph,
  edgedict::Dict{Int,Graphs.Edge{TreeClique}},
  vertdict::Dict{Int,TreeClique})

  # go through all source and target nodes
  for edge in edgedict
    if !haskey(vertdict, edge[2].source.index)
      vertdict[edge[2].source.index] = edge[2].source
    end
    if !haskey(vertdict, edge[2].target.index)
      vertdict[edge[2].target.index] = edge[2].target
    end
  end
  nothing
end

function edgelist2edgedict(edgelist::Array{Graphs.Edge{TreeClique},1})
  edgedict = Dict{Int,Graphs.Edge{TreeClique}}()
  for edge in edgelist
    edgedict[edge.index] = edge
  end
  return edgedict
end

#


text:
using SPADL

get_events_data(:wyscout)

# Teams


text:
# run this script in the `docs` project to fix the doctests if they are out of date
# carefully review the changes before committing!

using Documenter, KeywordSearch

DocMeta.setdocmeta!(KeywordSearch, :DocTestSetup, :(using KeywordSearch); recursive=true)

if success(`git diff --quiet`)
    doctest(KeywordSearch; fix=true)
else
    error("Git repo dirty; commit changes before fixing doctests.")
end

max_stars_repo_name: beacon-biosignals/KeywordSearch

max_issues_repo_name: beacon-biosignals/KeywordSearch

max_forks_repo_name: beacon-biosignals/KeywordSearch


text:
"""
    function upgrade_branches(
       base_network::String,
       output_location::String,
       fbase;
       indices=[],
       output_type="input"
       )
Creates a folder of network data for the network `base_network` with one line converted to LFAC, one for each index in `indices`, or if `indices` is empty, for every non-transformer branch in the network.
"""
function upgrade_branches(
      base_network::String,
      output_location::String,
      fbase;
      subnet_params=Dict(),
      interface_params=Dict(),
      indices=[],
      output_type="input",
      noloss=false
   )

   PowerModels.logger_config!("error")

   filename = split(base_network, "/")[end]
   println("$filename")
   filetype = split(lowercase(filename), '.')[end]
   base_name = filename[1:(end-length(filetype)-1)]
   (base_name, filetype) = splitext(filename)
   filetype = Unicode.normalize(filetype, casefold=true)
   fileext = filetype == ".m" ? filetype : ".json"
   if lowercase(output_type) == ".json"
      fileext = ".json"
   end
   # println(fileext)
   base_dir = base_network[1:end-length(filename)]

   # println("base_network: $base_network")
   # println("base_dir: $base_dir")

   base_subnet_df = CSV.read("$base_dir/subnetworks.csv", DataFrame, copycols=true)

   # Subnet parameters
   if !("var_f" in keys(subnet_params))
      subnet_params["var_f"] = false
   end
   if !("fmin" in keys(subnet_params))
      subnet_params["fmin"] = fbase
   end
   if !("fmax" in keys(subnet_params))
      subnet_params["fmax"] = fbase
   end
   if !("var_f_lfac" in keys(subnet_params))
      subnet_params["var_f_lfac"] = true
   end
   if !("fbase_lfac" in keys(subnet_params))
      subnet_params["fbase_lfac"] = fbase
   end
   if !("fmin_lfac" in keys(subnet_params))
      subnet_params["fmin_lfac"] = 0.1
   end
   if !("fmax_lfac" in keys(subnet_params))
      subnet_params["fmax_lfac"] = fbase
   end
   if "contingency_file" in names(base_subnet_df)
      ctg_filename = base_subnet_df.contingency_file[1]
      subnet_df = DataFrame(
         index=[1,2],
         file=["$(base_name)$(fileext)", "lfac_$(base_name)$(fileext)"],
         variable_f=[subnet_params["var_f"], subnet_params["var_f_lfac"]],
         f_base=[fbase, subnet_params["fbase_lfac"]],
         f_min=[subnet_params["fmin"], subnet_params["fmin_lfac"]],
         f_max=[subnet_params["fmax"], subnet_params["fmax_lfac"]],
         contingency_file=[ctg_filename, "lfac_$(ctg_filename)"],
         interfaces_file="interfaces_file" in names(base_subnet_df) ? [base_subnet_df.interfaces_file[1], base_subnet_df.interfaces_file[1]] : ["interfaces.csv","interfaces.csv"]
      )
   else
      subnet_df = DataFrame(
         index=[1,2],
         file=["$(base_name)$(fileext)", "lfac_$(base_name)$(fileext)"],
         variable_f=[subnet_params["var_f"], subnet_params["var_f_lfac"]],
         f_base=[fbase, subnet_params["fbase_lfac"]],
         f_min=[subnet_params["fmin"], subnet_params["fmin_lfac"]],
         f_max=[subnet_params["fmax"], subnet_params["fmax_lfac"]],
         interfaces_file="interfaces_file" in names(base_subnet_df) ? [base_subnet_df.interfaces_file[1], base_subnet_df.interfaces_file[1]] : ["interfaces.csv","interfaces.csv"]
      )
   end

   data = PowerModels.parse_file(base_network)

   ctg_br_dict = Dict{Int64, Any}()

   if "contingency_file" in names(subnet_df)
      ctg_list = parse_con_file(base_dir*"/"*subnet_df.contingency_file[1])

      for (ctg_idx,ctg) in enumerate(ctg_list)
         if ctg["component"] == "branch"
            circuit = parse(Int64, ctg["ckt"])

            br_idx = parse(Int64, sort(
               collect(keys(
                  filter(p->(
                     (last(p)["f_bus"]==ctg["i"] && last(p)["t_bus"]==ctg["j"])
                     || (last(p)["f_bus"]==ctg["j"] && last(p)["t_bus"]==ctg["i"])),
                     data["branch"]
                     )
                  )
               )
            )[circuit])

            # println("br_idx: $(br_idx)")
            ctg_br_dict[br_idx] = ctg

         end
      end
   end

   # Converter limits and loss parameters in p.u.
   if !("imax" in keys(interface_params))
      interface_params["imax"] = 40.0
   end
   if !("vmax" in keys(interface_params))
      interface_params["vmax"] = 1.1
   end
   if !("c1" in keys(interface_params))
      interface_params["c1"] = noloss ? 0 : 1e-4
   end
   if !("c2" in keys(interface_params))
      interface_params["c2"] = noloss ? 0 : 1e-4
   end
   if !("c3" in keys(interface_params))
      interface_params["c3"] = noloss ? 0 : 1e-4
   end
   if !("sw1" in keys(interface_params))
      interface_params["sw1"] = noloss ? 0 : 1e-4
   end
   if !("sw2" in keys(interface_params))
      interface_params["sw2"] = noloss ? 0 : 1e-4
   end
   if !("sw3" in keys(interface_params))
      interface_params["sw3"] = noloss ? 0 : 1e-4
   end
   if !("M" in keys(interface_params))
      interface_params["M"] = 0.9
   end
   if !("smax" in keys(interface_params))
      interface_params["smax"] = 50.0
   end
   if !("R" in keys(interface_params))
      interface_params["R"] = noloss ? 0 : 1.0e-5
   end
   if !("X" in keys(interface_params))
      interface_params["X"] = noloss ? 0 : 1.0e-3
   end
   if !("G" in keys(interface_params))
      interface_params["G"] = 0.0
   end
   if !("B" in keys(interface_params))
      interface_params["B"] = noloss ? 0 : 1.0e-4
   end
   if !("transformer" in keys(interface_params))
      interface_params["transformer"] = false
   end
   if !("tap" in keys(interface_params))
      interface_params["tap"] = 1.0
   end
   if !("shift" in keys(interface_params))
      interface_params["shift"] = 0.0
   end
   s_max = 100
   p_max = 0.0
   q_max = 0.0
   loss_param = 0.0
   lfac_vmin = 0.9
   lfac_vmax = 1.1
   # println(indices)
   for (i,(br_idx,branch)) in enumerate(data["branch"])
      lfac_ctg_br_dict = Dict{Int64, Any}()
      ctg_br_dict_tmp = deepcopy(ctg_br_dict)
      br_specified = ((length(indices) == 0) || (parse(Int64, br_idx) in indices)) # Defines whether or not the branch was specified as an input. Include all if empty set was specified.
      if (branch["br_status"]==1) && (!branch["transformer"]) && br_specified
         println("LFAC upgrade on branch $(branch["index"])")

         # iface_keys = ["index", "subnet_index", "bus"]
         # iface_types = [Int64, Int64, Int64]
         # # Add all other parameters in the interface_params dictionary
         # for (key,val) in keys(interface_params)
         #    push!(iface_keys, key)
         #    push!(iface_types, typeof(val))
         # end
         iface_df = DataFrame(index=Int64[],
            subnet_index=Int64[],
            bus=Int64[],
            imax=Float64[],
            vmax=Float64[],
            c1=Float64[],
            c2=Float64[],
            c3=Float64[],
            sw1=Float64[],
            sw2=Float64[],
            sw3=Float64[],
            M=Float64[],
            smax=Float64[],
            R=Float64[],
            X=Float64[],
            G=Float64[],
            B=Float64[],
            transformer=Bool[],
            tap=Float64[],
            shift=Float64[]
            )
         # println("sourceid: $(branch["source_id"])")
         if (length(branch["source_id"])>3) && (Unicode.normalize(branch["source_id"][4], casefold=true) == "sc") # is a series compensation branch
            seriescomp_bus = branch["sc_bus"]
            # println("seriescomp_bus: $seriescomp_bus")
            # Set two interface buses for base network
            subnet_index = 1
            base_f_bus = nothing
            base_t_bus = nothing
            next_f_bus = nothing
            next_t_bus = nothing
            lfac_next_branch = nothing
            if branch["f_bus"] == seriescomp_bus
               sc_f_bus = seriescomp_bus
               sc_t_bus = 2 # series compensation branch goes from seriescomp_bus to the t_bus in the LFAC network (always bus 2)
               base_t_bus = copy(branch["t_bus"])
               for (next_br_idx,next_branch) in data["branch"]
                  if next_branch["index"] != branch["index"] # exclude series compensation branch
                     if next_branch["f_bus"] == seriescomp_bus
                        base_f_bus = copy(next_branch["t_bus"])
                        lfac_next_branch = next_branch
                        next_f_bus = seriescomp_bus
                        next_t_bus = 1 # existing branch goes from seriescomp_bus to the f_bus in the LFAC network (bus 1)
                        println("Series compensation. Including connected branch $(lfac_next_branch["index"]) in LFAC network.")
                        break
                     elseif next_branch["t_bus"] == seriescomp_bus
                        base_f_bus = copy(next_branch["f_bus"])
                        lfac_next_branch = next_branch
                        next_t_bus = seriescomp_bus
                        next_f_bus = 1 # existing branch goes to seriescomp_bus from the f_bus in the LFAC network (bus 1)
                        println("Series compensation. Including connected branch $(lfac_next_branch["index"]) in LFAC network.")
                        break
                     end
                  end
               end
            elseif branch["t_bus"] == seriescomp_bus
               sc_t_bus = seriescomp_bus
               sc_f_bus = 1 # series compensation branch goes to seriescomp_bus from the f_bus in the LFAC network (always bus 1)
               base_f_bus = copy(branch["f_bus"])
               for (next_br_idx,next_branch) in data["branch"]
                  if next_branch["index"] != branch["index"] # exclude series compensation branch
                     if next_branch["f_bus"] == seriescomp_bus
                        base_t_bus = copy(next_branch["t_bus"])
                        lfac_next_branch = next_branch
                        next_f_bus = seriescomp_bus
                        next_t_bus = 2 # existing branch goes from seriescomp_bus to the t_bus in the LFAC network (bus 2)
                        println("Series compensation. Including connected branch $(lfac_next_branch["index"]) in LFAC network.")
                        break
                     elseif next_branch["t_bus"] == seriescomp_bus
                        base_t_bus = copy(next_branch["f_bus"])
                        lfac_next_branch = next_branch
                        next_t_bus = seriescomp_bus
                        next_f_bus = 2 # existing branch goes to seriescomp_bus from the t_bus in the LFAC network (bus 2)
                        println("Series compensation. Including connected branch $(lfac_next_branch["index"]) in LFAC network.")
                        break
                     end
                  end
               end
            else
               throw(ArgumentError("The series compensation bus $seriescomp_bus specified for branch $(branch["index"]) is not connected to the branch. f_bus=$f_bus and t_bus=$t_bus"))
            end

            lfac_buslookup = Dict(
               base_f_bus=>1,
               base_t_bus=>2,
               seriescomp_bus=>seriescomp_bus
            )

            # create output directory
            output_dir = "$output_location/br$(br_idx)_$(lfac_next_branch["index"])"
            if !isdir(output_dir)
               mkpath(output_dir)
            end

            # LFAC "from bus" interface
            iface_index = 1
            bus = base_f_bus
            interface_params["index"] = iface_index
            interface_params["subnet_index"] = subnet_index
            interface_params["bus"] = bus
            push!(iface_df, interface_params)
            f_area = copy(data["bus"]["$bus"]["area"])
            f_zone = copy(data["bus"]["$bus"]["zone"])
            f_kv_base = copy(data["bus"]["$bus"]["base_kv"])

            # LFAC "to bus" interface
            iface_index = 2
            bus = base_t_bus
            interface_params["index"] = iface_index
            interface_params["subnet_index"] = subnet_index
            interface_params["bus"] = bus
            push!(iface_df, interface_params)
            t_area = copy(data["bus"]["$bus"]["area"])
            t_zone = copy(data["bus"]["$bus"]["zone"])
            t_kv_base = copy(data["bus"]["$bus"]["base_kv"])

            for (bus_idx,bus) in data["bus"]
               bus["name"] = split(bus["name"], ";")[1]
            end
            println("series comp branch: $(branch["index"])")
            println("lfac next branch: $(lfac_next_branch["index"])")
            println("saving network at $output_dir/$(base_name)$(fileext)")
            # Disable branches and save network
            # Disabling branches is the only change to the base network.
            branch["br_status"] = 0
            lfac_next_branch["br_status"] = 0
            if fileext == ".m"
               io = open("$output_dir/$(base_name).m", "w");
               export_matpower(io, data)
               close(io)
            else
               stringnet = JSON.json(data)
               open("$output_dir/$(base_name)$(fileext)", "w") do f
                  write(f, stringnet)
               end
            end
            # Restore branch status
            branch["br_status"] = 1
            lfac_next_branch["br_status"] = 1

            # Create and save LFAC subnetwork
            lfac_case = Dict{String,Any}()
            lfac_case["name"] = data["name"]
            lfac_case["source_version"] = data["source_version"]
            lfac_case["baseMVA"] = data["baseMVA"]
            lfac_case["per_unit"] = true
            # Bus data
            lfac_case["bus"] = Dict{String, Any}(
               "1"=>Dict{String, Any}(
                  "bus_type"=>3,
                  "vmin"=>lfac_vmin,
                  "vmax"=>lfac_vmax,
                  "bus_i"=>1,
                  "index"=>1,
                  "va"=>0.0,
                  "vm"=>1.0,
                  "name"=>"Bus 1\tLF",
                  "area"=>f_area,
                  "zone"=>f_zone,
                  "base_kv"=>f_kv_base
               ),
               "2"=>Dict{String, Any}(
                  "bus_type"=>1,
                  "vmin"=>lfac_vmin,
                  "vmax"=>lfac_vmax,
                  "bus_i"=>2,
                  "index"=>2,
                  "va"=>0.0,
                  "vm"=>1.0,
                  "name"=>"Bus 2\tLF",
                  "area"=>t_area,
                  "zone"=>t_zone,
                  "base_kv"=>t_kv_base
               ),
               "$seriescomp_bus"=>Dict{String, Any}(
                  "bus_type"=>1,
                  "vmin"=>lfac_vmin,
                  "vmax"=>lfac_vmax,
                  "bus_i"=>seriescomp_bus,
                  "index"=>seriescomp_bus,
                  "va"=>0.0,
                  "vm"=>1.0,
                  "name"=>"SC Bus\tLF",
                  "area"=>data["bus"]["$seriescomp_bus"]["area"],
                  "zone"=>data["bus"]["$seriescomp_bus"]["zone"],
                  "base_kv"=>data["bus"]["$seriescomp_bus"]["base_kv"]
               )
            )
            # if "apply_ang_lim" in keys(data["bus"]["$seriescomp_bus"])
            #    lfac_case["bus"]["$seriescomp_bus"]["apply_ang_lim"] = data["bus"]["$seriescomp_bus"]["apply_ang_lim"]
            # end
            # Branch Data
            lfac_branch = deepcopy(branch)
            lfac_branch["f_bus"] = sc_f_bus
            lfac_branch["t_bus"] = sc_t_bus

            lfac_next_branch_copy = deepcopy(lfac_next_branch)
            lfac_next_branch_copy["f_bus"] = next_f_bus
            lfac_next_branch_copy["t_bus"] = next_t_bus
            # Change buses for alt angle limits to corresponding LFAC buses
            if "alt_ang_lim" in keys(lfac_branch)
               lfac_branch["alt_ang_lim"] = [lfac_buslookup[bus] for bus in lfac_branch["alt_ang_lim"]]
            end
            if "alt_ang_lim" in keys(lfac_next_branch_copy)
               lfac_next_branch_copy["alt_ang_lim"] = [lfac_buslookup[bus] for bus in lfac_next_branch_copy["alt_ang_lim"]]
            end

            lfac_case["branch"] = Dict{String, Any}(
               "$(lfac_branch["index"])"=>lfac_branch,
               "$(lfac_next_branch_copy["index"])"=>lfac_next_branch_copy
            )
            # Gen data
            lfac_case["gen"] = Dict{String, Any}(
               "1"=>Dict{String, Any}(
                  "gen_status"=>0,
                  "index"=>1,
                  "gen_bus"=>1,
                  "vg"=>1.0
               )
            )
            # Other data
            lfac_case["storage"] = Dict{String, Any}()
            lfac_case["switch"] = Dict{String, Any}()
            lfac_case["shunt"] = Dict{String, Any}()
            lfac_case["load"] = Dict{String, Any}()
            lfac_case["dcline"] = Dict{String, Any}()

            if fileext == ".m"
               io = open("$output_dir/lfac_$(base_name).m", "w");
               export_matpower(io, lfac_case)
               close(io)
            else
               stringnet = JSON.json(lfac_case)
               open("$output_dir/lfac_$(base_name)$(fileext)", "w") do f
                  write(f, stringnet)
               end
            end

         else # is not a series compensation branch

            # Check if branch is connected to a series compensation bus.
            # If so, skip the branch so that series compensated branches are not
            # split up in different subnetworks.
            sc_connected = false
            for (other_br_idx,other_branch) in data["branch"]
               # println("sourceid: $(other_branch["source_id"])")
               if (length(other_branch["source_id"])>3) && (Unicode.normalize(other_branch["source_id"][1], casefold=true) == "branch") && (Unicode.normalize(other_branch["source_id"][4], casefold=true) == "sc") # is a series compensation branch
                  # println("series compensation branch $(other_branch["index"])")
                  # println("other_branch: $other_branch")
                  if "sc_bus" ∉ keys(other_branch)
                     other_branch["sc_bus"] = other_branch["source_id"][2]
                  end
                  seriescomp_bus = other_branch["sc_bus"]
                  if (branch["f_bus"] == seriescomp_bus) || (branch["t_bus"] == seriescomp_bus)
                     sc_connected = true
                     break
                  end
               end
            end
            if sc_connected
               println("skipped branch $(branch["index"]) with series compensation")
               continue
            end
            # create output directory
            output_dir = "$output_location/br$(br_idx)"
            if !isdir(output_dir)
               mkpath(output_dir)
            end
            # Set two interface buses for base network
            subnet_index = 1
            iface_index = 1
            base_f_bus = branch["f_bus"]
            interface_params["index"] = iface_index
            interface_params["subnet_index"] = subnet_index
            interface_params["bus"] = base_f_bus
            push!(iface_df, interface_params)

            f_area = copy(data["bus"]["$base_f_bus"]["area"])
            f_zone = copy(data["bus"]["$base_f_bus"]["zone"])
            f_kv_base = copy(data["bus"]["$base_f_bus"]["base_kv"])
            iface_index = 2
            base_t_bus = branch["t_bus"]
            interface_params["index"] = iface_index
            interface_params["subnet_index"] = subnet_index
            interface_params["bus"] = base_t_bus
            push!(iface_df, interface_params)
            t_area = copy(data["bus"]["$base_t_bus"]["area"])
            t_zone = copy(data["bus"]["$base_t_bus"]["zone"])
            t_kv_base = copy(data["bus"]["$base_t_bus"]["base_kv"])
            for (bus_idx,bus) in data["bus"]
               if "name" in keys(bus)
                  bus["name"] = split(bus["name"], ";")[1]
               end
            end
            # Disable branch and save network
            branch["br_status"] = 0
            if fileext == ".m"
               io = open("$output_dir/$(base_name).m", "w");
               export_matpower(io, data)
               close(io)
            else
               stringnet = JSON.json(data)
               open("$output_dir/$(base_name)$(fileext)", "w") do f
                  write(f, stringnet)
               end
            end
            # Restore branch status
            branch["br_status"] = 1
            # Create and save LFAC subnetwork
            lfac_case = Dict{String,Any}()
            lfac_case["name"] = data["name"]
            lfac_case["source_version"] = data["source_version"]
            lfac_case["baseMVA"] = data["baseMVA"]
            lfac_case["per_unit"] = true
            # Bus data
            lfac_case["bus"] = Dict{String, Any}(
               "1"=>Dict{String, Any}(
                  "bus_type"=>3,
                  "vmin"=>lfac_vmin,
                  "vmax"=>lfac_vmax,
                  "bus_i"=>1,
                  "index"=>1,
                  "va"=>0.0,
                  "vm"=>1.0,
                  "name"=>"Bus 1\tLF",
                  "area"=>f_area,
                  "zone"=>f_zone,
                  "base_kv"=>f_kv_base
               ),
               "2"=>Dict{String, Any}(
                  "bus_type"=>1,
                  "vmin"=>lfac_vmin,
                  "vmax"=>lfac_vmax,
                  "bus_i"=>2,
                  "index"=>2,
                  "va"=>0.0,
                  "vm"=>1.0,
                  "name"=>"Bus 2\tLF",
                  "area"=>t_area,
                  "zone"=>t_zone,
                  "base_kv"=>t_kv_base
               )
            )
            # Branch Data
            lfac_branch = copy(branch)
            lfac_branch["f_bus"] = 1
            lfac_branch["t_bus"] = 2
            lfac_case["branch"] = Dict{String, Any}(
               "$(lfac_branch["index"])"=>lfac_branch
            )
            # Gen data
            lfac_case["gen"] = Dict{String, Any}(
               "1"=>Dict{String, Any}(
                  "gen_status"=>0,
                  "index"=>1,
                  "gen_bus"=>1,
                  "vg"=>1.0
               )
            )
            # Other data
            lfac_case["storage"] = Dict{String, Any}()
            lfac_case["switch"] = Dict{String, Any}()
            lfac_case["shunt"] = Dict{String, Any}()
            lfac_case["load"] = Dict{String, Any}()
            lfac_case["dcline"] = Dict{String, Any}()
            if fileext == ".m"
               io = open("$output_dir/lfac_$(base_name).m", "w");
               export_matpower(io, lfac_case)
               close(io)
            else
               stringnet = JSON.json(lfac_case)
               open("$output_dir/lfac_$(base_name)$(fileext)", "w") do f
                  write(f, stringnet)
               end
            end

            lfac_buslookup = Dict(
               base_f_bus=>1,
               base_t_bus=>2
            )

         end

         br_idx_int = parse(Int64,br_idx)
         if br_idx_int in keys(ctg_br_dict_tmp)
            bus_i = lfac_buslookup[ctg_br_dict_tmp[br_idx_int]["i"]]
            bus_j = lfac_buslookup[ctg_br_dict_tmp[br_idx_int]["j"]]
            lfac_ctg_br_dict[br_idx_int] = pop!(ctg_br_dict_tmp, br_idx_int)
            lfac_ctg_br_dict[br_idx_int]["i"] = bus_i
            lfac_ctg_br_dict[br_idx_int]["j"] = bus_j
            lfac_ctg_br_dict[br_idx_int]["ckt"] = 1
         end

         # Set two interface buses for LFAC line
         subnet_index = 2
         # LFAC "from bus" inerface
         iface_index = 1
         bus = 1
         interface_params["index"] = iface_index
         interface_params["subnet_index"] = subnet_index
         interface_params["bus"] = bus
         push!(iface_df, interface_params)
         # LFAC "to bus" inerface
         iface_index = 2
         bus = 2
         interface_params["index"] = iface_index
         interface_params["subnet_index"] = subnet_index
         interface_params["bus"] = bus
         push!(iface_df, interface_params)
         # Save interfaces file
         interfaces_filename = "interfaces_file" in names(base_subnet_df) ? base_subnet_df.interfaces_file[1] : "interfaces.csv"
         CSV.write("$output_dir/$interfaces_filename", iface_df)
         # Save subnetworks file
         CSV.write("$output_dir/subnetworks.csv", subnet_df)

         # println("output_dir: $output_dir")
         if "contingency_file" in names(subnet_df)
            write_branch_ctg_file(ctg_br_dict_tmp, output_dir*"/"*ctg_filename)
            write_branch_ctg_file(lfac_ctg_br_dict, output_dir*"/"*"lfac_"*ctg_filename)
         end
      end
   end
end

function enumerate_branches(
      base_network::String,
      output_location::String,
      fbase;
      subnet_params=Dict(),
      interface_params=Dict(),
      indices=[],
      output_type="input",
      noloss=false
   )
   println("The function enumerate_branches is deprecated. Please use upgrade_branches instead.")
   upgrade_branches(
      base_network,
      output_location,
      fbase;
      subnet_params,
      interface_params,
      indices,
      output_type,
      noloss
   )
end

function write_branch_ctg_file(ctg_dict, filename)
   open(filename, "w") do io
      for (idx,ctg) in ctg_dict
         output_str = "CONTINGENCY $(ctg["label"])\n$(uppercase(ctg["action"])) $(uppercase(ctg["component"])) FROM BUS $(ctg["i"]) TO BUS $(ctg["j"]) CIRCUIT $(ctg["ckt"])\nEND\n"
         write(io, output_str)
      end
      write(io,"END\n")
   end
end

function change_angle_lims(base_network, output_location, angle_lim_deg, indices=[]; output_type=".json")
   filename = split(base_network, "/")[end]
   prefix = base_network[1:(end-length(filename))]
   println("$filename")

   filetype = split(lowercase(filename), '.')[end]
   base_name = filename[1:(end-length(filetype)-1)]
   (base_name, filetype) = splitext(filename)
   filetype = Unicode.normalize(filetype, casefold=true)
   fileext = filetype == ".m" ? filetype : ".json"
   if lowercase(output_type) == ".json"
      fileext = ".json"
   end

   output_dir = "$output_location/"
   if !isdir(output_dir)
      mkpath(output_dir)
   end
   cp("$prefix/interfaces.csv", "$output_dir/interfaces.csv", force=true)
   cp("$prefix/subnetworks.csv", "$output_dir/subnetworks.csv", force=true)
   if isfile("$prefix/contingencies.con")
      cp("$prefix/contingencies.con", "$output_dir/contingencies.con", force=true)
   end

   data = PowerModels.parse_file(base_network)

   for (i,(br_idx,branch)) in enumerate(data["branch"])
      if ((length(indices) == 0) || (parse(Int64, br_idx) in indices))
         branch["angmin"] = -angle_lim_deg * pi/180.0
         branch["angmax"] = angle_lim_deg * pi/180.0
      end
   end
   if fileext == ".m"
      io = open("$output_dir/$(base_name).m", "w");
      export_matpower(io, data)
      close(io)
   else
      stringnet = JSON.json(data)
      open("$output_dir/$(base_name)$(fileext)", "w") do f
         write(f, stringnet)
      end
   end

   if isfile("$prefix/lfac_$filename")
      data = PowerModels.parse_file("$prefix/lfac_$filename")
      for (i,(br_idx,branch)) in enumerate(data["branch"])
         if ((length(indices) == 0) || (parse(Int64, br_idx) in indices))
            branch["angmin"] = -angle_lim_deg * pi/180.0
            branch["angmax"] = angle_lim_deg * pi/180.0
         end
      end
      if fileext == ".m"
         io = open("$output_dir/$(base_name).m", "w");
         export_matpower(io, data)
         close(io)
      else
         stringnet = JSON.json(data)
         open("$output_dir/$(base_name)$(fileext)", "w") do f
            write(f, stringnet)
         end
      end
   end
end

function change_lims(base_network, output_location, angle_lim_deg, indices=[])
   filename = split(base_network, "/")[end]
   prefix = base_network[1:(end-length(filename))]
   println("$filename")

   output_dir = "$output_location/"
   if !isdir(output_dir)
      mkpath(output_dir)
   end
   cp("$prefix/interfaces.csv", "$output_dir/interfaces.csv", force=true)
   cp("$prefix/subnetworks.csv", "$output_dir/subnetworks.csv", force=true)

   data = PowerModels.parse_file(base_network)

   for (i,(br_idx,branch)) in enumerate(data["branch"])
      if ((length(indices) == 0) || (parse(Int64, br_idx) in indices))
         branch["angmin"] = -angle_lim_deg * pi/180.0
         branch["angmax"] = angle_lim_deg * pi/180.0
      end
   end
   io = open("$output_dir/$filename", "w");
   export_matpower(io, data)
   close(io)

   if isfile("$prefix/lfac_$filename")
      data = PowerModels.parse_file("$prefix/lfac_$filename")
      for (i,(br_idx,branch)) in enumerate(data["branch"])
         if ((length(indices) == 0) || (parse(Int64, br_idx) in indices))
            branch["angmin"] = -angle_lim_deg * pi/180.0
            branch["angmax"] = angle_lim_deg * pi/180.0
         end
      end
      io = open("$output_dir/lfac_$filename", "w");
      export_matpower(io, data)
      close(io)
   end
end

function scale_load(
      scaling,
      areas,
      directory,
      objective,
      x_axis,
      y_axis;
      gen_areas=[],
      area_transfer=[],
      gen_zones=[],
      zone_transfer=[],
      plot_vert_line::Tuple=([],""),
      plot_horiz_line::Tuple=([],""),
      xlimits::Array{Any,1}=[],
      ylimits::Array{Any,1}=[],
      output_plot_label::Tuple{String,String}=("",""),
      scopf::Bool=false,
      contingency::Int64=0
   )

   scale_all_areas = false
   if length(areas) == 0
      scale_all_areas = true
   end
   # get load pd and qd from areas of interest and put them in dict_filt
   sn_data = read_sn_data(directory)
   dict_filt = Dict()
   for (subnet_idx,sn_subnet) in sn_data["sn"]
      println("number of loads: $(length(sn_subnet["load"]))")
      for (load_idx,load) in sn_subnet["load"]
         if (sn_subnet["bus"]["$(load["load_bus"])"]["area"] in areas) || (scale_all_areas)
            pd_keys = ["sn",subnet_idx,"load",load_idx,"pd"]
            pd = load["pd"]
            set_nested!(dict_filt, pd_keys, pd)
            qd_keys = ["sn",subnet_idx,"load",load_idx,"qd"]
            qd = load["qd"]
            set_nested!(dict_filt, qd_keys, qd)
         end
      end
   end

   # get the nested keys values
   (keys,values) = traverse_nested(dict_filt)

   # apply scaling
   # println("keys: $keys")
   # println("values: $values")
   values_scaled = collect(eachrow(values .* scaling'))
   # println("values_scaled: $values_scaled")
   params = (keys,values_scaled)
   println("number of scaled loads: $(length(values_scaled))")
   (results_dict, output_plot) = run_multiple_params(
      directory,
      objective,
      x_axis,
      y_axis,
      params,
      gen_areas,
      area_transfer,
      gen_zones,
      zone_transfer,
      plot_vert_line,
      plot_horiz_line,
      xlimits,
      ylimits,
      output_plot_label,
      scopf,
      contingency
   )
   return (results_dict, output_plot)
end


text:
filename = tempname()
# The mode "c" stands for creating a new file (clobber)
ds = NCDatasets.Dataset(filename,"c")

# define the dimension "lon" and "lat" with the size 10 and 11 resp.
ds.dim["lon"] = 10
ds.dim["lat"] = 11

v = NCDatasets.defVar(ds,"var_with_missing_data",UInt8,("lon","lat"))

data = trues(10,11)

# previously we got an infinite recursion
v[1:10,1:11] = data
@test Bool.(v[:,:]) ≈ data

v[:,:] = data
@test Bool.(v[:,:]) ≈ data

v[:] = data
@test Bool.(v[:,:]) ≈ data

NCDatasets.close(ds)


text:
#a1 = Dict(1=>"one", 2=>"two")

a1 = Dict{String, Int32}()

f = open( "test.bed", "r" )
n = countlines( f )  # 计算数据行数
seekstart( f )       # 返回文件起始位置

for i = 1:n
    s = split( readline(f), "\t")
    a = parse(Int64, s[2])
    b = parse(Int64, s[3])

    for pos = (a-125):(b+125)
        p_x = string(s[1], pos)
        a1[p_x] = 1
    end

end

println(length(a1))


text:
using GLWindow, ModernGL, GLAbstraction, Reactive, GeometryTypes, Colors, GLVisualize
using Gtk, ModernGL

resolution = (600, 500)
name = "GTK + GLVisualize"
parent = Gtk.Window(name, resolution..., true, true)
Gtk.visible(parent, true)
Gtk.setproperty!(parent, Symbol("is-focus"), false)
box = Gtk.Box(:v)
push!(parent, box)

sl = Gtk.Scale(false, linspace(0, 1, 100))
push!(box, sl)
gl_area = Gtk.GLArea()
Gtk.gl_area_set_required_version(gl_area, 3, 3)
push!(box, gl_area)
Gtk.setproperty!(box, :expand, gl_area, true)
adj = Gtk.Adjustment(sl)

Gtk.signal_connect(gl_area, "render") do gl_area, gdk_context
    val = getproperty(adj, :value, Float64)
    glClearColor(val, 0,0,1)
    glClear(GL_COLOR_BUFFER_BIT)
    glFlush()
    return false
end

Gtk.showall(parent)
# Gtk.gtk_main()

prerender=()->glDisable(GL_DEPTH_TEST), # draw over other items
postrender=()->glEnable(GL_DEPTH_TEST)

_view(visualize(
    points, :lines, indices = indices,
    prerender=()->glDisable(GL_DEPTH_TEST), # draw over other items
    postrender=()->glEnable(GL_DEPTH_TEST)
), w, camera = :fixed_pixel)


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using ValueShapes, Distributions

@testset "variate_transform" begin
    @testset "identity_vt" begin
        @test @inferred(BAT.IdentityVT(ScalarShape{Real}())) isa BAT.VariateTransform{<:ScalarShape,<:ScalarShape}
        @test @inferred(BAT.IdentityVT(ArrayShape{Int}(2))) isa BAT.VariateTransform{<:ArrayShape,<:ArrayShape}
        @test @inferred(BAT.IdentityVT(ArrayShape{Float32}(2, 3))) isa BAT.VariateTransform{<:ArrayShape,<:ArrayShape}

        ntshape = NamedTupleShape(a = ScalarShape{Real}(), b = ArrayShape{Int}(2))
        @test @inferred(BAT.IdentityVT(ntshape)) isa BAT.VariateTransform{<:ValueShapes.NamedTupleShape{(:a, :b)},<:ValueShapes.NamedTupleShape{(:a, :b)}}

        nttrafo = BAT.IdentityVT(ntshape)
        ntvalue = (a = 4.2, b = [5, 7])

        @test @inferred(varshape(nttrafo)) == ntshape
        @test @inferred(inv(nttrafo)) === nttrafo

        @test @inferred(BAT.apply_vartrafo(nttrafo, ntvalue, 0)).v === ntvalue
        @test BAT.apply_vartrafo(nttrafo, ntvalue, 0.79).ladj == 0.79
        @test ismissing(BAT.apply_vartrafo(nttrafo, ntvalue, missing).ladj)
        @test @inferred(BAT.apply_vartrafo(inv(nttrafo), ntvalue, 0.79)).v === ntvalue
        @test BAT.apply_vartrafo(inv(nttrafo), ntvalue, 0.79).ladj == 0.79
        @test ismissing(BAT.apply_vartrafo(inv(nttrafo), ntvalue, missing).ladj)

        @test @inferred(nttrafo ∘ nttrafo) === nttrafo

        arrayshape_trafo = @inferred(BAT.IdentityVT(ArrayShape{Real}(1)))
        hoeldertable_trafo = @inferred((arrayshape_trafo)(HoelderTable(), volcorr=Val(false)))

        @test hoeldertable_trafo.volcorr == @inferred(BAT.TDNoCorr())
        @test hoeldertable_trafo.trafo isa BAT.IdentityVT
    end
end