text:
using AlgebraicPetri
using OrdinaryDiffEq
using DiffEqFlux, Flux
using Catalyst # https://docs.juliahub.com/Catalyst/QGu8a/6.11.0/tutorials/parameter_estimation/
using Plots
using JSON
using CSV
using DataFrames

SIR = LabelledReactionNet{Float64, Float64}([:S=>1000.0, :I=>0.0, :R=>0.0],
                                            (:inf=>0.0008)=>((:S,:I)=>(:I,:I)),
                                            (:rec=>0.1)=>(:I=>:R));
SVIIR = LabelledReactionNet{Float64, Float64}([:S=>1000.0, :V=>0.0, :I_U=>1.0, :I_V=>0.0, :R=>0.0],
                                            (Symbol("β_{UU}")=>0.0008)=>((:S,:I_U)=>(:I_U,:I_U)),
                                            (Symbol("β_{UV}")=>0.0008)=>((:S,:I_V)=>(:I_U,:I_V)),
                                            (Symbol("β_{VU}")=>0.0008)=>((:V,:I_U)=>(:I_V,:I_U)),
                                            (Symbol("β_{VV}")=>0.0008)=>((:V,:I_V)=>(:I_V,:I_V)),
                                            (:γ_U=>0.1)=>(:I_U=>:R),
                                            (:γ_V=>0.1)=>(:I_V=>:R),
                                            (:ν=>0.0)=>(:S=>:V));
SIR_param_guess = [1e-8, 1e-4, 1e7, 1e5, 1e-10]
SVIIR_param_guess = [1e-8, 1e-8,1e-8, 1e-8,1/14, 1/14, 1e-4, 6e5, 1e1, 1e3, 1e-1, 1e-1]
model = SVIIR
param_guess = SVIIR_param_guess
calc_inf(s,p) = s[1]*s[3]*p[1] + s[1]*s[4]*p[2] + s[2]*s[3]*p[3] + s[2]*s[4]*p[4]

def_calc(s, p) = s[1]*s[2]*p[1]
function make_loss(pn, prob, times, data; calc_inf=def_calc, states=Dict(), rates=Dict())
    function loss(p)
        cur_p = exp.(p)
        u0 = exp.(p[(1:ns(pn)) .+ nt(pn)])
        for (k,v) in rates
            cur_p[k] = v
        end
        for (k,v) in states
            u0[k] = v
        end
        prob′ = remake(prob, p=cur_p, u0=u0, tspan=(0.0,150.0))
        sol = solve(prob′, Tsit5())
        sum(abs2, data .- [calc_inf(sol(t), cur_p) for t in times]), sol
    end
end

sviir_rxn = ReactionSystem(SVIIR)

prob = ODEProblem(sviir_rxn, concentrations(model), (0.0,150.0), vcat(rates(model), concentrations(model)))

state_data = CSV.read("test/georgia.csv", DataFrame)

times = 1:(length(state_data[:, :date])-1)
data = state_data[2:end, :cases] .- state_data[1:(end-1), :cases];

avg_range = 3
week_times = 1:(length(data) - avg_range * 2)
week_avg = map((avg_range+1):(length(data)-avg_range)) do i
    sum(data[(i-avg_range):(i+avg_range)]) ./ (2*avg_range+1)
end
plot(week_times, week_avg)

l_func = make_loss(model, prob, week_times, week_avg; calc_inf=calc_inf, rates=Dict(5=>1/14, 6=>1/14), states=Dict(2=>1e-9, 4=>1e-9));

p = DiffEqFlux.sciml_train(l_func,log.(param_guess),ADAM(0.1),maxiters = 1000)

n_solve = solve(remake(prob, p=exp.(p.u), u0=exp.(p.u[(1:ns(model)) .+ nt(model)])), Tsit5())
p_times = range(1,149, length=1000)
#plot(n_solve; labels=reshape(String.(snames(model)), (1,ns(model))))
plot(week_times, week_avg; linewidth=4, label="Georgia Data", yrange=(0,10000), yaxis="New daily infections (persons)", xaxis="Time (days)")
plot!(p_times, [calc_inf(n_solve(t), exp.(p.u)) for t in p_times]; linewidth=4, label="Estimated Data")

plot(n_solve; labels=reshape(String.(snames(model)), (1,ns(model))), linewidth=4, yaxis="Population (persons)", xaxis="Time (days)")


text:
module PICAnalysisTools

export apply_analytic, analytic_slice,
    track_particles, initial_cube,
    downsample_approx,
    mean_quantity, phase_space_mean

using PICDataStructures
using PICDataStructures: approx_target_size, dropgriddims
using RecursiveArrayTools
using Unitful
using FileTrees
using StaticArrays
using NumericalIntegration
using Transducers, ThreadsX
using Statistics
using ProgressLogging

include("analytic.jl")
include("tracking.jl")
include("statistics.jl")

end

ext: jl


text:
export PDDocument

mutable struct PDDocument
    texts::Vector{PDText}
    str::String
    char2text::Vector{Int}
    annotations::Vector
end

function PDDocument(texts::Vector{PDText})
    str = join(map(t -> t.str, texts))
    char2text = Int[]
    for i = 1:length(texts)
        t = texts[i]
        for _ = 1:sizeof(t.str)
            push!(char2text, i)
        end
    end
    PDDocument(texts, str, char2text, [])
end

function annotate!(doc::PDDocument, str::String, label::String, start::Int=1)
    r = findnext(str, doc.str, start)
    if r != nothing
        i = doc.char2text[first(r)]
        j = doc.char2text[last(r)]
        push!(doc.annotations, (i,j,label))
    end
    r
end

function toconll(doc::PDDocument)
    dict = Dict()
    annotations = map(doc.annotations) do (i,j,label)
        id = get!(dict, label, length(dict)+1)
        (i, j, label, id)
    end
    labels = map(doc.texts) do t
        ["O" for _ = 1:length(dict)]
    end
    for (i,j,label,id) in annotations
        if i == j
            labels[i][id] = "S-$label"
        else
            labels[i][id] = "B-$label"
            for k = i+1:j-1
                labels[k][id] = "I-$label"
            end
            labels[j][id] = "E-$label"
        end
    end
    lines = map(zip(doc.texts,labels)) do (t,l)
        strs = [t.str, string(t.page), string(t.fcoord), string(t.gcoord), l...]
        join(strs, "\t")
    end
    collect(lines)
end

function toxml(doc::PDDocument)

end


text:
# Convert slamindb Functions

function getmaxfactorid(conn, session::AbstractString)
  loadtx = transaction(conn)
  query =  "match (n:$(session):FACTOR)
            where not (n:NEWDATA)
            with id(n) as idn, n.exVertexId as nexvid
            order by nexvid desc limit 1
            return idn, nexvid"
  cph = loadtx(query, submit=true)
  exvid = 0
  if length(cph.results[1]["data"]) > 0
    neoid = cph.results[1]["data"][1]["row"][1]
    exvid = cph.results[1]["data"][1]["row"][2]
  end
  return exvid
end

"""
    getnewvertdict(conn, session)

Return a dictionary with frtend and mongo_keys json string information for :NEWDATA
elements in Neo4j database.
"""
function getnewvertdict(conn, session::AbstractString)

  loadtx = transaction(conn)
  query = "match (n:$(session))-[:DEPENDENCE]-(f:NEWDATA:$(session):FACTOR) where n.ready=1 or f.ready=1 return distinct n, f"
  cph = loadtx(query, submit=true)
  # loadresult = commit(loadtx)
  # @show cph.results[1]

  newvertdict = DataStructures.SortedDict{Int, Dict{Symbol, Dict{AbstractString,Any}}, Base.Order.ForwardOrdering}(Base.Order.ForwardOrdering())
  # mongokeydict = Dict{Int, Dict{AbstractString,Any}}()

  for val in cph.results[1]["data"]
    i = 0
    for elem in val["meta"]
      # @show elem["type"]    # @show rdict["type"]
      i+=1
      newvertdict[elem["id"]] = Dict{Symbol, Dict{AbstractString,Any}}()
      for (k,nv) in val["row"][i]
        if Symbol(k) == :frtend
          newvertdict[elem["id"]][Symbol(k)] = JSON.parse(nv)
        else
          newvertdict[elem["id"]][Symbol(k)] = Dict{AbstractString, Any}("val"=> nv)
        end
      end
      # rdict = JSON.parse(val["row"][i]["frtend"])
      # newvertdict[elem["id"]][:frtend] = rdict
      # if haskey(val["row"][i], "mongo_keys")
      #   # @show val["row"][i]["mongo_keys"]
      #   newvertdict[elem["id"]][:mongokeys] = JSON.parse(val["row"][i]["mongo_keys"])
      # end
      # # if uppercase(rdict["type"])=="POSE" || uppercase(rdict["type"])=="FACTOR"
      #   # npsym = Symbol(string("x",parse(Int, rdict["userid"])+1)) # TODO -- fix :x0 requirement
    end
    # println()
  end

  return newvertdict
end

"""
    parseMergeVertAttr(v, elem)

Parse elem dictionary according to thin ```frtend``` interface from other languages,
and merge contents into attributes of ```v::Graphs.ExVertex```.
"""
function parseMergeVertAttr!(v::Graphs.ExVertex, elem)
  mongos = Dict()
  for (k,va) in elem
    if k == :frtend
      v.attributes[string(k)] = JSON.json(va)
    elseif k == :mongo_keys
      mongos = JSON.parse(va["val"])
      v.attributes[string(k)] = va["val"]
    elseif k == :ready
      v.attributes[string(k)] = typeof(va["val"]) == Int ? va["val"] : parse(Int,va["val"])
    else
      warn("setting $(k) to $(typeof(va["val"]))")
      v.attributes[string(k)] = va["val"]  # this is replacing data incorrectly
    end
  end
  return mongos
end

"""
    mergeCloudVertex!(...)

Hodgepodge function to merge data in CloudVertex
"""
function mergeCloudVertex!{T <: AbstractString}(cg::CloudGraph,
      v::Graphs.ExVertex,
      neoNodeId::Int,
      alreadyexists::Bool,
      neo4jNode,
      existlbs::Vector{AbstractString},
      mongos;
      labels::Vector{T}=String[]  )
  #
  cv = CloudVertex()
  if alreadyexists
    # simply fetch existing cloudgraph if it exists, AGAIN
    cv = CloudGraphs.get_vertex(cg, neoNodeId)
    cgv = cloudVertex2ExVertex(cv)
    # NOTE, overwrite all values, assuming this is ONLY HAPPENING WITH VARIABLENODES
    # not checking to unpack Packed types
    for (ke,val) in cgv.attributes
      v.attributes[ke] = val
    end
  else
    cv = exVertex2CloudVertex( v )
    cv.neo4jNode = neo4jNode
    cv.neo4jNodeId = neoNodeId
    # can't fetch in cloudgraphs with thin interface: (fntend, mongo_keys, ready)
    cv.isValidNeoNodeId = true
    filter!(e->e!="NEWDATA",existlbs)
    cv.labels = union(existlbs, labels)
  end

  # want to check the mongo keys anyway, since they could have changed
  mergeBigDataElements!(cv.bigData.dataElements, mongos)
  return cv
end

function mergeBigDataElements!(bdes::Vector{BigDataElement}, mongos::Dict)
  if mongos != nothing
    for (k,oid) in mongos

      haskey = false
      for bd in bdes
        if bd.description == k
          warn("skipping bigDataElement $(bd.description), assumed to be repeat during merge.")
          haskey = true
        end
      end
      if !haskey
        println("transfer mongo oid $((k,oid))")
        # empty data, since we will call update_NeoBigData and not save_BigData
        bdei = CloudGraphs.BigDataElement(k, Vector{UInt8}(), string(oid))
        push!(bdes, bdei);
      end
    end
  # else
  #   println("no mongo")
  end
  nothing
end

function mergeValuesIntoCloudVert!{T <: AbstractString}(fgl::FactorGraph,
      neoNodeId::Int,
      elem,
      uidl,
      v::Graphs.ExVertex;
      labels::Vector{T}=String[]  )
  #

  # why am I getting a node again (because we don't have the labels here)?
  # TODO -- pass labels via elem dictionary
  neo4jNode = Neo4j.getnode(fgl.cg.neo4j.graph, neoNodeId)
  existlbs = Vector{AbstractString}(neo4jNode.metadata["labels"])

  alreadyexists = sum(existlbs .== "NEWDATA") == 0

  # parse dictionary of values retrieved from Neo4j
  mongos = parseMergeVertAttr!(v, elem)

  if !haskey(fgl.cgIDs, uidl)
    fgl.cgIDs[uidl] = neoNodeId
  else
    fgl.cgIDs[uidl] == neoNodeId ? nothing : error("trying to merge neo ids $(fgl.cgIDs[uidl]) and $(neoNodeId)")
  end

  # merge CloudVertex
  cv = mergeCloudVertex!(fgl.cg, v, neoNodeId, alreadyexists, neo4jNode, existlbs, mongos, labels=labels )

  # update merged vertex to database
  CloudGraphs.update_vertex!(fgl.cg, cv, true)
  # CloudGraphs.update_NeoBigData!(fgl.cg, cv)
  nothing
end

function parsePose2Pose2Constraint(lkl, elem)
  msm = split(elem["meas"], ' ')
  cov = zeros(3,3)
  cov[1,2], cov[1,3], cov[2,3] = parse(Float64, msm[5]), parse(Float64, msm[6]), parse(Float64, msm[8])
  cov += cov'
  cov[1,1], cov[2,2], cov[3,3] = parse(Float64, msm[4]), parse(Float64, msm[7]), parse(Float64, msm[9])
  zij = zeros(3,1)
  zij[:,1] = [parse(msm[1]);parse(msm[2]);parse(msm[3])]
  if length(lkl) >=4
    if lkl[4] == "STDEV"
      cov = cov^2
    end
  end
  return Pose2Pose2(zij, cov, [1.0])
end

