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

import LinearAlgebra: BlasInt
@static if isdefined(LinearAlgebra, :ARPACKException)
    import LinearAlgebra: ARPACKException
else
    struct ARPACKException <: Exception
        info::BlasInt
    end

    function Base.showerror(io::IO, ex::ARPACKException)
        print(io, "ARPACKException: ")
        if ex.info == -8
            print(io, "error return from calculation of a real Schur form.")
        elseif ex.info == -9
            print(io, "error return from calculation of eigenvectors.")
        elseif ex.info == -14
            print(io, string("did not find any eigenvalues to sufficient accuracy. ",
                "Try with a different starting vector or more Lanczos vectors ",
                "by increasing the value of ncv."))
        else
            print(io, "unspecified ARPACK error: $(ex.info)")
        end
    end
end

## aupd and eupd wrappers

function aupd_wrapper(T, matvecA!::Function, matvecB::Function, solveSI::Function, n::Integer,
                      sym::Bool, cmplx::Bool, bmat::String,
                      nev::Integer, ncv::Integer, which::String,
                      tol::Real, maxiter::Integer, mode::Integer, v0::Vector)
    lworkl = cmplx ? ncv * (3*ncv + 5) : (sym ? ncv * (ncv + 8) :  ncv * (3*ncv + 6) )
    TR = cmplx ? T.types[1] : T
    TOL = Ref{TR}(tol)

    v     = Matrix{T}(undef, n, ncv)
    workd = Vector{T}(undef, 3*n)
    workl = Vector{T}(undef, lworkl)
    rwork = cmplx ? Vector{TR}(undef, ncv) : Vector{TR}()

    if isempty(v0)
        resid = Vector{T}(undef, n)
        info  = Ref{BlasInt}(0)
    else
        resid = deepcopy(v0)
        info  = Ref{BlasInt}(1)
    end
    iparam = zeros(BlasInt, 11)
    ipntr  = zeros(BlasInt, (sym && !cmplx) ? 11 : 14)
    ido    = Ref{BlasInt}(0)

    iparam[1] = BlasInt(1)       # ishifts
    iparam[3] = BlasInt(maxiter) # maxiter
    iparam[7] = BlasInt(mode)    # mode

    zernm1 = 0:(n-1)

    while true
        if cmplx
            naupd(ido, bmat, n, which, nev, TOL, resid, ncv, v, n,
                  iparam, ipntr, workd, workl, lworkl, rwork, info)
        elseif sym
            saupd(ido, bmat, n, which, nev, TOL, resid, ncv, v, n,
                  iparam, ipntr, workd, workl, lworkl, info)
        else
            naupd(ido, bmat, n, which, nev, TOL, resid, ncv, v, n,
                  iparam, ipntr, workd, workl, lworkl, info)
        end
        if info[] != 0
            throw(ARPACKException(info[]))
        end

        x = view(workd, ipntr[1] .+ zernm1)
        y = view(workd, ipntr[2] .+ zernm1)
        if mode == 1  # corresponds to dsdrv1, dndrv1 or zndrv1
            if ido[] == 1
                matvecA!(y, x)
            elseif ido[] == 99
                break
            else
                throw(ARPACKException("unexpected behavior"))
            end
        elseif mode == 3 && bmat == "I" # corresponds to dsdrv2, dndrv2 or zndrv2
            if ido[] == -1 || ido[] == 1
                y[:] = solveSI(x)
            elseif ido[] == 99
                break
            else
                throw(ARPACKException("unexpected behavior"))
            end
        elseif mode == 2 # corresponds to dsdrv3, dndrv3 or zndrv3
            if ido[] == -1 || ido[] == 1
                matvecA!(y, x)
                if sym
                    x[:] = y    # overwrite as per Remark 5 in dsaupd.f
                end
                y[:] = solveSI(y)
            elseif ido[] == 2
                y[:] = matvecB(x)
            elseif ido[] == 99
                break
            else
                throw(ARPACKException("unexpected behavior"))
            end
        elseif mode == 3 && bmat == "G" # corresponds to dsdrv4, dndrv4 or zndrv4
            if ido[] == -1
                y[:] = solveSI(matvecB(x))
            elseif  ido[] == 1
                y[:] = solveSI(view(workd,ipntr[3] .+ zernm1))
            elseif ido[] == 2
                y[:] = matvecB(x)
            elseif ido[] == 99
                break
            else
                throw(ARPACKException("unexpected behavior"))
            end
        else
            throw(ArgumentError("ARPACK mode ($mode) not yet supported"))
        end
    end

    return (resid, v, n, iparam, ipntr, workd, workl, lworkl, rwork, TOL)
end

function eupd_wrapper(T, n::Integer, sym::Bool, cmplx::Bool, bmat::String,
                      nev::Integer, which::String, ritzvec::Bool,
                      TOL::Ref, resid, ncv::Integer, v, ldv, sigma, iparam, ipntr,
                      workd, workl, lworkl, rwork)
    howmny = "A"
    select = Vector{BlasInt}(undef, ncv)
    info   = Ref{BlasInt}(0)

    dmap = x -> abs.(x)
    if iparam[7] == 3 # shift-and-invert
        dmap = x -> abs.(1 ./ (x .- sigma))
    elseif which == "LR" || which == "LA" || which == "BE"
        dmap = real
    elseif which == "SR" || which == "SA"
        dmap = x -> -real(x)
    elseif which == "LI"
        dmap = imag
    elseif which == "SI"
        dmap = x -> -imag(x)
    end

    if cmplx
        d = Vector{T}(undef, nev+1)
        sigmar = Ref{T}(sigma)
        workev = Vector{T}(undef, 2ncv)
        neupd(ritzvec, howmny, select, d, v, ldv, sigmar, workev,
              bmat, n, which, nev, TOL, resid, ncv, v, ldv,
              iparam, ipntr, workd, workl, lworkl, rwork, info)
        if info[] != 0
            throw(ARPACKException(info[]))
        end

        p = sortperm(dmap(d[1:nev]), rev=true)
        return ritzvec ? (d[p], v[1:n, p],iparam[5],iparam[3],iparam[9],resid) : (d[p],iparam[5],iparam[3],iparam[9],resid)
    elseif sym
        d = Vector{T}(undef, nev)
        sigmar = Ref{T}(sigma)
        seupd(ritzvec, howmny, select, d, v, ldv, sigmar,
              bmat, n, which, nev, TOL, resid, ncv, v, ldv,
              iparam, ipntr, workd, workl, lworkl, info)
        if info[] != 0
            throw(ARPACKException(info[]))
        end

        p = sortperm(dmap(d), rev=true)
        return ritzvec ? (d[p], v[1:n, p],iparam[5],iparam[3],iparam[9],resid) : (d[p],iparam[5],iparam[3],iparam[9],resid)
    else
        dr = Vector{T}(undef, nev+1)
        di = Vector{T}(undef, nev+1)
        fill!(dr,NaN)
        fill!(di,NaN)
        sigmar = Ref{T}(real(sigma))
        sigmai = Ref{T}(imag(sigma))
        workev = Vector{T}(undef, 3*ncv)
        neupd(ritzvec, howmny, select, dr, di, v, ldv, sigmar, sigmai,
              workev, bmat, n, which, nev, TOL, resid, ncv, v, ldv,
              iparam, ipntr, workd, workl, lworkl, info)
        if info[] != 0
            throw(ARPACKException(info[]))
        end
        evec = complex.(Matrix{T}(undef, n, nev+1), Matrix{T}(undef, n, nev+1))

        j = 1
        while j <= nev
            if di[j] == 0
                evec[:,j] = v[:,j]
            else # For complex conjugate pairs
                evec[:,j]   = v[:,j] + im*v[:,j+1]
                evec[:,j+1] = v[:,j] - im*v[:,j+1]
                j += 1
            end
            j += 1
        end
        if j == nev+1 && !isnan(di[j])
            if di[j] == 0
                evec[:,j] = v[:,j]
                j += 1
            else
                throw(ARPACKException("unexpected behavior"))
            end
        end

        d = complex.(dr, di)

        if j == nev+1
            p = sortperm(dmap(d[1:nev]), rev=true)
        else
            p = sortperm(dmap(d), rev=true)
            p = p[1:nev]
        end

        return ritzvec ? (d[p], evec[1:n, p],iparam[5],iparam[3],iparam[9],resid) : (d[p],iparam[5],iparam[3],iparam[9],resid)
    end
end

for (T, saupd_name, seupd_name, naupd_name, neupd_name) in
    ((:Float64, :dsaupd_, :dseupd_, :dnaupd_, :dneupd_),
     (:Float32, :ssaupd_, :sseupd_, :snaupd_, :sneupd_))
    @eval begin
        function naupd(ido, bmat, n, evtype, nev, TOL::Ref{$T}, resid::Vector{$T}, ncv, v::Matrix{$T}, ldv,
                       iparam, ipntr, workd::Vector{$T}, workl::Vector{$T}, lworkl, info)
            ccall(($(string(naupd_name)), libarpack), Cvoid,
                  (Ref{BlasInt}, Ptr{UInt8}, Ref{BlasInt}, Ptr{UInt8}, Ref{BlasInt},
                   Ptr{$T}, Ptr{$T}, Ref{BlasInt}, Ptr{$T}, Ref{BlasInt},
                   Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ref{BlasInt}, Ref{BlasInt}, Clong, Clong),
                  ido, bmat, n, evtype, nev,
                  TOL, resid, ncv, v, ldv,
                  iparam, ipntr, workd, workl, lworkl, info, sizeof(bmat), sizeof(evtype))
        end

        function neupd(rvec, howmny, select, dr, di, z, ldz, sigmar, sigmai,
                  workev::Vector{$T}, bmat, n, evtype, nev, TOL::Ref{$T}, resid::Vector{$T}, ncv, v, ldv,
                  iparam, ipntr, workd::Vector{$T}, workl::Vector{$T}, lworkl, info)
            ccall(($(string(neupd_name)), libarpack), Cvoid,
                  (Ref{BlasInt}, Ptr{UInt8}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{$T}, Ref{BlasInt},
                   Ref{$T}, Ref{$T}, Ptr{$T}, Ptr{UInt8}, Ref{BlasInt}, Ptr{UInt8}, Ref{BlasInt},
                   Ptr{$T}, Ptr{$T}, Ref{BlasInt}, Ptr{$T}, Ref{BlasInt},
                   Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ref{BlasInt}, Ref{BlasInt}, Clong, Clong, Clong),
                  rvec, howmny, select, dr, di, z, ldz,
                  sigmar, sigmai, workev, bmat, n, evtype, nev,
                  TOL, resid, ncv, v, ldv,
                  iparam, ipntr, workd, workl, lworkl, info, sizeof(howmny), sizeof(bmat), sizeof(evtype))
        end

        function saupd(ido, bmat, n, which, nev, TOL::Ref{$T}, resid::Vector{$T}, ncv, v::Matrix{$T}, ldv,
                       iparam, ipntr, workd::Vector{$T}, workl::Vector{$T}, lworkl, info)
            ccall(($(string(saupd_name)), libarpack), Cvoid,
                  (Ref{BlasInt}, Ptr{UInt8}, Ref{BlasInt}, Ptr{UInt8}, Ref{BlasInt},
                   Ptr{$T}, Ptr{$T}, Ref{BlasInt}, Ptr{$T}, Ref{BlasInt},
                   Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ref{BlasInt}, Ref{BlasInt}, Clong, Clong),
                  ido, bmat, n, which, nev,
                  TOL, resid, ncv, v, ldv,
                  iparam, ipntr, workd, workl, lworkl, info, sizeof(bmat), sizeof(which))
        end

        function seupd(rvec, howmny, select, d, z, ldz, sigma,
                       bmat, n, evtype, nev, TOL::Ref{$T}, resid::Vector{$T}, ncv, v::Matrix{$T}, ldv,
                       iparam, ipntr, workd::Vector{$T}, workl::Vector{$T}, lworkl, info)
            ccall(($(string(seupd_name)), libarpack), Cvoid,
                  (Ref{BlasInt}, Ptr{UInt8}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ref{BlasInt},
                   Ptr{$T}, Ptr{UInt8}, Ref{BlasInt}, Ptr{UInt8}, Ref{BlasInt},
                   Ptr{$T}, Ptr{$T}, Ref{BlasInt}, Ptr{$T}, Ref{BlasInt},
                   Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ref{BlasInt}, Ref{BlasInt}, Clong, Clong, Clong),
                  rvec, howmny, select, d, z, ldz,
                  sigma, bmat, n, evtype, nev,
                  TOL, resid, ncv, v, ldv,
                  iparam, ipntr, workd, workl, lworkl, info, sizeof(howmny), sizeof(bmat), sizeof(evtype))
        end
    end
end

