text:
using OceanTurb, Printf

using OceanTurb.ModularKPP: DiagnosticPlumeModel

@use_pyplot_utils

usecmbright()

function makeplot!(axs, model)
    ModularKPP.update_state!(model)

    ΔT = model.state.plume.T - model.solution.T

    ϵ = CellField(model.grid)
    for i in eachindex(ϵ)
        @inbounds ϵ[i] = ModularKPP.entrainment(model.grid.zc[i], model)
    end

    for ax in axs
        sca(ax); cla()
    end

    sca(axs[1]); xlabel(L"\breve T - T"); ylabel(L"z \, \mathrm{(m)}")
    sca(axs[2]); xlabel(L"T")
    sca(axs[3]); xlabel(L"\breve W^2")
    sca(axs[4]); xlabel(L"\epsilon")

    sca(axs[1])
    plot(ΔT)
    removespines("top", "right")

    sca(axs[2])
    plot(model.solution.T, label="Environment \$ T \$")
    plot(model.state.plume.T, "--", label="Plume \$ \\breve T \$")
    removespines("top", "left", "right")

    sca(axs[3])
    plot(model.state.plume.W²)
    removespines("top", "left", "right")

    sca(axs[4])
    plot(ϵ)
    removespines("top", "left", "right")

    axs[2].tick_params(left=false, labelleft=false)
    axs[3].tick_params(left=false, labelleft=false)
    axs[4].tick_params(left=false, labelleft=false)

    return nothing
end

Qb = 1e-7
N² = 1e-5
Δt = 1minute

# The standard setup except with a plume model rather than a counter-gradient flux model.
model = ModularKPP.Model(;         grid = UniformGrid(N=64, H=256),
                                stepper = :BackwardEuler,
                           nonlocalflux = DiagnosticPlumeModel(Ca=0.1, Ce=0.4, Cα=1.0))

# Initial condition and fluxes
dTdz = N² / (model.constants.α * model.constants.g)
T₀(z) = 20 + dTdz*z
model.solution.T = T₀
OceanTurb.set!(model.state.plume.T, T₀)

Qθ = Qb / (model.constants.α * model.constants.g)

model.bcs.T.top = FluxBoundaryCondition(Qθ)
model.bcs.T.bottom = GradientBoundaryCondition(dTdz)

fig, axs = subplots(ncols=4, sharey=true)

run_until!(model, Δt, 24hour)
makeplot!(axs, model)


text:
using Documenter
using GeometricIntegrators

makedocs(;
    authors="Michael Kraus",
    sitename="Symplectic Runge-Kutta Methods for Degenerate Lagrangian Systems",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://michakraus.github.io/papers-srk-methods-for-degenerate-lagrangian-systems",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",

        "Lotka-Volterra 2d (singular Lagrangian)" => [
            "Symplectic Gauss-Legendre Runge-Kutta Methods" => "lotka-volterra-2d-singular/lotka-volterra-2d-singular-srk.md",
            "Gauss-Legendre Runge-Kutta Methods" => "lotka-volterra-2d-singular/lotka-volterra-2d-singular-firk.md",
            "Gauss-Legendre VPRK Methods" => "lotka-volterra-2d-singular/lotka-volterra-2d-singular-vprk-gauss.md",
            "Lobatto IIIA-IIIB VPRK Methods" => "lotka-volterra-2d-singular/lotka-volterra-2d-singular-vprk-lobatto-ab.md",
            "Lobatto IIIB-IIIA VPRK Methods" => "lotka-volterra-2d-singular/lotka-volterra-2d-singular-vprk-lobatto-ba.md",
            "Radau VPRK Methods" => "lotka-volterra-2d-singular/lotka-volterra-2d-singular-vprk-radau.md",
        ],

        "Lotka-Volterra 2d (symmetric Lagrangian)" => [
            "Symplectic Gauss-Legendre Runge-Kutta Methods" => "lotka-volterra-2d-symmetric/lotka-volterra-2d-symmetric-srk.md",
            "Gauss-Legendre VPRK Methods" => "lotka-volterra-2d-symmetric/lotka-volterra-2d-symmetric-vprk-gauss.md",
            "Gauss-Legendre Runge-Kutta Methods" => "lotka-volterra-2d-symmetric/lotka-volterra-2d-symmetric-firk.md",
            "Lobatto IIIA-IIIB VPRK Methods" => "lotka-volterra-2d-symmetric/lotka-volterra-2d-symmetric-vprk-lobatto-ab.md",
            "Lobatto IIIB-IIIA VPRK Methods" => "lotka-volterra-2d-symmetric/lotka-volterra-2d-symmetric-vprk-lobatto-ba.md",
            "Radau VPRK Methods" => "lotka-volterra-2d-symmetric/lotka-volterra-2d-symmetric-vprk-radau.md",
        ],
    ],
)

deploydocs(;
    repo="github.com/michakraus/paper-srk-methods-for-degenerate-lagrangian-systems",
    devbranch="main"
)

ext: jl


text:
using DataStructures: nlargest

struct GlobalPool{A}
    aggr::Symbol
    cluster::A
    function GlobalPool(aggr::Symbol, dim...)
        cluster = ones(Int64, dim)
        new{typeof(cluster)}(aggr, cluster)
    end
end

(g::GlobalPool)(X::AbstractArray) = pool(g.aggr, g.cluster, X)

struct LocalPool{A}
    aggr::Symbol
    cluster::A
end

function LocalPool(aggr::Symbol, cluster::AbstractArray)
    LocalPool{typeof(cluster)}(aggr, cluster)
end

(l::LocalPool)(X::AbstractArray) = pool(l.aggr, l.cluster, X)

struct TopKPool{T,S}
    A::AbstractMatrix{T}
    k::Int
    p::AbstractVector{S}
    Ã::AbstractMatrix{T}
end

function TopKPool(adj::AbstractMatrix, k::Int, in_channel::Integer; init=glorot_uniform)
    TopKPool(adj, k, init(in_channel), similar(adj, k, k))
end

