text:
"""
Basic continuous univariate parameter state type

# Constructors

## BasicContUnvParameterState{N<:Real}(value::N, <optional arguments>)

Construct a basic continuous univariate parameter state with some ``value``.

### Optional arguments:

* ``diagnostickeys::Vector{Symbol}=Symbol[]``: the diagnostic keys of the state.
* ``diagnosticvalues::Vector=Array(Any, length(diagnostickeys))``: the diagnostic values of the state.

## BasicContUnvParameterState{N<:Real}(<optional arguments>)

Construct a basic continuous univariate parameter state with an uninitialized ``value`` (``NaN``).

### Optional arguments:

* ``diagnostickeys::Vector{Symbol}=Symbol[]``: the diagnostic keys of the state.
* ``::Type{N}=Float64``: the element type of the state value.
* ``diagnosticvalues::Vector=Array(Any, length(diagnostickeys))``: the diagnostic values of the state.

# Examples

```julia
julia> state = BasicContUnvParameterState(-1.25, [:accept], [false])
Klara.BasicContUnvParameterState{Float64}(-1.25,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,Bool[false],[:accept])

julia> state.value
-1.25

julia> diagnostics(state)
Dict{Symbol,Bool} with 1 entry:
  :accept => false
```
"""
mutable struct BasicContUnvParameterState{N<:Real} <: ParameterState{Continuous, Univariate}
  "Vector value of basic continuous univariate parameter state"
  value::N
  "Value of log-likelihood at the state's value"
  loglikelihood::N
  "Value of log-prior at the state's value"
  logprior::N
  "Value of log-target at the state's value"
  logtarget::N
  "Value of gradient of log-likelihood at the state's value"
  gradloglikelihood::N
  "Value of gradient of log-prior at the state's value"
  gradlogprior::N
  "Value of gradient of log-target at the state's value"
  gradlogtarget::N
  "Value of metric tensor of log-likelihood at the state's value"
  tensorloglikelihood::N
  "Value of metric tensor of log-prior at the state's value"
  tensorlogprior::N
  "Value of metric tensor of log-target at the state's value"
  tensorlogtarget::N
  "Value of derivatives of metric tensor of log-likelihood at the state's value"
  dtensorloglikelihood::N
  "Value of derivatives of metric tensor of log-prior at the state's value"
  dtensorlogprior::N
  "Value of derivatives of metric tensor of log-target at the state's value"
  dtensorlogtarget::N
  "Diagnostic values associated with the sampling of the state"
  diagnosticvalues::Vector
  "Diagnostic keys associated with the sampling of the state"
  diagnostickeys::Vector{Symbol}
  diffmethods::Union{DiffMethods, Void}
  diffstate::Union{DiffState, Void}
end

function BasicContUnvParameterState(
  value::N,
  diagnostickeys::Vector{Symbol}=Symbol[],
  diffmethods::Union{DiffMethods, Void}=nothing,
  diffopts::Union{DiffOptions, Void}=nothing,
  diagnosticvalues::Vector=Array{Any}(length(diagnostickeys))
) where N<:Real
  v = convert(N, NaN)

  diffstate = DiffState()

  if diffopts != nothing
    if diffopts.order == 1
      for (i, field) in ((1, :resultll), (2, :resultlp), (3, :resultlt))
        if diffopts.targets[i]
          setfield!(diffstate, field, DiffResults.DiffResult(zero(value), zero(value)))
        end
      end
    else
      error("Derivative order must be 1, got order=$order")
    end
  end

  BasicContUnvParameterState{N}(
    value, v, v, v, v, v, v, v, v, v, v, v, v, diagnosticvalues, diagnostickeys, diffmethods, diffstate
  )
end

BasicContUnvParameterState(
  diagnostickeys::Vector{Symbol}=Symbol[],
  ::Type{N}=Float64,
  diffmethods::Union{DiffMethods, Void}=nothing,
  diffopts::Union{DiffOptions, Void}=nothing,
  diagnosticvalues::Vector=Array{Any}(length(diagnostickeys))
) where {N<:Real} =
  BasicContUnvParameterState(convert(N, NaN), diagnostickeys, diffmethods, diffopts, diagnosticvalues)

value_support{N<:Real}(::Type{BasicContUnvParameterState{N}}) = Continuous
value_support(::BasicContUnvParameterState) = Continuous

variate_form{N<:Real}(::Type{BasicContUnvParameterState{N}}) = Univariate
variate_form(::BasicContUnvParameterState) = Univariate

eltype{N<:Real}(::Type{BasicContUnvParameterState{N}}) = N
eltype{N<:Real}(::BasicContUnvParameterState{N}) = N

generate_empty(
  state::BasicContUnvParameterState,
  diffmethods::Union{DiffMethods, Void}=nothing,
  diffopts::Union{DiffOptions, Void}=nothing
) = BasicContUnvParameterState(state.diagnostickeys, eltype(state), diffmethods, diffopts)


text:
using Pkg; Pkg.activate(pwd()); Pkg.instantiate()
using Revise

@time using TumorGrowth

c = TumorGrowth.HyperCube(zeros(3), 10)
lattice = TumorGrowth.triangular_lattice(c; n=11)

sc = TumorGrowth.GLMakie.meshscatter(
    reshape(TumorGrowth.Point{3}.(lattice), prod(size(lattice))), markersize=0.1,
    )
display(sc)

lattice

lattice[1]
sc = plotting( DataFrame(position = reshape(lattice, prod(size(lattice)))),
    color = :black, markersize=0.1)
TumorGrowth.GLMakie.zlims!(sc, (0,4))

lattice

ext: jl

max_stars_repo_name: aangaji/TumorGrowth

max_issues_repo_name: aangaji/TumorGrowth

max_forks_repo_name: aangaji/TumorGrowth


text:
#                    N    Ne   limit
const ABC_bounds  = [10   0.1     0;
                     500  1.0   1000]

#                    N    K   η_max
const ECA_bounds  = [ 10  2    0.0;
                      500 10    4]

#                     N  F  CR
const DE_bounds   = [10  0  0;
                     500 2  1.0]

#                    N   C1   C2   omega
const PSO_bounds  = [10    0   0    0;
                     500   4   4    1.0]

function ABC_errors(Φ, benchmark)
    Errors_shared = SharedArray{Float64}(length(benchmark), NRUNS)
    Cost_shared = SharedArray{Int}(length(benchmark), NRUNS)

    Φ_ = SharedArray{Float64}([round(Φ[:N]), Φ[:Ne], round(Φ[:limit])])

    @sync @distributed for fn in 1:length(benchmark)
        N = Int(Φ_[1])
        Ne = round(Int, Φ_[2]*N)
        nruns = 0
        ff(x) = begin
            nruns += 1
            benchmark[fn].F(x)
        end
        for r = 1:NRUNS
            x, fx = ABC(ff, benchmark[fn].bounds_ul;
                                N = N,
                            limit = Int(Φ_[3]),
                               Ne = Ne,
                               No = N - Ne,
                                max_evals=max_NFEs,
                                termination=accuracy_termination,
                                )

            Errors_shared[fn, r] = fx < desired_accu ? 0.0 : fx
            Cost_shared[fn, r] = nruns
            nruns *= 0
        end
    end

    Errors = Matrix(Errors_shared)
    Cost = Matrix(Cost_shared)

    return  [Errors, Cost]

end

function ECA_errors(Φ, benchmark)
    Errors_shared = SharedArray{Float64}(length(benchmark), NRUNS)
    Cost_shared = SharedArray{Int}(length(benchmark), NRUNS)

    Φ_ = SharedArray{Float64}([ round(Φ[:N]), round(Φ[:K]), Φ[:eta]  ])

    @sync @distributed for fn in 1:length(benchmark)
        nruns = 0
        ff(x) = begin
            nruns += 1
            benchmark[fn].F(x)
        end
        for r = 1:NRUNS
            x, fx = eca(ff, D;
                                N = Int(Φ_[1]),
                                K = Int(Φ_[2]),
                                η_max = Φ_[3],
                                max_evals=max_NFEs,
                                termination=accuracy_termination,
                                limits  =  benchmark[fn].bounds_ul,
                                showResults=false)

            Errors_shared[fn, r] = fx < desired_accu ? 0.0 : fx
            Cost_shared[fn, r] = nruns
            nruns *= 0
        end
    end

    Errors = Matrix(Errors_shared)
    Cost = Matrix(Cost_shared)

    return  [Errors, Cost]

end

function DE_errors(Φ, benchmark)
    Errors_shared = SharedArray{Float64}(length(benchmark), NRUNS)
    Cost_shared = SharedArray{Int}(length(benchmark), NRUNS)
    Φ_ = SharedArray{Float64}([ round(Φ[:N]), Φ[:F], Φ[:CR] ])

    @sync @distributed for fn in 1:length(benchmark)
        nruns = 0
        ff(x) = begin
            nruns += 1
            benchmark[fn].F(x)
        end
        for r = 1:NRUNS
            x, fx = DE(ff, D;
                                N  = Int(Φ_[1]),
                                F  = Φ_[2],
                                CR = Φ_[3],
                                max_evals=max_NFEs,
                                termination=accuracy_termination,
                                limits  =  benchmark[fn].bounds_ul,
                                showResults=false)

            Errors_shared[fn, r] = fx < desired_accu ? 0.0 : fx
            Cost_shared[fn, r] = nruns
            nruns *= 0
        end
    end

    Errors = Matrix(Errors_shared)
    Cost = Matrix(Cost_shared)

    return  [Errors, Cost]

end

function PSO_errors(Φ, benchmark)
    Errors_shared = SharedArray{Float64}(length(benchmark), NRUNS)
    Cost_shared = SharedArray{Int}(length(benchmark), NRUNS)
    Φ_ = SharedArray{Float64}( [ round(Φ[:N]), Φ[:C1], Φ[:C2], Φ[:omega] ] )

    @sync @distributed for fn in 1:length(benchmark)
        nruns = 0
        ff(x) = begin
            nruns += 1
            benchmark[fn].F(x)
        end
        for r = 1:NRUNS
            x, fx = pso(ff, D;
                                N  = Int(Φ_[1]),
                                C1 = Φ_[2],
                                C2 = Φ_[3],
                                ω  = Φ_[4],
                                max_evals=max_NFEs,
                                termination=accuracy_termination,
                                limits  =  benchmark[fn].bounds_ul,
                                showResults=false)

            Errors_shared[fn, r] = fx < desired_accu ? 0.0 : fx
            Cost_shared[fn, r] = nruns
            nruns *= 0
        end
    end

    Errors = Matrix(Errors_shared)
    Cost = Matrix(Cost_shared)
    # Cost = Matrix(Cost_shared)

    return  [Errors, Cost]

end

function getErrors(Φ, benchmark, name)
    if name == "abc"
        return ABC_errors(Φ, benchmark)
    elseif name == "eca"
        return ECA_errors(Φ, benchmark)
    elseif name == "de"
        return DE_errors(Φ, benchmark)
    elseif name == "pso"
        return PSO_errors(Φ, benchmark)
    end
end


text:
using UnivariateFunctions: create_linear_interpolation, create_constant_interpolation_to_right, create_quadratic_spline
using UnivariateFunctions: create_constant_interpolation_to_left, years_from_global_base, evaluate
using Dates

tol = 10*eps()
const global_base_date = Date(2000,1,1)
StartDate = Date(2018, 7, 21)

x = StartDate .+ Dates.Day.(2 .* (1:1000 .- 1))

function ff(x::Date)
    days_between = years_from_global_base(x)
    return log(days_between) + sqrt(days_between)
end
y = ff.(x)

spline = create_linear_interpolation(x,y)
# Test if interpolating
all(abs.(evaluate.(Ref(spline), x) .- y) .< tol)

x_float = years_from_global_base.(x)
coefficient_in_first_interval = (y[2] - y[1])/(x_float[2] - x_float[1])
(coefficient_in_first_interval - spline.functions_[1].functions_[2].a_) < tol
other_coefficient = (y[11] - y[10])/(x_float[11] - x_float[10])
(other_coefficient - spline.functions_[10].functions_[2].a_) < tol

left_const = create_constant_interpolation_to_left(x,y)
right_const = create_constant_interpolation_to_right(x,y)
all(abs.(evaluate.(Ref(left_const), x_float .- tol) .- evaluate.(Ref(right_const), x_float)) .< tol)

# Making interpolation splies with date periods.
xx = Array{Month,1}([Month(1), Month(2), Month(3)])
xx2 = Array{DatePeriod,1}([Month(1), Month(2), Year(1)])
yy = [1,2,3]
abs(create_constant_interpolation_to_left(xx,yy)(Day(45)) - 2.0) < eps()
abs(evaluate(create_constant_interpolation_to_right(xx,yy), Day(15)) - 1.0) < eps()
abs(evaluate(create_linear_interpolation(xx,yy), Month(1)) - 1.0) < eps()
abs(evaluate(create_quadratic_spline(xx,yy), Month(2)) - 2.0) < eps()

abs(evaluate(create_constant_interpolation_to_left(xx2,yy), Day(45)) - 2.0) < eps()
abs(evaluate(create_constant_interpolation_to_right(xx2,yy), Day(15)) - 1.0) < eps()
abs(evaluate(create_linear_interpolation(xx2,yy), Month(1)) - 1.0) < eps()
abs(evaluate(create_quadratic_spline(xx2,yy), Month(2)) - 2.0) < eps()


text:
using LaTeX
import LaTeX: getrequirements
using Test, Dates

# Check that package dependency is working
@test haskey(
    getrequirements(Section("Code", [Code("""test_getrequirements()""")])),
    "texments")

@testset "Test a minimal document" begin

    doc = document([
        DocumentClass("article", ["11pt", "letterpaper"]),
        Date(2015, 11, 23),
        Title("Sample Article"),
        Author("Andi Andromeda"),
        Section("Code", [Code("""# minimal example""")])])
    @test split(doc, '\n')[1] == "\\documentclass[11pt,letterpaper]{article}"
    @test occursin("\\date{2015-11-23}", doc)
    @test occursin("\\author{Andi Andromeda}", doc)
    @test occursin("\\section", doc)  # articles have section as top-level
    @test occursin("\\usestyle{default}", doc)  # needed for code

    # Test a minimal report
    rep = report(Section("Test", "Lorum ipsum."))
    @test split(rep, '\n')[1] == "\\documentclass[11pt,a4paper]{report}"
    @test occursin("\\chapter", rep)  # reports have chapter as top-level

