text:
module Clockwork

export Clock
import Base: +, *, ^, -

immutable Clock
    hour::Int
    on_minute::Bool
end

Clock(hour::Integer) = Clock(Int(hour), false)

const glyph_for_clock = Dict{Clock, Symbol}()
for i=1:12
    s1 = Symbol(Char(0x1F550+(i-1)))
    s2 = Symbol(Char(0x1F550+(i-1)+12))
    t1 = Clock(i, false)
    t2 = Clock(i, true)
    @eval const $s1=$t1
    @eval const $s2=$t2
    glyph_for_clock[t1] = s1
    glyph_for_clock[t2] = s2
    @eval export $s1, $s2
end

for (op, carry) in [(:+, 1), (:-, -1)]
    @eval begin
        function $op(c1::Clock, c2::Clock)
            if c1.on_minute && c2.on_minute
                hour, on_minute = $op(c1.hour, c2.hour)+$carry, false
            elseif !c1.on_minute && !c2.on_minute
                hour, on_minute = $op(c1.hour, c2.hour), false
            else
                hour, on_minute = $op(c1.hour, c2.hour), true
            end
            Clock(mod1(hour,12), on_minute)
        end
    end
end

for (op, inner_op) in [(:*,:+), (:^, :*)]
    @eval begin
        function $op(c1::Clock, c2::Clock)
            c = c1
            for n=1:c2.hour-1
                c = $inner_op(c, c1)
            end
            c
        end
    end
end

for op in [:+, :*, :^, :-]
    @eval $op(n::Integer, c::Clock) = $op(Clock(n), c)
    @eval $op(c::Clock, n::Integer) = $op(c, Clock(n))
end

function Base.show(io::IO, c::Clock)
    print(io, glyph_for_clock[c])
end

function Base.convert(::Type{Dates.CompoundPeriod}, c::Clock)
    d = Dates.Hour(c.hour)
    if c.on_minute
        d += Dates.Minute(30)
    end
    Dates.CompoundPeriod(d)
end

function Base.convert(::Type{DateTime}, c::Clock)
    args = Any[Dates.Hour(c.hour)]
    if c.on_minute
        push!(args, Dates.Minute(30))
    end
    DateTime(args...)
end

Base.DateTime(c::Clock) = convert(DateTime, c)

_hour(d::Dates.DateTime) = Dates.hour(d)
_minute(d::Dates.DateTime) = Dates.minute(d)
_hour(d::Dates.Hour) = Dates.value(d)
_minute(d::Dates.Minute) = Dates.value(d)
_hour(x) = 0
_minute(x) = 0
for op in [:_hour, :_minute]
    @eval $op(d::Dates.CompoundPeriod) = sum(map($op, d.periods))
end

Base.convert(::Type{Clock}, d::Dates.AbstractTime) =
  Clock(mod1(_hour(d), 12), _minute(d)==30)

function Base.convert(::Type{Clock}, fmt::AbstractString)
    dt = DateTime(fmt,"H:M")
    convert(Clock, dt)
end

Base.zero(::Type{Clock}) = Clock(12, false)
Base.one(::Type{Clock}) = Clock(1, false)
Base.zero(::Clock) = zero(Clock)
Base.one(::Clock) = one(Clock)

Base.rand(rng, ::Type{Clock}) = Clock(rand(rng, 1:12), rand(rng, Bool))

end # module


text:
# Test spatial weights matrix
@testset "Choropleth Maps" begin

    guerry = sdataset("Guerry")
    boston = sdataset("Bostonhsg")

    @testset "EqualIntervals" begin
        mc = mapclassify(EqualIntervals(5), guerry.Litercy)

        @test length(mc) == 5
        @test counts(mc) == [20, 20, 21, 12, 12]
        lbound, ubound = bounds(mc)
        @test lbound == [12, 24.4, 36.8, 49.2, 61.6]
        @test ubound == [24.4, 36.8, 49.2, 61.6, 74]
        @test mc.lower == :open
        @test mc.upper == :closed
        @test levels(mc) == string.(1:5)
        @test assignments(mc)[1:3] == [3, 4, 1]
        @test assignments(mc)[83:85] == [1, 5, 3]
    end

    @testset "Quantiles" begin
        mc = mapclassify(Quantiles(5), guerry.Litercy)

        @test length(mc) == 5
        @test counts(mc) == [18, 19, 15, 16, 17]
        lbound, ubound = bounds(mc)
        @test lbound ≈ [12, 23, 31, 43, 54.4]
        @test ubound ≈ [23, 31, 43, 54.4, 74]
        @test mc.lower == :open
        @test mc.upper == :closed
    end

    @testset "CustomBreaks" begin
        mc = mapclassify(CustomBreaks([20, 30, 60]), guerry.Litercy)

        @test length(mc) == 4
        @test counts(mc) == [14, 18, 41, 12]
        lbound, ubound = bounds(mc)
        @test lbound ≈ [12, 20, 30, 60]
        @test ubound ≈ [20, 30, 60, 74]
        @test mc.lower == :open
        @test mc.upper == :closed
    end

    @testset "NaturalBreaks" begin
        mc = mapclassify(NaturalBreaks(5), guerry.Litercy)

        @test length(mc) == 5
        @test counts(mc) == [15, 24, 22, 15, 9]
        lbound, ubound = bounds(mc)
        @test lbound == [12, 21, 34, 49, 63.0]
        @test ubound == [21, 34, 49, 63, 74.0]
        @test mc.lower == :open
        @test mc.upper == :closed
    end

    @testset "BoxPlot" begin
        mc = mapclassify(BoxPlot(), guerry.Litercy)

        @test length(mc) == 6
        @test counts(mc) == [0, 23, 21, 20 ,21, 0]
        lbound, ubound = bounds(mc)
        @test lbound == [-Inf, -15.5, 25, 38, 52, 92.5]
        @test ubound == [-15.5, 25, 38, 52, 92.5, Inf]
        @test mc.lower == :open
        @test mc.upper == :closed
    end

    @testset "StdMean" begin
        mc = mapclassify(StdMean(), guerry.Litercy)

        @test length(mc) == 6
        @test counts(mc) == [0, 15, 30, 24, 16, 0]
        lbound, ubound = bounds(mc)
        @test lbound ≈ [-Inf, 4.272507, 21.706842, 39.141176, 56.575511, 74.009846] atol = 1e-5
        @test ubound ≈ [4.272507, 21.706842, 39.141176, 56.575511, 74.009846, Inf] atol = 1e-5
        @test mc.lower == :open
        @test mc.upper == :closed
    end

    @testset "Percentiles" begin
        mc = mapclassify(Percentiles(), guerry.Litercy)

        @test length(mc) == 6
        @test counts(mc) == [1, 9, 34, 32, 8, 1]
        lbound, ubound = bounds(mc)
        @test lbound ≈ [12, 12.84, 18, 38, 65.4, 73.16] atol = 1e-5
        @test ubound ≈ [12.84, 18, 38, 65.4, 73.16, 74] atol = 1e-5
        @test mc.lower == :open
        @test mc.upper == :closed
    end

    @testset "Closed-Open" begin
        mc = mapclassify(Quantiles(5), guerry.Litercy, lower = :closed, upper = :open)

        @test length(mc) == 5
        @test counts(mc) == [15, 17, 18, 18, 17]
        lbound, ubound = bounds(mc)
        @test lbound ≈ [12, 23, 31, 43, 54.4]
        @test ubound ≈ [23, 31, 43, 54.4, 74]
        @test mc.lower == :closed
        @test mc.upper == :open
    end

    @testset "Graduated Labels" begin
        mc = mapclassify(Quantiles(5), guerry.Litercy)
        @test maplabels(mc) ==
            ["[12.0, 23.0] (18)"
            "(23.0, 31.0] (19)"
            "(31.0, 43.0] (15)"
            "(43.0, 54.4] (16)"
            "(54.4, 74.0] (17)"]

        @test maplabels(mc, sep = " -- ") ==
            ["[12.0 -- 23.0] (18)"
            "(23.0 -- 31.0] (19)"
            "(31.0 -- 43.0] (15)"
            "(43.0 -- 54.4] (16)"
            "(54.4 -- 74.0] (17)"]

        @test maplabels(mc, counts = false) ==
            ["[12.0, 23.0]"
            "(23.0, 31.0]"
            "(31.0, 43.0]"
            "(43.0, 54.4]"
            "(54.4, 74.0]"]

        mc = mapclassify(StdMean(), guerry.Litercy)
        @test maplabels(mc, digits = 4) ==
            ["[-Inf, 4.2725] (0)"
            "(4.2725, 21.7068] (15)"
            "(21.7068, 39.1412] (30)"
            "(39.1412, 56.5755] (24)"
            "(56.5755, 74.0098] (16)"
            "(74.0098, Inf] (0)"]

        mc = mapclassify(Quantiles(5), guerry.Litercy, lower = :closed, upper = :open)
        @test maplabels(mc) ==
            ["[12.0, 23.0) (15)"
            "[23.0, 31.0) (17)"
            "[31.0, 43.0) (18)"
            "[43.0, 54.4) (18)"
            "[54.4, 74.0] (17)"]
    end

    @testset "Unique" begin
        mc = mapclassify(Unique(), guerry.Region)

        @test length(mc) == 5
        @test counts(mc) == [17, 17, 17, 17, 17]
        @test levels(mc) == ["E", "N", "C", "S", "W"]
    end

    @testset "Unique Labels" begin
        mc = mapclassify(Unique(), guerry.Region)
        @test maplabels(mc) ==
            ["E (17)"
            "N (17)"
            "C (17)"
            "S (17)"
            "W (17)"]

        @test maplabels(mc, counts = false) ==
            ["E"
            "N"
            "C"
            "S"
            "W"]
    end

    @testset "CoLocation" begin
        mc = mapclassify(CoLocation(StdMean()), [guerry.Litercy, guerry.Donatns])

        @test length(mc) == 7
        @test counts(mc) == [61, 0, 0, 19, 3, 2, 0]
        @test levels(mc) == ["Other", "1", "2", "3", "4", "5", "6"]

        @test maplabels(mc) ==
            ["Other (61)"
            "1 (0)"
            "2 (0)"
            "3 (19)"
            "4 (3)"
            "5 (2)"
            "6 (0)"]
        @test maplabels(mc, counts = false) ==
            ["Other"
            "1"
            "2"
            "3"
            "4"
            "5"
            "6"]
    end

    @testset "CoLocation String" begin
        LitHigh = repeat(["Low"], length(guerry.Litercy))
        LitHigh[guerry.Litercy .> mean(guerry.Litercy)] .= "High"
        DonHigh = repeat(["Low"], length(guerry.Donatns))
        DonHigh[guerry.Donatns .> mean(guerry.Donatns)] .= "High"

        mc = mapclassify(CoLocation(Unique()), [LitHigh, DonHigh])

        @test length(mc) == 3
        @test counts(mc) == [49, 10, 26]
        @test levels(mc) == ["Other", "High", "Low"]

        @test maplabels(mc) ==
            ["Other (49)"
            "High (10)"
            "Low (26)"]
        @test maplabels(mc, counts = false) ==
            ["Other"
            "High"
            "Low"]
    end

    @testset "Plot recipes" begin
        # Polygons
        @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), guerry.geometry, guerry.Litercy, EqualIntervals())
        # Points
        @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), boston.geometry, boston.MEDV, EqualIntervals())
        # Empty group in legend
        @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), boston.geometry, boston.MEDV, StdMean())
        # DataFrame and Symbol
        @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), guerry, :Litercy, EqualIntervals())
        # Paired palette with k <= 12
        @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), guerry, :Region, Unique())
        # Co-location
        @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), guerry.geometry, [guerry.Litercy, guerry.Donatns], CoLocation(StdMean()))
        # DataRrame Co-location
        @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), guerry, [:Litercy, :Donatns], CoLocation(StdMean()))
    end

end

ext: jl


text:
using GUI
using Test

include("googlefonts.jl")


text:
"""
    create_local_upstream(rootdir = tempname(); verbose = true, branch_name = "main")

Creates a local repo usable as a remote for testing.
Returns a tuple (;url, upstream_repo)
"""
function create_local_upstream(
        rootdir = tempname();
        verbose = true,
        branch_name = "main"
    )

    # home dirs
    upstream_repo = joinpath(rootdir, "upstream")
    _rm(upstream_repo)
    mkpath(upstream_repo)

    # create upstream
    verbose && @info("setting up upstream")
    url = _url_from_file(upstream_repo)
    _run("git -C $(upstream_repo) --bare init 2>&1"; verbose)
    verbose && println(read(joinpath(upstream_repo, "config"), String))
    verbose && println("\n")

    # make first commit
    cdir = tempname()
    mkpath(cdir)
    _run("git clone $(url) $(cdir) 2>&1"; verbose)
    dumpfile = joinpath(cdir, "README.md")
    write(dumpfile, "# TEST")
    _run("git -C $(cdir) add -A 2>&1"; verbose)
    commit_msg = "Test upstream created"
    user_name = get_global_config("user.name", "GitLink")
    user_email = get_global_config("user.email", "fake@email.com")
    _run("git -C $(cdir) -c user.name='$(user_name)' -c user.email='$(user_email)' commit -am '$(commit_msg)' 2>&1"; verbose)
    _run("git -C $(cdir) push 2>&1"; verbose)

    # rename branch
    curr_branch = _curr_branch(upstream_repo)
    if curr_branch != branch_name
        _run("git -C $(upstream_repo) branch -m $(curr_branch) $(branch_name) 2>&1"; verbose)
    end

    verbose && println("\n", "-"^60, "\n")
    _run("git -C $(upstream_repo) --no-pager log 2>&1"; verbose)

    _rm(cdir)

    return (;url, upstream_repo)
end


text:
push!(LOAD_PATH, string(pwd(),"/../../code/src"))
using ModifiedStokesSolver
using Base.Test

include("ModStokesFMM2D.jl")

alpha = 20.0

## define source and target points

ns = 500
nt = 600

srand(1)
stokeslet_str = 1 - 2*rand(2,ns)
src = 1 - 2*rand(2,ns)
fvec = 1 - 2*rand(2,ns)
nvec = 1 - 2*rand(2,ns)
targ = 1 - 2*rand(2,nt)

# Compute using the kernels from ModifiedStokesSolver
function mss_stokeslet_direct(src, targ, stokeslet_str, alpha)
    ns = size(src, 2)
    nt = size(targ, 2)
    u = zeros(2, nt)
    for i=1:ns
        for j=1:nt
            f = stokeslet_str[:, i]
            r = src[:,i]-targ[:,j]
            u[:,j] += stokeslet(r, f, alpha)
        end
    end
    return u
end
function mss_stresslet_direct(src, targ, fvec, nvec, alpha)
    ns = size(src, 2)
    nt = size(targ, 2)
    u = zeros(2, nt)
    for i=1:ns
        for j=1:nt
            r = src[:,i]-targ[:,j]
            u[:,j] += stresslet(r, fvec[:,i], nvec[:,i], alpha)
        end
    end
    return u
end

## Stokeslet
u1 = mss_stokeslet_direct(src, targ, stokeslet_str, alpha)
u2 = fmm_stokeslet_direct(src, targ, stokeslet_str, alpha)
u3 = fmm_stokeslet_targ(src, targ, stokeslet_str, alpha)

## Stresslet
d1 = mss_stresslet_direct(src, targ, fvec, nvec, alpha)
d2 = fmm_stresslet_direct(src, targ, fvec, nvec, alpha)
d3 = fmm_stresslet_targ(src, targ, fvec, nvec, alpha)

## Print + test
println(" * Stokeslet")
Gerrdirect = norm(u1[:]-u2[:]) / norm(u1[:], Inf)
Gerrfmm = norm(u1[:]-u3[:]) / norm(u1[:], Inf)
@show Gerrdirect
@show Gerrfmm

println(" * Stresslet")
Terrdirect = norm(d1[:]-d2[:]) / norm(d1[:], Inf)
Terrfmm = norm(d1[:]-d3[:]) / norm(d1[:], Inf)
@show Terrdirect
@show Terrfmm


text:
using Dates

