text:
using IRTools: IR,
               Variable,
               Pipe,
               xcall,
               var,
               prewalk,
               postwalk,
               blocks,
               predecessors,
               successors,
               argument!,
               arguments,
               branches,
               exprtype,
               insertafter!,
               finish,
               expand!,
               prune!,
               substitute!,
               substitute,
               block,
               block!,
               branch!,
               return!,
               stmt
using Base: @get!

@inline tuple_va(N, xs) = xs
@inline tuple_va(N, x, xs...) = (x, tuple_va(N, xs...)...)
@inline tuple_va(::Val{N}, ::Nothing) where {N} = ntuple(_ -> nothing, Val(N))

iscall(x, m::Module, n::Symbol) = isexpr(x, :call) && x.args[1] == GlobalRef(m, n)

gradindex(x, i) = x[i]
gradindex(::Nothing, i) = nothing
xgetindex(x, i...) = xcall(Base, :getindex, x, i...)
xgradindex(x, i) = xcall(Zygote, :gradindex, x, i)

normalise!(ir) = ir |> IRTools.merge_returns!

function instrument_new!(ir, v, ex)
    isexpr(ex, :new) ? (ir[v] = xcall(Zygote, :__new__, ex.args...)) :
    isexpr(ex, :splatnew) ? (ir[v] = xcall(Zygote, :__splatnew__, ex.args...)) : ex
end

# Hack to work around fragile constant prop through overloaded functions
unwrapquote(x) = x
unwrapquote(x::QuoteNode) = x.value

is_literal_getproperty(ex) =
    (
     iscall(ex, Base, :getproperty) ||
     iscall(ex, Core, :getfield) || iscall(ex, Base, :getfield)
    ) && ex.args[3] isa Union{QuoteNode,Integer}

function instrument_getproperty!(ir, v, ex)
    is_literal_getproperty(ex) ?
    (
     ir[v] = xcall(Zygote, :literal_getproperty, ex.args[2], Val(unwrapquote(ex.args[3])))
    ) :
    ex
end

is_literal_getindex(ex) =
    iscall(ex, Base, :getindex) &&
    length(ex.args) == 3 && ex.args[3] isa Union{Integer,QuoteNode}

function instrument_getindex!(ir, v, ex)
    is_literal_getindex(ex) ?
    (ir[v] = xcall(Zygote, :literal_getindex, ex.args[2], Val(unwrapquote(ex.args[3])))) :
    ex
end

is_literal_iterate(ex) =
    iscall(ex, Base, :indexed_iterate) &&
    length(ex.args) >= 3 && ex.args[3] isa Union{Integer,QuoteNode}

function instrument_iterate!(ir, v, ex)
    is_literal_iterate(ex) ?
    (
     ir[v] = xcall(
         Zygote,
         :literal_indexed_iterate,
         ex.args[2],
         Val(unwrapquote(ex.args[3])),
         ex.args[4:end]...,
     )
    ) :
    ex
end

function instrument_literals!(ir, v, ex)
    ex = instrument_getproperty!(ir, v, ex)
    ex = instrument_getindex!(ir, v, ex)
    ex = instrument_iterate!(ir, v, ex)
end

function istrackable(x)
    x isa GlobalRef && x.mod ∉ (Base, Core) || return false
    isconst(x.mod, x.name) || return true
    x = getfield(x.mod, x.name)
    !(x isa Type || sizeof(x) == 0)
end

function instrument_global!(ir, v, ex)
    if istrackable(ex)
        ir[v] = xcall(Zygote, :unwrap, QuoteNode(ex), ex)
    else
        ir[v] =
            prewalk(ex) do x
                istrackable(x) || return x
                insert!(
                    ir,
                    v,
                    stmt(xcall(Zygote, :unwrap, QuoteNode(x), x), type = exprtype(x)),
                )
            end
    end
end

function instrument(ir::IR)
    pr = Pipe(ir)
    for (v, st) in pr
        ex = st.expr
        isexpr(ex, :foreigncall) && continue
        isexpr(ex, :enter, :leave) && error("try/catch is not supported.")
        ex = instrument_new!(pr, v, ex)
        ex = instrument_literals!(pr, v, ex)
        ex = instrument_global!(pr, v, ex)
    end
    return finish(pr)
end

const BranchNumber = UInt8

function record_branches!(ir::IR)
    brs = Dict{Int,Variable}()
    for bb in blocks(ir)
        preds = predecessors(bb)
        length(preds) > 1 || continue
        brs[bb.id] = argument!(bb, BranchNumber(0), BranchNumber)
        i = length(arguments(bb))
        n = 0
        for aa in blocks(ir), br in branches(aa)
            br.block == bb.id && (arguments(br)[i] = BranchNumber(n += 1))
        end
    end
    return ir, brs
end

ignored_f(f) = f in (
    GlobalRef(Base, :not_int),
    GlobalRef(Core.Intrinsics, :not_int),
    GlobalRef(Core, :(===)),
    GlobalRef(Core, :apply_type),
    GlobalRef(Core, :typeof),
    GlobalRef(Core, :throw),
    GlobalRef(Base, :kwerr),
    GlobalRef(Core, :kwfunc),
    GlobalRef(Core, :isdefined),
)
ignored_f(ir, f) = ignored_f(f)
ignored_f(ir, f::Variable) = ignored_f(get(ir, f, nothing))

ignored(ir, ex) = isexpr(ex, :call) && ignored_f(ir, ex.args[1])
ignored(ir, ex::Variable) = ignored(ir, ir[ex])

# TODO: remove this once we don't mess with type inference
function _forward_type(Ts)
    usetyped || return Any
    all(T -> isconcretetype(T) || T <: DataType, Ts) || return Any
    T = Core.Compiler.return_type(_forward, Tuple{Context,Ts...})
    return T == Union{} ? Any : T
end

isvalidtype(jT, yT) = jT <: Tuple && length(jT.parameters) == 2 && jT.parameters[1] <: yT

function primal(ir::IR)
    pr = Pipe(ir)
    pbs = Dict{Variable,Variable}()
    argument!(pr, at = 1)
    cx = argument!(pr, Context, at = 2)
    for (v, st) in pr
        ex = st.expr
        if isexpr(ex, :call) && !ignored(ir, ex)
            yT = exprtype(ir, v)
            T = _forward_type(exprtype.((ir,), ex.args))
            if yT == Any || isvalidtype(T, yT)
                yJ = insert!(
                    pr,
                    v,
                    stmt(xcall(Zygote, :_forward, cx, ex.args...), line = ir[v].line),
                )
                pr[v] = xgetindex(yJ, 1)
                J = insertafter!(
                    pr,
                    v,
                    stmt(
                        xgetindex(yJ, 2),
                        type = T == Any ? Any : T.parameters[2],
                        line = ir[v].line,
                    ),
                )
                pbs[v] = substitute(pr, J)
            else
                yJ = insert!(pr, v, xcall(Zygote, :_forward, cx, ex.args...))
                y = insert!(pr, v, xgetindex(yJ, 1))
                J = insert!(pr, v, stmt(xgetindex(yJ, 2), line = ir[v].line))
                pr[v] = xcall(Zygote, :typeassert, y, yT)
                pbs[v] = substitute(pr, J)
            end
        end
    end
    pr = finish(pr)
    pr, brs = record_branches!(pr)
    return pr, brs, pbs
end

struct Primal
    ir::IR
    pr::IR
    varargs::Union{Int,Nothing}
    branches::Dict{Int,Variable}
    pullbacks::Dict{Variable,Variable}
end

function Primal(ir::IR; varargs = nothing)
    ir = instrument(normalise!(ir))
    pr, brs, pbs = primal(ir)
    Primal(expand!(ir), pr, varargs, brs, pbs)
end

# Backwards Pass

struct Alpha
    id::Int
end

Base.show(io::IO, x::Alpha) = print(io, "@", x.id)

alpha(x) = x
alpha(x::Variable) = Alpha(x.id)
Variable(a::Alpha) = Variable(a.id)

sig(b::IRTools.Block) =
    unique([arg for br in branches(b) for arg in br.args if arg isa Variable])
sig(pr::Primal) = Dict(b.id => sig(b) for b in blocks(pr.ir))

# TODO unreachables?
function adjointcfg(pr::Primal)
    ir = empty(pr.ir)
    return!(ir, nothing)
    for b in blocks(pr.ir)[2:end]
        block!(ir)
        preds = predecessors(b)
        rb = block(ir, b.id)
        for i = 1:length(preds)
            cond = i == length(preds) ? nothing :
                push!(rb, xcall(Base, :(!==), alpha(pr.branches[b.id]), BranchNumber(i)))
            branch!(rb, preds[i].id, unless = cond)
        end
        if !isempty(branches(b)) && branches(b)[end] == IRTools.unreachable
            branch!(rb, 0)
        end
    end
    sigs = sig(pr)
    for b in blocks(ir)[1:end-1], i = 1:length(sigs[b.id])
        argument!(b)
    end
    argument!(blocks(ir)[end])
    return ir, sigs
end

branchfor(ir, (from, to)) =
    get(filter(br -> br.block == to, branches(block(ir, from))), 1, nothing)

xaccum(ir) = nothing
xaccum(ir, x) = x
xaccum(ir, xs...) = push!(ir, xcall(Zygote, :accum, xs...))

function adjoint(pr::Primal)
    ir, sigs = adjointcfg(pr)
    for b in reverse(blocks(pr.ir))
        rb = block(ir, b.id)
        grads = Dict()
        grad(x, x̄) = push!(get!(grads, x, []), x̄)
        grad(x) = xaccum(rb, get(grads, x, [])...)
    # Backprop through (successor) branch arguments
        for i = 1:length(sigs[b.id])
            grad(sigs[b.id][i], arguments(rb)[i])
        end
    # Backprop through statements
        for v in reverse(keys(b))
            ex = b[v].expr
            if haskey(pr.pullbacks, v)
                g = push!(
                    rb,
                    stmt(
                        Expr(:call, alpha(pr.pullbacks[v]), grad(v)),
                        line = b[v].line,
                    ),
                )
                for (i, x) in enumerate(ex.args)
                    x isa Variable || continue
                    grad(x, push!(rb, stmt(xgradindex(g, i), line = b[v].line)))
                end
            elseif ex isa Core.PiNode
                grads[ex.val] = grads[v]
            elseif isexpr(ex, GlobalRef, :call, :isdefined, :inbounds, :meta)
            elseif isexpr(ex)
                push!(
                    rb,
                    stmt(
                        xcall(Base, :error, "Can't differentiate $(ex.head) expression"),
                        line = b[v].line,
                    ),
                )
            else # A literal value
                continue
            end
        end
        if b.id > 1 # Backprop through (predecessor) branch arguments
            gs = grad.(arguments(b))
            for br in branches(rb)
                br.block == 0 && continue
                br′ = branchfor(pr.ir, br.block => b.id)
                br′ == nothing && continue
                ins = br′.args
                for i = 1:length(br.args)
                    ā = [gs[j] for j = 1:length(ins) if ins[j] == sigs[br.block][i]]
                    br.args[i] = xaccum(rb, ā...)
                end
            end
        else # Backprop function arguments
            gs = [grad(arg) for arg in arguments(pr.ir)]
            Δ = push!(
                rb,
                pr.varargs == nothing ? xcall(Zygote, :tuple, gs...) :
                xcall(Zygote, :tuple_va, Val(pr.varargs), gs...),
            )
            branches(rb)[1].args[1] = Δ
        end
    end
    return ir
end

struct Adjoint
    primal::IR
    adjoint::IR
end

function Adjoint(ir::IR; varargs = nothing, normalise = true)
    pr = Primal(ir, varargs = varargs)
    adj = adjoint(pr) |> prune!
    if normalise
        permute!(adj, length(adj.blocks):-1:1)
        adj = IRTools.domorder!(adj) |> IRTools.renumber
    end
    Adjoint(pr.pr, adj)
end


text:
# # Custom potential

# We solve the 1D Gross-Pitaevskii equation with a custom potential.
# This is similar to [Gross-Pitaevskii equation in one dimension](@ref) and we
# show how to define local potentials attached to atoms, which allows for
# instance to compute forces.
using DFTK
using LinearAlgebra

# First, we define a new element which represents a nucleus generating a custom
# potential
struct ElementCustomPotential <: DFTK.Element
    pot_real::Function      # Real potential
    pot_fourier::Function   # Fourier potential
end

# We need to extend two methods to access the real and Fourier forms of
# the potential during the computations performed by DFTK
function DFTK.local_potential_fourier(el::ElementCustomPotential, q::Real)
    return el.pot_fourier(q)
end
function DFTK.local_potential_real(el::ElementCustomPotential, r::Real)
    return el.pot_real(r)
end

# We set up the lattice. For a 1D case we supply two zero lattice vectors
a = 10
lattice = a .* [[1 0 0.]; [0 0 0]; [0 0 0]];

# In this example, we want to generate two Gaussian potentials generated by
# two nuclei localized at positions ``x_1`` and ``x_2``, that are expressed in
# ``[0,1)`` in fractional coordinates. ``|x_1 - x_2|`` should be different from
# ``0.5`` to break symmetry and get nonzero forces.
x1 = 0.2
x2 = 0.8;
# We define the width of the Gaussian potential generated by one nucleus
L = 0.5;
# We set the potential in its real and Fourier forms
pot_real(x) = exp(-(x/L)^2)
pot_fourier(q::T) where {T <: Real} = exp(- (q*L)^2 / 4);
# And finally we build the elements and set their positions in the `atoms`
# array. Note that in this example `pot_real` is not required as all applications
# of local potentials are done in the Fourier space.
nucleus = ElementCustomPotential(pot_real, pot_fourier)
atoms = [nucleus => [x1*[1,0,0], x2*[1,0,0]]];

# Setup the Gross-Pitaevskii model
C = 1.0
α = 2;
n_electrons = 1  # Increase this for fun
terms = [Kinetic(),
         AtomicLocal(),
         PowerNonlinearity(C, α),
]
model = Model(lattice; atoms=atoms, n_electrons=n_electrons, terms=terms,
              spin_polarization=:spinless);  # use "spinless electrons"