for (T, TR, naupd_name, neupd_name) in
    ((:ComplexF64, :Float64, :znaupd_, :zneupd_),
     (:ComplexF32, :Float32, :cnaupd_, :cneupd_))
    @eval begin
        function naupd(ido, bmat, n, evtype, nev, TOL::Ref{$TR}, resid::Vector{$T}, ncv, v::Matrix{$T}, ldv,
                       iparam, ipntr, workd::Vector{$T}, workl::Vector{$T}, lworkl,
                       rwork::Vector{$TR}, info)
            ccall(($(string(naupd_name)), libarpack), Cvoid,
                  (Ref{BlasInt}, Ptr{UInt8}, Ref{BlasInt}, Ptr{UInt8}, Ref{BlasInt},
                   Ptr{$TR}, Ptr{$T}, Ref{BlasInt}, Ptr{$T}, Ref{BlasInt},
                   Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ref{BlasInt}, Ptr{$TR}, Ref{BlasInt}),
                  ido, bmat, n, evtype, nev,
                  TOL, resid, ncv, v, ldv,
                  iparam, ipntr, workd, workl, lworkl, rwork, info)
        end

        function neupd(rvec, howmny, select, d, z, ldz, sigma, workev::Vector{$T},
                       bmat, n, evtype, nev, TOL::Ref{$TR}, resid::Vector{$T}, ncv, v::Matrix{$T}, ldv,
                       iparam, ipntr, workd::Vector{$T}, workl::Vector{$T}, lworkl,
                       rwork::Vector{$TR}, info)
            ccall(($(string(neupd_name)), libarpack), Cvoid,
                  (Ref{BlasInt}, Ptr{UInt8}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ref{BlasInt},
                   Ptr{$T}, Ptr{$T}, Ptr{UInt8}, Ref{BlasInt}, Ptr{UInt8}, Ref{BlasInt},
                   Ptr{$TR}, Ptr{$T}, Ref{BlasInt}, Ptr{$T}, Ref{BlasInt},
                   Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ref{BlasInt}, Ptr{$TR}, Ref{BlasInt}),
                  rvec, howmny, select, d, z, ldz,
                  sigma, workev, bmat, n, evtype, nev,
                  TOL, resid, ncv, v, ldv,
                  iparam, ipntr, workd, workl, lworkl, rwork, info)
        end
    end
end


text:
# Bogumił Kamiński, 2019-03-25

# Reworked example from
# https://github.com/FluxML/model-zoo/blob/master/vision/mnist/mlp.jl

using Statistics
using Flux

const TRAIN = (class=UInt8[], image=Vector{Float64}[])
const TEST = (class=UInt8[], image=Vector{Float64}[])

for (filename, data) in [("mnist_train.int", TRAIN),
                         ("mnist_test.int", TEST)]
    open(filename) do f
        while !eof(f)
            c = read(f, UInt8)
            v = reshape(read(f, 28^2), 28, 28)
            push!(data.class, c)
            push!(data.image, vec(v[4:end-3, 4:end-3]) ./ 255)
        end
    end
end

Y = Flux.onehotbatch(TRAIN.class, 0:9);
X = reduce(hcat, TRAIN.image);

m = Chain(Dense(484, 64, relu),
          Dense(64, 64, relu),
          Dense(64, 64, relu),
          Dense(64, 10),
          softmax)

loss(x, y) = Flux.crossentropy(m(x), y)
accuracy(x, y) = mean(Flux.onecold(m(x)) .== Flux.onecold(y))
dataset = Iterators.repeated((X, Y), 200)
evalcb = () -> @show(loss(X, Y))
opt = ADAM()

Flux.train!(loss, params(m), dataset, opt, cb = Flux.throttle(evalcb, 10))
accuracy(X, Y)

tY = Flux.onehotbatch(TEST.class, 0:9);
tX = reduce(hcat, TEST.image);
accuracy(tX, tY)

ext: jl

max_stars_repo_path: 6. MNIST/mnist_3_dl.jl

max_issues_repo_path: 6. MNIST/mnist_3_dl.jl

max_forks_repo_path: 6. MNIST/mnist_3_dl.jl


text:
# Julia wrapper for header: src/dai.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0
# and then hand modified by Jason Knight <jason@jasonknight.us>

import Base: length, getindex, setindex!, copy, ==, in, show, isless, searchsortedlast, deepcopy_internal

export Var, label, states
export VarSet, insert!, erase!, labels, nrStates, calcLinearState, calcState, conditionalState, conditionalStateBoth, conditionalState2
export Factor, vars, entropy, embed, normalize!, p, vars_safe
export FactorGraph, numVars, numFactors, numEdges, setBackedFactor!, clearBackups!, restoreFactors!, readFromFile
export InfAlg, init!, run!, iterations, properties, marginal, belief

typealias _VarSet Ptr{Void}
typealias _PropertySet Ptr{Void}
typealias _Factor Ptr{Void}
typealias _FactorGraph Ptr{Void}
typealias _InfAlg Ptr{Void}

############################
# Var
############################

immutable Var
  label::Csize_t
  states::Csize_t
  function Var(label, states)
    assert(label>=0)
    assert(states>0)
    new(label,states)
  end
end

typealias _Var Ptr{Var}

function label(x::Var)
  int(x.label)
end
function states(x::Var)
  int(x.states)
end
function ==(v1::Var, v2::Var)
  v1.label == v2.label
  v1.states == v2.states
end
function show(io::IO, v::Var)
  print(io, "Var(label=$(label(v)), states=$(states(v)))")
end

############################
# VarSet
############################
type VarSet
  hdl::_VarSet
  function VarSet(hdl::_VarSet, own=true)
    v = new(hdl)
    own && finalizer(v, wrapdai_varset_delete)
    v
  end
end

function VarSet()
  VarSet( ccall( (:wrapdai_varset_create, libdai), _VarSet, ()) )
end

function VarSet(vars...)
  assert(eltype(vars) == Var)
  vsp = ccall( (:wrapdai_varset_create, libdai), _VarSet, ())
  for var in vars
    ccall( (:wrapdai_varset_insert, libdai), _VarSet, (_VarSet, _Var), vsp, &var)
  end
  VarSet(vsp)
end
function +(v1::Var, v2::Var)
  VarSet(v1,v2)
end

function wrapdai_varset_delete(x::VarSet)
  ccall( (:wrapdai_varset_delete, libdai), Void, (_VarSet,), x.hdl)
end
function copy(vs::VarSet)
  VarSet(ccall( (:wrapdai_varset_clone, libdai), _VarSet, (_VarSet,), vs.hdl))
end
function deepcopy_internal(vs::VarSet, stack::ObjectIdDict)
  haskey(stack,vs) && return vs
  return stack[vs] = copy(vs)
end

function insert!(vs::VarSet, v::Var)
  ccall( (:wrapdai_varset_insert, libdai), _VarSet, (_VarSet, _Var), vs.hdl, &v)
  vs
end
function labels(vs::VarSet)
  ptrs = pointer_to_array(ccall( (:wrapdai_varset_vars, libdai), _Var, (_VarSet,), vs.hdl),
    int(length(vs)))
  Int[label(x) for x in ptrs]
end
function searchsortedlast(vs::VarSet, v::Var)
  ccall( (:wrapdai_varset_searchsortedlast, libdai), Cint, (_VarSet,_Var), vs.hdl, &v)
end
function nrStates(vs::VarSet)
  ccall( (:wrapdai_varset_nrStates, libdai), Cint, (_VarSet,), vs.hdl)
end
function length(vs::VarSet)
  ccall( (:wrapdai_varset_size, libdai), Cint, (_VarSet,), vs.hdl)
end
function erase!(vs::VarSet, v::Var)
  ccall( (:wrapdai_varset_erase, libdai), _VarSet, (_VarSet, _Var), vs.hdl, &v)
  vs
end
function sub!(vs1::VarSet, vs2::VarSet)
  ccall( (:wrapdai_varset_remove, libdai), _VarSet, (_VarSet, _VarSet), vs1.hdl, vs2.hdl)
end
function add!(vs1::VarSet, vs2::VarSet)
  ccall( (:wrapdai_varset_addm, libdai), _VarSet, (_VarSet, _VarSet), vs1.hdl, vs2.hdl)
end
function +(vs1::VarSet, vs2::VarSet)
  VarSet(ccall( (:wrapdai_varset_add, libdai), _VarSet, (_VarSet, _VarSet), vs1.hdl, vs2.hdl))
end
function +(vs::VarSet, v::Var)
  VarSet(ccall( (:wrapdai_varset_add_one, libdai), _VarSet, (_VarSet, _Var), vs.hdl, &v))
end
function -(vs1::VarSet, vs2::VarSet)
  VarSet(ccall( (:wrapdai_varset_sub, libdai), _VarSet, (_VarSet, _VarSet), vs1.hdl, vs2.hdl))
end
function -(vs::VarSet, v::Var)
  VarSet(ccall( (:wrapdai_varset_sub_one, libdai), _VarSet, (_VarSet, _Var), vs.hdl, &v))
end

function calcLinearState(vs::VarSet, statevals)
  #for (i,v) in enumerate(vars(vs))
    #0 < statevals[i] <= states(v) || throw(BoundsError())
  #end
  #assert(length(vs) == length(statevals))
  tot = 0
  for i=1:length(statevals)
    tot += (statevals[i]-1)*2^(i-1)
  end
  tot+1
end

function pris_calcLinearState(vs::VarSet, statevals)
  for (i,v) in enumerate(vars(vs))
    0 < statevals[i] <= states(v) || throw(BoundsError())
  end
  assert(length(vs) == length(statevals))
  for i=1:length(statevals)
    statevals[i] -= 1
  end
  ret = 1+ccall( (:wrapdai_varset_calcLinearState, libdai), Csize_t,
    (_VarSet, Ptr{Csize_t}, _VarSet), vs.hdl, uint(statevals), C_NULL)
  for i=1:length(statevals)
    statevals[i] += 1
  end
  ret
end

function calcLinearState(vs::VarSet, statevals, orig::VarSet)
  for (i,v) in enumerate(vars(orig))
    if statevals[i] > states(v)
      @show vs, statevals, orig
      @show i,v
      @show statevals[i], states(v)
      throw(BoundsError())
    end
  end
  assert(length(orig) == length(statevals))
  for i=1:length(statevals)
    statevals[i] -= 1
  end
  ret = 1+ccall( (:wrapdai_varset_calcLinearState, libdai), Csize_t,
    (_VarSet, Ptr{Csize_t}, _VarSet), vs.hdl, uint(statevals), orig.hdl)
  for i=1:length(statevals)
    statevals[i] += 1
  end
  ret
end

function calcState(vs::VarSet, state)
  0 < state <= nrStates(vs) || throw(BoundsError())
  states = Array(Csize_t, length(vs))
  ccall( (:wrapdai_varset_calcState, libdai), Void, (_VarSet, Csize_t, Ptr{Csize_t}),
    vs.hdl, state-1, states)
  broadcast!(+,states,states,1)
  return states
end

function conditionalState(v::Var, pars::VarSet, vstate, parstate)
  #0 < vstate <= states(v) || throw(BoundsError())
  #0 < parstate <= nrStates(pars) || throw(BoundsError())
  tot = 0
  inserti = searchsortedlast(pars,v)
  correction = 0
  ps = parstate - 1
  for i=1:length(pars)
    if inserti < i && correction == 0
      tot += (vstate-1)*2^(i-1)
      correction = 1
    end
    stateval = 2^(i+correction-1) & (ps<<correction)
    tot += stateval#*2^(i+correction-1)
  end
  if correction == 0 # append variable
    tot += (vstate-1)*2^(length(pars))
  end
  tot+1
end

function conditionalStateBoth(v::Var, pars::VarSet, parstate)
  tot0 = 0
  tot1 = 0
  inserti = searchsortedlast(pars,v)
  correction = 0
  ps = parstate - 1
  for i=1:length(pars)
    if inserti < i && correction == 0
      tot1 += 2^(i-1)
      correction = 1
    end
    stateval = 2^(i+correction-1) & (ps<<correction)
    tot0 += stateval
    tot1 += stateval
  end
  if correction == 0 # append variable
    tot1 += 2^(length(pars))
  end
  tot0+1,tot1+1
end

function pris_conditionalState(v::Var, pars::VarSet, vstate, parstate)
  0 < vstate <= states(v) || throw(BoundsError())
  0 < parstate <= nrStates(pars) || throw(BoundsError())
  1 + ccall( (:wrapdai_varset_conditionalState, libdai), Csize_t,
    (_Var, _VarSet, Csize_t, Csize_t), &v, pars.hdl, vstate-1, parstate-1)
end

function conditionalState2(all::VarSet, v1::Var, v2::Var, pars::VarSet, vstate, parstate)
  #0 < vstate <= states(v1) || throw(BoundsError())
  #0 < parstate <= nrStates(pars) || throw(BoundsError())
  tot1 = 0
  tot2 = 0
  inserti1 = searchsortedlast(all,v1) #this is 'off by one' and accounted for below
  inserti2 = searchsortedlast(all,v2)
  ps = parstate - 1
  correction = 0
  for i=1:length(all)
    if inserti1 == i
      temp = (vstate-1)*2^(i-1)
      tot1 += temp
      tot2 += temp
      correction += 1
    elseif inserti2 == i
      # tot1 missing on purpose here
      tot2 += 2^(i-1)
      correction += 1
    else
      stateval = 2^(i-1) & (ps<<correction)
      tot1 += stateval
      tot2 += stateval
    end
  end
  tot1 += 1
  tot2 += 1
  #pris = pris_conditionalState2(all,v1,v2,pars,vstate,parstate)
  #if pris != (tot1,tot2)
    #raise(Exception())
    #@show v1
    #@show v2
    #@show pars
    #@show vstate
    #@show parstate
    #@show tot1+1,tot2+1
    #@show pris
  #end
  tot1,tot2
end

function pris_conditionalState2(all::VarSet, v1::Var, v2::Var, pars::VarSet, vstate, parstate)
  0 < vstate <= states(v1) || throw(BoundsError())
  0 < parstate <= nrStates(pars) || throw(BoundsError())
  i1 = 1 + ccall( (:wrapdai_varset_conditionalState2, libdai), Csize_t,
    (_Var, _Var, _VarSet, Csize_t, Csize_t, Csize_t),
    &v1, &v2, pars.hdl, vstate-1, 0, parstate-1)
  i2 = 1 + ccall( (:wrapdai_varset_conditionalState2, libdai), Csize_t,
    (_Var, _Var, _VarSet, Csize_t, Csize_t, Csize_t),
    &v1, &v2, pars.hdl, vstate-1, 1, parstate-1)
  return i1,i2
end

function in(v::Var, vs::VarSet)
  ccall( (:wrapdai_varset_contains, libdai), Bool, (_VarSet, _Var), vs.hdl, &v)