end


text:
module UnifiedChineseDictionary

end


text:
mutable struct GFGenerateState
    trace::DynamicDSLTrace
    constraints::ChoiceMap
    weight::Float64
    visitor::AddressVisitor
    params::Dict{Symbol,Any}
end

function GFGenerateState(gen_fn, args, constraints, params)
    trace = DynamicDSLTrace(gen_fn, args)
    GFGenerateState(trace, constraints, 0., AddressVisitor(), params)
end

function traceat(state::GFGenerateState, dist::Distribution{T},
              args, key) where {T}
    local retval::T

    # check that key was not already visited, and mark it as visited
    visit!(state.visitor, key)

    # check for constraints at this key
    constrained = has_value(state.constraints, key)
    !constrained && check_no_submap(state.constraints, key)

    # get return value
    if constrained
        retval = get_value(state.constraints, key)
    else
        retval = random(dist, args...)
    end

    # compute logpdf
    score = logpdf(dist, retval, args...)

    # add to the trace
    add_choice!(state.trace, key, retval, score)

    # increment weight
    if constrained
        state.weight += score
    end

    retval
end

function traceat(state::GFGenerateState, gen_fn::GenerativeFunction{T,U},
              args, key) where {T,U}
    local subtrace::U
    local retval::T

    # check key was not already visited, and mark it as visited
    visit!(state.visitor, key)

    # check for constraints at this key
    constraints = get_submap(state.constraints, key)

    # get subtrace
    (subtrace, weight) = generate(gen_fn, args, constraints)

    # add to the trace
    add_call!(state.trace, key, subtrace)

    # update weight
    state.weight += weight

    # get return value
    retval = get_retval(subtrace)

    retval
end

function splice(state::GFGenerateState, gen_fn::DynamicDSLFunction,
                args::Tuple)
    prev_params = state.params
    state.params = gen_fn.params
    retval = exec(gen_fn, state, args)
    state.params = prev_params
    retval
end

function generate(gen_fn::DynamicDSLFunction, args::Tuple,
                    constraints::ChoiceMap)
    state = GFGenerateState(gen_fn, args, constraints, gen_fn.params)
    retval = exec(gen_fn, state, args)
    set_retval!(state.trace, retval)
    (state.trace, state.weight)
end


text:
"""
    download_zip(url::AbstractString, name::AbstractString, dirpath::String; <keyword arguments>)

Download zipped dataset `name` from `url`, and place it to `dirpath`.

# Keyword arguments
- `usecache=true`: if `true` and `dirpath` already exists, it skips download. If `false`, it force download by removing `dirpath` (if exists).
- `unzip=true`: whether extract datasets from zip.
- `removezip=false`: whether remove original zip after extraction.
"""
function download_zip(
    url::AbstractString,
    name::AbstractString,
    dirpath::String;
    usecache = true,
    unzip = true,
    removezip = false,
)
    if usecache && isdir(dirpath)
        return
    end
    if isdir(dirpath)
        rm(dirpath, force = true, recursive = true)
    end

    # download zip
    mkpath(dirpath)
    zippath = joinpath(dirpath, "$(name).zip")

    io = open(zippath, "w")
    try
        HTTP.request("GET", url, response_stream = io)
    catch e
        if isfile(zippath)
            rm(zippath)
        end
        throw(e)
    finally
        close(io)
    end

    # unzip
    if unzip
        r = ZipFile.Reader(zippath)
        try
            for f in r.files
                if f.name[end] == '/'
                    continue
                end
                println("Extracting $(f.name)")
                io = open(joinpath(dirpath, basename(f.name)), "w")
                write(io, read(f, String))
            end
        finally
            close(r)
        end
    end

    # remove zip
    if removezip
        rm(zippath)
    end
end


text:
module test_bukdu_assoc

using Test
using Bukdu: Assoc

assoc = Assoc("A" => "B")
@test !isempty(assoc)

@test pairs(["A" => "B"]) == pairs(assoc)

params = Assoc("x" => "2")
@test haskey(params, :x)
@test get(params, :x, 0) == 2
@test something(params.x, "") == "2"
@test params.x == "2"
@test params[:x] == "2"
@test params["x"] == "2"

@test !haskey(params, :y)
@test get(params, :y, 3) == 3
@test something(params.y, "") == ""
@test params.y === nothing
@test params[:y] === nothing
@test params["y"] === nothing

params[:x] = "3"
@test params.x == "3"
params[:y] = "5"
@test params.y == "5"
@test length(params) == 2
@test keys(params) == ["x", "y"]
@test values(params) == ["3", "5"]

empty!(params)
@test !haskey(params, "x")

assoc = Assoc("a" => 2, "b" => [1,2])
@test !isempty(assoc)
@test assoc.a in assoc.b

# SubString is generated by url param value parsing
assoc = Assoc("z" => SubString("2", 1:1))
@test get(assoc, :z, 0) == 2

for (k, v) in assoc
    @test k == "z"
    @test v == "2"
end

end # module test_bukdu_assoc


text:
function sumofintervals(a::Array{Tuple{T,T}})::T where {T<:Number}
    # (x,false) = start point at x, (x,true) = end point at x
    points = reduce( append!, map( p -> [(p[1],false),(p[2],true)], a ), init=Tuple{T,Bool}[] )
    sort!( points )
    sum = 0
    in_interval = 0
    start_interval = 0
    for p in points
        if p[2] # true = end point
            in_interval -= 1
            if in_interval == 0
                sum += p[1] - start_interval
            end
        else # false = start point
            in_interval += 1
            if in_interval == 1
                start_interval = p[1]
            end
        end
    end
    sum
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
import Random
import Zygote
import RestrictedBoltzmannMachines as RBMs

using Test: @test, @testset, @test_broken, @inferred
using Statistics: mean, var
using SpecialFunctions: erfcx
using Distributions: truncated, Normal

for a = -10:10
    d = truncated(Normal(); lower=a)
    @test RBMs.tnmean(a) ≈ mean(d)
    @test RBMs.tnvar(a)  ≈ var(d)
    @test a ≤ RBMs.tnmean(a) < Inf
    @test 0 ≤ RBMs.tnvar(a) ≤ 1
end
@test 1e20 ≤ RBMs.tnmean(1e20) < Inf
@test_broken 0 ≤ RBMs.tnvar(1e80) ≤ 1

@testset "sqrt1half" begin
    @test (@inferred RBMs.sqrt1half(5)) ≈ 5.1925824035672520156
    @test (@inferred RBMs.sqrt1half(0)) == 1
    @test RBMs.sqrt1half(-1) == RBMs.sqrt1half(1) ≈ 1.6180339887498948482
    @test isnan(@inferred RBMs.sqrt1half(NaN))
    @test RBMs.sqrt1half(Inf) == RBMs.sqrt1half(-Inf) == Inf
    @test RBMs.sqrt1half(1e300) ≈ 1e300
end

@testset "randnt" begin
    @test (@inferred RBMs.randnt(0)) > 0
    @test (@inferred RBMs.randnt(1e300)) == 1e300
    @test (@inferred RBMs.randnt(Inf)) == Inf
    @test isnan(@inferred RBMs.randnt(NaN))
    @test (@inferred RBMs.randnt(floatmax(Float64))) == floatmax(Float64)
end

Random.seed!(18)

@inferred RBMs.randnt_half(1.0, 2.0)
@inferred RBMs.randnt_half(Float32(1.0), Float32(2.0))
@test RBMs.randnt_half(Float32(1.0), Float32(2.0)) isa Float32

# compare exact 1st and 2nd moments to Monte Carlo estimates
m1(μ,σ) = μ + σ * √(2/π) / erfcx(-μ/σ/√2)
m2(μ,σ) = μ^2 + σ^2 + μ * σ * √(2/π) / erfcx(-μ/σ/√2)

for μ = -1:1, σ = 1:2
    samples = [RBMs.randnt_half(μ,σ) for _ = 1:10^6]
    @test mean(samples.^1) ≈ m1(μ,σ) atol=1e-2
    @test mean(samples.^2) ≈ m2(μ,σ) atol=1e-2
end

# broadcasted versions
μ = 3randn(2,2); σ = 3rand(2,2)
dμ, dσ = Zygote.gradient(μ,σ) do μ,σ
    mean(m1.(μ,σ))
end


text:
struct IntSeq{T}
  bits::T
end

singleton(::Type{IntSeq{T}}, i) where {T} = ..
@post singleton(ret, ::Type{T}, i) = length(ret) == 1 && first(ret) == i

# @inline combine(a::Vector{T}, b::Vector{T}) where T =
@inline append(a::IntSeq{T}, b::IntSeq{T}) where T = vcat(a, b)
@inline append(a::Vector{T}, b::T) where T = vcat(a, T[b])
@inline singletonid(::Type{Vector{T}}, i::Vararg{T, N}) where {T, N} = T[i...]

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# Functions

export lorem_chars, lorem_words, lorem_sentences, lorem_paragraphs

function lorem_chars(volume)
    if volume > length(lorem_txt)
        throw(BoundsError)
    else
        return lorem_txt[1:volume]
    end
end

function lorem_words(volume)
    temp = split(lorem_txt)
    if volume > length(temp)
        throw(BoundsError)
    else
        return join(temp[1:volume], " ")
    end
end

function lorem_sentences(volume)
    temp = split(lorem_txt, "\.")
    if volume > length(temp)
        throw(BoundsError)
    else
        return join(temp[1:volume], "\.")
    end
end

function lorem_paragraphs(volume)
    temp = split(lorem_txt, "\n")
    if volume > length(temp)
        throw(BoundsError)
    else
        return join(temp[1:volume], "\n")
    end
end

max_stars_repo_name: ColinFay/LoremJulia

max_issues_repo_name: ColinFay/LoremJulia

max_forks_repo_name: ColinFay/LoremJulia


text:
using Test

global flog = open("log.testPrimitive", "a+")

#@testset "Primitive Variations Tests" begin
# modelからfalseを除外するのは適切なのかどうか?

 @testset "1 var Test" begin
  @test solveQ([],[[:a]], [[:a]]) == [[]]
  @test solveQ([],[[:a]], [[:b]]) == []
  @test solveQ([:x],[[:a]], [[:a]]) == [[:x]]
  @test solveQ([:x],[[:b]], [[:a]]) == []
  @test solveQ([:x],[[:x]], [[:a]]) == [[:a]]

  @test solveQ([:x],[[:a],[:b]], [[:a]]) == []
  @test solveQ([:x],[[:a],[:b]], [[:a],[:b]]) == [[:x]]

  @test solveQ([:x],[[:x]], [[:a],[:b]]) == [[:a],[:b]]
  @test solveQ([:x],[[:a],[:b]], [[:b],[:a]]) == [[:x]]

  @test solveQ([:x],[[:a :x],[:x :b]], [[:a :b],[:a :b]]) == []

 end

 @testset "2 vars Test" begin
  @test solveQ([:x,:y],[[:x :y]],[[:a :b]]) == [[:a,:b]]
  @test solveQ([:x,:y],[[:x :b],[:a :y]],[[:a :b]]) == [[:a,:b]]
  @test solveQ([:x,:y],[[:x :b],[:y :a]],[[:a :b]]) == []
  @test solveQ([:x,:y],[[:a :b]],[[:a :b]]) == [[:x,:y]]

# 変数をいれかえても結果は同じになる。
  @test solveQ([:x,:y],[[:x :b]],[[:a :b]]) == [[:a,:y]]
  @test solveQ([:x,:y],[[:y :b]],[[:a :b]]) == [[:x,:a]]

# 2 facts
  @test solveQ([:x,:y],[[:x :b],[:y :a]],[[:a :b],[:b :a]]) == [[:a,:b]]
  @test solveQ([:x,:y],[[:x :b],[:y :b]],[[:a :b]]) == [[:a,:a]]
# [:x :x]の強制力はないのか?
# xとyは独立にそれぞれ別のfactとunifyされるので、このような答えになる。
  @test solveQ([:x,:y],[[:x :x],[:y :y]],[[:a :a],[:b :b]]) == [[:a,:a],[:a,:b],[:b,:a],[:b,:b]]
# これが同一性を示しているように見えるのはたまたまのようだ
  @test solveQ([:x,:y],[[:x :x],[:y :y]],[[:a :a]]) == [[:a,:a]]
#これなら同一性になる
  @test solveQ([:x],[[:x :x]],[[:a :a],[:b :b]]) == [[:a],[:b]]
# ↑は、この世界での同一性を示すインスタンスの集合ということ

# Queryの間のandを考えたいときは、個別のsolveQ()を実行した結果を見るしかなさそう

 end

 @testset "2 gramatical const Test" begin
#これは確かに同一性だが、一般的な「同一性」は表現できていない
  @test solveQ([:x],[[:P :x :x]], [[:P :a :a]]) == [[:a]]
  @test solveQ([:x],[[:P :x :x]], [[:P :a :a],[:P :b :b]]) == [[:a],[:b]]

 end

 @testset "2 grammatical Test" begin
  @test solveQ([:x,:y],[[:x :is :kudamono],[:x :is :red]],[[:ringo :is :kudamono],[:ringo :is :aka],[:mikan :is :yellow]]) == []
  @test solveQ([:x,:y],[[:x :known :kudamono],[:x :is :red]],[[:mikan :known :kudamono],[:ringo :known :kudamono],[:ringo :is :aka],[:mikan :is :yellow]]) == []

 end

 @testset "3 3 vars Test" begin
  @test solveQ([:x,:y,:z],[[:x :y],[:y :z],[:z :x]], [[:a :b],[:b :c],[:c :a]]) == [[:a,:b,:c],[:b,:c,:a],[:c,:a,:b]]
 end