function (t::TopKPool)(X::AbstractArray)
    y = t.p' * X / norm(t.p)
    idx = topk_index(y, t.k)
    t.Ã .= view(t.A, idx, idx)
    X_ = view(X, :, idx) .* σ.(view(y, idx)')
    return X_
end


text:
module CosmoMMF

using LinearAlgebra
using FFTW, AbstractFFTs
using StaticArrays
using YAML
using JLD2
using Strided
using Images

include("util.jl")
include("filter.jl")
include("tagging.jl")
# include("../test/runtests.jl")

end


text:
module Verbose_print
    import InteractiveUtils
    abstract type Verbose_level end

    struct Verbose_3  <: Verbose_level
        fp::Union{Nothing,IOStream}
        Verbose_3() = new(nothing)
        Verbose_3(filename::String) = new(open(filename,"w"))
        Verbose_3(fp::IOStream) = new(fp)
    end

    struct Verbose_2  <: Verbose_level
        fp::Union{Nothing,IOStream}
        Verbose_2() = new(nothing)
        Verbose_2(filename::String) = new(open(filename,"w"))
        Verbose_2(fp::IOStream) = new(fp)
    end

    struct Verbose_1  <: Verbose_level
        fp::Union{Nothing,IOStream}
        Verbose_1() = new(nothing)
        Verbose_1(filename::String) = new(open(filename,"w"))
        Verbose_1(fp::IOStream) = new(fp)
    end

    function Base.flush(v::Verbose_level)
        if v.fp != nothing
            flush(v.fp)
        end
    end

    function InteractiveUtils.versioninfo(v::Verbose_level)
        InteractiveUtils.versioninfo()
        if v.fp != nothing
            InteractiveUtils.versioninfo(v.fp)
        end
    end

    function println_verbose1(v::Verbose_3,val...)
        println(val...)
        if v.fp != nothing
            println(v.fp,val...)
        end
    end

    function println_verbose2(v::Verbose_3,val...)
        println(val...)
        if v.fp != nothing
            println(v.fp,val...)
        end
    end

    function println_verbose3(v::Verbose_3,val...)
        println(val...)
        if v.fp != nothing
            println(v.fp,val...)
        end
    end

    function println_verbose1(v::Verbose_2,val...)
        println(val...)
        if v.fp != nothing
            println(v.fp,val...)
        end
    end

    function println_verbose2(v::Verbose_2,val...)
        println(val...)
        if v.fp != nothing
            println(v.fp,val...)
        end
    end

    function println_verbose3(v::Verbose_2,val...)
        return
    end

    function println_verbose1(v::Verbose_1,val...)
        println(val...)
        if v.fp != nothing
            println(v.fp,val...)
        end
    end

    function println_verbose2(v::Verbose_1,val...)
        return
    end

    function println_verbose3(v::Verbose_1,val...)
        return
    end

    function print_verbose1(v::Verbose_3,val...)
        print(val...)
        if v.fp != nothing
            print(v.fp,val...)
        end
    end

    function print_verbose2(v::Verbose_3,val...)
        print(val...)
        if v.fp != nothing
            print(v.fp,val...)
        end
    end

    function print_verbose3(v::Verbose_3,val...)
        print(val...)
        if v.fp != nothing
            print(v.fp,val...)
        end
    end

    function print_verbose1(v::Verbose_2,val...)
        print(val...)
        if v.fp != nothing
            print(v.fp,val...)
        end
    end

    function print_verbose2(v::Verbose_2,val...)
        print(val...)
        if v.fp != nothing
            print(v.fp,val...)
        end
    end

    function print_verbose3(v::Verbose_2,val...)
        return
    end

    function print_verbose1(v::Verbose_1,val...)
        print(val...)
        if v.fp != nothing
            print(v.fp,val...)
        end
    end

    function print_verbose2(v::Verbose_1,val...)
        return
    end

    function print_verbose3(v::Verbose_1,val...)
        return
    end
end


text:
export heatmap!, plot_percolation_prob, closeall

"display labeled square lattice "
function heatmap!(latt::TwoDLattice)
    if !latt.lattice_properties.islabeled
        label_components!(latt)
    end

    labeled_site = latt.lattice_properties.labeled_lattice_sites
    linsize = latt.lattice_config.linear_size
    maxlabel = maximum(labeled_site)
    colorchange = Dict(i=>j for (i,j) in zip(1:maxlabel, shuffle(1:maxlabel)))

    row, col = size(labeled_site)
    site = zeros(row, col)
    @simd for j in 1:col
        @simd for i in 1:row
            if labeled_site[i,j] != 0
                @inbounds site[i,j] = colorchange[labeled_site[i,j]]
            end
        end
    end

    heatmap(site, colorbar=false, aspect_ratio=:equal, c=:lighttest)
end

"""
plot_percolation_prob(latticetype::String, linsize::Int, ps, pinc, pf, nsample::Int)

plot percolation probability vs occupied probability.
"""
function plot_percolation_prob(latticetype::String, linsize::Int, ps, pinc, pf, nsample::Int)
    prob = ps:pinc:pf
    percoprob = zeros(length(prob))

    for (i,p) in enumerate(prob)
        nperco = 0
        for iter in 1:nsample
            if latticetype == "square"
                l = Square(linsize, p)
            elseif latticetype == "triangular"
                l = Triangular(linsize, p)
            else
                l = Honeycomb(linsize, p)
            end
            nperco += ispercolation!(l)
        end
        percoprob[i] = nperco / nsample
    end

    plot(prob, percoprob,
         xlabel="Concentration",
         ylabel="Percolation Probability",
         title="Percolation Probability",
         legend=false, marker=:circle)
end


text:
function plotKDE(fgl::FactorGraph, sym::Symbol;
      dims=nothing,
      title="",
      levels::Int=5,
      fill::Bool=false,
      api::DataLayerAPI=dlapi  )
  #
  p = getVertKDE(fgl,sym, api=api)
  # mmarg = length(marg) > 0 ? marg : collect(1:Ndim(p))
  # mp = marginal(p,mmarg)
  plotKDE(p, levels=levels, dims=dims, title=string(sym, "  ", title), fill=fill )
end
function plotKDE(fgl::FactorGraph, syms::Vector{Symbol};
      addt::Vector{BallTreeDensity}=BallTreeDensity[],
      dims=nothing,
      title=nothing,
      levels=3,
      api::DataLayerAPI=dlapi  )
  #
  COLORS = ["black";"red";"green";"blue";"cyan";"deepskyblue"]
  MP = BallTreeDensity[]
  LEG = String[]
  # mmarg = Int[]
  for sym in syms
    p = getVertKDE(fgl,sym, api=api)
    # mmarg = length(marg) > 0 ? marg : collect(1:Ndim(p))
    # mp = marginal(p,mmarg)
    push!(MP, p)
    push!(LEG, string(sym))
  end
  for p in addt
    # mp = marginal(p,mmarg)
    push!(MP, p)
    push!(LEG, "add")
  end
  plotKDE(MP,c=COLORS[1:length(MP)], levels=levels, dims=dims, legend=LEG, title=title)
end
# function plotKDE(fgl::FactorGraph, sym::Symbol;
#       marg=nothing,
#       levels::Int=5  )
#   #
#   warn("Depricated call, use dims=Int[..] as keyword instead.")
#   plotKDE(fgl, sym, dims=marg, levels=levels  )
# end
# function plotKDE(fgl::FactorGraph, syms::Vector{Symbol};
#       addt::Vector{BallTreeDensity}=BallTreeDensity[],
#       marg=nothing,
#       levels=3  )
#   #
#   warn("Depricated call, use dims=Int[..] as keyword instead.")
#   plotKDE(fgl, syms, dims=marg, addt=addt, levels=levels  )
# end

"""
    plotKDEofnc(fg,fsym)

plot absolute values of variables and measurement surrounding fsym factor.
"""
function plotKDEofnc(fgl::FactorGraph, fsym::Symbol;
    marg=nothing,
    N::Int=100  )
  #
  fnc = nothing
  if haskey(fgl.fIDs, fsym)
    fnc = getfnctype( fgl, fgl.fIDs[fsym] )
  else
    error("fIDs doesn't have $(fsym)")
  end
  p = kde!( getSample(fnc, N)[1]  )
  # mmarg = length(marg) > 0 ? marg : collect(1:Ndim(p))
  plotKDE(fgl, lsf(fgl, fsym), addt=[p], marg=marg)
end

"""
    plotKDEresiduals(fg,fsym)

Trye plot relative values of variables and measurement surrounding fsym factor.
"""
function plotKDEresiduals(fgl::FactorGraph,
      fsym::Symbol;
      N::Int=100,
      levels::Int=3,
      dims=nothing,
      fill=false,
      api::DataLayerAPI=localapi  )
  #
  COLORS = ["black";"red";"green";"blue";"cyan";"deepskyblue"]
  fnc = getfnctype( fgl, fgl.fIDs[fsym] )
  # @show sxi = lsf(fgl, fsym)[1]
  # @show sxj = lsf(fgl, fsym)[2]
  fct = getVert(fgl, fsym, nt=:fnc, api=api)
  @show sxi = getData(fct).fncargvID[1]
  @show sxj = getData(fct).fncargvID[2]
  xi = getVal(fgl, sxi, api=api)
  xj = getVal(fgl, sxj, api=api)
  measM = getSample(fnc, N)
  meas = length(measM) == 1 ? (0*measM[1], ) : (0*measM[1], measM[2])
  @show size(meas[1])
  d = size(measM[1],1)
  RES = zeros(d,N)
  for i in 1:N
    res = zeros(d)
    fnc(res, i, meas, xi, xj)
    RES[:,i] = res
    if length(measM) > 1
      if measM[2][i] == 0
        RES[:,i] = 0.5*randn(d)
      end
    end
  end
  pR = kde!(RES)
  pM = kde!(measM[1])
  fncvar = getfnctype(fct)
  plotKDE([pR;pM], c=COLORS[1:2], dims=dims,levels=3, legend=["residual";"model"], fill=fill, title=string(fsym, ", ", fncvar))
end

"""
    plotPriorsAtCliq(treel, lb, cllb)

Plot the product of priors and incoming upward messages for variable in clique.

plotPriorsAtCliq(tree, :x2, :x1[, marg=[1;2]])
"""
function plotPriorsAtCliq(tree::BayesTree, lb::Symbol, cllb::Symbol;
        dims::Vector{Int}=Int[],
        levels::Int=1,
        fill::Bool=false  )
  #
  COLORS = ["black";"red";"green";"blue";"cyan";"deepskyblue"]

  cliq = whichCliq(tree, cllb)
  cliqprs = cliq.attributes["debug"].priorprods[1]

  vidx = 1
  for lbel in cliqprs.lbls
    if lbel == lb
      break;
    else
      vidx += 1
    end
  end
  dims = length(dims)>0 ? dims : collect(1:size(cliqprs.prods[vidx].prev,1))

  arr = BallTreeDensity[]
  push!(arr, marginal(kde!(cliqprs.prods[vidx].prev), dims)  )
  push!(arr, marginal(kde!(cliqprs.prods[vidx].product), dims)  )
  len = length(cliqprs.prods[vidx].potentials)
  lg = String["p";"n"]
  i=0
  for pot in cliqprs.prods[vidx].potentials
    push!(arr, marginal(pot, dims)  )
    i+=1
    push!(lg, cliqprs.prods[vidx].potentialfac[i])
  end
  cc = COLORS[1:(len+2)]
  # @show length(arr), length(cc), length(lg)
  plotKDE(arr, c=cc, legend=lg, levels=levels, fill=fill );
end

"""
    plotUpMsgsAtCliq(treel, lb, cllb)

Draw the up pass belief of lb from clique cllb.

plotUpMsgsAtCliq(tree, :x2, :x1)
"""
function plotUpMsgsAtCliq(treel::BayesTree, lb::Symbol, cllb::Symbol;
      show::Bool=true,
      w=20cm, h=15cm,
      levels::Int=1,
      marg::Vector{Int}=Int[] )
  #
  cliq = whichCliq(treel, string(cllb))
  cliqoutmsg = cliq.attributes["debug"].outmsg
  lbls = cliq.attributes["debug"].outmsglbls

  bel = kde!(cliqoutmsg.p[lbls[lb]])
  bel = length(marg)==0 ? bel : marginal(bel, marg)

  plotKDE(bel)
end

function plotMCMC(treel::BayesTree, lbll::Symbol;
      delay::Int=200,
      show::Bool=true,
      w=20cm, h=15cm,
      levels::Int=1,
      dims=nothing  )
  #
  cliq = whichCliq(treel, string(lbll))
  cliqdbg = cliq.attributes["debug"]

  vidx = 1
  for lb in cliqdbg.mcmc[1].lbls
    if lb == lbll
      break;
    else
      vidx += 1
    end
  end

  tmpfilepath = joinpath(dirname(@__FILE__),"tmpimgs")
  ARR = BallTreeDensity[]
  COLORS = ["black";"red";"green";"blue";"cyan";"deepskyblue";"magenta"]
  for i in 1:length(cliqdbg.mcmc)
    ppr = kde!(cliqdbg.mcmc[i].prods[vidx].prev)
    ppp = kde!(cliqdbg.mcmc[i].prods[vidx].product)
    ARR = [ARR;ppr;ppr;cliqdbg.mcmc[i].prods[vidx].potentials]
  end
  rangeV = getKDERange(ARR)
  ppp = nothing
  for i in 1:length(cliqdbg.mcmc)
    ppr = kde!(cliqdbg.mcmc[i].prods[vidx].prev)
    ppp = kde!(cliqdbg.mcmc[i].prods[vidx].product)
    arr = [ppr;ppp;cliqdbg.mcmc[i].prods[vidx].potentials]
    len = length(cliqdbg.mcmc[i].prods[vidx].potentials)
    lg = String["p";"n";cliqdbg.mcmc[i].prods[vidx].potentialfac] #map(string, 1:len)]
    cc = plotKDE(arr, c=COLORS[1:(len+2)], legend=lg, levels=levels, fill=true, axis=rangeV, dims=dims );
    Gadfly.draw(PNG(joinpath(tmpfilepath,"$(string(lbll))mcmc$(i).png"),w,h),cc)
  end
  # draw initial and final result
  pp0 = kde!(cliqdbg.mcmc[1].prods[vidx].prev)
  i = 0
  cc = plotKDE([pp0], c=[COLORS[1]], legend=["0"], levels=levels, fill=true, axis=rangeV, dims=dims );
  Gadfly.draw(PNG(joinpath(tmpfilepath,"$(string(lbll))mcmc$(i).png"),w,h),cc)

  i = length(cliqdbg.mcmc)+1
  cc = plotKDE([pp0;ppp], c=COLORS[1:2], legend=["0";"n"], levels=levels, fill=true, axis=rangeV, dims=dims );
  Gadfly.draw(PNG(joinpath(tmpfilepath,"$(string(lbll))mcmc$(i).png"),w,h),cc)
  # generate output
  run(`convert -delay $(delay) $(tmpfilepath)/$(string(lbll))mcmc*.png $(tmpfilepath)/$(string(lbll))mcmc.gif`)
  !show ? nothing : (@async run(`eog $(tmpfilepath)/$(string(lbll))mcmc.gif`) )
  return "$(tmpfilepath)/$(string(lbll))mcmc.gif"
end

function drawOneMC!(cliqMC::CliqGibbsMC, minmax, mcmc=0; offs=2.0)

    if mcmc>minmax[4]
      minmax[4]=mcmc
    end

    i = 0.0
    for prod in cliqMC.prods
        prodVal = kde!(prod.product,"lcv")
        plotKDEProd!([prodVal;prod.potentials],minmax, h=-i*offs, mcmc=mcmc)
        i += 1.0
    end

end

function drawMCMCDebug(cliq; offs=2.0)
    println("$(cliq.attributes["label"])")
    cliqDbg = cliq.attributes["data"].debug
    MCd = 100.0/length(cliqDbg.mcmc)
    MCo=0.0
    minmax=[99999,-99999,-10,-99999.0]
    for mc in cliqDbg.mcmc
        drawOneMC!(mc, minmax, MCo, offs=offs)
        MCo+=MCd
    end

    n = 3
    y = linspace(minmax[1], minmax[2], n)
    x = linspace(minmax[3],minmax[4],n)
    xgrid = repmat(x',n,1)
    ygrid = repmat(y,1,n)
    z = zeros(n,n)

    for i in 1:length(cliqDbg.mcmc[1].prods)
      surf(xgrid,ygrid,z-(i-1)*offs,alpha=0.04, linewidth=0.0)
    end
end

function drawTreeUpwardMsgs(fgl::FactorGraph, bt::BayesTree; N=300)
    len = length(bt.cliques)-1
    vv = Array{Gadfly.Compose.Context,1}(len)
    #r = Array{RemoteRef,1}(len)
    i = 0
    for cliq in bt.cliques
        if cliq[1] == 1 println("No upward msgs from root."); continue; end
        @show cliq[2].attributes["label"]
        i+=1
        vv[i] = drawHorDens(fgl, cliq[2].attributes["data"].debug.outmsg.p, N)
    end
    #[r[j] = @spawn drawCliqueMsgs(bt.cliques[j+1]) for j in 1:len]
    #[vv[j] = fetch(r[j]) for j in 1:len]
    vv
end

function drawFrontalDens(fg::FactorGraph, bt::BayesTree;
                        N=300, gt=Union{})
    len = length(bt.cliques)
    vv = Array{Gadfly.Compose.Context,1}(len)
    i = 0
    for cliq in bt.cliques
        #@show cliq[2].attributes["label"]
        lenfr = length(cliq[2].attributes["data"].frontalIDs)

        p = Array{BallTreeDensity,1}(lenfr)
        j=0
        #pvals = Array{Array{Float64,2},1}(lenfr)
        gtvals = Dict{Int,Array{Float64,2}}()
        lbls = String[]

        for frid in cliq[2].attributes["data"].frontalIDs
            j+=1
            p[j] = getVertKDE(fg, frid) # getKDE(fg.v[frid])
            # p[j] = kde!(fg.v[frid].attributes["val"])

            #pvals[j] = fg.v[frid].attributes["val"]

            if gt!=Union{}
              gtvals[j] = gt[dlapi.getvertex(fg,frid).attributes["label"]] # fg.v[frid].
              #push!(gtvals, gt[fg.v[frid].attributes["label"]][1])
              #push!(gtvals, gt[fg.v[frid].attributes["label"]][2])
            end
            push!(lbls, dlapi.getvertex(fg,frid).attributes["label"]) # fg.v[frid].

        end

        #r = Array{RemoteRef,1}(lenfr)
        #[r[j] = @spawn kde!(pvals[j]) for j in 1:lenfr]
        #[p[j] = fetch(r[j]) for j in 1:lenfr]

        i+=1
        if length(gtvals) > 0
          #gtvals = reshape(gtvals,2,round(Int,length(gtvals)/2))'
          vv[i] = drawHorDens(p, N, gt=gtvals, lbls=lbls)
        else
          vv[i] = drawHorDens(p, N,lbls=lbls)
        end
    end
    #
    vv
end

# for some reason we still need msgPlots of right size in the global for these functions to work.
# precall drawTreeUpwardMsgs or drawFrontalDens to make this work properly TODO
function vstackedDensities(msgPlots)
    #msgPlots = f(fg, bt) # drawTreeUpwardMsgs
    evalstr = ""
    for i in 1:length(msgPlots)
        evalstr = string(evalstr, ",msgPlots[$(i)]")
    end
    eval(parse(string("vstack(",evalstr[2:end],")")))
end

function investigateMultidimKDE(p::BallTreeDensity, p0::BallTreeDensity)
    co = ["black"; "blue"]
    h = Union{}
    x = plotKDE([marginal(p,[1]); marginal(p0,[1])], c=co )
    y = plotKDE([marginal(p,[2]); marginal(p0,[2])], c=co )
    if p.bt.dims >= 3
      th = plotKDE([marginal(p,[3]); marginal(p0,[3])], c=co )
      h = hstack(x,y,th)
    else
      h = hstack(x,y)
    end

    return h
end

function investigateMultidimKDE(p::Array{BallTreeDensity,1})
    co = ["black"; "blue"; "green"; "red"; "magenta"; "cyan"; "cyan1"; "cyan2";
    "magenta"; "cyan"; "cyan1"; "cyan2"; "magenta"; "cyan"; "cyan1"; "cyan2"; "magenta";
    "cyan"; "cyan1"; "cyan2"; "magenta"; "cyan"; "cyan1"; "cyan2"]
    # compute all the marginals
    Pm = Array{Array{BallTreeDensity,1},1}()
    push!(Pm,stackMarginals(p,1)) #[marginal(p[1],[1]); marginal(p[2],[1])]
    push!(Pm,stackMarginals(p,2)) #[marginal(p[1],[2]); marginal(p[2],[2])]

    h = Union{}
    x = plotKDE(Pm[1], c=co )
    y = plotKDE(Pm[2], c=co )
    if p[1].bt.dims >= 3
      #Pm3 = [marginal(p[1],[3]); marginal(p[2],[3])]
      push!(Pm,stackMarginals(p,3)) # [marginal(p[1],[3]); marginal(p[2],[3])]
      th = plotKDE(Pm[3], c=co )
      h = hstack(x,y,th)
    else
      h = hstack(x,y)
    end
    return h
end

function investigateMultidimKDE(p::BallTreeDensity)
    x = plotKDE(marginal(p,[1]) )
    y = plotKDE(marginal(p,[2]) )
    if p.bt.dims >= 3
      th = plotKDE(marginal(p,[3]) )
      return hstack(x,y,th)
    end
    return hstack(x,y)
end

function vArrPotentials(potens::Dict{Int,EasyMessage})
  vv = Array{Gadfly.Compose.Context,1}(length(potens))
  i = 0
  oned=false
  for p in potens
      i+=1
      pb = kde!(p[2].pts, p[2].bws)
      if size(p[2].pts,1) > 3
        # vv[i] = plotKDE(pb)
        error("can't handle higher dimensional plots here yet")
      elseif size(p[2].pts,1) > 1
        vv[i] = investigateMultidimKDE(pb)
      else
        vv[i] = plotKDE(pb)
      end
  end
  return vv
end

function draw(em::EasyMessage;xlbl="X")
  p = Union{}
  if size(em.pts,1) == 1
    p=plotKDE(kde!(em),xlbl=xlbl)
  else
    p=plotKDE(kde!(em))
  end
  return p
end

function whosWith(cliq::Graphs.ExVertex)
  println("$(cliq.attributes["label"])")
  for pot in cliq.attributes["data"].potentials
      println("$(pot)")
  end
  nothing
end

function whosWith(bt::BayesTree, frt::String)
    whosWith(whichCliq(bt,frt))
end

function drawUpMsgAtCliq(fg::FactorGraph, cliq::Graphs.ExVertex)
    for id in keys(cliq.attributes["data"].debug.outmsg.p)
        print("$(dlapi.getvertex(fg,id).attributes["label"]), ") #fg.v[id].
    end
    println("")
    sleep(0.1)
    potens = cliq.attributes["data"].debug.outmsg.p
    vArrPotentials(potens)
end

function drawUpMsgAtCliq(fg::FactorGraph, bt::BayesTree, lbl::String)
    drawUpMsgAtCliq(fg, whichCliq(bt, lbl) )
end

# function drawDwnMsgAtCliq(fg::FactorGraph, cliq::Graphs.ExVertex)
function dwnMsgsAtCliq(fg::FactorGraph, cliq::Graphs.ExVertex)
  for id in keys(cliq.attributes["data"].debugDwn.outmsg.p)
      print("$(dlapi.getvertex(fg,id).label), ") # fg.v[id].
  end
  println("")
  sleep(0.1)
  potens = cliq.attributes["data"].debugDwn.outmsg.p
  potens
end

function dwnMsgsAtCliq(fg::FactorGraph, bt::BayesTree, lbl::String)
    dwnMsgsAtCliq(fg, whichCliq(bt, lbl) )
end

function drawPose2DMC!(plots::Array{Gadfly.Compose.Context,1}, cliqMC::CliqGibbsMC)

    for prod in cliqMC.prods
        prodVal = kde!(prod.product,"lcv") #cliqMC.prods[1]
        push!(plots, plotKDE([prodVal;prod.potentials]) )
    end
    vstackedPlots(plots)
end

function mcmcPose2D!(plots::Array{Gadfly.Compose.Context,1}, cliqDbg::DebugCliqMCMC, iter::Int=1)
    # for mc in cliqDbg.mcmc
    mc = cliqDbg.mcmc[iter]
    v = drawPose2DMC!(plots, mc)
    # end
    return v
end

function drawUpMCMCPose2D!(plots::Array{Gadfly.Compose.Context,1}, cliq::Graphs.ExVertex, iter::Int=1)
    whosWith(cliq)
    cliqDbg = cliq.attributes["data"].debug
    sleep(0.1)
    mcmcPose2D!(plots, cliqDbg, iter)
end

function drawUpMCMCPose2D!(plots::Array{Gadfly.Compose.Context,1}, bt::BayesTree, frt::String, iter::Int=1)
    drawUpMCMCPose2D!(plots, whichCliq(bt,frt), iter)
end

function drawDwnMCMCPose2D!(plots::Array{Gadfly.Compose.Context,1}, cliq::Graphs.ExVertex, iter::Int=1)
    whosWith(cliq)
    cliqDbg = cliq.attributes["data"].debugDwn
    sleep(0.1)
    mcmcPose2D!(plots, cliqDbg, iter)
end

function drawDwnMCMCPose2D!(plots::Array{Gadfly.Compose.Context,1}, bt::BayesTree, frt::String, iter::Int=1)
    drawDwnMCMCPose2D!(plots, whichCliq(bt,frt), iter)
end

function drawLbl(fgl::FactorGraph, lbl::Symbol)
    # v = dlapi.getvertex(fgl,lbl)
    # plotKDE(kde!(getVal(v)))
    plotKDE(getVertKDE(fgl,lbl))
end
drawLbl{T <: AbstractString}(fgl::FactorGraph, lbl::T) = drawLbl(fgl, Symbol(lbl))

function predCurrFactorBeliefs(fgl::FactorGraph, fc::Graphs.ExVertex)
  # TODO update to use ls and lsv functions
  prjcurvals = Dict{String, Array{BallTreeDensity,1}}()
  for v in dlapi.outneighbors(fgl, fc)
    pred = kde!(evalFactor2(fgl, fc, v.index))
    curr = kde!(getVal(v))
    prjcurvals[v.attributes["label"]] = [curr; pred]
  end
  return prjcurvals, collect(keys(prjcurvals))
end

function drawHorDens(fgl::FactorGraph, pDens::Dict{Int,EasyMessage}, N=200)
  p = BallTreeDensity[]
  lbls = String[]
  for pd in pDens
    push!(p, kde!(pd[2].pts,pd[2].bws))
    push!(lbls, dlapi.getvertex(fgl,pd[1]).attributes["label"])
  end
  @show lbls
  drawHorDens(p,N,lbls=lbls)
end

function drawHorBeliefsList(fgl::FactorGraph, lbls::Array{Symbol,1};
                        nhor::Int=-1,gt=nothing,N::Int=200, extend=0.1)
  len = length(lbls)
  pDens = BallTreeDensity[]
  for lb in lbls
    ptkde = getVertKDE(fgl,lb)
    push!(pDens, ptkde )
  end

  if nhor<1
    nhor = round(Int,sqrt(len))
  end
  vlen = ceil(Int, len/nhor)
  vv = Array{Gadfly.Compose.Context,1}(vlen)
  conslb = deepcopy(lbls)
  vidx = 0
  for i in 1:nhor:len
    pH = BallTreeDensity[]
    gtvals = Dict{Int,Array{Float64,2}}()
    labels = String[]
    for j in 0:(nhor-1)
      if i+j <= len
        push!(pH, pDens[i+j])
        push!(labels, string(lbls[i+j]))
        if gt != nothing gtvals[j+1] = gt[lbls[i+j]]  end
      end
    end
    vidx+=1
    if gt !=nothing
      vv[vidx] = KernelDensityEstimate.drawHorDens(pH, N=N, gt=gtvals, lbls=labels, extend=extend)
    else
      vv[vidx] = KernelDensityEstimate.drawHorDens(pH, N=N, lbls=labels, extend=extend)
    end
  end
  vv
end

function drawFactorBeliefs(fgl::FactorGraph, flbl::Symbol)
  if !haskey(fgl.fIDs, flbl)
    println("no key $(flbl)")
    return nothing
  end
  # for fc in fgl.f
    # if fc[2].attributes["label"] == flbl

    fc = fgl.g.vertices[fgl.fIDs[flbl]]  # fc = fgl.f[fgl.fIDs[flbl]]
      prjcurvals, lbls = predCurrFactorBeliefs(fgl, fc)
      if length(lbls) == 3
        return vstack(
        plotKDE(prjcurvals[lbls[1]]),
        plotKDE(prjcurvals[lbls[2]]),
        plotKDE(prjcurvals[lbls[3]]),
        )
      elseif length(lbls) == 2
        return vstack(
        plotKDE(prjcurvals[lbls[1]]),
        plotKDE(prjcurvals[lbls[2]]),
        )
      elseif length(lbls) == 1
        return plotKDE(prjcurvals[lbls[1]])
      end

    # end
  # end
  nothing
end
drawFactorBeliefs{T <: AbstractString}(fgl::FactorGraph, flbl::T) = drawFactorBeliefs(fgl, Symbol(flbl))

function getColorsByLength(len::Int)
  COLORS = ["black";"red";"green";"blue";"deepskyblue";"yellow";"magenta"]
  morecyan = ["cyan" for i in (length(COLORS)+1):len]
  retc = [COLORS; morecyan]
  return retc[1:len]
end

"""
    plotLocalProduct(fgl::FactorGraph, lbl::Symbol; N::Int=100, dims::Vector{Int}=Int[])

Plot the proposal belief from neighboring factors to `lbl` in the factor graph (ignoring Bayes tree representation),
and show with new product approximation for reference.
"""
function plotLocalProduct(fgl::FactorGraph, lbl::Symbol; N::Int=100, dims::Vector{Int}=Int[], api::DataLayerAPI=dlapi, levels::Int=1)
  warn("not showing partial constraints, but included in the product")
  arr = Array{BallTreeDensity,1}()
  lbls = String[]
  push!(arr, getVertKDE(fgl, lbl, api=api))
  push!(lbls, "curr")
  pp, parr, partials, lb = localProduct(fgl, lbl, N=N, api=api)
  if pp != parr[1]
    push!(arr,pp)
    push!(lbls, "prod")
    for a in parr
      push!(arr, a)
    end
    lbls = union(lbls, lb)
  end
  dims = length(dims) > 0 ? dims : collect(1:Ndim(pp))
  colors = getColorsByLength(length(arr))
  return plotKDE(arr, dims=dims, levels=levels, c=colors, legend=lbls, title=string("Local product, ",lbl))
end

"""
    plotLocalProduct{T <: AbstractString}(fgl::FactorGraph, lbl::T; N::Int=100, dims::Vector{Int}=Int[])

Plot the proposal belief from neighboring factors to `lbl` in the factor graph (ignoring Bayes tree representation),
and show with new product approximation for reference. String version is obsolete and will be deprecated.
"""
plotLocalProduct{T <: AbstractString}(fgl::FactorGraph, lbl::T; N::Int=100, dims::Vector{Int}=Int[]) = plotLocalProduct(fgl, Symbol(lbl), N=N, dims=dims)

function saveplot(pl;name="pl",frt=:png,w=25cm,h=25cm,nw=false,fill=true)
  if frt==:png
    Gadfly.draw(PNG(string(name,".png"),w,h),pl)
    # if fill run(`composite $(name).png plB.png $(name).png`) end
    if !nw run(`eog $(name).png`) end
  end
  if frt==:pdf
    Gadfly.draw(PDF(string(name,".pdf"),w,h),pl)
    if !nw run(`evince $(name).pdf`) end
  end
  nothing
end

function animateVertexBelief(FGL::Array{FactorGraph,1}, lbl; nw=false)
  len = length(FGL)
  [saveplot(plotLocalProduct(FG[i],lbl),h=15cm,w=30cm,name="gifs/pl$(i)",nw=true) for i=1:len];
  run(`convert -delay 100 gifs/pl*.png result.gif`)
  if !nw run(`eog result.gif`) end
  nothing
end

function fixRotWrapErr!(RT::Array{Float64,1})

  for i in 1:length(RT)
    if RT[i] > pi
      RT[i] = abs(RT[i]-2.0*pi)
    end
  end
  nothing
end

function asyncUniComp(fgl::FactorGraph, isamdict::Dict{Int,Array{Float64,1}})
  r,rt,lb = computeGraphResiduals(fgl,isamdict);
  fixRotWrapErr!(rt)
  return [sqrt(Base.mean(r.^2));maximum(abs(r));sqrt(Base.mean(rt.^2));maximum(rt)]
end

function unimodalCompare(FGL::Array{FactorGraph,1},isamdict::Dict{Int,Array{Float64,1}})
  len = length(FGL)
  RMS = Float64[]
  MAX = Float64[]
  RMSth = Float64[]
  MAXth = Float64[]

  rr = Future[] #RemoteRef[]

  for fgl in FGL
    push!(rr, remotecall(uppA(),asyncUniComp, fgl, isamdict))
  end

  for r in rr
    err = fetch(r)
    push!(RMS, err[1])
    push!(MAX, err[2])
    push!(RMSth, err[3])
    push!(MAXth, err[4])
  end

  x=0:(len-1)
  df1 = DataFrame(x=x, y=RMS, label="rms")
  df2 = DataFrame(x=x, y=MAX, label="max")
  df3 = DataFrame(x=x, y=RMSth*180.0/pi, label="rmsth")
  df4 = DataFrame(x=x, y=MAXth*180.0/pi, label="maxth")
  df = vcat(df1, df2)
  dfth = vcat(df3,df4)

  return df,dfth
end

function asyncAnalyzeSolution(fgl::FactorGraph, sym::Symbol)
  lbl = string(sym)
  pp, arr, partials = localProduct(fgl, lbl)
  lpm = getKDEMax(pp)
  em = getKDEMax(getVertKDE(fgl,lbl))
  err1 = norm(lpm[1:2]-em[1:2])
  err2 = 0.0
  if lbl[1]=='x'
    err2 = abs(lpm[3]-em[3])
  end
  return [err1;err2]
end

function analyzeSolution(FGL::Array{FactorGraph,1},fggt=Union{})
  len = length(FGL)
  RMS = Float64[]
  MAX = Float64[]
  RMSth = Float64[]
  MAXth = Float64[]
  for fgl in FGL
    xLB, lLB = ls(fgl)
    ERR = Float64[]
    ERRth = Float64[]
    ALB = [xLB;lLB]
    rr = Future[] #RemoteRef[]
    for lbl in ALB
      push!(rr, remotecall(uppA(),asyncAnalyzeSolution, fgl, lbl))
      # err
      # push!(ERR, err[1])
      # if lbl[1]=='x'
      #   push!(ERRth, err[2])
      # end
    end

    idx = 1
    for r in rr
      err = fetch(r)
      push!(ERR, err[1])
      if ALB[idx][1]=='x'
        push!(ERRth, err[2])
      end
      idx += 1
    end
    push!(RMS, sqrt(Base.mean(ERR.^2)))
    push!(MAX, maximum(abs(ERR)))
    push!(RMSth, sqrt(Base.mean(ERRth.^2)))
    push!(MAXth, maximum(ERRth))
  end

  x=0:(len-1)
  df1 = DataFrame(x=x, y=RMS, label="rms")
  df2 = DataFrame(x=x, y=MAX, label="max")
  df3 = DataFrame(x=x, y=RMSth*180.0/pi, label="rmsth")
  df4 = DataFrame(x=x, y=MAXth*180.0/pi, label="maxth")
  df = vcat(df1, df2)
  dfth = vcat(df3,df4)
  return df,dfth
end
# discrete_color_manual(colors...; levels=nothing,order=nothing) is deprecated, use color_discrete_manual(colors...; levels=levels,order=order) instead.

function drawAnalysis(df,dfth)
  return vstack(
  Gadfly.plot(df, x="x", y="y", color="label", Geom.line,
         Scale.discrete_color_manual("red","black")),
  Gadfly.plot(dfth, x="x", y="y", color="label", Geom.line,
        Scale.discrete_color_manual("red","black"))
        )
end

function getAllFGsKDEs(fgD::Array{FactorGraph,1}, vertid::Int)
  ret = Array{BallTreeDensity,1}()
  for i in 1:length(fgD)
    push!(ret, getVertKDE(fgD[i],vertid) )
  end
  return ret
end

function drawAllPose2DBeliefs(plots::Array{Gadfly.Compose.Context,1}, fgD::Array{FactorGraph,1})
    ids = sort(collect(keys(fgD[1].v)))
    co = ["black"; "blue"; "green"; "red"; "magenta"; "cyan"; "cyan1"; "cyan2"]
    println(co[1:length(fgD)])
    for i in ids
        @show dlapi.getvertex(fgD[1],i).attributes["label"] #fgD[1].v[i].
            kdes = getAllFGsKDEs(fgD, i)
            push!(plots, plotKDE(  kdes  )) # [kde!(getVal(V)); kde!(getVal(V0))]
    end
    vstackedPlots(plots)
end

# legacy function -- use the array version instead
function drawAllPose2DBeliefs(plots::Array{Gadfly.Compose.Context,1}, fgD::FactorGraph, fgD0=Union{})
  println("WARNING: drawAllPose2DBeliefs -- legacy function -- use the array version instead.")
  if fgD0 != Union{}
    drawAllPose2DBeliefs(plots, [fgD;fgD0])
  else
    drawAllPose2DBeliefs(plots, [fgD])
  end
end

function drawComicStripLM(fgD::Array{FactorGraph,1})
    comicA = Array{Gadfly.Plot,1}()
    for fgd in fgD
        cv = drawPosesLandms(fgd)
        # cv = drawPoses(fgd)
        push!(comicA,cv)
    end
    hstack(comicA)
end

function drawComicStrip(fgD::Array{FactorGraph,1})
    comicA = Array{Gadfly.Plot,1}()
    for fgd in fgD
        cv = drawPoses(fgd)
        push!(comicA,cv)
    end
    hstack(comicA)
end

function compositeComic(fnc::Function, fgGT, fgA::Array{FactorGraph,1})
    v = Union{}
    @show length(fgA)
    if length(fgA) == 2
        Gadfly.set_default_plot_size(25cm, 10cm)
        v = fnc([fgA[1:2];fgGT])
    elseif length(fgA) == 3
        Gadfly.set_default_plot_size(25cm, 20cm)
        v = vstack(fnc(fgA[1:2])
        ,fnc([fgA[3];fgGT])    )
    elseif length(fgA) == 4
        Gadfly.set_default_plot_size(25cm, 20cm)
        v = vstack(fnc(fgA[1:3])
        ,fnc([fgA[4];fgGT])    )
    elseif length(fgA) == 7
        Gadfly.set_default_plot_size(25cm, 25cm)
        v = vstack(fnc(fgA[1:3])
        ,fnc(fgA[4:6])
        ,fnc([fgA[7];fgGT])    )
    elseif length(fgA) == 10
        Gadfly.set_default_plot_size(25cm, 25cm)
        v = vstack(fnc(fgA[1:3])
        ,fnc(fgA[4:6])
        ,fnc(fgA[7:9])
        ,fnc([fgA[10];fgGT])    )
    elseif length(fgA) == 13
        Gadfly.set_default_plot_size(25cm, 30cm)
        v = vstack(fnc(fgA[1:3])
        ,fnc(fgA[4:6])
        ,fnc(fgA[7:9])
        ,fnc(fgA[10:12])
        ,fnc([fgA[13];fgGT])    )
    end
    v
end

function compositeComic(fnc::Function, fgA::Array{FactorGraph,1})
    v = Union{}
    @show length(fgA)
    if length(fgA) == 2
        Gadfly.set_default_plot_size(25cm, 10cm)
        v = fnc(fgA[1:2])
    elseif length(fgA) == 3
        Gadfly.set_default_plot_size(25cm, 20cm)
        v = fnc(fgA[1:3])
    elseif length(fgA) == 4
        Gadfly.set_default_plot_size(25cm, 25cm)
        v = vstack(fnc(fgA[1:2])
        ,fnc(fgA[3:4]))
    elseif length(fgA) == 6
        Gadfly.set_default_plot_size(25cm, 25cm)
        v = vstack(fnc(fgA[1:3])
        ,fnc(fgA[4:6])  )
    elseif length(fgA) == 9
        Gadfly.set_default_plot_size(25cm, 25cm)
        v = vstack(fnc(fgA[1:3])
        ,fnc(fgA[4:6])
        ,fnc(fgA[7:9])    )
    elseif length(fgA) == 12
        Gadfly.set_default_plot_size(25cm, 30cm)
        v = vstack(fnc(fgA[1:3])
        ,fnc(fgA[4:6])
        ,fnc(fgA[7:9])
        ,fnc(fgA[10:12])    )
    end
    v
end


text:
println("------------------------------------")
println("|     Fields and vector spaces     |")
println("------------------------------------")
ti = time()
@testset TimedTestSet "Fields" begin
    @test isa(ℝ, Field)
    @test isa(ℂ, Field)
    @test eval(Meta.parse(sprint(show, ℝ))) == ℝ
    @test eval(Meta.parse(sprint(show, ℂ))) == ℂ
    @test ℝ ⊆ ℝ
    @test ℝ ⊆ ℂ
    @test ℂ ⊆ ℂ
    @test !(ℂ ⊆ ℝ)

    for T in (Int8, Int16, Int32, Int64, BigInt)
        @test one(T) ∈ ℝ
        @test one(Rational{T}) ∈ ℝ
        @test !(one(Complex{T}) ∈ ℝ)
        @test !(one(Complex{Rational{T}}) ∈ ℝ)
        @test one(T) ∈ ℂ
        @test one(Rational{T}) ∈ ℂ
        @test one(Complex{T}) ∈ ℂ
        @test one(Complex{Rational{T}} ∈ ℂ)

        @test T ⊆ ℝ
        @test Rational{T} ⊆ ℝ
        @test !(Complex{T} ⊆ ℝ)
        @test !(Complex{Rational{T}} ⊆ ℝ)
        @test T ⊆ ℂ
        @test Rational{T} ⊆ ℂ
        @test Complex{T} ⊆ ℂ
        @test Complex{Rational{T}} ⊆ ℂ
    end
    for T in (Float32, Float64, BigFloat)
        @test one(T) ∈ ℝ
        @test !(one(Complex{T}) ∈ ℝ)
        @test one(T) ∈ ℂ
        @test one(Complex{T} ∈ ℂ)

        @test T ⊆ ℝ
        @test !(Complex{T} ⊆ ℝ)
        @test T ⊆ ℂ
        @test Complex{T} ⊆ ℂ
    end
end

@testset TimedTestSet "ElementarySpace: CartesianSpace" begin
    d = 2
    V = ℝ^d
    @test eval(Meta.parse(sprint(show, V))) == V
    @test eval(Meta.parse(sprint(show, typeof(V)))) == typeof(V)
    @test isa(V, VectorSpace)
    @test isa(V, ElementarySpace)
    @test isa(V, InnerProductSpace)
    @test isa(V, EuclideanSpace)
    @test isa(V, CartesianSpace)
    @test !isdual(V)
    @test !isdual(V')
    @test V == CartesianSpace(Trivial() => d) == CartesianSpace(Dict(Trivial() => d))
    @test @inferred(hash(V)) == hash(deepcopy(V))
    @test V == @inferred(dual(V)) == @inferred(conj(V)) == @inferred(adjoint(V))
    @test field(V) == ℝ
    @test @inferred(sectortype(V)) == Trivial
    @test @inferred(sectors(V)) == (Trivial(),)
    @test @inferred(TensorKit.hassector(V, Trivial()))
    @test @inferred(dim(V)) == d == @inferred(dim(V, Trivial()))
    @test dim(@inferred(typeof(V)())) == 0
    @test @inferred(TensorKit.axes(V)) == Base.OneTo(d)
    @test V == ℝ[d] == ℝ[](d) == typeof(V)(d)
    W = @inferred ℝ[1]
    @test @inferred(oneunit(V)) == W == oneunit(typeof(V))
    @test @inferred(⊕(V,V)) == ℝ^(2d)
    @test @inferred(⊕(V,oneunit(V))) == ℝ^(d+1)
    @test @inferred(⊕(V,V,V,V)) == ℝ^(4d)
    @test @inferred(fuse(V,V)) == ℝ^(d^2)
    @test @inferred(fuse(V,V',V,V')) == ℝ^(d^4)
    @test @inferred(flip(V)) == V'
    @test @inferred(min(V, ℝ^3)) == V
    @test @inferred(max(V', ℝ^3)) == ℝ^3
end

@testset TimedTestSet "ElementarySpace: ComplexSpace" begin
    d = 2
    V = ℂ^d
    @test eval(Meta.parse(sprint(show, V))) == V
    @test eval(Meta.parse(sprint(show, V'))) == V'
    @test eval(Meta.parse(sprint(show, typeof(V)))) == typeof(V)
    @test isa(V, VectorSpace)
    @test isa(V, ElementarySpace)
    @test isa(V, InnerProductSpace)
    @test isa(V, EuclideanSpace)
    @test isa(V, ComplexSpace)
    @test !isdual(V)
    @test isdual(V')
    @test V == ComplexSpace(Trivial() => d) == ComplexSpace(Dict(Trivial() => d))
    @test @inferred(hash(V)) == hash(deepcopy(V)) != hash(V')
    @test @inferred(dual(V)) == @inferred(conj(V)) == @inferred(adjoint(V)) != V
    @test @inferred(field(V)) == ℂ
    @test @inferred(sectortype(V)) == Trivial
    @test @inferred(sectors(V)) == (Trivial(),)
    @test @inferred(TensorKit.hassector(V, Trivial()))
    @test @inferred(dim(V)) == d == @inferred(dim(V, Trivial()))
    @test @inferred(TensorKit.axes(V)) == Base.OneTo(d)
    @test V == ℂ[d] == ℂ[](d) == typeof(V)(d)
    W = @inferred ℂ[1]
    @test @inferred(oneunit(V)) == W == oneunit(typeof(V))
    @test @inferred(⊕(V, V)) == ℂ^(2d)
    @test_throws SpaceMismatch (⊕(V, V'))
    @test_throws MethodError (⊕(ℝ^d, ℂ^d))
    @test_throws MethodError (⊗(ℝ^d, ℂ^d))
    @test @inferred(⊕(V,V)) == ℂ^(2d)
    @test @inferred(⊕(V,oneunit(V))) == ℂ^(d+1)
    @test @inferred(⊕(V,V,V,V)) == ℂ^(4d)
    @test @inferred(fuse(V,V)) == ℂ^(d^2)
    @test @inferred(fuse(V,V',V,V')) == ℂ^(d^4)
    @test @inferred(flip(V)) == V'
    @test @inferred(min(V, ℂ^3)) == V
    @test @inferred(max(V', ℂ[3]')) == ℂ[3]'
end

@testset TimedTestSet "ElementarySpace: GeneralSpace" begin
    d = 2
    V = GeneralSpace{ℂ}(d)
    @test eval(Meta.parse(sprint(show, V))) == V
    @test eval(Meta.parse(sprint(show, dual(V)))) == dual(V)
    @test eval(Meta.parse(sprint(show, conj(V)))) == conj(V)
    @test eval(Meta.parse(sprint(show, typeof(V)))) == typeof(V)
    @test !isdual(V)
    @test isdual(V')
    @test !isdual(conj(V))
    @test isdual(conj(V'))
    @test !TensorKit.isconj(V)
    @test !TensorKit.isconj(V')
    @test TensorKit.isconj(conj(V))
    @test TensorKit.isconj(conj(V'))
    @test isa(V, VectorSpace)
    @test isa(V, ElementarySpace)
    @test !isa(V, InnerProductSpace)
    @test !isa(V, EuclideanSpace)
    @test @inferred(hash(V)) == hash(deepcopy(V)) != hash(V')
    @test @inferred(dual(V)) != @inferred(conj(V)) != V
    @test @inferred(field(V)) == ℂ
    @test @inferred(sectortype(V)) == Trivial
    @test @inferred(TensorKit.hassector(V, Trivial()))
    @test @inferred(dim(V)) == d == @inferred(dim(V, Trivial()))
    @test @inferred(TensorKit.axes(V)) == Base.OneTo(d)
end

@testset TimedTestSet "ElementarySpace: RepresentationSpace{$G}" for G in (ℤ₂, ℤ₃, ℤ₄, U₁, CU₁, SU₂, FibonacciAnyon, ℤ₃ × ℤ₄, U₁ × SU₂, SU₂ × SU₂, ℤ₂ × FibonacciAnyon × FibonacciAnyon)
    if Base.IteratorSize(values(G)) === Base.IsInfinite()
        set = unique(vcat(one(G), [randsector(G) for k = 1:10]))
        gen = (c=>2 for c in set)
    else
        gen = (values(G)[k]=>(k+1) for k in 1:length(values(G)))
    end
    V = RepresentationSpace(gen)
    @test eval(Meta.parse(sprint(show, V))) == V
    @test eval(Meta.parse(sprint(show, V'))) == V'
    @test V' == RepresentationSpace(gen; dual = true)
    @test V == @inferred RepresentationSpace(gen...)
    @test V' == @inferred RepresentationSpace(gen...; dual = true)
    @test V == @inferred RepresentationSpace(tuple(gen...))
    @test V' == @inferred RepresentationSpace(tuple(gen...); dual = true)
    @test V == @inferred RepresentationSpace(Dict(gen))
    @test V' == @inferred RepresentationSpace(Dict(gen); dual = true)
    @test V == @inferred RepresentationSpace{G}(gen)
    @test V' == @inferred RepresentationSpace{G}(gen; dual = true)
    @test V == @inferred RepresentationSpace{G}(gen...)
    @test V' == @inferred RepresentationSpace{G}(gen...; dual = true)
    @test V == @inferred RepresentationSpace{G}(Dict(gen))
    @test V' == @inferred RepresentationSpace{G}(Dict(gen); dual = true)
    @test V == @inferred typeof(V)(c=>dim(V,c) for c in sectors(V))
    if G isa ZNIrrep
        @test V == @inferred typeof(V)(V.dims)
        @test V' == @inferred typeof(V)(V.dims; dual = true)
    end
    @test @inferred(hash(V)) == hash(deepcopy(V)) != hash(V')
    @test V == RepresentationSpace(reverse(collect(gen))...)
    @test eval(Meta.parse(sprint(show, V))) == V
    @test eval(Meta.parse(sprint(show, typeof(V)))) == typeof(V)
    # space with no sectors
    @test dim(@inferred(typeof(V)())) == 0
    # space with a single sector
    W = @inferred ℂ[one(G)=>1]
    @test W == RepresentationSpace(one(G)=>1, randsector(G) => 0)
    @test @inferred(oneunit(V)) == W == oneunit(typeof(V))
    # randsector never returns trivial sector, so this cannot error
    @test_throws ArgumentError RepresentationSpace(one(G)=>1, randsector(G) => 0, one(G)=>3)
    @test eval(Meta.parse(sprint(show, W))) == W
    @test isa(V, VectorSpace)
    @test isa(V, ElementarySpace)
    @test isa(V, InnerProductSpace)
    @test isa(V, EuclideanSpace)
    @test isa(V, RepresentationSpace)
    @test isa(V, RepresentationSpace{G})
    @test isa(V, Base.IteratorSize(values(G)) == Base.IsInfinite() ?
                    TensorKit.GenericRepresentationSpace{G} :
                    TensorKit.FiniteRepresentationSpace{G})
    @test @inferred(dual(V)) == @inferred(conj(V)) == @inferred(adjoint(V)) != V
    @test @inferred(field(V)) == ℂ
    @test @inferred(sectortype(V)) == G
    slist = @inferred sectors(V)
    @test @inferred(TensorKit.hassector(V, first(slist)))
    @test @inferred(dim(V)) == sum((@inferred(dim(s)*dim(V, s))) for s in slist)
    if hasfusiontensor(G)
        @test @inferred(TensorKit.axes(V)) == Base.OneTo(dim(V))
    end
    @test @inferred(⊕(V,V)) == RepresentationSpace{G}(c=>2dim(V,c) for c in sectors(V))
    @test @inferred(⊕(V,V,V,V)) == RepresentationSpace{G}(c=>4dim(V,c) for c in sectors(V))
    @test @inferred(⊕(V,oneunit(V))) ==
            RepresentationSpace{G}(c=>isone(c)+dim(V,c) for c in sectors(V))
    @test @inferred(fuse(V,oneunit(V))) == V
    d = Dict{G,Int}()
    for a in sectors(V), b in sectors(V)
        for c in a ⊗ b
            d[c] = get(d, c, 0) + dim(V, a)*dim(V, b)*Nsymbol(a,b,c)
        end
    end
    @test @inferred(fuse(V,V)) == RepresentationSpace(d)
    @test @inferred(flip(V)) ==
            RepresentationSpace{G}(conj(c)=>dim(V,c) for c in sectors(V))'
    @test @inferred(⊕(V, V)) == @inferred max(V, ⊕(V, V))
    @test V == @inferred min(V, ⊕(V, V))
    @test min(V, RepresentationSpace(one(G)=>3)) == RepresentationSpace(one(G)=>2)
    @test_throws SpaceMismatch (⊕(V, V'))
end

@testset TimedTestSet "ProductSpace{ℂ}" begin
    V1, V2, V3, V4 = ℂ[1], ℂ[2], ℂ[3], ℂ[4]
    P = @inferred ProductSpace(V1, V2, V3, V4)
    @test eval(Meta.parse(sprint(show, P))) == P
    @test eval(Meta.parse(sprint(show, typeof(P)))) == typeof(P)
    @test isa(P, VectorSpace)
    @test isa(P, CompositeSpace)
    @test spacetype(P) == ComplexSpace
    @test sectortype(P) == Trivial
    @test @inferred(hash(P)) == hash(deepcopy(P)) != hash(P')
    @test P == deepcopy(P)
    @test P == typeof(P)(P...)
    @test @inferred(dual(P)) == P'
    @test @inferred(field(P)) == ℂ
    @test @inferred(*(V1, V2, V3, V4)) == P
    @test @inferred(⊗(V1, V2, V3, V4)) == P
    @test @inferred(⊗(V1, V2⊗V3⊗V4)) == P
    @test @inferred(⊗(V1⊗V2, V3⊗V4)) == P
    @test @inferred(⊗(V1, V2, V3⊗V4)) == P
    @test @inferred(⊗(V1, V2⊗V3, V4)) == P
    @test @inferred(⊗(P)) == P
    @test @inferred(⊗(V1)) == ProductSpace(V1)
    @test eval(Meta.parse(sprint(show, ⊗(V1)))) == ⊗(V1)
    @test @inferred(one(V1)) == @inferred(one(typeof(V1))) ==
                @inferred(one(P)) == @inferred(one(typeof(P))) ==
                ProductSpace{ComplexSpace}(())
    @test eval(Meta.parse(sprint(show, one(P)))) == one(P)
    @test @inferred(⊗(one(P), P)) == P
    @test @inferred(⊗(P, one(P))) == P
    @test @inferred(⊗(one(P), one(P))) == one(P)
    @test @inferred(adjoint(P)) == dual(P) == V4' ⊗ V3' ⊗ V2' ⊗ V1'
    @test @inferred(dims(P)) == map(dim, (V1, V2, V3, V4))
    @test @inferred(dim(P)) == prod(dim, (V1, V2, V3, V4))
    @test @inferred(dim(P, 2)) == dim(V2)
    @test @inferred(sectors(P)) ==
            (mapreduce(sectors, (a, b)->tuple(a..., b...), (V1, V2, V3, V4)),)
    cube(x) = x^3
    @test @inferred(cube(V1)) == V1 ⊗ V1 ⊗ V1
    N = 3
    @test V1^N == V1 ⊗ V1 ⊗ V1
    @test P^2 == P ⊗ P
    @test @inferred(dims(P, first(sectors(P)))) == dims(P)
    @test @inferred(blocksectors(P)) == (Trivial(),)
    @test @inferred(blockdim(P, first(blocksectors(P)))) == dim(P)
    @test Base.IteratorEltype(P) == Base.IteratorEltype(typeof(P)) ==
                                    Base.IteratorEltype(P.spaces)
    @test Base.IteratorSize(P) == Base.IteratorSize(typeof(P)) ==
                                    Base.IteratorSize(P.spaces)
    @test Base.eltype(P) == Base.eltype(typeof(P)) == typeof(V1)
    @test eltype(collect(P)) == typeof(V1)
    @test collect(P) == [V1, V2, V3, V4]
end

@testset TimedTestSet "ProductSpace{SU₂Space}" begin
    V1, V2, V3, V4, V5 = SU₂Space(0=>3, 1//2=>1), SU₂Space(0=>2, 1=>1),
                            SU₂Space(1//2=>1, 1=>1)', SU₂Space(0=>2, 1//2=>2),
                            SU₂Space(0=>1, 1//2=>1, 3//2=>1)'
    W = TensorKit.HomSpace(V1 ⊗ V2, V3 ⊗ V4 ⊗ V5)
    @test W == (V3 ⊗ V4 ⊗ V5 → V1 ⊗ V2)
    @test W == (V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5)
    @test W' == (V1 ⊗ V2 → V3 ⊗ V4 ⊗ V5)
    @test eval(Meta.parse(sprint(show, W))) == W
    @test eval(Meta.parse(sprint(show, typeof(W)))) == typeof(W)
    @test spacetype(W) == SU₂Space
    @test sectortype(W) == SU₂
    @test W[1] == V1
    @test W[2] == V2
    @test W[3] == V3'
    @test W[4] == V4'
    @test W[5] == V5'
    @test @inferred(hash(W)) == hash(deepcopy(W)) != hash(W')
    @test W == deepcopy(W)
end

@testset TimedTestSet "HomSpace" begin
    V1, V2, V3 = SU₂Space(0=>3, 1//2=>1), SU₂Space(0=>2, 1=>1), SU₂Space(1//2=>1, 1=>1)'
    P = @inferred ProductSpace(V1, V2, V3)
    @test eval(Meta.parse(sprint(show, P))) == P
    @test eval(Meta.parse(sprint(show, typeof(P)))) == typeof(P)
    @test isa(P, VectorSpace)
    @test isa(P, CompositeSpace)
    @test spacetype(P) == SU₂Space
    @test sectortype(P) == SU₂
    @test @inferred(hash(P)) == hash(deepcopy(P)) != hash(P')
    @test @inferred(dual(P)) == P'
    @test @inferred(field(P)) == ℂ
    @test @inferred(*(V1, V2, V3)) == P
    @test @inferred(⊗(V1, V2, V3)) == P
    @test @inferred(adjoint(P)) == dual(P) == V3' ⊗ V2' ⊗ V1'
    @test @inferred(⊗(V1)) == ProductSpace(V1)
    @test @inferred(one(V1)) == @inferred(one(typeof(V1))) ==
                @inferred(one(P)) == @inferred(one(typeof(P))) ==
                ProductSpace{ComplexSpace}(())
    @test @inferred(dims(P)) == map(dim, (V1, V2, V3))
    @test @inferred(dim(P)) == prod(dim, (V1, V2, V3))
    for s in @inferred(sectors(P))
        @test hassector(P, s)
        @test @inferred(dims(P, s)) == dim.((V1, V2, V3), s)
    end
    @test sum(dim(c)*blockdim(P, c) for c in @inferred(blocksectors(P))) == dim(P)
end
tf = time()
printstyled("Finished vector space tests in ",
            string(round(tf-ti; sigdigits=3)),
            " seconds."; bold = true, color = Base.info_color())
println()


text:
#abstract type AbstractTarget end
mutable struct Target1D
    size::Float64
    position::Float64
end

"""
    TargetRect(a,b,position)
defines a rectangular target by specifying the lenght a, the height b and the position (range)
"""
mutable struct TargetRect <:AbstractTarget
    a::Float64
    b::Float64
    position::Union{Array{Float64,1},Nothing}
    ρ::Union{Float64,Nothing}
end

"""
    TargetCirc(ρ,position)
defines a circular target by specifying the radius ρ and the position (range)
"""
mutable struct TargetCirc <:AbstractTarget
    ρ::Float64 #radius
    position::Float64
end
"""
    TargetSphe(ρ,position)
defines a spherical target by specifying the radius ρ and the position (range)
"""
mutable struct TargetSphe <:AbstractTarget
    ρ::Float64 #radius
    position::Float64
end

mutable struct Fuze
    offset::Float64
    height::Float64
end


text:
"""
Test Sundials:
This case study a three bus system with 1 machine (One d- One q-: 4th order model), a VSM of 19 states and an infinite source. The connection between buses 2 and 3 is modeled using dynamic lines.
The perturbation trips two of the three circuits of line between buses 1 and 2, triplicating its impedance.
The objective of the test is to assess if the model can exploit Sundials' multiple linear solvers
"""

##################################################
############### LOAD DATA ########################
##################################################

include(joinpath(dirname(@__FILE__), "data_tests/test11.jl"))

##################################################
############### SOLVE PROBLEM ####################
##################################################
solver_list = [
    :Dense,
    #:Band Requires jac_upper, jac_lower
    :LapackDense,
    #:LapackBand Requires jac_upper, jac_lower,
    :GMRES,
    :BCG,
    :PCG,
    :TFQMR,
    :KLU,
]

#time span
tspan = (0.0, 40.0)

#Define Fault: Change of YBus
Ybus_change = NetworkSwitch(
    1.0, #change at t = 1.0
    Ybus_fault,
) #New YBus

function test_sundials(solver)
    path = (joinpath(pwd(), "test-sundials"))
    !isdir(path) && mkdir(path)
    try
        #Define Simulation Problem
        sim = Simulation!(
            ImplicitModel,
            path,
            threebus_sys, #system
            tspan, #time span
            Ybus_change, #Type of Fault
        )

        #Obtain small signal results for initial conditions
        small_sig = small_signal_analysis(sim)

        #Solve problem
        @info "$(solver)" @time execute!(sim, IDA(linear_solver = solver);)

        #Obtain data for voltages
        series = get_voltage_magnitude_series(sim, 102)
        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test10_x0_init
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end
        @test (diff[1] < 1e-3)
        @test sim.solution.retcode == :Success
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
    return
end

@testset "Sundials Tests" begin
    for name in solver_list
        @testset "$(name)" begin
            test_sundials(name)
        end
    end
end


text:
import CLBLAS
using OpenCL: cl

const clblas = CLBLAS
clblas.setup()

alpha = cl.cl_float(5)
beta = cl.cl_float(2)
A = float32([1.0 2.0 3.0; 4.0 5.0 6.0])  #rand(Float32, 2, 3)
B = float32([1.0 2.0; 3.0 4.0; 5.0 6.0]) #rand(Float32, 3, 2)
C = float32([1.0 2.0; 3.0 4.0])#rand(Float32, 2, 2)

data = (5.0*A*B) + (2.0*C)

future = clblas.clblasSgemm(UInt32(0), UInt32(0), alpha, A, B, beta, C)

result = fetch(future)
println(result)
println("------------")
println(data)
@assert(result == data)
info("success!")

clblas.teardown()


text:
using Ipopt
using ForwardDiff
using CuArrays, CUDAnative
using TimerOutputs
using SparseDiffTools
using SparseArrays
using acopf
timeroutput = TimerOutput()

function test(Pg0, Qg0, Vm0, Va0, timeroutput, opfdata, arraytype; max_iter=100)

  if arraytype == CuArray
    T = CuVector
  elseif arraytype == Array
    T = Vector
  else
    error("Unkown array type $arraytype.")
  end

  nPg = size(Pg0,1) ; nQg = size(Qg0,1) ; nVm = size(Vm0,1) ; nVa = size(Va0,1)
  n = nPg + nQg + nVm + nVa
  nbus = length(opfdata.buses)
  nline = length(opfdata.lines)
  ngen = length(opfdata.generators)
  m = 2 * nbus + 2 * nline

  function get_jac_sparsity()
    println("Computing coloring...")
    x = Vector{Float64}(ones(Float64, n))
    Pg = Vector{Float64}(zeros(Float64, nPg))
    Qg = Vector{Float64}(zeros(Float64, nQg))
    Va = Vector{Float64}(zeros(Float64, nVa))
    Vm = Vector{Float64}(zeros(Float64, nVm))
    rbalconst = Vector{Float64}(zeros(Float64, nbus))
    ibalconst = Vector{Float64}(zeros(Float64, nbus))
    limitsto = Vector{Float64}(zeros(Float64, nline))
    limitsfrom = Vector{Float64}(zeros(Float64, nline))
    function constraints(x)
      Pg = x[1:nPg]
      Qg = x[nPg+1:nPg+nQg]
      Va = x[nPg+nQg+1:nPg+nQg+nVa]
      Vm = x[nPg+nQg+nVa+1:end]
      acopf.update_arrays!(arrays, Pg, Qg, Va, Vm, timeroutput)
      acopf.constraints(rbalconst, ibalconst, limitsto, limitsfrom, arrays, timeroutput)
      y[1:nbus] = rbalconst[:]
      y[nbus+1:2*nbus] = ibalconst[:]
      y[2*nbus+1:2*nbus+nline] = limitsto[:]
      y[2*nbus+nline+1:end] = limitsfrom[:]
      return y
    end
    cux = Vector{Float64}(x)
    cfg = ForwardDiff.JacobianConfig(constraints, cux)
    rbalconst   = Vector{eltype(cfg)}(undef, nbus)
    ibalconst   = Vector{eltype(cfg)}(undef, nbus)
    limitsto    = Vector{eltype(cfg)}(undef, nline)
    limitsfrom  = Vector{eltype(cfg)}(undef, nline)
    y = Vector{eltype(cfg)}(undef, 2*nbus+2*nline)
    diffcux = Vector{eltype(cfg)}(x)
    Pg = diffcux[1:nPg]
    Qg = diffcux[nPg+1:nPg+nQg]
    Va = diffcux[nPg+nQg+1:nPg+nQg+nVa]
    Vm = diffcux[nPg+nQg+nVa+1:end]
    arrays = acopf.create_arrays(Pg, Qg, Va, Vm, opfdata, timeroutput, Array)
    fjac = cux -> ForwardDiff.jacobian(constraints, cux, cfg)
    jac = sparse(fjac(cux))
    colors = unique(matrix_colors(jac))
    println("Done.")
    return colors
  end
  cuPg = T{Float64}(zeros(Float64, nPg))
  cuQg = T{Float64}(zeros(Float64, nQg))
  cuVa = T{Float64}(zeros(Float64, nVa))
  cuVm = T{Float64}(zeros(Float64, nVm))
  curbalconst = T{Float64}(zeros(Float64, nbus))
  cuibalconst = T{Float64}(zeros(Float64, nbus))
  culimitsto = T{Float64}(zeros(Float64, nline))
  culimitsfrom = T{Float64}(zeros(Float64, nline))
  # println("Number of Jacobian colors: ", size(get_jac_sparsity(),1))
  # return
  io = open("log.out", "w+")
  function myprint(name,var)
      return
      write(io, name * ": ")
      for v in var
        if v != 0.0
          write(io, string(v) * " ")
        end
      end
      write(io,"\n")
  end
  function eval_f(x::Vector{Float64})
    @timeit timeroutput "eval_f" begin
    cuPg[:] = x[1:nPg]
    cuQg[:] = x[nPg+1:nPg+nQg]
    cuVa[:] = x[nPg+nQg+1:nPg+nQg+nVa]
    cuVm[:] = x[nPg+nQg+nVa+1:end]
    arrays = acopf.create_arrays(cuPg, cuQg, cuVa, cuVm, opfdata, timeroutput, arraytype)
    obj = acopf.objective(arrays, timeroutput)
    myprint("fx", x)
    myprint("fy", obj)
    end
    return obj
  end

  function eval_g(x::Vector{Float64}, g::Vector{Float64})
    @timeit timeroutput "eval_g" begin
    cuPg[:] = x[1:nPg]
    cuQg[:] = x[nPg+1:nPg+nQg]
    cuVa[:] = x[nPg+nQg+1:nPg+nQg+nVa]
    cuVm[:] = x[nPg+nQg+nVa+1:end]
    arrays = acopf.create_arrays(cuPg, cuQg, cuVa, cuVm, opfdata, timeroutput, arraytype)
    acopf.constraints(curbalconst, cuibalconst, culimitsto, culimitsfrom, arrays, timeroutput)
    g[1:nbus] = curbalconst[:]
    g[nbus+1:2*nbus] = cuibalconst[:]
    g[2*nbus+1:2*nbus+nline] = culimitsto[:]
    g[2*nbus+nline+1:end] = culimitsfrom[:]
    myprint("gx", x)
    myprint("g", g)
    end
  end

  function eval_grad_f(x::Vector{Float64}, grad_f::Vector{Float64})
    @timeit timeroutput "eval_grad_f" begin
    function objective(x)
      Pg = x[1:nPg]
      Qg = x[nPg+1:nPg+nQg]
      Va = x[nPg+nQg+1:nPg+nQg+nVa]
      Vm = x[nPg+nQg+nVa+1:end]
      arrays = acopf.create_arrays(Pg, Qg, Va, Vm, opfdata, timeroutput, arraytype)
      return acopf.objective(arrays, timeroutput)
    end
    cux = T{Float64}(x)
    g = cux -> ForwardDiff.gradient(objective, cux)
    grad_f[:] = g(cux)
    myprint("grad_x", x)
    myprint("grad_f", grad_f)
    end
    return
  end

  function eval_jac_g(x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, values::Vector{Float64})
    @timeit timeroutput "eval_jac_g" begin
    if mode == :Structure
      idx = 1
      for c in 1:m #number of constraints
        for i in 1:n # number of variables
          rows[idx] = c ; cols[idx] = i
          idx += 1
        end
      end
    else
      function constraints(x)
        Pg = x[1:nPg]
        Qg = x[nPg+1:nPg+nQg]
        Va = x[nPg+nQg+1:nPg+nQg+nVa]
        Vm = x[nPg+nQg+nVa+1:end]
        acopf.update_arrays!(arrays, Pg, Qg, Va, Vm, timeroutput)
        acopf.constraints(rbalconst, ibalconst, limitsto, limitsfrom, arrays, timeroutput)
        y[1:nbus] = rbalconst[:]
        y[nbus+1:2*nbus] = ibalconst[:]
        y[2*nbus+1:2*nbus+nline] = limitsto[:]
        y[2*nbus+nline+1:end] = limitsfrom[:]
        return y
      end
      cux = T{Float64}(x)
      cfg = ForwardDiff.JacobianConfig(constraints, cux)
      rbalconst   = T{eltype(cfg)}(undef, nbus)
      ibalconst   = T{eltype(cfg)}(undef, nbus)
      limitsto    = T{eltype(cfg)}(undef, nline)
      limitsfrom  = T{eltype(cfg)}(undef, nline)
      y = T{eltype(cfg)}(undef, 2*nbus+2*nline)
      diffcux = T{eltype(cfg)}(x)
      Pg = diffcux[1:nPg]
      Qg = diffcux[nPg+1:nPg+nQg]
      Va = diffcux[nPg+nQg+1:nPg+nQg+nVa]
      Vm = diffcux[nPg+nQg+nVa+1:end]
      arrays = acopf.create_arrays(Pg, Qg, Va, Vm, opfdata, timeroutput, arraytype)
      fjac = cux -> ForwardDiff.jacobian(constraints, cux, cfg)
      jac = fjac(cux)
      k = 1
      for i in 1:m
        for j in 1:n
          values[k] = jac[i,j]
          k += 1
        end
      end
      myprint("jac_x", x)
      myprint("jac_g", values)
    end
    end
    return
  end

  function eval_h(x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, obj_factor::Float64, lambda::Vector{Float64}, values::Vector{Float64})
    @timeit timeroutput "eval_h" begin
    if mode == :Structure
      idx = 1
      for row = 1:n
        for col = 1:row
          rows[idx] = row
          cols[idx] = col
          idx += 1
        end
      end
    else
      function objective(x)
        Pg = x[1:nPg]
        Qg = x[nPg+1:nPg+nQg]
        Va = x[nPg+nQg+1:nPg+nQg+nVa]
        Vm = x[nPg+nQg+nVa+1:end]
        arrays = acopf.create_arrays(Pg, Qg, Va, Vm, opfdata, timeroutput, arraytype)
        return acopf.objective(arrays, timeroutput)
      end

      @timeit timeroutput "moving to GPU" begin
      cux = T{Float64}(x)
      end
      h = cux -> ForwardDiff.hessian(objective, cux)
      objhess = h(cux)
      @timeit timeroutput "moving from GPU objective" begin
      k = 1
      for i in 1:n
        for j in 1:i
          values[k] = obj_factor * objhess[i,j]
          k += 1
        end
      end
      end
      function constraints(x)
        Pg = x[1:nPg]
        Qg = x[nPg+1:nPg+nQg]
        Va = x[nPg+nQg+1:nPg+nQg+nVa]
        Vm = x[nPg+nQg+nVa+1:end]
        acopf.update_arrays!(arrays, Pg, Qg, Va, Vm, timeroutput)
        acopf.constraints(rbalconst, ibalconst, limitsto, limitsfrom, arrays, timeroutput)
        y[1:nbus] = rbalconst[:]
        y[nbus+1:2*nbus] = ibalconst[:]
        y[2*nbus+1:2*nbus+nline] = limitsto[:]
        y[2*nbus+nline+1:end] = limitsfrom[:]
        return y
      end
      @timeit timeroutput "remaining" begin
      cfg1 = ForwardDiff.JacobianConfig(nothing, x, ForwardDiff.Chunk(x))
      cfg = ForwardDiff.JacobianConfig(nothing, cfg1.duals, ForwardDiff.Chunk(x))
      rbalconst   = T{eltype(cfg)}(undef, nbus)
      ibalconst   = T{eltype(cfg)}(undef, nbus)
      limitsto    = T{eltype(cfg)}(undef, nline)
      limitsfrom  = T{eltype(cfg)}(undef, nline)
      y = T{eltype(cfg)}(undef, 2*nbus+2*nline)
      diffcux = T{eltype(cfg)}(x)
      Pg = diffcux[1:nPg]
      Qg = diffcux[nPg+1:nPg+nQg]
      Va = diffcux[nPg+nQg+1:nPg+nQg+nVa]
      Vm = diffcux[nPg+nQg+nVa+1:end]
      end
      @timeit timeroutput "compute hessian" begin
      arrays = acopf.create_arrays(Pg, Qg, Va, Vm, opfdata, timeroutput, arraytype)
      fhess = x -> ForwardDiff.jacobian(x -> ForwardDiff.jacobian(constraints, x), x)
      hess = fhess(cux)
      end
      @timeit timeroutput "reshape" begin
      hess = reshape(hess, m, n, n)
      end

      @timeit timeroutput "moving from GPU constraint" begin
      for l in 1:m
        k = 1
        for i in 1:n
          for j in 1:i
            values[k] += lambda[l] * hess[l,i,j]
            k += 1
          end
        end
      end
      end
      myprint("h_x", x)
      myprint("h", values)
    end
    end
    return
  end
  x_L = Vector{Float64}(undef,n)
  x_U = Vector{Float64}(undef,n)
  for (i,g) in enumerate(opfdata.generators)
    x_L[i] = g.Pmin
    x_U[i] = g.Pmax
    x_L[nPg+i] = g.Qmin
    x_U[nPg+i] = g.Qmax
  end
  for (i,b) in enumerate(opfdata.buses)
    x_L[nPg + nQg + i] = -Inf
    x_U[nPg + nQg + i] = Inf
    x_L[nPg + nQg + nVa + i] = b.Vmin
    x_U[nPg + nQg + nVa + i] = b.Vmax
  end
  x_L[nPg+nQg+1] = 0.0
  x_U[nPg+nQg+1] = 0.0

  g_L = Vector{Float64}(undef, m)
  g_U = Vector{Float64}(undef, m)
  for i in 1:2*nbus g_L[i] = 0.0 end
  for i in 1:2*nbus g_U[i] = 0.0 end
  for i in 2*nbus+1:m g_L[i] = -Inf end
  for i in 2*nbus+1:m g_U[i] = 0.0 end

  myprint("x_L", x_L)
  myprint("x_U", x_U)

  myprint("g_L", g_L)
  myprint("g_U", g_U)

  idx = 1
  for row = 1:n
    for col = 1:row
      idx += 1
    end
  end
  prob = createProblem(n, x_L, x_U, m, g_L, g_U, m*n, idx-1,
                      eval_f, eval_g,eval_grad_f, eval_jac_g, eval_h)
  prob.x = 	[0.0 for i in 1:n]

  # Callback
  function intermediate(alg_mod::Int, iter_count::Int,
  obj_value::Float64, inf_pr::Float64, inf_du::Float64, mu::Float64,
  d_norm::Float64, regularization_size::Float64, alpha_du::Float64, alpha_pr::Float64,
  ls_trials::Int)
  return iter_count < max_iter  # Interrupts after one iteration.
  end

  setIntermediateCallback(prob, intermediate)

  solvestat = solveProblem(prob)
  close(io)
  return prob.obj_val

end


text:
export domaintype,
    codomaintype,
    coefficienttype,
    prectype

"""
A `Dictionary{S,T}` is an ordered family of functions, in which each function
maps a variable of type `S` to a variable of type `T`. The dictionary can be
thought of as an array, where the elements are functions that map `S` to `T`.

A `Dictionary{S,T}` has domain type `S` and codomain type `T`. The domain type
corresponds to the type of a domain in the `DomainSets.jl` package, and it is the
type of the expected argument to the elements of the dictionary. The
codomain type is the type of the output.

Each dictionary is ordered via its index set: the ordering is determined by the
iterator of the index set. A dictionary `d` can be indexed in several ways:
- the linear index is a positive natural number between `1` and `length(d)`
- the native index is an index that more closely corresponds to the conventional
  mathematical notation of the dictionary, if that differs from linear indexing.
  For example, polynomials may have degree ranging from `0` to `length(d)-1`.
  Fourier series may have negative frequencies.

Some dictionaries define other types of indices, such as multilinear and product
indices. See `CompositeDictionary` and `ProductDictionary`. There is always a
bijection between the different types of index sets supported by a dictionary.

A dictionary is the basic building block of the package. They can support several
operations, such as interpolation, evaluation, differentiation and others. This
functionality is made available using a specific interface. See the documentation
of the features for a description of that interface and the syntax.
"""
abstract type Dictionary{S,T}
end

# Useful abstraction for special cases
const Dictionary1d{S <: Number,T} = Dictionary{S,T}
# Warning: these are shaky definitions of multivariate functions
const Dictionary2d{S <: Number,T} = Union{Dictionary{NTuple{2,S},T},Dictionary{SVector{2,S},T}}
const Dictionary3d{S <: Number,T} = Union{Dictionary{NTuple{3,S},T},Dictionary{SVector{3,S},T}}
const Dictionary4d{S <: Number,T} = Union{Dictionary{NTuple{4,S},T},Dictionary{SVector{4,S},T}}

"The type of the elements of the domain of the dictionary."
domaintype(::Type{<:Dictionary{S,T}}) where {S,T} = S
domaintype(dict::Dictionary) = domaintype(typeof(dict))

"The type of the elements of the codomain of the dictionary."
codomaintype(::Type{<:Dictionary{S,T}}) where {S,T} = T
codomaintype(dict::Dictionary) = codomaintype(typeof(dict))

"The type of the expansion coefficients in a dictionary."
# By default we set it equal to the codomaintype
coefficienttype(D::Type{<:Dictionary{S,T}}) where {S,T} = codomaintype(D)
coefficienttype(dict::Dictionary) = coefficienttype(typeof(dict))

prectype(D::Type{<:Dictionary}) = prectype(domaintype(D), codomaintype(D))
numtype(D::Type{<:Dictionary}) = numtype(domaintype(D), codomaintype(D))

# The dimension of a function set is the dimension of its domain type
dimension(dict::Dictionary) = dimension(domaintype(dict))

dimension(dict::Dictionary, i) = dimension(component(dict, i))

"Are the functions in the dictionary real-valued?"
isreal(d::Dictionary) = isreal(codomaintype(d))

"Is the dictionary a (truncation of a) basis?"
isbasis(d::Dictionary) = false

"Is the dictionary orthogonal (with respect to the given measure)?"
isorthogonal(d::Dictionary) = hasmeasure(d) && isorthogonal(d, measure(d))
isorthogonal(d::Dictionary, μ::Measure) = isorthonormal(d, μ)

"Is the dictionary orthonormal (with respect to the given measure)?"
isorthonormal(d::Dictionary) = hasmeasure(d) && isorthonormal(d, measure(d))
isorthonormal(d::Dictionary, μ::Measure) = false

isorthogonal(d::Dictionary, grid::AbstractGrid) = isorthogonal(d, discretemeasure(grid))
isorthonormal(d::Dictionary, grid::AbstractGrid) = isorthonormal(d, discretemeasure(grid))

"Is the dictionary biorthogonal (with respect to the given measure)?"
isbiorthogonal(d::Dictionary) = hasmeasure(d) && isbiorthogonal(d, measure(d))
isbiorthogonal(d::Dictionary, μ::Measure) = isorthogonal(d, μ)

size(d::Dictionary, j) = size(d)[j]

"""
The length and size of a dictionary may not be enough to uniquely determine the
size of the dictionary, if it has additional internal structure. The output of
dimensions is such that `resize(dict, dimensions(dict)) == dict`.
"""
dimensions(d::Dictionary1d) = length(d)
dimensions(d::Dictionary) = size(d)

length(d::Dictionary) = prod(size(d))

firstindex(d::Dictionary) = first(eachindex(d))
lastindex(d::Dictionary) = last(eachindex(d))

#############################
# Domain and codomain type
#############################

similar(s::Dictionary, ::Type{T}) where {T} = similar(s, T, size(s))

similar(s::Dictionary, size::Int...) = similar(s, domaintype(s), size...)
similar(s::Dictionary, dims::Base.Dims) = similar(s, domaintype(s), dims...)

similar(s::Dictionary, ::Type{T}, dims::Base.Dims) where {T} = similar(s, T, dims...)

# This is a default routine that ony appies when nothing changes.
function similar(dict::Dictionary{S,T}, ::Type{S}, n::Int) where {S,T}
    @assert n == length(dict)
    dict
end

resize(s::Dictionary, dims...) = similar(s, domaintype(s), dims...)

widen(d::Dictionary) = similar(d, widen(domaintype(d)))

"Return a set of zero coefficients in the native format of the set."
zeros(s::Dictionary) = zeros(coefficienttype(s), s)
ones(s::Dictionary) = ones(coefficienttype(s), s)

# By default we assume that the native format corresponds to an array of the
# same size as the set. This is not true, e.g., for multidicts.
zeros(::Type{T}, s::Dictionary) where {T} = zeros(T, size(s))
ones(::Type{T}, s::Dictionary) where {T} = ones(T, size(s))

"What is the type of the coefficient vector of the dictionary?"
containertype(d::Dictionary) = typeof(zeros(d))

"Transform `a` to a coefficient vector for the given dictionary."
tocoefficientformat(a, d::Dictionary) = reshape(a, size(d))

function rand(dict::Dictionary)
    c = zeros(dict)
    T = coefficienttype(dict)
    for i in eachindex(c)
        c[i] = rand(T)
    end
    c
end

###########
# Indexing
###########

IndexStyle(d::Dictionary) = IndexLinear()

"""
Dictionaries are ordered lists. Their ordering is defined by the way their
index sets are ordered.

The `ordering` of a dictionary returns a list-like object that can be indexed
with integers between `1` and `length(dict)`. This operation returns the
corresponding native index. This defines the ordering of the native index set
of the dictionary.
"""
ordering(dict::Dictionary) = Base.OneTo(length(dict))

eachindex(d::Dictionary) = eachindex(IndexStyle(d), d)
eachindex(::IndexLinear, d::Dictionary) = axes1(d)
eachindex(::IndexCartesian, d::Dictionary) = CartesianIndices(axes(d))

axes1(d::Dictionary) = Base.OneTo(length(d))
axes(d::Dictionary) = map(Base.OneTo, size(d))

"Compute the native index corresponding to the given index."
native_index(dict::Dictionary, idx) = _native_index(dict, idx)
# We redirect to a fallback _native_index in case the concrete dictionary
# did not implement native_index.  We explicitly convert a linear index using the ordering.
# Anything else we throw an error because the index looks invalid
_native_index(dict::Dictionary, idx::Int) = ordering(dict)[idx]
_native_index(dict::Dictionary, idx::NativeIndex) = idx
_native_index(dict::Dictionary, idx::AbstractShiftedIndex) = idx
_native_index(dict::Dictionary, idx) = throw(ArgumentError("invalid index: $idx"))

"Compute the linear index corresponding to the given index."
linear_index(dict::Dictionary, idx) = _linear_index(dict, idx)
# We can accept an integer unchanged, anything else we pass to the ordering
_linear_index(dict::Dictionary, idx::Int) = idx
_linear_index(dict::Dictionary, idxn) = ordering(dict)[idxn]

##################################################
## Conversion between coefficient representations
##################################################

"""
Convert the set of coefficients in the native format of the dictionary
to a linear list in a vector.
"""
function linearize_coefficients(dict::Dictionary, coef_native)
    coef_linear = zeros(eltype(coef_native), length(dict))
    linearize_coefficients!(dict, coef_linear, coef_native)
end

linearize_coefficients!(dict::Dictionary, coef_linear::Vector, coef_native) =
    copyto!(coef_linear, coef_native)
# Note that copyto! is defined in util/common.jl

"""
Convert a linear set of coefficients back to the native representation of the dictionary.
"""
function delinearize_coefficients(dict::Dictionary, coef_linear::Vector)
    coef_native = zeros(eltype(coef_linear), dict)
    delinearize_coefficients!(dict, coef_native, coef_linear)
end

delinearize_coefficients!(dict::Dictionary, coef_native, coef_linear::Vector) =
    copyto!(coef_native, coef_linear)

"Promote the given coefficients to the native representation of the dictionary."
native_coefficients(dict::Dictionary, coef) = _native_coefficients(dict, coef)
# TODO: we create an unnecessary copy here if the native type is a vector
_native_coefficients(dict::Dictionary, coef::Vector) = delinearize_coefficients(dict, coef)
_native_coefficients(dict::Dictionary, coef) = coef

# Sets have a native size and a linear size. However, there is not necessarily a
# bijection between the two. You can always convert a native size to a linear size,
# but the other direction can be done in general only approximately.
# For example, a 2D tensor product set can only support sizes of the form n1 * n2. Its native size may be
# (n1,n2) and its linear size n1*n2, but not any integer n maps to a native size tuple.
# By convention, we denote a native size variable by size_n.
"Compute the native size best corresponding to the given linear size."
approximate_native_size(d::Dictionary, size_l) = size_l

"Compute the linear size corresponding to the given native size."
linear_size(d::Dictionary, size_n) = size_n

"Suggest a suitable size, close to `n`, to resize the given dictionary."
approx_length(d::Dictionary, n::Int) = n
approx_length(d::Dictionary, n::Real) = approx_length(d, round(Int,n))

###############################
## Properties of function sets
###############################

"""
Does the dictionary implement a differentiation operator?
An optional second argument may specify an exact order of the derivative.
"""
hasderivative(Φ::Dictionary) = false

function hasderivative(Φ::Dictionary, order)
    # We have to be correct for zero order derivatives, because even if the
    # dictionary itself does not support derivatives, it could be part of a
    # composition where other elements do
    if orderiszero(order)
        true
    else
        order == 1 ? hasderivative(Φ) : false
    end
end

"Does the dictionary implement an antiderivative?"
hasantiderivative(Φ::Dictionary) = false

function hasantiderivative(Φ::Dictionary, order)
    if orderiszero(order)
        true
    else
        order == 1 ? hasantiderivative(Φ) : false
    end
end

"Does the dictionary have an associated interpolation grid?"
hasinterpolationgrid(d::Dictionary) = false

"Does the dictionary have a transform associated with some space?"
hastransform(d1::Dictionary, d2) = false

# Convenience functions: default grid, and conversion from grid to space
hastransform(d::Dictionary) = hasinterpolationgrid(d) && hastransform(d, interpolation_grid(d))
hastransform(d::Dictionary, grid::AbstractGrid) =
    hastransform(d, GridBasis{coefficienttype(d)}(grid))

"Does the dictionary support extension and restriction operators?"
hasextension(d::Dictionary) = false

# A concrete Dictionary may also override extension_set and restriction_set
# The default is simply to resize.
extension_set(d::Dictionary, n) = resize(d, n)
restriction_set(d::Dictionary, n) = resize(d, n)

###############################
## Iterating over dictionaries
###############################

# Default iterator over sets of functions: based on underlying index iterator.
function iterate(d::Dictionary)
    iter = eachindex(d)
    first_item, first_state = iterate(iter)
    (d[first_item], (iter, (first_item, first_state)))
end

function iterate(d::Dictionary, state)
    iter, iter_tuple = state
    iter_item, iter_state = iter_tuple
    next_tuple = iterate(iter, iter_state)
    if next_tuple != nothing
        next_item, next_state = next_tuple
        (d[next_item], (iter,next_tuple))
    end
end

include("dict_evaluation.jl")
include("dict_moments.jl")


text:
"""
    getAccumulatedTradeSeries(arr)
generates the profit and loss series by a given array containing profit and losses
be aware: data must be sorted by date upfront!
@todo: generalize type definition
"""
function getAccumulatedProfitLossSeries(profitLossSeriesArray::Array{Float64,1})
    accumulatedProfitAndLossSeries=accumulate(+, profitLossArray);
    return accumulatedProfitAndLossSeries
end

max_stars_repo_name: RainerNeu/tradeStatistics

max_issues_repo_name: RainerNeu/tradeStatistics

max_forks_repo_name: RainerNeu/tradeStatistics


text:
richprint(io::IO, dom::Vector{Node}, @nospecialize(x)) = richprint_default(io, dom, x)

function richprint_default(io::IO, dom::Vector{Node}, @nospecialize(x))
    t = typeof(x)
    isbitstype(t) || (Base.show_circular(io, x) && return)
    let id = startNode(dom, io, x, "")
        if which(Base.show, (typeof(io), typeof(x))).sig === Tuple{typeof(Base.show), IO, ANY}
            # override Base.default_show behavior
            nf = nfields(t)
            nb = sizeof(x)
            if nf != 0 || nb == 0
                Base.show(io, t)
                print(io, '(')
                recur_io = IOContext(io, :SHOWN_SET => x)
                for i in 1:nf
                    nextChild(dom, io, id)
                    f = fieldname(t, i)
                    if !isdefined(x, f)
                        print(io, undef_ref_str)
                    else
                        richprint(recur_io, dom, getfield(x, f))
                    end
                    if i < nf
                        print(io, ", ")
                    end
                end
                print(io, ')')
            else
                Base.show_default(io, x)
            end
        else
            Base.show(io, x)
        end
        endNode(dom, io, id)
    end
end

max_stars_repo_name: vtjnash/Firelight


text:
# ---
# title: 881. Boats to Save People
# id: problem881
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Two Pointers, Greedy
# link: <https://leetcode.com/problems/boats-to-save-people/description/>
# hidden: true
# ---
#
# The `i`-th person has weight `people[i]`, and each boat can carry a maximum
# weight of `limit`.
#
# Each boat carries at most 2 people at the same time, provided the sum of the
# weight of those people is at most `limit`.
#
# Return the minimum number of boats to carry every given person.  (It is
# guaranteed each person can be carried by a boat.)
#
#
#
# **Example 1:**
#
#
#
#     Input: people = [1,2], limit = 3
#     Output: 1
#     Explanation: 1 boat (1, 2)
#
#
# **Example 2:**
#
#
#
#     Input: people = [3,2,2,1], limit = 3
#     Output: 3
#     **Explanation** : 3 boats (1, 2), (2) and (3)
#
#
# **Example 3:**
#
#
#
#     Input: people = [3,5,3,4], limit = 5
#     Output: 4
#     **Explanation** : 4 boats (3), (3), (4), (5)
#
# **Note** :
#
#   * `1 <= people.length <= 50000`
#   * `1 <= people[i] <= limit <= 30000`
#
#
## @lc code=start
using LeetCode

function num_rescue_boats(people::Vector{Int}, limit::Int)
    sort!(people)
    a, b = 1, length(people)
    res = 0
    while a <= b
        res += 1
        a == b && break
        b -= 1
        (people[a] + people[b] <= limit) && (a += 1)
    end
    res
end
## @lc code=end


text:
module ProximityEnrichment

using GenomicFeatures, Statistics, StatsBase, HypothesisTests, Distributions

export proxenrich

"""
    intervaldist(intA, intB)

    intA        Genomic Interval A
    intB        Genomic Interval B

    returns distance between two intervals, 0 if intersecting.
"""
function intervaldist(intA, intB)
    (seqname(intA) != seqname(intB)) && return typemax(Int)

    ΔAB = leftposition(intB) - rightposition(intA)
    ΔBA = leftposition(intA) - rightposition(intB)

    (ΔBA ≤ 0) && (ΔAB ≤ 0) && return 0

    min(abs(ΔAB), abs(ΔBA))
end

"""
    closestinterval(intervals, i)

    intervals       Set of genomic locations
    i               Single genomic location

    returns distance between i and member of intervals closest to i
"""
function closestinterval(intervals, i)

    ind = searchsorted(intervals, i)

    Δ = typemax(Int)
    (1 ≤ ind.start ≤ length(intervals)) && (Δ = min(Δ, intervaldist(intervals[ind.start], i)))
    (1 ≤ ind.stop  ≤ length(intervals)) && (Δ = min(Δ, intervaldist(intervals[ind.stop],  i)))

    ifelse(Δ == typemax(Int), -1, Δ)
end

"""
    closest_tss_peak(gene_intervals, peak_intervals)

    gene_intervals   Intervals describing TSS
    peak_intervals   Intervals describing peak positions

    returns vector of distances
"""
closest_tss_peak(gene_intervals, peak_intervals) = [closestinterval(peak_intervals, iv) for iv in gene_intervals]

"""
    getx(d, xp)

    Aligns distances between peaks and genes that exist with positions to calculate fisher

    returns new_x and indexing of xp
"""
function getx(d, xp)
    ud = sort!(unique(d[d .> xp[1]]))
    dud = diff(ud)
    δ = step(xp)

    ind = findfirst(d -> d < δ, dud)

    if isnothing(ind)
        return xp, ind
    else
        xpi = Int(floor((ud[ind] - xp[1])/δ)) + 2
        return [ud[1:ind] ; xp[xpi:end]], ind
    end
end

"""
    fishertest(indA, indB)

    indA   logical vector selecting items in list A
    indB   logical vector selecting items in list B

    returns pvalue of right tail and oddsratio
"""
function fishertest(indA, indB; tail=:right)
    a = sum( indA  .&  indB )
    b = sum( indA  .& .!indB )
    c = sum(.!indA .&  indB )
    d = sum(.!indA .& .!indB )

    if iszero(a) || iszero(b) || iszero(c) || iszero(d)
       return 1.0, 1.0
    end
    or = (a/c)/(b/d)

    pvalue(FisherExactTest(a, b, c, d), tail=tail), or
end

"""
    hypertest(indA, indB)

    indA   logical vector selecting items in list A
    indB   logical vector selecting items in list B

    returns pvalue of right tail and oddsratio
"""
function hypertest(indA, indB; tail=:right)

    ### Hypergeometric distribution is parametised by
    ## success in population (number of positives in B)
    ## failures in population (numnber of negatives in B)
    ## Number of trials (total positives in A)
    ## pvalue is ccdf of k where k is the number of positves in A and B
    s = sum(indB)          ## sucesses
    f = length(indB) - s   ## failures
    n = sum(indA)          ## trials

    k = sum( indA  .&  indB )
    b = sum( indA  .& .!indB )
    c = sum(.!indA .&  indB )
    d = sum(.!indA .& .!indB )

    if iszero(k) || iszero(b) || iszero(c) || iszero(d)
       return 1.0, 1.0
    end

    if tail == :right
        hp = ccdf(Hypergeometric(s, f, n), k - 1)
    elseif tail == :left
        hp = cdf(Hypergeometric(s, f, n), k)
    end

    or = (k/c)/(b/d)
    hp, or
end

"""
    proxenrich(xp, genes, peaks, geneind=trues(length(genes)), peakind=trues(length(peaks)); trans = x -> log10(x + 1), label="")

    xp          points to calculate enrichment at
    genes       Vector of Intervals describing TSS coordinates
    peaks       Vector of Intervals describing peaks
    geneind     logical vector selecting the foreground, i.e. genes of interest
    peakind     logical vector selecting the peaks to

    Calculates right tail pvalues using hypergeometric or FisherExactTest

    returns named tuple of relevant stats

"""
function proxenrich(xp, genes, peaks, geneind=trues(length(genes)), peakind=trues(length(peaks)); trans = x -> log10(x + 1), label="", testfun=hypertest, tail=:right)
    ### Calculate distance closest peak to each tss

    !issorted(genes) && error("Genes not sorted")
    !issorted(peaks) && error("Peaks not sorted")
    db = closest_tss_peak(genes, peaks[peakind])
    validdist = db .!= -1 ### closest_tss_peak returns -1 if no peaks on same chromosome
    Δb = trans.(db)     ## background
    Δg = Δb[geneind]    ## foreground

    ### Calculate subset of xp to avoid calculating pvalues at distances below/above closest/furtherest gene
    xg, xgi = getx(Δg, xp)

    ### enrichment over background, not necessary for fisher calculation, but useful to plot
    ## Remove Calculation of ecdf for now
    # fc = ecdf(Δg)(xg)./ecdf(Δb)(xg)

    ### Fisher tests for number of genes within δ against gene list
    count = [sum(geneind .& (Δb .≤ δ) .& validdist) for δ ∈ xg]
    ft    = [testfun(geneind, (Δb .≤ δ) .& validdist, tail=tail) for δ ∈ xg]

    pvalue = first.(ft)
    or     = last.(ft)
    (xg=xg, xgi=xgi, pvalue=pvalue, or=or, n=sum(geneind), np=size(peaks[peakind]), count=count, label=label)
end

end # module

max_stars_repo_name: owensnick/ProximityEnrichment

max_issues_repo_name: owensnick/ProximityEnrichment

max_forks_repo_name: owensnick/ProximityEnrichment


text:
"""
StagedPolicy

Policy for finite horizon POMDP solvers.
Wraps multiple Policies, one for each stage.
"""
struct StagedPolicy{M<:Union{POMDP,MDP}, P<:Policy}<:Policy
    pomdp::M
    staged_policies::Array{P, 1}

    function StagedPolicy(m::M, p::Vector{P}) where {M<:Union{POMDP,MDP}, P<:Policy}
        HorizonLength(m) == InfiniteHorizon() ? throw(ArgumentError("Argument pomdp should be
                             valid Finite Horizon POMDP with methods from FiniteHorizonPOMDPs.jl
                             interface implemented. If you are completely sure that you implemented
                             all of them, you should also check if you have defined
                             HorizonLength(::Type{<:MyFHMDP})")) : new{M, P}(m, p)
     end
end

### functional methods ###
value(p::StagedPolicy, b, t) = t <= horizon(p.pomdp) ? value(p.staged_policies[t], b) : 0.

action(p::StagedPolicy, b, t) = t <= horizon(p.pomdp) ? action(p.staged_policies[t], b) : ordered_actions(p.pomdp)[1]

# Given that, for example, POMDPSimulators.jl does not work with action(p, b, t),
# we need to somehow introduce the stage variable
function action(p::StagedPolicy, b::DiscreteBelief)
    t = -1
    for i in 1:length(b.b)
        if b.b[i] > 0.0
            if t < 0
                t = stage(p.pomdp, b.state_list[i])
            elseif t != stage(p.pomdp, b.state_list[i])
                error("Not all states in the belief had the same stage")
            end
        end
    end
    return action(p, b, t)
end

updater(p::StagedPolicy) = DiscreteUpdater(p.pomdp)

"""
Return an iterator of alpha vector-action pairs of given stage in the policy.
"""
alphapairs(p::StagedPolicy, t::Int64) = alphapairs(p.staged_policies[t])

"""
Return the alpha vectors of given stage
"""
alphavectors(p::StagedPolicy, t::Int64) = alphavectors(p.staged_policies[t])

actionvalues(p::StagedPolicy, b, t::Int64) = actionvalues(p.staged_policies[t], b)

function Base.push!(p::StagedPolicy{<:Union{POMDP,MDP},<:Policy}, alpha::Vector{Float64}, a, t)
    push!(p.staged_policies[t].alphas, alpha)
    push!(p.staged_policies[t].action_map, a)
end


text:
using ITensors,
      Test

include("util.jl")

function basicRandomMPO(sites; dim=4)
  M = MPO(sites)
  N = length(M)
  links = [Index(dim,"n=$(n-1),Link") for n=1:N+1]
  for n=1:N
    M[n] = randomITensor(links[n],sites[n],sites[n]',links[n+1])
  end
  M[1] *= delta(links[1])
  M[N] *= delta(links[N+1])
  return M
end

@testset "MPO Basics" begin
  N = 6
  sites = [Index(2,"Site,n=$n") for n=1:N]
  @test length(MPO()) == 0
  O = MPO(sites)
  @test length(O) == N

  str = split(sprint(show, O), '\n')
  @test str[1] == "MPO"
  @test length(str) == length(O) + 2

  O[1] = ITensor(sites[1], prime(sites[1]))
  @test hasind(O[1],sites[1])
  @test hasind(O[1],prime(sites[1]))
  P = copy(O)
  @test hasind(P[1],sites[1])
  @test hasind(P[1],prime(sites[1]))
  # test constructor from Vector{ITensor}
  K = randomMPO(sites)
  @test ITensors.data(MPO(copy(ITensors.data(K)))) == ITensors.data(K)

  @testset "orthogonalize!" begin
    phi = randomMPS(sites)
    K = randomMPO(sites)
    orthogonalize!(phi, 1)
    orthogonalize!(K, 1)
    orig_inner = ⋅(phi, K, phi)
    orthogonalize!(phi, div(N, 2))
    orthogonalize!(K, div(N, 2))
    @test ⋅(phi, K, phi) ≈ orig_inner
  end

  @testset "inner <y|A|x>" begin
    phi = randomMPS(sites)
    K = randomMPO(sites)
    @test maxlinkdim(K) == 1
    psi = randomMPS(sites)
    phidag = dag(phi)
    prime!(phidag)
    phiKpsi = phidag[1]*K[1]*psi[1]
    for j = 2:N
      phiKpsi *= phidag[j]*K[j]*psi[j]
    end
    @test phiKpsi[] ≈ inner(phi,K,psi)

    badsites = [Index(2,"Site") for n=1:N+1]
    badpsi = randomMPS(badsites)
    @test_throws DimensionMismatch inner(phi,K,badpsi)

    # make bigger random MPO...
    for link_dim in 2:5
        mpo_tensors  = ITensor[ITensor() for ii in 1:N]
        mps_tensors  = ITensor[ITensor() for ii in 1:N]
        mps_tensors2 = ITensor[ITensor() for ii in 1:N]
        mpo_link_inds = [Index(link_dim, "r$ii,Link") for ii in 1:N-1]
        mps_link_inds = [Index(link_dim, "r$ii,Link") for ii in 1:N-1]
        mpo_tensors[1] = randomITensor(mpo_link_inds[1], sites[1], sites[1]')
        mps_tensors[1] = randomITensor(mps_link_inds[1], sites[1])
        mps_tensors2[1] = randomITensor(mps_link_inds[1], sites[1])
        for ii in 2:N-1
            mpo_tensors[ii] = randomITensor(mpo_link_inds[ii], mpo_link_inds[ii-1], sites[ii], sites[ii]')
            mps_tensors[ii] = randomITensor(mps_link_inds[ii], mps_link_inds[ii-1], sites[ii])
            mps_tensors2[ii] = randomITensor(mps_link_inds[ii], mps_link_inds[ii-1], sites[ii])
        end
        mpo_tensors[N] = randomITensor(mpo_link_inds[N-1], sites[N], sites[N]')
        mps_tensors[N] = randomITensor(mps_link_inds[N-1], sites[N])
        mps_tensors2[N] = randomITensor(mps_link_inds[N-1], sites[N])
        K   = MPO(mpo_tensors, 0, N+1)
        psi = MPS(mps_tensors, 0, N+1)
        phi = MPS(mps_tensors2, 0, N+1)
        orthogonalize!(psi, 1; maxdim=link_dim)
        orthogonalize!(K, 1; maxdim=link_dim)
        orthogonalize!(phi, 1; normalize=true, maxdim=link_dim)
        phidag = dag(phi)
        prime!(phidag)
        phiKpsi = phidag[1]*K[1]*psi[1]
        for j = 2:N
          phiKpsi *= phidag[j]*K[j]*psi[j]
        end
        @test scalar(phiKpsi) ≈ inner(phi,K,psi)
    end
  end

  @testset "inner <By|A|x>" begin
    phi = makeRandomMPS(sites)

    K = makeRandomMPO(sites,chi=2)
    J = makeRandomMPO(sites,chi=2)

    psi = makeRandomMPS(sites)
    phidag = dag(phi)
    prime!(phidag, 2)
    Jdag = dag(J)
    prime!(Jdag)
    for j ∈ eachindex(Jdag)
      swapprime!(Jdag[j],2,3)
      swapprime!(Jdag[j],1,2)
      swapprime!(Jdag[j],3,1)
    end

    phiJdagKpsi = phidag[1]*Jdag[1]*K[1]*psi[1]
    for j ∈ eachindex(psi)[2:end]
      phiJdagKpsi = phiJdagKpsi*phidag[j]*Jdag[j]*K[j]*psi[j]
    end

    @test phiJdagKpsi[] ≈ inner(J,phi,K,psi)

    badsites = [Index(2,"Site") for n=1:N+1]
    badpsi = randomMPS(badsites)
    @test_throws DimensionMismatch inner(J,phi,K,badpsi)
  end

  @testset "error_contract" begin
    phi = makeRandomMPS(sites)
    K = makeRandomMPO(sites,chi=2)

    psi = makeRandomMPS(sites)

    dist = sqrt(abs(1 + (inner(phi,phi) - 2*real(inner(phi,K,psi)))
                        /inner(K,psi,K,psi)))
    @test dist ≈ error_contract(phi,K,psi)

    badsites = [Index(2,"Site") for n=1:N+1]
    badpsi = randomMPS(badsites)
    # Apply K to phi and check that error_contract is close to 0.
    Kphi = contract(K, phi; method="naive", cutoff=1E-8)
    @test error_contract(Kphi, K, phi) ≈ 0. atol=1e-4

    @test_throws DimensionMismatch contract(K,badpsi;method="naive", cutoff=1E-8)
    @test_throws DimensionMismatch error_contract(phi,K,badpsi)
  end

  @testset "contract" begin
    phi = randomMPS(sites)
    K   = randomMPO(sites)
    @test maxlinkdim(K) == 1
    psi = randomMPS(sites)
    psi_out = contract(K, psi; maxdim=1)
    @test inner(phi,psi_out) ≈ inner(phi,K,psi)
    @test_throws ArgumentError contract(K, psi, method="fakemethod")

    badsites = [Index(2,"Site") for n=1:N+1]
    badpsi = randomMPS(badsites)
    @test_throws DimensionMismatch contract(K,badpsi)

    # make bigger random MPO...
    for link_dim in 2:5
      mpo_tensors  = ITensor[ITensor() for ii in 1:N]
      mps_tensors  = ITensor[ITensor() for ii in 1:N]
      mps_tensors2 = ITensor[ITensor() for ii in 1:N]
      mpo_link_inds = [Index(link_dim, "r$ii,Link") for ii in 1:N-1]
      mps_link_inds = [Index(link_dim, "r$ii,Link") for ii in 1:N-1]
      mpo_tensors[1] = randomITensor(mpo_link_inds[1], sites[1], sites[1]')
      mps_tensors[1] = randomITensor(mps_link_inds[1], sites[1])
      mps_tensors2[1] = randomITensor(mps_link_inds[1], sites[1])
      for ii in 2:N-1
        mpo_tensors[ii] = randomITensor(mpo_link_inds[ii], mpo_link_inds[ii-1], sites[ii], sites[ii]')
        mps_tensors[ii] = randomITensor(mps_link_inds[ii], mps_link_inds[ii-1], sites[ii])
        mps_tensors2[ii] = randomITensor(mps_link_inds[ii], mps_link_inds[ii-1], sites[ii])
      end
      mpo_tensors[N] = randomITensor(mpo_link_inds[N-1], sites[N], sites[N]')
      mps_tensors[N] = randomITensor(mps_link_inds[N-1], sites[N])
      mps_tensors2[N] = randomITensor(mps_link_inds[N-1], sites[N])
      K   = MPO(mpo_tensors, 0, N+1)
      psi = MPS(mps_tensors, 0, N+1)
      phi = MPS(mps_tensors2, 0, N+1)
      orthogonalize!(psi, 1; maxdim=link_dim)
      orthogonalize!(K, 1; maxdim=link_dim)
      orthogonalize!(phi, 1; normalize=true, maxdim=link_dim)
      psi_out = contract(deepcopy(K), deepcopy(psi); maxdim=10*link_dim, cutoff=0.0)
      @test inner(phi, psi_out) ≈ inner(phi, K, psi)
    end
  end

  @testset "add(::MPO, ::MPO)" begin
    shsites = siteinds("S=1/2", N)
    K = randomMPO(shsites)
    L = randomMPO(shsites)
    M = add(K, L)
    @test length(M) == N
    psi = randomMPS(shsites)
    k_psi = contract(K, psi, maxdim=1)
    l_psi = contract(L, psi, maxdim=1)
    @test inner(psi, k_psi + l_psi) ≈ ⋅(psi, M, psi) atol=5e-3
    @test inner(psi, sum([k_psi, l_psi])) ≈ dot(psi, M, psi) atol=5e-3
    for dim in 2:4
        shsites = siteinds("S=1/2",N)
        K = basicRandomMPO(shsites; dim=dim)
        L = basicRandomMPO(shsites; dim=dim)
        M = K + L
        @test length(M) == N
        psi = randomMPS(shsites)
        k_psi = contract(K, psi)
        l_psi = contract(L, psi)
        @test inner(psi, k_psi + l_psi) ≈ dot(psi, M, psi) atol=5e-3
        @test inner(psi, sum([k_psi, l_psi])) ≈ inner(psi, M, psi) atol=5e-3
        psi = randomMPS(shsites)
        M = add(K, L; cutoff=1E-9)
        k_psi = contract(K, psi)
        l_psi = contract(L, psi)
        @test inner(psi, k_psi + l_psi) ≈ inner(psi, M, psi) atol=5e-3
    end
  end

  @testset "contract(::MPO, ::MPO)" begin
    psi = randomMPS(sites)
    K = randomMPO(sites)
    L = randomMPO(sites)
    @test maxlinkdim(K) == 1
    @test maxlinkdim(L) == 1
    KL = contract(prime(K), L, maxdim=1)
    psi_kl_out = contract(prime(K), contract(L, psi, maxdim=1), maxdim=1)
    @test inner(psi,KL,psi) ≈ inner(psi, psi_kl_out) atol=5e-3

    # where both K and L have differently labelled sites
    othersitesk = [Index(2,"Site,aaa") for n=1:N]
    othersitesl = [Index(2,"Site,bbb") for n=1:N]
    K = randomMPO(sites)
    L = randomMPO(sites)
    for ii in 1:N
      replaceind!(K[ii], sites[ii]', othersitesk[ii])
      replaceind!(L[ii], sites[ii]', othersitesl[ii])
    end
    KL = contract(K, L, maxdim=1)
    psik = randomMPS(othersitesk)
    psil = randomMPS(othersitesl)
    psi_kl_out = contract(K, contract(L, psil, maxdim=1), maxdim=1)
    @test inner(psik,KL,psil) ≈ inner(psik, psi_kl_out) atol=5e-3

    badsites = [Index(2,"Site") for n=1:N+1]
    badL = randomMPO(badsites)
    @test_throws DimensionMismatch contract(K,badL)
  end

  @testset "*(::MPO, ::MPO)" begin
    psi = randomMPS(sites)
    K = randomMPO(sites)
    L = randomMPO(sites)
    @test maxlinkdim(K) == 1
    @test maxlinkdim(L) == 1
    KL = *(prime(K), L, maxdim=1)
    psi_kl_out = *(prime(K), *(L, psi, maxdim=1), maxdim=1)
    @test ⋅(psi, KL, psi) ≈ dot(psi, psi_kl_out) atol=5e-3

    # where both K and L have differently labelled sites
    othersitesk = [Index(2,"Site,aaa") for n=1:N]
    othersitesl = [Index(2,"Site,bbb") for n=1:N]
    K = randomMPO(sites)
    L = randomMPO(sites)
    for ii in 1:N
      replaceind!(K[ii], sites[ii]', othersitesk[ii])
      replaceind!(L[ii], sites[ii]', othersitesl[ii])
    end
    KL = *(K, L, maxdim=1)
    psik = randomMPS(othersitesk)
    psil = randomMPS(othersitesl)
    psi_kl_out = *(K, *(L, psil, maxdim=1), maxdim=1)
    @test dot(psik, KL, psil) ≈ psik ⋅ psi_kl_out atol=5e-3

    badsites = [Index(2,"Site") for n=1:N+1]
    badL = randomMPO(badsites)
    @test_throws DimensionMismatch K * badL
  end

  sites = siteinds("S=1/2",N)
  O = MPO(sites,"Sz")
  @test length(O) == N # just make sure this works

  @test_throws ArgumentError randomMPO(sites, 2)
  @test isnothing(linkind(MPO(fill(ITensor(), N), 0, N + 1), 1))
end

nothing


text:
"""
    ADI.Metrics

This module provides code for analyzing the results from ADI in a way that is interpretable statistically. Some of the key functionalities are signal-to-noise, significance, the receiver operating characteristic, and the contrast curve.
"""
module Metrics

export detectionmap,
       snr,
       significance,
       noise,
       contrast_curve,
       throughput,
       stimmap,
       slimmap,
       stim_threshold

include("contrast.jl")
include("snr.jl")
include("stim.jl")
include("ensemble.jl")

end