end
function ==(vs1::VarSet, vs2::VarSet)
  ccall( (:wrapdai_varset_isequal, libdai), Bool, (_VarSet, _VarSet), vs1.hdl, vs2.hdl)
end
function isless(vs1::VarSet, vs2::VarSet)
  ccall( (:wrapdai_varset_isless, libdai), Bool, (_VarSet, _VarSet), vs1.hdl, vs2.hdl)
end
function vars(vs::VarSet)
  pointer_to_array(ccall( (:wrapdai_varset_vars, libdai), _Var, (_VarSet,), vs.hdl),
    int(length(vs)),false)
end
function show(io::IO, vs::VarSet)
  if length(vs) == 0
    print(io, "VarSet()")
  else
    println(io, "VarSet(")
    for v=vars(vs)
      show(io, v)
      println(io, ",")
    end
    print(io,")")
  end
end
############################
# Factor
############################

type Factor
  hdl::_Factor
  function Factor(hdl::_Factor)
    v = new(hdl)
    finalizer(v, wrapdai_factor_delete)
    v
  end
end
function Factor(v::Var...)
  if length(v) == 1
    return Factor(ccall( (:wrapdai_factor_create_var, libdai), _Factor, (_Var,), &v[1]))
  else
    return Factor(Varset(v))
  end
end
function Factor(vs::VarSet)
  Factor(ccall( (:wrapdai_factor_create_varset, libdai), _Factor, (_VarSet,), vs.hdl))
end
function Factor(vs::VarSet, vals::Vector{Float64})
  assert(length(vals) == nrStates(vs))
  Factor(ccall( (:wrapdai_factor_create_varset_vals, libdai), _Factor, (_VarSet, Ptr{Cdouble}), vs.hdl, vals))
end

function wrapdai_factor_delete(fac::Factor)
  ccall( (:wrapdai_factor_delete, libdai), Void, (_Factor,), fac.hdl)
end
function length(fac::Factor)
  ccall( (:wrapdai_factor_numvars, libdai), Cint, (_Factor,), fac.hdl)
end
function labels(fac::Factor)
  vsp = ccall( (:wrapdai_factor_vars_unsafe, libdai), _VarSet, (_Factor,), fac.hdl)
  lvsp = ccall( (:wrapdai_varset_size, libdai), Cint, (_VarSet,), vsp)
  ptrs = pointer_to_array(ccall( (:wrapdai_varset_vars, libdai), _Var, (_VarSet,), vsp),
    int(lvsp))
  Int[label(x) for x in ptrs]
end
function setindex!(fac::Factor, val, index)
  0 < index <= nrStates(fac) || throw(BoundsError())
  ccall( (:wrapdai_factor_set, libdai), Void, (_Factor, Csize_t, Cdouble), fac.hdl, index-1, val)
end
function getindex(fac::Factor, index)
  0 < index <= nrStates(fac) || throw(BoundsError())
  ccall( (:wrapdai_factor_get, libdai), Cdouble, (_Factor, Csize_t), fac.hdl, index-1)
end
function in(var::Var,fac::Factor)
  var in vars(fac)
end
function copy(fac::Factor)
  Factor(ccall( (:wrapdai_factor_clone, libdai), _Factor, (_Factor,), fac.hdl))
end
function deepcopy_internal(fac::Factor, stack::ObjectIdDict)
  haskey(stack,fac) && return fac
  return stack[fac] = copy(fac)
end
function vars(fac::Factor)
  VarSet(ccall( (:wrapdai_factor_vars_unsafe, libdai), _VarSet, (_Factor,), fac.hdl), false)
end
function vars_safe(fac::Factor)
  VarSet(ccall( (:wrapdai_factor_vars, libdai), _VarSet, (_Factor,), fac.hdl), true)
end
function nrStates(fac::Factor)
  int(ccall( (:wrapdai_factor_nrStates, libdai), Csize_t, (_Factor,), fac.hdl))
end
function entropy(fac::Factor)
  ccall( (:wrapdai_factor_entropy, libdai), Cdouble, (_Factor,), fac.hdl)
end
function marginal(fac::Factor, vs::VarSet)
  vsp = ccall( (:wrapdai_factor_vars_unsafe, libdai), _VarSet, (_Factor,), fac.hdl)
  assert(ccall( (:wrapdai_varset_isless, libdai), Bool, (_VarSet, _VarSet), vs.hdl, vsp))
  Factor(ccall( (:wrapdai_factor_marginal, libdai), _Factor, (_Factor, _VarSet), fac.hdl, vs.hdl))
end
function embed(fac::Factor, vs::VarSet)
  #vsp = ccall( (:wrapdai_factor_vars_unsafe, libdai), _VarSet, (_Factor,), fac.hdl)
  #assert(ccall( (:wrapdai_varset_isless, libdai), Bool, (_VarSet, _VarSet), vsp, vs.hdl))
  Factor(ccall( (:wrapdai_factor_embed, libdai), _Factor, (_Factor, _VarSet), fac.hdl, vs.hdl))
end
#function embed(fac::Factor, v::Var)
  ##vsp = ccall( (:wrapdai_factor_vars_unsafe, libdai), _VarSet, (_Factor,), fac.hdl)
  ##assert(ccall( (:wrapdai_varset_isless, libdai), Bool, (_VarSet, _VarSet), vsp, vs.hdl))
  #Factor(ccall( (:wrapdai_factor_embed_one, libdai), _Factor, (_Factor, _Var), fac.hdl, &v))
#end

function embed(fac::Factor, v::Var)
  #create blank factor with right nodes
  oldvars = vars(fac)
  inserti = searchsortedlast(oldvars,v)
  newfac = Factor(oldvars+v)
  numoldvars = length(oldvars)
  for origstate=1:nrStates(fac)
    temp = fac[origstate]
    tot0 = 0
    tot1 = 0
    correction = 0
    ps = origstate - 1
    for i=1:numoldvars
      if inserti < i && correction == 0
        #tot0 += 0*2^(i-1)
        tot1 += 2^(i-1)
        correction = 1
      end
      stateval = 2^(i+correction-1) & (ps<<correction)
      tot0 += stateval
      tot1 += stateval
    end
    if correction == 0 # append variable
      #tot0 += (vstate-1)*2^(length(oldvars))
      tot1 += 2^numoldvars
    end
    newfac[tot0+1] = temp
    newfac[tot1+1] = temp
  end
  newfac
end

function normalize!(fac::Factor)
  ccall( (:wrapdai_factor_normalize, libdai), Cdouble, (_Factor,), fac.hdl)
end
function ==(fac1::Factor, fac2::Factor)
  return ccall( (:wrapdai_factor_isequal, libdai), Bool, (_Factor, _Factor), fac1.hdl, fac2.hdl)
end
# the below is just a test, I should really copy the memory
function p(fac::Factor)
  return pointer_to_array(ccall((:wrapdai_factor_p, libdai), Ptr{Cdouble}, (_Factor,), fac.hdl), (nrStates(fac),), false)
end
function show(io::IO, fac::Factor)
  if length(fac) == 0
    print(io, "Factor()")
  else
    println(io, "Factor(")
    show(io,vars(fac))
    println(io, ",")
    for i=1:nrStates(fac)
      println(io, "$i => $(fac[i])")
    end
    print(io,")")
  end
end
############################
# FactorGraph
############################

type FactorGraph
  hdl::_FactorGraph
  function FactorGraph(hdl::_FactorGraph)
    v = new(hdl)
    finalizer(v, wrapdai_fg_delete)
    v
  end
end

function FactorGraph(factors...)
  assert(eltype(factors)==Factor)
  numfacs = length(factors)
  facs = [x.hdl for x in factors]
  FactorGraph(ccall( (:wrapdai_fg_create_facs, libdai), _FactorGraph, (_Factor, Cint), facs, numfacs))
end

function wrapdai_fg_delete(fg::FactorGraph)
  ccall( (:wrapdai_fg_delete, libdai), Void, (_FactorGraph,), fg.hdl)
end

function vars(fg::FactorGraph)
  pointer_to_array((ccall( (:wrapdai_fg_vars, libdai), _Var, (_FactorGraph,), fg.hdl)),(numVars(fg),),false)
end
function copy(fg::FactorGraph)
  FactorGraph(ccall( (:wrapdai_fg_clone, libdai), _FactorGraph, (_FactorGraph,), fg.hdl))
end
function deepcopy_internal(fg::FactorGraph, stack::ObjectIdDict)
  haskey(stack,fg) && return fg
  return stack[fg] = copy(fg)
end
function numVars(fg::FactorGraph)
  int(ccall( (:wrapdai_fg_nrVars, libdai), Csize_t, (_FactorGraph,), fg.hdl))
end
function numFactors(fg::FactorGraph)
  int(ccall( (:wrapdai_fg_nrFactors, libdai), Csize_t, (_FactorGraph,), fg.hdl))
end
function numEdges(fg::FactorGraph)
  int(ccall( (:wrapdai_fg_nrEdges, libdai), Csize_t, (_FactorGraph,), fg.hdl))
end
function printBPG(fg::FactorGraph)
  ccall( (:wrapdai_fg_printBPG, libdai), Void, (_FactorGraph,), fg.hdl)
end
function getindex(fg::FactorGraph, ind)
  0 < ind <= numFactors(fg) || throw(BoundsError())
  Factor(ccall( (:wrapdai_fg_factor, libdai), _Factor, (_FactorGraph, Cint), fg.hdl, ind-1))
end
function setindex!(fg::FactorGraph, fac::Factor, ind)
  0 < ind <= numFactors(fg) || throw(BoundsError())
  ccall( (:wrapdai_fg_setFactor, libdai), Void,
    (_FactorGraph, Cint, _Factor), fg.hdl, ind-1, fac.hdl)
end
function ==(fg1::FactorGraph, fg2::FactorGraph)
  nf1 = numFactors(fg1)
  nf2 = numFactors(fg2)
  nf1 == nf2 || return false
  for i = 1:nf1
    fg1[i] == fg2[i] || return false
  end
  return true
end
function setBackedFactor!(fg::FactorGraph, ind, fac::Factor)
  0 < ind <= numFactors(fg) || throw(BoundsError())
  ccall( (:wrapdai_fg_setFactor_backup, libdai), Void,
    (_FactorGraph, Cint, _Factor), fg.hdl, ind-1, fac.hdl)
end
#function getVar(fg::FactorGraph, ind)
  #0 < ind <= numVars(fg) || throw(BoundsError())
  #Var(ccall( (:wrapdai_fg_var, libdai), _Var, (_FactorGraph, Csize_t), fg.hdl, ind-1))
#end
function clearBackups!(fg::FactorGraph)
  ccall( (:wrapdai_fg_clearBackups, libdai), Void, (_FactorGraph,), fg.hdl)
end
function restoreFactors!(fg::FactorGraph)
  ccall( (:wrapdai_fg_restoreFactors, libdai), Void, (_FactorGraph,), fg.hdl)
end
function readFromFile(fg::FactorGraph, text)
  FactorGraph(ccall( (:wrapdai_fg_readFromFile, libdai),
    Void, (_FactorGraph, Ptr{Uint8}), fg.hdl, text))
end
function show(io::IO, fg::FactorGraph)
  if numVars(fg) == 0
    print(io, "FactorGraph()")
  else
    println(io, "FactorGraph(")
    for v=vars(fg)
      show(io, v)
      println(io, ",")
    end
    println(io, "Factors: ")
    for i=1:numFactors(fg)
      print(io, "$i: ")
      show(io, fg[i])
      println(io,",")
    end
    print(io,")")
  end
end
############################
# JunctionTree
############################

type InfAlg
  hdl::_InfAlg
  function InfAlg(hdl::_InfAlg)
    v = new(hdl)
    finalizer(v, wrapdai_infalg_delete)
    v
  end
end
function InfAlg(fg::FactorGraph, desc="")
  #   OLD:
  #   full = "JTREE[updates=HUGIN,verbose=0]"
  #   full = "BP[inference=SUMPROD,verbose=0,updates=SEQRND,logdomain=0,tol=1e-9,maxiter=10000,damping=0.0]"

  # From libDAI tests/aliases
  # full = "BP[inference=SUMPROD,updates=SEQMAX,logdomain=0,tol=1e-9,maxiter=10000,damping=0.0]"
  # full = "FBP[inference=SUMPROD,updates=SEQMAX,logdomain=0,tol=1e-9,maxiter=10000,damping=0.0]"
  # full = "TRWBP[updates=SEQFIX,tol=1e-9,maxiter=10000,logdomain=0,nrtrees=0]"
  full = "JTREE[inference=SUMPROD,updates=HUGIN]"
  # full = "JTREE[inference=SUMPROD,updates=SHSH]"
  # full = "MF[tol=1e-9,maxiter=10000,damping=0.0,init=UNIFORM,updates=NAIVE]"
  # full = "TREEEP[type=ORG,tol=1e-9,maxiter=10000]"
  # full = "MR[updates=FULL,inits=RESPPROP,tol=1e-9]"
  # full = "HAK[doubleloop=0,clusters=MIN,init=UNIFORM,tol=1e-9,maxiter=10000]"
  # full = "HAK[doubleloop=0,clusters=LOOP,init=UNIFORM,loopdepth=4,tol=1e-9,maxiter=10000]"
  # full = "HAK[doubleloop=1,clusters=LOOP,init=UNIFORM,loopdepth=6,tol=1e-9,maxiter=10000]"
  # full = "LC[cavity=FULL,reinit=1,updates=SEQRND,maxiter=10000,cavainame=BP,cavaiopts=[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0],tol=1e-9]"
  # full = "LC[cavity=FULL,reinit=1,updates=SEQFIX,maxiter=10000,cavainame=TREEEP,cavaiopts=[type=ORG,tol=1e-9,maxiter=10000],tol=1e-9]"
  # full = "GIBBS[maxiter=10000,burnin=100,restart=10000]"
  # full = "DECMAP[ianame=BP,iaopts=[inference=MAXPROD,updates=SEQRND,logdomain=1,tol=1e-9,maxiter=10000,damping=0.1,verbose=0],reinit=1,verbose=0]"
  # full = "GLC[rgntype=SINGLE,cavity=UNIFORM,updates=SEQRND,maxiter=10000,tol=1e-9,cavainame=BP,cavaiopts=[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0],inainame=EXACT,inaiopts=[],tol=1e-9,verbose=0]"
  InfAlg(ccall( (:wrapdai_newInfAlgFromString, libdai), _InfAlg, (Ptr{Uint8}, _FactorGraph), length(desc) == 0 ? full : desc, fg.hdl))