# We discretize using a moderate Ecut and run a SCF algorithm to compute forces
# afterwards. As there is no ionic charge associated to `nucleus` we have to specify
# a starting density and we choose to start from a zero density.
Ecut = 500
basis = PlaneWaveBasis(model, Ecut, kgrid=(1, 1, 1))
ρ = zeros(complex(eltype(basis)), basis.fft_size)
scfres = self_consistent_field(basis, tol=1e-8, ρ=from_fourier(basis, ρ))
scfres.energies
# Computing the forces can then be done as usual:
hcat(forces(scfres)...)

# Extract the converged total local potential
tot_local_pot = DFTK.total_local_potential(scfres.ham)[:, 1, 1]; # use only dimension 1

# Extract other quantities before plotting them
ρ = real(scfres.ρ.real)[:, 1, 1]  # converged density
ψ_fourier = scfres.ψ[1][:, 1];    # first kpoint, all G components, first eigenvector
ψ = G_to_r(basis, basis.kpoints[1], ψ_fourier)[:, 1, 1]
ψ /= (ψ[div(end, 2)] / abs(ψ[div(end, 2)]));

using Plots
x = a * vec(first.(DFTK.r_vectors(basis)))
p = plot(x, real.(ψ), label="real(ψ)")
plot!(p, x, imag.(ψ), label="imag(ψ)")
plot!(p, x, ρ, label="ρ")
plot!(p, x, tot_local_pot, label="tot local pot")

ext: jl


text:
include("ImpedanceAnalyzer_Functions.jl")

function MagPhaseBoxDotImpedanceStats_AgainstKnown(ZMagVec::Vector,ZPhaseVec::Vector,KnownImpedance::Complex,newPlot,labelVal,CVal)

    ZMagNormalized = ZMagVec ./ abs(KnownImpedance)

    KnownRealVec = repeat([real(KnownImpedance)], Int(length(ZMagVec)))
    KnownImagVec = repeat([(imag(KnownImpedance))], Int(length(ZMagVec)))
    KnownMagVec = repeat([abs(KnownImpedance)], Int(length(ZMagVec)))
    KnownPhaseVec = repeat([(angle(KnownImpedance)*180/pi)], Int(length(ZMagVec)))
    if newPlot
        Tmp = display(violin([log10.(KnownMagVec) KnownPhaseVec],[ZMagNormalized ZPhaseVec], layout=2,label = [labelVal*" Mag" labelVal*" Phase"],color=CVal,legend=[true false],xlabel=["Log₁₀(|Z|)" "∠Z [deg]"]))
        # display(StatsPlots.scatter(ZMagVec,KnownMagVec))
    else
        display(violin!([log10.(KnownMagVec) KnownPhaseVec],[ZMagNormalized ZPhaseVec], layout=2,label = [labelVal*" Mag" labelVal*" Phase"],color=CVal,legend=[true false]))
        # display(StatsPlots.scatter!(ZMagVec,KnownMagVec))
    end

end

function PlotManyElements(PathArray,KnownImpedanceVec,NameVec,ColorVec)

    for i in 1:length(KnownImpedanceVec)

        StatsData = UnpackStatsRead(PathArray[i])
        ZPhaseVec = StatsData[2]
        ZMagVec = StatsData[3]

        KnownImpedance = KnownImpedanceVec[i]
        if i==1
            MagPhaseBoxDotImpedanceStats_AgainstKnown(ZMagVec,ZPhaseVec,Complex(KnownImpedance),true,NameVec[i],ColorVec[i])
        else
            MagPhaseBoxDotImpedanceStats_AgainstKnown(ZMagVec,ZPhaseVec,Complex(KnownImpedance),false,NameVec[i],ColorVec[i])
        end

    end

end

function PlotMeanErrorStd(PathArray,KnownImpedanceVec,NameVec)
    ZStdMagVec =zeros(length(PathArray))
    ZMeanMagVec =zeros(length(PathArray))
    for i in 1:length(KnownImpedanceVec)

        StatsData = UnpackStatsRead(PathArray[i])
        ZMeanMagVec[i] = StatsData[1]["Mean_ZMag"]
        ZStdMagVec[i] = StatsData[1]["Std_ZMag"]
    end

    Error = (abs.(KnownImpedanceVec) .- ZMeanMagVec) ./ abs.(KnownImpedanceVec)
    ZStdMagVecNormalized = ZStdMagVec ./ abs.(KnownImpedanceVec)
    scatter([abs.(KnownImpedanceVec) abs.(KnownImpedanceVec)],[Error ZStdMagVecNormalized], layout = (2,1),xscale = :log10,xlabel=["" "|Z|"],ylabel = ["Error[UL]" "Normalized Std. Dev |Z|"],legend = false)

end

PathArray = [
"C:\\Users\\MPI\\Documents\\GitHub\\Simple_Impedance_Analyzer\\Julia Data Processing\\Accuracy_Testing_Data\\2022_01_08_17_39_02_R100m.tdms"

"C:\\Users\\MPI\\Documents\\GitHub\\Simple_Impedance_Analyzer\\Julia Data Processing\\Accuracy_Testing_Data\\2022_01_08_10_00_22_R1Ohms.tdms"

"C:\\Users\\MPI\\Documents\\GitHub\\Simple_Impedance_Analyzer\\Julia Data Processing\\Accuracy_Testing_Data\\2022_01_08_17_49_59_R10m.tdms"

"C:\\Users\\MPI\\Documents\\GitHub\\Simple_Impedance_Analyzer\\Julia Data Processing\\Accuracy_Testing_Data\\2022_01_08_10_10_25_R100kOhms.tdms"

"C:\\Users\\MPI\\Documents\\GitHub\\Simple_Impedance_Analyzer\\Julia Data Processing\\Accuracy_Testing_Data\\2022_01_08_10_15_25_R10kOhms.tdms"

"C:\\Users\\MPI\\Documents\\GitHub\\Simple_Impedance_Analyzer\\Julia Data Processing\\Accuracy_Testing_Data\\2022_01_08_10_20_21_R264Ohms.tdms"]

KnownImpedanceVec = [PolarDeg2Complex(0.0996,.75), PolarDeg2Complex(1.0075,.9) ,PolarDeg2Complex(0.00983,2), PolarDeg2Complex(99800,0.01), PolarDeg2Complex(9838,0.01) , PolarDeg2Complex(264.8,0)]

NameVec = ["100mΩ" , "1Ω" , "10mΩ" , "100kΩ" ,"10kΩ" , "264Ω"]
ColorVec = ["blue" , "red" , "green" , "purple" ,"yellow" , "cyan"]

PathArray = [
"C:\\Users\\MPI\\Documents\\GitHub\\Simple_Impedance_Analyzer\\Julia Data Processing\\Accuracy_Testing_Data\\2022_01_08_10_30_11_C10nF.tdms"

"C:\\Users\\MPI\\Documents\\GitHub\\Simple_Impedance_Analyzer\\Julia Data Processing\\Accuracy_Testing_Data\\2022_01_08_10_35_17_C3n9F.tdms"

"C:\\Users\\MPI\\Documents\\GitHub\\Simple_Impedance_Analyzer\\Julia Data Processing\\Accuracy_Testing_Data\\2022_01_08_10_40_13_C22nF.tdms"

"C:\\Users\\MPI\\Documents\\GitHub\\Simple_Impedance_Analyzer\\Julia Data Processing\\Accuracy_Testing_Data\\2022_01_08_10_45_46_C80nF.tdms"]

KnownImpedanceVec = [PolarDeg2Complex(1651,-89.98) ,PolarDeg2Complex(4458,-89.11), PolarDeg2Complex(709.5,-89.98), PolarDeg2Complex(195.95,-89.45)]

NameVec = ["10nF" , "3.9nF" , "22nF" , "80nF"]
ColorVec = ["blue" , "red" , "green" , "purple"]


text:
function myreverse(phrase)

end


text:
using Dates, Luxor
french_months = ["janvier", "février", "mars", "avril","mai", "juin","juillet", "août", "septembre", "octobre","novembre", "décembre"];
french_monts_abbrev=["janv","févr","mars","avril","mai","juin","juil","août","sept","oct","nov","déc"];
french_days=["lundi","mardi","mercredi","jeudi","vendredi","samedi","dimanche"];
Dates.LOCALES["french"] = Dates.DateLocale(french_months,french_monts_abbrev,french_days, [""]);
t=now()

function drawneedle(l)

    goldblend = Luxor.blend(Luxor.Point(-200, 0), Luxor.Point(200, 0))
    Luxor.addstop(goldblend, 0.0,  "gold4")
    Luxor.addstop(goldblend, 0.5,  "gold3")
    p0  = Luxor.Point(1,1)
    pp1 = Luxor.Point(5,l/3)
    pp2 = Luxor.Point(0,l)
    pp3 = Luxor.Point(-5,l/3)
    Luxor.setline(3)
    ptlist = [p0,pp1,pp2,pp3]
    ptlist2 = [p0,pp1+(2, 0),pp2,pp3-(2, 0)]
    #Luxor.sethue("ivory")
    Luxor.setblend(goldblend)
    Luxor.poly(ptlist, :fill, close=true)
    Luxor.sethue("black")
    Luxor.poly(ptlist2, :stroke, close=true)

end
function drawclock(radius, action=:none)
    t  = Dates.now()
    txtday = Dates.dayname(t ; locale="french")
    txtmonth = Dates.monthname(t ; locale="french")
    s  = Dates.second(t)
    m  = Dates.minute(t)
    h  = Dates.hour(t)
    p0 = Luxor.Point(0,0)
    Luxor.origin()
    Luxor.setopacity(0.1)
    #Luxor.background("lavender")
    bgblend = Luxor.blend(Luxor.Point(0, 0), 5, Luxor.Point(0, 0), 150)
    Luxor.addstop(bgblend, 0.0,  "seagreen")
    Luxor.addstop(bgblend, 0.20, "seagreen1")
    Luxor.addstop(bgblend, 0.30, "aquamarine1")
    Luxor.addstop(bgblend, 0.60,  "aquamarine2")
    Luxor.addstop(bgblend, 0.70,  "aquamarine3")
    Luxor.addstop(bgblend, 0.85,  "aquamarine3")
    setblend(bgblend)
    #box(O, 400, 100, :fill)
    Luxor.circle(0,0, radius+60, :fill)
    Luxor.setopacity(0.8)
    Luxor.sethue("black")
    Luxor.setline(1)
    Luxor.circle(0,0, radius+60, :stroke)
    #draw 12 lines
    ht = 12
    for theta in range(0, step=pi/6, length=12)
        Luxor.@layer begin
            Luxor.rotate(theta)
            Luxor.translate(0, -radius-40)
            Luxor.sethue("black")
            Luxor.setline(3)
            Luxor.line(p0,Luxor.Point(0,10), :stroke)
            Luxor.translate(0, -2)
            fontface("Agenda-Black")
            fontsize(radius/5)
            Luxor.text("$ht" ,0,0,halign=:center)
            if ht >11
                ht = 0
            end
            ht += 1
        end
    end
    # draw quarters
    for theta in range(0, step=pi/2, length=4)
        starsize = radius/10
        Luxor.@layer begin
            Luxor.rotate(theta)
            Luxor.translate(0, -radius-20)
            Luxor.sethue("goldenrod4")
            Luxor.star(0,0, starsize,6, 0.5, 0, :fill)
        end
    end
    # draw seconds
    Luxor.rotate(0)
    for theta in range(0, step=pi/30, length=s+1)
            Luxor.origin()
            Luxor.rotate(theta)
            Luxor.translate(0, -radius-60)
            Luxor.sethue("red")
            Luxor.circle(0,0, 2, :fill)
    end

    #draw drawneedle for minutes
    Luxor.origin()
    Luxor.rotate(pi+(m*(pi/30)))
    drawneedle(radius+25)

    #draw drawneedle for hour
    deltah = m / 90
    Luxor.origin()
    Luxor.rotate(pi+((deltah + h) * (pi/6)))
    drawneedle(radius-15)

    Luxor.origin()
    Luxor.sethue("black")
    Luxor.circle(0,0, 10, :fill)
    Luxor.sethue("darkblue")
    fontface("Agenda-Black")
    fontsize(20)
    Luxor.text(Dates.format(t, "E d U";locale="french"),0,-20 ,halign=:center)
    Luxor.text(Dates.format(t, "yyyy"),0, -40 ,halign=:center)
    fontsize(30)
    Luxor.text(Dates.format(t, " HH:MM:SS"),-10,50,halign=:center)
end

max_stars_repo_name: nodrygo/jAnalogAlarmClock

max_issues_repo_name: nodrygo/jAnalogAlarmClock

max_forks_repo_name: nodrygo/jAnalogAlarmClock


text:
# using BoxCoxTrans
import YeoJohnsonTrans as YJT
using Tullio

function YJT.transform(X::AbstractArray, λs::Vector)
    Y = deepcopy(X)
    for j in size(X, 1)
        Y[j, :] = YJT.transform(X[j, :], λs[j])
    end
    return Y
end

"""
Yeo-Johnson transformation according to a matrix X
"""
function yeojohnson(X::AbstractArray)
    λs = Array{Float64, 1}(undef, size(X, 1))
    for (i, x) ∈ eachrow(X) |> enumerate
        λs[i], details = YJT.lambda(x; interval=(-10.0, 10.0))
    end
    return Y -> YJT.transform(Y, λs)
end
export yeojohnson


text:
#
println("Dc) Test of the PhotoIonization module with ASF from an internally generated initial- and final-state multiplet.")
#
setDefaults("print summary: open", "zzz-PhotoIonization.sum")
setDefaults("method: continuum, asymptotic Coulomb")  ## setDefaults("method: continuum, Galerkin")
setDefaults("method: normalization, pure Coulomb")   ## setDefaults("method: normalization, pure Coulomb")    setDefaults("method: normalization, pure sine")

wa = Atomic.Computation("xx",  Nuclear.Model(20.); grid=JAC.Radial.Grid("grid: by given parameters"; rnt = 2.0e-5, h = 5.0e-2, hp = 2.0e-2, NoPoints = 800),
                        initialConfigs=[Configuration("[Ne] 3s^2 3p^6")],
                        finalConfigs  =[Configuration("[Ne] 3s^2 3p^5"), Configuration("[Ne] 3s 3p^6") ],
                        process = JAC.Photo,
                        processSettings=PhotoIonization.Settings([E1, M1], [JAC.UseCoulomb, JAC.UseBabushkin], [50., 60., 70., 80., 90., 100.], false, true, true, true,
                                        false, Tuple{Int64,Int64}[], ExpStokes(1., 0., 0.)) )