this_year = year(today(Dates.today())
this_month = month(Dates.today())
this_day = day(Dates.today())

function read_ohlc(stock::String, fm::Int, fd::Int, fy::Int, tm::Int, td::Int, ty::Int, period::String)

# take care of yahoo's 0 indexing for month
    fm -= 1
    tm -= 1

    ydata     = readlines(`curl -s "http://ichart.finance.yahoo.com/table.csv?s=$stock&a=$fm&b=$fd&c=$fy&d=$tm&e=$td&f=$ty&g=$period"`)
    numstring = ydata[2:end]
    sa  = split(numstring[1], ',')'

    for i in 2:length(numstring)
        sa  = [sa ; split(numstring[i], ',')']
    end

    time_conversion = map(x -> parse("yyyy-MM-dd", x), convert(Array{UTF16String}, sa[:,1]))

    ohlc = DataFrame(quote
        Date  = $time_conversion
        Open  = float($sa[:,2])
        High  = float($sa[:,3])
        Low   = float($sa[:,4])
        Close = float($sa[:,5])
        Vol   =   int($sa[:,6])
        Adj   = float($sa[:,7])
    end)

    flipud(ohlc)

end

# monthly data

window = 3 # years

monthly_ohlc_over_window(stock::String) = read_ohlc(stock::String, this_month, this_day, this_year - window, this_month,  this_day, this_year, "m")

# TODO make the portfolio list unlimited in length
function returns{T <: Union(ASCIIString, UTF8String)}(list::Array{T})

    str1 = list[1]
    str2 = list[2]
    str3 = list[3]
    str4 = list[4]
    str5 = list[5]

    s1 = monthly_ohlc_over_window(list[1])
    s2 = monthly_ohlc_over_window(list[2])
    s3 = monthly_ohlc_over_window(list[3])
    s4 = monthly_ohlc_over_window(list[4])
    s5 = monthly_ohlc_over_window(list[5])

    r1 = log_return(s1["Adj"])
    r2 = log_return(s2["Adj"])
    r3 = log_return(s3["Adj"])
    r4 = log_return(s4["Adj"])
    r5 = log_return(s5["Adj"])

    ohlc = DataFrame(quote
        $str1 = $r1
        $str2 = $r2
        $str3 = $r3
        $str4 = $r4
        $str5 = $r5
    end)
    ohlc[2:end, :]
end

# function read_monthlies(filename::String)

#     ohlc  = read_table(filename)
#     time_conversion =  Calendar.parse("MMM-yy", ohlc[:,1])

#     within!(ohlc, quote
#         Date = $(time_conversion)
#     end)

# # flipud(df)
# end


text:
using Distributions
using Test

const n_samples = 5_000_001

mu = [1.0, 2.0, 3.0]
C = [4. -2. -1.; -2. 5. -1.; -1. -1. 6.]
h = mu
J = C

for d in [
    Dirichlet(3, 2.0),
    Dirichlet([2.0, 1.0, 3.0]),
    IsoNormal(mu, 2.0),
    DiagNormal(mu, [1.5, 2.0, 2.5]),
    MvNormal(mu, C),
    IsoNormalCanon(h, 2.0),
    DiagNormalCanon(h, [1.5, 2.0, 1.2]),
    MvNormalCanon(h, J)]

    println(d)
    dmean = mean(d)
    dcov = cov(d)
    dent = entropy(d)

    x = rand(d, n_samples)
    xmean = vec(mean(x, 2))
    z = x .- xmean
    xcov = (z * z') * (1 / n_samples)

    lp = logpdf(d, x)
    xent = -mean(lp)

    println("expected mean  = $dmean")
    println("empirical mean = $xmean")
    println("--> abs.dev = $(maximum(abs(dmean - xmean)))")

    println("expected cov  = $dcov")
    println("empirical cov = $xcov")
    println("--> abs.dev = $(maximum(abs(dcov - xcov)))")

    println("expected entropy = $dent")
    println("empirical entropy = $xent")
    println("--> abs.dev = $(abs(dent - xent))")

    println()
end


text:
# =============================
# Structs to define a sequence
# =============================

abstract type AbstractSequence end

# "result" of the algorithm
struct Sequence{ST<:AbstractStop} <: AbstractSequence
    hex_id::String
    stops::Vector{ST}
end

name(S::Sequence) = S.hex_id
stop_names(S::Sequence) = name.(S.stops)
stops(S::Sequence) = S.stops

function Base.:∈(stop_id::String, S::Sequence)
    return stop_id ∈ stop_names(S)
end

function Base.push!(S::Sequence, stop)
    push!(S.stops, stop)
end

function Base.length(S::Sequence)
    length(S.stops)
end

# ===================================
# Sequencing algorithms interface
# ===================================

abstract type AbstractSequencingAlgorithm end

# sort a route given by sequence
# we assume that the stops in R and in Q are equal if they
# have the same name
function apply_sequence!(R::AbstractRoute, Q::Sequence)
    S = R.stops
    aux = similar(S)

    for (i, qi) in enumerate(Q.stops)
        idx = findfirst(x -> name(x) == name(qi), S)
        aux[i] = S[idx]
    end
    S .= aux
    return S
end

function _min_zone(R::AbstractRoute; initial_zone=zone(R.stops[2]), tol=1e-4)

    # split stops by zone
    points = point_cloud(R)
    zone_names = zones(R, remove_station=true)
    num_zones = length(zone_names)

    # overapproximate with convex polygons
    P = [overapproximate(ConvexHullArray(pi), HPolygon, 1e-3) for pi in points]

    # preallocate output array
    out = Vector{String}() # FIXME use Vector{Zone} ??
    push!(out, initial_zone)

    # compute cost matrix
    M = Matrix{Float64}(undef, num_zones, num_zones)
    for i in 1:num_zones
        for j in 1:num_zones
            M[i, j] = hausdorff_distance(P[i], P[j], ε=tol)
        end
    end
    T = NamedArray(M, (zone_names, zone_names))

    # for each zone, find the next one that we didn't visit and that is closer
    x = initial_zone    #  convert(String, R.stops[2].zone_id)

    # iterate over zones
    while length(out) < num_zones
        s = sort(view(T, x, :))

        i = 2 # ignoro la posicion actual (tiempo = 0)
        y = s.dicts[1].keys[i]
        while (y ∈ out) # iterar hasta encontrar uno que no haya visitado
            i += 1
            i > length(s.dicts[1]) && return out
            y = s.dicts[1].keys[i]
        end
        push!(out, y)
        x = y
    end

    return out
end


text:
"DSL for Financial Contracts."
module FinancialDSL

import BusinessDays
import InterestRates
import ForwardDiff
import OptimizingIR
import Distributions
import Roots

const CAL_BRL = BusinessDays.Brazil()
BusinessDays.initcache(CAL_BRL)

include("BS.jl")
include("currencies/currencies.jl")
include("marketdata.jl")
include("core/core.jl")

end # module


text:
struct nqueenReward <: SeaPearl.AbstractReward end

function SeaPearl.set_reward!(::Type{SeaPearl.StepPhase}, lh::SeaPearl.LearnedHeuristic{SR,nqueenReward,O}, model::SeaPearl.CPModel, symbol::Union{Nothing,Symbol}) where {SR <: SeaPearl.AbstractStateRepresentation,O <: SeaPearl.ActionOutput}
    lh.current_reward -= 1
    nothing
end

function SeaPearl.set_reward!(::Type{SeaPearl.DecisionPhase}, lh::SeaPearl.LearnedHeuristic{SR,nqueenReward,O}, model::SeaPearl.CPModel) where {SR <: SeaPearl.AbstractStateRepresentation,O <: SeaPearl.ActionOutput}
    lh.current_reward -= 1
    nothing
end

function SeaPearl.set_reward!(::Type{SeaPearl.EndingPhase}, lh::SeaPearl.LearnedHeuristic{SR,nqueenReward,O}, model::SeaPearl.CPModel, symbol::Union{Nothing,Symbol}) where {SR <: SeaPearl.AbstractStateRepresentation,O <: SeaPearl.ActionOutput}
    lh.current_reward += 100 / model.statistics.numberOfNodes
    nothing
end

max_stars_repo_name: corail-research/SeaPearlZoo

max_issues_repo_name: corail-research/SeaPearlZoo

max_forks_repo_name: corail-research/SeaPearlZoo


text:
using Test, Random

println("Not testing solution to Exercise 2")

#=
@testset "Testing solution to Exercise 2" begin
@testset "Running ex2.jl" begin
   include("../ex2.jl")
end;

@testset "Testing that variables exist" begin
   @test @isdefined(response_2a)
   @test @isdefined(response_2b)
   @test @isdefined(response_2c)
   @test @isdefined(response_2d)
   @test @isdefined(response_2e)
   @test @isdefined(response_2f)
   @test @isdefined(response_2g)
   @test @isdefined(response_2h)
   @test @isdefined(response_2i)
   @test @isdefined(response_2j)
   @test @isdefined(response_2k)
   @test @isdefined(response_2l)
   @test @isdefined(response_2m)
end;

@testset "Testing that variables are not missing" begin
   @test !ismissing(response_2a)
   @test !ismissing(response_2b)
   @test !ismissing(response_2c)
   @test !ismissing(response_2d)
   @test !ismissing(response_2e)
   @test !ismissing(response_2f)
   @test !ismissing(response_2g)
   @test !ismissing(response_2h)
   @test !ismissing(response_2i)
   @test !ismissing(response_2j)
   @test !ismissing(response_2k)
   @test !ismissing(response_2l)
   @test !ismissing(response_2m)
end;

end; # Exercise 2
=#


text:
@testset "Test conversion from JSONModel to StandardModel" begin

    model_path = download_data_file(
        "http://bigg.ucsd.edu/static/models/e_coli_core.json",
        joinpath("data", "e_coli_core.json"),
        "7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
    )

    jsonmodel = load_model(model_path)
    stdmodel = convert(StandardModel, jsonmodel)

    # test if same reaction ids
    @test issetequal(reactions(jsonmodel), reactions(stdmodel))
    @test issetequal(metabolites(jsonmodel), metabolites(stdmodel))
    @test issetequal(genes(jsonmodel), genes(stdmodel))
    # not the best tests since it is possible that error could cancel each other out:
    @test sum(stoichiometry(jsonmodel)) == sum(stoichiometry(stdmodel))
    jlbs, jubs = bounds(jsonmodel)
    slbs, subs = bounds(jsonmodel)
    @test sum(jlbs) == sum(slbs)
    @test sum(jubs) == sum(subs)
end


text:
Base.@kwdef mutable struct Benchmark
    b::Vector{Float64} # benchmark weights
    c::Vector{Float64} # Climate costs per $ of revenue
    σ::Vector{Float64} # Volatility
    ρ::Matrix # Correlation matrix
    Σ::Union{Matrix,Nothing} = nothing # Covariance Matrix
end

ρ = [
    1 0.8 0.7 0.6 0.7 0.5 0.7 0.6
    0.8 1 0.75 0.65 0.5 0.6 0.5 0.65
    0.7 0.75 1 0.8 0.7 0.7 0.7 0.7
    0.6 0.65 0.8 1 0.85 0.8 0.75 0.75
    0.7 0.5 0.7 0.85 1 0.6 0.8 0.65
    0.5 0.6 0.7 0.8 0.6 1 0.5 0.7
    0.7 0.5 0.7 0.75 0.8 0.5 1 0.8
    0.6 0.65 0.7 0.75 0.65 0.7 0.8 1
]

using LinearAlgebra
function get_cov_from_corr!(bench::Benchmark)::Benchmark
    D = diagm(bench.σ) # diagonal matrix of volatiltiy
    bench.Σ = D * bench.ρ * D # multiplying columns and rows of correlation matrix by volatiltiy
    return bench
end

our_benchmark = Benchmark(b = [0.23, 0.19, 0.17, 0.13, 0.09, 0.08, 0.06, 0.05],
                        c= [-1.2, 0.8, 2.75, 1.60, -2.75, -1.3, 0.90, -1.70],
                        σ = [0.22, 0.20, 0.25, 0.18, 0.35, 0.23, 0.13, 0.29],
                        ρ = ρ)

get_cov_from_corr!(our_benchmark)

function get_benchmark_score(bench::Benchmark)::Float64
    return round(bench.b' * bench.c; digits= 2)
end

get_benchmark_score(our_benchmark)

Base.@kwdef mutable struct TiltedPortfolio
    x::Vector{Float64} # Optimal weights
end

# portfolio = TiltedPortfolio(x = [0.1, 0.1, 0.3, 0.3, 0.2])

function get_excess_score(portfolio::TiltedPortfolio, bench::Benchmark)::Float64
    return round((portfolio.x - bench.b)' * bench.c;digits = 5)
end

function tracking_error_volatility(portfolio::TiltedPortfolio, bench::Benchmark)::Float64
    return sqrt((portfolio.x - bench.b)' * bench.Σ * (portfolio.x - bench.b))
end

# get_excess_score(portfolio, our_benchmark)

# tracking_error_volatility(portfolio, our_benchmark)

using JuMP, COSMO # for the optimization

# tilt the benchmark according to gamma
function portfolio_tilting(bench::Benchmark, γ::Float64)::TiltedPortfolio
    n = length(bench.c) # number of assets
    model = JuMP.Model(COSMO.Optimizer)
    # the optimal weights we want to find
    @variable(model, x[1:n])
    @objective(model, Min, 1/2 * x' * bench.Σ * x - x' * (γ * bench.c + bench.Σ * bench.b))
    @constraint(model, zeros(n) .<= x .<= ones(n))
    @constraint(model, ones(n)' * x == 1)
    JuMP.optimize!(model)
    portfolio = TiltedPortfolio(x = JuMP.value.(x))
    return portfolio
end

# simulate for gammas between 0 and 10
using Plots
gammas = [i for i in 0.0:0.001:0.012]
excess_scores = zeros(length(gammas))
excess_te = zeros(length(gammas))

for i in 1:length(gammas)
    new_portfolio = portfolio_tilting(our_benchmark, gammas[i])
    excess_scores[i] = get_excess_score(new_portfolio, our_benchmark)
    excess_te[i] = tracking_error_volatility(new_portfolio, our_benchmark)
end

plot(excess_te, excess_scores, xlabel = "Tracking error volatility",
                                ylabel = "Excess ESG score",
                                label = "")

max_stars_repo_name: tlorans/SustainableFinance

max_issues_repo_name: tlorans/SustainableFinance

max_forks_repo_name: tlorans/SustainableFinance


text:
pickN = Function[]

""" `bracket_N(ϕ;  kwargs...)`
   - ϕ is a linemodel
   - kwargs from `bracket`
      - minimization in the interval [a,b]
      - solution to satisfy   α ≤ ϕ'(t\\*) ≤ β and ϕ(t\\*)<ϕ(a)
      - stp : stopping object; atol and rtol are bypassed and the conditions above apply

driver to call `bracket` with `pick_inN`, the plain Newton pick in scheme
"""
function bracket_N(ϕ;  kwargs...)
    reset!(ϕ)
    return bracket(ϕ,  pick_in = pick_inN; kwargs...)
end
push!(pickN, bracket_N)

""" `bracket_N2(ϕ, a, b, α, β, stp;  kwargs...)`
   - ϕ is a linemodel
   - kwargs from `bracket`
      - minimization in the interval [a,b]
      - solution to satisfy   α ≤ ϕ'(t\\*) ≤ β and ϕ(t\\*)<ϕ(a)
      - stp : stopping object; atol and rtol are bypassed and the conditions above apply

driver to call `bracket` with `pick_inN2`, the variant Newton pick in scheme
"""
function bracket_N2(ϕ;  kwargs...)
    reset!(ϕ)
    return bracket(ϕ,  pick_in = pick_inN2; kwargs...)
end
push!(pickN, bracket_N2)

""" bracket_N3(h, a, b, α, β, stp;  kwargs...)
   - ϕ is a linemodel
   - kwargs from `bracket`
      - minimization in the interval [a,b]
      - solution to satisfy   α ≤ ϕ'(t\\*) ≤ β and ϕ(t\\*)<ϕ(a)
      - stp : stopping object; atol and rtol are bypassed and the conditions above apply

driver to call `bracket` with `pick_inN3`, the accelerated Newton pick in scheme
"""
function bracket_N3(ϕ;  kwargs...)
    reset!(ϕ)
    resetN3!()
    return bracket(ϕ,  pick_in = pick_inN3; kwargs...)
end
push!(pickN, bracket_N3)

picks = Function[]

""" driver to call bracket with pick_ins, the plain secant pick in scheme """
function bracket_s(ϕ;  kwargs...)
    reset!(ϕ)
    resetSec!()
    return bracket(ϕ,  pick_in = pick_ins; kwargs...)
end
push!(picks, bracket_s)

""" driver to call bracket with pick_ins2, the accelerated secant pick in scheme
using 2 secant iterations on the cubic interpolant
"""
function bracket_s2(ϕ;  kwargs...)
    reset!(h)
    resetSec!()
    return bracket(ϕ,  pick_in = pick_ins2; kwargs...)
end
push!(picks, bracket_s2)

""" driver to call bracket with pick_ins2N, the accelerated secant pick in scheme
using 2 Newton iterations on the cubic interpolant
"""
function bracket_s2N(ϕ;  kwargs...)
    reset!(ϕ)
    resetSec!()
    return bracket(ϕ,  pick_in = pick_ins2N; kwargs...)
end
push!(picks, bracket_s2N)

""" driver to call bracket with pick_ins2Opt, the accelerated secant pick in scheme
using 2 Newton iterations on the cubic interpolant (variant)
"""
function bracket_s2Opt(ϕ;  kwargs...)
    reset!(ϕ)
    resetSec!()
    return bracket(ϕ,  pick_in = pick_ins2Opt; kwargs...)
end
push!(picks, bracket_s2Opt)

""" driver to call bracket with pick_insB, the basic bissection pick in scheme
"""

function bracket_B(ϕ;  kwargs...)
    reset!(ϕ)
    return bracket(ϕ,  pick_in = pick_inB; kwargs...)
end

export pickN, picks

for n in nameof.(pickN ∪ picks )
    @eval export $n
end

export bracket_B


text:
# # 120: Differing species sets in regions, 1D
# ([source code](SOURCE_URL))

module Example120_ThreeRegions1D

using Printf
using VoronoiFVM
using ExtendableGrids
using GridVisualize

function main(;n=30,Plotter=nothing,plot_grid=false, verbose=false,unknown_storage=:sparse,tend=10)
    h=3.0/(n-1)
    X=collect(0:h:3.0)
    grid=VoronoiFVM.Grid(X)
    cellmask!(grid,[0.0],[1.0],1)
    cellmask!(grid,[1.0],[2.1],2)
    cellmask!(grid,[1.9],[3.0],3)

    subgrid1=subgrid(grid,[1])
    subgrid2=subgrid(grid,[1,2,3])
    subgrid3=subgrid(grid,[3])

    if plot_grid
        plotgrid(grid,Plotter=Plotter)
        return
    end

    eps=[1,1,1]
    k=[1,1,1]

    physics=VoronoiFVM.Physics(
    num_species=3,
    reaction=function(f,u,node)
        if node.region==1
            f[1]=k[1]*u[1]
            f[2]=-k[1]*u[1]
        elseif node.region==3
            f[2]=k[3]*u[2]
            f[3]=-k[3]*u[2]
        else
            f[1]=0
        end
    end,

        flux=function(f,u0,edge)
            u=unknowns(edge,u0)

        if edge.region==1
            f[1]=eps[1]*(u[1,1]-u[1,2])
            f[2]=eps[2]*(u[2,1]-u[2,2])
        elseif edge.region==2
            f[2]=eps[2]*(u[2,1]-u[2,2])
        elseif edge.region==3
            f[2]=eps[2]*(u[2,1]-u[2,2])
            f[3]=eps[3]*(u[3,1]-u[3,2])
        end
    end,

    source=function(f,node)
        if node.region==1
           f[1]=1.0e-4*(3.0-node[1])
        end
    end,

    storage=function(f,u,node)
        f.=u
    end
    )

    sys=VoronoiFVM.System(grid,physics,unknown_storage=unknown_storage)

    enable_species!(sys,1,[1])
    enable_species!(sys,2,[1,2,3])
    enable_species!(sys,3,[3])

    boundary_dirichlet!(sys,3,2,0.0)

    inival=unknowns(sys)
    U=unknowns(sys)
    inival.=0

    control=VoronoiFVM.NewtonControl()
    control.verbose=verbose
    tstep=0.01
    time=0.0
    istep=0
    testval=0
    p=GridVisualizer(Plotter=Plotter,layout=(1,1))
    while time<tend
        time=time+tstep
        solve!(U,inival,sys,control=control,tstep=tstep)
        inival.=U
        if verbose
            @printf("time=%g\n",time)
        end
        tstep*=1.1
        istep=istep+1
        testval=U[2,5]

        U1=view(U[1,:],subgrid1)
        U2=view(U[2,:],subgrid2)
        U3=view(U[3,:],subgrid3)

        scalarplot!(p[1,1],subgrid1, U1,label="spec1", color=(0.5,0,0),xlimits=(0,3),flimits=(0,1e-3),
                  title=@sprintf("three regions t=%.3g",time))
        scalarplot!(p[1,1],subgrid2, U2,label="spec2", color=(0.0,0.5,0),clear=false)
        scalarplot!(p[1,1],subgrid3, U3,label="spec3", color=(0.0,0.0,0.5),clear=false,show=true)
        sleep(1.0e-2)
    end
    return testval
end

function test()
    testval=0.0005954993329548969
    main(unknown_storage=:sparse) ≈ testval &&
        main(unknown_storage=:dense) ≈ testval
end

end


text:
"""
    struct GeneralSpace{𝕜} <: ElementarySpace{𝕜}

A finite-dimensional space over an arbitrary field `𝕜` without additional structure. It is
thus characterized by its dimension, and whether or not it is the dual and/or conjugate
space. For a real field `𝕜`, the space and its conjugate are the same.
"""
struct GeneralSpace{𝕜} <: ElementarySpace{𝕜}
    d::Int
    dual::Bool
    conj::Bool
    function GeneralSpace{𝕜}(d::Int, dual::Bool, conj::Bool) where {𝕜}
        d >= 0 ||
            throw(ArgumentError("Dimension of a vector space should be bigger than zero"))
        if 𝕜 isa Field
            new{𝕜}(Int(d), dual, (𝕜 ⊆ ℝ) ? false : conj)
        else
            throw(ArgumentError("Unrecognised scalar field: $𝕜"))
        end
    end
end
GeneralSpace{𝕜}(d::Int; dual::Bool = false, conj::Bool = false) where {𝕜} =
    GeneralSpace{𝕜}(d, dual, conj)

dim(V::GeneralSpace) = V.d
isdual(V::GeneralSpace) = V.dual
isconj(V::GeneralSpace) = V.conj

Base.axes(V::GeneralSpace) = Base.OneTo(dim(V))

dual(V::GeneralSpace{𝕜}) where {𝕜} =
    GeneralSpace{𝕜}(dim(V), !isdual(V), isconj(V))
Base.conj(V::GeneralSpace{𝕜}) where {𝕜} =
    GeneralSpace{𝕜}(dim(V), isdual(V), !isconj(V))

function Base.show(io::IO, V::GeneralSpace{𝕜}) where {𝕜}
    if isconj(V)
        print(io, "conj(")
    end
    print(io, "GeneralSpace{", 𝕜, "}(", dim(V), ")")
    if isdual(V)
        print(io, "'")
    end
    if isconj(V)
        print(io, ")")
    end
end
Base.show(io::IO, ::Type{GeneralSpace}) = print(io, "GeneralSpace")


text:
using DReal

import DReal: init_dreal!, Context, Var, add!, set_precision!, is_satisfiable, model

begin
# Blank slate
reset_ctx!()

# ## Example
# ## =======
init_dreal!()
ctx = Context(qf_nra)
set_precision!(ctx, 0.00001)
x = Var(ctx, Float64, "x", -3.141592, 3.141592)
y = Var(ctx, Float64, "y", -3.141592, 3.141592)
add!(sin(x) == cos(y))
res = is_satisfiable(ctx)
if res
  println("x = $(model(ctx,x))")
  println("y = $(model(ctx,y))")
end
end


text:
"""
In-Memory implementation of Kernel Spectral Clustering
"""
function kscnet()
  # TODO:
  return "KSCnet.jl not yet implemented."
end


text:
@testset "Conversion from and to SBML model" begin
    json_model = download_data_file(
        "http://bigg.ucsd.edu/static/models/iJO1366.json",
        joinpath("data", "iJO1366.json"),
        "9376a93f62ad430719f23e612154dd94c67e0d7c9545ed9d17a4d0c347672313",
    )

    jm = load_json_model(json_model)
    sm = convert(StandardModel, jm)
    jm2 = convert(JSONModel, sm)

    @test Set(reactions(jm)) == Set(reactions(sm))
    @test Set(reactions(jm)) == Set(reactions(jm2))
end


text:
"""
$SIGNATURES

Initialize an empty in-memory DistributedFactorGraph `::DistributedFactorGraph` object.
"""
function initfg(dfg::T=InMemDFGType(solverParams=SolverParams());
                                sessionname="NA",
                                robotname="",
                                username="",
                                cloudgraph=nothing) where T <: AbstractDFG
#
  #
  return dfg
end

#init an empty fg with a provided type and SolverParams
function initfg(::Type{T};solverParams=SolverParams(),
                      sessionname="NA",
                      robotname="",
                      username="",
                      cloudgraph=nothing) where T <: AbstractDFG
  #
  return T(solverParams=solverParams)
end

function initfg(::Type{T},solverParams::SolverParams;
                      sessionname="NA",
                      robotname="",
                      username="",
                      cloudgraph=nothing) where T <: AbstractDFG
  #
  return T{SolverParams}(solverParams=solverParams)
end

# Should deprecate in favor of TensorCast.jl
reshapeVec2Mat(vec::Vector, rows::Int) = reshape(vec, rows, round(Int,length(vec)/rows))

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

"""
    $(SIGNATURES)

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

"""
    $(SIGNATURES)

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

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

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

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

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

"""
    $SIGNATURES

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

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

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

function setValKDE!(vd::VariableNodeData,
                    val::Array{Float64,2},
                    setinit::Bool=true,
                    inferdim::Float64=0.0)::Nothing
  # recover variableType information
  sty = getVariableType(vd)
  p = AMP.manikde!(val, getManifolds(sty))
  setValKDE!(vd, p, setinit, inferdim)
  nothing
end

function setValKDE!(v::DFGVariable,
                    val::Array{Float64,2},
                    bws::Array{Float64,2},
                    setinit::Bool=true,
                    inferdim::Float64=0;
                    solveKey::Symbol=:default)::Nothing
  # recover variableType information
  setValKDE!(getSolverData(v, solveKey), val, bws[:,1], setinit, inferdim )

  nothing
end

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

"""
    $SIGNATURES

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

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

"""
    $SIGNATURES

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

"""
    $SIGNATURES

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

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

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

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

function DefaultNodeDataParametric( dodims::Int,
                                    dims::Int,
                                    variableType::InferenceVariable;
                                    initialized::Bool=true,
                                    dontmargin::Bool=false)::VariableNodeData

  # this should be the only function allocating memory for the node points
  if false && initialized
    error("not implemented yet")
    # pN = AMP.manikde!(randn(dims, N), variableType.manifolds);
    #
    # sp = Int[0;] #round.(Int,range(dodims,stop=dodims+dims-1,length=dims))
    # gbw = getBW(pN)[:,1]
    # gbw2 = Array{Float64}(undef, length(gbw),1)
    # gbw2[:,1] = gbw[:]
    # pNpts = getPoints(pN)
    # #initval, stdev
    # return VariableNodeData(pNpts,
    #                         gbw2, Symbol[], sp,
    #                         dims, false, :_null, Symbol[], variableType, true, 0.0, false, dontmargin)
  else
    sp = round.(Int,range(dodims,stop=dodims+dims-1,length=dims))
    return VariableNodeData(zeros(dims, 1),
                            zeros(dims,1), Symbol[], sp,
                            dims, false, :_null, Symbol[], variableType, false, 0.0, false, dontmargin, 0, 0, :parametric)
  end

end

function setDefaultNodeDataParametric!(v::DFGVariable, variableType::InferenceVariable; kwargs...)
  vnd = DefaultNodeDataParametric(0, variableType |> getDimension, variableType; kwargs...)
  setSolverData!(v, vnd, :parametric)
  return nothing
end

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

      pN = AMP.manikde!(randn(dims, N), getManifolds(varType));

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

"""
    $SIGNATURES

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

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

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

# get instance from variableType
_variableType(varType::InferenceVariable) = varType
_variableType(varType::Type{<:InferenceVariable}) = varType()

"""
$(SIGNATURES)

Add a variable node `label::Symbol` to `dfg::AbstractDFG`, as `varType<:InferenceVariable`.

Notes
-----
- keyword `nanosecondtime` is experimental and intended as the whole subsection portion -- i.e. accurateTime = (timestamp MOD second) + Nanosecond

Example
-------

```julia
fg = initfg()
addVariable!(fg, :x0, Pose2)
```
"""
function addVariable!(dfg::AbstractDFG,
                      label::Symbol,
                      varTypeU::Union{T, Type{T}};
                      N::Int=getSolverParams(dfg).N,
                      solvable::Int=1,
                      timestamp::Union{DateTime,ZonedDateTime}=now(localzone()),
                      nanosecondtime::Union{Nanosecond,Int64,Nothing}=nothing,
                      dontmargin::Bool=false,
                      labels::Union{Vector{Symbol},Nothing}=nothing,
                      tags::Vector{Symbol}=Symbol[],
                      smalldata=Dict{Symbol, DFG.SmallDataTypes}(),
                      checkduplicates::Bool=true,
                      initsolvekeys::Vector{Symbol}=getSolverParams(dfg).algorithms ) where T<:InferenceVariable
  #
  varType = _variableType(varTypeU)
  # TODO Remove deprecation in v0.16
  if :ut in fieldnames(T)
    Base.depwarn("Field `ut` (microseconds) for variable type ($T) has been deprecated please use DFGVariable.nstime, kwarg: nanosecondtime", :addVariable!)
    if isnothing(nanosecondtime)
      varType.ut == -9999999999 && error("please define a time for type $(T), use FGVariable.nstime, kwarg: nanosecondtime")
      nanosecondtime = Nanosecond(varType.ut*1000)
    else
      @warn "Nanosecond time has been specified as $nanosecondtime, ignoring `ut` field value: $(varType.ut)."
    end
  elseif isnothing(nanosecondtime)
    nanosecondtime = Nanosecond(0)
  end

  # deprecate in v0.16
  labels isa Vector ? (union!(tags, labels); @warn("labels is deprecated, use tags instead")) : nothing

  union!(tags, [:VARIABLE])
  v = DFGVariable(label, varType; tags=Set(tags), smallData=smalldata, solvable=solvable, timestamp=timestamp, nstime=Nanosecond(nanosecondtime))

  (:default in initsolvekeys) &&
    setDefaultNodeData!(v, 0, N, getDimension(varType), initialized=false, varType=varType, dontmargin=dontmargin) # dodims

  (:parametric in initsolvekeys) &&
    setDefaultNodeDataParametric!(v, varType, initialized=false, dontmargin=dontmargin)

  DFG.addVariable!(dfg, v)

  return v
end

"""
    $(SIGNATURES)

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

"""
    $(SIGNATURES)

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

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

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

  # get solvefor manifolds
  manis = length(Xi)==0 || sfidx==0 ? (:null,) : getManifolds(Xi[sfidx])

  # FIXME, forcing maxlen to N results in errors (see test/testVariousNSolveSize.jl) see #105
  # maxlen = N == 0 ? maxlen : N
  return maxlen, sfidx, manis
end

function parseusermultihypo(multihypo::Nothing, nullhypo::Float64)
  verts = Symbol[]
  mh = nothing
  return mh, nullhypo
end
function parseusermultihypo(multihypo::Vector{Float64}, nullhypo::Float64)
  mh = nothing
  if 0 < length(multihypo)
    multihypo2 = multihypo
    multihypo2[1-1e-10 .< multihypo] .= 0.0
    # check that terms sum to full probability
    @assert abs(sum(multihypo2) % 1) < 1e-10  || 1-1e-10 < sum(multihypo2) % 1 "ensure multihypo sums to a (or nearly, 1e-10) interger, see #1086"
    # check that only one variable broken into fractions
    @assert sum(multihypo2[1e-10 .< multihypo2]) ≈ 1
    # force normalize something that is now known to be close
    multihypo2 ./= sum(multihypo2)

    mh = Categorical(Float64[multihypo2...] )
  end
  return mh, nullhypo
end

# import IncrementalInference: prepgenericconvolution, convert

"""
    $SIGNATURES

Function to calculate measurement dimension from factor sampling.

Notes
- Will not work in all situations, but good enough so far.
  - # TODO standardize via domain or manifold definition...??
"""
function calcZDim(cf::CalcFactor{T}) where {T <: FunctorInferenceType}
  #
  # zdim = T != GenericMarginal ? size(getSample(usrfnc, 2)[1],1) : 0
  zdim = if T != GenericMarginal
    # vnds = Xi # (x->getSolverData(x)).(Xi)
    # NOTE try to make sure we get matrix back (not a vector)
    smpls = sampleFactor(cf, 2)[1]
    size(smpls,1)
  else
    0
  end
  return zdim
end

function prepgenericconvolution(Xi::Vector{<:DFGVariable},
                                usrfnc::T;
                                multihypo::Union{Nothing, Distributions.Categorical}=nothing,
                                nullhypo::Real=0.0,
                                threadmodel=MultiThreaded  ) where {T <: FunctorInferenceType}
  #
  ARR = Array{Array{Float64,2},1}()
  maxlen, sfidx, manis = prepareparamsarray!(ARR, Xi, nothing, 0) # Nothing for init.
  fldnms = fieldnames(T) # typeof(usrfnc)

  # standard factor metadata
  sflbl = 0==length(Xi) ? :null : getLabel(Xi[end])
  fmd = FactorMetadata(Xi, getLabel.(Xi), ARR, sflbl, nothing)
  cf = CalcFactor( usrfnc, fmd, 0, 1, (Matrix{Float64}(undef,0,0),), ARR)

  zdim = calcZDim(cf)
  # zdim = T != GenericMarginal ? size(getSample(usrfnc, 2)[1],1) : 0
  certainhypo = multihypo !== nothing ? collect(1:length(multihypo.p))[multihypo.p .== 0.0] : collect(1:length(Xi))

  ccw = CommonConvWrapper(
          usrfnc,
          zeros(1,0),
          zdim,
          ARR,
          fmd,
          specialzDim = sum(fldnms .== :zDim) >= 1,
          partial = sum(fldnms .== :partial) >= 1,
          hypotheses=multihypo,
          certainhypo=certainhypo,
          nullhypo=nullhypo,
          threadmodel=threadmodel
        )
  #
  return ccw
end

# TODO perhaps consolidate with constructor?
"""
$SIGNATURES

Generate the default factor data for a new DFGFactor.
"""
function getDefaultFactorData(
      dfg::AbstractDFG,
      Xi::Vector{<:DFGVariable},
      usrfnc::T;
      multihypo::Vector{<:Real}=Float64[],
      nullhypo::Float64=0.0,
      threadmodel=SingleThreaded,
      eliminated::Bool = false,
      potentialused::Bool = false,
      edgeIDs = Int[],
      solveInProgress = 0 ) where T <: FunctorInferenceType
  #

  # prepare multihypo particulars
  # storeMH::Vector{Float64} = multihypo == nothing ? Float64[] : [multihypo...]
  mhcat, nh = parseusermultihypo(multihypo, nullhypo)

  # allocate temporary state for convolutional operations (not stored)
  ccw = prepgenericconvolution(Xi, usrfnc, multihypo=mhcat, nullhypo=nh, threadmodel=threadmodel)

  # and the factor data itself
  return FunctionNodeData{CommonConvWrapper{T}}(eliminated, potentialused, edgeIDs, ccw, multihypo, ccw.certainhypo, nullhypo, solveInProgress)

end

"""
    $SIGNATURES

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

"""
    $SIGNATURES

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

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

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

Related

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

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

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

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

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

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

"""
    $(SIGNATURES)

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

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

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

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

  didinit = true

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

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

"""
    $(TYPEDSIGNATURES)

Method to manually initialize a variable using a set of points.

Notes
- Disable automated graphinit on `addFactor!(fg, ...; graphinit=false)
  - any un-initialized variables will automatically be initialized by `solveTree!`

Example:

```julia
# some variable is added to fg
addVariable!(fg, :somepoint3, ContinuousEuclid{2})

# data is organized as (row,col) == (dimension, samples)
pts = randn(2,100)
initManual!(fg, :somepoint3, pts)

# manifold management should be done automatically.
# note upgrades are coming to consolidate with Manifolds.jl, see RoME #244

## it is also possible to initManual! by using existing factors, e.g.
initManual!(fg, :x3, [:x2x3f1])
```

DevNotes
- TODO better document graphinit and treeinit.
"""
function initManual!( variable::DFGVariable,
                      ptsArr::BallTreeDensity)
  #
  setValKDE!(variable, ptsArr, true)
  return nothing
end
function initManual!( dfg::AbstractDFG,
                      label::Symbol,
                      belief::BallTreeDensity)
  #
  variable = getVariable(dfg, label)
  initManual!(variable, belief)
  return nothing
end
function initManual!( dfg::AbstractDFG,
                      label::Symbol,
                      usefcts::Vector{Symbol} )
  #
  @info "initManual! $label"
  pts = predictbelief(dfg, label, usefcts)[1]
  vert = getVariable(dfg, label)
  Xpre = AMP.manikde!(pts, getVariableType(vert) |> getManifolds )
  setValKDE!(vert, Xpre, true)
  return nothing
end

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

const initVariableManual! = initManual!

"""
    $SIGNATURES

Set solveKey values of `dest::AbstractDFG` according to `initKey::Symbol=:graphinit` values.

Notes
- Some flexibility for using two DFGs and different key values, see Examples and code for details.
- Can also be specific with `varList::Vector{Symbol}`.
- Returns `dest` graph.
- Uses the supersolve mechanism.

Examples
```julia
resetInitialValues!(fg)
resetInitialValues!(fg1,fg2)  # into 1 from 2
resetInitialValues!(fg1,fg1,:myotherinit)  # use different init value into solveKey :default
resetInitialValues!(fg1,fg1,:graphinit, :mysolver) # not into solveKey=:default but :mysolver
resetInitialValues!(fg1, varList=[:x1;:l3])  # Specific variables only

# Into `fgNew` object, leaving `fg` untouched
fgNew = deepcopy(fg)
resetInitialValues!(fgNew,fg)
```

Related

initManual!, graphinit (keyword)
"""
function resetInitialValues!(dest::AbstractDFG,
                             src::AbstractDFG=dest,
                             initKey::Symbol=:graphinit,
                             solveKey::Symbol=:default;
                             varList::AbstractVector{Symbol}=ls(dest))
  #
  for vs in varList
    vnd = getSolverData(getVariable(src, vs), initKey)
    # guess we definitely want to use copy to preserve the initKey memory
    updateVariableSolverData!(dest,vs,vnd,solveKey,true; warn_if_absent=false)
  end
  return dest
end
const resetInitValues! = resetInitialValues!

"""
    $SIGNATURES

Ensure that no variables set as `solvable=1` are floating free without any connected `solvable=1` factors.  If any found, then set those 'free' variable's `solvable=solvableFallback` (default `0`).

Related

ensureAllInitialized!
"""
function ensureSolvable!(dfg::AbstractDFG; solvableTarget::Int=1, solvableFallback::Int=0)
  # workaround in case isolated variables occur
  solvVars = ls(dfg, solvable=solvableTarget)
  varHasFact = (x->length(ls(dfg,x, solvable=solvableTarget))==0).(solvVars)
  blankVars = solvVars[findall(varHasFact)]
  if 0 < length(blankVars)
    @warn("solveTree! dissallows solvable variables without any connected solvable factors -- forcing solvable=0 on $(blankVars)")
    (x->setSolvable!(dfg, x, solvableFallback)).(blankVars)
  end
  return blankVars
end

"""
    $SIGNATURES

Perform `graphinit` over all variables with `solvable=1` (default).

Related

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

function assembleFactorName(dfg::AbstractDFG,
                            Xi::Vector{<:DFGVariable} )
  #

  existingFactorLabels = listFactors(dfg)
  existingFactorLabelDict = Dict(existingFactorLabels .=> existingFactorLabels)
  namestring = ""
  for vert in Xi #f.Xi
    namestring = string(namestring,vert.label)
  end
  opt = getSolverParams(dfg)
  for i in 1:opt.maxincidence
    tempnm = string(namestring, "f$i")
    if !haskey(existingFactorLabelDict, Symbol(tempnm))
      namestring = tempnm
      break
    end
    i != opt.maxincidence ? nothing : error("Artificial restriction to not connect more than $(opt.maxincidence) factors to a variable (bad for sparsity), try setting getSolverParams(fg).maxincidence=1000 to adjust this restriction.")
  end
  return Symbol(namestring)
end

"""
    $(SIGNATURES)

Add factor with user defined type <: FunctorInferenceType to the factor graph
object. Define whether the automatic initialization of variables should be
performed.  Use order sensitive `multihypo` keyword argument to define if any
variables are related to data association uncertainty.
"""
function addFactor!(dfg::AbstractDFG,
                    Xi::Vector{<:DFGVariable},
                    usrfnc::R;
                    multihypo::Vector{Float64}=Float64[],
                    nullhypo::Float64=0.0,
                    solvable::Int=1,
                    tags::Vector{Symbol}=Symbol[],
                    timestamp::Union{DateTime,ZonedDateTime}=now(localzone()),
                    graphinit::Bool=getSolverParams(dfg).graphinit,
                    threadmodel=SingleThreaded,
                    suppressChecks::Bool=false  ) where
                      {R <: FunctorInferenceType}
  #
  # depcrecation

  varOrderLabels = [v.label for v=Xi]
  namestring = assembleFactorName(dfg, Xi)
  solverData = getDefaultFactorData(dfg, Xi, deepcopy(usrfnc), multihypo=multihypo, nullhypo=nullhypo, threadmodel=threadmodel)
  newFactor = DFGFactor(Symbol(namestring),
                        varOrderLabels,
                        solverData;
                        tags=Set(union(tags, [:FACTOR])),
                        solvable=solvable,
                        timestamp=timestamp)
  #

  success = DFG.addFactor!(dfg, newFactor)

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

  return newFactor
end

function addFactor!(dfg::AbstractDFG,
                    xisyms::Vector{Symbol},
                    usrfnc::FunctorInferenceType;
                    multihypo::Vector{<:Real}=Float64[],
                    nullhypo::Float64=0.0,
                    solvable::Int=1,
                    timestamp::Union{DateTime,ZonedDateTime}=now(localzone()),
                    tags::Vector{Symbol}=Symbol[],
                    graphinit::Bool=getSolverParams(dfg).graphinit,
                    threadmodel=SingleThreaded,
                    suppressChecks::Bool=false  )
  #
  # depcrecation

  # basic sanity check for unary vs n-ary
  if !suppressChecks && length(xisyms) == 1 && !(usrfnc isa AbstractPrior) && !(usrfnc isa Mixture)
    @warn("Listing only one variable $xisyms for non-unary factor type $(typeof(usrfnc))")
  end

  variables = getVariable.(dfg, xisyms)
  # verts = map(vid -> DFG.getVariable(dfg, vid), xisyms)
  addFactor!(dfg, variables, usrfnc, multihypo=multihypo, nullhypo=nullhypo, solvable=solvable, tags=tags, graphinit=graphinit, threadmodel=threadmodel, timestamp=timestamp )
end

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

"""
    $SIGNATURES

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

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

Future
- TODO: `A` should be sparse data structure (when we exceed 10'000 var dims)
- TODO: Incidence matrix is rectagular and adjacency is the square.
"""
function getEliminationOrder( dfg::G;
                              ordering::Symbol=:qr,
                              solvable::Int=1,
                              constraints::Vector{Symbol}=Symbol[]) where G <: AbstractDFG
  #
  @assert 0 == length(constraints) || ordering == :ccolamd "Must use ordering=:ccolamd when trying to use constraints"
  # Get the sparse adjacency matrix, variable, and factor labels
  adjMat, permuteds, permutedsf = DFG.getBiadjacencyMatrix(dfg, solvable=solvable)
  # adjMat, permuteds, permutedsf = DFG.getAdjacencyMatrixSparse(dfg, solvable=solvable)

  # Create dense adjacency matrix

  p = Int[]
  if ordering==:chol
    # hack for dense matrix....
    A = adjMat
    p = cholesky(Matrix(A'A),Val(true)).piv
    @warn "check that cholesky ordering is not reversed -- basically how much fill in (separator size) are you seeing???  Long skinny chains in tree is bad."
  elseif ordering==:qr
    # hack for dense matrix....
    A = Array(adjMat)
    # this is the default
    q,r,p = qr(A, Val(true))
    p .= p |> reverse
  elseif ordering==:ccolamd
    cons = zeros(SuiteSparse_long, length(adjMat.colptr) - 1)
    cons[findall(x->x in constraints, permuteds)] .= 1
    p = Ccolamd.ccolamd(adjMat, cons)
    @warn "Ccolamd is experimental in IIF at this point in time."
  else
    prtslperr("getEliminationOrder -- cannot do the requested ordering $(ordering)")
  end

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

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

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

function addChainRuleMarginal!( dfg::AbstractDFG,
                                Si::Vector{Symbol} )
  #

  lbls = String[]
  genmarg = GenericMarginal()
  Xi = map(v -> DFG.getVariable(dfg, v), Si)
  # @info "adding marginal to"
  # for x in Xi
  #   @info "x.index=",x.index
  # end
  addFactor!( dfg, Xi, genmarg, graphinit=false, suppressChecks=true )
  nothing
end

function rmVarFromMarg( dfg::AbstractDFG,
                        fromvert::DFGVariable,
                        gm::Vector{DFGFactor}  )
  #

  @debug " - Removing $(fromvert.label)"
  for m in gm
    @debug "Looking at $(m.label)"
    for n in DFG.getNeighbors(dfg, m) #x1, x2
      if n == fromvert.label # n.label ==? x1
        @debug "   - Breaking link $(m.label)->$(fromvert.label)..."
        @debug "     - Original links: $(DFG.ls(dfg, m))"
        remvars = setdiff(DFG.ls(dfg, m), [fromvert.label])
        @debug "     - New links: $remvars"

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

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

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

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

      if typeof(_getCCW(fct)) == CommonConvWrapper{GenericMarginal}
        push!(gm, fct)
      end
    end

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

    # mark variable
    getSolverData(vert).eliminated = true

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

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

  end
  return nothing
end

"""
    $(SIGNATURES)

Get KernelDensityEstimate kde estimate stored in variable node.
"""
getBelief(vnd::VariableNodeData) = AMP.manikde!(getVal(vnd), getBW(vnd)[:,1], getVariableType(vnd) |> getManifolds)
getBelief(v::DFGVariable, solvekey::Symbol=:default) = getKDE(getSolverData(v, solvekey))
getBelief(dfg::AbstractDFG, lbl::Symbol, solvekey::Symbol=:default) = getKDE(getVariable(dfg, lbl), solvekey)

const getKDE = getBelief

#


text:
using StaticArrays: StaticArray, SMatrix, SVector
using LinearAlgebra

function seed(v::SVector{N}) where N
    SMatrix{N,N,eltype(v)}(I)
end

function _seed(v, ij)
    i, j = Tuple(ij)
    vi = v[i]
    return (i==j) ? one(vi) : zero(vi)
end

function seed(v)
    vv = vec(v)
    ax = axes(vv, 1)
    return _seed.(Ref(vv), CartesianIndices((ax, ax)))
end

###
### Derivative object:
###

"""
    D(f)

`D(f)(x) * v` computes ``\\frac{df}{dx}(x) ⋅ v``
"""
struct D{T,F}
    f::F
    x::T

    D(f) = new{Nothing,typeof(f)}(f, nothing)
    (dd::D{<:Nothing,F})(x::T) where {T,F} = new{T,F}(dd.f, x)
    D(dd::D) = D(DI(dd.f))
end

"""
    DI(f)

`DI(f)(x)` is a convenience function to materialize the derivative, gradient or
Jacobian of `f` at `x`.
"""
DI(f) = x->D(f)(x) * mul_identity(x)

mul_identity(x::AbstractArray) = I
mul_identity(x) = one(x)

# WARNING: It assume that the number type is commutative
Base.:*(v::Number, dd::D) = dd * v
function Base.:*(dd::D{<:Number}, v::Number)
    ps, derivative_dual = dualrun() do
        dual = Dual(dd.x, v)
        partials(dual), dd.f(dual)
    end
    derivative = map(derivative_dual) do a
        pa = partials(a)
        pa isa Zero ? zero(ps) : pa
    end
    return derivative
end

unwrap_adj(x::Union{Transpose,Adjoint}) = unwrap_adj(parent(x))
unwrap_adj(x) = x

function Base.:*(dd::D{<:AbstractArray}, V::Union{AbstractArray,UniformScaling})
    checkinput(dd.x, V)
    xx_partial = V isa UniformScaling ? seed(dd.x) :
                    V isa AbstractVector ? V :
                    V'
    duals = dualrun() do
        dualarray = DualArray(dd.x, xx_partial)
        dd.f(dualarray)
    end

    J_dual = unwrap_adj(duals)
    J_sz = xx_partial isa AbstractVector ?
               (length(J_dual), ) :
               (length(J_dual), size(xx_partial, 2))

    if J_dual isa AbstractArray # Jacobian
        if J_dual isa DualArray
            ps = allpartials(J_dual)
            return isjvp(J_dual) ? ps : ps'
        elseif J_dual isa StaticArray
            return extract_diffresult(J_dual, J_sz)
        else
            # `f: R^n -> R^m` so the Jacobian is `m × n`
            J = similar(J_dual, valtype(eltype(J_dual)), J_sz)
            extract_diffresult!(J, J_dual)
            return J
        end
    else # gradient
        ps = partials(J_dual)
        if V isa AbstractVector
            return ps isa Zero ? zero(J_dual) : first(ps)
        else
            rhs = V isa UniformScaling ? dd.x : @view(V[1, :])
            return ps isa Zero ? (false * rhs)' : ps'
        end
    end
end

@noinline throwinput() = throw(ArgumentError("In `D(f)(x) * v`, `x` and `v` must have the same axes."))
checkinput(x, v) = (all(i->axes(x, i) == axes(v, i), 1:ndims(x)) || throwinput(); nothing)
checkinput(_, ::UniformScaling) = nothing

function extract_diffresult(xs, mn::Tuple)
    isjvp = mn isa Tuple{<:Integer}
    m, n = isjvp ? (first(mn), 1) : mn
    tup = mapreduce((x,y)->tuple(x..., y...), xs.data) do x
        ps = partials(x)
        if ps isa Zero
            ntuple(_->false, n)
        else
            ps.data
        end
    end
    return isjvp ? SVector{m}(tup) : SMatrix{n, m}(tup)'
end

function extract_diffresult!(J, ds::AbstractArray)
    @inbounds for (d, i) in zip(ds, axes(J, 1))
        if partials(d) isa Zero
            J[i, :] .= false
        else
            J[i, :] .= partials(d)
        end
    end
    return nothing
end

@inline Cassette.overdub(ctx::TaggedCtx, f::typeof(extract_diffresult!), args...) = f(args...)
@inline Cassette.overdub(ctx::TaggedCtx, f::typeof(extract_diffresult), args...) = f(args...)

# pretty printing
function Base.show(io::IO, dd::D{T}) where T
    print(io, "D(", nameof(dd.f), ')')
    if !(T <: Nothing)
        print(io, '(')
        show(IOContext(io, :compact => true, :limit => true), dd.x)
        print(io, ')')
    end
    return nothing
end


text:
using LinearAlgebra

"""
vout = eceftoseof(vin, satpos, satvel, vintype)

ECEFで表されるベクトルVinを、SEOFに座標変換する

# Argments
 - `vin`：座標変換を行う対象となるベクトル @ECEF
 - `satpos`：衛星位置ベクトル @ECEF
 - `satvel`：衛星速度ベクトル @ECEF
 - `vintype`：vinの種類（省略可、省略時vintype="dir"）　vintype="pos"->位置ベクトル vintype="dir"->方向ベクトル

# Return
 - `vout`：座標変換後のベクトル @SEOF

"""

function eceftoseof(vin, satpos, satvel, vintype="dir")
    # SEOFのz軸計算
    seof_z = -1.0 * satpos / norm(satpos)
    # SEOFのx軸計算
    normvel = satvel / norm(satvel)
    cθ = dot(seof_z, normvel)
    seof_x = normvel - cθ*seof_z
    # SEOFのy軸計算
    seof_y = cross(vec(seof_z), vec(seof_x))
    seof_y = seof_y / norm(seof_y)

    # ECEF -> SEOF()の座標変換行列の作成
    M_eceftoseof = seof_x'
    M_eceftoseof = hcat(M_eceftoseof, reshape(seof_y, 3, 1))
    M_eceftoseof = hcat(M_eceftoseof, seof_z')

    # 座標変換(回転)
    vout = M_eceftoseof \ vin'
    if vintype == "pos"
        satpos_seof = M_eceftoseof \ satpos'
        # 座標変換（原点の移動)
        vout = vout' .- satpos_seof'
    elseif vintype == "dir"
        vout = vout'
    else
        throw(ArgumentError("Enter \"dir\" or \"pos\" for vintype"))
    end

    return vout

end

"""
vout = eceftoseof(vin, satpos, satvel, vintype)

ECEFで表されるベクトルVinを、SEOFに座標変換する

# Argments
 - `vin`：座標変換を行う対象となるベクトル @ECEF
 - `satpos`：衛星位置ベクトル @ECEF
 - `satvel`：衛星速度ベクトル @ECEF
 - `vintype`：vinの種類（省略可、省略時vintype="dir"）　vintype="pos"->位置ベクトル vintype="dir"->方向ベクトル

# Return
 - `vout`：座標変換後のベクトル @SEOF

"""

function seoftoscsf(v, qua)
    vout = zeros(1,3)
    println(v)
    println(qua)
    vout_q = qua * vec(v) / qua

    vout[1]=vout_q.q1
    vout[2]=vout_q.q2
    vout[3]=vout_q.q3

    return vout
end

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
#import Pkg
#Pkg.add("Gadfly")

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

using PyPlot
using DelimitedFiles

printstyled("STARTING THE TEST...\n",bold=true,color=:light_cyan)

printstyled("Loading the pictures... ",bold=true,color=:bold)
m = ImgtoFloat("test.png")
mi = ImgtoFloat("test_ID.png")
printstyled("DONE\n",bold=false,color=:yellow)

printstyled("Creating Lighting Environnement... ",bold=true,color=:bold)
sl = SceneLights()
printstyled("DONE\n",bold=false,color=:yellow)

printstyled("Setting the layers... ",bold=true,color=:bold)
spwnLayer!(sl,"land",5.0)
spwnLayer!(sl,"rocks",-0.5)
printstyled("DONE\n",bold=false,color=:yellow)

printstyled("Spawning the lights... ",bold=true,color=:bold)
spwnLight!(sl,"jose",[2,2],"central",20.0,0.0,[1.0,1.0,0.0])
spwnLight!(sl,"mob",[5,18],"land",20.0,0.0,[1.0,0.0,1.0])
printstyled("DONE\n",bold=false,color=:yellow)

printstyled("Spawning the objects... ",bold=true,color=:bold)
spwnObject!(sl,"carreau",[4,4],"central","test.png","test_ID.png")
printstyled("DONE\n",bold=false,color=:yellow)

askSceneLights(sl)

printstyled("Processing... ",bold=true,color=:bold)
pp,dd = Ray_ID(sl,"carreau")
printstyled("DONE\n",bold=false,color=:yellow)

printstyled("Displaying the pictures... ",bold=true,color=:bold)
subplot(1,3,1)
imshow(m)
subplot(1,3,2)
imshow(mi)
subplot(1,3,3)
imshow(pp)
scatter(sl.Light["jose"].Coord[1]-sl.Object["carreau"].Coord[1],sl.Light["jose"].Coord[2]-sl.Object["carreau"].Coord[2])
scatter(sl.Light["mob"].Coord[1]-sl.Object["carreau"].Coord[1],sl.Light["mob"].Coord[2]-sl.Object["carreau"].Coord[2])
printstyled("DONE\n",bold=false,color=:yellow)

printstyled("ENDING THE TEST...\n",bold=true,color=:light_magenta)

#How to display
#bring light


text:
unfliprotate(img) = rotl90(img[:, end:-1:1])
fliprotate(img) = rotr90(img)[:, end:-1:1]
# tmix=frames=15:weights="1 1 1 1 1 1 1 1 1 1 1 1 1 1 1"
function getimg(file, t)
    imgraw = FFMPEG.ffmpeg() do exe
        read(`$exe -loglevel 8 -ss $t -i $file -vf yadif=1,scale=sar"*"iw:ih -vframes 1 -f image2pipe -`)
    end
    return fliprotate(ImageMagick.load_(imgraw))
end

function plotrawpoi(poi::POI, file)
    ind = Node(1)
    img = lift(ind) do i
        getimg(poi.video, time(poi, i))
    end
    timestamp = lift(ind) do i
        string("time stamp: ", round(time(poi)[i], digits = 3), " sec")
    end
    spot = lift(ind) do i
        space(poi, i)
    end
    fig = Figure()
    ax = fig[1, 1] = Axis(fig, aspect = DataAspect(), yreversed = true, xlabel = "X (pixel)", ylabel = "Y (pixel)", title = timestamp)
    # scene, layout = layoutscene(0)
    # ax = layout[1,1] = LAxis(scene, aspect = DataAspect(), yreversed = true, xlabel = "X (pixel)", ylabel = "Y (pixel)", title = timestamp)
    image!(ax, img)
    # layout[0, :] = LText(scene, timestamp)
    scatter!(ax, spot, color = :transparent, strokecolor = :red, strokewidth = 3)
    tightlimits!(ax)
    if length(time(poi)) == 1
        recordit(file, fig)
    else
        recordit(poi, file, ind, fig)
        # recordmovie(poi, file, ind, scene)
    end
end

recordit(file, scene) = save("$file.png", scene)

function recordit(poi, file, ind, scene)
    for (j,i) in enumerate(round.(Int, range(1, stop = length(time(poi)) - 1, length = 5)))
        ind[] = i
        save("$file$j.png", scene)
    end
end

function recordmovie(poi, file, ind, scene)
    GLMakie.Makie.inline!(true)
    record(scene, "$file.mp4", eachindex(time(poi))) do i
        ind[] = i
    end
    GLMakie.Makie.inline!(false)
end

decompose(imgw) = (UnitRange.(axes(imgw)), parent(imgw))
function plotcalibration(calibration::Calibration, calib, file)
    fig = Figure()
    ax = fig[1, 1] = Axis(fig, aspect = DataAspect(), yreversed = true, xlabel = "X (pixel)", ylabel = "Y (pixel)", title = "Raw", backgroundcolor = :black)
    # scene, layout = layoutscene(0)
    # ax = layout[1,1] = LAxis(scene, aspect = DataAspect(), yreversed = true, xlabel = "X (pixel)", ylabel = "Y (pixel)", title = "Raw", backgroundcolor = :black)
    img = getimg(calibration.video, calibration.extrinsic)
    image!(ax, img)
    tightlimits!(ax)
    ax = fig[1, 2] = Axis(fig, aspect = DataAspect(), yreversed = true, xlabel = "X (cm)", ylabel = "Y (cm)", title = "Calibrated", backgroundcolor = :black)
    # ax = layout[1,2] = LAxis(scene, aspect = DataAspect(), yreversed = true, xlabel = "X (cm)", ylabel = "Y (cm)", title = "Calibrated", backgroundcolor = :black)
    indices, imgw = decompose(calibrate(calib, img))
    image!(ax, indices..., imgw)
    tightlimits!(ax)
    fig[2,1:2] = Label(fig, "ϵ (± pixel) = $(round.(calib.ϵ, digits = 2))")
    save("$file.png", fig)
end

function _label_position_distances(ps)
    p1, p2 = ps
    Δ = p2 - p1
    d = norm(Δ)
    label = string(round(Int, d))
    pos = p1 + d/2*LinearAlgebra.normalize(Δ)
    return label, pos
end

function _plotpoic(ax, indices, img, poic, expected_locations)
    image!(ax, indices..., img)
    ps = [Pair(only(p1), only(p2)) for (p1, p2) in combinations(space.(values(poic)), 2)]
    lps = _label_position_distances.(ps)
    # textlayer = textlayer!(ax)
    annotations(first.(lps), Point2f0.(last.(lps)), color = :yellow, align = (:center, :center), textsize = 12)#, rotation = π/6)
    linesegments!(ax, ps, color = :white)
    opts = (marker = '+', color = :red, strokewidth = 0, markersize = 6)
    scatter!(ax, only.(space.(values(poic))); opts..., color = :green)
    labels = String.(keys(poic))
    pos = [Point(only(space(xy)) .- (0, 5)) for xy in values(poic)]
    annotations!(labels, pos, color = :yellow, align = (:center, :bottom), textsize = 12)#, rotation = π/6)
    if !isempty(expected_locations)
        scatter!(ax, Point.(values(expected_locations)); opts...)
    end
end

function plotcalibratedpoi(pois, calib, file, expected_locations, calib2)
    fig = Figure()
    # scene, layout = layoutscene(0)
    # poic = Dict(k => calibrate(calib, v) for (k,v) in poi)
    x = first(values(pois))
    img = getimg(x.video, time(x))
    imgww = calibrate(calib, img)
    indices, imgw = decompose(imgww)
    ax1 = fig[1, 1] = Axis(fig, aspect = DataAspect(), yreversed = true, xlabel = "X (cm)", ylabel = "Y (cm)", title = "Calibrated", backgroundcolor = :black)
    # ax1 = layout[1,1] = LAxis(scene, aspect = DataAspect(), yreversed = true, xlabel = "X (cm)", ylabel = "Y (cm)", title = "Calibrated", backgroundcolor = :black)
    _plotpoic(ax1, indices, imgw, pois, expected_locations)
    for v in values(pois)
        # map!(calib2, space(v), space(v))
        calibrate!(v, calib2)# = map!(xy -> calibrate(c, xy), space(poi), space(poi))
    end
    # imgw = warp(OffsetArray(imgw, indices...), inv(calib2))
    # imgw = warp(imgw, inv(calib2), indices)
    # indices = parent.(axes(imgw))
    # imgw = parent(imgw)
    #
    # @show calib2.tform
    # @show typeof(calib2.tform)
    # @show typeof(imgww)
    # @show axes(imgww)
    # @show calib2.itform
    # @show typeof(calib2.itform)
    indices, imgw = decompose(calibrate(calib2, imgww))
    ax2 = fig[1, 2] = Axis(fig, aspect = DataAspect(), yreversed = true, xlabel = "X (cm)", ylabel = "Y (cm)", title = "& corrected", backgroundcolor = :black)
    # ax2 = layout[1,2] = LAxis(scene, aspect = DataAspect(), yreversed = true, xlabel = "X (cm)", ylabel = "Y (cm)", title = "& corrected", backgroundcolor = :black)
    _plotpoic(ax2, indices, imgw, pois, expected_locations)
    linkaxes!(ax1, ax2)
    mx = minimum(first.(only.(space.(values(pois))))) - 100
    Mx = maximum(first.(only.(space.(values(pois))))) + 100
    my = minimum(last.(only.(space.(values(pois))))) - 100
    My = maximum(last.(only.(space.(values(pois))))) + 100
    limits!(ax1, mx, Mx, My, my)
    hideydecorations!(ax2)
    hidexdecorations!(ax1, ticklabels = false, ticks = false, grid = false)
    hidexdecorations!(ax2, ticklabels = false, ticks = false, grid = false)
    fig[2,1:2] = Label(fig, "X (cm)")
    save(joinpath(file, "calibrated POIs.png"), fig)
end

# function textlayer!(ax::Axis)
#     pxa = lift(AbstractPlotting.zero_origin, ax.scene.px_area)
#     Scene(ax.scene, pxa, raw = true, camera = campixel!)
# end
#
# function AbstractPlotting.annotations!(textlayer::Figure, ax::Axis, texts, positions; kwargs...)
#     positions = positions isa Observable ? positions : Observable(positions)
#     screenpositions = lift(positions, ax.scene.camera.projectionview, ax.scene.camera.pixel_space) do positions, pv, pspace
#         p4s = to_ndim.(Vec4f0, to_ndim.(Vec3f0, positions, 0.0), 1.0)
#         p1m1s = [pv *  p for p in p4s]
#         projected = [inv(pspace) * p1m1 for p1m1 in p1m1s]
#         pdisplay = [Point2(p[1:2]...) for p in projected]
#     end
#     annotations!(textlayer, texts, screenpositions; kwargs...)
# end


text:
function create_simple_discrete_bandit(num_actions::Int; noise=1.0)
    A = 1:num_actions
    function simple_reward(a, noise)
        return a + randn() * noise
    end
    r = a -> simple_reward(a, noise)
    meta = Dict{Symbol, Any}()
    meta[:minreward] = -Inf
    meta[:maxreward] = Inf
    meta[:minobjective] = 1
    meta[:maxobjective] = num_actions
    b = BanditProblem(A, r, meta)

    return b
end

function create_simple_discrete_contextualbandit(num_actions::Int; noise_x=1.0, noise_r=1.0)
    X = [-2noise_x 2noise_x]
    A = 1:num_actions

    function sample_context(noise)
        return clamp(randn() * noise, -2*noise, 2*noise)
    end

    function simple_reward(x, a, noise)
        return x*a + randn() * noise
    end

    d = ()->sample_context(noise_x)
    r = (x,a) -> simple_reward(x,a, noise_r)
    meta = Dict{Symbol, Any}()
    meta[:minreward] = -Inf
    meta[:maxreward] = Inf
    meta[:minobjective] = -(√2 / π)
    meta[:maxobjective] = num_actions * (√2 / π)
    b = ContextualBanditProblem(X,A,d,r,meta)
    return b
end


text:
# Command: Find Elements
"""
    Elements(session::Session, location_strategy::AbstractString, value::AbstractString)::Vector{Element}

Find Elements

See also: [`Element`](@ref)
"""
function Elements(
    session::Session,
    location_strategy::AbstractString,
    value::AbstractString,
)::Vector{Element}
    location_strategy ∈
    ["css selector", "link text", "partial link text", "tag name", "xpath"] ||
    throw(ArgumentError("""location_strategy must be "css selector", "link text", "partial link text", "tag name" or "xpath"."""))
    @unpack addr, id = session
    response = HTTP.post(
        "$addr/session/$id/elements",
        [("Content-Type" => "application/json; charset=utf-8")],
        JSON3.write(Dict("using" => location_strategy, "value" => value)),
    )
    @assert response.status == 200
    Element.(session, first(values(node)) for node in JSON3.read(response.body).value)
end