end

function wrapdai_infalg_delete(ia::InfAlg)
  ccall( (:wrapdai_ia_delete, libdai), Void, (_InfAlg,), ia.hdl)
end

function copy(ia::InfAlg)
  InfAlg(ccall( (:wrapdai_ia_clone, libdai), _InfAlg, (_InfAlg,), ia.hdl))
end

function deepcopy_internal(ia::InfAlg, stack::ObjectIdDict)
  haskey(stack,ia) && return ia
  return stack[ia] = copy(ia)
end

function init!(ia::InfAlg)
  ccall( (:wrapdai_ia_init, libdai), Void, (_InfAlg,), ia.hdl)
end
function run!(ia::InfAlg)
  ccall( (:wrapdai_ia_run, libdai), Void, (_InfAlg,), ia.hdl)
end
function iterations(ia::InfAlg)
  ccall( (:wrapdai_ia_iterations, libdai), Csize_t, (_InfAlg,), ia.hdl)
end
function properties(ia::InfAlg)
  ptr = ccall( (:wrapdai_ia_printProperties, libdai), Ptr{Uint8}, (_InfAlg,), ia.hdl)
  bytestring(ptr)
end
function marginal(ia::InfAlg, vs::VarSet, reinit=false)
  Factor(ccall( (:wrapdai_ia_calcMarginal, libdai), _Factor, (_InfAlg, _VarSet, Cuchar), ia.hdl, vs.hdl, reinit))
end
function belief(ia::InfAlg, vs::VarSet)
  Factor(ccall( (:wrapdai_ia_belief, libdai), _Factor, (_InfAlg, _VarSet), ia.hdl, vs.hdl))
end

############################
# PropertySet
############################
function wrapdai_ps_create(name)
  ccall( (:wrapdai_ps_create, libdai), _PropertySet, (Ptr{Uint8},), name)
end
function wrapdai_ps_delete(ps::_PropertySet)
  ccall( (:wrapdai_ps_delete, libdai), _PropertySet, (_PropertySet,), ps)
end


text:
reload("Wkhtmltox")

module Try
end

module Try
using Wkhtmltox
using Base.Test

using Base.Markdown

md"""
list :
  - **abcs** : dvfqsdf
  - **qsdfs** qsdfsdf
"""

#### conversion to pdf ####

# exfile = joinpath(dirname(@__FILE__),"../examples/index.html")
exfile = joinpath(dirname(@__FILE__),"../examples/example.html")
# outfile = tempname() * ".pdf"
outfile = "c:/temp/img.png"

img_init(1)

is = ImgSettings("in" => exfile,
                 "out" => outfile,
                 "fmt" => "png")

conv = ImgConverter(is)

run(conv)

conv = nothing
img_deinit()

#### conversion to image ####


text:
export plot_mixture_by_component, plot_mixture_sequential

using PyPlot

function plot_mixture_by_component(d::AbstractMixtureModel, x; lw=1, show_legend=false)
    @assert length(d) == 1
    K = length(probs(d))
    cmap = ColorMap("bwr")
    for (i, j) in enumerate(sort(collect(1:K), by = (x -> probs(d)[x]), rev=true))
        plot(x, probs(d)[j] * pdf(d.components[j], x), color=cmap.o(i/K*0.45),
        lw=lw,
        label=@sprintf("component %d: p=%.1f", i, probs(d)[i]))
    end
    plot(x, pdf(d, x), lw=1, ls=":", color="red", alpha=0.7, label="overall")
    if show_legend
        legend()
    end
end

function plot_mixture_sequential(d::AbstractMixtureModel, x; lw=1, rescaled=true)
    @assert length(d) == 1
    K = length(probs(d))
    cmap = ColorMap("Blues")
    subplot(2,1,1)
    for k in 1:K
        if probs(d)[k] < 1e-3
            continue
        end
        priors = probs(d)[1:k]
        w_sum = sum(priors)
        priors = priors / w_sum
        sub_d = MixtureModel(d.components[1:k], priors)
        if rescaled
            plot(x, pdf(sub_d, x), lw=lw, color=cmap.o(w_sum))
        else
            plot(x, w_sum * pdf(sub_d, x), lw=lw, color=cmap.o(w_sum))
        end
    end
    xlabel("x"); ylabel("density")
    subplot(2,1,2)
    w = zeros(K)
    for k in 1:K
        w[k] = k==1 ? probs(d)[k] : probs(d)[k] + w[k-1]
        plot(k, w[k], ".", ms=10, color=cmap.o(w[k]))
    end
    plot(collect(1:K), w, "--", color="blue", alpha=0.5)
    xlabel("component added"); ylabel("sum of p")
end

"plot fitness score vs. iterations"
function plot_loss(iter_vec, loss_vec::Vector{Float64}; color="blue")
    plot(iter_vec, loss_vec, ".-", color=color)
    xlabel("iteration"); ylabel("cross entropy"); grid()
end

"plot fitness score vs. iterations, with best possible fitness score"
function plot_loss{T<:Integer}(iter_vec::AbstractVector{T}, loss_vec::Vector{Float64}, true_d::Distribution; color="blue")
    plot_loss(iter_vec, loss_vec, color=color)
    axhline(y=self_entropy(true_d), ls="--", color="red")
end

function plot_loss(loss_vec::Vector{Float64}, true_d::Distribution; color="blue")
    plot_loss(1:length(loss_vec), loss_vec, true_d, color=color)
end

"generate a meshgrid and evaluate f(.) on each point"
function ndgrid_eval{T}(v1::AbstractVector{T}, v2::AbstractVector{T}, f::Function)
    m, n = length(v1), length(v2)
    v1 = reshape(v1, m, 1)
    v2 = reshape(v2, 1, n)
    X = repmat(v1, 1, n)
    Y = repmat(v2, m, 1)
    Z = zeros(m, n)
    for i in 1:m
        for j in 1:n
            Z[i, j] = f([X[i,j], Y[i,j]])
        end
    end
    return (X, Y, Z)
end

"plot density in 3D"
function plot_density_3D{T}(d::Distribution, x::AbstractVector{T}=linspace(-3,3,128), y::AbstractVector{T}=linspace(-3,3,128))
    X, Y, Z = ndgrid_eval(x, y, z -> pdf(d, z))
    cmap = ColorMap("bwr")
    subplot(2,2,1, projection="3d")
    plot_surface(X, Y, Z, alpha=0.3)
    xlabel("x"); ylabel("y");
    subplot(2,2,4)
    contour(X, Y, Z); grid()
    xlabel("x"); ylabel("y")
    subplot(2,2,2)
    contour(Z, Y, X); grid()
    ylabel("y")
    subplot(2,2,3)
    contour(X, Z, Y); grid()
    xlabel("x")

end

"plot and compare 2D distributions"
function plot_compare_3D{T}(true_d::Distribution, approx_d::Distribution,
                            x::AbstractVector{T}, y::AbstractVector{T}; n=1000)
    X, Y, Z1 = ndgrid_eval(x, y, z -> pdf(true_d, z))
    X, Y, Z2 = ndgrid_eval(x, y, z -> pdf(approx_d, z))

    figure()
    plot_surface(X, Y, Z1, rstride=8, cstride=8, alpha=0.3, color="blue")
    plot_surface(X, Y, Z2, rstride=8, cstride=8, alpha=0.3, color="green")
    xlabel("x"); ylabel("y"); zlabel("density")

    figure()
    samples_true = rand(true_d, n)
    samples_approx = rand(approx_d, n)
    plot(samples_true[1,:][:], samples_true[2,:][:], color="blue", alpha=0.3, ".", label="true")
    plot(samples_approx[1,:][:], samples_approx[2,:][:], color="green", alpha=0.3, ".", label="approx")
    legend()
end

"read libsvm format"
function read_libsvm_format(filename::AbstractString, n::Int, p::Int)
    y = Vector{Int}(n)
    X = zeros(n, p)
    open(filename) do fr
        iter = 0
        for line in eachline(fr)
            iter += 1
            x = zeros(p)
            for (i,s) in enumerate(split(rstrip(line)))
                if i==1
                    y[iter] = parse(Int, s)
                else
                    j, v = split(s, ":")
                    j = parse(Int, j)
                    v = parse(Float64, v)
                    X[iter, j] = v
                end
            end
            if iter==n
                break
            end
        end
    end
    @show size(X), size(y)
    return X, y
end

max_stars_repo_name: richardkwo/boosting-VB

max_issues_repo_name: richardkwo/boosting-VB

max_forks_repo_name: richardkwo/boosting-VB


text:
module Time

using Dates

using ..Mattermost

function time(req)
    return OutgoingWebhookResponse(Dates.format(Dates.now(Dates.UTC), Dates.ISODateTimeFormat))
end

end


text:
"""
   Initialize(petsclib)

Initialize the `petsclib`.

Manual: [`PetscInitializeNoArguments`](https://petsc.org/release/docs/manualpages/Sys/PetscInitializeNoArguments.html)
"""
function Initialized end
@for_libpetsc function Initialized(::$UnionPetscLib)
    r_flag = Ref{PetscBool}()
    @chk ccall(
        (:PetscInitialized, $petsc_library),
        PetscErrorCode,
        (Ptr{PetscBool},),
        r_flag,
    )
    return r_flag[] == PETSC_TRUE
end

"""
   Initialized(petsclib)

Check if `petsclib` is initialized

Manual: [`PetscInitialized`](https://petsc.org/release/docs/manualpages/Sys/PetscInitialized.html)
"""
function Initialize end
@for_libpetsc function Initialize(::$UnionPetscLib)
    if !Initialized($petsclib)
        MPI.Initialized() || MPI.Init()
        @chk ccall(
            (:PetscInitializeNoArguments, $petsc_library),
            PetscErrorCode,
            (),
        )

        # disable signal handler
        @chk ccall((:PetscPopSignalHandler, $petsc_library), PetscErrorCode, ())

        atexit(() -> Finalize($PetscLib))
    end
    return nothing
end

"""
   Finalize(petsclib)

Check if `petsclib` is initialized

Manual: [`PetscFinalize`](https://petsc.org/release/docs/manualpages/Sys/PetscFinalize.html)
"""
function Finalize end
@for_libpetsc function Finalize(::$UnionPetscLib)
    if !Finalized($petsclib)
        @chk ccall((:PetscFinalize, $petsc_library), PetscErrorCode, ())
    end
    return nothing
end

"""
   Finalized(petsclib)

Check if `petsclib` is finalized

Manual: [`PetscFinalized`](https://petsc.org/release/docs/manualpages/Sys/PetscFinalized.html)
"""
function Finalized end
@for_libpetsc function Finalized(::$UnionPetscLib)
    r_flag = Ref{PetscBool}()
    @chk ccall(
        (:PetscFinalized, $petsc_library),
        PetscErrorCode,
        (Ptr{PetscBool},),
        r_flag,
    )
    return r_flag[] == PETSC_TRUE
end


text:
Base.@kwdef struct SpringParam
    ls::Float64 # Spring length
    ms::Float64 # Mass
    ks::Float64 # Spring Strength
    ns::Int     # Mass number
end
getvalues(p::SpringParam) = p.ls,p.ms,p.ks,p.ns

Base.@kwdef struct InitParam
    λ::Float64      # Wavelength of the initial deformation
    shift::Float64  # Maximal initial node displacement
    pos::Float64    # Central position of the node displacement
end
getvalues(p::InitParam) = p.λ,p.shift,p.pos

Base.@kwdef struct AnimParam
    δt::Float64      # Time step for the dynamic simulation
    nδt::Int         # Number of timesteps
    nδtperframe::Int # Number of timesteps per animation frame
end
getvalues(p::AnimParam) = p.δt,p.nδt,p.nδtperframe


text:
module FactoredEikonalFastMarching
export EikonalParam
export EikonalTempMemory
export getEikonalParam,getEikonalTempMemory
export solveFastMarchingUpwindGrad

using jInv.Mesh

EikIdxType  = Int32
EikTtype	= Float64 # here we choose the floating type for the calculation
export EikIdxType
const imposeMonotonicity = false;

"""
	type EikonalParam

	Stores results of a Factored Eikonal forward problem

	| nabla (T0(x) x T1(x)) | = m(x)

	where T0 is the analytical solution for m=1 and T1 is the factor
	to be computed.

	Fields:

	Mesh
	kappaSquared::Array{Float64}    - squared slowness model
	src         ::Array{Int64,1}    - index of source location
	HO          ::Bool              - flag for higher-order discretization
	T1          ::Array             - factor of Eikonal
	ordering    ::Array{EikIdxType} - order of points as were determined by
	                                  Fast Marching (used in sensitivity computation)
	OP          ::Array{Int8}       - type of forward/backward operator in each
	                                  voxel (used in sensitivity computation)

    Constructor: getEikonalParam
"""
mutable struct EikonalParam
    Mesh        :: RegularMesh # Mesh
	kappaSquared:: Array{Float64} # Squared slowness model
	src         :: Array{Int64,1}
	HO          :: Bool
	T1 			:: Array
	ordering	:: Array{EikIdxType,1}
	OP		    :: Array{Int8};