wb = perform(wa)
setDefaults("print summary: close", "")

#==
                        processSettings=PhotoIonization.Settings([E1, M1], [JAC.UseCoulomb, JAC.UseBabushkin], [3000., 3200.], false, true, true, true,
                                        true, Tuple{Int64,Int64}[(1,1)]) )
==#


text:
@with_kw struct SimpleLightDark <: POMDPs.POMDP{Int,Int,Float64}
    discount::Float64       = 0.95
    correct_r::Float64      = 100.0
    incorrect_r::Float64    = -100.0
    light_loc::Int          = 10
    radius::Int             = 60
end
POMDPs.discount(p::SimpleLightDark) = p.discount
POMDPs.isterminal(p::SimpleLightDark, s::Number) = !(s in -p.radius:p.radius)

const ACTIONS = [-10, -1, 0, 1, 10]
POMDPs.actions(p::SimpleLightDark) = ACTIONS
POMDPs.n_actions(p::SimpleLightDark) = length(actions(p))
const ACTION_INDS = Dict(a=>i for (i,a) in enumerate(actions(SimpleLightDark())))
POMDPs.actionindex(p::SimpleLightDark, a::Int) = ACTION_INDS[a]

POMDPs.states(p::SimpleLightDark) = -p.radius:p.radius + 1
POMDPs.n_states(p::SimpleLightDark) = length(states(p))
POMDPs.stateindex(p::SimpleLightDark, s::Int) = s+p.radius+1

function POMDPs.transition(p::SimpleLightDark, s::Int, a::Int)
    if a == 0
        return SparseCat(SVector(p.radius+1), SVector(1.0))
    else
        return SparseCat(SVector(clamp(s+a, -p.radius, p.radius)), SVector(1.0))
    end
end

POMDPs.observation(p::SimpleLightDark, sp) = Normal(sp, abs(sp - p.light_loc) + 0.0001)

function POMDPs.reward(p::SimpleLightDark, s, a)
    if a == 0
        return s == 0 ? p.correct_r : p.incorrect_r
    else
        return -1.0
    end
end

function POMDPs.initialstate_distribution(p::SimpleLightDark)
    ps = ones(2*div(p.radius,2)+1)
    ps /= length(ps)
    return SparseCat(div(-p.radius,2):div(p.radius,2), ps)
end

@with_kw struct DSimpleLightDark <: POMDPs.POMDP{Int, Int, Int}
    sld::SimpleLightDark = SimpleLightDark()
    binsize::Float64     = 1.0
end

POMDPs.generate_o(p::DSimpleLightDark, sp, rng::AbstractRNG) = floor(Int, rand(rng, observation(p.sld, sp))/p.binsize)
POMDPs.generate_o(p::DSimpleLightDark, a, sp, rng::AbstractRNG) = generate_o(p, sp, rng)
function POMDPs.generate_sor(p::DSimpleLightDark, s, a, rng::AbstractRNG)
    sp = generate_s(p, s, a, rng)
    o = generate_o(p, sp, rng)
    r = reward(p, s, a, sp)
    return sp, o, r
end

function POMDPModelTools.obs_weight(p::DSimpleLightDark, sp::Int, o::Int)
    cod = observation(p.sld, sp)
    return cdf(cod, (o+1)*p.binsize) - cdf(cod, o*p.binsize)
end

POMDPs.discount(p::DSimpleLightDark) = discount(p.sld)
POMDPs.isterminal(p::DSimpleLightDark, s::Number) = isterminal(p.sld, s)
POMDPs.actions(p::DSimpleLightDark) = actions(p.sld)
POMDPs.n_actions(p::DSimpleLightDark) = n_actions(p.sld)
POMDPs.actionindex(p::DSimpleLightDark, a::Int) = actionindex(p.sld, a)
POMDPs.states(p::DSimpleLightDark) = states(p.sld)
POMDPs.n_states(p::DSimpleLightDark) = n_states(p.sld)
POMDPs.stateindex(p::DSimpleLightDark, s::Int) = stateindex(p.sld, s)
POMDPs.transition(p::DSimpleLightDark, s::Int, a::Int) = transition(p.sld, s, a)
POMDPs.reward(p::DSimpleLightDark, s, a) = reward(p.sld, s, a)
POMDPs.initialstate_distribution(p::DSimpleLightDark) = initial_state_distribution(p.sld)

struct LDHeuristic <: Policy
    p::SimpleLightDark
    q::AlphaVectorPolicy{SimpleLightDark, Int}
    std_thresh::Float64
end

struct LDHSolver <: Solver
    q::QMDPSolver
    std_thresh::Float64
end

LDHSolver(;std_thresh::Float64=0.1, kwargs...) = LDHSolver(QMDPSolver(;kwargs...), std_thresh)

POMDPs.solve(sol::LDHSolver, pomdp::SimpleLightDark) = LDHeuristic(pomdp, solve(sol.q, pomdp), sol.std_thresh)

POMDPs.action(p::LDHeuristic, s::Int) = action(p.q, s)
Random.seed!(p::LDHeuristic, s) = p

function POMDPs.action(p::LDHeuristic, b::AbstractParticleBelief)
    s = std(particles(b))
    if s <= p.std_thresh
        return action(p.q, b)
    else
        m = mean(particles(b))
        ll = p.p.light_loc
        if m == ll
            return -1*Int(sign(ll))
        elseif abs(m-ll) >= 10
            return -10*Int(sign(m-ll))
        else
            return -Int(sign(m-ll))
        end
    end
end

struct LDSide <: Solver end

mutable struct LDSidePolicy{LD} <: Policy
    q::AlphaVectorPolicy{LD, Int}
end

POMDPs.solve(solver::LDSide, pomdp::Union{SimpleLightDark,DSimpleLightDark}) = LDSidePolicy(solve(QMDPSolver(), pomdp))
Random.seed!(p::LDSidePolicy, s) = p

function POMDPs.action(p::LDSidePolicy, b)
    if pdf(b, mode(b)) > 0.9
        return action(p.q, b)
    else
        return 10
    end
end


text:
@info "Loading test packages..."

using nnlearn, BGHMM, HMMBase, BioSequences, Distributions, Distributed, Random, Serialization, Test
import StatsFuns: logsumexp

@info "Beginning tests..."

Random.seed!(786)
O=1000;S=50

@testset "PWM source prior setup, PWM source initialisation and manipulation functions" begin
    #test dirichlet prior estimation from wm inputs
    wm_input = [.0 .2 .3 .5; .0 .2 .3 .5]
    est_dirichlet_vec = nnlearn.estimate_dirichlet_prior_on_wm(wm_input,4.0)
    @test typeof(est_dirichlet_vec) == Vector{Dirichlet{Float64}}
    for pos in 1:length(est_dirichlet_vec)
        @test isapprox(est_dirichlet_vec[pos].alpha, [0.,.8,1.2,2.0])
    end

    bad_input = wm_input .* 2
    @test_throws DomainError nnlearn.estimate_dirichlet_prior_on_wm(bad_input,4.0)

    wm_input = [.1 .2 .3 .4; .1 .2 .3 .4]
    est_dirichlet_vec = nnlearn.estimate_dirichlet_prior_on_wm(wm_input,4.0)
    @test typeof(est_dirichlet_vec) == Vector{Dirichlet{Float64}}
    for pos in 1:length(est_dirichlet_vec)
        @test est_dirichlet_vec[pos].alpha == [.4,.8,1.2,1.6]
    end

    length_range = 2:2

    #test informative/uninformative source prior vector assembly
    test_priors = nnlearn.assemble_source_priors(2, [wm_input], 4.0, length_range)
    @test length(test_priors)  == 2
    for pos in 1:length(test_priors[1])
        @test test_priors[1][pos].alpha == [.4,.8,1.2,1.6]
    end
    for pos in 1:length(test_priors[2])
        @test test_priors[2][pos].alpha == ones(4)/4
    end

    #test source wm initialisation from priors
    test_sources = nnlearn.init_logPWM_sources(test_priors, length_range)
    for source in test_sources
        for pos in 1:size(source[1])[1]
            @test isprobvec(exp.(source[1][pos,:]))
        end
    end

    #test that wm_shift is returning good shifted probvecs
    rando_dist=Dirichlet([.25,.25,.25,.25])
    for i in 1:1000
        wm=rand(rando_dist)
        new_wm=nnlearn.wm_shift(wm,Weibull(1.5,.1))
        @test isprobvec(new_wm)
        @test wm!=new_wm
    end

    #test that legal new sources are generated by permute_source_weights
    permuted_weight_sources=deepcopy(test_sources)
    permuted_weight_sources[1]=nnlearn.permute_source_weights(permuted_weight_sources[1],1.,Weibull(1.5,.1))
    permuted_weight_sources[2]=nnlearn.permute_source_weights(permuted_weight_sources[2],1.,Weibull(1.5,.1))
    @test permuted_weight_sources != test_sources
    for (s,source) in enumerate(permuted_weight_sources)
        for pos in 1:size(source[1],1)
            @test isprobvec(exp.(source[1][pos,:]))
            @test source[1][pos,:] != test_sources[s][1][pos,:]
        end
    end

    #test that get_length_params returns legal length shifts
    lls=1:10
    pr=1:5
    for i in 1:1000
        srcl=rand(1:10)
        sign, permute_length=nnlearn.get_length_params(srcl, lls, pr)
        @test pr[1]<=permute_length<=pr[end]
        @test sign==-1 || sign==1
        @test lls[1]<=(srcl+(sign*permute_length))<=lls[end]
    end

    permuted_length_sources=deepcopy(test_sources)
    permuted_length_sources[1]=nnlearn.permute_source_length(permuted_length_sources[1],test_priors[1],1:3,1:10)
    permuted_length_sources[2]=nnlearn.permute_source_length(permuted_length_sources[2],test_priors[2],1:3,1:10)
    for (s, source) in enumerate(permuted_length_sources)
        @test size(source[1],1) != size(test_sources[1][1],1)
        @test 1<=size(source[1],1)<=3
    end

    info_test_wm=[1. 0. 0. 0.
    .94 .02 .02 .02
    .82 .06 .06 .06
    .7 .1 .1 .1
    .67 .11 .11 .11
    .52 .16 .16 .16
    .4 .2 .2 .2
    .25 .25 .25 .25]

    #test eroding sources by finding most informational position and cutting off when information drops below threshold
    infovec=nnlearn.get_pwm_info(log.(info_test_wm))
    @test infovec==[2.0, 1.5774573308022544, 1.0346297041419121, 0.6432203505529603, 0.5620360019822908, 0.2403724636586433, 0.07807190511263773, 0.0]

    erosion_test_source=(log.([.25 .25 .25 .25
                         .2 .4 .2 .2
                         .7 .1 .1 .1
                         .06 .06 .06 .82
                         .7 .1 .1 .1
                         .25 .25 .25 .25]),1)

    infovec=nnlearn.get_pwm_info(erosion_test_source[1])
    start_idx, end_idx = nnlearn.get_erosion_idxs(infovec, .25, 2:8)
    @test start_idx==3
    @test end_idx==5

    eroded_pwm,eroded_prior_idx=nnlearn.erode_source(erosion_test_source,2:8,.25)
    for pos in 1:size(eroded_pwm,1)
        @test isprobvec(exp.(eroded_pwm[pos,:]))
    end
    @test eroded_prior_idx==3
    @test isapprox(exp.(eroded_pwm),[.7 .1 .1 .1
    .06 .06 .06 .82
    .7 .1 .1 .1])
end

@testset "Mix matrix initialisation and manipulation functions" begin
    #test mix matrix init
    prior_mix_test=nnlearn.init_mix_matrix((trues(2,10),0.0),2, 20)
    @test all(prior_mix_test[:,1:10])
    @test !any(prior_mix_test[:,11:20])

    @test sum(nnlearn.init_mix_matrix((falses(0,0),1.0), O, S)) == O*S
    @test sum(nnlearn.init_mix_matrix((falses(0,0),0.0), O, S)) == 0
    @test 0 < sum(nnlearn.init_mix_matrix((falses(0,0),0.5), O, S)) < O*S

    #test mix matrix decorrelation
    empty_mixvec=falses(O)
    one_mix=nnlearn.mixvec_decorrelate(empty_mixvec,1)
    @test sum(one_mix)==1

    empty_mix = falses(O,S)
    new_mix,clean=nnlearn.mix_matrix_decorrelate(empty_mix, 500)
    @test 0 < sum(new_mix) <= 500
    @test !all(clean)

    full_mix = trues(O,S)
    less_full_mix,clean=nnlearn.mix_matrix_decorrelate(full_mix, 500)

    @test O*S-sum(less_full_mix) <= 500
    @test !all(clean)

    #test matrix similarity and dissimilarity functions
    test_mix=falses(O,S)
    test_mix[1:Int(floor(O/2)),:].=true
    test_idx=3
    compare_mix=deepcopy(test_mix)
    compare_mix[:,test_idx] .= .!compare_mix[:,test_idx]
    @test nnlearn.most_dissimilar(test_mix,compare_mix)==test_idx

    src_mixvec=falses(O)
    src_mixvec[Int(ceil(O/2)):end].=true
    @test nnlearn.most_similar(src_mixvec,compare_mix)==test_idx
end