function recoverConstraintType(cgl::CloudGraph,
            elem;
            mongokeys::Dict=Dict(),
            N::Int=200  )
  #
  lkl = split(elem["lklh"], ' ')
  if lkl[1]=="PR2"
    msm = split(elem["meas"], ' ')
    cov = zeros(3,3)
    cov[1,2], cov[1,3], cov[2,3] = parse(Float64, msm[5]), parse(Float64, msm[6]), parse(Float64, msm[8])
    cov += cov'
    cov[1,1], cov[2,2], cov[3,3] = parse(Float64, msm[4]), parse(Float64, msm[7]), parse(Float64, msm[9])
    zi = zeros(3,1)
    zi[:,1] = [parse(msm[1]);parse(msm[2]);parse(msm[3])]
    return PriorPose2(MvNormal(zi, cov^2))
  elseif lkl[1]=="PTPR2"
    msm = split(elem["meas"], ' ')
    cov = zeros(2,2)
    cov[1,2] = parse(Float64, msm[4])^2
    cov += cov'
    cov[1,1], cov[2,2] = parse(Float64, msm[3])^2, parse(Float64, msm[5])^2
    zi = zeros(2)
    zi[1:2] = [parse(msm[1]);parse(msm[2])]
    return PriorPoint2(MvNormal(zi, cov))
  elseif lkl[1]=="PP2"
    return parsePose2Pose2Constraint(lkl, elem)
    # msm = split(elem["meas"], ' ')
    # cov = zeros(3,3)
    # cov[1,2], cov[1,3], cov[2,3] = parse(Float64, msm[5]), parse(Float64, msm[6]), parse(Float64, msm[8])
    # cov += cov'
    # cov[1,1], cov[2,2], cov[3,3] = parse(Float64, msm[4]), parse(Float64, msm[7]), parse(Float64, msm[9])
    # zij = zeros(3,1)
    # zij[:,1] = [parse(msm[1]);parse(msm[2]);parse(msm[3])]
    # return Pose2Pose2(zij, cov^2, [1.0])
  elseif lkl[1]=="BR"
    msm = split(elem["meas"], ' ')
    return Pose2DPoint2DBearingRange{Normal, Normal}(
                  Normal(parse(msm[1]), parse(Float64, msm[3]) ),
                  Normal(parse(msm[2]),parse(Float64, msm[5]) )  )
  elseif lkl[1]=="rangeBearingMEAS"
    @show rangekey = mongokeys["range"]
    rngs = bin2arr(CloudGraphs.read_MongoData(cgl, rangekey))
    # rngs = Vector{Float64}(elem["range"] )
    @show bearingkey = mongokeys["bearing"]
    bearing = bin2arr(CloudGraphs.read_MongoData(cgl, bearingkey))
    # bearing = Vector{Float64}(elem["bearing"] )
    @show size(rngs), size(bearing)
    prange = resample(kde!(rngs),N)
    pbear = resample(kde!(bearing),N)
    # warn("temporary return")
    # return prange, pbear
    return Pose2DPoint2DBearingRangeDensity(pbear, prange)
  elseif lkl[1]=="rangeMEAS"
    @show rangekey = mongokeys["range"]
    rngs = bin2arr(CloudGraphs.read_MongoData(cgl, rangekey))
    # rngs = Vector{Float64}(elem["range"] )
    @show size(rngs)
    prange = resample(kde!(rngs),N)
    # warn("temporary return")
    # return prange
    return Pose2DPoint2DRangeDensity(prange)
  else
    return error("Don't know how to convert $(lkl[1]) to a factor")
  end
end

function populatenewvariablenodes!(fgl::FactorGraph, newvertdict::SortedDict; N::Int=100)

  for (neoNodeId,elem) in newvertdict
    # @show neoNodeId
    nlbsym = Symbol()
    uidl = 0
    labels = AbstractString["$(fgl.sessionname)"]
    initvals = Array{Float64,2}()
    if elem[:frtend]["t"] == "P"
      uidl = elem[:frtend]["uid"]+1 # TODO -- remove +1 and allow :x0, :l0
      nlbsym = Symbol(string('x', uidl))
      initvals = 0.1*randn(3,N) # TODO -- fix init to proper values
      # v = addNode!(fgl, nlbsym, , 0.01*eye(3), N=N, ready=0, uid=uidl,api=localapi)
      push!(labels,"POSE")
    elseif elem[:frtend]["t"] == "L"
      warn("using hack counter for LANDMARKS uid +200000")
      uidl = elem[:frtend]["tag_id"]+200000 # TODO complete hack
      nlbsym = Symbol(string('l', uidl))
      initvals = 0.1*randn(2,N) # Make sure autoinit still works properly
      push!(labels,"LANDMARK")
    end
    if !haskey(fgl.IDs, nlbsym) && nlbsym != Symbol()
      # @show nlbsym, size(initvals)
      # TODO remove initstdev, deprecated
      v = addNode!(fgl, nlbsym, initvals, 0.01*eye(size(initvals,2)), N=N, ready=0, uid=uidl,api=localapi)
      mergeValuesIntoCloudVert!(fgl, neoNodeId, elem, uidl, v, labels=labels)
      println()
    end
  end
  println("done populating new variables")
  nothing
end

function populatenewfactornodes!(fgl::FactorGraph, newvertdict::SortedDict, maxfuid::Int)
  foffset = 100000
  fuid = maxfuid >= foffset ? maxfuid : maxfuid+foffset # offset factor values
  warn("using hack counter for FACTOR uid starting at $(fuid)")
  # neoNodeId = 63363
  # elem = newvertdict[neoNodeId]
  for (neoNodeId,elem) in newvertdict
    if elem[:frtend]["t"] == "F"
      # @show neoNodeId
      if !haskey(elem,:ready)
        # warn("missing ready field")
        continue
      end
      if Int(elem[:ready]["val"]) != 1
        warn("ready/val field not equal to 1")
        continue
      end
      # verts relating to this factor
      verts = Vector{Graphs.ExVertex}()
      i=0
      for bf in split(elem[:frtend]["btwn"], ' ')
        i+=1
        # uid = 0
        uid = parse(Int,bf)+1
        if ((elem[:frtend]["lklh"][1:2] == "BR" ||
            elem[:frtend]["lklh"] == "rangeBearingMEAS" ||
            elem[:frtend]["lklh"] == "rangeMEAS" ) && i==2 ) ||
            ( elem[:frtend]["lklh"] == "PTPR2 G 2 STDEV" && i==1 )
          # detect bearing range factors being added between pose and landmark
          warn("using hack counter for LANDMARKS uid +$(2*foffset)")
          uid = parse(Int,bf)+ 2*foffset # complete hack
        end
        push!(verts, fgl.g.vertices[uid])
      end
      # the factor type
      usrfnc = !haskey(elem,:mongo_keys) ? recoverConstraintType(fgl.cg, elem[:frtend]) : recoverConstraintType(fgl.cg, elem[:frtend], mongokeys=JSON.parse(elem[:mongo_keys]["val"]) )
      fuid += 1
      vert = addFactor!(fgl, verts, usrfnc, ready=0, api=localapi, uid=fuid, autoinit=true)
      println("at populatenewfactornodes!, btwn="*elem[:frtend]["btwn"])
      mergeValuesIntoCloudVert!(fgl, neoNodeId, elem, fuid, vert, labels=["FACTOR";"$(fgl.sessionname)"])

      # TODO -- upgrade to add multple variables
      # for vv in verts
      vv = verts[end]
      @show vv.label, Base.mean(vv.attributes["data"].val, 2)
      dlapi.updatevertex!(fgl, vv, updateMAPest=false)
      # end
    end
  end
  nothing
end

"""
    updatenewverts!(fgl::FactorGraph; N::Int)

Convert vertices of session in Neo4j DB with Caesar.jl's required data elements
in preparation for MM-iSAMCloudSolve process.
"""
function updatenewverts!(fgl::FactorGraph; N::Int=100)
  sortedvd = getnewvertdict(fgl.cg.neo4j.connection, fgl.sessionname)
  populatenewvariablenodes!(fgl, sortedvd, N=N)
  maxfuid = getmaxfactorid(fgl.cg.neo4j.connection, fgl.sessionname)
  populatenewfactornodes!(fgl, sortedvd, maxfuid)
  nothing
end

"""
    resetentireremotesession(conn, session)

match (n:\$(session))
remove n.backendset, n.ready, n.data, n.bigData, n.label, n.packedType, n.exVertexId, n.shape, n.width
set n :NEWDATA
return n
"""
function resetentireremotesession(conn, session::AbstractString; segment::AbstractString="")
  loadtx = transaction(conn)
  query = segment == "" ? "match (n:$(session)) " : "match (n:$(session):$(segment)) "
  query = query*"where exists(n.frtend)
           remove n.backendset, n.ready, n.data, n.bigData, n.label, n.packedType, n.exVertexId, n.shape, n.width, n.MAP_est
           set n :NEWDATA"
  cph = loadtx(query, submit=true)
  loadresult = commit(loadtx)
  nothing
end

#


text:
"""
P are the pointers
I contains indexes into `strides` (dynamic) and `X` (static)
X static strides
"""
struct GroupedStridedPointers{P,C,B,R,I,X,O}
    ptrs::P
    strides::X
    offsets::O
end

@inline function GroupedStridedPointers{P,C,B,R,I}(ptrs::P, strides::X, offsets::O) where {P,C,B,R,I,X,O}
    GroupedStridedPointers{P,C,B,R,I,X,O}(ptrs, strides, offsets)
end

"""
G is a tuple(tuple((A_ind,A's dim),(A_ind,A's dim)), ())
it gives the groups.
"""
@inline function grouped_strided_pointer(A::Tuple{Vararg{<:AbstractArray,N}}, ::Val{G}) where {N,G}
    grouped_strided_pointer(
        map(memory_reference, A),
        map(contiguous_axis, A),
        map(contiguous_batch_size, A),
        map(stride_rank, A),
        map(bytestrides, A),
        map(offsets, A),
        map(ArrayInterface.dense_dims, A),
        Val{G}()
    )
end

@generated function grouped_strided_pointer(
    ptrs::P, contig_axis::C, batch_sz::B, r::R, x::X, o::O, d::D, ::Val{()}
) where {P,C,B,R,X,O,D}
    # no groups
    N = length(P.parameters)
    q = Expr(:block, Expr(:meta, :inline))
    i = 0
    Ct = Expr(:tuple); Bt = Expr(:tuple); Rt = Expr(:tuple); Xt = Expr(:tuple); Ot = Expr(:tuple); It = Expr(:tuple)
    for n in 1:N
        push!(Ct.args, C.parameters[n].parameters[1])
        push!(Bt.args, B.parameters[n].parameters[1])
        push!(Rt.args, R.parameters[n].parameters[1])
        Xₙ = X.parameters[n];# Oₙ = O.parameters[n];
        Itt = Expr(:tuple)
        if !iszero(length(Xₙ.parameters))
            xₙ = Symbol(:x_,n)
            oₙ = Symbol(:o_,n)
            push!(q.args, Expr(:(=), xₙ, Expr(:ref, :x, n)))
            push!(q.args, Expr(:(=), oₙ, Expr(:ref, :o, n)))
            for j ∈ 1:length(Xₙ.parameters)
                push!(Xt.args, Expr(:ref, xₙ, j))
                push!(Ot.args, Expr(:ref, oₙ, j))
                push!(Itt.args, (i += 1))
            end
        end
        push!(It.args, Itt)
    end
    push!(q.args, :(GroupedStridedPointers{$P,$Ct,$Bt,$Rt,$It}($Xt, $Ot)))
    q
end

function ordered_rank_and_sort(R::NTuple{N,Int}) where {N}
    sp = ntuple(zero, Val{N}())
    r = ntuple(n -> sum(R[n] .≥ R), Val{N}())
    @inbounds for n in 1:N
        sp = Base.setindex(sp, n, r[n])
    end
    (r, sp)::NTuple{2,NTuple{N,Int}}
end

function matching_values(Xₙ, j, Xₚ, k)::Bool
    xₙⱼ = Xₙ[j]
    xₙⱼ === Int && return false
    xₚₖ = Xₚ[k]
    xₚₖ === Int && return false
    xₙⱼ === xₚₖ
end

function check_match(pm, Rₙ, Sₙ, Xₙ, Dₙ, j, Rₚ, Sₚ, Xₚ, Dₚ, k)::Bool
    # if strides match, we're done
    # @show Xₙ, j, Xₚ, k
    matching_values(Xₙ, j, Xₚ, k) && return true
    Rₙⱼ = Rₙ[j]
    Rₚₖ = Rₚ[k]
    if Rₚₖ > 1 && Rₙⱼ > 1
        # Otherwise, we check preceding axis
        nind = Sₙ[Rₙⱼ-1]
        pind = Sₚ[Rₚₖ-1]
        # matching_values(Xₙ, nind, Xₚ, pind) && return true
        # for them being of equal size,
        # @show pm[nind, pind], Dₙ[nind], Dₚ[pind], nind, pind
        if pm[nind, pind] && Dₙ[nind] && Dₚ[pind]
            return check_match(pm, Rₙ, Sₙ, Xₙ, Dₙ, nind, Rₚ, Sₚ, Xₚ, Dₚ, pind)
        end
    # elseif isone(Rₚₖ) && isone(Rₙⱼ)
    end
    false
end

@generated function grouped_strided_pointer(
    ptrs::P, contig_axis::C, batch_sz::B, r::R, x::X, o::O, d::D, ::Val{G}
) where {P,C,B,R,X,O,D,G}
    N = length(P.parameters)
    # We search for stride matches here
    # first we go over the groups, looking for matches
    m = Matrix{Matrix{Bool}}(undef, N, N)
    for g ∈ G
        # each `g` is a tuple of (array ind, paired dim)
        ng = length(g)
        # @show g
        for i ∈ 1:ng
            ai, di = g[i]
            # @show ai, di
            idim = length(X.parameters[ai].parameters)
            for j ∈ i+1:ng
                aj, dj = g[j]
                # @show aj, dj
                jdim = length(X.parameters[aj].parameters)
                pm = if isassigned(m, LinearIndices(m)[ai,aj])
                    m[ai, aj]
                elseif ai < aj
                    m[ai,aj] = m[aj,ai] = fill(false, jdim, idim)
                    # m[aj,ai] = fill(false, idim, jdim)
                elseif ai > aj
                    m[ai,aj] = m[aj,ai] = fill(false, idim, jdim)
                    # m[ai,aj] = fill(false, jdim, idim)
                end
                if ai < aj
                    pm[dj, di] = true
                elseif ai > aj
                    pm[di, dj] = true
                end
            end
        end
    end
    # @show m G
    # Here we check for static size equivalence info
    # m contains information on matching sizes
    q = Expr(:block, Expr(:meta, :inline))
    i = 0
    Ct = Expr(:tuple); Bt = Expr(:tuple); Rt = Expr(:tuple); Xt = Expr(:tuple); Ot = Expr(:tuple); It = Expr(:tuple)
    for n in 1:N
        push!(Ct.args, C.parameters[n].parameters[1])
        push!(Bt.args, B.parameters[n].parameters[1])
        Rₙ = R.parameters[n].parameters[1]
        push!(Rt.args, Rₙ)
        Xₙ = X.parameters[n].parameters
        Itt = Expr(:tuple)
        ndim = length(Xₙ)
        if ndim > 0
            Rₙ, Sₙ = ordered_rank_and_sort(Rₙ)
            Dₙ = D.parameters[n].parameters[1]
            xₙ = Symbol(:x_,n)
            oₙ = Symbol(:o_,n)
            xₙ_oₙ_not_extracted = true
            for j in 1:ndim
                # Here, we now check if we actually need to add info, or if we can find it
                match = false
                Rₙⱼ = (Rₙ[j])::Int
                # nprev_max = isone(Rₙⱼ) ? 0 : n - 1
                for nprev in 1:n - 1
                    # @show n, nprev
                    isassigned(m, LinearIndices(m)[n,nprev]) || continue
                    pm = m[n, nprev] # just accessing lower triangle, but matrix is symmetric
                    # @show pm
                    # pm is ndim_n x ndim_nprev
                    # we need to search back through stride ranks...
                    for k ∈ axes(pm,1)
                        Rₚ, Sₚ = ordered_rank_and_sort(R.parameters[nprev].parameters[1])
                        Dₚ = D.parameters[nprev].parameters[1]
                        Xₚ = X.parameters[nprev].parameters
                        match |= check_match(pm, Rₙ, Sₙ, Xₙ, Dₙ, j, Rₚ, Sₚ, Xₚ, Dₚ, k)
                        if match
                            # must push nprev's ind
                            push!(Itt.args, It.args[nprev].args[k])
                            break
                        end
                    end
                    match && break
                end
                if !match
                    if xₙ_oₙ_not_extracted
                        push!(q.args, Expr(:(=), xₙ, Expr(:ref, :x, n)))
                        push!(q.args, Expr(:(=), oₙ, Expr(:ref, :o, n)))
                        xₙ_oₙ_not_extracted = false
                    end
                    push!(Xt.args, Expr(:ref, xₙ, j))
                    push!(Ot.args, Expr(:ref, oₙ, j))
                    push!(Itt.args, (i += 1))
                end
            end
        end
        push!(It.args, Itt)
    end
    push!(q.args, :(GroupedStridedPointers{$P,$Ct,$Bt,$Rt,$It}(ptrs, $Xt, $Ot)))
    q