end

"""
mutable struct EikonalTempMemory
This struct is used for having a reusable memory allocation for consecutive runs of eikonal solves.
Done:: An array of booleans for knowing which variables are known.
J	:: An array of indices for the heap.
V   :: An array of values for the heap.
"""
mutable struct EikonalTempMemory
	Done		:: Array{Bool};
	J  			:: Array{Int64,1};
	V			:: Array{Float64,1};
end

function getEikonalParam(Mesh:: RegularMesh,kappaSquared::Array{Float64},src:: Array{Int64,1},HO:: Bool)
	if prod(Mesh.n.+1) > (2^31-1)
		error("Eikonal is using Int32 for storing the sensitivities and N is too large. Change to EikIdxType in Eikoonal.jl to Int64");
	end
	pEik = EikonalParam(Mesh,kappaSquared,src,HO,[],[],[]);
	return pEik;
end

function getEikonalTempMemory(n::Array{Int64,1})
ntup = tuple(n...);
dim = length(n);
N = prod(n);
mem = EikonalTempMemory(zeros(Bool,tuple(n.+4...)),zeros(Int64,N),zeros(Float64,N));
return mem;
end

import jInv.Utils.clear!
function clear!(param::EikonalParam)
param.kappaSquared = [];
param.T1 = [];
param.ordering = [];
param.OP = [];
return param;
end

function solveFastMarchingUpwindGrad(pEik::EikonalParam, mem::EikonalTempMemory)
if minimum(pEik.kappaSquared) < 1e-16
	error("Fast Marching::kappaSquared is negative somewhere");
end

dim = pEik.Mesh.dim;
if dim==3
	factoredFastMarching3D(pEik,mem);
elseif dim==2
	factoredFastMarching(pEik,mem);
end
return;
end

include("heap.jl");
include("helpFuncs2D.jl");
include("helpFuncs3D.jl");
include("factoredFastMarching2D.jl");
include("factoredFastMarching3D.jl");
include("getUpdatedFactoredVal2D.jl");
include("getUpdatedFactoredVal3D.jl");
include("getAnalyticEikonalSolution2D.jl");
include("getAnalyticEikonalSolution3D.jl");
include("getAnalyticalMediums.jl")
include("EikSensitivities.jl");

end


text:
# provides compatibility with Distributions.jl
# maintained by @zsunberg

import Distributions
import Distributions: Distribution, UnivariateDistribution, MultivariateDistribution, Matrixvariate
import Distributions: Categorical, MvNormal

@generated function rand(rng::AbstractRNG, d::Distribution)
    Core.println("""
         WARNING: You are using a $d distribution from Distributions.jl. This will work, but simulations will not be repeatable because Distributions.jl does not support non-global random number generators (i.e. there is no `rand(rng::AbstractRNG, d::$d)`)

         For more information, or help Distributions.jl to support this, see https://github.com/JuliaStats/Distributions.jl/issues/436

         To disable this warning, use `rand(rng::AbstractRNG, d::$d) = rand(d)` (you will still not be able to repeat simulations).
         """)
    return quote
        rand(d)
    end
end

rand(rng::AbstractRNG, d::UnivariateDistribution) = Distributions.quantile(d, rand(rng)) # fallback

function Base.iterate(d::Categorical)
    iterate(1:Distributions.ncategories(d))
end
function Base.iterate(d::Categorical, state)
    iterate(1:Distributions.ncategories(d), state)
end

support(d::Categorical) = 1:Distributions.ncategories(d)

# sampletype(d::UnivariateDistribution) = eltype(d)
# sampletype(d::MultivariateDistribution) = Vector{eltype(d)}
# sampletype(d::Distribution{Matrixvariate}) = Matrix{eltype(d)}

# for MvNormal - this should be removed once Distributions.jl PR #597 is in a tagged Distributions.jl release
rand(rng::AbstractRNG, d::MvNormal) = _rand!(rng, d, Vector{eltype(d)}(undef, length(d)))

function _rand!(rng::AbstractRNG, d::MvNormal, x::VecOrMat)
    Distributions.add!(Distributions.unwhiten!(d.Σ, randn!(rng, x)), d.μ)
end


text:
using QuickPOMDPs
using Test

using POMDPs
using POMDPPolicies
using POMDPSimulators
using BeliefUpdaters
using POMDPModelTools
using Random
using Compose

include("discrete_explicit.jl")
include("quick.jl")


text:
# 'cstrt_bieri_ngon' provides the vertices coordinates (a,b) of a small
# n-gon for n = 2m and m >= 3
function cstrt_bieri_ngon(n)
if mod(n,2) == 0 && n>=6
    # initialization
    s = (2*sqrt(114)-7)/22;
    t = (84*s^2-272*s+175)/(4*(22*s+7));
    c = (7792*s^4+16096*s^3+2568*s^2-6248*s+223)*pi^2/(768*(22*s+7))-(226*s^2+84*s*t-22*t^2-542*s-136*t+303)/(2*(22*s+7));
    u = s*pi/n + t*pi/n^2 - c*pi/n^3;
    v = (pi/2-u)/(n/2-1);
    d = asin(sin(u)+sin(u+3*v/2)/(2*cos(v/2)))-u-v;
    x = zeros(Int(n/2));
    x[1] = u;
    x[2] = v+d;
    x[3] = v-d;
    if n >= 8
        x[4:end] .= v;
    end
    # construction
    a = zeros(n-1); b = zeros(n-1);
    a[Int(n/2-1)] = sin(x[1]); a[Int(n/2+1)] = -a[Int(n/2-1)];
    b[Int(n/2-1)] = cos(x[1]); b[Int(n/2+1)] = b[Int(n/2-1)];
    for i = 2:Int(n/2-1)
        a[Int(mod(i*(n-1)/2-mod(i,2)/2,n))] = a[Int(mod((i-1)*(n-1)/2-mod(i-1,2)/2,n))]-(-1)^i*sin(sum(x[1:i]));
        a[Int(n-mod(i*(n-1)/2-mod(i,2)/2,n))] = -a[Int(mod(i*(n-1)/2-mod(i,2)/2,n))];
        b[Int(mod(i*(n-1)/2-mod(i,2)/2,n))] = b[Int(mod((i-1)*(n-1)/2-mod(i-1,2)/2,n))]-(-1)^i*cos(sum(x[1:i]));
        b[Int(n-mod(i*(n-1)/2-mod(i,2)/2,n))] = b[Int(mod(i*(n-1)/2-mod(i,2)/2,n))];
    end
    a[Int(n/2)] = 0; b[Int(n/2)] = 1;
    return(a,b)
end
end

max_stars_repo_name: cbingane/optigon

max_issues_repo_name: cbingane/optigon

max_forks_repo_name: cbingane/optigon


text:
using Test
using NNlib
using Zygote
using NNlibCUDA
using ForwardDiff: Dual
using CUDA
CUDA.allowscalar(false)

include("test_utils.jl")

if CUDA.has_cuda()
    include("activations.jl")
    include("batchedmul.jl")
    include("upsample.jl")
    include("conv.jl")
    include("pooling.jl")
    include("softmax.jl")
    include("batchnorm.jl")
end


text:
module QCGANS
using Compat

using Yao
using Yao.Zoo
using Yao.Blocks

import Yao.Registers: tracedist
import Base: gradient

export QCGAN, qcgan, p0g, p0t, loss
"""
Quantum Circuit GAN.

Reference:
    Benedetti, M., Grant, E., Wossnig, L., & Severini, S. (2018). Adversarial quantum circuit learning for pure state approximation, 1–14.
"""
struct QCGAN{N, GT<:MatrixBlock{N}, DT<:MatrixBlock, OT<:MatrixBlock, CT<:AbstractBlock}
    target::DefaultRegister
    generator::GT
    discriminator::DT
    reg0::DefaultRegister
    witness_op::OT

    circuit::CT
    grots::Vector{RotationGate}
    drots::Vector{RotationGate}
end

function QCGAN(target::DefaultRegister, gen::MatrixBlock, dis::MatrixBlock)
    N = nqubits(target)
    c = sequence(gen, addbit(1), dis)
    witness_op = put(N+1, (N+1)=>P0)
    QCGAN{N, typeof(gen), typeof(dis), typeof(witness_op), typeof(c)}(target, gen, dis, zero_state(N), witness_op, c, collect_rotblocks(gen), collect_rotblocks(dis))
end

function qcgan(target::DefaultRegister, depth_gen::Int, depth_disc::Int, pairs1::Vector{<:Pair}, pairs2::Vector{<:Pair})
    n = nqubits(target)
    generator = diff_circuit(n, depth_gen, pairs1)
    discriminator = diff_circuit(n+1, depth_disc, pairs2)
    return QCGAN(target, generator, discriminator)
end

p0g(qcg::QCGAN) = ()->expect(qcg.witness_op, apply!(copy(qcg.reg0), qcg.circuit)) |> real
p0t(qcg::QCGAN) = ()->expect(qcg.witness_op, apply!(copy(qcg.target), qcg.circuit[2:end])) |> real
loss(qcg::QCGAN) = ()->(qcg |> p0t)() - (qcg |> p0g)()

tracedist(qcg::QCGAN) = tracedist(qcg.target, apply!(copy(qcg.reg0), qcg.generator))[]

function gradient(qcg::QCGAN)
    ggrad_g = opgrad(qcg |> p0g, qcg.grots)
    dgrad_g = opgrad(qcg |> p0g, qcg.drots)
    dgrad_t = opgrad(qcg |> p0t, qcg.drots)
    -ggrad_g, dgrad_t - dgrad_g
end

end

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
module CalcNumberLoadBearing

include("./../../constraints/geometric/constants.jl")
using JuMP
using .Constants

export calc_number_lb

function calc_number_lb(m)

    @expression(m, num_of_boxes_on_top_i[i=1:length(boxes)],
        sum(m[:if_simple_lb_var][i, j] for j in 1:length(boxes) if i!=j)
    )

    @constraint(m, c_num_of_boxes_on_top_i[i=1:length(boxes)-1],
        m[:num_of_boxes_on_top_i][i] <= lb_number_of_boxes
    )

    return m
end

end


text:
function RLCore.Experiment(
    ::Val{:JuliaRL},
    ::Val{:SAC},
    ::Val{:Pendulum},
    ::Nothing;
    save_dir = nothing,
    seed = 123,
)
    if isnothing(save_dir)
        t = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
        save_dir = joinpath(pwd(), "checkpoints", "JuliaRL_SAC_Pendulum_$(t)")
    end

    lg = TBLogger(joinpath(save_dir, "tb_log"), min_level = Logging.Info)
    rng = StableRNG(seed)
    inner_env = PendulumEnv(T = Float32, rng = rng)
    A = action_space(inner_env)
    low = A.left
    high = A.right
    ns = length(state(inner_env))
    na = 1

    env = ActionTransformedEnv(
        inner_env;
        action_mapping = x -> low + (x[1] + 1) * 0.5 * (high - low),
    )
    init = glorot_uniform(rng)

    create_policy_net() = NeuralNetworkApproximator(
        model = GaussianNetwork(
            pre = Chain(
                Dense(ns, 30, relu),
                Dense(30, 30, relu),
            ),
            μ = Chain(Dense(30, na, initW = init)),
            logσ = Chain(Dense(30, na, x -> clamp.(x, typeof(x)(-10), typeof(x)(2)), initW = init)),
        ),
        optimizer = ADAM(0.003),
    )

    create_q_net() = NeuralNetworkApproximator(
        model = Chain(
            Dense(ns + na, 30, relu; initW = init),
            Dense(30, 30, relu; initW = init),
            Dense(30, 1; initW = init),
        ),
        optimizer = ADAM(0.003),
    )

    agent = Agent(
        policy = SACPolicy(
            policy = create_policy_net(),
            qnetwork1 = create_q_net(),
            qnetwork2 = create_q_net(),
            target_qnetwork1 = create_q_net(),
            target_qnetwork2 = create_q_net(),
            γ = 0.99f0,
            ρ = 0.995f0,
            α = 0.2f0,
            batch_size = 64,
            start_steps = 1000,
            start_policy = RandomPolicy(Space([-1.0..1.0 for _ in 1:na]); rng = rng),
            update_after = 1000,
            update_every = 1,
            rng = rng,
        ),
        trajectory = CircularArraySARTTrajectory(
            capacity = 10000,
            state = Vector{Float32} => (ns,),
            action = Vector{Float32} => (na,),
        ),
    )

    stop_condition = StopAfterStep(10_000)
    total_reward_per_episode = TotalRewardPerEpisode()
    time_per_step = TimePerStep()
    hook = ComposedHook(
        total_reward_per_episode,
        time_per_step,
        DoEveryNStep() do t, agent, env
            with_logger(lg) do
                @info(
                    "training",
                    reward_term = agent.policy.reward_term,
                    entropy_term = agent.policy.entropy_term,
                )
                if is_terminated(env)
                    @info "training" reward = total_reward_per_episode.reward log_step_increment = 0
                end
            end
        end,
    )

    Experiment(agent, env, stop_condition, hook, "# Play Pendulum with SAC")
end


text:
struct BDM3DRefSpace{T} <: RefSpace{T,12} end