@testset "Observation setup and model scoring functions" begin
    #test a trivial scoring example
    obs=[BioSequences.DNASequence("AAAAA")]
    order_seqs = BGHMM.get_order_n_seqs(obs, 0)
    coded_seqs = BGHMM.code_seqs(order_seqs)
    obs=Array(transpose(coded_seqs))

    source_pwm = zeros(1,4)
    source_pwm[1,1] = 1
    log_pwm = log.(source_pwm)

    source_stop=5

    @test nnlearn.score_source(obs[1:5,1], log_pwm, source_stop) == [0.0 -Inf; 0.0 -Inf; 0.0 -Inf; 0.0 -Inf; 0.0 -Inf]

    #make sure revcomp_pwm is reversing pwms across both dimensions
    revcomp_test_pwm = zeros(2,4)
    revcomp_test_pwm[1,1] = 1
    revcomp_test_pwm[2,3] = 1
    log_revcomp_test_pwm = log.(revcomp_test_pwm)
    @test nnlearn.revcomp_pwm(log_revcomp_test_pwm) == [-Inf 0 -Inf -Inf
                                                        -Inf -Inf -Inf 0]

    #test a more complicated scoring example
    obs=[BioSequences.DNASequence("ATGATGATGATG")]
    order_seqs = BGHMM.get_order_n_seqs(obs, 0)
    coded_seqs = BGHMM.code_seqs(order_seqs)
    obs=Array(transpose(coded_seqs))

    source_pwm = [.7 .1 .1 .1
                  .1 .1 .1 .7
                  .1 .1 .7 .1]
    log_pwm = log.(source_pwm)

    source_start = 1
    source_stop = 10

    @test isapprox(exp.(nnlearn.score_source(obs[1:12,1], log_pwm, source_stop)),
        [.7^3 .1^3; .1^3 .1^3; .1^3 (.1*.7^2); .7^3 .1^3; .1^3 .1^3; .1^3 (.1*.7^2); .7^3 .1^3; .1^3 .1^3; .1^3 (.1*.7^2); .7^3 .1^3])

    #test scoring of multiple obs and sources
    obs=[BioSequences.DNASequence("ATGATGATGATG")
         BioSequences.DNASequence("TGATGATGATGA")]
    order_seqs = BGHMM.get_order_n_seqs(obs, 0)
    coded_seqs = BGHMM.code_seqs(order_seqs)
    obs=Array(transpose(coded_seqs))

    source_pwm_2 = [.6 .1 .1 .2
                    .2 .1 .1 .6
                    .1 .2 .6 .1]

    target_o1_s1 = [.7^3 .1^3; .1^3 .1^3; .1^3 (.1*.7^2); .7^3 .1^3; .1^3 .1^3; .1^3 (.1*.7^2); .7^3 .1^3; .1^3 .1^3; .1^3 (.1*.7^2); .7^3 .1^3]
    target_o1_s2 = [.6^3 (.2*.1^2); (.2*.1^2) (.2*.1^2); (.2*.1^2) (.2*.6^2); .6^3 (.2*.1^2); (.2*.1^2) (.2*.1^2); (.2*.1^2) (.2*.6^2); .6^3 (.2*.1^2); (.2*.1^2) (.2*.1^2); (.2*.1^2) (.2*.6^2); .6^3 (.2*.1^2)]
    target_o2_s1 = [.1^3 .1^3; .1^3 (.1*.7^2); .7^3 .1^3; .1^3 .1^3; .1^3 (.1*.7^2); .7^3 .1^3; .1^3 .1^3; .1^3 (.1*.7^2); .7^3 .1^3; .1^3 .1^3]
    target_o2_s2 = [(.2*.1^2) (.2*.1^2); (.2*.1^2) (.2*.6^2); .6^3 (.2*.1^2); (.2*.1^2) (.2*.1^2); (.2*.1^2) (.2*.6^2); .6^3 (.2*.1^2); (.2*.1^2) (.2*.1^2); (.2*.1^2) (.2*.6^2); .6^3 (.2*.1^2); (.2*.1^2) (.2*.1^2)]

    sources = [(log.(source_pwm), 1),(log.(source_pwm_2), 1)]
    dbl1_srcs=[(log.(source_pwm), 1),(log.(source_pwm), 1)]
    source_wmls = [size(source[1])[1] for source in sources]

    position_start = 1
    offsets=[0,0]
    mix_matrix = trues(2,2)

    score_mat1 = nnlearn.score_obs_sources(sources, Vector(obs[:,1]), 12, source_wmls)

    @test isapprox(exp.(score_mat1[1]),target_o1_s1)
    @test isapprox(exp.(score_mat1[2]),target_o1_s2)

    score_mat2 = nnlearn.score_obs_sources(sources, Vector(obs[:,2]), 12, source_wmls)

    @test isapprox(exp.(score_mat2[1]),target_o2_s1)
    @test isapprox(exp.(score_mat2[2]),target_o2_s2)

    #test score weaving and IPM likelihood calculations
    o=1
    bg_scores = log.(fill(.5, (12,2)))
    log_motif_expectation = log(0.5 / size(bg_scores)[1])
    obs_source_indices = findall(mix_matrix[o,:])
    obs_cardinality = length(obs_source_indices)
    penalty_sum = sum(exp.(fill(log_motif_expectation,obs_cardinality)))
    cardinality_penalty=log(1.0-penalty_sum)

    lh_target = -8.87035766177774

    o1_lh = nnlearn.weave_scores(12, view(bg_scores,:,1), score_mat1, findall(mix_matrix[1,:]), source_wmls, log_motif_expectation, cardinality_penalty)
    @test isapprox(lh_target,o1_lh)
    o2_lh = nnlearn.weave_scores(12, view(bg_scores,:,2), score_mat2, findall(mix_matrix[2,:]), source_wmls, log_motif_expectation, cardinality_penalty)

    lh,cache = nnlearn.IPM_likelihood(sources, obs, [12,12], bg_scores, mix_matrix, true,true)
    @test o1_lh==cache[1]
    @test o2_lh==cache[2]
    @test isapprox(nnlearn.lps(o1_lh,o2_lh),lh)

    #test source penalization
    dbl_score_mat= nnlearn.score_obs_sources(dbl1_srcs, Vector(obs[:,1]),12,source_wmls)

    naive=nnlearn.weave_scores(12, view(bg_scores,:,1), Vector{Matrix{Float64}}(), Vector{Int64}(), Vector{Int64}(), log_motif_expectation,cardinality_penalty)

    single=nnlearn.weave_scores(12, view(bg_scores,:,1), [dbl_score_mat[1]], [1], [source_wmls[1]], log_motif_expectation, cardinality_penalty)

    double=nnlearn.weave_scores(12, view(bg_scores,:,1), dbl_score_mat, [1,2], source_wmls, log_motif_expectation, cardinality_penalty)

    triple=nnlearn.weave_scores(12, view(bg_scores,:,1), [dbl_score_mat[1] for i in 1:3], [1,2,3], [3 for i in 1:3], log_motif_expectation, cardinality_penalty)

    @test (single-naive) > (double-single) > (triple-double)

    naive_target=-16.635532333438686
    naive = nnlearn.IPM_likelihood(sources, obs, [findfirst(iszero,obs[:,o])-1 for o in 1:size(obs)[2]], bg_scores,falses(size(obs)[2],length(sources)))
    @test naive==naive_target

    #test IPM_likelihood clean vector and cache calculations
    baselh,basecache = nnlearn.IPM_likelihood(sources, obs, [12,12], bg_scores, mix_matrix, true,true)

    clean=[true,false]

    unchangedlh,unchangedcache=nnlearn.IPM_likelihood(sources, obs, [12,12], bg_scores, mix_matrix, true, true, basecache, clean)

    changed_lh,changedcache=nnlearn.IPM_likelihood(sources, obs, [12,12], bg_scores, BitMatrix([true true; false false]), true, true, basecache, clean)

    indep_lh, indepcache=nnlearn.IPM_likelihood(sources,obs,[12,12], bg_scores, BitMatrix([true true; false false]), true, true)

    @test baselh==unchangedlh!=changed_lh==indep_lh
    @test basecache==unchangedcache!=changedcache==indepcache
end

@testset "Model permutation functions" begin
    source_pwm = [.7 .1 .1 .1
    .1 .1 .1 .7
    .1 .1 .7 .1]

    source_pwm_2 = [.6 .1 .1 .2
    .2 .1 .1 .6
    .1 .2 .6 .1]

    pwm_to_erode = [.25 .25 .25 .25
                    .97 .01 .01 .01
                    .01 .01 .01 .97
                    .01 .01 .97 .01
                    .25 .25 .25 .25]

    src_length_limits=2:5

    source_priors = nnlearn.assemble_source_priors(3, [source_pwm, source_pwm_2], 4.0, src_length_limits)
    mix_prior=0.2

    bg_scores = log.(fill(.25, (12,2)))
    obs=[BioSequences.DNASequence("ATGATGATGATG")
    BioSequences.DNASequence("TGATGATGATGA")]
    order_seqs = BGHMM.get_order_n_seqs(obs, 0)
    coded_seqs = BGHMM.code_seqs(order_seqs)
    obs=Array(transpose(coded_seqs))
    obsl=[findfirst(iszero,obs[:,o])-1 for o in 1:size(obs)[2]]

    test_model = nnlearn.ICA_PWM_model("test", source_priors, (falses(0,0),mix_prior), bg_scores, obs, src_length_limits)

    duplicate_sources=[(source_pwm,1) for i in 1:3]
    cons_check, cons_idxs = nnlearn.consolidate_check(duplicate_sources)
    @test !cons_check
    @test cons_idxs == [[2,3],[1,3],[1,2]]

    mix=falses(10,3)
    mix[1:2,1].=true
    mix[9:10,2].=true

    con_obs=hcat(obs,obs,obs,obs,obs)
    con_obsl=[12 for i in 1:10]
    con_bg=log.(fill(.5, (12,10)))

    cons_model=nnlearn.consolidate_srcs(cons_idxs, duplicate_sources, mix, con_obs, con_obsl, con_bg, source_priors, test_model.informed_sources, src_length_limits)
    @test cons_model.sources[1]==duplicate_sources[1]
    @test cons_model.sources[2]!=duplicate_sources[2]
    @test cons_model.sources[3]!=duplicate_sources[3]
    @test cons_model.mix_matrix[:,1]==[true, true, false, false, false, false, false, false, true, true]
    @test "FM from consolidate" in cons_model.flags
    @test "nofit" in cons_model.flags

    ps_model= nnlearn.permute_source(test_model, test_model.log_Li, obs, obsl, bg_scores, source_priors, 1000,.3,.5)
    @test ps_model.log_Li > test_model.log_Li
    @test ps_model.sources != test_model.sources
    @test ps_model.mix_matrix == test_model.mix_matrix
    @test "PS from test" in ps_model.flags

    pm_model= nnlearn.permute_mix(test_model, test_model.log_Li,obs, obsl, bg_scores, source_priors, 1000)
    @test pm_model.log_Li > test_model.log_Li
    @test pm_model.sources == test_model.sources
    @test pm_model.mix_matrix != test_model.mix_matrix
    @test "PM from test" in pm_model.flags

    badmix_lh=nnlearn.IPM_likelihood(test_model.sources,obs,obsl, bg_scores, trues(2,3))
    badmix=nnlearn.ICA_PWM_model("badmix",test_model.sources, test_model.informed_sources, test_model.source_length_limits, trues(2,3), badmix_lh,[""])
    psfm_model=nnlearn.perm_src_fit_mix(badmix, badmix.log_Li,obs, obsl, bg_scores, source_priors, 1000)
    @test psfm_model.log_Li > badmix.log_Li
    @test psfm_model.sources != badmix.sources
    @test psfm_model.mix_matrix != badmix.mix_matrix
    @test "PSFM from badmix" in psfm_model.flags

    fm_model=nnlearn.fit_mix(test_model,obs,obsl,bg_scores)
    @test fm_model.log_Li > test_model.log_Li
    @test fm_model.sources == test_model.sources
    @test fm_model.mix_matrix != test_model.mix_matrix
    @test "FM from test" in fm_model.flags
    @test "nofit" in fm_model.flags

    post_fm_psfm=nnlearn.perm_src_fit_mix(fm_model, test_model.log_Li, obs, obsl, bg_scores, source_priors, 1000)
    @test "PSFM from candidate" in post_fm_psfm.flags
    @test "nofit" in post_fm_psfm.flags

    rd_model=nnlearn.random_decorrelate(test_model, test_model.log_Li, obs, obsl, bg_scores, source_priors, 1000)
    @test rd_model.log_Li > test_model.log_Li
    @test rd_model.sources != test_model.sources
    @test rd_model.mix_matrix != test_model.mix_matrix
    @test "RD from test" in rd_model.flags

    rs_model=nnlearn.reinit_src(test_model, test_model.log_Li, obs, obsl, bg_scores, source_priors, (falses(0,0),.5), 1000, false)
    @test rs_model.log_Li > test_model.log_Li
    @test rs_model.sources != test_model.sources
    @test rs_model.mix_matrix != test_model.mix_matrix
    @test "RS from test" in rs_model.flags

    erosion_sources=[(log.(source_pwm),1),(log.(source_pwm_2),1),(log.(pwm_to_erode),1)]

    eroded_mix=trues(2,3)

    erosion_lh=nnlearn.IPM_likelihood(erosion_sources,obs,obsl, bg_scores, eroded_mix)

    erosion_model=nnlearn.ICA_PWM_model("erode", erosion_sources, test_model.informed_sources, test_model.source_length_limits,eroded_mix, erosion_lh, [""])

    eroded_model=nnlearn.erode_model(erosion_model, erosion_model.log_Li, obs, obsl, bg_scores, source_priors, 1000)
    @test eroded_model.log_Li > erosion_model.log_Li
    @test eroded_model.sources != erosion_model.sources
    @test eroded_model.mix_matrix == erosion_model.mix_matrix
    @test "EM from erode" in eroded_model.flags
    @test eroded_model.sources[1]==erosion_model.sources[1]
    @test eroded_model.sources[2]==erosion_model.sources[2]
    @test eroded_model.sources[3]!=erosion_model.sources[3]
    @test size(eroded_model.sources[3][1],1)==3

    path=randstring()
    test_record = nnlearn.Model_Record(path, rs_model.log_Li)
    serialize(path, rs_model)

    dm_model=nnlearn.distance_merge([test_record], test_model, test_model.log_Li, obs, obsl, bg_scores, source_priors, 1000)
    @test dm_model.log_Li > test_model.log_Li
    @test dm_model.sources != test_model.sources
    @test dm_model.mix_matrix != test_model.mix_matrix
    @test "FM from consolidate" in dm_model.flags
    #@test "DM from test" in dm_model.flags

    sm_model=nnlearn.similarity_merge([test_record], test_model, test_model.log_Li, obs, obsl, bg_scores, source_priors, 1000)
    @test sm_model.log_Li > test_model.log_Li
    @test sm_model.sources != test_model.sources
    @test sm_model.mix_matrix != test_model.mix_matrix
    #@test "SM from test" in sm_model.flags

    librarian=addprocs(1)
    @everywhere import nnlearn

    ddm_model=nnlearn.distance_merge(librarian[1], [test_record], test_model, test_model.log_Li, obs, obsl, bg_scores, source_priors, 1000)
    @test ddm_model.log_Li > test_model.log_Li
    @test ddm_model.sources != test_model.sources
    @test ddm_model.mix_matrix != test_model.mix_matrix
    @test ("DM from test" in ddm_model.flags) || ("FM from consolidate" in ddm_model.flags)

    dsm_model=nnlearn.similarity_merge(librarian[1], [test_record], rd_model, test_model.log_Li, obs, obsl, bg_scores, source_priors, 1000)
    @test dsm_model.log_Li > test_model.log_Li
    @test dsm_model.sources != test_model.sources
    @test dsm_model.mix_matrix != test_model.mix_matrix
    #@test "SM from test" in dsm_model.flags

    rmprocs(librarian)
    rm(path)