end

@generated function stridedpointers(gsp::GroupedStridedPointers{P,C,B,R,I,X,O}) where {P,C,B,R,X,O,I}
    t = Expr(:tuple)
    for i ∈ eachindex(I)
        Iᵢ = I[i]
        Nᵢ = length(Iᵢ)
        p = Expr(:ref, :ptrs, i)
        # curly = Expr(:curly, :StridedPointer, )
        x = Expr(:tuple); o = Expr(:tuple)
        for j ∈ Iᵢ
            push!(x.args, Expr(:ref, :strds, j))
            push!(o.args, Expr(:ref, :offs, j))
        end
        push!(t.args, Expr(:call, :stridedpointer, p, :(Contiguous{$(C[i])}()), :(ContiguousBatch{$(B[i])}()), :(StrideRank{$(R[i])}()), x, o))
    end
    Expr(:block, Expr(:meta,:inline), :(ptrs = gsp.ptrs), :(strds = gsp.strides), :(offs = gsp.offsets), t)
end


text:
@testset "Ml function tests" begin show_testset_section() #Scope for test data

p = exp.((2*pi*im/7) .* collect(0:6))

#=
                  1.0 + 0.0im
   0.6234898018587336 + 0.7818314824680298im
 -0.22252093395631434 + 0.9749279121818236im
   -0.900968867902419 + 0.43388373911755823im
  -0.9009688679024191 - 0.433883739117558im
  -0.2225209339563146 - 0.9749279121818236im
   0.6234898018587334 - 0.7818314824680299im
=#

@testset "cplxpair()" begin show_testset_description()
	rtol = 1e-12

	pA = exp.((2*pi*im/7) .* collect(0:6))
	pA = RSDeltaSigmaPort.cplxpair(pA)

	p = -0.900968867902419 + 0.43388373911755823im
		@test pA[1] ≈ conj(p) rtol = rtol
		@test pA[2] ≈ p       rtol = rtol
	p = -0.22252093395631434 + 0.9749279121818236im
		@test pA[3] ≈ conj(p) rtol = rtol
		@test pA[4] ≈ p       rtol = rtol
	p = 0.6234898018587336 + 0.7818314824680298im
		@test pA[5] ≈ conj(p) rtol = rtol
		@test pA[6] ≈ p       rtol = rtol
	p = 1.0 + 0.0im
		@test pA[7] ≈ p rtol = rtol
end
@testset "poly()" begin show_testset_description()
	@test RSDeltaSigmaPort.poly([3,4,5]) == [1, -12, 47, -60]
end
@testset "interp()" begin show_testset_description()
	abstol = 1e-15
	v = sin.(0:.1:2π)
	OSR = 32
	@test all(abs.(RSDeltaSigmaPort.interp(v, OSR)[1:OSR:end] .- v[:]) .< abstol)
end

end


text:
@testset "Moments" begin
    begin

        # Fake results, with mean 4.0, standard deviation 2.0 and skew 0.0
        fields = [2.0, 4.0, 6.0]
        results_freq = [FrequencySimulationResult(reshape([complex(f)],1,1),[SVector(0.0)],[0.0]) for f in fields]
        results_time = [TimeSimulationResult(reshape([f],1,1),[SVector(0.0)],[0.0]) for f in fields]

        moments_freq = statistical_moments(results_freq, 3)
        @test moments_freq[1][1,1] ≈ 4.0 &&
              moments_freq[2][1,1] ≈ 2.0 &&
              moments_freq[3][1,1] ≈ 0.0

        moments_time = statistical_moments(results_time, 3)
        @test moments_time[1][1,1] ≈ 4.0 &&
              moments_time[2][1,1] ≈ 2.0 &&
              moments_time[3][1,1] ≈ 0.0
    end
end


text:
using IterTools
using LinearAlgebra
using Logging
using Primes

include("ODE.jl")

#------------------------------------------------------------------------------

function eval_frac(f, point)
    if applicable(numerator, f)
        return divexact(evaluate(numerator(f), point), evaluate(denominator(f), point))
    end
    return evaluate(f, point)
end

#------------------------------------------------------------------------------

function differentiate_solution(
        ode::ODE{P},
        params::Dict{P, T},
        ic::Dict{P, T},
        inputs::Dict{P, Array{T, 1}},
        prec::Int
    ) where {T <: Generic.FieldElem, P <: MPolyElem{T}}
    """
    Input: the same as for power_series_solutions
    Output: a tuple consisting of the power series solution and
            a dictionary of the form (u, v) => power series, where u is a state variable
            v is a state of parameter, and the power series is the partial derivative of
            the function u w.r.t. v evaluated at the solution
    """
    @debug "Computing the power series solution of the system"
    ps_sol = power_series_solution(ode, params, ic, inputs, prec)
    ps_ring = parent(Iterators.peel(values(ps_sol))[1])
    for p in ode.parameters
        ps_sol[p] = ps_ring(params[p])
    end
    sol_vec = [ps_sol[v] for v in gens(ode.poly_ring)]

    @debug "Building the variational system at the solution"
    # Y' = AY + B
    vars = vcat(ode.x_vars, ode.parameters)
    SA = MatrixSpace(ps_ring, length(ode.x_vars), length(ode.x_vars))
    A = SA([
        eval_frac(derivative(ode.equations[vars[i]], vars[j]), sol_vec)
        for i in 1:length(ode.x_vars), j in 1:length(ode.x_vars)
    ])
    SB = MatrixSpace(ps_ring, length(ode.x_vars), length(vars))
    B = zero(SB)
    for i in 1:length(ode.x_vars)
        for j in (length(ode.x_vars) + 1):length(vars)
            B[i, j] = eval_frac(derivative(ode.equations[vars[i]], vars[j]), sol_vec)
        end
    end
    # TODO: make use of one() function (problems modulo prime)
    initial_condition = zero(MatrixSpace(base_ring(ode.poly_ring), length(ode.x_vars), length(vars)))
    for i in 1:length(ode.x_vars)
        initial_condition[i, i] = 1
    end

    @debug "Solving the variational system and forming the output"
    sol_var_system = ps_matrix_linear_de(A, B, initial_condition, prec)
    return (
        ps_sol,
        Dict(
            (vars[i], vars[j]) => sol_var_system[i, j]
            for i in 1:length(ode.x_vars), j in 1:length(vars)
        )
    )
end

#------------------------------------------------------------------------------

function differentiate_output(
        ode::ODE{P},
        outputs::Array{P, 1},
        params::Dict{P, T},
        ic::Dict{P, T},
        inputs::Dict{P, Array{T, 1}},
        prec::Int
    ) where {T <: Generic.FieldElem, P <: MPolyElem{T}}
    """
    Similar to differentiate_solution but computes partial derivatives of a prescribed outputs
    returns a list of dictionaries var => dy / dvar
    """
    @debug "Computing partial derivatives of the solution"
    ps_sol, sol_diff = differentiate_solution(ode, params, ic, inputs, prec)
    ps_ring = parent(Iterators.peel(values(ps_sol))[1])
    for p in ode.parameters
        ps_sol[p] = ps_ring(params[p])
    end
    sol_vec = [ps_sol[v] for v in gens(ode.poly_ring)]

    @debug "Evaluating the partial derivatives of the outputs"
    result = []
    for y in outputs
        push!(result, Dict())
        for x in ode.x_vars
            result[end][x] = sum([eval_frac(derivative(y, xx), sol_vec) * sol_diff[(xx, x)] for xx in ode.x_vars])
        end
        for p in ode.parameters
            result[end][p] = sum([eval_frac(derivative(y, xx), sol_vec) * sol_diff[(xx, p)] for xx in ode.x_vars])
            result[end][p] += eval_frac(derivative(y, p), sol_vec)
        end
    end

    return result
end

#------------------------------------------------------------------------------

function get_degree_and_coeffsize(f::MPolyElem{Nemo.fmpq})
    """
    for f being polynomial over QQ returns a tuple
    (degree, max_coef_size)
    """
    if length(f) == 0
        return (0, 1)
    end
    max_coef = 1
    for c in coeffs(f)
        max_coef = max(max_coef, 2 * height_bits(c))
    end
    return (total_degree(f), max_coef)
end

function get_degree_and_coeffsize(f::Generic.Frac{<: MPolyElem{Nemo.fmpq}})
    """
    for f being rational function over QQ returns a tuple
    (degree, max_coef_size)
    """
    num_deg, num_coef = get_degree_and_coeffsize(numerator(f))
    den_deg, den_coef = get_degree_and_coeffsize(denominator(f))
    return (max(num_deg, den_deg), max(num_coef, den_coef))
end

#------------------------------------------------------------------------------

function compute_defect(ode::ODE{P}, outputs::Array{P, 1}, p::Float64 = 0.99) where P <: MPolyElem{Nemo.fmpq}
    """
    Computed the identifiability defect (Definition 2.6) of an ode system with respect
    to the given outputs with probability at least p
    """

    # Computing the prime using Theorem 1.1 from https://doi.org/10.1006/jsco.2002.0532
    @debug "Computing the prime number"
    d, h = 1, 1
    for f in vcat(collect(values(ode.equations)), outputs)
        df, hf = get_degree_and_coeffsize(f)
        d = max(d, df)
        h = max(h, hf)
    end
    mu = ceil(1 / (1 - sqrt(p)))
    n = length(ode.x_vars)
    m = length(outputs)
    r = length(ode.u_vars)
    ell = length(ode.parameters)
    D = 4 * (n + ell)^2 * (n + m) * d
    Dprime = D * (2 * log(n + ell + r + 1) + log(mu * D)) + 4 * (n + ell)^2 * ((n + m) * h + log(2 * n * D))
    prime = Primes.nextprime(Int(ceil(2 * mu * Dprime)))
    @debug "The prime is $prime"
    prime = 2^31 - 1
    F = GF(prime)

    @debug "Reducing the system modulo prime"
    ode_red = reduce_ode_mod_p(ode, prime)
    outputs_red = [_reduce_poly_mod_p(poly, prime) for poly in outputs]
    prec = length(ode.x_vars) + length(ode.parameters)
    params_vals = Dict(p => F(rand(1:prime)) for p in ode_red.parameters)
    ic = Dict(x => F(rand(1:prime)) for x in ode_red.x_vars)
    inputs = Dict{typeof(outputs_red[1]), Array{Nemo.gfp_elem, 1}}(u => [F(rand(1:prime)) for i in 1:prec] for u in ode_red.u_vars)

    @debug "Computing the output derivatives"
    output_derivatives = differentiate_output(ode_red, outputs_red, params_vals, ic, inputs, prec)

    @debug "Building the matrices"
    Jac_x_param = zero(MatrixSpace(F, length(outputs_red) * prec, prec))
    xs_params = vcat(ode_red.x_vars, ode_red.parameters)
    for i in 1:length(outputs_red)
        for j in 1:prec
            for k in 1:length(xs_params)
                Jac_x_param[(i - 1) * prec + j, k] = coeff(output_derivatives[i][xs_params[k]], j - 1)
            end
        end
    end
    Jac_x = Jac_x_param[:, 1:length(ode_red.x_vars)]

    @debug "Computing the ranks"
    return LinearAlgebra.rank(Jac_x) + ell - LinearAlgebra.rank(Jac_x_param)
end

#------------------------------------------------------------------------------

function bound_number_experiments(ode::ODE{P}, outputs::Array{P, 1}, p::Float64=0.99) where P <: MPolyElem{Nemo.fmpq}
    """
    Input:
        - ode, an ODE object representing an ODE system over Q
        - outputs, rational function being the outputs of the system
    Output: a named tuple with field loc and glob giving upper bounds to the
    NumExpLoc and NumExpGlob, respectively (for definitions, see the paper)
    """
    defect_prev = length(ode.parameters)
    for i in 1:(length(ode.parameters) + 1)
        ode_replicated, outputs_replicated = generate_replica(ode, outputs, i)
        defect_cur = compute_defect(ode_replicated, outputs_replicated, 1. - (1. - p) / length(ode.parameters))
        if defect_prev == defect_cur
            return (loc = i - 1, glob = i)
        end
        defect_prev = defect_cur
    end
    throw(Core.ErrorException("Too many steps for defect to stabilize. Something is wrong here"))
end

#------------------------------------------------------------------------------

max_stars_repo_name: pogudingleb/ExperimentsBound

max_issues_repo_name: pogudingleb/ExperimentsBound

max_forks_repo_name: pogudingleb/ExperimentsBound


text:
# Structured matrices

#####
##### `Diagonal`
#####

rrule(::Type{<:Diagonal}, d::AbstractVector) = Diagonal(d), Rule(diag)

rrule(::typeof(diag), A::AbstractMatrix) = diag(A), Rule(Diagonal)

#####
##### `Symmetric`
#####

rrule(::Type{<:Symmetric}, A::AbstractMatrix) = Symmetric(A), Rule(_symmetric_back)

_symmetric_back(ΔΩ) = UpperTriangular(ΔΩ) + LowerTriangular(ΔΩ)' - Diagonal(ΔΩ)
_symmetric_back(ΔΩ::Union{Diagonal,UpperTriangular}) = ΔΩ

#####
##### `Adjoint`
#####

# TODO: Deal with complex-valued arrays as well
rrule(::Type{<:Adjoint}, A::AbstractMatrix{<:Real}) = Adjoint(A), Rule(adjoint)
rrule(::Type{<:Adjoint}, A::AbstractVector{<:Real}) = Adjoint(A), Rule(vec∘adjoint)

rrule(::typeof(adjoint), A::AbstractMatrix{<:Real}) = adjoint(A), Rule(adjoint)
rrule(::typeof(adjoint), A::AbstractVector{<:Real}) = adjoint(A), Rule(vec∘adjoint)

#####
##### `Transpose`
#####

rrule(::Type{<:Transpose}, A::AbstractMatrix) = Transpose(A), Rule(transpose)
rrule(::Type{<:Transpose}, A::AbstractVector) = Transpose(A), Rule(vec∘transpose)

rrule(::typeof(transpose), A::AbstractMatrix) = transpose(A), Rule(transpose)
rrule(::typeof(transpose), A::AbstractVector) = transpose(A), Rule(vec∘transpose)