function (f::BDM3DRefSpace)(p)

    u,v,w = parametric(p)

    tu = tangents(p,1)
    tv = tangents(p,2)
    tw = tangents(p,3)

    j = jacobian(p)
    d=1/j

    return SVector((

        (value= 2*(-v*tu+v*tv)/j , divergence= 2*d),
        (value= 2*(-w*tu+w*tw)/j , divergence= 2*d),
        (value= 2*((u+v+w-1)*tu)/j , divergence= 2*d),

        (value= 2*(-w*tv+w*tw)/j , divergence= 2*d),
        (value= 2*((u+v+w-1)*tv)/j , divergence= 2*d),
        (value= 2*(u*tu-u*tv)/j , divergence= 2*d),

        (value= 2*((u+v+w-1)*tw)/j , divergence= 2*d),
        (value= 2*(u*tu-u*tw)/j , divergence= 2*d),
        (value= 2*(v*tv-v*tw)/j  , divergence= 2*d),

        (value= 2*(u*tu)/j      , divergence= 2*d),
        (value= 2*(v*tv)/j       , divergence= 2*d),
        (value= 2*(w*tw)/j       , divergence= 2*d)

    ))
end


text:
#using Revise
using CmdStan
using Plots
using StatsPlots

#%% Input Data #################################################################
y = [607, 583, 521, 494, 369, 782, 570, 678, 467, 620, 425, 395, 346, 361, 310,
300, 382, 294, 315, 323, 421, 339, 398, 328, 335, 291, 329, 310, 294, 321, 286,
349, 279, 268, 293, 310, 259, 241, 243, 272, 247, 275, 220, 245, 268, 357, 273,
301, 322, 276, 401, 368, 149, 507, 411, 362, 358, 355, 362, 324, 332, 268, 259,
274, 248, 254, 242, 286, 276, 237, 259, 251, 239, 247, 260, 237, 206, 242, 361,
267, 245, 331, 357, 284, 263, 244, 317, 225, 254, 253, 251, 314, 239, 248, 250,
200, 256, 233, 427, 391, 331, 395, 337, 392, 352, 381, 330, 368, 381, 316, 335,
316, 302, 375, 361, 330, 351, 186, 221, 278, 244, 218, 126, 269, 238, 194, 384,
154, 555, 387, 317, 365, 357, 390, 320, 316, 297, 354, 266, 279, 327, 285, 258,
267, 226, 237, 264, 510, 490, 458, 425, 522, 927, 555, 550, 516, 548, 560, 545,
633, 496, 498, 223, 222, 309, 244, 207, 258, 255, 281, 258, 226, 257, 263, 266,
238, 249, 340, 247, 216, 241, 239, 226, 273, 235, 251, 290, 473, 416, 451, 475,
406, 349, 401, 334, 446, 401, 252, 266, 210, 228, 250, 265, 236, 289, 244, 327,
274, 223, 327, 307, 338, 345, 381, 369, 445, 296, 303, 326, 321, 309, 307, 319,
288, 299, 284, 278, 310, 282, 275, 372, 295, 306, 303, 285, 316, 294, 284, 324,
264, 278, 369, 254, 306, 237, 439, 287, 285, 261, 299, 311, 265, 292, 282, 271,
268, 270, 259, 269, 249, 261, 425, 291, 291, 441, 222, 347, 244, 232, 272, 264,
190, 219, 317, 232, 256, 185, 210, 213, 202, 226, 250, 238, 252, 233, 221, 220,
287, 267, 264, 273, 304, 294, 236, 200, 219, 276, 287, 365, 438, 420, 396, 359,
405, 397, 383, 360, 387, 429, 358, 459, 371, 368, 452, 358, 371]

ind = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8,
 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,
 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11,
  11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
   12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14,
   14, 14, 14, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 15, 16, 16, 16, 16,
   16, 17, 17, 17, 17, 17, 18, 18, 18, 18, 18, 19, 19, 19, 19, 19, 20, 20, 20,
   20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 21, 21,
   21, 21, 21, 21, 21, 21, 21, 21, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22,
   22, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 24, 24, 24, 24, 24, 24, 24, 24,
   24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 25, 25, 25, 25, 25, 25, 25,
   25, 25, 25, 25, 25, 25, 26, 26, 26, 26, 26, 27, 27, 27, 27, 27, 28, 28, 28,
   28, 28, 28, 28, 28, 28, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 30, 30, 30,
   30, 30, 30, 31, 31, 31, 31, 31, 32, 32, 32, 32, 32, 33, 34, 34, 34, 34, 34,
   34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34];

#%% Functions ##################################################################
## Jespers HDI:
function hdi(theta_samp,alpha=0.05)
    cred_mass = 1.0-alpha
    ci = zeros(2)
    if length(size(theta_samp))>1   # multidimensional theta
        K,N = size(theta_samp)
        cis = zeros(2,K)
        for k in 1:K    # for each dimension

            ts = theta_samp[k,:]
            sind = sortperm(ts)     # get series of indices that sorts theta values
            sts = ts[sind]  # sorted theta values

            # Shifting a 95% index block from far left to far right and use the
            # samllest interval that can be found during the shift
            N = length(sind)    # number of sample points
            length_ci = Inf     # start with inf HDI interval and shrink
            for i in 1:Int(floor(N*alpha))  # iterate over the 5% lowest values (low to high)
                i2 = Int(floor(N*cred_mass)+i)  # the 5% highest values (low to high), one 95% block away on the other side
                prop_ci = [sts[i],sts[i2]]      # proposed interval
                length_prop_ci = prop_ci[2]-prop_ci[1]
                if length_prop_ci < length_ci
                    ci = prop_ci
                    length_ci = ci[2]-ci[1]
                end
            end
            cis[:,k] = ci

        end
        return cis
    else
        N = length(theta_samp)

        ts = theta_samp
        sind = sortperm(ts)
        sts = ts[sind]

        N = length(sind)
        length_ci = Inf
        for i in 1:Int(floor(N*alpha))
            i2 = Int(floor(N*cred_mass)+i)
            prop_ci = [sts[i],sts[i2]]
            length_prop_ci = prop_ci[2]-prop_ci[1]
            if length_prop_ci < length_ci
                ci = prop_ci
                length_ci = ci[2]-ci[1]
            end
        end
        return ci
    end
end

# A recursive application of the HDI algorithm to find the mode:
function find_mode(sample_pts;runs=2)
    # re-apply the hdi algorithm, until there is only two points or so left
    # so the ci interval becomes smaller and smaller
    # in the end, return an average over the interval
    ci = zeros(2)
    i_left = 0
    i_right = 0

    # Make 50% credibility intervals
    alpha = 0.5;
    cred_mass = 1.0 - alpha;

    # sort the pts:
    sind = sortperm(sample_pts)
    sts = sample_pts[sind]
    N = length(sind)

    if N > runs
        # HDI algorithm, but return indices too:
        length_ci = Inf
        for i in 1:Int(floor(N*alpha))
            i2 = Int(floor(N*cred_mass)+i)
            prop_ci = [sts[i], sts[i2]]
            length_prop_ci = prop_ci[2] - prop_ci[1]
            if length_prop_ci < length_ci
                i_left  = i
                i_right = i2
                ci = prop_ci
                length_ci = ci[2] - ci[1]
            end
        end

        # Recursion:
        ci = find_mode(sts[i_left:i_right],runs=runs)

    else
        ci = [sts[1],sts[end]]
    end

    return ci
end

#%% Logarithmise and Transform Data ############################################
J = maximum(ind)    # number of individuals
I = length(y)       # number of observations

# Go to log-space:
logy = log.(y)

logMean = mean(logy)
logStd = std(logy)

# Mean-centre and scale each element:
zlogy = (logy .- logMean) ./ logStd;

################################################################################
############################# STANS MCMC #######################################

#%% Stan Setup #################################################################

# If run from command line externally:
#projDir = dirname(@__FILE__)

# If run from Jupyter/Hydrogen, maybe change to suit you:
projDir= "/lhome/johhub/Desktop/ABDA/A5"
#projDir= "/lhome/johhub/Desktop/ABDA/A5"
tmpDir = projDir*"/tmp"

noOfChains = 4
N = 10^6 / noOfChains   # more than 10^6 samples make the histograms thin
keepchains = false
burnIn = 10^4

modelString = "
data {
    int<lower=1> J;         // number of individuals for which data was observed
    int<lower=1> I;         // observation number in the total vector
    int<lower=1> ID[I];     // individual ID in the total vector
    real zlogy[I];           // reaction time measurements in the total vector
}
parameters {
    real theta[J];          // mean for each individual
    real<lower=0.000001> sigma;    // same std for all individuals
    real mu;                // mean for the group
    real<lower=0.000001> tau;      // std of the group (only on group = all individuals)
}
transformed parameters { // no transformed variables to use
// how to transform here?
}
model {
    for (i in 1:I)
        zlogy[i] ~ normal(theta[ID[i]],sigma);
    for (j in 1:J)
        theta[j] ~ normal(mu,tau);
    // no prior for sigma and tau is equivalent to a uniform prior
    // however this is an <<improper>> prior and can lead to problems
}
generated quantities {
    real zlogy_pred;
    real theta_pred;
    theta_pred = normal_rng(mu,tau);
    zlogy_pred = normal_rng(theta_pred,sigma);
}
";

#%%### 2) Data for the Stan model, note that variable names must
# correspond to defined model in stan-code
observedData = Dict("I" => I,
                    "J" => J,
                    "zlogy" => zlogy,
                    "ID" => ind);

### 3) Chain specs
myModel = Stanmodel(
                Sample(save_warmup=false,
                       num_warmup=burnIn,
                       num_samples=N,
                       thin=1),   # thin: Period between saved samples
                name = "reactionTime-V1",
                model = modelString,
                printsummary = false,
                tmpdir = tmpDir,
                nchains = noOfChains);    # number of chains (4 default)

#%%### 4) Run CmdStan:
# rc:     return code (0 if all is fine)
# chn:    chain results
# cnames: vector of variable names
rc, chn, cnames = stan(myModel,
                       observedData,
                       projDir,
                       diagnostics = false,
                       CmdStanDir = CMDSTAN_HOME);

#%% Transform back (undo mean-centering and scaling and go to non-log space)
#
# we assumed mean-centered and scaled log-normal distributions for theta
# by our "tricks" above, i.e. μ_y, σ_y fixed
#                    z_logy   ~  Normal(θ, σ^2)
#       (log(y) - m_y) / σ_y  ~  Normal(θ, σ^2)
#       (log(y) - m_y) / σ_y  =  θ + σ * x
# where y is the observation and x the predictor, then
#       log(y) = (θ + σ * x) * σ_y + m_y
#       log(y) = θ*σ_y + σ*σ_y*x  + m_y
#       log(y) = θ*σ_y + m_y + σ*σ_y*x
#       log(y) = θ_trans + σ_trans * x
#       log(y) ~ Normal(θ_trans, σ_trans^2)
#           y  = exp(θ_trans + σ_trans * x)
#           y  ~ LogNormal(θ_trans, σ_trans^2)
#         E[y] = exp(θ_trans + σ_trans^2 / 2)
#
# The same is true for the prior distribution yielding theta:
# ()μ and τ come from hyperdistributions)
#                         θ ~ Normal(μ,τ^2)
#     (θ_trans - m_y) / σ_y = μ + τ * xi
#                   θ_trans = μ*σ_y + m_y + τ*σ_y*xi
#                   θ_trans = μ_trans + τ_trans*xi
#                   θ_trans ~ Normal(μ_trans, τ_trans^2)
#
#                    log(y) = μ_trans + τ_trans*xi + σ_trans*x
# xi and x are both ~ N(0,1)
# a sum of two distributions, i.e. here N(0,τ_trans^2) + N(0,σ_trans^2), are a
# convolution of the two, thus their sum is: N(0,τ_trans^2 + σ_trans^2)
# assuming independence, and therefore
#                    log(y) = μ_trans + sqrt(τ_trans^2 + σ_trans^2)*x
#                    log(y) ~ Normal(μ_trans, (σ_trans^2 + τ_trans^2))
#                      E[y] = exp(μ_trans + (σ_trans^2 + τ_trans^2) / 2)
#                      E[y] = exp(μ_trans + σ_trans^2 / 2 + τ_trans^2 / 2)

#%% check the names and positions of the vars (might change with naming)
chn.value[1,:,1]

#%%
# Originals
# indices in chains:
# 5 = mu
# 7 = sigma
# 9 = tau
# 10 = theta1
# 43 = theta34
# 44 = theta_pred
# 46 = logy_pred

# Throw all the chain samples in a single vetor (vertical concatenation):
θ = 1.0 * chn.value[:,10:(J+10-1),1]
μ = 1.0 * chn.value[:,5,1]
σ = 1.0 * chn.value[:,7,1]
τ = 1.0 * chn.value[:,9,1]
logy_pred = 1.0 * chn.value[:,46,1]
for i in 2:noOfChains
    global θ = vcat(θ, 1.0 * chn.value[:,10:(J+10-1),i])
    global μ = vcat(μ, 1.0 * chn.value[:,5,i])
    global σ = vcat(σ, 1.0 * chn.value[:,7,i])
    global τ = vcat(τ, 1.0 * chn.value[:,9,i])
    global logy_pred = vcat(logy_pred, 1.0 * chn.value[:,46,i])
end

##### RE-READ RESULTS IF CHAIN HAS ALREADY RUN
# only for re-use purposes, please ignore!
# using CSV
#
# N = 10^6
#
# ϕ = Array{Float64,1}(undef,N)
# θ = Array{Float64,2}(undef,N,J)
# μ = Array{Float64,1}(undef,N)
# σ = Array{Float64,1}(undef,N)
# τ = Array{Float64,1}(undef,N)
#
# for i in 1:noOfChains
#     # subindices
#     from = N_chain*(i-1) + 1
#     to = N_chain*i
#
#     # read into a DataFrame:
#     chain_csv = CSV.read(tmpDir*"/reactionTime-V1_samples_$i.csv"; comment="#", normalizenames=true)
#
#     θ_df = chain_csv[:,r"theta"]    # regex
#     for j in 1:J
#         θ[from:to,j] = θ_df[:,j]  # need to bring on array form, from data frame
#     end
#
#     μ[from:to,:] = chain_csv[:, r"mu"][:,1]
#     σ[from:to,:] = chain_csv[:, r"sigma"][:,1]
#     τ[from:to,:] = chain_csv[:, r"tau"][:,1]
# end
######