end

@testset "Ensemble assembly and nested sampling functions" begin
    ensembledir = randstring()
    spensembledir = randstring()
    distdir = randstring()

    source_pwm = [.7 .1 .1 .1
    .1 .1 .1 .7
    .1 .1 .7 .1]

    source_pwm_2 = [.6 .1 .1 .2
    .2 .1 .1 .6
    .1 .2 .6 .1]

    src_length_limits=2:12
    no_sources=4

    source_priors = nnlearn.assemble_source_priors(no_sources, [source_pwm, source_pwm_2], 4.0, src_length_limits)
    mix_prior=.5

    bg_scores = log.(fill(.1, (30,27)))
    obs=[BioSequences.DNASequence("CCGTTGACGATGTGATGAATAATGAAAGAA")
    BioSequences.DNASequence("CCCCGATGATGACCGTTGACCAGATGGATG")
    BioSequences.DNASequence("CCCCGATGATGACCCCGATTTTGAAAAAAA")
    BioSequences.DNASequence("TCATCATGCTGATGATGAATCAGATGAAAG")
    BioSequences.DNASequence("TGATGAATCTGACCCAGATGCCGATTTTGA")
    BioSequences.DNASequence("CCCCGATTTTGATCAGGATGAATAAAGAAA")
    BioSequences.DNASequence("TCATGGGCTGATGAACCGTTGACGATGAAA")
    BioSequences.DNASequence("TCATCCTGCTGACCCCGATTTCAGTGAAAA")
    BioSequences.DNASequence("TGATGAATAAAGTCATCCTGCATGTGAAAA")
    BioSequences.DNASequence("CCGTTGACGATGTGATGAATGATAAAGAAA")
    BioSequences.DNASequence("CCCCGATGATGACCGATGTTGACGATGAAA")
    BioSequences.DNASequence("CCCCGATGATGAATGCCCCGATTTTGAAAA")
    BioSequences.DNASequence("TCATCATGCTGATGATGAATAAAGAAAAAA")
    BioSequences.DNASequence("TGATGAATCTGACCCCGATCAGTTTGAAAA")
    BioSequences.DNASequence("CCCCGATTTTGATCAGATGGATGAATAAAG")
    BioSequences.DNASequence("TCATGGGCTGAACCGTTGACAGCGATGAAA")
    BioSequences.DNASequence("TCATCCTGCTCAGGACCCCGATTTTATGGA")
    BioSequences.DNASequence("TGATGAATCAGAAAGTCATCCTGCATGTGA")
    BioSequences.DNASequence("CCGTTGACCAGGATGTGATGAATAAAGAAA")
    BioSequences.DNASequence("CCCCGATGATGCAGACCGTTGACGATGAAA")
    BioSequences.DNASequence("CCCCGATGACAGTGACCCAGCCGATTTTGA")
    BioSequences.DNASequence("TCATCATGCTGAATGTGATGAATAAAAAAA")
    BioSequences.DNASequence("TGATGATGAATCTGAATGCCCCGATTTTGA")
    BioSequences.DNASequence("CCCCGATATGTTTGATGACAGTGAATAAAG")
    BioSequences.DNASequence("TCATGATGGGCTGAACCGTTGACGATGAAA")
    BioSequences.DNASequence("TCACAGTCCTGCTGACCCCGATTATGTTGA")
    BioSequences.DNASequence("TGATGATGAATAAAGTCATGATCCTGCTGA")
    ]

    order_seqs = BGHMM.get_order_n_seqs(obs, 0)
    coded_seqs = BGHMM.code_seqs(order_seqs)
    obs=Array(transpose(coded_seqs))
    position_start=1;offsets=[0,0]

    ensemble = nnlearn.Bayes_IPM_ensemble(ensembledir, 150, source_priors, (falses(0,0),mix_prior), bg_scores, obs, src_length_limits)
    ensemble = nnlearn.Bayes_IPM_ensemble(ensembledir, 200, source_priors, (falses(0,0),mix_prior), bg_scores, obs, src_length_limits) #test resumption

    sp_ensemble = nnlearn.Bayes_IPM_ensemble(spensembledir, 200, source_priors, (falses(0,0),mix_prior), bg_scores, obs, src_length_limits)

    @test length(ensemble.models) == 200
    for model in ensemble.models
        @test -1950 < model.log_Li < -1200
    end

    @test length(sp_ensemble.models) == 200
    for model in sp_ensemble.models
        @test -1950 < model.log_Li < -1200
    end

    assembler=addprocs(1)

    @everywhere using nnlearn

    dist_ensemble=nnlearn.Bayes_IPM_ensemble(assembler, distdir, 150, source_priors, (falses(0,0),mix_prior), bg_scores, obs, src_length_limits)
    dist_ensemble=nnlearn.Bayes_IPM_ensemble(assembler, distdir, 200, source_priors, (falses(0,0),mix_prior), bg_scores, obs, src_length_limits) #test resumption

    @test length(dist_ensemble.models) == 200
    for model in ensemble.models
        @test -1950 < model.log_Li < -1200
    end

    rmprocs(assembler)
    rm(distdir, recursive=true)

    permute_limit = 900
    job_sets=[
([
    ("PS", (no_sources)),
    ("PM", (no_sources)),
    ("PSFM", (no_sources)),
    ("PSFM", (no_sources, .8, 1.)),
    ("FM", ()),
    ("DM", (no_sources)),
    ("SM", (no_sources)),
    ("RD", (no_sources)),
    ("RI", (no_sources)),
    ("EM", (no_sources))
],[.025, .025, .025, .025, .775, .025, .025, .025, .025, .025]),
([
    ("PS", (no_sources)),
    ("PM", (no_sources)),
    ("PSFM", (no_sources)),
    ("PSFM", (no_sources, .8, 1.)),
    ("FM", ()),
    ("DM", (no_sources)),
    ("SM", (no_sources)),
    ("RD", (no_sources)),
    ("RI", (no_sources)),
    ("EM", (no_sources))
],[.15, .15, .15, .05, .05, 0.10, 0.10, 0.10, 0.05, .10]),
]
    job_limit=8
    job_set_thresh=[-Inf,ensemble.naive_lh]
    param_set=(job_sets,job_set_thresh,job_limit)

    @info "Testing threaded convergence..."
    sp_logZ = nnlearn.ns_converge!(sp_ensemble, param_set, permute_limit, 10., model_display=4)

    @test length(sp_ensemble.models) == 200
    @test length(sp_ensemble.log_Li) == length(sp_ensemble.log_Xi) == length(sp_ensemble.log_wi) == length(sp_ensemble.log_Liwi) == length(sp_ensemble.log_Zi) == length(sp_ensemble.Hi) == sp_ensemble.model_counter-200
    for i in 1:length(sp_ensemble.log_Li)-1
        @test sp_ensemble.log_Li[i] <= sp_ensemble.log_Li[i+1]
    end
    for i in 1:length(sp_ensemble.log_Zi)-1
        @test sp_ensemble.log_Zi[i] <= sp_ensemble.log_Zi[i+1]
    end
    @test sp_logZ > -1400.0

    @info "Testing multiprocess convergence..."
    @info "Spawning worker pool..."
    librarians=addprocs(1, topology=:master_worker)
    worker_pool=addprocs(2, topology=:master_worker)
    @everywhere using nnlearn,Random
    @everywhere Random.seed!(1)

    ####CONVERGE############
    final_logZ = nnlearn.ns_converge!(ensemble, [param_set,param_set], permute_limit, librarians, worker_pool, 10., backup=(true,250), model_display=4)

    rmprocs(worker_pool)
    rmprocs(librarians)

    @test length(ensemble.models) == 200
    @test length(ensemble.log_Li) == length(ensemble.log_Xi) == length(ensemble.log_wi) == length(ensemble.log_Liwi) == length(ensemble.log_Zi) == length(ensemble.Hi) == ensemble.model_counter-200
    for i in 1:length(ensemble.log_Li)-1
        @test ensemble.log_Li[i] <= ensemble.log_Li[i+1]
    end
    for i in 1:length(ensemble.log_Zi)-1
        @test ensemble.log_Zi[i] <= ensemble.log_Zi[i+1]
    end
    @test typeof(final_logZ) == Float64
    @test sp_logZ > -1400.0

    @info "Tests complete!"

    rm(ensembledir, recursive=true)
    rm(spensembledir, recursive=true)
end


text:
export getDiamondModel,getBallModel ;

function getDiamondModel(n::Array{Int64})
n_tup = tuple(n...);
n = collect(n_tup);
u = zeros(Float32,n_tup);
X,Y,Z = meshgrid(-div(n[1],2):div(n[1],2),-div(n[2],2):div(n[2],2),-div(n[3],2):div(n[3],2))
I = (abs.(X)+abs.(Y)+abs.(Z)) .<= div(n[1],4);
u[I] .= 1.0;
return u;
end

function getBallModel(n::Array{Int64})
n_tup = tuple(n...);
n = collect(n_tup);
u = zeros(Float32,n_tup);
X,Y,Z = meshgrid(-div(n[1],2):div(n[1],2),-div(n[2],2):div(n[2],2),-div(n[3],2):div(n[3],2))
I = (sqrt.(X.^2+Y.^2 + Z.^2)) .<= div(n[1],4);
u[I] .= 1.0;
return u;
end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
using Test
using StaticArrays
using LinearAlgebra
using TrajOptPlots
using MeshCat
using GeometryTypes
const TO = TrajectoryOptimization

model_qv = Dynamics.FreeBody{UnitQuaternion{Float64,VectorPart},Float64}()
model_qe = Dynamics.FreeBody{UnitQuaternion{Float64,ExponentialMap},Float64}()
model_qp = Dynamics.FreeBody{UnitQuaternion{Float64,MRPMap},Float64}()
model_qg = Dynamics.FreeBody{UnitQuaternion{Float64,CayleyMap},Float64}()
model_p = Dynamics.FreeBody{MRP{Float64},Float64}()
model_e = Dynamics.FreeBody{RPY{Float64},Float64}()

x,u = rand(model_p)

model_p isa Dynamics.FreeBody{<:Rotation}
size(model_e) == (12,6)
size(model_qv) == (13,6)

# Set up Trajectory Optimization Problem
model = model_qv

# Params
N = 101
tf = 10.
x0 = zeros(model)[1]

r0 = @SVector zeros(3)
q0 = zero(UnitQuaternion)
v0 = @SVector zeros(3)
ω0 = @SVector zeros(3)
x0 = Dynamics.build_state(model, r0, q0, v0, ω0)
@test x0[4] == 1

# Objective
Q = Diagonal(@SVector fill(0.1, 12))
R = Diagonal(@SVector fill(0.05, 6))
Qf = Diagonal(@SVector fill(10.0, 12))

rf = 0*@SVector [10,10,10]
qf = UnitQuaternion(RPY(pi/4, 0, 0))
xf = Dynamics.build_state(model, rf, qf, v0, ω0)
Gf = Dynamics.state_diff_jacobian(model, xf)

Q = Gf*Q*Gf'
Qf = Gf*Qf*Gf'
obj = LQRObjective(Q, R, Qf, xf, N, checks=false)

# Build Problem
prob = Problem(model, obj, xf, tf, x0=x0)

# Solve the problem
ilqr = iLQRSolver(prob)
ilqr.opts.verbose = true
U0 = [@SVector zeros(6) for k = 1:N-1]
initial_controls!(ilqr, U0)
solve!(ilqr)
q = states(ilqr)[N][4:7]
q'xf[4:7]

vis = Visualizer(); open(vis);
dims = [1,1,3]
robot = HyperRectangle(Vec((-1 .* dims ./2)...), Vec(dims...)); setobject!(vis["robot"], robot);
visualize!(vis, model, get_trajectory(ilqr))

Z = get_trajectory(ilqr)
states(ilqr)[end]
TO.stage_cost(obj[N], Z[end])
controls(ilqr)
ilqr.∇F
ilqr.K

# Step through the solve
solver = ilqr
initial_controls!(solver, U0)
rollout!(solver)
cost(solver)
Z = solver.Z
TO.state_diff_jacobian!(solver.G, solver.model, Z)

solver.G
controls(ilqr)
states(ilqr)[end]

# Visualizer
using MeshCat, GeometryTypes, CoordinateTransformations
using TrajOptPlots
vis = Visualizer()
open(vis)
robot = HyperRectangle(Vec((-1 .* dims ./2)...), Vec(dims...))
setobject!(vis["robot"], robot)

visualize!(vis, model, get_trajectory(ilqr))
Z = get_trajectory(ilqr)
X = states(Z)
states(ilqr)[2]
q = orientation(model, X[end])
settransform!(vis["robot"], LinearMap(Quat(q)))


text:
mutable struct TimeStepWizard{T}
              cfl :: T
    diffusive_cfl :: T
       max_change :: T
       min_change :: T
           max_Δt :: T
           min_Δt :: T
               Δt :: T
end

"""
    TimeStepWizard(cfl=0.1, max_change=2.0, min_change=0.5, max_Δt=Inf, min_Δt=0.0, Δt=0.01)

A type for calculating adaptive time steps based on capping the CFL number at `cfl`.

On calling `update_Δt!(wizard, model)`, the `TimeStepWizard` computes a time-step such that
``cfl = max(u/Δx, v/Δy, w/Δz) Δt``, where ``max(u/Δx, v/Δy, w/Δz)`` is the maximum ratio
between model velocity and along-velocity grid spacing anywhere on the model grid. The new
`Δt` is constrained to change by a multiplicative factor no more than `max_change` or no
less than `min_change` from the previous `Δt`, and to be no greater in absolute magnitude
than `max_Δt` and no less than `min_Δt`.
"""
TimeStepWizard(; cfl=0.1, diffusive_cfl=Inf, max_change=2.0, min_change=0.5, max_Δt=Inf, min_Δt=0.0, Δt=0.01) =
        TimeStepWizard{typeof(Δt)}(cfl, diffusive_cfl, max_change, min_change, max_Δt, min_Δt, Δt)