kb = [[:H :は :魔法少女],[:M :は :魔法少女],[:H :持つ :盾],[:M :持つ :弓],[:弓 :は :武器],[:盾 :は :防具]]
 @testset "outof range" begin
  @test solveQ([:x,:y,:z,:w],[[:x :は :y],[:x :持つ :z],[:z :は :w]],kb) == [[:H,:魔法少女,:盾,:防具],[:M,:魔法少女,:弓,:武器]]

 end

 @testset "END OF testPrimitive" begin end

#end
close(flog)

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
module Normalize

using LinearAlgebra
using Statistics, StatsBase
using Random, Distributions

include("../src/util.jl")
using .Utility

# XXX: return coefficients (?)
function random_polynomial(order, scale)
    a = randn(order+1, order+1)
    return (u,v) -> sum(scale^(p₁+p₂)*a[p₁+1,p₂+1]*(u-.5)^p₁*(v-.5)^p₂ for p₁ ∈ 0:order for p₂ ∈ 0:order if (p₁ + p₂) ≤ order)
end

function rescale(x)
    x = x .- minimum(x,dims=2)
    x = x ./ mean(x,dims=2)

    return x
end

function random_surface(dim, len, order, scale)
    u, v = rand(len), rand(len)
    Λ = [ random_polynomial(order, scale) for d ∈ 1:dim ]
    w = hcat((λ.(u,v) for λ ∈ Λ)...)'
    w = rescale(w)

    return w, u, v
end

function random_mean(w, ngene)
    ncell = size(w,2)

    g = rand(Gamma(5,5), ngene) .+ 1
    c = rand(Gamma(5,5), ncell) .+ 1
    μ = Diagonal(g)*(rand(ngene,size(w,1))/size(w,1))*w*Diagonal(c)

    return μ, (gene=g, cell=c)
end

function random_sample(μ)
    ngene, ncell = size(μ)

    ϕ = rand(Gamma(10,10), ngene) .+ 1
    sample = Array{Int,2}(undef, ngene, ncell)
    for i in 1:ngene
        r = ϕ[i]
        for j in 1:ncell
            p = ϕ[i] / (μ[i,j] .+ ϕ[i])
            sample[i,j] = rand(NegativeBinomial(r,p))
        end
    end

    return sample, ϕ
end

function generate(ngene, ncell, rank; order=10)
    w, u, v   = random_surface(rank, ncell, order, 2)
    μ, scale  = random_mean(w, ngene)
    sample, ϕ = random_sample(μ)

    return (
        sample = sample,
        mean   = μ,
        ϕ      = ϕ,
        scale  = scale,
        w      = w,
        u      = u,
        v      = v,
    )
end

function main()
    model = generate(8000,2000,25; order=15)
    σ² = model.sample.*(model.sample .+ model.ϕ) ./ (1 .+ model.ϕ)
    u, v = sinkhorn(σ²)

    N = Diagonal(sqrt.(u))*model.sample*Diagonal(sqrt.(v))
    F = svd(N)

    return model, (norm=N, svd=F), u, v
end

end

if abspath(PROGRAM_FILE) == @__FILE__
    Normalize.main()
end

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
function checkcompilers(app)

cpustatus = Sys.which(app.cpucompiler);
if cpustatus != nothing
    display("Using " * app.cpucompiler * " compiler for CPU source code");
else
    display(app.cpucompiler * " compiler is not available on your system.");
    cpustatus = Sys.which("g++");
    if cpustatus != nothing
        display("However, g++ compiler is available on your system. It will be used for compiling CPU source code.");
        app.cpucompiler = "g++";
    else
        error("C++ compiler is not available on your system.");
    end
end

if app.mpiprocs>1
    mpistatus = Sys.which(app.mpicompiler);
    if mpistatus != nothing
        display("Using " * app.mpicompiler * " compiler for MPI source code");
    else
        display(app.mpicompiler * " compiler is not available on your system.");
        mpistatus = Sys.which("mpicxx");
        if mpistatus != nothing
            display("However, mpicxx compiler is available on your system. It will be used for compiling MPI source code.");
            app.mpicompiler = "mpicxx";
            app.mpirun = "mpirun";
        else
            error("MPI compiler is not available on your system.");
        end
    end
end

if app.platform == "gpu"
    gpustatus = Sys.which(app.gpucompiler);
    if gpustatus != nothing
        display("Using " * app.gpucompiler * " compiler for GPU source code");
    else
        display(app.gpucompiler * " compiler is not available on your system.");
        gpustatus = Sys.which("nvcc");
        if gpustatus != nothing
            display("However, nvcc compiler is available on your system. It will be used for compiling GPU source code.");
            app.gpucompiler = "nvcc";
        else
            error("GPU compiler is not available on your system.");
        end
    end
end

mdir = pwd();
ii = findlast(app.codename, mdir);
coredir = mdir[1:ii[end]] * "/lib";

cpulib = 0;
gpulib = 0;
if Sys.isapple()
    if (isfile(coredir * "/Mac/commonCore.a")) && (isfile(coredir * "/Mac/opuCore.a"))
        cpulib = 1;
    end
    if (isfile(coredir * "/Mac/gpuCore.a"))
        gpulib = 1;
    end
elseif Sys.isunix()
    if (isfile(coredir * "/Linux/commonCore.a")) && (isfile(coredir * "/Linux/opuCore.a"))
        cpulib = 1;
    end
    if (isfile(coredir * "/Linux/gpuCore.a"))
        gpulib = 1;
    end
elseif Sys.iswindows()
    if (isfile(coredir * "/Windows/commonCore.a")) && (isfile(coredir * "/Windows/opuCore.a"))
        cpulib = 1;
    end
    if (isfile(coredir * "/Windows/gpuCore.a"))
        gpulib = 1;
    end
end

if cpulib==0
    display("Generating CPU core libraries.");
    genlib(app.cpucompiler, "", coredir);
end

if gpulib==0
    if (length(app.gpucompiler)>0) && (app.platform == "gpu")
        display("Generating GPU core library.");
        genlib("", app.gpucompiler, coredir);
    end
end

return cpulib, gpulib

end


text:
import Base: +, isoperator

nodes(x) = 1
function nodes(x::Expr)
  n = 1
  for y in x.args
    n += nodes(y)
  end
  return n
end

struct Insert
  idx::Vector{Int}
  tree
end

cost(i::Insert) = 1

function Base.show(io::IO, i::Insert)
  print(io, "insert ")
  show(io, i.tree)
end

struct Delete
  idx::Vector{Int}
  tree
end

cost(d::Delete) = nodes(d.tree)

function Base.show(io::IO, d::Delete)
  print(io, "delete ")
  show(io, d.tree)
end

struct Replace
  idx::Vector{Int}
  old
  new
end

cost(r::Replace) = nodes(r.old)

function Base.show(io::IO, r::Replace)
  print(io, "replace ")
  show(io, r.old)
  print(io, " => ")
  show(io, r.new)
end

struct Patch
  cost::Int
  ps::Vector{Any}
end

Patch(ps) = Patch(isempty(ps) ? 0 : sum(cost.(ps)), ps)

@forward Patch.ps Base.isempty, Base.length

a::Patch + b::Patch = Patch(a.cost + b.cost, vcat(a.ps, b.ps))

function Base.show(io::IO, p::Patch)
  println(io, "Expression patch:")
  for p in p.ps
    println(io, p)
  end
end

best(a::Patch) = a
best(a::Patch, b::Patch, c...) = best(a.cost < b.cost ? a : b, c...)

isopchange(x1, x2) =
  isexpr(x1, :call) &&
  x1.args[1] isa Symbol &&
  x2.args[1] isa Symbol &&
  x1.args[1] != x2.args[1] &&
  (isoperator(x1.args[1]) || isoperator(x2.args[1]))

shouldreplace(x1, x2) = x1 != x2

shouldreplace(x1::Expr, x2::Expr) = x1.head != x2.head || isopchange(x1, x2)

canshift(head, i) = !(head == :call && i == 1)

function nleft(xs)
  h = zeros(Int, length(xs) + 1)
  for i = length(xs):-1:1
    h[i] = h[i+1] + nodes(xs[i])
  end
  return h
end

htable(x, y) = Broadcast.broadcasted((x, y) -> x > y ? x-y : y > x ? 1 : 0, nleft(x), nleft(y)')

function fdiff(head, ii, f1, f2)
  m, n = length(f1)+1, length(f2)+1
  ps = Dict{Tuple{Int,Int},Patch}()
  q = PriorityQueue{Tuple{Int,Int},Int}()
  complete = falses(m, n)
  h = htable(f1, f2)
  function visit!(i, j, p)
    haskey(ps, (i, j)) && (p.cost >= ps[(i,j)].cost) && return
    ps[(i,j)] = p
    q[(i,j)] = p.cost + h[i,j]
  end
  visit!(1, 1, Patch([]))
  while !isempty(q)
    (i, j) = dequeue!(q)
    (i,j) == (m,n) && return ps[(i,j)]
    complete[i,j] = true
    p = ps[(i, j)]
    if canshift(head, i)
      i < m && !complete[i+1,j] && visit!(i+1, j, p + Patch([Delete([ii...,i], f1[i])]))
      j < n && !complete[i,j+1] && visit!(i, j+1, p + Patch([Insert([ii...,i], f2[j])]))
    end
    i < m && j < n && !complete[i+1,j+1] && visit!(i+1, j+1, p + diff([ii...,i],f1[i],f2[j]))
  end
end

diff(ii, x1, x2) =
  (x1,x2) isa NTuple{2,Union{LineNumberNode,Nothing}} ? Patch([]) :
  shouldreplace(x1, x2)  ? Patch([Replace(ii, x1, x2)]) :
  x1 == x2               ? Patch([]) :
  fdiff(x1.head, ii, x1.args, x2.args)

diff(a, b) = diff([], a, b)


text:
export getVariableOrder

"""
    $SIGNATURES

Get the frontal variable IDs `::Int` for a given clique in a Bayes (Junction) tree.
"""
getCliqFrontalVarIds(cliqdata::Union{BayesTreeNodeData, PackedBayesTreeNodeData})::Vector{Symbol} = cliqdata.frontalIDs
getCliqFrontalVarIds(cliq::TreeClique)::Vector{Symbol} = getCliqFrontalVarIds(getData(cliq))

"""
    $SIGNATURES

Get the frontal variable IDs `::Int` for a given clique in a Bayes (Junction) tree.
"""
getFrontals(cliqd::Union{TreeClique,BayesTreeNodeData})::Vector{Symbol} = getCliqFrontalVarIds(cliqd)

"""
    $SIGNATURES

Create a new clique.
"""
function addClique!(bt::BayesTree, dfg::G, varID::Symbol, condIDs::Array{Symbol}=Symbol[])::ExVertex where G <: AbstractDFG
  bt.btid += 1
  clq = Graphs.add_vertex!(bt.bt, ExVertex(bt.btid,string("Clique",bt.btid)))
  bt.cliques[bt.btid] = clq

  setLabel(clq, "")

  # Specific data container
  setData!(clq, emptyBTNodeData())

  appendClique!(bt, bt.btid, dfg, varID, condIDs)
  return clq
end

"""
    $SIGNATURES

Generate the label for particular clique (used by graphviz for visualization).
"""
function makeCliqueLabel(dfg::G, bt::BayesTree, clqID::Int)::String where G <: AbstractDFG
  clq = bt.cliques[clqID]
  flbl = ""
  clbl = ""
  for fr in getData(clq).frontalIDs
    flbl = string(flbl, DFG.getVariable(dfg,fr).label, ",")
  end
  for sepr in getData(clq).separatorIDs
    clbl = string(clbl, DFG.getVariable(dfg,sepr).label, ",")
  end
  setLabel(clq, string(flbl, ": ", clbl))
end

"""
    $SIGNATURES

Add the separator for the newly created clique.
"""
function appendSeparatorToClique(bt::BayesTree, clqID::Int, seprIDs::Array{Symbol,1})
  union!(getData(bt.cliques[clqID]).separatorIDs, seprIDs)
  nothing
end

"""
    $SIGNATURES

Add a new frontal variable to clique.

DevNotes
- TODO, define what "conditionals" are CLEARLY!!
"""
function appendClique!(bt::BayesTree,
                       clqID::Int,
                       dfg::AbstractDFG,
                       varID::Symbol,
                       seprIDs::Array{Symbol,1}=Symbol[] )::Nothing
  #
  clq = bt.cliques[clqID]
  var = DFG.getVariable(dfg, varID)

  # add frontal variable
  push!(getData(clq).frontalIDs, varID)

  # total dictionary of frontals for easy access
  bt.frontals[varID] = clqID

  # TODO - confirm this, append to cliq separator ??
  # @info "going for: appendSeparatorToClique on (clqID, seprIDs)=($clqID, $seprIDs)"
  appendSeparatorToClique(bt, clqID, seprIDs)
  makeCliqueLabel(dfg, bt, clqID)
  return nothing
end

"""
    $SIGNATURES

Instantiate a new child clique in the tree.
"""
function newChildClique!(bt::BayesTree, dfg::AbstractDFG, CpID::Int, varID::Symbol, Sepj::Array{Symbol,1})
  # physically create the new clique
  chclq = addClique!(bt, dfg, varID, Sepj)
  parent = bt.cliques[CpID]
  # Staying with Graphs.jl for tree in first stage
  edge = Graphs.make_edge(bt.bt, parent, chclq)
  Graphs.add_edge!(bt.bt, edge)

  return chclq
end

"""
    $SIGNATURES

Post order tree traversal and build potential functions
"""
function findCliqueFromFrontal(bt::BayesTree, frtlID::Int)
  for cliqPair in bt.cliques
    id = cliqPair[1]
    cliq = cliqPair[2]
    for frtl in getFrontals(cliq)
      if frtl == frtlID
        return cliq
      end
    end
  end
  error("Clique with desired frontal ID not found")
end

"""
    $SIGNATURES

Find parent clique Cp that containts the first eliminated variable of Sj as frontal.
"""
function identifyFirstEliminatedSeparator(dfg::AbstractDFG,
                                          elimorder::Vector{Symbol},
                                          firvert::DFGVariable,
                                          Sj=solverData(firvert).separator)::DFGVariable
  #
  firstelim = 99999999999
  for s in Sj
    temp = something(findfirst(isequal(s), elimorder), 0) # findfirst(p, s)
    if (temp < firstelim)
      firstelim = temp
    end
  end
  DFG.getVariable(dfg, elimorder[firstelim])
end

"""
    $SIGNATURES

Eliminate a variable and add to tree cliques accordingly.

Dev Notes
- `p` should be elimination order.
- `var` is next variable to be added to the tree.
- TODO, make sure this works for disjoint graphs
  - Check laplacian, check eigen == 1 is disjoint sets

References
Kaess et al.: Bayes Tree, WAFR, 2010, [Alg. 2]
Kaess et al.: iSAM2, IJRR, 2011, [Alg. 3]
Fourie, D.: mmisam, PhD thesis, 2017. [Chpt. 5]
"""
function newPotential(tree::BayesTree, dfg::G, var::Symbol, elimorder::Array{Symbol,1}) where G <: AbstractDFG
    firvert = DFG.getVariable(dfg,var)
    # no parent
    if (length(solverData(firvert).separator) == 0)
      # if (length(tree.cliques) == 0)
        # create new root
        addClique!(tree, dfg, var)
      # else
      #   # add to root
      #   @warn "root append clique is happening"
      #   appendClique!(tree, 1, dfg, var)
      # end
    else
      # find parent clique Cp that containts the first eliminated variable of Sj as frontal
      Sj = solverData(firvert).separator
      felbl = identifyFirstEliminatedSeparator(dfg, elimorder, firvert, Sj).label
      # get clique id of first eliminated frontal
      CpID = tree.frontals[felbl]
      # look to add this conditional to the tree
      cliq = tree.cliques[CpID]
      # clique of the first eliminated frontal
      unFC = union(getCliqFrontalVarIds(cliq), getCliqSeparatorVarIds(cliq))
      # if the separator of this new variable is identical to the (entire) clique of the firstly eliminated frontal.
      if (sort(unFC) == sort(Sj))
        # just add new variable as frontal to this clique
        # insert conditional (p(var|sepr)) into clique CpID -- i.e. just adding a frontal
        # @info "adding new frontal $var to existing clique $CpID"
        appendClique!(tree, CpID, dfg, var)
      else
        # a new child clique is required here (this becomes parent)
        # @info "adding new child clique with parent separator."
        newChildClique!(tree, dfg, CpID, var, Sj)
      end
    end
end

"""
    $SIGNATURES

Build the whole tree in batch format.
"""
function buildTree!(tree::BayesTree, dfg::AbstractDFG, elimorder::Array{Symbol,1})
  revorder = reverse(elimorder,dims=1) # flipdim(p, 1), fixing #499
  # prevVar = 0
  for var in revorder
    @info "Adding $var to tree..."
    newPotential(tree, dfg, var, elimorder)
    prevVar = var
  end

  return tree
end

"""
    $SIGNATURES

Open view to see the graphviz exported Bayes tree, assuming default location and
viewer app.  See keyword arguments for more details.
"""
function showTree(;filepath::String="/tmp/caesar/bt.pdf",
                   viewerapp::String="evince"  )
  #
  try
    @async run(`$(viewerapp) $(filepath)`)
  catch ex
    @warn "not able to show via $(viewerapp) $(filepath)"
    @show ex
    @show stacktrace()
  end
end

"""
    $SIGNATURES

Draw the Bayes (Junction) tree by means of `.dot` files and `pdf` reader.

Notes
- Uses system install of graphviz.org.
- Can also use Linux tool `xdot`.
"""
function drawTree(treel::BayesTree;
                  show::Bool=false,                  # must remain false for stability and automated use in solver
                  filepath::String="/tmp/caesar/bt.pdf",
                  viewerapp::String="evince",
                  imgs::Bool=false )
  #
  fext = filepath[(end-2):end] #split(filepath, '.')[end]
  fpwoext = filepath[1:(end-4)]# split(filepath, '.')[end-1]
  mkpath(dirname(fpwoext))

  # modify a deepcopy
  btc = deepcopy(treel)
  for (cid, cliq) in btc.cliques
    if imgs
      firstlabel = split(getLabel(cliq),',')[1]
      spyCliqMat(cliq, suppressprint=true) |> exportimg("/tmp/$firstlabel.png")
      cliq.attributes["image"] = "/tmp/$firstlabel.png"
      setLabel(cliq, "")
    end
    delete!(cliq.attributes, "data")
  end

  fid = IOStream("")
  try
    fid = open("$(fpwoext).dot","w+")
    write(fid,to_dot(btc.bt))
    close(fid)
    run(`dot $(fpwoext).dot -T $(fext) -o $(filepath)`)
  catch ex
    @warn ex
    @show stacktrace()
  finally
    close(fid)
  end

  show ? showTree(viewerapp=viewerapp, filepath=filepath) : nothing
end

"""
    $SIGNATURES

Draw the Bayes (junction) tree with LaTeX labels by means of `.dot` and `.tex`
files.

Notes
- Uses system install of graphviz.org.
- Uses external python `dot2tex` tool (`pip install dot2tex`).

Related:

drawTree
"""
function generateTexTree(treel::BayesTree;
                         filepath::String="/tmp/caesar/bayes/bt")
    btc = deepcopy(treel)
    for (cid, cliq) in btc.cliques
        label = getLabel(cliq)

        # Get frontals and separator, and split into elements.
        frt, sep = split(label,':')
        efrt = split(frt,',')
        esep = split(sep,',')

        # Transform frontals into latex.
        newfrontals = ""
        for l in efrt
            # Split into symbol and subindex (letter and number).
            parts = split(l, r"[^a-z0-9]+|(?<=[a-z])(?=[0-9])|(?<=[0-9])(?=[a-z])")
            if size(parts)[1] == 2
                newfrontals = string(newfrontals, "\\bm{", parts[1], "}_{", parts[2], "}, ")
            elseif size(parts)[1] == 3
                newfrontals = string(newfrontals, "\\bm{", parts[2], "}_{", parts[3], "}, ")
            end
        end
        # Remove the trailing comma.
        newfrontals = newfrontals[1:end-2]

        # Transform separator into latex.
        newseparator = ""
        if length(sep) > 1
            for l in esep
                # Split into symbol and subindex.
                parts = split(l, r"[^a-z0-9]+|(?<=[a-z])(?=[0-9])|(?<=[0-9])(?=[a-z])")
                if size(parts)[1] == 2
                    newseparator = string(newseparator, "\\bm{", parts[1], "}_{", parts[2], "}, ")
                elseif size(parts)[1] == 3
                    newseparator = string(newseparator, "\\bm{", parts[2], "}_{", parts[3], "}, ")
                end
            end
        end
        # Remove the trailing comma.
        newseparator = newseparator[1:end-2]
        # Create full label and replace the old one.
        newlabel = string(newfrontals, ":", newseparator)
        setLabel(cliq, newlabel)
    end

    # Use new labels to produce `.dot` and `.tex` files.
    fid = IOStream("")
    try
        mkpath(joinpath((split(filepath,'/')[1:(end-1)])...) )
        fid = open("$(filepath).dot","w+")
        write(fid, to_dot(btc.bt))
        close(fid)
        # All in one command.
        run(`dot2tex -tmath --preproc $(filepath).dot -o $(filepath)proc.dot`)
        run(`dot2tex $(filepath)proc.dot -o $(filepath).tex`)
    catch ex
        @warn ex
        @show stacktrace()
    finally
        close(fid)
    end

    return btc
end

"""
    $SIGNATURES

Build Bayes/Junction/Elimination tree from a given variable ordering.
"""
function buildTreeFromOrdering!(dfg::G,
                                p::Vector{Symbol};
                                drawbayesnet::Bool=false,
                                maxparallel::Int=50,
                                solvable::Int=1 ) where G <: InMemoryDFGTypes
  #
  t0 =time_ns()
  println()
  fge = deepcopy(dfg)
  @info "Building Bayes net..."
  buildBayesNet!(fge, p, maxparallel=maxparallel, solvable=solvable)

  @info "Staring the Bayes tree construction from Bayes net"
  tree = emptyBayesTree()
  tree.variableOrder = p
  buildTree!(tree, fge, p)

  if drawbayesnet
    println("Bayes Net")
    sleep(0.1)
    fid = open("bn.dot","w+")
    write(fid,to_dot(fge.bn))
    close(fid)
  end

  println("Find potential functions for each clique")
  cliq = tree.cliques[1] # start at the root
  buildCliquePotentials(dfg, tree, cliq, solvable=solvable); # fg does not have the marginals as fge does

  tree.buildTime = (time_ns()-t0)/1e9
  return tree
end

"""
    $SIGNATURES

Build Bayes/Junction/Elimination tree from a given variable ordering.
"""
function buildTreeFromOrdering!(dfg::DFG.AbstractDFG,
                                p::Vector{Symbol};
                                drawbayesnet::Bool=false,
                                maxparallel::Int=50  )
  #
  println()

  @info "Copying to a local DFG"
  fge = InMemDFGType(params=getSolverParams(dfg))#GraphsDFG{SolverParams}(params=SolverParams())
  DistributedFactorGraphs._copyIntoGraph!(dfg, fge, union(getVariableIds(dfg), getFactorIds(dfg)), true)

  println("Building Bayes net from cloud...")
  buildBayesNet!(fge, p, maxparallel=maxparallel)

  tree = emptyBayesTree()
  tree.variableOrder = p
  buildTree!(tree, fge, p)

  if drawbayesnet
    println("Bayes Net")
    sleep(0.1)
    fid = open("bn.dot","w+")
    write(fid,to_dot(fge.bn))
    close(fid)
  end

  println("Find potential functions for each clique")
  cliq = tree.cliques[1] # start at the root
  buildCliquePotentials(dfg, tree, cliq); # fg does not have the marginals as fge does

  return tree
end

"""
    $SIGNATURES

Build Bayes/Junction/Elimination tree.

Notes
- Default to free qr factorization for variable elimination order.
"""
function prepBatchTree!(dfg::AbstractDFG;
                        ordering::Symbol=:qr,
                        variableOrder::Union{Nothing, Vector{Symbol}}=nothing,
                        drawpdf::Bool=false,
                        show::Bool=false,
                        filepath::String="/tmp/caesar/bt.pdf",
                        viewerapp::String="evince",
                        imgs::Bool=false,
                        drawbayesnet::Bool=false,
                        maxparallel::Int=50 )
  #
  p = variableOrder != nothing ? variableOrder : getEliminationOrder(dfg, ordering=ordering)

  # for debuggin , its useful to have the variable ordering
  if drawpdf
    ispath(getLogPath(dfg)) ? nothing : Base.mkpath(getLogPath(dfg))
    open(joinLogPath(dfg,"variableOrder.txt"), "a") do io
      writedlm(io, string.(reshape(p,1,:)), ',')
    end
  end

  tree = buildTreeFromOrdering!(dfg, p, drawbayesnet=drawbayesnet, maxparallel=maxparallel)

  # GraphViz.Graph(to_dot(tree.bt))
  # Michael reference -- x2->x1, x2->x3, x2->x4, x2->l1, x4->x3, l1->x3, l1->x4
  #Michael reference 3sig -- x2l1x4x3    x1|x2
  println("Bayes Tree")
  if drawpdf
    drawTree(tree, show=show, filepath=filepath, viewerapp=viewerapp, imgs=imgs)
  end

  return tree
end

"""
    $SIGNATURES

Partial reset of basic data fields in `::VariableNodeData` of `::FunctionNode` structures.
"""
function resetData!(vdata::VariableNodeData)::Nothing
  vdata.eliminated = false
  vdata.BayesNetOutVertIDs = Symbol[]
  # vdata.BayesNetVertID = :_null # TODO dont use nothing, see DFG issue #16
  vdata.separator = Symbol[]
  nothing
end

function resetData!(vdata::FunctionNodeData)::Nothing
  vdata.eliminated = false
  vdata.potentialused = false
  nothing
end

"""
    $SIGNATURES

Wipe data from `dfg` object so that a completely fresh Bayes/Junction/Elimination tree
can be constructed.
"""
function resetFactorGraphNewTree!(dfg::AbstractDFG)::Nothing
  for v in DFG.getVariables(dfg)
    resetData!(solverData(v))
  end
  for f in DFG.getFactors(dfg)
    resetData!(solverData(f))
  end
  nothing
end

"""
    $SIGNATURES

Reset factor graph and build a new tree from the provided variable ordering `p`.
"""
function resetBuildTreeFromOrder!(fgl::AbstractDFG, p::Vector{Symbol})::BayesTree
  resetFactorGraphNewTree!(fgl)
  return buildTreeFromOrdering!(fgl, p)
end

"""
    $(SIGNATURES)

Build a completely new Bayes (Junction) tree, after first wiping clean all
temporary state in fg from a possibly pre-existing tree.

Related:

buildTreeFromOrdering!
"""
function wipeBuildNewTree!(dfg::G;
                           ordering::Symbol=:qr,
                           drawpdf::Bool=false,
                           show::Bool=false,
                           filepath::String="/tmp/caesar/bt.pdf",
                           viewerapp::String="evince",
                           imgs::Bool=false,
                           maxparallel::Int=50,
                           variableOrder::Union{Nothing, Vector{Symbol}}=nothing  )::BayesTree where G <: AbstractDFG
  #
  resetFactorGraphNewTree!(dfg);
  return prepBatchTree!(dfg, variableOrder=variableOrder, ordering=ordering, drawpdf=drawpdf, show=show, filepath=filepath, viewerapp=viewerapp, imgs=imgs, maxparallel=maxparallel);
end

"""
    $(SIGNATURES)

Return the TreeClique node object that represents a clique in the Bayes
(Junction) tree, as defined by one of the frontal variables `frt<:AbstractString`.

Notes
- Frontal variables only occur once in a clique per tree, therefore is a unique identifier.

Related:

getCliq, getTreeAllFrontalSyms
"""
getCliq(bt::BayesTree, frt::Symbol) = bt.cliques[bt.frontals[frt]]

"""
    $(SIGNATURES)

Return the TreeClique node object that represents a clique in the Bayes
(Junction) tree, as defined by one of the frontal variables `frt<:AbstractString`.

Notes
- Frontal variables only occur once in a clique per tree, therefore is a unique identifier.

Related:

getCliq, getTreeAllFrontalSyms
"""
whichCliq(bt::BayesTree, frt::Symbol) = getCliq(bt, frt)
whichCliq(bt::BayesTree, frt::T) where {T <: AbstractString} = whichCliq(bt, Symbol(string(frt)))

"""
    $SIGNATURES

Return boolean on whether the frontal variable `frt::Symbol` exists somewhere in the `::BayesTree`.
"""
hasCliq(bt::BayesTree, frt::Symbol)::Bool = haskey(bt.frontals, frt)

"""
    $SIGNATURES

Return depth in tree as `::Int`, with root as depth=0.

Related

getCliq
"""
function getCliqDepth(tree, cliq)::Int
  prnt = getParent(tree, cliq)
  if length(prnt) == 0
    return 0
  end
  return getCliqDepth(tree, prnt[1]) + 1
end
getCliqDepth(tree::BayesTree, sym::Symbol)::Int = getCliqDepth(tree, getCliq(tree, sym))

"""
    $(SIGNATURES)

Set the upward passing message for Bayes (Junction) tree clique `cliql`.

Dev Notes
- TODO setUpMsg! should also set inferred dimension
"""
function setUpMsg!(cliql::ExVertex, msgs::Dict{Symbol, BallTreeDensity})
  @error "setUpMsg!, use inferred dimension version instead"
  getData(cliql).upMsg = msgs
end

function setUpMsg!(cliql::ExVertex, msgs::TempBeliefMsg) #Dict{Symbol, Tuple{BallTreeDensity, Float64}}
  # ms = Dict{Symbol, BallTreeDensity}()
  # for (id, val) in msgs
  #   ms[id] = val[1]
  # end
  getData(cliql).upMsg = msgs #ms
  nothing
end

"""
    $(SIGNATURES)

Return the last up message stored in `cliq` of Bayes (Junction) tree.
"""
getUpMsgs(cliql::TreeClique) = getData(cliql).upMsg
getUpMsgs(btl::BayesTree, sym::Symbol) = getUpMsgs(whichCliq(btl, sym))

"""
    $(SIGNATURES)

Return the last up message stored in `cliq` of Bayes (Junction) tree.
"""
getCliqMsgsUp(cliql::TreeClique) = upMsg(cliql)
getCliqMsgsUp(treel::BayesTree, frt::Symbol) = getCliqMsgsUp(getCliq(treel, frt))

"""
    $(SIGNATURES)

Set the downward passing message for Bayes (Junction) tree clique `cliql`.
"""
function setDwnMsg!(cliql::ExVertex, msgs::TempBeliefMsg) #Dict{Symbol, BallTreeDensity}
  getData(cliql).dwnMsg = msgs
end

"""
    $(SIGNATURES)

Return the last down message stored in `cliq` of Bayes (Junction) tree.
"""
getDwnMsgs(cliql::TreeClique) = getData(cliql).dwnMsg
getDwnMsgs(btl::BayesTree, sym::Symbol) = getDwnMsgs(whichCliq(btl, sym))

"""
    $(SIGNATURES)

Return the last down message stored in `cliq` of Bayes (Junction) tree.
"""
getCliqMsgsDown(cliql::TreeClique) = getDwnMsgs(cliql)

function appendUseFcts!(usefcts,
                        lblid::Symbol,
                        fct::DFGFactor )
                        # fid::Symbol )
  #
  union!(usefcts, Symbol(fct.label))
  # for tp in usefcts
  #   if tp == fct.index
  #     return nothing
  #   end
  # end
  # tpl = fct.label
  # push!(usefcts, tpl )
  nothing
end

"""
    $SIGNATURES

Return list of factors which depend only on variables in variable list in factor
graph -- i.e. among variables.

Notes
-----
* `unused::Bool=true` will disregard factors already used -- i.e. disregard where `potentialused=true`
"""
function getFactorsAmongVariablesOnly(dfg::G,
                                      varlist::Vector{Symbol};
                                      unused::Bool=true  ) where G <: AbstractDFG
  # collect all factors attached to variables
  prefcts = Symbol[]
  for var in varlist
    union!(prefcts, DFG.ls(dfg, var))
  end

  almostfcts = Symbol[]
  if unused
    # now check if those factors have already been added
    for fct in prefcts
      vert = DFG.getFactor(dfg, fct)
      if !solverData(vert).potentialused
        push!(almostfcts, fct)
      end
    end
  else
    almostfcts = prefcts
  end

  # Select factors that have all variables in this clique var list
  usefcts = Symbol[]
  for fct in almostfcts
    if length(setdiff(DFG.getNeighbors(dfg, fct), varlist)) == 0
      push!(usefcts, fct)
    end
  end

  return usefcts
end

"""
    $SIGNATURES

Get all factors connected to frontal variables

Dev Notes
- Why not just do this, `ffcs = union( map(x->ls(fgl, x), frtl) )`?
"""
function getCliqFactorsFromFrontals(fgl::G,
                                    cliq::TreeClique,
                                    varlist::Vector{Symbol};
                                    inseparator::Bool=true,
                                    unused::Bool=true,
                                    solvable::Int=1  ) where {G <: AbstractDFG}
  #
  frtls = getData(cliq).frontalIDs
  seprs = getData(cliq).separatorIDs
  allids = [frtls;seprs]

  usefcts = Symbol[]
  for frsym in frtls
    # usefcts = Int[]
    for fctid in ls(fgl, frsym)
        fct = getFactor(fgl, fctid)
        if !unused || !solverData(fct).potentialused
            loutn = ls(fgl, fctid, solvable=solvable)
            # deal with unary factors
            if length(loutn)==1
                union!(usefcts, Symbol[Symbol(fct.label);])
                # appendUseFcts!(usefcts, loutn, fct) # , frsym)
                solverData(fct).potentialused = true
            end
            # deal with n-ary factors
            for sep in loutn
                if sep == frsym
                    continue # skip the frsym itself
                end
                insep = sep in allids
                if !inseparator || insep
                    union!(usefcts, Symbol[Symbol(fct.label);])
                    solverData(fct).potentialused = true
                    if !insep
                        @info "cliq=$(cliq.index) adding factor that is not in separator, $sep"
                    end
                end
            end
        end
    end
  end

  return usefcts
end

"""
    $SIGNATURES

Return `::Bool` on whether factor is a partial constraint.
"""
isPartial(fcf::T) where {T <: FunctorInferenceType} = :partial in fieldnames(T)
function isPartial(fct::DFGFactor)  #fct::TreeClique
  fcf = solverData(fct).fnc.usrfnc!
  isPartial(fcf)
end

"""
    $SIGNATURES

Determine and set the potentials for a particular `cliq` in the Bayes (Junction) tree.
"""
function setCliqPotentials!(dfg::G,
                            bt::BayesTree,
                            cliq::TreeClique;
                            solvable::Int=1  ) where G <: AbstractDFG
  #
  varlist = getCliqVarIdsAll(cliq)

  @info "using all factors connected to frontals and attached to separator"
  fctsyms = getFactorsAmongVariablesOnly(dfg, varlist, unused=true )
  # filter only factors connected to frontals (for upward)
  frtfcts = union(map(x->ls(dfg, x), getCliqFrontalVarIds(cliq))...)
  fctsyms = intersect(fctsyms, frtfcts)

  getData(cliq).potentials = fctsyms
  getData(cliq).partialpotential = Vector{Bool}(undef, length(fctsyms))
  fcts = map(x->getFactor(dfg, x), fctsyms)
  getData(cliq).partialpotential = map(x->isPartial(x), fcts)
  for fct in fcts
    solverData(fct).potentialused = true
  end

  @info "finding all frontals for down WIP"
  ffctsyms = getCliqFactorsFromFrontals(dfg, cliq, Symbol[], inseparator=false, unused=false, solvable=solvable)
  # fnsyms = getCliqVarsWithFrontalNeighbors(csmc.dfg, csmc.cliq)
  getData(cliq).dwnPotentials = ffctsyms
  getData(cliq).dwnPartialPotential = map(x->isPartial(getFactor(dfg,x)), ffctsyms )

  nothing
end

function getCliquePotentials!(dfg::G,
                            bt::BayesTree,
                            cliq::TreeClique  ) where G <: AbstractDFG
  #
  @warn "getCliquePotentials! deprecated, use getCliqPotentials instead."
  getCliqPotentials(dfg, bt, cliq)
end

function cliqPotentialIDs(cliq::TreeClique)
  potIDs = Symbol[]
  for idfct in getData(cliq).potentials
    push!(potIDs, idfct)
  end
  return potIDs
end

"""
    $SIGNATURES

Collect and returl all child clique separator variables.
"""
function collectSeparators(bt::BayesTree, cliq::TreeClique)::Vector{Symbol}
  allseps = Symbol[]
  for child in out_neighbors(cliq, bt.bt)#tree
      allseps = [allseps; getData(child).separatorIDs]
  end
  return allseps
end

"""
    $SIGNATURES

Return boolean matrix of factor by variable (row by column) associations within
clique, corresponds to order presented by `getCliqFactorIds` and `getCliqAllVarIds`.
"""
function getCliqAssocMat(cliq::TreeClique)
  getData(cliq).cliqAssocMat
end

"""
    $SIGNATURES

Return boolean matrix of upward message singletons (i.e. marginal priors) from
child cliques.  Variable order corresponds to `getCliqAllVarIds`.
"""
getCliqMsgMat(cliq::TreeClique) = getData(cliq).cliqMsgMat

"""
    $SIGNATURES

Return boolean matrix of factor variable associations for a clique, optionally
including (`showmsg::Bool=true`) the upward message singletons.  Variable order
corresponds to `getCliqAllVarIds`.
"""
function getCliqMat(cliq::TreeClique; showmsg::Bool=true)
  assocMat = getCliqAssocMat(cliq)
  msgMat = getCliqMsgMat(cliq)
  mat = showmsg ? [assocMat;msgMat] : assocMat
  return mat
end

"""
    $SIGNATURES

Get `cliq` separator (a.k.a. conditional) variable ids`::Symbol`.
"""
getCliqSeparatorVarIds(cliqdata::Union{BayesTreeNodeData, PackedBayesTreeNodeData})::Vector{Symbol} = cliqdata.separatorIDs
getCliqSeparatorVarIds(cliq::TreeClique)::Vector{Symbol} = getCliqSeparatorVarIds(getData(cliq))

"""
    $SIGNATURES

Get `cliq` potentials (factors) ids`::Int`.
"""
getCliqFactorIds(cliqdata::BayesTreeNodeData)::Vector{Symbol} = cliqdata.potentials
getCliqFactorIds(cliq::TreeClique)::Vector{Symbol} = getCliqFactorIds(getData(cliq))

"""
    $SIGNATURES

Get all `cliq` variable ids`::Symbol`.

Related

getCliqVarIdsAll, getCliqAllFactIds, getCliqVarsWithFrontalNeighbors
"""
function getCliqAllVarIds(cliq::TreeClique)::Vector{Symbol}
  frtl = getCliqFrontalVarIds(cliq)
  cond = getCliqSeparatorVarIds(cliq)
  union(frtl,cond)
end

"""
    $SIGNATURES

Return all variables (including frontal factor neighbors).

Dev Notes
- TODO needs to be refactored and optimized.

Related

getCliqAllVarIds
"""
function getCliqVarsWithFrontalNeighbors(fgl::G,
                                         cliq::TreeClique;
                                         solvable::Int=1  ) where {G <: AbstractDFG}
  #
  frtl = getCliqFrontalVarIds(cliq)
  cond = getCliqSeparatorVarIds(cliq)
  syms = Symbol[]
  union!(syms,Symbol.(frtl))
  union!(syms,Symbol.(cond))

  # TODO Can we trust factors are frontal connected?
  ffcs = union( map(x->ls(fgl, x, solvable=solvable), frtl)... )
  # @show ffcs = getData(cliq).potentials
  neig = union( map(x->ls(fgl, x, solvable=solvable), ffcs)... )
  union!(syms, Symbol.(neig))
  return syms
end

"""
    $SIGNATURES

Get all `cliq` variable ids`::Symbol`.

Related

getCliqAllVarIds, getCliqAllFactIds
"""
getCliqVarIdsAll(cliq::TreeClique)::Vector{Symbol} = getCliqAllVarIds(cliq::TreeClique)

"""
    $SIGNATURES

Get all `cliq` factor ids`::Symbol`.

DEPRECATED, use getCliqFactorIdsAll instead.

Related

getCliqVarIdsAll, getCliqFactors
"""
getCliqFactorIdsAll(cliqd::BayesTreeNodeData) = cliqd.potentials
getCliqFactorIdsAll(cliq::TreeClique) = getCliqFactorIdsAll(getData(cliq))
getCliqFactorIdsAll(treel::BayesTree, frtl::Symbol) = getCliqFactorIdsAll(getCliq(treel, frtl))

const getCliqFactors = getCliqFactorIdsAll

"""
    $SIGNATURES

Get all `cliq` factor ids`::Symbol`.

DEPRECATED, use getCliqFactorIdsAll instead.

Related

getCliqVarIdsAll
"""
function getCliqAllFactIds(cliqd::BayesTreeNodeData)
    @warn "getCliqAllFactIds deprecated, use getCliqFactorIdsAll instead."
    return getCliqFactorIdsAll(cliqd)
end

function getCliqAllFactIds(cliq::TreeClique)
    @warn "getCliqAllFactIds deprecated, use getCliqFactorIdsAll instead."
    return getCliqFactorIdsAll(getData(cliq))
end

"""
    $SIGNATURES

Get all `cliq` variable labels as `::Symbol`.
"""
function getCliqAllVarSyms(dfg::G, cliq::TreeClique)::Vector{Symbol} where G <: AbstractDFG
  # Symbol[getSym(dfg, varid) for varid in getCliqAllVarIds(cliq)]
  @warn "deprecated getCliqAllVarSyms, use getCliqAllVarIds instead."
  getCliqAllVarIds(cliq) # not doing all frontals
end

"""
    $SIGNATURES

Return dictionary of down messages consisting of all frontal and separator beliefs of this clique.

Notes:
- Fetches numerical results from `subdfg` as dictated in `cliq`.
"""
function getCliqDownMsgsAfterDownSolve(subdfg::G, cliq::TreeClique)::TempBeliefMsg where G <: AbstractDFG
  # Dict{Symbol, BallTreeDensity}
  # where the return msgs are contained
  container = TempBeliefMsg() # Dict{Symbol,BallTreeDensity}()

  # go through all msgs one by one
  for sym in getCliqAllVarIds(cliq)
    container[sym] = (getKDE(subdfg, sym), getVariableInferredDim(subdfg, sym))
  end

  # return the result
  return container
end

"""
    $SIGNATURES

Get variable ids`::Int` with prior factors associated with this `cliq`.

Notes:
- does not include any singleton messages from upward or downward message passing.
"""
function getCliqVarIdsPriors(cliq::TreeClique,
                             allids::Vector{Symbol}=getCliqAllVarIds(cliq),
                             partials::Bool=true  )::Vector{Symbol}
  # get ids with prior factors associated with this cliq
  amat = getCliqAssocMat(cliq)
  prfcts = sum(amat, dims=2) .== 1

  # remove partials as per request
  !partials ? nothing : (prfcts .&= getData(cliq).partialpotential)

  # return variable ids in `mask`
  mask = sum(amat[prfcts[:],:], dims=1)[:] .> 0
  return allids[mask]
end

"""
    $SIGNATURES

Return the number of factors associated with each variable in `cliq`.
"""
getCliqNumAssocFactorsPerVar(cliq::TreeClique)::Vector{Int} = sum(getCliqAssocMat(cliq), dims=1)[:]

"""
    $SIGNATURES

Get `cliq` variable IDs with singleton factors -- i.e. both in clique priors and up messages.
"""
function getCliqVarSingletons(cliq::TreeClique,
                              allids::Vector{Symbol}=getCliqAllVarIds(cliq),
                              partials::Bool=true  )::Vector{Symbol}
  # get incoming upward messages (known singletons)
  mask = sum(getCliqMsgMat(cliq),dims=1)[:] .>= 1
  upmsgids = allids[mask]

  # get ids with prior factors associated with this cliq
  prids = getCliqVarIdsPriors(cliq, getCliqAllVarIds(cliq), partials)

  # return union of both lists
  return union(upmsgids, prids)
end

"""
    $SIGNATURES

Get each clique subgraph association matrix.
"""
function compCliqAssocMatrices!(dfg::G, bt::BayesTree, cliq::TreeClique) where G <: AbstractDFG
  frtl = getCliqFrontalVarIds(cliq)
  cond = getCliqSeparatorVarIds(cliq)
  inmsgIDs = collectSeparators(bt, cliq)
  potIDs = cliqPotentialIDs(cliq)
  # Construct associations matrix here
  # matrix has variables are columns, and messages/constraints as rows
  cols = [frtl;cond]
  getData(cliq).inmsgIDs = inmsgIDs
  getData(cliq).potIDs = potIDs
  cliqAssocMat = Array{Bool,2}(undef, length(potIDs), length(cols))
  cliqMsgMat = Array{Bool,2}(undef, length(inmsgIDs), length(cols))
  fill!(cliqAssocMat, false)
  fill!(cliqMsgMat, false)
  for j in 1:length(cols)
    for i in 1:length(inmsgIDs)
      if cols[j] == inmsgIDs[i]
        cliqMsgMat[i,j] = true
      end
    end
    for i in 1:length(potIDs)
      idfct = getData(cliq).potentials[i]
      if idfct == potIDs[i] # sanity check on clique potentials ordering
        # TODO int and symbol compare is no good
        for vertidx in solverData(DFG.getFactor(dfg, idfct)).fncargvID
          if vertidx == cols[j]
            cliqAssocMat[i,j] = true
          end
        end
      else
        prtslperr("compCliqAssocMatrices! -- potential ID ordering was lost")
      end
    end
  end
  getData(cliq).cliqAssocMat = cliqAssocMat
  getData(cliq).cliqMsgMat = cliqMsgMat
  nothing
end

function countSkips(bt::BayesTree)
  skps = 0
  for cliq in bt.cliques
    m = getCliqMat(cliq[2])
    mi = map(Int,m)
    skps += sum(map(Int,sum(mi, dims=1) .== 1))
  end
  return skps
end

function skipThroughMsgsIDs(cliq::TreeClique)
  cliqdata = getData(cliq)
  numfrtl1 = floor(Int,length(cliqdata.frontalIDs)+1)
  condAssocMat = cliqdata.cliqAssocMat[:,numfrtl1:end]
  condMsgMat = cliqdata.cliqMsgMat[:,numfrtl1:end]
  mat = [condAssocMat;condMsgMat];
  mab = sum(map(Int,mat),dims=1) .== 1
  mabM = sum(map(Int,condMsgMat),dims=1) .== 1
  mab = mab .& mabM
  # rang = 1:size(condMsgMat,2)
  msgidx = cliqdata.separatorIDs[vec(collect(mab))]
  return msgidx
end

function directPriorMsgIDs(cliq::TreeClique)
  frtl = getData(cliq).frontalIDs
  sepr = getData(cliq).separatorIDs
  cols = [frtl;sepr]
  mat = getCliqMat(cliq, showmsg=true)
  singr = sum(map(Int,mat),dims=2) .== 1
  rerows = collect(1:length(singr))
  b = vec(collect(singr))
  rerows2 = rerows[b]
  sumsrAc = sum(map(Int,mat[rerows2,:]),dims=1)
  sumc = sum(map(Int,mat),dims=1)
  pmSkipCols = (sumsrAc - sumc) .== 0
  return cols[vec(collect(pmSkipCols))]
end

function directFrtlMsgIDs(cliq::TreeClique)
  numfrtl = length(getData(cliq).frontalIDs)
  frntAssocMat = getData(cliq).cliqAssocMat[:,1:numfrtl]
  frtlMsgMat = getData(cliq).cliqMsgMat[:,1:numfrtl]
  mat = [frntAssocMat; frtlMsgMat];
  mab = sum(map(Int,mat),dims=1) .== 1
  mabM = sum(map(Int,frtlMsgMat),dims=1) .== 1
  mab = mab .& mabM
  return getData(cliq).frontalIDs[vec(collect(mab))]
end

function directAssignmentIDs(cliq::TreeClique)
  # NOTE -- old version been included in iterated variable stack
  assocMat = getData(cliq).cliqAssocMat
  msgMat = getData(cliq).cliqMsgMat
  mat = [assocMat;msgMat];
  mab = sum(map(Int,mat),dims=1) .== 1
  mabA = sum(map(Int,assocMat),dims=1) .== 1
  mab = mab .& mabA
  # TODO -- use proper accessor methods
  frtl = getData(cliq).frontalIDs
  sepr = getData(cliq).separatorIDs
  cols = [frtl;sepr]
  return cols[vec(collect(mab))]
  # also calculate how which are conditionals
end

function mcmcIterationIDs(cliq::TreeClique)
  mat = getCliqMat(cliq)
  # assocMat = getData(cliq).cliqAssocMat
  # msgMat = getData(cliq).cliqMsgMat
  # mat = [assocMat;msgMat];

  sum(sum(map(Int,mat),dims=1)) == 0 ? error("mcmcIterationIDs -- unaccounted variables") : nothing
  mab = 1 .< sum(map(Int,mat),dims=1)
  cols = getCliqAllVarIds(cliq)

  # must also include "direct variables" connected through projection only
  directvars = directAssignmentIDs(cliq)
  usset = union(directvars, cols[vec(collect(mab))])
  # NOTE -- fix direct vs itervar issue, DirectVarIDs against Iters should also Iter
  # NOTE -- using direct then mcmcIter ordering to prioritize non-msg vars first
  return setdiff(usset, getData(cliq).directPriorMsgIDs)
end

function getCliqMatVarIdx(cliq::TreeClique, varid::Symbol, allids=getCliqAllVarIds(cliq) )
  len = length(allids)
  [1:len;][allids .== varid][1]
end

"""
    $SIGNATURES

Determine and return order list of variable ids required for minibatch Gibbs iteration inside `cliq`.

Notes
- Singleton factors (priors and up messages) back of the list
- least number of associated factor variables earlier in list
- Same as getCliqVarSolveOrderUp
"""
function mcmcIterationIdsOrdered(cliq::TreeClique)
  # get unordered iter list
  alliter = mcmcIterationIDs(cliq)

  # get all singletons
  allsings = getCliqVarSingletons(cliq)
  singletonvars = intersect(alliter, allsings)

  # get all non-singleton iters
  nonsinglvars = setdiff(alliter, singletonvars)

  # sort nonsingletons ascending number of factors
  mat = getCliqMat(cliq)
  lenfcts = sum(mat, dims=1)
  nonslen = zeros(length(nonsinglvars))
  for i in 1:length(nonsinglvars)
    varid = nonsinglvars[i]
    varidx = getCliqMatVarIdx(cliq, varid)
    nonslen[i] = lenfcts[varidx]
  end
  p = sortperm(nonslen)
  ascnons = nonsinglvars[p]

  # sort singleton vars ascending number of factors
  singslen = zeros(length(singletonvars))
  for i in 1:length(singletonvars)
    varid = singletonvars[i]
    varidx = getCliqMatVarIdx(cliq, varid)
    singslen[i] = lenfcts[varidx]
  end
  p = sortperm(singslen)
  ascsing = singletonvars[p]

  return [ascnons; ascsing]
end

"""
    $SIGNATURES

Determine and return order list of variable ids required for minibatch Gibbs iteration inside `cliq`.

Notes
- Singleton factors (priors and up messages) back of the list
- least number of associated factor variables earlier in list
- Same as mcmcIterationIdsOrdered
"""
function getCliqVarSolveOrderUp(cliq::TreeClique)
  return mcmcIterationIdsOrdered(cliq)
end

"""
    $(SIGNATURES)

Prepare the variable IDs for nested clique Gibbs mini-batch calculations, by assembing these clique data fields:
- `directPriorMsgIDs`
- `directvarIDs`
- `itervarIDs`
- `msgskipIDs`
- `directFrtlMsgIDs`

"""
function setCliqMCIDs!(cliq::TreeClique)
  getData(cliq).directPriorMsgIDs = directPriorMsgIDs(cliq)

  # NOTE -- directvarIDs are combined into itervarIDs
  getData(cliq).directvarIDs = directAssignmentIDs(cliq)
  # TODO find itervarIDs that have upward child singleton messages and update them last in iter list
  getData(cliq).itervarIDs = mcmcIterationIdsOrdered(cliq)  #mcmcIterationIDs(cliq)

  getData(cliq).msgskipIDs = skipThroughMsgsIDs(cliq)
  getData(cliq).directFrtlMsgIDs = directFrtlMsgIDs(cliq)

  # TODO add initialization sequence var id list too

  nothing
end

# post order tree traversal and build potential functions
function buildCliquePotentials(dfg::G, bt::BayesTree, cliq::TreeClique; solvable::Int=1) where G <: AbstractDFG
    for child in out_neighbors(cliq, bt.bt)#tree
        buildCliquePotentials(dfg, bt, child)
    end
    @info "Get potentials $(getLabel(cliq))"
    setCliqPotentials!(dfg, bt, cliq, solvable=solvable)

    compCliqAssocMatrices!(dfg, bt, cliq)
    setCliqMCIDs!(cliq)

    nothing
end

"""
    $(SIGNATURES)

Return a vector of child cliques to `cliq`.
"""
function childCliqs(treel::BayesTree, cliq::TreeClique)
    childcliqs = Vector{TreeClique}(undef, 0)
    for cl in Graphs.out_neighbors(cliq, treel.bt)
        push!(childcliqs, cl)
    end
    return childcliqs
end
function childCliqs(treel::BayesTree, frtsym::Symbol)
  childCliqs(treel,  whichCliq(treel, frtsym))
end

"""
    $(SIGNATURES)

Return a vector of child cliques to `cliq`.
"""
getChildren(treel::BayesTree, frtsym::Symbol) = childCliqs(treel, frtsym)
getChildren(treel::BayesTree, cliq::TreeClique) = childCliqs(treel, cliq)

"""
    $SIGNATURES

Return a vector of all siblings to a clique, which defaults to not `inclusive` the calling `cliq`.
"""
function getCliqSiblings(treel::BayesTree, cliq::TreeClique, inclusive::Bool=false)::Vector{TreeClique}
  prnt = getParent(treel, cliq)
  if length(prnt) > 0
    allch = getChildren(treel, prnt[1])
  end
  if inclusive
    return allch
  end
  sibs = TreeClique[]
  for ch in allch
    if ch.index != cliq.index
      push!(sibs, ch)
    end
  end
  return sibs
end

"""
    $(SIGNATURES)

Return `cliq`'s parent clique.
"""
function parentCliq(treel::BayesTree, cliq::TreeClique)
    Graphs.in_neighbors(cliq, treel.bt)
end
function parentCliq(treel::BayesTree, frtsym::Symbol)
  parentCliq(treel,  whichCliq(treel, frtsym))
end

"""
    $(SIGNATURES)

Return `cliq`'s parent clique.
"""
getParent(treel::BayesTree, afrontal::Union{Symbol, TreeClique}) = parentCliq(treel, afrontal)

"""
    $SIGNATURES

Return one symbol (a frontal variable) from each clique in the `::BayesTree`.

Notes
- Frontal variables only occur once in a clique per tree, therefore is a unique identifier.

Related:

whichCliq, printCliqHistorySummary
"""
function getTreeAllFrontalSyms(fgl::G, tree::BayesTree) where G <: AbstractDFG
  cliqs = tree.cliques
  syms = Vector{Symbol}(undef, length(cliqs))
  for (id,cliq) in cliqs
    syms[id] = getCliqFrontalVarIds(cliq)[1]
  end
  return syms
end

"""
    $SIGNATURES

Get the `::Condition` variable for a clique, likely used for delaying state transitions in
state machine solver.
"""
getSolveCondition(cliq::TreeClique) = getData(cliq).solveCondition

"""
    $SIGNATURES

Return dictionary of all up belief messages currently in a Bayes `tree`.

Related

getUpMsgs
"""
function getTreeCliqUpMsgsAll(tree::BayesTree)::Dict{Int,TempBeliefMsg}
  allUpMsgs = Dict{Int,TempBeliefMsg}()
  for (idx,cliq) in tree.cliques
    msgs = getUpMsgs(cliq)
    allUpMsgs[cliq.index] = TempBeliefMsg()
    for (lbl,msg) in msgs
      # TODO capture the inferred dimension as part of the upward propagation
      allUpMsgs[cliq.index][lbl] = msg
    end
  end
  return allUpMsgs
end

"""
    $SIGNATURES

Convert tree up messages dictionary to a new dictionary relative to variables specific messages and their depth in the tree

Notes
- Return data in `TempUpMsgPlotting` format:
    Dict{Symbol,   -- is for variable label
     Vector{       -- multiple msgs for the same variable
      Symbol,      -- Clique index
      Int,         -- Depth in tree
      BTD          -- Belief estimate
      inferredDim  -- Information count
     }
"""
function stackCliqUpMsgsByVariable(tree::BayesTree,
                                   tmpmsgs::Dict{Int, TempBeliefMsg}  )::TempUpMsgPlotting
  #
  # start of the return data structure
  stack = TempUpMsgPlotting()

  # look at all the clique level data
  for (cidx,tmpmsg) in tmpmsgs
    # look at all variables up msg from each clique
    for (sym,msgdim) in tmpmsg
      # create a new object for a particular variable if hasnt been seen before
      if !haskey(stack,sym)
        stack[sym] = Vector{Tuple{Symbol, Int, BallTreeDensity, Float64}}()
      end
      # assemble metadata
      cliq = tree.cliques[cidx]
      frt = getCliqFrontalVarIds(cliq)[1]
      # add this belief msg and meta data to vector of variable entry
      push!(stack[sym], (frt, getCliqDepth(tree, cliq),msgdim[1], msgdim[2]))
    end
  end

  return stack
end

"""
    $SIGNATURES

Return the variable order stored in a tree object.
"""
getVariableOrder(treel::BayesTree)::Vector{Symbol} = treel.variableOrder

getEliminationOrder(treel::BayesTree) = treel.variableOrder

"""
    $SIGNATURES

EXPERIMENTAL, Save a Bayes (Junction) tree object to file.

Notes
- Converts and saves to JLD2 format a set of `PackedBayesTreeNodeData` objects.
- IIF issue #481

Related

IIF.loadTree, DFG.saveDFG, DFG.loadDFG, JLD2.@save, JLD2.@load
"""
function saveTree(treel::BayesTree,
                  filepath=joinpath("/tmp","caesar","savetree.jld2") )
  #
  savetree = deepcopy(treel)
  for i in 1:length(savetree.cliques)
    if  getData(savetree.cliques[i]) isa BayesTreeNodeData
      setData!(savetree.cliques[i], convert(PackedBayesTreeNodeData, getData(savetree.cliques[i])))
    end
  end

  JLD2.@save filepath savetree
  return filepath
end

function saveTree(treeArr::Vector{BayesTree},
                  filepath=joinpath("/tmp","caesar","savetrees.jld2") )
  #
  savetree = deepcopy(treeArr)
  for savtre in savetree, i in 1:length(savtre.cliques)
    if getData(savtre.cliques[i]) isa BayesTreeNodeData
      setData!(savtre.cliques[i], convert(PackedBayesTreeNodeData, getData(savtre.cliques[i])))
    end
  end

  JLD2.@save filepath savetree
  return filepath
end

"""
    $SIGNATURES

EXPERIMENTAL, Save a Bayes (Junction) tree object to file.

Notes
- Converts and saves to JLD2 format a set of `PackedBayesTreeNodeData` objects.
- IIF issue #481

Related

IIF.saveTree, DFG.saveDFG, DFG.loadDFG, JLD2.@save, JLD2.@load
"""
function loadTree(filepath=joinpath("/tmp","caesar","savetree.jld2"))
  data = @load filepath savetree

  # convert back to a type that which could not be serialized by JLD2
  if savetree isa Vector
      for savtre in savetree, i in 1:length(savtre.cliques)
        if getData(savtre.cliques[i]) isa PackedBayesTreeNodeData
          setData!(savtre.cliques[i], convert(BayesTreeNodeData, getData(savtre.cliques[i])))
        end
      end
  else
    for i in 1:length(savetree.cliques)
      if getData(savetree.cliques[i]) isa PackedBayesTreeNodeData
        setData!(savetree.cliques[i], convert(BayesTreeNodeData, getData(savetree.cliques[i])))
      end
    end
  end

  # return loaded and converted tree
  return savetree
end


text:
##################################################
# Visualization
##################################################

using RecipesBase
using Colors: RGB, HSV
using Statistics

elementcolor(::Hydrogen) = RGB(1., 1., 1.)
elementcolor(::Carbon) = RGB(0., 0., 0.)
elementcolor(::Nitrogen) = RGB(0., 0., 1.)
elementcolor(::Oxygen) = RGB(1., 0., 0.)
volume(::Hydrogen) = 5
volume(::Carbon) = 7
volume(::Nitrogen) = 8
volume(::Oxygen) = 9

@recipe function f(mol::Molecule; unitradius=15, local_axes=false)
    xs = [position(atom)[1] for atom in mol.atoms]
    ys = [position(atom)[2] for atom in mol.atoms]
    zs = [position(atom)[3] for atom in mol.atoms]
    max_distance = max(maximum(xs) - minimum(xs),
                       maximum(ys) - minimum(ys),
                       maximum(zs) - minimum(zs))
    max_distance += sqrt(max_distance)

    markercolor --> elementcolor.(mol.atoms)
    markersize --> unitradius .* ((3/4π) .* volume.(mol.atoms)) .^ 1/3
    @series begin
        segments = [Float64[] for _ in 1:3]
        linewidth --> 2
        linecolor --> :black
        seriestype := :path3d
        primary := false
        markershape := :none

        for bond in mol.bonds
            p_i = position(bond.atom1)
            p_j = position(bond.atom2)
            push!.(segments, p_i, p_j, NaN)
        end
        (segments[1], segments[2], segments[3])
    end

    if local_axes
        for (j, col) in zip(1:3, (:green, :blue, :red))
            @series begin
                segments = [Float64[] for _ in 1:3]
                linewidth --> 5
                linecolor --> col
                seriestype := :path3d
                primary := false
                markershape := :none

                for atom in mol.atoms
                    # for (p0, U) in zip(axes.origins, axes.axes)
                    p_i = position(atom)
                    p_j = p_i + atom.axes[:, j]
                    push!.(segments, p_i, p_j, NaN)
                end
                (segments[1], segments[2], segments[3])
            end
        end
    end

    xlims = mean(xs) - max_distance/2, mean(xs) + max_distance/2
    ylims = mean(ys) - max_distance/2, mean(ys) + max_distance/2
    zlims = mean(zs) - max_distance/2, mean(zs) + max_distance/2
    xticks = floor(xlims[1]):ceil(xlims[2])
    yticks = floor(ylims[1]):ceil(ylims[2])
    zticks = floor(zlims[1]):ceil(zlims[2])

    xlims --> xlims
    xticks --> (xticks, string.(xticks))
    ylims --> ylims
    yticks --> (yticks, string.(yticks))
    zlims --> zlims
    zticks --> (zticks, string.(zticks))

    (xs, ys, zs)
end

# Projection

@recipe function f(mol::Molecule, weights::Vector; unitradius=15, local_axes=false, pscale = 10)
    xs = [position(atom)[1] for atom in mol.atoms]
    ys = [position(atom)[2] for atom in mol.atoms]
    zs = [position(atom)[3] for atom in mol.atoms]
    max_distance = max(maximum(xs) - minimum(xs),
                       maximum(ys) - minimum(ys),
                       maximum(zs) - minimum(zs))
    max_distance += sqrt(max_distance)

    xlims = mean(xs) - max_distance/2, mean(xs) + max_distance/2
    ylims = mean(ys) - max_distance/2, mean(ys) + max_distance/2
    zlims = mean(zs) - max_distance/2, mean(zs) + max_distance/2
    xticks = floor(xlims[1]):ceil(xlims[2])
    yticks = floor(ylims[1]):ceil(ylims[2])
    zticks = floor(zlims[1]):ceil(zlims[2])

    xlims --> xlims
    xticks --> (xticks, string.(xticks))
    ylims --> ylims
    yticks --> (yticks, string.(yticks))
    zlims --> zlims
    zticks --> (zticks, string.(zticks))

    @series begin
        xs = [position(atom)[1] for atom in mol.atoms]
        ys = [position(atom)[2] for atom in mol.atoms]
        zs = [position(atom)[3] for atom in mol.atoms]
        markercolor --> elementcolor.(mol.atoms)
        markersize --> unitradius .* ((3/4π) .* volume.(mol.atoms)) .^ 1/3
        seriesalpha --> 0.1

        (xs, ys, zs)
    end

    @series begin
        segments = [Float64[] for _ in 1:3]
        linewidth --> 2
        linecolor --> :black
        seriestype := :path3d
        primary := false
        markershape := :none

        for bond in mol.bonds
            p_i = position(bond.atom1)
            p_j = position(bond.atom2)
            push!.(segments, p_i, p_j, NaN)
        end
        (segments[1], segments[2], segments[3])
    end

    if local_axes
        for (j, col) in zip(1:3, (:green, :blue, :red))
            @series begin
                segments = [Float64[] for _ in 1:3]
                linewidth --> 5
                linecolor --> col
                seriestype := :path3d
                primary := false
                markershape := :none

                for atom in mol.atoms
                    # for (p0, U) in zip(axes.origins, axes.axes)
                    p_i = position(atom)
                    p_j = p_i + atom.axes[:, j]
                    push!.(segments, p_i, p_j, NaN)
                end
                (segments[1], segments[2], segments[3])
            end
        end
    end

    @series begin
        xs, ys, zs = [Float64[] for _ in 1:3]
        markercolor := markercolor = []#HSV.(180/π .* angle.(weights), 1, 1)
        markersize := markersize = []#unitradius .* abs.(weights)

        i = 1
        for atom in mol.atoms
            x, y, z = position(atom)
            axes = unitradius .* atom.axes ./ pscale
            for orbital in atom.orbitals
                if orbital isa Orbital"*s"
                    push!(xs, x)
                    push!(ys, y)
                    push!(zs, z)
                    push!(markersize, unitradius*abs(weights[i]))
                    push!(markercolor, HSV(180/π*angle(weights[i]), 1, 1))
                end
                for (idx, O) in ((1, Orbital"*p_x"), (2, Orbital"*p_y"), (3, Orbital"*p_z"))
                    if orbital isa O
                        a = axes .* abs(weights[i])
                        push!(xs, x - a[1, idx], x + a[1, idx])
                        push!(ys, y - a[2, idx], y + a[2, idx])
                        push!(zs, z - a[3, idx], z + a[3, idx])
                        for c in 0:1
                            push!(markersize, unitradius*abs(weights[i]))
                            push!(markercolor, HSV(180c + 180/π*angle(weights[i]), 1, 1))
                        end
                    end
                end
                i += 1
            end
        end

        (xs, ys, zs)
    end
end

# With arrows

@recipe function f(g::Function, mol::Molecule, arrows::Vector; unitradius=15)
    xs = [position(atom)[1] for atom in mol.atoms]
    ys = [position(atom)[2] for atom in mol.atoms]
    zs = [position(atom)[3] for atom in mol.atoms]
    max_distance = max(maximum(xs) - minimum(xs),
                       maximum(ys) - minimum(ys),
                       maximum(zs) - minimum(zs))
    max_distance += sqrt(max_distance)

    markercolor --> elementcolor.(mol.atoms)
    markersize --> unitradius .* ((3/4π) .* volume.(mol.atoms)) .^ 1/3
    @series begin
        segments = [Float64[] for _ in 1:3]
        linewidth --> 2
        linecolor --> :black
        seriestype := :path3d
        primary := false
        markershape := :none

        for bond in mol.bonds
            p_i = position(bond.atom1)
            p_j = position(bond.atom2)
            push!.(segments, p_i, p_j, NaN)
        end
        (segments[1], segments[2], segments[3])
    end

    for (j, col) in zip(1:3, (:green, :blue, :red))
        @series begin
            segments = [Float64[] for _ in 1:3]
            linewidth --> 5
            linecolor --> col
            seriestype := :path3d
            primary := false
            markershape := :none

            for (atom, arrow) in zip(filter(g, mol.atoms), arrows)
                # for (p0, U) in zip(axes.origins, axes.axes)
                p_i = position(atom)
                p_j = p_i + arrow[j]
                push!.(segments, p_i, p_j, NaN)
            end
            (segments[1], segments[2], segments[3])
        end
    end

    xlims = mean(xs) - max_distance/2, mean(xs) + max_distance/2
    ylims = mean(ys) - max_distance/2, mean(ys) + max_distance/2
    zlims = mean(zs) - max_distance/2, mean(zs) + max_distance/2
    xticks = floor(xlims[1]):ceil(xlims[2])
    yticks = floor(ylims[1]):ceil(ylims[2])
    zticks = floor(zlims[1]):ceil(zlims[2])

    xlims --> xlims
    xticks --> (xticks, string.(xticks))
    ylims --> ylims
    yticks --> (yticks, string.(yticks))
    zlims --> zlims
    zticks --> (zticks, string.(zticks))

    (xs, ys, zs)
end


text:
module Neo4j2Julia
    include("jolt/Neo4jconnection.jl")
    include("jolt/jolt_parsing.jl")
    include("plotting/plot.jl")

    using MetaGraphs
    using LightGraphs

    export create_graph
    export plot

    """
        create_graph(...)

    Creates a MetaGraph connecting given a Neo4j query and authentication

    # Arguments
    - `usr::String`: the username for the Neo4j database (defualt: neo4j)
    - `passw::String`: the password for the Neo4j database (defualt: neo4j)
    - `URI::String`: the URI for connecting the Neo4j database (defualt: localhost)
    - `port::String`: the port where the database is accesible (given URI) (defualt: 7474)
    - `query::String`: the query to create the graph from (defualt: all database)

    # Returns
    - `MetaGraph`: the graph

    # Throws
    - `IOError[...](ECONNREFUSED)[...]`: If there was no possible conection to the database
    """
    function create_graph(
        usr::String = "neo4j",
        passw::String = "neo4j",
        URI::String = "localhost",
        port::String = "7474",
        query::String = "MATCH (n) OPTIONAL MATCH (n)-[r]-() RETURN n, r;")

    # Check whether the database is accesible or not
    check_accessibility(URI, port)

    # Now access the data the query asks for
    splitted_response = get_data(usr, passw, URI, port, query)

    g = MetaDiGraph()
    set_indexing_prop!(g, :id) # the "database" indexes by the id property

    jolt_parse(g, splitted_response)
    end
end


text:
module MasterEquationSolve

using PyCall
function __init__()
    global qt = pyimport("qutip")
    global sp = pyimport("scipy.sparse")
end

include("qutip_interface.jl")
include("liouvillian.jl")
include("mesolve.jl")
include("main.jl")

end

# this has to be outside the module so __init__() has been called already
if abspath(PROGRAM_FILE) == @__FILE__
    MasterEquationSolve.julia_main()
end


text:
# # Platform 9¾
# In this example we will explore Platform 9 and three quarters of King's cross station, London
#
# At the beginning of an academic year of Hogwarts School of Witchcraft and Wizardry a
# large number of students choose to travel using the Hogwards Express.
# The train departs from platform 9 and three quarters, King's Cross Station, London.
# To reach this platform wizards and witches must run into an arch between the platforms 9 and
# 10 to be teleported the 9¾ platform.
# However due a certain number of mishaps/accidents over the years the Ministry of Magic
# has made new rules.
#
# Only those heading into the wall at a predetermined speed of 4.3 m/s may be teleported.
# Too slow and the wizrd might be folling around and too fast might result in a minor injury
# or worse... Being spotted by a muggle which the ministry will then have obliviate.
# PS : The ministry has now enforced the metric system for all official work. Hence the speed
# is mentioned in meters per second.

using PiCraft

# First let us make the platforms (Muggle and Wizard world)
# Let's first make the muggle platform (9 & 10) and then place a platform
# perpendicular to that beneath

p = getTile()
setBlocks(p .+ (1, 0, 1), p .+ (12, 1, 40), Block(43, 5)) # Make a cobblestone platform
setBlocks(p .+ (5, 2, 1), p .+ (8, 12, 40), Block(45))    # For Brick arches

# Now we want to make arches in the brick column.
# For now semi-circular arches will do
# The equation for a general circle is `x^2 + y^2 = r^2` where `r` is the radius
# So in general if we have a point: (a, b, c) and we wish to make a circle around it we might
# do something like this:- (Say in the x-y plane)
# ```
# for x = -r:sign(r):r, y = -r:sign(r):r
#     if x^2 + y^2 <= r^2
#         setBlock(p .+ (x, y, 0), Block(1,0))
#     end
# end
# ```
# However, since we have to make a semicircle we should only iterate y from 0 to r

r = 9
c1 = p .+ (5, 2, 11) # Center for 1st arch
for y = 0:r, z = -r:r
    if y^2 + z^2 < r^2
        setBlocks(c1 .+ (0, y, z), c1 .+ (3, y, z), Block(0))
    end
end

c2 = p .+ (5, 2, 40 - 10)
for y = 0:r, z = -r:r
    if y^2 + z^2 < r^2
        setBlocks(c2 .+ (0, y, z), c2 .+ (3, y, z), Block(0))
    end
end

# ![Platform.gif](../assets/img/Platform9_75/Platform.gif)

# Now we must now choose a reasonable condition to teleport
# One way can be: If player velocity in `z` direction and position is encosed within
# x: p[1] + 5:8 , y: p[2] + 2, z: 20.3 then teleport to desired location.
#
# Now how do we calculate velocity ?
# Velocity is the rate at which an object's position changes w.r.t time.
# We can compute velocity as follows:
# ```
# p1 = getPos()
# sleep(0.05)
# p2 = getPos()
# vel = (p2 - p1)./0.05 # velocity = distance/time
# ```
# This gives us the velocity component of each direction.
# To get the true instantaneous velocity we have to take the norm of this.
# `v = norm([velocity...]) # We effectively does v = sqrt(vel_x^2 + vel_y^2 + vel_z^2)`
# However we only care about a wizard's z-component of velocity which will be
# vel[3]

function getVel()
    p1 = getPos()
    sleep(0.05)
    vel = getPos() .- p1
    return vel[3]/0.05
end

t = time()
while time() - t < 10
    post("Velocity = "*string.(getVel())*" m/s")
end

# ![getVel.png](../assets/img/Platform9_75/getVel.png)

# Since platform 9¾ remains closed for most of the year we can add a time constraint
# to automatically close it.

t = time()
post("Platform 9¾ open!")
while time() - t < 30 # Platform closes in 30 seconds
    v = getVel()
    pos = getPos()
    if v >= 4.20 && v <= 4.32
        if pos[1] < p[1] + 9 && pos[1] > p[1] + 4 && pos[2] == p[2] + 2 && pos[3] > p[3] + 19.3
            setPos(p)
        end
    end
end
post("Platform 9¾ closed!")

# Now, We wish to make an identical copy of this planform.

m = copyModel(p .+ (1, 0, 1), p .+ (12, 12, 40))
buildModel(rotate(m, dir = :l, dim = :y), p .+ (1, 15, 1))
# This built a copy of the platform which is perpendicular and above it.
# Place the new plaform anywhere you desire.
# Then set `p` to the coordinates of the arch.


text:
function stable_test_method(foo::UInt32)
    foo+foo
end

#Used to ensure check_method gets right method
function stable_test_method(foo::Number)
    bar = "unstable"
    if foo > 10
        bar = foo
    end
    bar
end

function unstable_variables(foo::UInt8)
    bar = 1
    if foo > 8
        bar = foo
    end
    "result:$bar"
end

function unstable_return(foo::UInt32)
    if foo > 1
        1
    else
        1.0
    end
end

function unstable_combo(foo::UInt32)
    bar = 1
    if foo > 1
        bar = foo
    end
    bar
end

@testset "Stability Analysis" begin
    #check_function
    @test Tuple{Any, StabilityReport}[((UInt32,), StabilityReport())] == check_function(stable_test_method, ((UInt32,),))

    @test Tuple{Any, StabilityReport}[((UInt32,), StabilityReport()), ((Float64,), StabilityReport(Tuple{Symbol, Type}[(:bar, Union{String, Float64}), (:return, Union{String, Float64})]))] == check_function(stable_test_method, ((UInt32,), (Float64,)))

    #check_method
    @test StabilityReport() == check_method(stable_test_method, (UInt32,))
    @test StabilityReport() == check_method(unstable_return, (UInt32,), Dict(:return=>Number))
    @test StabilityReport() == check_method(unstable_variables, (UInt8,), Dict(:bar=>Integer))

    @test StabilityReport(Tuple{Symbol, Type}[(:bar, Union{UInt8, Int})]) == check_method(unstable_variables, (UInt8,))
    @test StabilityReport(Tuple{Symbol, Type}[(:return, Union{Int, Float64})]) == check_method(unstable_return, (UInt32,))

    @test StabilityReport(Tuple{Symbol, Type}[(:bar, Union{UInt32, Int})]) == check_method(unstable_combo, (UInt32,), Dict(:return=>Number))
    @test StabilityReport(Tuple{Symbol, Type}[(:return, Union{Int, UInt32})]) == check_method(unstable_combo, (UInt32,), Dict(:bar=>Integer))
    @test StabilityReport(Tuple{Symbol, Type}[(:bar, Union{UInt32, Int}), (:return, Union{Int, UInt32})]) == check_method(unstable_combo, (UInt32,))

    #is_stable(::StabilityReport)
    @test is_stable(StabilityReport())
    @test is_stable(StabilityReport(Vector{Tuple{Symbol, Type}}(undef, 0)))
    @test !is_stable(StabilityReport(Tuple{Symbol, Type}[(:x, Number)]))
    @test !is_stable(StabilityReport(Tuple{Symbol, Type}[(:x, Number),
                                                         (:y, Array),
                                                         (:z, AbstractArray{UInt8, 1})]))
    @test !is_stable(StabilityReport(Tuple{Symbol, Type}[(:return, Number)]))
    @test !is_stable(StabilityReport(Tuple{Symbol, Type}[(:x, Number), (:return, Number)]))

    #is_stable(::Vector{StabilityReport})
    @test is_stable([StabilityReport()])
    @test is_stable([StabilityReport(), StabilityReport(), StabilityReport()])
    @test !is_stable([StabilityReport(Tuple{Symbol, Type}[(:return, Number)])])
    @test !is_stable([StabilityReport(), StabilityReport(Tuple{Symbol, Type}[(:return, Number)])])
    @test !is_stable([StabilityReport(Tuple{Symbol, Type}[(:return, Number)]), StabilityReport()])

    #is_stable(::Vector{Tuple{Any, StabilityReport}})
    @test is_stable([((), StabilityReport())])
    @test is_stable([((), StabilityReport()), ((Int32,), StabilityReport()), ((Float64,),StabilityReport())])
    @test !is_stable([((String,),StabilityReport(Tuple{Symbol, Type}[(:return, Number)]))])
    @test !is_stable([((String,),StabilityReport()), ((), StabilityReport(Tuple{Symbol, Type}[(:return, Number)]))])
    @test !is_stable([((Tuple, String),StabilityReport(Tuple{Symbol, Type}[(:return, Number)])), ((),StabilityReport())])
end


text:
using DrWatson
@quickactivate "MEngProject"
using MEngProject, CUDA, DifferentialEquations, PyPlot, NNlib,  ImageFiltering, Images, MEngProject.LaminartKernels, MEngProject.LaminartInitFunc, MEngProject.Utils, BenchmarkTools, Test

using OrdinaryDiffEq, ParameterizedFunctions, LSODA, Sundials, DiffEqDevTools, Noise

batch = 1

files = readdir(datadir("img"))

tspan = (0.0f0, 200f0)

batch_ = string(batch,"_",rand(1000:9999))
mkdir(plotsdir(string("batch",batch_)))
for file in files[2:end]

	p = LaminartInitFunc.parameterInit_conv_gpu(datadir("img",file), Parameters.parameters_f32);

	u0 = cu(reshape(zeros(Float32, p.dim_i, p.dim_j*(5*2+2)), p.dim_i, p.dim_j, 5*2+2,1))

	arr1 = u0[:, :, 1:2,:]
	arr2 = u0[:, :, 1:1,:];

	f = LaminartFunc.LamFunction(
		similar(arr1), #x
		similar(arr1), #m
		similar(arr1), #s
		similar(arr2), #x_lgn,
		similar(arr1), #C,
		similar(arr1), #H_z,
		similar(arr1), # dy_temp,
		similar(arr1), # dm_temp,
		similar(arr1), # dz_temp,
		similar(arr1), # ds_temp,
		similar(arr2), # dv_temp,
		similar(arr1), # H_z_temp,
		similar(arr2), #  V_temp_1,
		similar(arr2), #  V_temp_2,
		similar(arr1), #  Q_temp,
		similar(arr1), #   P_temp
    )

	prob = ODEProblem(f, u0, tspan, p)
# 	@benchmark sol = solve(prob)
	sol = solve(prob)

	for t ∈ [25,50,100,200]
	#for t ∈ [25,50,100,200,400,800]
# 			for t ∈ [25,50,100]

		axMax = findmax(sol(t)[:,:,:,1])[1]

		for k ∈ 1:2:10
			fig, ax = plt.subplots()
			k2=k+1

			im = ax.imshow(sol(t)[:,:,k,1], cmap=matplotlib.cm.PRGn,
						   vmax=axMax, vmin=-axMax)
			im2 = ax.imshow(sol(t)[:,:,k+1,1], cmap=matplotlib.cm.RdBu_r,
						   vmax=axMax, vmin=-axMax, alpha=0.5)

			cbar = fig.colorbar(im2,  shrink=0.9, ax=ax)
			cbar.ax.set_xlabel("\$k=$k2\$")
					cbar = fig.colorbar(im,  shrink=0.9, ax=ax)
			cbar.ax.set_xlabel("\$k=$k\$")
			layer=Utils.layers[k]
				plt.title("Layer: $layer, \$t=$t\$")
				plt.axis("off")
				fig.tight_layout()
			plt.savefig(plotsdir(string("batch",batch_),string(file,"_",t,"_",Utils.la[k],".png")))
		end

		k=11
		fig, ax = plt.subplots()
		k2=k+1

		im = ax.imshow(sol[:,:,k,1,t], cmap=matplotlib.cm.PRGn,
					   vmax=axMax, vmin=-axMax)
		im2 = ax.imshow(sol[:,:,k+1,1,t], cmap=matplotlib.cm.RdBu_r,
					   vmax=axMax, vmin=-axMax, alpha=0.5)

		cbar = fig.colorbar(im2,  shrink=0.9, ax=ax)
		cbar.ax.set_xlabel("\$v^-\$")
				cbar = fig.colorbar(im,  shrink=0.9, ax=ax)
		cbar.ax.set_xlabel("\$v^+\$")

		layer=Utils.layers[k]
			plt.title("Layer: $layer, \$t=$t\$")
			plt.axis("off")
			fig.tight_layout()

		plt.savefig(plotsdir(string("batch",batch_),string(file,"_",t,"_",Utils.la[k],".png")))
	end

# time plot
	fig, axs = plt.subplots()
	lines = ["b--","g--","r-.","c-.","r:","c:","m-","y-","m:","y:","b-.","g-."]
	for k ∈ 1:12
		layer=Utils.layers_1[k]
		axs.plot(sol[findmax(sol[:,:,k,1,end])[2][1],findmax(sol[:,:,k,1,end])[2][2],k,1,:],lines[k], label="$layer")
	end
	axs.set_xlabel("Time")
	axs.set_ylabel("Activation")
	plt.legend()
	fig.tight_layout()
	plt.savefig(plotsdir(string("batch",batch_),string(file,"_time.png")))
end