#%%
# Un-scale and un-mean-centre:
θ_trans = θ .* logStd .+ logMean
#θ_trans = θ
μ_trans = μ .* logStd .+ logMean
#μ_trans = μ
σ_trans = σ .* logStd
#σ_trans = σ
τ_trans = τ .* logStd
#τ_trans = τ
logy_trans = logy_pred .* logStd .+ logMean
#logy_trans = logy_pred

# Get into non-log space:
θ_trans_unLog = exp.(θ_trans .+ 0.5 .* repeat(σ_trans,1,J).^2);
μ_trans_unLog = exp.(μ_trans .+ 0.5 .* σ_trans.^2 .+ 0.5 .* τ_trans.^2);
logy_trans_unLog = exp.(logy_trans);

### Exporting
#using CSV, DataFrames
#exp_τ = Array{Float64,2}(undef,length(τ),2)
#exp_τ[:,1] = τ
#exp_τ[:,2] = τ_trans
#CSV.write(projDir*"/export_tau_A5.csv",  DataFrame(exp_τ), writeheader=false)

#exp_μ = Array{Float64,2}(undef,length(μ),2)
#exp_μ[:,1] = μ
#exp_τ[:,2] = μ_trans
#CSV.write(projDir*"/export_mu_A5.csv",  DataFrame(exp_τ), writeheader=false)

################################################################################
############################# TASKS ############################################

######################## Task A-1 (the Dude)####################################
#%% Mode, Mean, HDIs
dude = θ_trans_unLog[:,4]
ω = mean(find_mode(dude))
μ_bar = mean(dude);
left,right = hdi(dude);

#%% Plotting
histogram(dude, bins=100, normalize=:pdf, label="MCMC", alpha=0.3, linealpha=0.1)   # Comes from StatsPlots now
dPlot = density!(dude,linewidth=3,label="density estimate")
dCurve = filter(!isnan,dPlot.series_list[1].plotattributes[:y])
dTopPoint = maximum(dCurve)
plot!([(left,0),(left,dTopPoint/2)], linewidth=3, color="green", label="HDI",
       annotations = (left, dTopPoint/2, text("$(Int(round(left)))",:green,:bottom)))
plot!([(right,0),(right,dTopPoint/2)],linewidth=3,color="green",label="HDI",
       annotations = (right, dTopPoint/2, text("$(Int(round(right)))",:green,:bottom)))
plot!([(μ_bar,0),(μ_bar,dTopPoint)], linewidth=3, color="blue", label="mean",
       annotations = (μ_bar, 0, text("$(Int(round(μ_bar)))",:blue,:top)))
plot!([(ω,0),(ω,dTopPoint)], linewidth=3, color="red", label="mode",
       annotations = (ω, dTopPoint, text("$(Int(round(ω)))",:red,:bottom)))

Plots.savefig("/home/johhub/Desktop/ABDA/A5/figs/A1-Dude-Stan.pdf")

######################## Task A-2-a-i (expectation new individual) #############
#%% Mode, Mean, HDIs
samplePts = μ_trans_unLog
ω = mean(find_mode(samplePts))
μ_bar = mean(samplePts);
left,right = hdi(samplePts);

#%% Plotting
histogram(samplePts, bins=100, normalize=:pdf, label="MCMC", alpha=0.3, linealpha=0.1)   # Comes from StatsPlots now
dPlot = density!(samplePts,linewidth=3,label="density estimate")
dCurve = filter(!isnan,dPlot.series_list[1].plotattributes[:y])
dTopPoint = maximum(dCurve)
plot!([(left,0),(left,dTopPoint/2)], linewidth=3, color="green", label="HDI",
       annotations = (left, dTopPoint/2, text("$(Int(round(left)))",:green,:bottom)))
plot!([(right,0),(right,dTopPoint/2)],linewidth=3,color="green",label="HDI",
       annotations = (right, dTopPoint/2, text("$(Int(round(right)))",:green,:bottom)))
plot!([(μ_bar,0),(μ_bar,dTopPoint)], linewidth=3, color="blue", label="mean",
       annotations = (μ_bar, 0, text("$(Int(round(μ_bar)))",:blue,:top)))
plot!([(ω,0),(ω,dTopPoint)], linewidth=3, color="red", label="mode",
       annotations = (ω, dTopPoint, text("$(Int(round(ω)))",:red,:bottom)))

Plots.savefig("/home/johhub/Desktop/ABDA/A5/figs/A2-Group-Stan.pdf")

########################  Task A-2-a-ii (new individual prediction)  ###########
samplePts = logy_trans_unLog
ω = mean(find_mode(samplePts))
μ_bar = mean(samplePts);
med = median(samplePts);
left,right = hdi(samplePts);

histogram(samplePts, bins=100, normalize=:pdf, label="Simulation", alpha=0.3, linealpha=0.1)   # Comes from StatsPlots now
dPlot = density!(samplePts,linewidth=3,label="density estimate")
dCurve = filter(!isnan,dPlot.series_list[1].plotattributes[:y])
dTopPoint = maximum(dCurve)
plot!([(left,0),(left,dTopPoint/2)], linewidth=3, color="green", label="HDI",
       annotations = (left, dTopPoint/2, text("$(Int(round(left)))",:green,:bottom)))
plot!([(right,0),(right,dTopPoint/2)],linewidth=3,color="green",label="HDI",
       annotations = (right, dTopPoint/2, text("$(Int(round(right)))",:green,:bottom)))
plot!([(μ_bar,0),(μ_bar,dTopPoint)], linewidth=3, color="blue", label="mean",
       annotations = (μ_bar, 0, text("$(Int(round(μ_bar)))",:blue,:top)))
plot!([(ω,0),(ω,dTopPoint)], linewidth=3, color="red", label="mode",
       annotations = (ω, dTopPoint, text("$(Int(round(ω)))",:red,:bottom)))
plot!([(med,0),(med,dTopPoint)], linewidth=3, color="purple", label="median",
       annotations = (med, 0, text("$(Int(round(med)))",:purple,:bottom)))

Plots.savefig("/home/johhub/Desktop/ABDA/A5/figs/A2-Pred-Stan.pdf")

########################  Task A-2-b (compare to website statistics)  ##########
#%% Taken from the homepage:
median_web = 273.0
mean_web = 284.0
println("Δ Median = ",median_web - med)  # using median from previous task
println("Δ Mean = ",mean_web - μ_bar)    # using mean from previous task

########################  Task A-3 #############################################
# Compare hierarchical theta to individual theta using sample means

#%% Logarithmic sample means
θ_means_log = zeros(J)
for j in 1:J
    θ_means_log[j] = mean(logy[ind .== j])
end

#%% Plot into one figure
# Get ordered indices:
sortIdx = sortperm(θ_means_log)
# Limits for the figure
myXlims = (minimum(θ_trans),maximum(θ_trans))
# Initialise the subplots
StatsPlots.plot(layout=(J, 1),size = (1000, 1500))
# Plot each theta:
for i in 1:(J-1)
    j = sortIdx[i]
    # Sampled thetas and their mean:
    histogram!(θ_trans[:,j], bins=100, normalize=:pdf, legend=false, alpha=0.3, linealpha=0.0,
            ann=(myXlims[1]+.05,4,"ind $j:"),ticks=nothing, yaxis=false, subplot=i, xlims=myXlims)
    vline!([mean(θ_trans[:,j])],linewidth=3, color="black", subplot=i, legend=false)

    # The (log) sample means of the initial data:
    vline!([θ_means_log[j]],linewidth=3, color="red", subplot=i, legend=false)
end

# The last one separately so I can see it in Hydrogen:
j = sortIdx[J]
histogram!(θ_trans[:,j], bins=100, normalize=:pdf, legend=false, alpha=0.3, linealpha=0.0,
            ann=(myXlims[1]+.05,10,"ind $J:"),ticks=nothing, yaxis=false, subplot=J, xlims=myXlims)
vline!([mean(θ_trans[:,j])],linewidth=3, color="black", subplot=J, legend=false)
vline!([θ_means_log[j]],linewidth=3, color="red", subplot=J, legend=false)

Plots.savefig("/home/johhub/Desktop/ABDA/A5/figs/A2-Comp-MLE-Slice-All.pdf")

# --------------- Old Code (please ignore) -------------------------
# Get a color map:
#curColor = get_color_palette(:auto, plot_color(:white), J)
#for n = 1:Int(ceil(J/5))
#    plot()
#    for j in ((n-1)*5+1):min((n*5),J)
#        histogram!(θ_trans[:,j], bins=100, normalize=:pdf,legend=false,alpha=0.2, linealpha=0.0, color=curColor[j])
#        vline!([θ_means_log[j]], linewidth=3, color=curColor[j])
#    end
#    Plots.savefig("/home/johhub/Desktop/ABDA/A5/figs/A2-Comp-MLE-Stan-$n.pdf")
#end


text:
########################################################################################
# Solution
########################################################################################
abstract type AbstractProbODESolution{T,N,S} <: DiffEqBase.AbstractODESolution{T,N,S} end
struct ProbODESolution{T,N,uType,puType,probType,uType2,DType,xType,tType,sType,pType,P,A,S,IType,DE} <: AbstractProbODESolution{T,N,uType}
    u::uType                # array of non-probabilistic function values
    pu::puType              # array of Gaussians
    p::probType             # ODE posterior
    u_analytic::uType2
    errors::DType
    x::xType
    t::tType
    sigmas::sType
    proposals::pType
    prob::P
    alg::A
    solver::S
    dense::Bool
    interp::IType
    retcode::Symbol
    destats::DE
end
function solution_new_retcode(sol::ProbODESolution{T,N}, retcode) where {T,N}
    ProbODESolution{T, N, typeof(sol.u), typeof(sol.pu), typeof(sol.p),
                    typeof(sol.u_analytic), typeof(sol.errors), typeof(sol.x),
                    typeof(sol.t), typeof(sol.sigmas), typeof(sol.proposals),
                    typeof(sol.prob), typeof(sol.alg), typeof(sol.solver),
                    typeof(sol.interp), typeof(sol.destats)}(
                        sol.u, sol.pu, sol.p, sol.u_analytic, sol.errors, sol.x, sol.t,
                        sol.sigmas, sol.proposals, sol.prob, sol.alg, sol.solver, sol.dense,
                        sol.interp, retcode, sol.destats)
end

function DiffEqBase.build_solution(
    prob::DiffEqBase.AbstractODEProblem,
    alg::ODEFilter,
    t, x, sigmas,
    proposals, solver;
    retcode = :Default,
    destats = nothing,
    timeseries_errors = length(x)>2,
    dense = true, dense_errors = dense,
    calculate_error = true,
    kwargs...)
    @unpack d, q, E0 = solver.constants

    @assert length(t) == length(x) == (length(sigmas)+1)

    x = StructArray(x)
    pu = StructArray(map(x -> E0 * x, x))
    u = pu.μ
    p = GaussianODEFilterPosterior(t, x, sigmas, solver)
    interp = p

    T = eltype(eltype(u))
    N = length((size(prob.u0)..., length(u)))

    u_analytic = nothing
    Vector{typeof(prob.u0)}()
    if DiffEqBase.has_analytic(prob.f)
        u_analytic = Vector{typeof(prob.u0)}()
        errors = Dict{Symbol,real(eltype(prob.u0))}()

        sol = ProbODESolution{
            T, N, typeof(u), typeof(pu), typeof(p),
            typeof(u_analytic), typeof(errors),
            typeof(x), typeof(t), typeof(sigmas), typeof(proposals),
            typeof(prob), typeof(alg), typeof(solver), typeof(interp), typeof(destats)}(
                u, pu, p, u_analytic, errors, x, t, sigmas, proposals, prob, alg, solver, dense, interp, retcode, destats)

        if calculate_error
            DiffEqBase.calculate_solution_errors!(sol;
                timeseries_errors=timeseries_errors,
                dense_errors=dense_errors)
        end

        return sol
    else
        return ProbODESolution{
            T, N, typeof(u), typeof(pu), typeof(p), Nothing, Nothing,
            typeof(x), typeof(t), typeof(sigmas), typeof(proposals),
            typeof(prob), typeof(alg), typeof(solver), typeof(interp), typeof(destats)}(
                u, pu, p, nothing, nothing, x, t, sigmas, proposals, prob, alg, solver, dense, interp, retcode, destats)
    end
end

function DiffEqBase.build_solution(sol::ProbODESolution{T,N}, u_analytic, errors) where {T,N}
    return ProbODESolution{
        T, N, typeof(sol.u), typeof(sol.pu), typeof(sol.p), typeof(u_analytic), typeof(errors),
        typeof(sol.x), typeof(sol.t), typeof(sol.sigmas), typeof(sol.proposals),
        typeof(sol.prob), typeof(sol.alg), typeof(sol.solver), typeof(sol.interp), typeof(sol.destats)}(
            sol.u, sol.pu, sol.p, u_analytic, errors, sol.x, sol.t, sol.sigmas, sol.proposals, sol.prob, sol.alg, sol.solver, sol.dense, sol.interp, sol.retcode, sol.destats)
end

########################################################################################
# Dense Output
########################################################################################
abstract type AbstractFilteringPosterior end
struct GaussianFilteringPosterior <: AbstractFilteringPosterior
    t
    x
    sigmas
    solver
end
function (posterior::GaussianFilteringPosterior)(tval::Real)
    @unpack t, x, sigmas, solver = posterior

    @unpack A!, Q!, d, q, E0, Precond, InvPrecond = solver.constants
    @unpack Ah, Qh = solver.cache
    if tval < t[1]
        error("Invalid t<t0")
    end
    if tval in t
        idx = sum(t .<= tval)
        @assert t[idx] == tval
        return x[idx]
    end

    idx = sum(t .<= tval)
    prev_t = t[idx]
    prev_rv = x[idx]
    σ² = sigmas[minimum((idx, end))]

    # Extrapolate
    h1 = tval - prev_t
    P, PI = Precond(h1), InvPrecond(h1)
    A!(Ah, h1)
    Q!(Qh, h1)
    Qh .*= σ²
    goal_pred = predict(P * prev_rv, Ah, Qh, PI)
    goal_pred = PI * goal_pred

    if !solver.smooth || tval >= t[end]
        return goal_pred
    end

    next_t = t[idx+1]
    next_smoothed = x[idx+1]

    # Smooth
    h2 = next_t - tval
    P, PI = Precond(h2), InvPrecond(h2)
    goal_pred = P * goal_pred
    next_smoothed = P * next_smoothed
    A!(Ah, h2)
    Q!(Qh, h2)
    Qh .*= σ²

    goal_smoothed, _ = smooth(goal_pred, next_smoothed, Ah, Qh)

    return PI * goal_smoothed

end
(posterior::GaussianFilteringPosterior)(tvals::AbstractVector) = StructArray(posterior.(tvals))
Base.getindex(post::GaussianFilteringPosterior, key) = post.x[key]

abstract type AbstractODEFilterPosterior <: DiffEqBase.AbstractDiffEqInterpolation end
struct GaussianODEFilterPosterior <: AbstractODEFilterPosterior
    filtering_posterior::GaussianFilteringPosterior
    proj
end

function GaussianODEFilterPosterior(t, x, sigmas, solver)
    GaussianODEFilterPosterior(
        GaussianFilteringPosterior(t, x, sigmas, solver),
        solver.constants.E0)
end
DiffEqBase.interp_summary(interp::GaussianODEFilterPosterior) = "Gaussian ODE Filter Posterior"
Base.getindex(post::GaussianODEFilterPosterior, key) = post.proj * post.filtering_posterior[key]

function (ode_posterior::GaussianODEFilterPosterior)(t::Real)
    return ode_posterior.proj * ode_posterior.filtering_posterior(t)
end
function (ode_posterior::GaussianODEFilterPosterior)(tvec::AbstractVector)
    return StructArray([ode_posterior.proj * ode_posterior.filtering_posterior(t)
                        for t in tvec])
end

(sol::ProbODESolution)(t::Real) = sol.p(t).μ
(sol::ProbODESolution)(t::AbstractVector) = DiffEqBase.DiffEqArray(sol.(t), t)

########################################################################################
# Plotting
########################################################################################
@recipe function f(sol::AbstractProbODESolution; c=1.96)
    times = range(sol.t[1], sol.t[end], length=1000)
    # values = stack(sol.u)
    pus = StructArray(sol.p.(times))
    values = stack(pus.μ)
    # vars = stack(diag.(sol.pu.Σ))
    vars = stack(diag.(pus.Σ))
    stds = sqrt.(vars)
    ribbon --> c * stds
    xguide --> "t"
    yguide --> "y(t)"
    return times, values
end

########################################################################################
# Sampling from a solution
########################################################################################
function get_zero_cross_indices(C)
    bad_idx = []
    D = size(C)[1]
    for i in 1:D
        if all(C[i, :] .< eps(eltype(C))) && all(C[:, i] .< eps(eltype(C)))
            push!(bad_idx, i)
        end
    end
    return bad_idx
end

"""Helper function to sample from our covariances, which often have a "cross" of zeros
For the 0-cov entries the outcome of the sampling is deterministic!"""
function _rand(x::Gaussian, n::Int=1)
    m, C = x.μ, x.Σ

    try
        chol = cholesky(Symmetric(C))
        sample = m .+ chol.L*randn(length(m), n)
        return sample
    catch e
        @warn "Cholesky failed; Try to sample more manually" e
        bad_idx = get_zero_cross_indices(C)

        D = length(m)
        @assert all(C[bad_idx, :] .< eps(eltype(C)))
        @assert all(C[:, bad_idx] .< eps(eltype(C)))
        reduced_idx = setdiff(1:D, bad_idx)
        reduced_x = Gaussian(m[reduced_idx], Symmetric(C[reduced_idx, reduced_idx]))
        reduced_sample = rand(reduced_x)

        sample = reduced_sample
        for i in bad_idx
            insert!(sample, i, m[i])
        end
        @assert sample[bad_idx] == m[bad_idx]

        return sample
    end
end

function sample_back(x_curr::Gaussian, x_next_sample::AbstractVector, Ah::AbstractMatrix, Qh::AbstractMatrix, PI=I)
    m_p, P_p = Ah*x_curr.μ, Ah*x_curr.Σ*Ah' + Qh
    P_p_inv = inv(Symmetric(P_p))
    Gain = x_curr.Σ * Ah' * P_p_inv

    m = x_curr.μ + Gain * (x_next_sample - m_p)

    P = ((I - Gain*Ah) * x_curr.Σ * (I - Gain*Ah)'
         + Gain * Qh * Gain')

    assert_nonnegative_diagonal(P)
    return Gaussian(m, P)
end

function sample(sol::ProbODESolution, n::Int=1)
    sample(sol.t, sol.x, sol.solver, n)
end
function sample(ts, xs, solver, n::Int=1)

    @unpack A!, Q!, d, q, E0, Precond, InvPrecond = solver.constants
    @unpack Ah, Qh = solver.cache
    dim = d*(q+1)

    x = xs[end]
    sample = _rand(x, n)
    @assert size(sample) == (dim, n)

    sample_path = zeros(length(ts), dim, n)
    sample_path[end, :, :] .= sample
    # @info "final value and samples" x.μ sample sample_path[end, :]

    for i in length(xs)-1:-1:1
        dt = ts[i+1] - ts[i]

        i_sigma = sum(solver.times .<= ts[i])
        σ² = solver.sol.sigmas[i_sigma]

        A!(Ah, dt)
        Q!(Qh, dt)
        Qh .*= σ²
        P, PI = Precond(dt), InvPrecond(dt)

        for j in 1:n
            sample_p = P*sample_path[i+1, :, j]
            x_prev_p = P*xs[i]

            prev_sample_p = sample_back(x_prev_p, sample_p, Ah, Qh, PI)

            # sample_path[i, :, j] .= PI*prev_sample_p.μ
            sample_path[i, :, j] .= PI*_rand(prev_sample_p)[:]
        end
    end

    return sample_path[:, 1:d, :]
end
function dense_sample(sol::ProbODESolution, n::Int=1)
    times = range(sol.t[1], sol.t[end], length=1000)
    states = StructArray(sol.p.filtering_posterior(times))

    sample(times, states, sol.solver, n), times
end

max_stars_repo_name: nathanaelbosch/capos

max_issues_repo_name: nathanaelbosch/capos

max_forks_repo_name: nathanaelbosch/capos


text:
# Julia Implementation of BLEU and Smooth BLEU score
# ref: https://github.com/tensorflow/nmt/blob/master/nmt/scripts/bleu.py#L56

# Example: bleu_score([["apple is apple"]], ["apple is appl"])

"""
    get_ngrams(segment, max_order)

Extracts all n-grams upto a given maximum order from an input segment. Returns the counter containing all n-grams upto max_order in segment
with a count of how many times each n-gram occurred.

# Arguments
 - `segment`: text segment from which n-grams will be extracted.
 - `max_order`: maximum length in tokens of the n-grams returned by this methods.

"""
function get_ngrams(segment, max_order)
    ngrams_count = OrderedDict()
    for order in 1:max_order
        for i in 1: (length(segment) - order+1)
            ngram = tuple(segment[i:i+order-1]...)
            if (ngram) in keys(ngrams_count)
                ngrams_count[ngram] += 1
            else
                ngrams_count[ngram] = 1
            end
        end
    end
    return ngrams_count
end

"""
    bleu_score(reference_corpus, translation_corpus; max_order=4, smooth=false)

Computes BLEU score of translated segments against one or more references. Returns the `BLEU score`, `n-gram precisions`, `brevity penalty`,
geometric mean of n-gram precisions, translation_length and  reference_length

# Arguments
 - `reference_corpus`: list of lists of references for each translation. Each reference should be tokenized into a list of tokens.
 - `translation_corpus`: list of translations to score. Each translation should be tokenized into a list of tokens.
 - `max_order`: maximum n-gram order to use when computing BLEU score.
 - `smooth=false`: whether or not to apply. Lin et al. 2004 smoothing.

"""
function bleu_score(reference_corpus, translation_corpus; max_order=4, smooth=false)
    matches_by_order = zeros(max_order)
    possible_matches_by_order = zeros(max_order)
    reference_length = 0
    translation_length = 0
    for (references, translation) in zip(reference_corpus, translation_corpus)
        reference_length += min([length(r) for r in references]...)
        translation_length += length(translation)
        merged_ref_ngram_counts = OrderedDict()
        for reference in references
          ref_ngrams = get_ngrams(reference, max_order)
          keys_union = union(keys(merged_ref_ngram_counts), keys(ref_ngrams))
          for key in keys_union
              try (b[key])
                  try (ref_ngrams[key])
                      merged_ref_ngram_counts[key] = max(merged_ref_ngram_counts[key], ref_ngrams[i])
                  catch error
                      continue
                  end
              catch error
                  merged_ref_ngram_counts[key] = ref_ngrams[key]
              end
           end
        end
        # print(length(merged_ref_ngram_counts),"\n")
        translation_ngram_counts = get_ngrams(translation, max_order)
        overlap = OrderedDict()
        keys_union = intersect(keys(merged_ref_ngram_counts), keys(translation_ngram_counts))
        for key in keys_union
                 overlap[key] = min(translation_ngram_counts[key],  merged_ref_ngram_counts[key])
        end
        for key in overlap
            matches_by_order[length(key[1])] += key[2]
        end
        for order in 1:max_order
            possible_matches = length(translation) - order + 1
            if possible_matches > 0
                possible_matches_by_order[order] += possible_matches
            end
        end
    end
    precisions = zeros(max_order)
    for i in 1:max_order
        if smooth
            precisions[i] = (matches_by_order[i] + 1.0) / (possible_matches_by_order[i] + 1.0)
        else
            if possible_matches_by_order[i]>0
                precisions[i] = (float(matches_by_order[i]) / possible_matches_by_order[i])
            end
        end
    end

    geo_mean = 0.0
    if min(precisions...) > 0
       p_log_sum = sum(log.(precisions)) / max_order
       geo_mean = exp(p_log_sum)
    end

    ratio =  translation_length / reference_length
    bp = 1.0
    if ratio <1.0
       bp = exp(1 - 1 /ratio)
    end

    bleu = geo_mean * bp
    return bleu, precisions, bp, geo_mean, translation_length, reference_length
end


text:
#ДАНО: Где-то на неограниченном со всех сторон поле и без внутренних перегородок имеется единственный маркер. Робот - в произвольной клетке поля.
#РЕЗУЛЬТАТ: Робот - в клетке с тем маркером.
function find_marker(r)
    steps_max=1
    side=Nord
    while ismarker(r)==0
        for _ in 1:2
            find_marker!(r,side,steps_max)
            side=next(side)
        end
        steps_max+=1
    end
end
function find_marker!(r,side,steps_max)
    for _ in 1:steps_max
        if ismarker(r)
            return nothing
        end
        move!(r,side)
    end
end
next(side)=HorizonSide(mod(Int(side)+1,4))
find_marker(r)


text:
using JLD2, ProgressMeter, DSP, Statistics
@load "dat/anonomous_data.jld"
Δ = 1/1.28
fitstart = 85
fitend = length(waveseries)-1

##
include("src/exampleanalysis/PaperAnalysis.jl")

## compute spectrogram, errorgram and mean direction for threshold choice
errorgram, Ω, spectime = PaperAnalysis.compute_errorgram(waveseries, wavetime)
wavedir = mod2pi.(PaperAnalysis.mt_wavedirection(permutedims(reduce(vcat, waveseries)), 2304, 0).power) # uses same Ω as above

##
"Compute the lower cutoff threshold from the errorgram."
function errorgramthreshold(R, Ω, tol, window)
    first_good = findprev(j->abs(mean(R[i] for i in j-window:j+window))>tol, window+1:length(R)-window, length(R)÷2)
    return isnothing(first_good) ? error("Failed to find cutoff") : Ω[first_good + window]
end

"Compute the threshold to remove the swell. Hard coded for the specific swell in this example."
function wavedir_threshold(dir,Ω)
    thresh = findprev(dir.>5,length(dir)÷2)
    return isnothing(thresh) ? NaN : Ω[thresh]
end

## choose threholds
swellend = fitend-50 # final record which actually contains swell
errorthreshold = [errorgramthreshold(a,Ω, 2, 20) for a in eachcol(errorgram[:,fitstart:fitend])]
swellthreshold = [wavedir_threshold(w,Ω) for w in eachcol(wavedir[:,fitstart:swellend])]
threshold = max.(errorthreshold,[swellthreshold;zeros(length(errorthreshold)-length(swellthreshold))])

## initialise fits
ini = [PaperAnalysis.InitialiseData(PaperAnalysis.WindSea,waveseries[i],Δ,threshold[k]) for (k,i) ∈ enumerate(fitstart:fitend)]

## fit models
fits = @showprogress [PaperAnalysis.FittedModel(i) for i ∈ ini]
@save "dat/wavefits.jld" fits fitstart fitend


text:
function pcap_equation(cur_reactor::AbstractReactor)
  cur_G_T = K_PC(cur_reactor)

  cur_G_T *= cur_reactor.sigma_v

  cur_equation = Equation(1, 0, -2, cur_G_T)

  cur_equation
end