"""
    update_Δt!(wizard, model)

Compute `wizard.Δt` given the velocities and diffusivities of `model`, and the parameters
of `wizard`.
"""
function update_Δt!(wizard, model)

    Δt = min(
             wizard.cfl * cell_advection_timescale(model),          # advective
             wizard.diffusive_cfl * cell_diffusion_timescale(model) # diffusive
            )

    # Put the kibosh on if needed
    Δt = min(wizard.max_change * wizard.Δt, Δt)
    Δt = max(wizard.min_change * wizard.Δt, Δt)
    Δt = clamp(Δt, wizard.min_Δt, wizard.max_Δt)

    wizard.Δt = Δt

    return nothing
end

(c::CFL{<:TimeStepWizard})(model) = c.Δt.Δt / c.timescale(model)


text:
import MetNets: av, va

# MetState interface
av(out::EPOut) = out.av
va(out::EPOut) = out.va

# EPOut
μ(out::EPOut) = out.μ
σ(out::EPOut) = out.σ


text:
# Julia wrapper for header: /usr/local/ffmpeg/2.4/include/libavutil/avutil.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    avutil_version,
    avutil_configuration,
    avutil_license,
    av_get_media_type_string,
    av_get_picture_type_char,
    av_x_if_null,
    av_int_list_length_for_size,
    av_fopen_utf8,
    av_get_time_base_q

function avutil_version()
    ccall((:avutil_version,libavutil),Uint32,())
end

function avutil_configuration()
    ccall((:avutil_configuration,libavutil),Ptr{Uint8},())
end

function avutil_license()
    ccall((:avutil_license,libavutil),Ptr{Uint8},())
end

function av_get_media_type_string(media_type::AVMediaType)
    ccall((:av_get_media_type_string,libavutil),Ptr{Uint8},(AVMediaType,),media_type)
end

function av_get_picture_type_char(pict_type::AVPictureType)
    ccall((:av_get_picture_type_char,libavutil),Uint8,(AVPictureType,),pict_type)
end

function av_x_if_null(p,x)
    ccall((:av_x_if_null,libavutil),Ptr{Void},(Ptr{Void},Ptr{Void}),p,x)
end

function av_int_list_length_for_size(elsize::Integer,list,term::Uint64)
    ccall((:av_int_list_length_for_size,libavutil),Uint32,(Uint32,Ptr{Void},Uint64),elsize,list,term)
end

function av_fopen_utf8(path,mode)
    ccall((:av_fopen_utf8,libavutil),Ptr{FILE},(Ptr{Uint8},Ptr{Uint8}),path,mode)
end

function av_get_time_base_q()
    ccall((:av_get_time_base_q,libavutil),AVRational,())
end


text:
# This file was generated, do not modify it. # hide
X = coerce(X, autotype(X, rules=(:discrete_to_continuous,)))

dtr_model = DTR()
dtr = machine(dtr_model, X, y)

fit!(dtr, rows=train)

ypred = MLJ.predict(dtr, rows=test)
round(rms(ypred, y[test]), sigdigits=3)

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
module BirdClef

include("allmodules.jl")

using .Data

const dataset = "birdclef-2022"

function main(data_path=joinpath(@__DIR__,"data"))
    @info "Looking for dataset using " data_path dataset
    if !has_dataset(data_path,dataset)
        @warn "You do not have the required kaggle dataset."
        @info "Attemping to download kaggle dataset"
        download_dataset(data_path,dataset)
    else
        @info "Found dataset!"
    end
end

export main,Data

end

max_stars_repo_name: BirdClef/birdclef

max_issues_repo_name: BirdClef/birdclef

max_forks_repo_name: BirdClef/birdclef


text:
module Arblib

using Arb_jll
import LinearAlgebra

if VERSION >= v"1.5.0-DEV.639"
    import Base: contains
end

export Mag,
    MagRef,
    Arf,
    ArfRef,
    Arb,
    ArbRef,
    Acb,
    AcbRef,
    ArbVector,
    AcbVector,
    ArbMatrix,
    AcbMatrix,
    ArbRefVector,
    AcbRefVector,
    ArbRefMatrix,
    AcbRefMatrix,
    ArbPoly,
    ArbSeries,
    AcbPoly,
    AcbSeries,
    ref

macro libarb(function_name)
    return (:($function_name), libarb)
end

macro libflint(function_name)
    return (:($function_name), Arb_jll.libflint)
end

include("arb_types.jl")
include("rounding.jl")
include("types.jl")
include("arbcall.jl")
include("manual_overrides.jl")

include("precision.jl")

include("setters.jl")
include("constructors.jl")
include("predicates.jl")
include("show.jl")
include("arithmetic.jl")
include("rand.jl")

include("ref.jl")
include("vector.jl")
include("matrix.jl")
include("array_common.jl")
include("eigen.jl")
include("poly.jl")
include("calc_integrate.jl")

include("arbcalls/mag.jl")
include("arbcalls/arf.jl")
include("arbcalls/arb.jl")
include("arbcalls/acb.jl")
include("arbcalls/arb_poly.jl")
include("arbcalls/acb_poly.jl")
include("arbcalls/arb_mat.jl")
include("arbcalls/acb_mat.jl")
include("arbcalls/acb_hypgeom.jl")
include("arbcalls/arb_hypgeom.jl")
include("arbcalls/acb_elliptic.jl")
include("arbcalls/acb_dirichlet.jl")
include("arbcalls/eigen.jl")

end # module


text:
module TestBasic

using ConicHulls, ConicHulls.Common, ConicHulls.Dets, ConicHulls.Hulls
using ConicHulls.RefHull
using ConicHulls.Hulls.dominates
using ConicHulls.Insertion.find_dominated_facet

function printhull(hull::ConicHull)
    println("ConicHull:")
    for facet in hull.facets
        println("  ", facet)
    end
end

function test(NC)
    initial_gs = eye(Int, NC)
    hull = create_hull(NC, initial_gs)
    NC, F, G = nconic(hull), ftype(hull), gtype(hull)

    gs, fs = Hulls.create_simplex(NC, F, G, initial_gs)
    # @show gs
    # @show fs

    @assert det(gs...) > 0
    @assert det(gs[[2, 1, (3:NC)...]]...) == -det(gs...)
    for k=1:NC
        @assert !dominates(gs[1], fs[1])
    end

    gs = copy(hull.generators)
    verify(hull)
    verify_hull(hull, gs)

    for facet in hull.facets
        x = sum([g.x for g in facet.generators])
        # Check that conic combination of generators is in the hull
        g_in = G(x) # todo: better / safer way to create new generator?
        @assert find_dominated_facet(hull, g_in) === nothing
        # Check that conic combination of the facet generators minus
        # the opposing generator dominates the facet (and no other)
        opposite_g = first(setdiff(Set(hull.generators), Set(facet.generators)))
        x_out = x - opposite_g.x
        g_out = G(x_out)
        # @show x_out
        @assert find_dominated_facet(hull, g_out) === facet
    end

    x = [-1, fill(1, NC-1)...]
#    printhull(hull)
#    println(g)
    g = add!(hull, x)
    @assert g != nothing
    verify(hull)

    push!(gs, g)
    verify_hull(hull, gs)
#    printhull(hull)
end

for k=1:2
    for NC in 3:5
        @time test(NC)
    end
end

end # module


text:
using Colors
using Images
using Boltzmann
using MNIST
using ImageView
using Base.Test

function run_mnist()
    X, y = testdata()  # test data is smaller, no need to downsample
    normalize_samples!(X)
    binarize!(X;threshold=0.01)

    TrainSet = X[:,1:9000]
    ValidSet = X[:,9001:end]
    HiddenUnits = 100;
    Epochs = 10;
    Gibbs = 1;
    LearnRate = 0.05
    MonitorEvery=1

    rbm = BernoulliRBM(28*28, HiddenUnits, (28,28); momentum=0.0)

    # Attempt with L2 weight decay
    finalrbm,monitor = fit(rbm, TrainSet;n_iter=Epochs,
                          weight_decay="l2",
                          decay_magnitude=0.001,
                          lr=LearnRate,
                          validation=ValidSet,
                          NormalizationApproxIter=Gibbs,
                          monitor_every=MonitorEvery,
                          monitor_vis=true)

    info("==========monitor==============")
    println(monitor)
    info("===============================")

    WriteMonitorChartPDF(finalrbm,monitor,X,"testmonitor.pdf")
    SaveMonitorHDF5(monitor,"testmonitor.h5")
end

run_mnist()


text:
#=
MathOptInterface wrapper of Hypatia solver
=#

const SV = MOI.SingleVariable
const SAF = MOI.ScalarAffineFunction
const VV = MOI.VectorOfVariables
const VAF = MOI.VectorAffineFunction

export Optimizer

"""
$(TYPEDEF)

A MathOptInterface optimizer type for Hypatia.
"""
mutable struct Optimizer{T <: Real} <: MOI.AbstractOptimizer
    use_dense_model::Bool # make the model use dense A and G data instead of sparse

    solver::Solvers.Solver{T} # Hypatia solver object
    model::Models.Model{T} # Hypatia model object

    # result data
    x::Vector{T}
    s::Vector{T}
    y::Vector{T}
    z::Vector{T}

    # data for transforming certificates
    obj_sense::MOI.OptimizationSense
    num_eq_constrs::Int
    constr_offset_eq::Vector{Int}
    constr_prim_eq::Vector{T}
    constr_offset_cone::Vector{Int}
    constr_prim_cone::Vector{T}
    nonpos_idxs::UnitRange{Int}
    interval_idxs::UnitRange{Int}
    interval_scales::Vector{T}
    moi_other_cones_start::Int
    moi_other_cones::Vector{MOI.AbstractVectorSet}

    function Optimizer{T}(;
        use_dense_model::Bool = true, # TODO should depend on the size and sparsity of A, G in the model
        solver_options... # TODO allow passing in a solver?
        ) where {T <: Real}
        opt = new{T}()
        opt.use_dense_model = use_dense_model
        if !haskey(solver_options, :syssolver)
            # choose default system solver based on use_dense_model
            sstype = (use_dense_model ? Solvers.QRCholDenseSystemSolver :
                Solvers.SymIndefSparseSystemSolver)
            solver_options = (solver_options..., syssolver = sstype{T}())
        end
        if !haskey(solver_options, :preprocess)
            # only preprocess if using dense model # TODO maybe should preprocess if sparse
            solver_options = (solver_options..., preprocess = use_dense_model)
        end
        if !haskey(solver_options, :reduce)
            # only reduce if using dense model
            solver_options = (solver_options..., reduce = use_dense_model)
        end
        if !haskey(solver_options, :init_use_indirect)
            # only use indirect if not using dense model
            solver_options = (solver_options...,
                init_use_indirect = !use_dense_model)
        end
        opt.solver = Solvers.Solver{T}(; solver_options...)
        return opt
    end
end

Optimizer(; options...) = Optimizer{Float64}(; options...) # default to Float64

MOI.get(::Optimizer, ::MOI.SolverName) = "Hypatia"
MOI.get(opt::Optimizer, ::MOI.RawSolver) = opt.solver

MOI.is_empty(opt::Optimizer) = (opt.solver.status == Solvers.NotLoaded)

MOI.empty!(opt::Optimizer) = (opt.solver.status = Solvers.NotLoaded)

MOI.supports(
    ::Optimizer{T},
    ::Union{
        MOI.ObjectiveSense,
        MOI.ObjectiveFunction{SV},
        MOI.ObjectiveFunction{SAF{T}},
        },
    ) where {T <: Real} = true

MOI.supports_constraint(
    ::Optimizer{T},
    ::Type{<:Union{SV, SAF{T}}},
    ::Type{<:Union{
        MOI.EqualTo{T},
        MOI.GreaterThan{T},
        MOI.LessThan{T},
        MOI.Interval{T},
        }},
    ) where {T <: Real} = true

MOI.supports_constraint(
    ::Optimizer{T},
    ::Type{<:Union{VV, VAF{T}}},
    ::Type{<:SupportedCones{T}},
    ) where {T <: Real} = true