#####
##### Triangular matrices
#####

rrule(::Type{<:UpperTriangular}, A::AbstractMatrix) = UpperTriangular(A), Rule(Matrix)

rrule(::Type{<:LowerTriangular}, A::AbstractMatrix) = LowerTriangular(A), Rule(Matrix)


text:
module HeroIcons

export HeroIcon

function __init__()
    if !isdir(joinpath(@__DIR__, "..", "deps", "heroicons-0.4.2"))
        error("HeroIcons should be rebuilt.")
    end
end

struct HeroIcon
    name::String
    style::Symbol
    data::String
    css::String
    function HeroIcon(name::String; style::Symbol = :outline, css = "height:24px;width:24px;")
        style in [:outline, :solid] || error("style must be :outline or :solid")
        data = read(joinpath(@__DIR__, "..", "deps", "heroicons-0.4.2", "optimized", "$style", "$name.svg"), String)
        new(name, style, data, css)
    end
end
function Base.show(io::IO, o::HeroIcon)
    println(io, "HeroIcon: \"$(o.name)\"")
    println(io, "  - style: ", o.style)
    println(io, "  - css:   ", o.css)
end

function Base.show(io::IO, ::MIME"text/html", o::HeroIcon)
    show(io, MIME"text/html"(), HTML("""
$(o.data)
"""))
end

end


text:
module GLInfo

using ModernGL

###############################################################################

using GLLists
using GLDebug
using GLExtender

###############################################################################

export createcontextinfo
export get_glsl_version_string

###############################################################################

global GLSL_VERSION = ""

###############################################################################

"""
TODO
"""
function createcontextinfo()
  global GLSL_VERSION
  glsl = split(unsafe_string(@GLCHECK glGetString(GL_SHADING_LANGUAGE_VERSION)), ['.', ' '])
  if length(glsl) >= 2
    glsl = VersionNumber(parse(Int, glsl[1]), parse(Int, glsl[2]))
    GLSL_VERSION = string(glsl.major) * rpad(string(glsl.minor),2,"0")
  else
    error("Unexpected version number string. Please report this bug! GLSL version string: $(glsl)")
  end

  glv = split(unsafe_string(@GLCHECK glGetString(GL_VERSION)), ['.', ' '])
  if length(glv) >= 2
    glv = VersionNumber(parse(Int, glv[1]), parse(Int, glv[2]))
  else
    error("Unexpected version number string. Please report this bug! OpenGL version string: $(glv)")
  end
  dict = Dict{Symbol,Any}(
      :glsl_version   => glsl,
      :gl_version     => glv,
      :gl_vendor      => unsafe_string(@GLCHECK glGetString(GL_VENDOR)),
      :gl_renderer  => unsafe_string(@GLCHECK glGetString(GL_RENDERER)),
      :gl_extensions => getExtensions(), #split(unsafe_string(glGetString(GL_EXTENSIONS))),
  )
end
export createcontextinfo

function get_glsl_version_string()
  if isempty(GLSL_VERSION)
    error("couldn't get GLSL version, GLUTils not initialized, or context not created?")
  end
  return "#version $(GLSL_VERSION)\n"
end
export get_glsl_version_string

"""
TODO
"""
function getInfoLog(obj::GLuint)
  # Return the info log for obj, whether it be a shader or a program.
  isShader = @GLCHECK glIsShader(obj)
  getiv = isShader == GL_TRUE ? glGetShaderiv : glGetProgramiv
  getInfo = isShader == GL_TRUE ? glGetShaderInfoLog : glGetProgramInfoLog
  # Get the maximum possible length for the descriptive error message
  len = GLint[0]
  @GLCHECK getiv(obj, GL_INFO_LOG_LENGTH, len)
  maxlength = len[]
  # TODO: Create a macro that turns the following into the above:
  # maxlength = @glPointer getiv(obj, GL_INFO_LOG_LENGTH, GLint)
  # Return the text of the message if there is any
  if maxlength > 0
    buffer = zeros(GLchar, maxlength)
    sizei = GLsizei[0]
    @GLCHECK getInfo(obj, maxlength, sizei, buffer)
    len = sizei[]
    unsafe_string(pointer(buffer), len)
  else
    ""
  end
end
export getInfoLog

end # GLInfo

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
# using FluidFluidQueues

##
# import Pkg
# Pkg.activate(".")
using DiscretisedFluidQueues, SparseArrays, LinearAlgebra, StableRNGs

# include("../src/operators.jl")
# include("../src/partition.jl")
##

using Test

T = [-2.5 2 0.5; 1 -2 1; 1 2 -3]
C = [0.0; 2.0; -6.0]

model = BoundedFluidQueue(T,C,6.0)
am = augment_model(model)
nodes = 0.0:1.0:6.0
order = 5
mesh = DGMesh(nodes,order)
dq = DiscretisedFluidQueue(am,mesh)
dq_notam = DiscretisedFluidQueue(model,mesh)
rates_lwr = [1.0; -3.0]
rates_upr = -[1.0; 2.0]
ffq_rates = Matrix(transpose([
    -1.0 0.0 3.0
    -1.0 2.0 0.0
    0.0 -2.0 3.0
    0.0 2.0 -3.0
    1.0 -2.0 0.0
    1.0 0.0 -3.0
]))
# ffq_rates = [ones(1,length(nodes)-1);-2*ones(1,length(nodes)-1);-3*ones(1,length(nodes)-1)]
_rates = Rates(ffq_rates,rates_lwr,rates_upr)
ffq = augment_model(FluidFluidQueue(dq_notam,_rates))
aug_rates = ffq.rates.interior
B = build_full_generator(dq)
ffq_B = FluidFluidGenerator(ffq)

for i in 1:n_phases(model)
    for j in 1:n_phases(model)
        @test B[(:,i,:),(:,j,:)]==ffq_B[(:,i),(:,j)]
        @test B[(:,i,:),(:,j,:)]==ffq_B[(:,i,:),(:,j,:)]
    end
end
R = RatesOperator(ffq)
for i in 1:n_phases(am)
    @show i
    lwr_idx = sum(DiscretisedFluidQueues._has_left_boundary(am.S)[1:i])
    upr_idx = sum(DiscretisedFluidQueues._has_right_boundary(am.S)[1:i])
    n₋ = DiscretisedFluidQueues._has_left_boundary(am.S,i) ? Int(rates_lwr[lwr_idx]==0.0) : 0
    n₊ = DiscretisedFluidQueues._has_right_boundary(am.S,i) ? Int(rates_upr[upr_idx]==0.0) : 0
    @test R[(Zero,i)]==zeros(sum(aug_rates[i,:].==0.0)*n_bases_per_cell(mesh)+n₊+n₋)
    n₋ = DiscretisedFluidQueues._has_left_boundary(am.S,i) ? Int(rates_lwr[lwr_idx]>0.0) : 0
    n₊ = DiscretisedFluidQueues._has_right_boundary(am.S,i) ? Int(rates_upr[upr_idx]>0.0) : 0
    @test R[(Plus,i)]==(1/max(1,i-1))*ones(sum(aug_rates[i,:].>0.0)*n_bases_per_cell(mesh)+n₊+n₋)
    n₋ = DiscretisedFluidQueues._has_left_boundary(am.S,i) ? Int(rates_lwr[lwr_idx]<0.0) : 0
    n₊ = DiscretisedFluidQueues._has_right_boundary(am.S,i) ? Int(rates_upr[upr_idx]<0.0) : 0
    @test R[(Minus,i)]==(1/max(1,i-1))*ones(sum(aug_rates[i,:].<0.0)*n_bases_per_cell(mesh)+n₊+n₋)
end
p_idx = Set(FluidFluidQueues._build_index_from_tuple_sign(ffq,(Plus,:,:)))
m_idx = Set(FluidFluidQueues._build_index_from_tuple_sign(ffq,(Minus,:,:)))
z_idx = Set(FluidFluidQueues._build_index_from_tuple_sign(ffq,(Zero,:,:)))
np_idx = Set(FluidFluidQueues._build_index_from_tuple_sign(ffq,(NotPlus,:,:)))
nm_idx = Set(FluidFluidQueues._build_index_from_tuple_sign(ffq,(NotMinus,:,:)))
nz_idx = Set(FluidFluidQueues._build_index_from_tuple_sign(ffq,(NotZero,:,:)))
@test intersect(p_idx,m_idx)==Set{Int64}()
@test intersect(z_idx,m_idx)==Set{Int64}()
@test intersect(z_idx,p_idx)==Set{Int64}()
@test intersect(nz_idx,z_idx)==Set{Int64}()
@test intersect(np_idx,p_idx)==Set{Int64}()
@test intersect(nm_idx,m_idx)==Set{Int64}()
@test intersect(nm_idx,p_idx)==p_idx
@test intersect(nm_idx,z_idx)==z_idx
@test intersect(np_idx,m_idx)==m_idx
@test intersect(np_idx,z_idx)==z_idx
@test intersect(nz_idx,p_idx)==p_idx
@test intersect(nz_idx,m_idx)==m_idx
@test union(p_idx,m_idx,z_idx)==Set(1:length(R))

@test R[Zero]==R[FluidFluidQueues._build_index_from_tuple_sign(ffq,(Zero,:,:))]
@test R[Plus]==R[FluidFluidQueues._build_index_from_tuple_sign(ffq,(Plus,:,:))]
@test R[Minus]==R[FluidFluidQueues._build_index_from_tuple_sign(ffq,(Minus,:,:))]
@test R[NotZero]==R[FluidFluidQueues._build_index_from_tuple_sign(ffq,(NotZero,:,:))]
@test R[NotPlus]==R[FluidFluidQueues._build_index_from_tuple_sign(ffq,(NotPlus,:,:))]
@test R[NotMinus]==R[FluidFluidQueues._build_index_from_tuple_sign(ffq,(NotMinus,:,:))]

D = InOutGenerator(ffq_B,0.0)
psi = build_psi(D)
xi = build_xi(ffq_B,psi)
marginal, point_mass, K = build_limit_dist_operators(ffq_B,D,R,psi,xi)

macro resetrng1()
    return :(StableRNG(04012022))
end
macro resetrng2()
    return :(StableRNG(204012022))
end
rng1 = @resetrng1()
n_sim = 100
XYvals = 6*rand(rng1,n_sim)
phi0 = rand(rng1,1:4,n_sim)
interior_rates = repeat(C,1,length(nodes)-1)
ffq = augment_model(FluidFluidQueue(dq_notam,Rates(interior_rates,[0.0;0.0],[0.0;0.0])))
ffq_notam = FluidFluidQueue(dq_notam,Rates(interior_rates,[0.0;0.0],[0.0;0.0]))
rng2 = @resetrng2()
SFM, YSims = simulate(ffq,FluidFluidQueues.n_jumps_y(10),
    (φ=phi0, X=XYvals, Y=XYvals),
    rng2,
)
rng2 = @resetrng2()
SFM_notam, YSims_notam = simulate(ffq_notam,FluidFluidQueues.n_jumps_y(10),
    (φ=max.(1,phi0.-1), X=XYvals, Y=XYvals),
    rng2,
)

@test SFM.X ≈ YSims
@test SFM_notam.X ≈ YSims_notam
@test SFM_notam.t ≈ SFM.t
@test SFM_notam.φ ≈ max.(1,SFM.φ.-1)


text: include("./dataCode/estimateModels/ApplicAndTestsOnData/estimateDirBinGasDataAndForecastTest.jl")

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
function getGatWhichFisAchieved(res, F,glim)
	# res[4][i,1] is the F values array
	# res[3][i] is the correspoding g value
	# print("gettin g")
	# print(res)
	g=Inf
	for i in 2:size(res[4], 1)
		thisRunValue = res[4][i,1]
		if thisRunValue<=F && res[3][i]<=glim
			print(i)
			print(thisRunValue)
			# If F is really close (but lower than) to res[4][i-1,1], multiplier should be small positive
			# If F is really close (but higher than) to res[4][i,1], multiplier should be close to 1
			g = res[3][i-1] + (res[3][i]-res[3][i-1])*((res[4][i-1,1]-F)/(res[4][i-1,1]-res[4][i,1]))
			break
		end
	end
	print("found g")
	print(g)
	g
end

function getF(res, glim)
	currentBest = Inf
	prevRunValue = Inf
	prevG = Inf
	for i in 1:size(res[4], 1)
		thisRunValue = res[4][i,1]
		if thisRunValue<currentBest && res[3][i]<=glim
			currentBest=thisRunValue
		end
		if res[3][i]>glim
			midValue  = thisRunValue + ((res[3][i]  - glim)/(res[3][i]-prevG))* (prevRunValue-thisRunValue)
			if midValue<currentBest
				currentBest=midValue
			end
			break
		end
		prevRunValue=thisRunValue
		prevG=res[3][i]
	end
	currentBest
end

function getFoverTen(res, glim)
	glim = glim/10
	currentBest = Inf
	prevRunValue = Inf
	prevG = Inf
	for i in 1:size(res[4], 1)
		thisRunValue = res[4][i,1]
		if thisRunValue<currentBest && res[3][i]<=glim
			currentBest=thisRunValue
		end
		if res[3][i]>glim
			midValue  = thisRunValue + ((res[3][i]  - glim)/(res[3][i]-prevG))* (prevRunValue-thisRunValue)
			if midValue<currentBest
				currentBest=midValue
			end
			break
		end
		prevRunValue=thisRunValue
		prevG=res[3][i]
	end
	currentBest
end

function getFtrain(res, glim)
	currentBest = Inf
	prevRunValue = Inf
	prevG = Inf
	for i in 1:size(res[4], 1)
		thisRunValue = res[4][i,end]
		if thisRunValue<currentBest && res[3][i]<=glim
			currentBest=thisRunValue
		end
		if res[3][i]>glim
			midValue  = thisRunValue + ((res[3][i]  - glim)/(res[3][i]-prevG))* (prevRunValue-thisRunValue)
			if midValue<currentBest
				currentBest=midValue
			end
			break
		end
		prevRunValue=thisRunValue
		prevG=res[3][i]
	end
	currentBest
end

function getFfinal(res, glim)
	currentBest = Inf
	prevG = Inf
	prevRunValue = Inf
	for i in 1:size(res[4], 1)
		# println("i = $i")
		thisRunValue = res[4][i,1]
		# println("thisRunValue = $thisRunValue")
		# println("res[3][i] = $(res[3][i])")
		if res[3][i]<=glim && ~isnan(thisRunValue)
			currentBest=thisRunValue
		end
		if res[3][i]>glim
			if  ~isnan(thisRunValue) && ~isinf(thisRunValue)
				currentBest = thisRunValue + ((res[3][i]  - glim)/(res[3][i]-prevG))* (prevRunValue-thisRunValue)
			end
			break
		end
		prevRunValue=thisRunValue
		prevG=res[3][i]
	end
	min(currentBest,res[4][1,1])
end

function getFtrainfinal(res, glim)
	currentBest = Inf
	prevG = Inf
	prevRunValue = Inf
	for i in 1:size(res[4], 1)
		thisRunValue = res[4][i,end]
		if res[3][i]<=glim && ~isnan(thisRunValue)
			currentBest=thisRunValue
		end
		if res[3][i]>glim
			if  ~isnan(thisRunValue) && ~isinf(thisRunValue)
				currentBest = thisRunValue + ((res[3][i]  - glim)/(res[3][i]-prevG))* (prevRunValue-thisRunValue)
			end
			break
		end
		prevRunValue=thisRunValue
		prevG=res[3][i]
	end
	min(currentBest,res[4][1,end])
end

function getPCC(res, glim)
	currentBest = Inf
	for i in 1:size(res[4], 1)
		thisRunValue = -res[4][i,2]
		if thisRunValue<currentBest && res[3][i]<=glim
			currentBest=thisRunValue
		end
	end
	currentBest
end

function getMCC(res, glim)
	currentBest = Inf
	for i in 1:size(res[4], 1)
		thisRunValue = -res[4][i,3]
		if thisRunValue<currentBest && res[3][i]<=glim
			currentBest=thisRunValue
		end
	end
	currentBest
end

# FIND THE LOWEST VALUE!!! Therefore, valueComputer must return the negative of pcc and mcc
function findBestStepsizeFactor(alg::Function, getThisRunValue::Function, bestPossible::Float64, thisString, returnResultIfExists::Function; outputLevel::Int64=0)

	mid = 50

	outputLevel>0 && println("Starting findBestStepsizeFactor, with $(2*mid-1) possible stepsizes, looking for $thisString")

	function getNextI(values)

		outputLevel>2 && println(values)

		# Weird cases
		for i in values
			if i<=bestPossible
				return false
			end
		end

		# Find minimizer
		(bestVal, bestI) = findmin(values)

		possibleUp=0
		possibleDown=0
		prevUpValue=0
		prevDownValue=0
		checkUp=true
		for i in bestI:99
			outputLevel>2 && println("  Looking at i = $i values[i]=$(values[i])")
			if isnan(values[i])
				outputLevel>1 && println(" Might look up! Prev = $(values[i-1]) ")
				possibleUp=i
				prevUpValue= values[i-1]
				break
			elseif i<=97 && values[i]<=values[i+1]<=values[i+2]
				outputLevel>1 && println(" $(values[i]) $(values[i+1]) $(values[i+2])")
				checkUp = false
				break
			elseif i==98 && values[i]<values[i+1]
				outputLevel>1 && println(" $(values[i]) $(values[i+1])")
				checkUp = false
				break
			elseif i==99
				outputLevel>1 && println(" $(values[i])")
				checkUp = false
				break
			end
		end

		checkDown=true
		for i in bestI:-1:1
			outputLevel>2 && println("  Looking at i = $i values[i]=$(values[i])")
			if isnan(values[i])
				outputLevel>1 && println(" Might look down! Prev = $(values[i+1]) ")
				possibleDown=i
				prevDownValue= values[i+1]
				break
			elseif i>=3 && values[i-2]>=values[i-1]>values[i]
				outputLevel>1 && println(" $(values[i-2]) $(values[i-1]) $(values[i])")
				checkDown = false
				break
			elseif i==2 && values[i-1]>values[i]
				outputLevel>1 && println(" $(values[i-1]) $(values[i])")
				checkDown = false
				break
			elseif i==1
				outputLevel>1 && println(" $(values[i])")
				checkDown = false
				break
			end
		end

		if checkUp==false && checkDown==false
			return false
		elseif checkUp==false
			return possibleDown
		elseif checkDown==false
			return possibleUp
		else
			if prevDownValue<prevUpValue
				return possibleDown
			else
				return possibleUp
			end
		end
	end

	# The middle one corresponds to 1 (or to 2^0) !!
	values=ones(2*mid-1)*NaN

	# Always start with -20,-10,0,10

	if returnResultIfExists(-20)==false
		theValueS = getThisRunValue(alg(-20))
		values[mid-20] = theValueS
		outputLevel>1 && println(" For stepsizePower -20 the value is $theValueS")
	end
	if returnResultIfExists(-10)==false
		theValueA = getThisRunValue(alg(-10))
		values[mid-10] = theValueA
		outputLevel>1 && println(" For stepsizePower -10 the value is $theValueA")
	end
	if returnResultIfExists(0)==false
		theValueB = getThisRunValue(alg(0))
		values[mid] = theValueB
		outputLevel>1 && println(" For stepsizePower 0 the value is $theValueB")
	end
	# if returnResultIfExists(10)==false
	# 	theValueC = getThisRunValue(alg(10))
	# 	values[mid+10] = theValueC
	# 	outputLevel>1 && println(" For stepsizePower 10 the value is $theValueC")
	# end

	bestI=0
	currentBest=Inf
	for i in 1:99
		existingResult = returnResultIfExists(i-mid)
		if existingResult != false
			values[i] = getThisRunValue(existingResult)
		end
		if values[i]<currentBest
			bestI = i
			currentBest = values[i]
		end
	end

	# New routine
	nextI = getNextI(values)
	while nextI != false
		values[nextI] = getThisRunValue(alg(nextI-mid))
		nextI = getNextI(values)
	end

	(bestVal, bestI) = findmin(values)
	outputLevel>0 && println("Best stepsize power found: $(bestI-mid)")
	outputLevel>0 && println("Best value found: $bestVal")
	((bestI-mid), bestVal)
end


text:
module cbox_

using IUP

export
	cbox

include("../src/libiup_h.jl")

#using IProfile
#@iprofile begin

#------------------------example html/examples/C/cbox.c ----------------------------------
function cbox()
	IupOpen()	#Initializes IUP

	IupControlsOpen()
	func_1()
	IupShowXY(IupGetHandle("dlg"),IUP_CENTER,IUP_CENTER)

	IupMainLoop()                               # Initializes IUP main loop
	IupClose()                                  # And close it when ready

	# Don't understand. The profiler reports minuscule times (max = 0.03 s)
	# but the whole thing takes about ~3.4 s (first run)
	#@iprofile report
end

# -----------------------------------------------------------------------------------------
function func_1 ()
	img = IupImage(32,32, img_bits1())
	IupSetHandle ("img1", img)
	IupSetAttribute (img, "0", "0 0 0")
	IupSetAttribute (img, "1", "BGCOLOR")
	IupSetAttribute (img, "2", "255 0 0")

	img = IupImage(32,32, img_bits2())
	IupSetHandle ("img2", img)
	IupSetAttribute (img, "0", "0 0 0")
	IupSetAttribute (img, "1", "0 255 0")
	IupSetAttribute (img, "2", "BGCOLOR")
	IupSetAttribute (img, "3", "255 0 0")

	_frm_1 = IupFrame(
				IupVbox(
					IupSetAttributes(IupButton("Button Text"), "CINDEX=1"),
					IupSetAttributes(IupButton(), "IMAGE=img1,CINDEX=2"),
					IupSetAttributes(IupButton(), "IMAGE=img1,IMPRESS=img2,CINDEX=3")
				))
	IupSetAttribute(_frm_1,"TITLE","IupButton")
	IupSetAttribute(_frm_1,"CX","10")
	IupSetAttribute(_frm_1,"CY","180")

	_frm_2 = IupFrame(
				IupVbox(
					IupSetAttributes(IupLabel("Label Text"), "CINDEX=1"),
					IupSetAttributes(IupLabel(), "SEPARATOR=HORIZONTAL,CINDEX=2"),
					IupSetAttributes(IupLabel(), "IMAGE=img1,CINDEX=3")
				))
	IupSetAttribute(_frm_2,"TITLE","IupLabel");
	IupSetAttribute(_frm_2,"CX","200");
	IupSetAttribute(_frm_2,"CY","250");

	_frm_3 = IupFrame(
				IupVbox(
					IupSetAttributes(IupToggle("Toggle Text"), "VALUE=ON,CINDEX=1"),
					IupSetAttributes(IupToggle(), "IMAGE=img1,IMPRESS=img2,CINDEX=2"),
					IupSetAttributes(IupFrame(IupRadio(IupVbox(
					IupSetAttributes(IupToggle("Toggle Text"), "CINDEX=3"),
					IupSetAttributes(IupToggle("Toggle Text"), "CINDEX=4")
					))), "TITLE=IupRadio")
				))
	IupSetAttribute(_frm_3,"TITLE","IupToggle");
	IupSetAttribute(_frm_3,"CX","400");
	IupSetAttribute(_frm_3,"CY","250");

	_text_1 = IupText()
	IupSetAttribute(_text_1,"VALUE","IupText Text");
	IupSetAttribute(_text_1,"SIZE","80x");
	IupSetAttribute(_text_1,"CINDEX","1");
	IupSetAttribute(_text_1,"CX","10");
	IupSetAttribute(_text_1,"CY","100");

	_ml_1 = IupMultiLine()
	IupSetAttribute(_ml_1,"VALUE","IupMultiline Text\nSecond Line\nThird Line");
	IupSetAttribute(_ml_1,"SIZE","80x60");
	IupSetAttribute(_ml_1,"CINDEX","1");
	IupSetAttribute(_ml_1,"CX","200");
	IupSetAttribute(_ml_1,"CY","100");

	_list_1 = IupList()
	IupSetAttribute(_list_1,"VALUE","1");
	IupSetAttribute(_list_1,"1","Item 1 Text");
	IupSetAttribute(_list_1,"2","Item 2 Text");
	IupSetAttribute(_list_1,"3","Item 3 Text");
	IupSetAttribute(_list_1,"CINDEX","1");
	IupSetAttribute(_list_1,"CX","10");
	IupSetAttribute(_list_1,"CY","10");

	_list_2 = IupList()
	IupSetAttribute(_list_2,"DROPDOWN","YES");
	IupSetAttribute(_list_2,"VALUE","2");
	IupSetAttribute(_list_2,"1","Item 1 Text");
	IupSetAttribute(_list_2,"2","Item 2 Text");
	IupSetAttribute(_list_2,"3","Item 3 Text");
	IupSetAttribute(_list_2,"CINDEX","2");
	IupSetAttribute(_list_2,"CX","200");
	IupSetAttribute(_list_2,"CY","10");

	_list_3 = IupList()
	IupSetAttribute(_list_3,"EDITBOX","YES");
	IupSetAttribute(_list_3,"VALUE","3");
	IupSetAttribute(_list_3,"1","Item 1 Text");
	IupSetAttribute(_list_3,"2","Item 2 Text");
	IupSetAttribute(_list_3,"3","Item 3 Text");
	IupSetAttribute(_list_3,"CINDEX","3");
	IupSetAttribute(_list_3,"CX","400");
	IupSetAttribute(_list_3,"CY","10");

	_cnv_1 = IupCanvas()
	IupSetAttribute(_cnv_1,"RASTERSIZE","100x100");
	IupSetAttribute(_cnv_1,"POSX","0");
	IupSetAttribute(_cnv_1,"POSY","0");
	IupSetAttribute(_cnv_1,"BGCOLOR","128 255 0");
	IupSetAttribute(_cnv_1,"CX","400")
	IupSetAttribute(_cnv_1,"CY","150");

	_ctrl_1 = IupVal()
	IupSetAttribute(_ctrl_1,"CX","600");
	IupSetAttribute(_ctrl_1,"CY","200");

	_cbox = IupCbox(_text_1,_ml_1,_list_1,_list_2,_list_3,_cnv_1,_ctrl_1,
				createtree(),
				create_mat(),
				_frm_1,_frm_2,_frm_3
			)
	IupSetAttribute(_cbox,"SIZE","480x200");

	hbox = IupSetAttributes(IupHbox(_cbox), "MARGIN=10x10");

	dlg = IupDialog(hbox);
	IupSetHandle("dlg",dlg);
	IupSetAttribute(dlg,"TITLE","Cbox Test");
end

# -----------------------------------------------------------------------------------------
function createtree()
	tree = IupTree()
	IupSetAttributes(tree, "FONT=COURIER_NORMAL_10,
						NAME=Figures,
						ADDBRANCH=3D,
						ADDBRANCH=2D,
						ADDLEAF1=trapeze,
						ADDBRANCH1=parallelogram,
						ADDLEAF2=diamond,
						ADDLEAF2=square,
						ADDBRANCH4=triangle,
						ADDLEAF5=scalenus,
						ADDLEAF5=isosceles,
						ADDLEAF5=equilateral,
						RASTERSIZE=180x180,
						VALUE=6,
						CTRL=ON,
						SHIFT=ON,
						CX=600,
						CY=10,
						ADDEXPANDED=NO")
	return tree
end

# -----------------------------------------------------------------------------------------
function create_mat()
	mat = IupMatrix()

	IupSetAttribute(mat, "NUMCOL", "1")
	IupSetAttribute(mat, "NUMLIN", "3")
	IupSetAttribute(mat, "NUMCOL_VISIBLE", "1")
	IupSetAttribute(mat, "NUMLIN_VISIBLE", "3")
	IupSetAttribute(mat, "EXPAND", "NO")
	IupSetAttribute(mat, "SCROLLBAR", "NO")

	IupSetAttribute(mat, "0:0", "Inflation")
	IupSetAttribute(mat, "1:0", "Medicine ")
	IupSetAttribute(mat, "2:0", "Food")
	IupSetAttribute(mat, "3:0", "Energy")
	IupSetAttribute(mat, "0:1", "January 2000")
	IupSetAttribute(mat, "1:1", "5.6")
	IupSetAttribute(mat, "2:1", "2.2")
	IupSetAttribute(mat, "3:1", "7.2")

	IupSetAttribute(mat,"BGCOLOR","255 255 255")
	IupSetAttribute(mat,"BGCOLOR1:0","255 128 0")
	IupSetAttribute(mat,"BGCOLOR2:1","255 128 0")
	IupSetAttribute(mat,"FGCOLOR2:0","255 0 128")
	IupSetAttribute(mat,"FGCOLOR1:1","255 0 128")

	IupSetAttribute(mat,"CX","600")
	IupSetAttribute(mat,"CY","250")

	return mat
end

# -----------------------------------------------------------------------------------------
function img_bits2()
	img = uint8([
	 2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2
	,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,2
	,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,2,2
	,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,2,2,2
	,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,2,2,2,2
	,2,2,2,2,2,2,2,2,2,2,3,3,3,3,1,1,3,3,3,3,3,3,3,3,3,3,2,2,2,2,2,2
	,2,2,2,2,2,2,2,2,2,3,3,3,3,3,1,1,3,3,3,3,3,3,3,3,3,2,2,2,2,2,2,2
	,2,2,2,2,2,2,2,2,3,3,3,3,3,3,1,1,3,3,3,3,3,3,3,3,2,2,2,2,2,2,2,2
	,3,3,3,3,3,3,3,3,3,3,3,3,3,3,1,1,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3
	,3,3,3,3,3,3,3,3,3,3,3,3,3,3,1,1,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3
	,3,3,3,0,3,3,3,3,3,3,3,3,3,3,1,1,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3
	,3,3,3,0,3,3,3,3,3,3,3,3,3,3,1,1,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3
	,3,3,3,0,3,0,3,0,3,3,0,3,3,3,1,1,0,3,3,3,0,0,3,0,3,3,0,0,0,3,3,3
	,3,3,3,0,3,0,0,3,0,0,3,0,3,0,1,1,3,0,3,0,3,3,0,0,3,0,3,3,3,0,3,3
	,3,3,3,0,3,0,3,3,0,3,3,0,3,3,1,1,3,0,3,0,3,3,3,0,3,0,3,3,3,0,3,3
	,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1
	,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1
	,3,3,3,0,3,0,3,3,0,3,3,0,3,0,1,1,3,0,3,0,3,3,0,0,3,0,3,3,3,0,3,3
	,3,3,3,0,3,0,3,3,0,3,3,0,3,3,1,1,0,0,3,3,0,0,3,0,3,3,0,0,0,3,3,3
	,3,3,3,3,3,3,3,3,3,3,3,3,3,3,1,1,3,3,3,3,3,3,3,0,3,3,3,3,3,3,3,3
	,3,3,3,3,3,3,3,3,3,3,3,3,3,3,1,1,3,3,3,0,3,3,3,0,3,3,3,3,3,3,3,3
	,3,3,3,3,3,3,3,3,3,3,3,3,3,3,1,1,3,3,3,3,0,0,0,3,3,3,3,3,3,3,3,3
	,3,3,3,3,3,3,3,3,3,3,3,3,3,3,1,1,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3
	,3,3,3,3,3,3,3,3,3,3,3,3,3,3,1,1,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3
	,2,2,2,2,2,2,2,3,3,3,3,3,3,3,1,1,3,3,3,3,3,2,2,2,2,2,2,2,2,2,2,2
	,2,2,2,2,2,2,3,3,3,3,3,3,3,3,1,1,3,3,3,3,2,2,2,2,2,2,2,2,2,2,2,2
	,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,2,2,2,2,2,2,2,2,2,2,2,2
	,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2
	,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2
	,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2
	,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2
	,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2
	])
end

function img_bits1()
	img = uint8([
	 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1
	,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1
	,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1
	,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1
	,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1
	,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1
	,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1
	,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1
	,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2
	,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2
	,2,2,2,0,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2
	,2,2,2,0,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2
	,2,2,2,0,2,0,2,0,2,2,0,2,2,2,0,0,0,2,2,2,0,0,2,0,2,2,0,0,0,2,2,2
	,2,2,2,0,2,0,0,2,0,0,2,0,2,0,2,2,2,0,2,0,2,2,0,0,2,0,2,2,2,0,2,2
	,2,2,2,0,2,0,2,2,0,2,2,0,2,2,2,2,2,0,2,0,2,2,2,0,2,0,2,2,2,0,2,2
	,2,2,2,0,2,0,2,2,0,2,2,0,2,2,0,0,0,0,2,0,2,2,2,0,2,0,0,0,0,0,2,2
	,2,2,2,0,2,0,2,2,0,2,2,0,2,0,2,2,2,0,2,0,2,2,2,0,2,0,2,2,2,2,2,2
	,2,2,2,0,2,0,2,2,0,2,2,0,2,0,2,2,2,0,2,0,2,2,0,0,2,0,2,2,2,0,2,2
	,2,2,2,0,2,0,2,2,0,2,2,0,2,2,0,0,0,0,2,2,0,0,2,0,2,2,0,0,0,2,2,2
	,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,0,2,2,2,2,2,2,2,2
	,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,0,2,2,2,0,2,2,2,2,2,2,2,2
	,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,0,0,0,2,2,2,2,2,2,2,2,2
	,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2
	,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2
	,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1
	,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1
	,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1
	,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1
	,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1
	,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1
	,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1
	,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1
	])
end

#end     # @iprofile begin

end 	# module


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

export threadid, nthreads, @threads

"""
    Threads.threadid()

Get the ID number of the current thread of execution. The master thread has ID `1`.
"""
threadid() = Int(ccall(:jl_threadid, Int16, ())+1)

# Inclusive upper bound on threadid()
"""
    Threads.nthreads()

Get the number of threads available to the Julia process. This is the inclusive upper bound
on `threadid()`.
"""
nthreads() = Int(unsafe_load(cglobal(:jl_n_threads, Cint)))

function _threadsfor(iter,lbody)
    lidx = iter.args[1]         # index
    range = iter.args[2]
    quote
        local threadsfor_fun
        let range = $(esc(range))
        function threadsfor_fun(onethread=false)
            r = range # Load into local variable
            lenr = length(r)
            # divide loop iterations among threads
            if onethread
                tid = 1
                len, rem = lenr, 0
            else
                tid = threadid()
                len, rem = divrem(lenr, nthreads())
            end
            # not enough iterations for all the threads?
            if len == 0
                if tid > rem
                    return
                end
                len, rem = 1, 0
            end
            # compute this thread's iterations
            f = firstindex(r) + ((tid-1) * len)
            l = f + len - 1
            # distribute remaining iterations evenly
            if rem > 0
                if tid <= rem
                    f = f + (tid-1)
                    l = l + tid
                else
                    f = f + rem
                    l = l + rem
                end
            end
            # run this thread's iterations
            for i = f:l
                local $(esc(lidx)) = @inbounds r[i]
                $(esc(lbody))
            end
        end
        end
        if threadid() != 1
            # only thread 1 can enter/exit _threadedregion
            Base.invokelatest(threadsfor_fun, true)
        else
            ccall(:jl_threading_run, Cvoid, (Any,), threadsfor_fun)
        end
        nothing
    end
end

"""
    Threads.@threads

A macro to parallelize a for-loop to run with multiple threads. This spawns `nthreads()`
number of threads, splits the iteration space amongst them, and iterates in parallel.
A barrier is placed at the end of the loop which waits for all the threads to finish
execution, and the loop returns.
"""
macro threads(args...)
    na = length(args)
    if na != 1
        throw(ArgumentError("wrong number of arguments in @threads"))
    end
    ex = args[1]
    if !isa(ex, Expr)
        throw(ArgumentError("need an expression argument to @threads"))
    end
    if ex.head === :for
        if ex.args[1] isa Expr && ex.args[1].head === :(=)
            return _threadsfor(ex.args[1], ex.args[2])
        else
            throw(ArgumentError("nested outer loops are not currently supported by @threads"))
        end
    else
        throw(ArgumentError("unrecognized argument to @threads"))
    end
end

"""
    Threads.@spawn expr

Create and run a [`Task`](@ref) on any available thread. To wait for the task to
finish, call [`wait`](@ref) on the result of this macro, or call [`fetch`](@ref)
to wait and then obtain its return value.

!!! note
    This feature is currently considered experimental.

!!! compat "Julia 1.3"
    This macro is available as of Julia 1.3.
"""
macro spawn(expr)
    thunk = esc(:(()->($expr)))
    var = esc(Base.sync_varname)
    quote
        local task = Task($thunk)
        task.sticky = false
        if $(Expr(:islocal, var))
            push!($var, task)
        end
        schedule(task)
        task
    end
end


text:
@testset "DensityInterface" begin
    using DensityInterface

    d_uv_continous = Normal(-1.5, 2.3)
    d_uv_discrete = Poisson(4.7)
    d_mv = MvNormal([2.3 0.4; 0.4 1.2])
    d_av = Distributions.MatrixReshaped(MvNormal(rand(10)), 2, 5)

    @testset "Distribution" begin
        for d in (d_uv_continous, d_uv_discrete, d_mv, d_av)
            x = rand(d)
            ref_logd_at_x = logpdf(d, x)
            DensityInterface.test_density_interface(d, x, ref_logd_at_x)

            # Stricter than required by test_density_interface:
            @test logfuncdensity(logdensityof(d)) === d
        end

        for di_func in (logdensityof, densityof)
            @test_throws ArgumentError di_func(d_uv_continous, [rand(d_uv_continous) for i in 1:3])
            @test_throws ArgumentError di_func(d_mv, hcat([rand(d_mv) for i in 1:3]...))
            @test_throws ArgumentError di_func(d_av, [rand(d_av) for i in 1:3])
        end
    end
end


text:
module NeuroEvolution

using Cambrian
using YAML
using JSON
using Random
using Formatting

println("NE.jl v0.1")

include("utils.jl")
include("config.jl")

include("gene.jl")
include("NEAT_network.jl")
include("individual.jl")
include("process.jl")

include("mutation.jl")
include("crossover.jl")

include("species.jl")
include("populate.jl")

include("HyperNEAT.jl")

include("loader.jl")

include("evolution.jl")
end # module


text:
# Bearing and Range constraints for 2D

# better to use bearingrange with [uniform bearing], numerical solving issue on 1D
mutable struct Pose2DPoint2DRange <: IncrementalInference.FunctorPairwise
    Zij::Vector{Float64} # bearing and range hypotheses as columns
    Cov::Float64
    W::Vector{Float64}
    Pose2DPoint2DRange() = new()
    Pose2DPoint2DRange(x...) = new(x[1],x[2],x[3])
end
function (pp2r::Pose2DPoint2DRange)(res::Array{Float64},
      idx::Int,
      meas::Tuple{Array{Float64,2}, Array{Float64,1}}, # from getSample
      xi::Array{Float64,2},
      lm::Array{Float64,2}  )
  #
  # TODO -- still need to add multi-hypotheses support here
  # this is the noisy range
  z = pp2r.Zij[1]+meas[1][1,idx]
  # theta = meas[2]
  # @show size(lm), size(xi), size(meas), size(meas[1]), size(meas[2])
  XX = lm[1,idx] - (z*cos(meas[2][idx]) + xi[1,idx])
  YY = lm[2,idx] - (z*sin(meas[2][idx]) + xi[2,idx])
  res[1] = XX^2 + YY^2
  nothing
end
function getSample(pp2::Pose2DPoint2DRange, N::Int=1)
  return (pp2.Cov*randn(1,N),  2*pi*rand(N))
end

#-------------------------------------------------------------------------------
# bearing and range available

mutable struct Pose2DPoint2DBearingRange{B <: Distributions.Distribution, R <: Distributions.Distribution} <: IncrementalInference.FunctorPairwise
    # Zij::Array{Float64,2} # bearing and range hypotheses as columns
    # Cov::Array{Float64,2}
    # W::Array{Float64,1}
    bearing::B
    range::R
    Pose2DPoint2DBearingRange{B,R}() where {B,R} = new{B,R}()
    Pose2DPoint2DBearingRange(x1::B,x2::R) where {B,R} = new{B,R}(x1,x2)
    Pose2DPoint2DBearingRange{B,R}(x1::B,x2::R) where {B,R} = new{B,R}(x1,x2)
end
function getSample(pp2br::Pose2DPoint2DBearingRange, N::Int=1)
  b = rand(pp2br.bearing, N)
  r = rand(pp2br.range, N)
  return ([b';r'], )
end
# define the conditional probability constraint
function (pp2br::Pose2DPoint2DBearingRange)(res::Array{Float64},
        idx::Int,
        meas::Tuple{Array{Float64,2}},
        xi::Array{Float64,2},
        lm::Array{Float64,2} )
  #
  # @show size(lm), size(xi), size(meas[1])
  # @show idx
  # @show meas[1][:,idx]
  # @show xi[:, idx]
  # @show lm[:,idx]
  res[1] = lm[1,idx] - (meas[1][2,idx]*cos(meas[1][1,idx]+xi[3,idx]) + xi[1,idx])
  res[2] = lm[2,idx] - (meas[1][2,idx]*sin(meas[1][1,idx]+xi[3,idx]) + xi[2,idx])
  nothing
end

# Support for database based solving

passTypeThrough(d::FunctionNodeData{Pose2DPoint2DRange}) = d

mutable struct PackedPose2DPoint2DBearingRange <: IncrementalInference.PackedInferenceType
    bmu::Float64 # 0rotations, 1translation in each column
    bsig::Float64
    rmu::Float64
    rsig::Float64
    PackedPose2DPoint2DBearingRange() = new()
    PackedPose2DPoint2DBearingRange(x1, x2, x3, x4) = new(x1, x2, x3, x4)
end
function convert(::Type{Pose2DPoint2DBearingRange{Normal{T}, Normal{T}}}, d::PackedPose2DPoint2DBearingRange) where T <: Real
  return Pose2DPoint2DBearingRange{Distributions.Normal{T}, Distributions.Normal{T}}(Distributions.Normal{T}(d.bmu,d.bsig), Distributions.Normal{T}(d.rmu, d.rsig))
end
function convert(::Type{PackedPose2DPoint2DBearingRange}, d::Pose2DPoint2DBearingRange{Normal{T}, Normal{T}}) where T
  return PackedPose2DPoint2DBearingRange(d.bearing.μ, d.bearing.σ,
                                         d.range.μ,   d.range.σ )
end

#-------------------------------------------------------------------------------
# bearing only available

# this factor type is still a work in progress
mutable struct Pose2DPoint2DBearing{B <: Distributions.Distribution} <: IncrementalInference.FunctorPairwise
    bearing::B
    Pose2DPoint2DBearing{B}() where {B} = new{B}()
    Pose2DPoint2DBearing(x1::B) where {B} = new{B}(x1)
    Pose2DPoint2DBearing{B}(x1::B) where {B} = new{B}(x1)
end
function getSample(pp2br::Pose2DPoint2DBearing, N::Int=1)
  return (rand(pp2br.bearing, N), )
end
# define the conditional probability constraint
function (pp2br::Pose2DPoint2DBearing)(res::Array{Float64},
        idx::Int,
        meas::Tuple,
        xi::Array{Float64,2},
        lm::Array{Float64,2} )
  #
  res[1] = meas[1][idx] - atan2(lm[2,idx]-xi[2,idx], lm[1,idx]-xi[1,idx])
  nothing
end

# to be deprecated
function pack3(xL1, xL2, p1, p2, p3, xF3)
    warn("RoME.BearingRange2D:pack3 to be deprecated")
    X = zeros(3)
    X[p1] = xL1
    X[p2] = xL2
    X[p3] = xF3
    return X
end

function bearrang!(residual::Array{Float64,1}, Z::Array{Float64,1}, X::Array{Float64,1}, L::Array{Float64,1})
  wTb = SE2(X)
  bTl = wTb\[L[1:2];1.0]
  b = atan2(bTl[2],bTl[1])
  residual[1] = Z[2]-norm(bTl[1:2])
  residual[2] = Z[1]-b
  nothing
end

# should be collapsed to use only numericRootGenericRandomizedFnc
function solveLandm(Zbr::Array{Float64,1}, par::Array{Float64,1}, init::Array{Float64,1})
    return numericRoot(bearrang!, Zbr, par, init)
    # return (nlsolve(   (l, res) -> bearrang!(res, Zbr, par, l), init )).zero
end

# old numeric residual function for pose 2 to pose 2 constraint function.
function solvePose2(Zbr::Array{Float64,1}, par::Array{Float64,1}, init::Array{Float64,1})
    # TODO -- rework to ominus oplus and residual type method
    p = collect(1:3);
    shuffle!(p);
    p1 = p.==1; p2 = p.==2; p3 = p.==3
    #@show init, par
    r = nlsolve(    (res, x) -> bearrang!(res, Zbr,  pack3(x[1], x[2], p1, p2, p3, init[p3]), par),
                    [init[p1];init[p2]] )
    return pack3(r.zero[1], r.zero[2], p1, p2, p3, init[p3]);
end

function solveSetSeps(fnc::Function, Zbr::Array{Float64,1}, CovZ::Array{Float64,2},
                      pars::Array{Float64,2}, inits::Array{Float64,2})
    out = zeros(size(inits))
    for i in 1:size(pars,2)
        ent = 1.0*[CovZ[1,1]*randn(); CovZ[2,2]*randn()]
        out[:,i] = fnc((Zbr+ent), vec(pars[:,i]), vec(inits[:,i]) )
    end
    return out
end

# Xid is the one you want to get back
function evalPotential(brpho::Pose2DPoint2DBearingRange, Xi::Array{Graphs.ExVertex,1}, Xid::Int; N::Int=100)
    # TODO -- add null hypothesis here, might even be done one layer higher in call stack
    val = Array{Float64,2}()
    ini = Array{Graphs.ExVertex,1}()
    par = Array{Graphs.ExVertex,1}()
    oth = Array{Graphs.ExVertex,1}()
    ff::Function = +
    nullhyp = 0.0
    # mmodes = 1 < length(brpho.W)
    # implicit equation portion -- multi-dependent function
    if Xid == Xi[1].index # brpho. ## find the pose

        ff = solvePose2
        # ini = brpho.Xi[1]
        par = Xi[2:end]
        for j in 1:(length(Xi)-1)
            push!(ini, Xi[1])
        end
        #println("Xid == brpho.Xi[1].index inits=", size(inits), ", par=", size(pars))
    elseif Xid == Xi[2].index # find landmark
        if length(Xi) > 2
            nullhyp = 0.5 # should be variable weight
            oth = Xi[3]
        end
        ff = solveLandm
        ini = Xi[2]
        par = Xi[1]
    elseif length(Xi) > 2
        nullhyp = 0.5 # should be variable weight
        if Xid == Xi[3].index # find second mode landmark
            ff = solveLandm
            ini = Xi[3]
            oth = Xi[2]
            par = Xi[1]
        end
    end
    if ff == +
        error("Bad evalPotential Pose2Point2DBearingRange")
    end

    # Gamma = Categorical(brpho.W)

    inits = getVal(ini)
    pars = getVal(par)
    others =  length(Xi) > 2 && Xid != Xi[1].index ? getVal(oth) : Union{}
    # add null hypothesis case
    len = length(Xi) > 2 && Xid != Xi[1].index ? size(others,2) : 0
    # gamma = mmodes ? rand(Gamma) : 1
    numnh = floor(Int, 2*nullhyp*len) # this doubles the value count for null cases
    nhvals = zeros(size(inits,1),numnh)
    for i in 1:numnh
        idx = floor(Int,len*rand()+1)
        # nhvals[:,i] = inits[:,idx] # WRONG! this should be the other landmark, not current landmark
        nhvals[:,i] = others[:,idx]
    end

    val = solveSetSeps(ff, vec(brpho.Zij[:,1]), brpho.Cov, pars, inits)

    return [val';nhvals']'
end

function evalPotentialNew(brpho::Pose2DPoint2DBearingRange, Xi::Array{Graphs.ExVertex,1}, Xid::Int; N::Int=100)
    # TODO -- add null hypothesis here, might even be done one layer higher in call stack
    val = Array{Float64,2}()
    ini = Array{Graphs.ExVertex,1}()
    par = Array{Graphs.ExVertex,1}()
    ff::Function = +
    mmodes = 1 < length(brpho.W)
    x = getVal(brpho.Xi[1])
    l1 = getVal(brpho.Xi[2])
    l2 = mmodes ? getVal(brpho.Xi[3]) : nothing
    nPts = size(x,2)

    pars = Array{Float64,2}(size(x))

    # discrete option, marginalized out before message is sent
    Gamma = mmodes ? rand(Categorical(brpho.W),nPts) : ones(Int,nPts)
    L1s = Gamma .== 1
    nl1s = sum(map(Int, L1s))
    nl2s = nPts-nl1s
    # implicit equation portion -- multi-dependent function

    if Xid == brpho.Xi[1].index # find the pose
        ff = solvePose2
        # par = brpho.Xi[2:end]
        push!(ini, brpho.Xi[1])
        for j in 1:nPts
            pars[:,j] = L1s[j] ? l1[:,j] : l2[:,j]
        end
    elseif Xid == brpho.Xi[2].index # find landmark
        ff = solveLandm
        pars = x
        push!(ini,brpho.Xi[2])
    elseif mmodes
        if Xid == brpho.Xi[3].index # find second mode landmark
            ff = solveLandm
            push!(ini,brpho.Xi[3])
            pars = x
        end
    end
    if ff == +
        error("Bad evalPotential Pose2Point2DBearingRange")
    end

    inits = getVal(ini)

    L1 = sample(getKDE(brpho.Xi[2]),numsampls) #?? # TODO you where here

    len = size(inits,2)
    numnh = floor(Int, 2*nullhyp*len) # this doubles the value count for null cases
    nhvals = zeros(size(inits,1),numnh)
    for i in 1:numnh
        idx = floor(Int,len*rand()+1)
        nhvals[:,i] = inits[:,idx]
    end

    val = solveSetSeps(ff, vec(brpho.Zij[:,1]), brpho.Cov, pars, inits)

    return [val';nhvals']'
end

# Solve for Xid, given values from vertices [Xi] and measurement rho
function evalPotential(rho::Pose2DPoint2DRange, Xi::Array{Graphs.ExVertex,1}, Xid::Int; N::Int=100)
  fromX, ret = nothing, nothing
  if Xi[1].index == Xid
    fromX = getVal( Xi[2] )
    ret = deepcopy(getVal( Xi[1] )) # carry pose yaw row over if required
    ret[3,:] = 2*pi*rand(size(fromX,2))-pi
  elseif Xi[2].index == Xid
    fromX = getVal( Xi[1] )
    ret = deepcopy(getVal( Xi[2] )) # carry pose yaw row over if required
  end
  r,c = size(fromX)
  theta = 2*pi*rand(c)
  noisy = rho.Cov*randn(c) + rho.Zij[1]

  for i in 1:c
    ret[1,i] = noisy[i]*cos(theta[i]) + fromX[1,i]
    ret[2,i] = noisy[i]*sin(theta[i]) + fromX[2,i]
  end

  return ret
end

# no longer needed
# function convert(::Type{PackedFunctionNodeData{PackedPose2DPoint2DBearingRange}}, d::FunctionNodeData{Pose2DPoint2DBearingRange})
#   return PackedFunctionNodeData{PackedPose2DPoint2DBearingRange}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
#           string(d.frommodule), convert(PackedPose2DPoint2DBearingRange, d.fnc))
# end
# function convert(::Type{FunctionNodeData{Pose2DPoint2DBearingRange}}, d::PackedFunctionNodeData{PackedPose2DPoint2DBearingRange})
#   return FunctionNodeData{Pose2DPoint2DBearingRange}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
#           Symbol(d.frommodule), convert(Pose2DPoint2DBearingRange, d.fnc))
# end
# function FNDencode(d::FunctionNodeData{Pose2DPoint2DBearingRange})
#   return convert(PackedFunctionNodeData{PackedPose2DPoint2DBearingRange}, d)
# end
# function FNDdecode(d::PackedFunctionNodeData{PackedPose2DPoint2DBearingRange})
#   return convert(FunctionNodeData{Pose2DPoint2DBearingRange}, d)
# end

# ------------------------------------------------------


text:
module TestStructArrays

using Mutabilities
using Mutabilities: isfrozen, isfrozenindex, isfrozenvalue
using StructArrays
using StaticArrays
using Test

@testset "isfrozen" begin
    @test !isfrozen(StructArray(a = [1]))
    @test !isfrozenindex(StructArray(a = [1]))
    @test !isfrozenvalue(StructArray(a = [1]))

    @test isfrozen(StructArray(a = 1:1))
    @test isfrozenindex(StructArray(a = 1:1))
    @test isfrozenvalue(StructArray(a = 1:1))

    @test !isfrozen(StructArray(a = [1 2]))
    @test isfrozenindex(StructArray(a = [1 2]))
    @test !isfrozenvalue(StructArray(a = [1 2]))
end

@testset "melt" begin
    @test melt(StructArray(a = 1:3)).a isa Vector{Int}
    @test melt(StructArray(a = SA[1, 2])).a isa Vector{Int}
    @test meltvalue(StructArray(a = SA[1, 2])).a isa MArray
    @test meltvalue(StructArray(a = MVector(1, 2))).a isa MArray
    @test meltindex(StructArray(a = SA[1, 2])).a isa Vector{Int}
    @test meltindex(StructArray(a = MVector(1, 2))).a isa Vector{Int}

    x = StructArray{Complex{Int}}(re = SA[1, 2], im = SA[3, 4])
    @test meltindex(x) isa StructArray{Complex{Int}}
    @test meltindex(x).re isa Vector{Int}
    @test meltvalue(x) isa StructArray{Complex{Int}}
    @test meltvalue(x).re isa MArray
    @test typeof(melt(x)) == typeof(meltindex(x))
end

@testset "freeze" begin
    @test freeze(StructArray(a = [1:3;])).a == 1:3
    @test isfrozen(freeze(StructArray(a = [1:3;])).a)
    @test freeze(StructArray(a = SA[1, 2])) === StructArray(a = SA[1, 2])
end

end  # module


text:
module TestMonteCarloPi
using LiterateTest.AssertAsTest: @assert
test() = include(joinpath(@__DIR__, "../examples/monte_carlo_pi.jl"))
end  # module


text:
using Test
using Statistics

include("../src/UncValue.jl")
using ..UncValue

a = Value(3.1415, 0.0012)
b = Value(2.7182818, 3.4e-6)
c = Value(36458.246, 25.64)

@testset "constructor" begin
    @test a.x == 3.1415
    @test a.σ == 0.0012

    @test_throws ErrorException Value(4.1, -2.4)

    d = Value(5.248)
    @test d.x == 5.248
    @test d.σ == 0

    d = Value(8.0, 1)
    @test typeof(d.x) == typeof(d.σ)

    d = Value(8, 0.2)
    @test typeof(d.x) == typeof(d.σ)
end

@testset "precision" begin
    @test precision(a) == -4
    @test precision(b) == -7
    @test precision(c) == 0
end

@testset "addition" begin
    @test typeof(a + 2) <: Value
    @test typeof(2 + a) <: Value
    s = a + b
    @test typeof(s) <: Value
    @test s.x == (a.x + b.x)
    @test s.σ ≈ hypot(a.σ, b.σ)

    n = -a
    @test typeof(n) <: Value
    @test n.x == (-a.x)
    @test n.σ == a.σ
end

@testset "substraction" begin
    @test typeof(a - 2) <: Value
    @test typeof(2 - a) <: Value
    s = a - b
    @test typeof(s) <: Value
    @test s.x ≈ (a.x - b.x)
    @test s.σ ≈ hypot(a.σ, b.σ)
end

@testset "product" begin
    @test typeof(a * 2) <: Value
    @test typeof(2 * a) <: Value
    s = a * b
    @test typeof(s) <: Value
    @test s.x ≈ (a.x * b.x)
    @test s.σ ≈ hypot(a.σ * b.x, a.x * b.σ)
end

@testset "division" begin
    @test typeof(a / 2) <: Value
    @test typeof(2 / a) <: Value
    s = a / b
    @test typeof(s) <: Value
    @test s.x ≈ (a.x / b.x)
    @test s.σ ≈ hypot(a.σ / b.x, a.x * b.σ / b.x^2)
end

@testset "left division" begin
    s = a \ 2
    @test typeof(s) <: Value
    @test s.x ≈ (a.x \ 2)
    @test s.σ ≈ a.x^2 \ 2 * a.σ

    s = 2 \ a
    @test typeof(s) <: Value
    @test s.x ≈ (2 \ a.x)
    @test s.σ ≈ 2 \ a.σ

    s = a \ b
    @test typeof(s) <: Value
    @test s.x ≈ (a.x \ b.x)
    @test s.σ ≈ hypot(a.x \ b.σ, a.x^2 \ b.x * a.σ)
end

@testset "integer division" begin
    j = Value{Int32}(2, 1)
    k = Value{Int32}(7, 2)

    s = j ÷ k.x
    @test typeof(s) <: Value
    @test s.x ≈ j.x ÷ k.x
    @test s.σ ≈ j.σ ÷ k.x

    s = k.x ÷ j
    @test typeof(s) <: Value
    @test s.x ≈ k.x ÷ j.x
    @test s.σ ≈ k.x * j.σ ÷ j.x^2

    s = k ÷ j
    @test typeof(s) <: Value
    @test s.x ≈ k.x ÷ j.x
    @test s.σ ≈ hypot(k.x * j.σ ÷ j.x^2, j.σ ÷ k.x)

    s = fld(j, k.x)
    @test typeof(s) <: Value
    @test s.x ≈ fld(j.x, k.x)
    @test s.σ ≈ fld(j.σ, k.x)

    s = fld(k.x, j)
    @test typeof(s) <: Value
    @test s.x ≈ fld(k.x, j.x)
    @test s.σ ≈ fld(k.x * j.σ, j.x^2)

    s = fld(k, j)
    @test typeof(s) <: Value
    @test s.x ≈ fld(k.x, j.x)
    @test s.σ ≈ floor(hypot(fld(k.x * j.σ, j.x^2), fld(j.σ, k.x)))

    s = cld(j, k.x)
    @test typeof(s) <: Value
    @test s.x ≈ cld(j.x, k.x)
    @test s.σ ≈ cld(j.σ, k.x)

    s = cld(k.x, j)
    @test typeof(s) <: Value
    @test s.x ≈ cld(k.x, j.x)
    @test s.σ ≈ cld(k.x * j.σ, j.x^2)

    s = cld(k, j)
    @test typeof(s) <: Value
    @test s.x ≈ cld(k.x, j.x)
    @test s.σ ≈ ceil(hypot(cld(k.x * j.σ, j.x^2), cld(j.σ, k.x)))
end

@testset "power" begin
    @test typeof(a^3) <: Value
    @test typeof(3^a) <: Value
    s = a^b
    @test typeof(s) <: Value
    @test s.x ≈ a.x^b.x
    @test s.σ ≈ hypot(b.x * a.x^(b.x-1) * a.σ, a.x^b.x * log(abs(a.x)) * b.σ)
end

@testset "equality" begin
    @test a == 3.1415
    @test (a == 3.141) == false
    @test a != 3.141

    @test 3.1415 == a
    @test (3.141 == a) == false
    @test 3.141 != a

    @test (a == b) == false
    @test a != b
end

@testset "inequality" begin
    @test a < 3.1416
    @test (a < 3.141) == false

    @test 3.1416 > a
    @test (3.141 > a) == false

    @test a <= 3.1415
    @test a <= 3.1416
    @test (a <= 3.141) == false

    @test 3.1416 >= a
    @test 3.1415 >= a
    @test (3.141 >= a) == false

    @test a > 3.141
    @test (a > 3.1416) == false

    @test 3.141 < a
    @test (3.1416 < a) == false

    @test a >= 3.1415
    @test a >= 3.1414
    @test (a >= 3.1416) == false

    @test 3.1414 <= a
    @test 3.1415 <= a
    @test (3.146 <= a) == false

    @test a > b
    @test b < a

    @test a >= a
    @test a >= b
    @test (b >= a) == false
    @test a <= a
    @test b <= a
    @test (a <= b) == false
end

@testset "abs" begin
    s = abs(a)
    @test typeof(s) <: Value
    @test s.x == a.x
    @test s.σ == a.σ

    s = abs(-a)
    @test s.x == a.x
    @test s.σ == a.σ
end

@testset "abs2" begin
    s = abs2(a)
    @test typeof(s) <: Value
    @test s.x == a.x^2
    @test s.σ == 2 * a.x * a.σ

    s = abs2(-a)
    @test s.x == a.x^2
    @test s.σ == 2 * a.x * a.σ
end

@testset "sin/cos/tan" begin
    s = sin(b)
    @test typeof(s) <: Value
    @test s.x ≈ sin(b.x)
    @test s.σ ≈ abs(cos(b.x)) * b.σ

    c = cos(b)
    @test typeof(c) <: Value
    @test c.x ≈ cos(b.x)
    @test c.σ ≈ abs(sin(b.x)) * b.σ

    t = tan(b)
    @test typeof(t) <: Value
    @test t.x ≈ tan(b.x)
    @test t.σ ≈ sec(b.x)^2 * b.σ
end

@testset "sind/cosd/tand" begin
    s = sind(b)
    @test typeof(s) <: Value
    @test s.x ≈ sind(b.x)
    @test s.σ ≈ abs(cosd(b.x)) * b.σ

    c = cosd(b)
    @test typeof(c) <: Value
    @test c.x ≈ cosd(b.x)
    @test c.σ ≈ abs(sind(b.x)) * b.σ

    t = tand(b)
    @test typeof(t) <: Value
    @test t.x ≈ tand(b.x)
    @test t.σ ≈ secd(b.x)^2 * b.σ
end

@testset "sinh/cosh/tanh" begin
    s = sinh(b)
    @test typeof(s) <: Value
    @test s.x ≈ sinh(b.x)
    @test s.σ ≈ cosh(b.x) * b.σ

    c = cosh(b)
    @test typeof(c) <: Value
    @test c.x ≈ cosh(b.x)
    @test c.σ ≈ abs(sinh(b.x)) * b.σ

    t = tanh(b)
    @test typeof(t) <: Value
    @test t.x ≈ tanh(b.x)
    @test t.σ ≈ sech(b.x)^2 * b.σ
end

@testset "sinpi/cospi" begin
    s = sinpi(b)
    @test typeof(s) <: Value
    @test s.x ≈ sinpi(b.x)
    @test s.σ ≈ π * abs(cospi(b.x)) * b.σ

    c = cospi(b)
    @test typeof(c) <: Value
    @test c.x ≈ cospi(b.x)
    @test c.σ ≈ π * abs(sinpi(b.x)) * b.σ
end

@testset "asin/acos/atan" begin
    d = Value(0.2435, 0.0658)

    s = asin(d)
    @test typeof(s) <: Value
    @test s.x ≈ asin(d.x)
    @test s.σ ≈ d.σ / sqrt(1 - d.x^2)

    c = acos(d)
    @test typeof(c) <: Value
    @test c.x ≈ acos(d.x)
    @test c.σ ≈ d.σ / sqrt(1 - d.x^2)

    t = atan(d)
    @test typeof(t) <: Value
    @test t.x ≈ atan(d.x)
    @test t.σ ≈ d.σ / (1 + d.x^2)

    t = atan(a, b)
    @test t.x ≈ atan(a.x, b.x)
    @test t.σ ≈ hypot(b.x * a.σ, a.x * b.σ) / hypot(a.x, b.x)^2

    t = atan(a, 0.1)
    @test t.x ≈ atan(a.x, 0.1)
    @test t.σ ≈ 0.1 * a.σ / hypot(a.x, 0.1)^2

    t = atan(0.1, a)
    @test t.x ≈ atan(0.1, a.x)
    @test t.σ ≈ 0.1 * a.σ / hypot(a.x, 0.1)^2
end

@testset "asind/acosd/atand" begin
    d = Value(0.2435, 0.0658)

    s = asind(d)
    @test typeof(s) <: Value
    @test s.x ≈ asind(d.x)
    @test s.σ ≈ d.σ / sqrt(1 - d.x^2)

    c = acosd(d)
    @test typeof(c) <: Value
    @test c.x ≈ acosd(d.x)
    @test c.σ ≈ d.σ / sqrt(1 - d.x^2)

    t = atand(d)
    @test typeof(t) <: Value
    @test t.x ≈ atand(d.x)
    @test t.σ ≈ d.σ / (1 + d.x^2)

    t = atand(a, b)
    @test t.x ≈ atand(a.x, b.x)
    @test t.σ ≈ hypot(b.x * a.σ, a.x * b.σ) / hypot(a.x, b.x)^2

    t = atand(a, 0.1)
    @test t.x ≈ atand(a.x, 0.1)
    @test t.σ ≈ 0.1 * a.σ / hypot(a.x, 0.1)^2

    t = atand(0.1, a)
    @test t.x ≈ atand(0.1, a.x)
    @test t.σ ≈ 0.1 * a.σ / hypot(a.x, 0.1)^2
end

@testset "asinh/acosh/atanh" begin
    d = Value(1.2435, 0.0658)

    s = asinh(d)
    @test typeof(s) <: Value
    @test s.x ≈ asinh(d.x)
    @test s.σ ≈ d.σ / sqrt(1 + d.x^2)

    c = acosh(d)
    @test typeof(c) <: Value
    @test c.x ≈ acosh(d.x)
    @test c.σ ≈ d.σ / sqrt(d.x^2 - 1)

    d -= 1
    t = atanh(d)
    @test typeof(t) <: Value
    @test t.x ≈ atanh(d.x)
    @test t.σ ≈ d.σ / (1 - d.x^2)
end

@testset "csc/sec/cot" begin
    s = csc(b)
    @test typeof(s) <: Value
    @test s.x ≈ csc(b.x)
    @test s.σ ≈ abs(cot(b.x) * csc(b.x)) * b.σ

    c = sec(b)
    @test typeof(c) <: Value
    @test c.x ≈ sec(b.x)
    @test c.σ ≈ abs(tan(b.x) * sec(b.x)) * b.σ

    t = cot(b)
    @test typeof(t) <: Value
    @test t.x ≈ cot(b.x)
    @test t.σ ≈ csc(b.x)^2 * b.σ
end

@testset "cscd/secd/cotd" begin
    s = cscd(b)
    @test typeof(s) <: Value
    @test s.x ≈ cscd(b.x)
    @test s.σ ≈ abs(cotd(b.x) * cscd(b.x)) * b.σ

    c = secd(b)
    @test typeof(c) <: Value
    @test c.x ≈ secd(b.x)
    @test c.σ ≈ abs(tand(b.x) * secd(b.x)) * b.σ

    t = cotd(b)
    @test typeof(t) <: Value
    @test t.x ≈ cotd(b.x)
    @test t.σ ≈ cscd(b.x)^2 * b.σ
end

@testset "csch/sech/coth" begin
    s = csch(b)
    @test typeof(s) <: Value
    @test s.x ≈ csch(b.x)
    @test s.σ ≈ abs(coth(b.x) * csch(b.x)) * b.σ

    c = sech(b)
    @test typeof(c) <: Value
    @test c.x ≈ sech(b.x)
    @test c.σ ≈ abs(tanh(b.x) * sech(b.x)) * b.σ

    t = coth(b)
    @test typeof(t) <: Value
    @test t.x ≈ coth(b.x)
    @test t.σ ≈ csch(b.x)^2 * b.σ
end

@testset "acsc/asec/acot" begin
    s = acsc(b)
    @test typeof(s) <: Value
    @test s.x ≈ acsc(b.x)
    @test s.σ ≈ b.σ / sqrt(b.x^2 - 1)

    c = asec(b)
    @test typeof(c) <: Value
    @test c.x ≈ asec(b.x)
    @test c.σ ≈ b.σ / sqrt(b.x^2 - 1)

    t = acot(b)
    @test typeof(t) <: Value
    @test t.x ≈ acot(b.x)
    @test t.σ ≈ b.σ / (1 + b.x^2)
end

@testset "acscd/asecd/acotd" begin
    s = acscd(b)
    @test typeof(s) <: Value
    @test s.x ≈ acscd(b.x)
    @test s.σ ≈ b.σ / sqrt(b.x^2 - 1)

    c = asecd(b)
    @test typeof(c) <: Value
    @test c.x ≈ asecd(b.x)
    @test c.σ ≈ b.σ / sqrt(b.x^2 - 1)

    t = acotd(b)
    @test typeof(t) <: Value
    @test t.x ≈ acotd(b.x)
    @test t.σ ≈ b.σ / (1 + b.x^2)
end

@testset "acsch/asech/acoth" begin
    d = Value(0.12435, 0.0658)
    s = acsch(b)
    @test typeof(s) <: Value
    @test s.x ≈ acsch(b.x)
    @test s.σ ≈ b.σ / sqrt(b.x^2 + 1)

    c = asech(d)
    @test typeof(c) <: Value
    @test c.x ≈ asech(d.x)
    @test c.σ ≈ d.σ / sqrt(1 - d.x^2)

    t = acoth(b)
    @test typeof(t) <: Value
    @test t.x ≈ acoth(b.x)
    @test t.σ ≈ b.σ / abs(1 - b.x^2)
end

@testset "deg2rad" begin
    d = deg2rad(a)
    @test typeof(d) <: Value
    @test d.x == deg2rad(a.x)
    @test d.σ == deg2rad(a.σ)

    d = rad2deg(a)
    @test typeof(d) <: Value
    @test d.x == rad2deg(a.x)
    @test d.σ == rad2deg(a.σ)
end

@testset "exp" begin
    r = exp(a)
    @test typeof(r) <: Value
    @test r.x == exp(a.x)
    @test r.σ ≈ r.x * a.σ

    r = exp2(a)
    @test typeof(r) <: Value
    @test r.x == exp2(a.x)
    @test r.σ ≈ r.x * a.σ * log(2)

    r = exp10(a)
    @test typeof(r) <: Value
    @test r.x == exp10(a.x)
    @test r.σ ≈ r.x * a.σ * log(10)
end

@testset "log" begin
    r = log(a)
    @test typeof(r) <: Value
    @test r.x == log(a.x)
    @test r.σ ≈ a.σ / a.x

    r = log2(a)
    @test typeof(r) <: Value
    @test r.x == log2(a.x)
    @test r.σ ≈ a.σ / (a.x * log(2))

    r = log10(a)
    @test typeof(r) <: Value
    @test r.x == log10(a.x)
    @test r.σ ≈ a.σ / (a.x * log(10))

    r = log(3, a)
    @test typeof(r) <: Value
    @test r.x == log(3, a.x)
    @test r.σ ≈ a.σ / (a.x * log(3))

    r = log1p(a)
    @test typeof(r) <: Value
    @test r.x == log1p(a.x)
    @test r.σ ≈ a.σ / (a.x + 1)
end

@testset "roots" begin
    r = sqrt(a)
    @test typeof(r) <: Value
    @test r.x == sqrt(a.x)
    @test r.σ == a.σ / (2 * sqrt(a.x))

    r = cbrt(a)
    @test typeof(r) <: Value
    @test r.x == cbrt(a.x)
    @test r.σ == a.σ / (3 * cbrt(a.x)^2)

    h = hypot(a, 3)
    @test typeof(h) <: Value
    @test h.x == hypot(a.x, 3)
    @test h.σ == a.σ * a.x / h

    h = hypot(3, a)
    @test typeof(h) <: Value
    @test h.x == hypot(3, a.x)
    @test h.σ == a.σ * a.x / h

    h = hypot(a, b)
    @test typeof(h) <: Value
    @test h.x == hypot(a.x, b.x)
    @test h.σ == hypot(a.σ * a.x, b.σ * b.x) / h
end

@testset "sign" begin
    @test sign(a) == sign(a.x)
    @test sign(Value(-2.4, 2.4)) == sign(-2.4)

    @test signbit(a) == signbit(a.x)
    @test signbit(Value(-2.4, 2.4)) == signbit(-2.4)
end

@testset "inv" begin
    r = inv(a)
    @test typeof(r) <: Value
    @test r.x == inv(a.x)
    @test r.σ == a.σ / a.x^2
end

@testset "approx" begin
    @test a ≈ 3.141
    @test (a ≈ 3.145) == false
    @test isapprox(a, 3.145; significance=3)

    @test 3.141 ≈ a
    @test (3.145 ≈ a) == false
    @test isapprox(3.145, a; significance=3)

    @test (a ≈ b) == false
end

@testset "cmp" begin
    @test cmp(a, 3) == 1
    @test cmp(3, a) == -1
    @test cmp(a, a) == 0
    @test cmp(a, b) == 1
    @test cmp(b, a) == -1
end

@testset "isless" begin
    @test isless(a, b) == false
    @test isless(b, a)
    @test isless(a, 3) == false
    @test isless(3, a)
end

@testset "clamp" begin
    @test clamp(3, a) == a.x - a.σ
    @test clamp(4, a) == a.x + a.σ
    @test clamp(3.1414, a) == 3.1414
end

@testset "min/max" begin
    @test min(a, b) == min(b, a) == b
    @test max(a, b) == max(b, a) == a
    @test min(a, 3) == min(3, a) == 3
    @test max(a, 4) == max(4, a) == 4

    @test typeof(min(a, 4)) <: Value
    @test typeof(max(a, 3)) <: Value

    @test max(a) == a.x + a.σ
    @test min(a) == a.x - a.σ

    @test min(a, b, c, -2, 10) == -2
    @test min(a, b, c, 34, 13, 9) == b
    @test max(a, b, c, 2, 10) == c
    @test max(a, 4, b, 2, 1) == 4
end

@testset "val/unc" begin
    @test val(a.x) == a.x
    @test unc(a.x) == 0

    @test val(a) == a.x
    @test unc(a) == a.σ

    A = fill(Value(a.x, a.σ), (5, 2, 6))

    @test eltype(A) <: Value
    @test eltype(val(A)) <: Real
    @test mean(val(A)) ≈ a.x
    @test mean(unc(A)) ≈ a.σ
end

@testset "set_unc" begin
    A = zeros(3, 7, 2)

    uncA = set_unc(A, 0.2)
    @test eltype(uncA) <: Value
    @test mean(unc(uncA)) ≈ 0.2

    uncA = set_unc(A, fill(0.3, (3, 7, 2)))
    @test eltype(uncA) <: Value
    @test mean(unc(uncA)) ≈ 0.3

    A = fill(Value(a.x, a.σ), (5, 2, 6))
    uncA = set_unc(A, 0.2)
    @test eltype(uncA) <: Value
    @test mean(unc(uncA)) ≈ 0.2

    uncA = set_unc(A, fill(0.3, (5, 2, 6)))
    @test eltype(uncA) <: Value
    @test mean(unc(uncA)) ≈ 0.3
end