# build representation as min c'x s.t. A*x = b, h - G*x in K
function MOI.copy_to(
    opt::Optimizer{T},
    src::MOI.ModelLike;
    copy_names::Bool = false,
    ) where {T <: Real}
    @assert !copy_names
    idx_map = MOI.Utilities.IndexMap()

    # variables
    n = MOI.get(src, MOI.NumberOfVariables()) # columns of A
    j = 0
    for vj in MOI.get(src, MOI.ListOfVariableIndices()) # MOI.VariableIndex
        j += 1
        idx_map[vj] = MOI.VariableIndex(j)
    end
    @assert j == n

    # objective function
    F = MOI.get(src, MOI.ObjectiveFunctionType())
    if F == SV
        obj = SAF{T}(MOI.get(src, MOI.ObjectiveFunction{F}()))
    elseif F == SAF{T}
        obj = MOI.get(src, MOI.ObjectiveFunction{F}())
    else
        error("function type $F not supported")
    end
    (Jc, Vc) = (Int[], T[])
    for t in obj.terms
        push!(Jc, idx_map[t.variable_index].value)
        push!(Vc, t.coefficient)
    end
    obj_offset = obj.constant
    if MOI.get(src, MOI.ObjectiveSense()) == MOI.MAX_SENSE
        Vc .*= -1
        obj_offset *= -1
    end
    opt.obj_sense = MOI.get(src, MOI.ObjectiveSense())
    model_c = Vector(sparsevec(Jc, Vc, n))

    # constraints
    get_src_cons(F, S) = MOI.get(src, MOI.ListOfConstraintIndices{F, S}())
    get_con_fun(con_idx) = MOI.get(src, MOI.ConstraintFunction(), con_idx)
    get_con_set(con_idx) = MOI.get(src, MOI.ConstraintSet(), con_idx)
    i = 0 # MOI constraint index

    # equality constraints
    p = 0 # rows of A (equality constraint matrix)
    (IA, JA, VA) = (Int[], Int[], T[])
    (Ib, Vb) = (Int[], T[])
    (Icpe, Vcpe) = (Int[], T[]) # constraint set constants for opt.constr_prim_eq
    constr_offset_eq = Vector{Int}()

    for F in (SV, SAF{T}), ci in get_src_cons(F, MOI.EqualTo{T})
        i += 1
        idx_map[ci] = MOI.ConstraintIndex{F, MOI.EqualTo{T}}(i)
        push!(constr_offset_eq, p)
        p += 1
        fi = get_con_fun(ci)
        si = get_con_set(ci)
        if F == SV
            push!(IA, p)
            push!(JA, idx_map[fi.variable].value)
            push!(VA, -1)
            push!(Vb, -si.value)
        else
            for vt in fi.terms
                push!(IA, p)
                push!(JA, idx_map[vt.variable_index].value)
                push!(VA, -vt.coefficient)
            end
            push!(Vb, fi.constant - si.value)
        end
        push!(Ib, p)
        push!(Icpe, p)
        push!(Vcpe, si.value)
    end

    for F in (VV, VAF{T}), ci in get_src_cons(F, MOI.Zeros)
        i += 1
        idx_map[ci] = MOI.ConstraintIndex{F, MOI.Zeros}(i)
        push!(constr_offset_eq, p)
        fi = get_con_fun(ci)
        dim = MOI.output_dimension(fi)
        if F == VV
            append!(IA, (p + 1):(p + dim))
            append!(JA, idx_map[vi].value for vi in fi.variables)
            append!(VA, -ones(T, dim))
        else
            for vt in fi.terms
                push!(IA, p + vt.output_index)
                push!(JA, idx_map[vt.scalar_term.variable_index].value)
                push!(VA, -vt.scalar_term.coefficient)
            end
            append!(Ib, (p + 1):(p + dim))
            append!(Vb, fi.constants)
        end
        p += dim
    end

    push!(constr_offset_eq, p)
    model_A = dropzeros!(sparse(IA, JA, VA, p, n))
    model_b = Vector(sparsevec(Ib, Vb, p))
    opt.num_eq_constrs = i
    opt.constr_prim_eq = Vector(sparsevec(Icpe, Vcpe, p))
    opt.constr_offset_eq = constr_offset_eq

    # conic constraints
    q = 0 # rows of G (cone constraint matrix)
    (IG, JG, VG) = (Int[], Int[], T[])
    (Ih, Vh) = (Int[], T[])
    (Icpc, Vcpc) = (Int[], T[]) # constraint set constants for opt.constr_prim_eq
    constr_offset_cone = Vector{Int}()
    cones = Cones.Cone{T}[]

    # build up one nonnegative cone from LP constraints
    nonneg_start = q

    # nonnegative-like constraints
    for F in (SV, SAF{T}), ci in get_src_cons(F, MOI.GreaterThan{T})
        i += 1
        idx_map[ci] = MOI.ConstraintIndex{F, MOI.GreaterThan{T}}(i)
        push!(constr_offset_cone, q)
        q += 1
        fi = get_con_fun(ci)
        si = get_con_set(ci)
        if F == SV
            push!(IG, q)
            push!(JG, idx_map[fi.variable].value)
            push!(VG, -1)
            push!(Vh, -si.lower)
        else
            for vt in fi.terms
                push!(IG, q)
                push!(JG, idx_map[vt.variable_index].value)
                push!(VG, -vt.coefficient)
            end
            push!(Vh, fi.constant - si.lower)
        end
        push!(Ih, q)
        push!(Vcpc, si.lower)
        push!(Icpc, q)
    end

    for F in (VV, VAF{T}), ci in get_src_cons(F, MOI.Nonnegatives)
        i += 1
        idx_map[ci] = MOI.ConstraintIndex{F, MOI.Nonnegatives}(i)
        push!(constr_offset_cone, q)
        fi = get_con_fun(ci)
        if F == VV
            for vj in fi.variables
                q += 1
                push!(IG, q)
                push!(JG, idx_map[vj].value)
                push!(VG, -1)
            end
        else
            dim = MOI.output_dimension(fi)
            for vt in fi.terms
                push!(IG, q + vt.output_index)
                push!(JG, idx_map[vt.scalar_term.variable_index].value)
                push!(VG, -vt.scalar_term.coefficient)
            end
            append!(Ih, (q + 1):(q + dim))
            append!(Vh, fi.constants)
            q += dim
        end
    end

    # nonpositive-like constraints
    nonpos_start = q

    for F in (SV, SAF{T}), ci in get_src_cons(F, MOI.LessThan{T})
        i += 1
        idx_map[ci] = MOI.ConstraintIndex{F, MOI.LessThan{T}}(i)
        push!(constr_offset_cone, q)
        q += 1
        fi = get_con_fun(ci)
        si = get_con_set(ci)
        if F == SV
            push!(IG, q)
            push!(JG, idx_map[fi.variable].value)
            push!(VG, 1)
            push!(Vh, si.upper)
        else
            for vt in fi.terms
                push!(IG, q)
                push!(JG, idx_map[vt.variable_index].value)
                push!(VG, vt.coefficient)
            end
            push!(Vh, -fi.constant + si.upper)
        end
        push!(Ih, q)
        push!(Vcpc, si.upper)
        push!(Icpc, q)
    end

    for F in (VV, VAF{T}), ci in get_src_cons(F, MOI.Nonpositives)
        i += 1
        idx_map[ci] = MOI.ConstraintIndex{F, MOI.Nonpositives}(i)
        push!(constr_offset_cone, q)
        fi = get_con_fun(ci)
        if F == VV
            for vj in fi.variables
                q += 1
                push!(IG, q)
                push!(JG, idx_map[vj].value)
                push!(VG, 1)
            end
        else
            dim = MOI.output_dimension(fi)
            for vt in fi.terms
                push!(IG, q + vt.output_index)
                push!(JG, idx_map[vt.scalar_term.variable_index].value)
                push!(VG, vt.scalar_term.coefficient)
            end
            append!(Ih, (q + 1):(q + dim))
            append!(Vh, -fi.constants)
            q += dim
        end
    end

    # single nonnegative cone
    opt.nonpos_idxs = (nonpos_start + 1):q
    if q > nonneg_start
        push!(cones, Cones.Nonnegative{T}(q - nonneg_start))
    end

    # build up one L_infinity norm cone from two-sided interval constraints
    interval_start = q
    SV_intvl = get_src_cons(SV, MOI.Interval{T})
    SAF_intvl = get_src_cons(SAF{T}, MOI.Interval{T})
    num_intervals = length(SV_intvl) + length(SAF_intvl)
    interval_scales = zeros(T, num_intervals)

    if num_intervals > 0
        i += 1
        push!(constr_offset_cone, q)
        q += 1
        push!(Ih, q)
        push!(Vh, 1)
    end

    interval_count = 0
    for (F, Cs) in ((SV, SV_intvl), (SAF{T}, SAF_intvl)), ci in Cs
        i += 1
        idx_map[ci] = MOI.ConstraintIndex{F, MOI.Interval{T}}(i)
        push!(constr_offset_cone, q)
        q += 1

        upper = get_con_set(ci).upper
        lower = get_con_set(ci).lower
        @assert isfinite(upper) && isfinite(lower)
        @assert upper > lower
        mid = (upper + lower) / T(2)
        scal = 2 / (upper - lower)

        fi = get_con_fun(ci)
        if F == SV
            push!(IG, q)
            push!(JG, idx_map[fi.variable].value)
            push!(VG, -scal)
            push!(Vh, -mid * scal)
        else
            for vt in fi.terms
                push!(IG, q)
                push!(JG, idx_map[vt.variable_index].value)
                push!(VG, -vt.coefficient * scal)
            end
            push!(Vh, (fi.constant - mid) * scal)
        end
        push!(Ih, q)
        push!(Vcpc, mid)
        push!(Icpc, q)
        interval_count += 1
        interval_scales[interval_count] = scal
    end

    opt.interval_idxs = (interval_start + 2):q
    opt.interval_scales = interval_scales
    if q > interval_start
        # exists at least one interval-type constraint
        push!(cones, Cones.EpiNormInf{T, T}(q - interval_start))
    end

    # non-LP conic constraints
    moi_other_cones_start = i + 1
    moi_other_cones = MOI.AbstractVectorSet[]

    for (F, S) in MOI.get(src, MOI.ListOfConstraints())
        if S <: LinearCones{T}
            continue # already copied these constraints
        end
        @assert S <: SupportedCones{T}

        for ci in get_src_cons(F, S)
            i += 1
            idx_map[ci] = MOI.ConstraintIndex{F, S}(i)
            push!(constr_offset_cone, q)
            fi = get_con_fun(ci)
            si = get_con_set(ci)
            push!(moi_other_cones, si)
            dim = MOI.output_dimension(fi)

            if F == VV
                JGi = (idx_map[vj].value for vj in fi.variables)
                IGi = permute_affine(si, 1:dim)
                VGi = rescale_affine(si, fill(-one(T), dim))
            else
                JGi = (idx_map[vt.scalar_term.variable_index].value
                    for vt in fi.terms)
                IGi = permute_affine(si, [vt.output_index for vt in fi.terms])
                VGi = rescale_affine(si, [-vt.scalar_term.coefficient
                    for vt in fi.terms], IGi)
                Ihi = permute_affine(si, 1:dim)
                Vhi = rescale_affine(si, fi.constants)
                Ihi = q .+ Ihi
                append!(Ih, Ihi)
                append!(Vh, Vhi)
            end

            IGi = q .+ IGi
            append!(IG, IGi)
            append!(JG, JGi)
            append!(VG, VGi)
            push!(cones, cone_from_moi(T, si))
            q += dim
        end
    end

    push!(constr_offset_cone, q)

    # finalize model
    model_G = dropzeros!(sparse(IG, JG, VG, q, n))
    model_h = Vector(sparsevec(Ih, Vh, q))

    opt.model = Models.Model{T}(model_c, model_A, model_b, model_G, model_h,
        cones; obj_offset = obj_offset)
    if opt.use_dense_model # convert A and G to dense
        Models.densify!(opt.model)
    end

    opt.constr_offset_cone = constr_offset_cone
    opt.constr_prim_cone = Vector(sparsevec(Icpc, Vcpc, q))
    opt.moi_other_cones_start = moi_other_cones_start
    opt.moi_other_cones = moi_other_cones

    return idx_map
end

function MOI.optimize!(opt::Optimizer{T}) where {T <: Real}
    # build and solve the model
    model = opt.model
    solver = opt.solver

    Solvers.load(solver, model)
    Solvers.solve(solver)

    status = Solvers.get_status(solver)
    primal_obj = Solvers.get_primal_obj(solver)
    dual_obj = Solvers.get_dual_obj(solver)
    opt.x = Solvers.get_x(solver)
    opt.y = Solvers.get_y(solver)
    opt.s = Solvers.get_s(solver)
    opt.z = Solvers.get_z(solver)

    # transform solution for MOI conventions
    opt.constr_prim_eq += model.b - model.A * opt.x
    opt.s[opt.nonpos_idxs] .*= -1
    opt.z[opt.nonpos_idxs] .*= -1
    opt.s[opt.interval_idxs] ./= opt.interval_scales
    i = opt.moi_other_cones_start - opt.num_eq_constrs
    for cone in opt.moi_other_cones
        if needs_untransform(cone)
            os = opt.constr_offset_cone
            idxs = (os[i] + 1):os[i + 1]
            @views untransform_affine(cone, opt.s[idxs])
            @views untransform_affine(cone, opt.z[idxs])
        end
        i += 1
    end
    opt.constr_prim_cone .+= opt.s
    opt.z[opt.interval_idxs] .*= opt.interval_scales

    return
end

MOI.supports(::Optimizer, ::MOI.Silent) = true
MOI.set(opt::Optimizer, ::MOI.Silent, value::Bool) = (opt.solver.verbose = !value)
MOI.get(opt::Optimizer, ::MOI.Silent) = !opt.solver.verbose

MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
MOI.set(opt::Optimizer, ::MOI.TimeLimitSec, value::Real) =
    (opt.solver.time_limit = value)
MOI.set(opt::Optimizer, ::MOI.TimeLimitSec, ::Nothing) =
    (opt.solver.time_limit = Inf)

function MOI.get(opt::Optimizer, ::MOI.TimeLimitSec)
    isfinite(opt.solver.time_limit) && return opt.solver.time_limit
    return
end

function MOI.get(opt::Optimizer, ::MOI.SolveTime)
    if opt.solver.status in (:NotLoaded, :Loaded)
        error("solve has not been called")
    end
    return Solvers.get_solve_time(opt.solver)
end

MOI.get(opt::Optimizer, ::MOI.RawStatusString) = string(opt.solver.status)

MOI.get(opt::Optimizer, ::MOI.BarrierIterations) = opt.solver.num_iters

function MOI.get(opt::Optimizer, ::MOI.TerminationStatus)
    status = opt.solver.status
    if status in (Solvers.NotLoaded, Solvers.Loaded)
        return MOI.OPTIMIZE_NOT_CALLED
    elseif status == Solvers.Optimal
        return MOI.OPTIMAL
    elseif status == Solvers.PrimalInfeasible
        return MOI.INFEASIBLE
    elseif status == Solvers.DualInfeasible
        return MOI.DUAL_INFEASIBLE
    elseif status == Solvers.SlowProgress
        return MOI.SLOW_PROGRESS
    elseif status == Solvers.IterationLimit
        return MOI.ITERATION_LIMIT
    elseif status == Solvers.TimeLimit
        return MOI.TIME_LIMIT
    else
        @warn("Hypatia status $(opt.solver.status) not handled")
        return MOI.OTHER_ERROR
    end
end

function MOI.get(opt::Optimizer, ::MOI.PrimalStatus)
    status = opt.solver.status
    if status == Solvers.Optimal
        return MOI.FEASIBLE_POINT
    elseif status == Solvers.PrimalInfeasible
        return MOI.INFEASIBLE_POINT
    elseif status == Solvers.DualInfeasible
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif status == Solvers.IllPosed
        return MOI.OTHER_RESULT_STATUS
    else
        return MOI.UNKNOWN_RESULT_STATUS
    end
end

function MOI.get(opt::Optimizer, ::MOI.DualStatus)
    status = opt.solver.status
    if status == Solvers.Optimal
        return MOI.FEASIBLE_POINT
    elseif status == Solvers.PrimalInfeasible
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif status == Solvers.DualInfeasible
        return MOI.INFEASIBLE_POINT
    elseif status == Solvers.IllPosed
        return MOI.OTHER_RESULT_STATUS
    else
        return MOI.UNKNOWN_RESULT_STATUS
    end
end

function MOI.get(opt::Optimizer, ::MOI.ObjectiveValue)
    raw_obj_val = Solvers.get_primal_obj(opt.solver)
    return ((opt.obj_sense == MOI.MAX_SENSE) ? -1 : 1) * raw_obj_val
end

function MOI.get(
    opt::Optimizer,
    ::Union{MOI.DualObjectiveValue, MOI.ObjectiveBound},
    )
    raw_dual_obj_val = Solvers.get_dual_obj(opt.solver)
    return ((opt.obj_sense == MOI.MAX_SENSE) ? -1 : 1) * raw_dual_obj_val
end

MOI.get(opt::Optimizer, ::MOI.ResultCount) = 1

MOI.get(
    opt::Optimizer,
    ::MOI.VariablePrimal,
    vi::MOI.VariableIndex,
    ) = opt.x[vi.value]

MOI.get(
    opt::Optimizer,
    a::MOI.VariablePrimal,
    vi::Vector{MOI.VariableIndex},
    ) = MOI.get.(opt, a, vi)

function MOI.get(
    opt::Optimizer,
    ::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{F, S},
    ) where {F <: MOI.AbstractFunction, S <: MOI.AbstractScalarSet}
    # scalar set
    i = ci.value
    if i <= opt.num_eq_constrs
        # constraint is an equality
        return opt.y[opt.constr_offset_eq[i] + 1]
    else
        # constraint is conic
        i -= opt.num_eq_constrs
        return opt.z[opt.constr_offset_cone[i] + 1]
    end
end

function MOI.get(
    opt::Optimizer,
    ::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{F, S},
    ) where {F <: MOI.AbstractFunction, S <: MOI.AbstractVectorSet}
    # vector set
    i = ci.value
    if i <= opt.num_eq_constrs
        # constraint is an equality
        os = opt.constr_offset_eq
        return opt.y[(os[i] + 1):os[i + 1]]
    else
        # constraint is conic
        i -= opt.num_eq_constrs
        os = opt.constr_offset_cone
        return opt.z[(os[i] + 1):os[i + 1]]
    end
end

MOI.get(
    opt::Optimizer,
    a::MOI.ConstraintDual,
    ci::Vector{MOI.ConstraintIndex},
    ) = MOI.get.(opt, a, ci)

function MOI.get(
    opt::Optimizer,
    ::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex{F, S},
    ) where {F <: MOI.AbstractFunction, S <: MOI.AbstractScalarSet}
    # scalar set
    i = ci.value
    if i <= opt.num_eq_constrs
        # constraint is an equality
        return opt.constr_prim_eq[opt.constr_offset_eq[i] + 1]
    else
        # constraint is conic
        i -= opt.num_eq_constrs
        return opt.constr_prim_cone[opt.constr_offset_cone[i] + 1]
    end
end

function MOI.get(
    opt::Optimizer,
    ::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex{F, S},
    ) where {F <: MOI.AbstractFunction, S <: MOI.AbstractVectorSet}
    # vector set
    i = ci.value
    if i <= opt.num_eq_constrs
        # constraint is an equality
        os = opt.constr_offset_eq
        return opt.constr_prim_eq[(os[i] + 1):os[i + 1]]
    else
        # constraint is conic
        i -= opt.num_eq_constrs
        os = opt.constr_offset_cone
        return opt.constr_prim_cone[(os[i] + 1):os[i + 1]]
    end
end

MOI.get(
    opt::Optimizer,
    a::MOI.ConstraintPrimal,
    ci::Vector{MOI.ConstraintIndex},
    ) = MOI.get.(opt, a, ci)


text:
"""
    B, T = brun(H)

Brun's integer-relations alg implemented as a matrix decomposition. Takes as
input the matrix `H` and returns a reduced basis `B` and `T`, a unimodular
transformation matrix such that `B = H*T`. Brun reduction is often done with
`pinv(H)` as input to yield `B = pinv(H)*T`.

See V. Brun, "En generalisation av kjedebrøken I," Skr. Vid
ensk. Selsk. Kristiana, Mat. Nat. Klasse, 1919.  See
https://archive.org/stream/skrifterutgitavv201chri#page/300/mode/2up

Follows code from D. Wuebben, D. Seethaler, J. Jalden, and G. Matz, "Lattice
Reduction - A Survey with Applications in Wireless Communications" IEEE
Signal Processing Magazine, March 2011

# Examples
```jldoctest
julia> H=[1 2; 3 4]; B,T=brun(H); T
2×2 Array{Int64,2}:
  3  -1
 -2   1

```
"""
function brun(H::AbstractArray{Td,2}) where Td

    # get size
    K, M = size(H)

    if rank(H) < K
        error("Input basis does not have full row rank")
    end

    Ti= getIntType(Td)

    grammian = H'*H
    metric_vec = sqrt.(abs.(diag(grammian')))
    mu  = grammian[:,1] # initialize mu with an arbitrary column of the Grammian
    T = Matrix{Ti}(I, K, K)   # unimodular transformation matrix
    B = copy(H)

    doBrun = true
    while doBrun

        # Calculation of Brun's update value
        zw,s = findmax(mu)
        zw = copy(mu)
        zw[s] = -typemax(zw[1])
        zw,t = findmax(zw)
        r = round(mu[s]/mu[t])

        # basis update
        updated_vec = B[:,s] .- r'*B[:,t]

        # if metric can be improved
        updated_norm = norm(updated_vec)
        if updated_norm < metric_vec[s]

            # update mu, unimodular and reduced matrices, metric
            mu[s] -= r*mu[t]
            T[:,s] .-= r'*view(T,:,t)  #T[:,s] -= r'*T[:,t]
            B[:,s] .= updated_vec
            metric_vec[s] = updated_norm

        else
            doBrun = false
        end

    end #

    return B, T
end


text:
"""
Original publication:
Tax, David MJ, and Robert PW Duin. "Support vector data description." Machine learning 54.1 (2004): 45-66.
"""
mutable struct VanillaSVDD <: SVDDClassifier
    state::ModelState

    # model parameters
    C::Float64
    kernel_fct::Kernel

    # training data
    data::Array{Float64,2}
    K::Array{Float64,2}
    adjust_K::Bool
    K_adjusted::Array{Float64,2}

    # fitted values
    alpha_values::Vector{Float64}
    const_term::Float64
    R::Float64

    function VanillaSVDD(data)
        m = new()
        m.C = 1.0
        m.state = model_created
        m.data = data
        m.adjust_K = false
        m.const_term = -Inf
        m.R = -Inf
        return m
    end
    VanillaSVDD(data, pools) = VanillaSVDD(data)
end

get_model_params(model::VanillaSVDD) = Dict(:C => model.C)

is_valid_param_value(model::VanillaSVDD, x::Type{Val{:C}}, v) = 0 <= v <= 1

set_pools!(model::VanillaSVDD, pools::Vector{Symbol}) = nothing
set_pools!(model::VanillaSVDD, pools::Dict{Symbol, Vector{Int}}) = nothing

function set_C!(model::VanillaSVDD, C::Number)
    @assert 0 <= C <= 1
    model.C = C
    return nothing
end

function solve!(model::VanillaSVDD, solver::JuMP.OptimizerFactory)
    debug(LOGGER, "[SOLVE] Setting up QP for VanillaSVDD with $(is_K_adjusted(model) ? "adjusted" : "non-adjusted") kernel matrix.")
    QP = Model(solver)
    K = is_K_adjusted(model) ? model.K_adjusted : model.K

    @variable(QP, 0 <= α[1:size(K,1)] <= model.C)
    @objective(QP, Max, sum(α[i]*K[i,i] for i in eachindex(α)) -
                        sum(α[i]*α[j] * K[i,j] for i in eachindex(α) for j in eachindex(α)))
    @constraint(QP, sum(α) == 1)
    debug(LOGGER, "[SOLVE] Solving QP with $(typeof(solver))...")
    JuMP.optimize!(QP)
    status = JuMP.termination_status(QP)
    debug(LOGGER, "[SOLVE] Finished with status: $(status).")
    model.alpha_values = JuMP.value.(α)
    return status
end

function get_support_vectors(model::VanillaSVDD)
    findall((model.alpha_values .> OPT_PRECISION) .& (model.alpha_values .< (model.C - OPT_PRECISION)))
end

get_alpha_prime(model::VanillaSVDD) = model.alpha_values


text:
function constructdevice!(m::JuMP.AbstractModel, netinjection::BalanceNamedTuple, category::Type{B}, category_formulation::Type{D}, system_formulation::Type{S}, sys::PowerSystems.PowerSystem; args...) where {B <: PowerSystems.Branch, D <: AbstractBranchForm, S <: AbstractFlowForm}

    flowvariables(m, system_formulation, sys.branches, sys.time_periods)

    thermalflowlimits(m, system_formulation, sys.branches, sys.time_periods)

end

function constructdevice!(m::JuMP.AbstractModel, netinjection::BalanceNamedTuple, category::Type{B}, category_formulation::Type{PiLine}, system_formulation::Type{StandardPTDF}, sys::PowerSystems.PowerSystem; args...) where {B <: PowerSystems.Branch}

    PTDF = [a.second for a in args if a.first == :PTDF][1]

    ac_branches_set = [b.name for b in sys.branches if isa(b,category)]

    isempty(PTDF) ? @error("NO PTDF matrix has been provided") : (size(PTDF.axes[1])[1] != length(ac_branches_set) ? @error("PTDF size is inconsistent") : true)

    flowvariables(m, system_formulation, sys.branches, sys.time_periods)

    thermalflowlimits(m, system_formulation, sys.branches, sys.time_periods)

    dc_networkflow(m, netinjection, PTDF)

end


text:
using Dates
using ArgParse

include("../btg.jl")
include("../datasets/load_abalone.jl")

s = ArgParseSettings()
@add_arg_table! s begin
    "--fast"
        help = "use fast or not"
        action = :store_true
end
parsed_args = parse_args(ARGS, s)

ind = 1:50
posx = 1:7
posc = 1:3
x = data[ind, posx]
Fx = data[ind, posc]
y = float(target[ind])

#prediction data
pind = 1:1
Fx0 = reshape(data[pind, posc], 1, length(posc))
x0 = reshape(data[pind, posx], 1, length(posx))

#parameter setting
rangeθ = [10.0 1000]
rangeλ = [0. 3.] #we will always used 1 range scale for lambda

trainingData1 = trainingData(x, Fx, y) #training data used for testing various functions
d = getDimension(trainingData1); n = getNumPts(trainingData1); p = getCovDimension(trainingData1)
btg1 = btg(trainingData1, rangeθ, rangeλ; quadtype = ["SparseGrid", "Gaussian"])
(pdf, cdf, dpdf) = solve(btg1); #initialize training_buffer_dicts, solve once so can use fast techiques to extrapolate submatrix determinants, etc.

"""
leave-one-out-training-data
"""
function lootd(td::AbstractTrainingData, i::Int64)
    x = getPosition(td)
    Fx = getCovariates(td)
    z = getLabel(td)
    x_minus_i = x[[1:i-1;i+1:end], :]
    Fx_minus_i = Fx[[1:i-1;i+1:end], :]
    z_minus_i = z[[1:i-1;i+1:end]]
    x_i = x[i:i, :]
    Fx_i = Fx[i:i, :]
    z_i = z[i:i, :]
    return trainingData(x_minus_i, Fx_minus_i, z_minus_i), x_i, Fx_i, z_i
end
###
###
### plotting
###
###
nrow=6; ncol=8
PyPlot.close("all") #close existing windows
plt, axs = PyPlot.subplots(nrow, ncol)
PyPlot.suptitle("Cross Validation $(parsed_args["fast"])", fontsize=10)
#figure(1)
before = Dates.now()
for j = 1:nrow*ncol
    ind1 = Int64(ceil(j/ncol))
    ind2 = Int64(j - ncol*(floor((j-.1)/ncol)))
    println("iteration $j")
    if j <nrow*ncol
        if parsed_args["fast"]
            (pdf, cdf, dpdf) = solve(btg1, validate = j);
            b1 = y0 -> pdf(x0, Fx0, y0);
            c1 = y0 -> cdf(x0, Fx0, y0);
            (x, y) = plt_data(b1, .01, 1.2, 100)
            (xc, yc) = plt_data(c1, .01, 1.2, 100)
        else
            (trainingdata_minus_i, x_i, Fx_i, z_i) = lootd(trainingData1, j)
            btg2 = btg(trainingdata_minus_i, rangeθ, rangeλ; quadtype = ["SparseGrid", "Gaussian"])
            (pdf_minus_i, cdf_minus_i, dpdf_minus_i) = solve(btg2)
            b1 = y -> pdf_minus_i(x_i, Fx_i, y)
            c1 = y -> cdf_minus_i(x_i, Fx_i, y)
            (x, y) = plt_data(b1, .01, 1.2, 100)
            (xc, yc) = plt_data(c1, .01, 1.2, 100)
        end
        axs[ind1, ind2].plot(x, y, color = "red", linewidth = 1.0, linestyle = "-")
        axs[ind1, ind2].plot(xc, yc, color = "orange", linewidth = 1.0, linestyle = "-")
        axs[ind1, ind2].axvline(x =  getLabel(btg1.trainingData)[j])
    else #do annotation in last box
        after = Dates.now()
        elapsedmin = round(((after - before) / Millisecond(1000))/60, digits=5)
        qtype =  btg1.quadType
        axs[ind1, ind2].annotate(
        L"$\theta: $" * "$rangeθ" * "\n" *
        L"$\lambda: $" * "$rangeλ" * "\n" *
        "$qtype \n" *
        "ind: " * "$ind \n" *
        "x: " * "$posx " *", cov: " * "$posc\n" *
        "time: " * "$elapsedmin",
        xy=[.03;.03],
        fontsize = 5.0)
        # *
        #L"$\lambda: $" * "$rangeλ" * ".\n" *
        #"quad: " * "$btg1.quadType",
    end
end
#for ax in axs
#    ax.set(xlabel="x-label", ylabel="y-label")
#end
# Hide x labels and tick labels for top plots and y ticks for right plots.
for ax in axs
    ax.label_outer()
end

PyPlot.savefig("figure/test_validation5_$(parsed_args["fast"]).pdf")

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg