text:
# The references to special cases are to Table of Integrals, Series, and Products, § 9.121, followed by NIST's DLMF.

"""
Compute the Gauss hypergeometric function `₂F₁(a, b; c; z)`.
"""
function _₂F₁(a, b, c, z; method::Symbol = :general)
    if real(b) < real(a)
        return _₂F₁(b, a, c, z; method = method) # ensure a ≤ b
    elseif isequal(a, c) # 1. 15.4.6
        return exp(-b*log1p(-z))
    elseif isequal(b, c) # 1. 15.4.6
        return exp(-a*log1p(-z))
    elseif isequal(c, 0.5)
        if isequal(a+b, 0) # 31. 15.4.11 & 15.4.12
            return cosnasinsqrt(2b, z)
        elseif isequal(a+b, 1) # 32. 15.4.13 & 15.4.14
            return cosnasinsqrt(1-2b, z)*exp(-log1p(-z)/2)
        elseif isequal(b-a, 0.5) # 15.4.7 & 15.4.8
            return expnlog1pcoshatanhsqrt(-2a, z)
        end
    elseif isequal(c, 1.5)
        if abeqcd(a, b, 0.5) # 13. 15.4.4 & 15.4.5
            return sqrtasinsqrt(z)
        elseif abeqcd(a, b, 1) # 14.
            return sqrtasinsqrt(z)*exp(-log1p(-z)/2)
        elseif abeqcd(a, b, 0.5, 1) # 15. 15.4.2 & 15.4.3
            return sqrtatanhsqrt(z)
        elseif isequal(a+b, 1) # 29. 15.4.15 & 15.4.16
            return sinnasinsqrt(1-2b, z)
        elseif isequal(a+b, 2) # 30.
            return sinnasinsqrt(2-2b, z)*exp(-log1p(-z)/2)
        elseif isequal(b-a, 0.5) # 4. 15.4.9 & 15.4.10
            return expnlog1psinhatanhsqrt(1-2a, z)
        end
    elseif isequal(c, 2)
        if abeqcd(a, b, 1) # 6. 15.4.1
            return log1pover(-z)
        elseif a ∈ ℤ && b == 1 # 5.
            return expm1nlog1p(1-a, -z)
        elseif a == 1 && b ∈ ℤ # 5.
            return expm1nlog1p(1-b, -z)
        end
    elseif isequal(c, 4)
        if abeqcd(a, b, 2)
            return logandpoly(z)
        end
    elseif isequal(c, 2.5) && abeqcd(a, b, 1, 1.5)
        return speciallog(z)
    end
    if method == :positive
        return _₂F₁positive(a, b, c, z)
    else#if method == :general
        return _₂F₁general(a, b, c, z) # catch-all
    end
end

"""
Compute the Gauss hypergeometric function `₂F₁(a, b; c; z)` with positive parameters a, b, and c and argument 0 ≤ z ≤ 1. Useful for statisticians.
"""
function _₂F₁positive(a, b, c, z)
    @assert a > 0 && b > 0 && c > 0 && 0 ≤ z ≤ 1
    if z == 1
        return _₂F₁argument_unity(a, b, c, float(z))
    else
        return _₂F₁maclaurin(a, b, c, z)
    end
end

"""
Compute the Gauss hypergeometric function `₂F₁(a, b; c; z)` with general parameters a, b, and c.
This polyalgorithm is designed based on the paper

N. Michel and M. V. Stoitsov, Fast computation of the Gauss hypergeometric function with all its parameters complex with application to the Pöschl–Teller–Ginocchio potential wave functions, Comp. Phys. Commun., 178:535–551, 2008.
"""
function _₂F₁general(a, b, c, z)
    if z == 1
        return _₂F₁argument_unity(a, b, c, float(z))
    elseif real(b) < real(a)
        return _₂F₁general(b, a, c, z)
    elseif !isalmostwellpoised(a, b, c)
        return exp((c-a-b)*log1p(-z))*_₂F₁general(c-a, c-b, c, z)
    elseif abs(z) ≤ ρ || -a ∈ ℕ₀ || -b ∈ ℕ₀
        return _₂F₁maclaurin(a, b, c, z)
    elseif abs(z/(z-1)) ≤ ρ
        return exp(-a*log1p(-z))_₂F₁maclaurin(a, c-b, c, z/(z-1))
    elseif abs(inv(z)) ≤ ρ
        return _₂F₁Inf(a, b, c, z)
    elseif abs(1-inv(z)) ≤ ρ
        return exp(-a*log1p(-z))*_₂F₁Inf(a, c-b, c, z/(z-1))
    elseif abs(1-z) ≤ ρ
        return _₂F₁one(a, b, c, z)
    elseif abs(inv(1-z)) ≤ ρ
        return exp(-a*log1p(-z))*_₂F₁one(a, c-b, c, z/(z-1))
    else
        return pFqweniger([a, b], [c], z) # _₂F₁taylor(a, b, c, z)
    end
end

"""
Compute the Gauss hypergeometric function `₂F₁(a, b; c; z)` with general parameters a, b, and c.
This polyalgorithm is designed based on the review

J. W. Pearson, S. Olver and M. A. Porter, Numerical methos for the computation of the confluent and Gauss hypergeometric functions, arXiv:1407.7786, 2014.
"""
function _₂F₁general2(a, b, c, z)
    T = promote_type(typeof(a), typeof(b), typeof(c), typeof(z))
    if z == 1
        return _₂F₁argument_unity(a, b, c, float(z))
    elseif real(b) < real(a)
        return _₂F₁general(b, a, c, z)
    elseif !isalmostwellpoised(a, b, c)
        return exp((c-a-b)*log1p(-z))*_₂F₁general2(c-a, c-b, c, z)
    elseif abs(z) ≤ ρ || -a ∈ ℕ₀ || -b ∈ ℕ₀
        return _₂F₁maclaurin(a, b, c, z)
    elseif abs(z / (z - 1)) ≤ ρ && absarg(1 - z) < convert(real(T), π) # 15.8.1
        w = z/(z-1)
        return _₂F₁maclaurin(a, c-b, c, w)*exp(-a*log1p(-z))
    elseif abs(inv(z)) ≤ ρ && absarg(-z) < convert(real(T), π)
        w = inv(z)
        if isapprox(a, b) # 15.8.8
            return gamma(c)/gamma(a)/gamma(c-a)*(-w)^a*_₂F₁logsumalt(a, c-a, z, w)
        elseif a-b ∉ ℤ # 15.8.2
            return gamma(c)*((-w)^a*gamma(b-a)/gamma(b)/gamma(c-a)*_₂F₁maclaurin(a, a-c+1, a-b+1, w)+(-w)^b*gamma(a-b)/gamma(a)/gamma(c-b)*_₂F₁maclaurin(b, b-c+1, b-a+1, w))
        end # TODO: full 15.8.8
    elseif abs(inv(1-z)) ≤ ρ && absarg(-z) < convert(real(T), π)
        w = inv(1-z)
        if isapprox(a, b) # 15.8.9
            return gamma(c)*exp(-a*log1p(-z))/gamma(a)/gamma(c-b)*_₂F₁logsum(a, c-a, z, w, 1)
        elseif a-b ∉ ℤ # 15.8.3
            return gamma(c)*(exp(-a*log1p(-z))*gamma(b-a)/gamma(b)/gamma(c-a)*_₂F₁maclaurin(a, c-b, a-b+1, w)+exp(-b*log1p(-z))*gamma(a-b)/gamma(a)/gamma(c-b)*_₂F₁maclaurin(b, c-a, b-a+1, w))
        end # TODO: full 15.8.9
    elseif abs(1-z) ≤ ρ && absarg(z) < convert(real(T), π) && absarg(1-z) < convert(real(T), π)
        w = 1-z
        if isapprox(c, a + b) # 15.8.10
            return gamma(c)/gamma(a)/gamma(b)*_₂F₁logsum(a, b, z, w, -1)
        elseif c - a - b ∉ ℤ # 15.8.4
            return gamma(c)*(gamma(c-a-b)/gamma(c-a)/gamma(c-b)*_₂F₁maclaurin(a, b, a+b-c+1, w)+exp((c-a-b)*log1p(-z))*gamma(a+b-c)/gamma(a)/gamma(b)*_₂F₁maclaurin(c-a, c-b, c-a-b+1, w))
        end # TODO: full 15.8.10
    elseif abs(1-inv(z)) ≤ ρ && absarg(z) < convert(real(T), π) && absarg(1-z) < convert(real(T), π)
        w = 1-inv(z)
        if isapprox(c, a + b) # 15.8.11
            return gamma(c)*z^(-a)/gamma(a)*_₂F₁logsumalt(a, b, z, w)
        elseif c - a - b ∉ ℤ # 15.8.5
            return gamma(c)*(z^(-a)*gamma(c-a-b)/gamma(c-a)/gamma(c-b)*_₂F₁maclaurin(a, a-c+1, a+b-c+1, w)+z^(a-c)*(1-z)^(c-a-b)*gamma(a+b-c)/gamma(a)/gamma(b)*_₂F₁maclaurin(c-a, 1-a, c-a-b+1, w))
        end # TODO: full 15.8.11
    elseif abs(z-0.5) > 0.5
        if isapprox(a, b) && !isapprox(c, a+0.5)
            return gamma(c)/gamma(a)/gamma(c-a)*(0.5-z)^(-a)*_₂F₁continuationalt(a, c, 0.5, z)
        elseif a-b ∉ ℤ
            return gamma(c)*(gamma(b-a)/gamma(b)/gamma(c-a)*(0.5-z)^(-a)*_₂F₁continuation(a, a+b, c, 0.5, z) + gamma(a-b)/gamma(a)/gamma(c-b)*(0.5-z)^(-b)*_₂F₁continuation(b, a+b, c, 0.5, z))
        end
    end
    return pFqweniger([a, b], [c], z)
end

# Special case of (-x)^a*_₂F₁ to handle LogNumber correctly in RiemannHilbert.jl
function mxa_₂F₁(a, b, c, z)
    if isequal(c, 2)
        if abeqcd(a, b, 1) # 6. 15.4.1
            return log1p(-z)
        end
    elseif isequal(c, 4)
        if abeqcd(a, b, 2)
            return 6*(-2 + (1-2/z)*log1p(-z))
        end
    end
    return (-z)^a*_₂F₁(a, b, c, z)
end

function _₂F₁argument_unity(a, b, c, z)
    T = promote_type(typeof(a), typeof(b), typeof(c), typeof(z))
    a, b, c = T(a), T(b), T(c)
    if iswellpoised(a, b, c)
        return exp(loggamma(c)-loggamma(c-a)+loggamma(c-a-b)-loggamma(c-b))
    elseif isalmostwellpoised(a, b, c)
        if a + b == c
            f = loggamma(c)-loggamma(a)-loggamma(b)
            if isreal(f)
                return T(Inf)
            else
                return T(Inf)*sign(exp(f))
            end
        else
            return T(NaN)
        end
    else # !isalmostwellpoised(a, b, c)
        f = loggamma(c)-loggamma(a)+gamma(a+b-c)-loggamma(b)
        if isreal(f)
            return T(Inf)
        else
            return T(Inf)*sign(exp(f))
        end
    end
end


text:
using POMCP
using Base.Test

using POMDPModels
using POMDPs
using POMDPToolbox
using ParticleFilters

using NBInclude

pc = ParticleCollection([1,2])

solver = POMCPSolver(estimate_value=RolloutEstimator(FeedWhenCrying()),
                    eps=0.01,
                    c=10.0,
                    tree_queries=50,
                    rng=MersenneTwister(2))

test_solver(solver, BabyPOMDP())

solver = POMCPSolver()
test_solver(solver, BabyPOMDP())

# test FORollout
solver = POMCPSolver(estimate_value=FORollout(FeedWhenCrying()),
                    eps=0.01,
                    c=10.0,
                    tree_queries=50,
                    rng=MersenneTwister(2))

test_solver(solver, BabyPOMDP())

# test for particle depletion
solver = POMCPSolver(estimate_value=RolloutEstimator(FeedWhenCrying()),
                     eps=0.01,
                     c=10.0,
                     tree_queries=5,
                     rng=MersenneTwister(2))

@test_throws ErrorException test_solver(solver, BabyPOMDP(), max_steps=100)

# test for DPW
solver = POMCPDPWSolver(tree_queries=100)

test_solver(solver,BabyPOMDP())

# test for enable/disable action pw
solver = POMCPDPWSolver(tree_queries=100,
                     eps=0.01,
                     c=10.0,
                     enable_action_pw=false,
                     rng=MersenneTwister(2))
test_solver(solver, BabyPOMDP())

test_solver(solver, LightDark1D(), updater=SIRParticleFilter(LightDark1D(), 100))
# test_solver(solver, TigerPOMDP(), updater=DiscreteUpdater

include("visualization.jl")
nbinclude("../notebooks/Display_Tree.ipynb")

nbinclude("../notebooks/Basic_Usage.ipynb")

nbinclude("../notebooks/Minimal_Example.ipynb")

include("Belief_and_Particle_Filter_Options.jl")

ext: jl


text:
using Immerse, Colors

hfig = figure()
x = range(0, stop=4pi, length=101)
display(plot(x=x, y=sin.(x), Geom.line(tag=:line)))

val = Immerse.setproperty!((hfig,:line), rand(1:5), :linewidth)
@test Immerse.getproperty((hfig,:line), :linewidth) == [val*mm]
val = Immerse.setproperty!((hfig,:line), RGB(rand(),rand(),rand()), :stroke)
@test Immerse.getproperty((hfig,:line), :stroke) == [coloralpha(val)]
Immerse.setproperty!((hfig,:line), false, :visible)
@test Immerse.getproperty((hfig,:line), :visible) == [false]
Immerse.setproperty!((hfig,:line), true, :visible)
@test Immerse.getproperty((hfig,:line), :visible) == [true]


text:
∑ = sum
using StatsFuns:logistic
sigmoid = logistic
using CSV:File
using DataFrames
using Printf:@sprintf

# DataTypes
using LinearAlgebra:Transpose
Numeric = Union{Int64,Float64}
NumericV = Union{Array{Int64,1},Array{Float64,1}}
NumericM = Union{Array{Int64,2},Array{Float64,2},Transpose{Float64,Array{Float64,2}},Transpose{Int64,Array{Int64,2}}}
arrayNumericV = Union{Array{Array{Int64,1},1},Array{Array{Float64,1},1}}
TupleNumericM = Tuple{NumericM, NumericM};

macro name(arg)
    string(arg)
end

function print_arr_dims(array::NumericV)::String
    arrayName = @name(array)
    outStr = @sprintf "FYI: %s is a numeric vector " arrayName
    outStr *= @sprintf "with dimensions %i × 1." size(array)[1]
    println(outStr)
end

function print_arr_dims(array::NumericM)::String
    arrayName = @name(array)
    outStr = @sprintf "FYI: %s is a numeric matrix " arrayName
    #outStr = @sprintf "with dimensions %i × %i." size(array)[1] size(array)[2]
    println(outStr)
end

max_stars_repo_name: luis-perez-one/courseStandfordML

max_issues_repo_name: luis-perez-one/courseStandfordML

max_forks_repo_name: luis-perez-one/courseStandfordML


text:
using Plots
using Formatting

function logistic(r,x)
    r*x*(1.0-x)
end

function logisticmap(;x0=0.21,r=4,n=10)
    results = [(x0,logistic(r,x0))];
    #N = 100; Note this line should not be here
    for i in 1:(n-1) # Note change of case N -> n.  Why?
      xold = results[end][2]
      xnew = logistic(r,xold)
      append!(results,[(xold,xold)])
      append!(results,[(xold,xnew)])
  end
  results
end

function vary(dname;x0_vals=[],r_vals=[],n_vals=[])
 if !isdir(dname)
  mkdir(dname)
 end
 fno = 1
 for x0 in x0_vals
  for r in r_vals
   for n in n_vals
   fname = format(dname*"/frame{:05d}.png",fno)
   fno = fno + 1
   simres = logisticmap(x0=x0, r=r, n=n)
   xvals = [x for (x, y) in simres]
   yvals = [y for (x, y) in simres]
   plot(xvals,yvals,legend=false,xaxis=false,yaxis=false,xlim=(0,1),ylim=(0,1));
   savefig(fname)
  end
 end
end

# To make frames varying r, specify x0_vals as an array with a single value
# along with whatever values of r you want to loop through
vary("frames_r";x0_vals=[0.5],r_vals=0:0.1:5)

vary("loop";x0_vals=[(r-1)/r],r_vals=[4])
# To make frames varying x0, specify r_vals as an array with a single value
# along with whatever values of x0 you want to loop through
vary("frames_x0";x0_vals=0:0.05:1,r_vals=[4.0])

for x0 in 0:0.0001:1
    global fno
    fname = format("framesxl/frame{:05d}.png",fno)
    fno = fno + 1
    simres = logisticmap(x0=x0, r=4, n=100)
    xvals = [x for (x, y) in simres]
    yvals = [y for (x, y) in simres]
    plot(xvals,yvals,legend=false,xaxis=false,yaxis=false,xlim=(0,1),ylim=(0,1));
    savefig(fname)
end

#if !isdir("frames_r")
 # mkdir("frames_r")
#end
function graphpngs(x0, delta_r, dname)
    #dictname=Dict()
    #dictname["dname"] = dname
    if !isdir(dname)
        mkdir(dname)
    end
    fno = 1
    global fno
    for r in 0:delta_r:5
        fname = format(dname*"/frame{:05d}.png",fno)
        fno = fno + 1
        simres = logisticmap(x0=x0, r=r, n=100)
        xvals = [x for (x, y) in simres]
        yvals = [y for (x, y) in simres]
        plot(xvals,yvals,legend=false,xaxis=false,yaxis=false,xlim=(0,1),ylim=(0,1));
        savefig(fname)
    end
end

function varyx0(delta_x0, r, dname)
    if !isdir(dname)
        mkdir(dname)
    end
    fno = 1
    global fno
    for x0 in 0:delta_x0:1
        fname = format(dname*"/frame{:05d}.png",fno)
        fno = fno + 1
        simres = logisticmap(x0=x0, r=r, n=100)
        xvals = [x for (x, y) in simres]
        yvals = [y for (x, y) in simres]
        plot(xvals,yvals,legend=false,xaxis=false,yaxis=false,xlim=(0,1),ylim=(0,1));
        savefig(fname)
    end
end

function varyr(x0, delta_r, dname)
     if !isdir(dname)
         mkdir(dname)
     end
     fno = 1
     global fno
     for r in 0:delta_r:5
         fname = format(dname*"/frame{:05d}.png",fno)
         fno = fno + 1
         simres = logisticmap(x0=x0, r=r, n=100)
         xvals = [x for (x, y) in simres]
         yvals = [y for (x, y) in simres]
         plot(xvals,yvals,legend=false,xaxis=false,yaxis=false,xlim=(0,1),ylim=(0,1));
         savefig(fname)
     end
end

function varyxr(x0, r, dname, v)
    if v==1
        varyx0(x0, r, dname)
    elseif v==2
        varyr(x0, r, dname)
    else
        println("To vary x0 set v to 1, to vary r set v to 2")
    end
end

function varyxr2(delta_x0, delta_r, dname, v)
    if v==1 #varys x0
        if !isdir(dname)
            mkdir(dname)
        end
        fno = 1
        global fno
        for x0 in 0:delta_x0:1
            fname = format(dname*"/frame{:05d}.png",fno)
            fno = fno + 1
            r = delta_r
            simres = logisticmap(x0=x0, r=r, n=100)
            xvals = [x for (x, y) in simres]
            yvals = [y for (x, y) in simres]
            plot(xvals,yvals,legend=false,xaxis=false,yaxis=false,xlim=(0,1),ylim=(0,1));
            savefig(fname)
        end
    elseif v==2 #varys r
        if !isdir(dname)
            mkdir(dname)
        end
        fno = 1
        global fno
        for r in 0:delta_r:5
            fname = format(dname*"/frame{:05d}.png",fno)
            fno = fno + 1
            x0=delta_x0
            simres = logisticmap(x0=x0, r=r, n=100)
            xvals = [x for (x, y) in simres]
            yvals = [y for (x, y) in simres]
            plot(xvals,yvals,legend=false,xaxis=false,yaxis=false,xlim=(0,1),ylim=(0,1));
            savefig(fname)
        end
    else
        println("To vary x0 set v to 1, to vary r set v to 2.")
    end
end

function graphpngs(x0, delta_r, dname)
    #dictname=Dict()
    #dictname["dname"] = dname
    if !isdir(dname)
        mkdir(dname)
    end
    fno = 1
    global fno
    for r in 0:delta_r:5
        fname = format(dname*"/frame{:05d}.png",fno)
        fno = fno + 1
        simres = logisticmap(x0=x0, r=r, n=100)
        xvals = [x for (x, y) in simres]
        yvals = [y for (x, y) in simres]
        plot(xvals,yvals,legend=false,xaxis=false,yaxis=false,xlim=(0,1),ylim=(0,1));
        savefig(fname)
    end
end

if !isdir("frames_time")
  mkdir("frames_time")
end

@time begin
    fno = 1 # Frame counter
    for x0 in 0:1:1.0
    #for r in 0:0.01:5
        # Make a filename (within new directory) that is numbered
        # padded to 5 digits with zeros (i.e. 1 -> 00001)
        global fno
        fname = format("frames_time/frame{:05d}.png",fno)
        fno = fno + 1

        simres = logisticmap(x0=x0,r=4,n=100)
        #simres = logisticmap(x0=0.2,r=r,n=100)
        xvals = [x for (x, y) in simres]
        yvals = [y for (x, y) in simres]
        plot(xvals,yvals,legend=false,xaxis=false,yaxis=false,xlim=(0,1),ylim=(0,1));
        savefig(fname)
    end
end

# NOTES
# These are cobweb plots
# https://en.wikipedia.org/wiki/Cobweb_plot

# The above code generates cobweb plots, fixing r and varying x0
# for all x0 from 0 to 1.0 in steps of 0.1

# Can you instead fix x0 and vary r from 0 to 5?

# Can you check how long it takes your computer to run through
# one iteration of the loop above?  How many frames could you
# make in an hour? - 18000  Try decreasing the step size for varying x0
# so that your computer is busy making plots for one hour.
# Do the same for the r version.

# To make a .gif from the .png files in the frames directory
# Need to first install ffmpeg, add to system path, then navigate to frames directory in command line
# and run ffmpeg like this:
# ffmpeg -i frame%05d.png -vf "fps=12,split[s0][s1];[s0]palettegen[p];[s1][p]paletteuse" -loop 0 output.gif

#

#xvals = []
#yvals = []

#for (x, y) in simres
    #append!(xvals, x)
    #append!(yvals, y)
#end


text:
push!(LOAD_PATH,"../src/")

using Documenter, Plasm

makedocs(
	format = :html,
	sitename = "Plasm.jl",
	pages = [
		"Home" => "index.md",
		"Visualization" => "larview.md",
		"Numbering" => "numbering.md",
		"Glossary" => "glossary.md"
	]
)


text:
const GAEVF = JuMP.GenericAffExpr{Float64, VariableRef}
const GQEVF = JuMP.GenericQuadExpr{Float64, VariableRef}

function moi_tests(
    model::DecisionModel,
    params::Bool,
    vars::Int,
    interval::Int,
    lessthan::Int,
    greaterthan::Int,
    equalto::Int,
    binary::Bool,
)
    JuMPmodel = PSI.get_jump_model(model)
    @test (:ParameterJuMP in keys(JuMPmodel.ext)) == params
    @test JuMP.num_variables(JuMPmodel) == vars
    @test JuMP.num_constraints(JuMPmodel, GAEVF, MOI.Interval{Float64}) == interval
    @test JuMP.num_constraints(JuMPmodel, GAEVF, MOI.LessThan{Float64}) == lessthan
    @test JuMP.num_constraints(JuMPmodel, GAEVF, MOI.GreaterThan{Float64}) == greaterthan
    @test JuMP.num_constraints(JuMPmodel, GAEVF, MOI.EqualTo{Float64}) == equalto
    @test ((JuMP.VariableRef, MOI.ZeroOne) in JuMP.list_of_constraint_types(JuMPmodel)) ==
          binary

    return
end

function psi_constraint_test(
    model::DecisionModel,
    constraint_keys::Vector{<:PSI.ConstraintKey},
)
    constraints = PSI.get_constraints(model)
    for con in constraint_keys
        @test get(constraints, con, nothing) !== nothing
    end
    return
end

function psi_aux_variable_test(
    model::DecisionModel,
    constraint_keys::Vector{<:PSI.AuxVarKey},
)
    op_container = PSI.get_optimization_container(model)
    vars = PSI.get_aux_variables(op_container)
    for key in constraint_keys
        @test get(vars, key, nothing) !== nothing
    end
    return
end

function psi_checkbinvar_test(
    model::DecisionModel,
    bin_variable_keys::Vector{<:PSI.VariableKey},
)
    container = PSI.get_optimization_container(model)
    for variable in bin_variable_keys
        for v in PSI.get_variable(container, variable)
            @test JuMP.is_binary(v)
        end
    end
    return
end

function psi_checkobjfun_test(model::DecisionModel, exp_type)
    model = PSI.get_jump_model(model)
    @test JuMP.objective_function_type(model) == exp_type
    return
end

function moi_lbvalue_test(model::DecisionModel, con_key::PSI.ConstraintKey, value::Number)
    for con in PSI.get_constraints(model)[con_key]
        @test JuMP.constraint_object(con).set.lower == value
    end
    return
end

function psi_checksolve_test(model::DecisionModel, status)
    model = PSI.get_jump_model(model)
    JuMP.optimize!(model)
    @test termination_status(model) in status
end

function psi_checksolve_test(model::DecisionModel, status, expected_result, tol=0.0)
    res = solve!(model)
    model = PSI.get_jump_model(model)
    @test termination_status(model) in status
    obj_value = JuMP.objective_value(model)
    @test isapprox(obj_value, expected_result, atol=tol)
end

function psi_ptdf_lmps(res::ProblemResults, ptdf)
    cp_duals = read_dual(res, PSI.ConstraintKey(CopperPlateBalanceConstraint, PSY.System))
    λ = Matrix{Float64}(cp_duals[:, propertynames(cp_duals) .!= :DateTime])

    flow_duals = read_dual(res, PSI.ConstraintKey(NetworkFlowConstraint, PSY.Line))
    μ = Matrix{Float64}(flow_duals[:, ptdf.axes[1]])

    buses = get_components(Bus, get_system(res))
    lmps = OrderedDict()
    for bus in buses
        lmps[get_name(bus)] = μ * ptdf[:, get_number(bus)]
    end
    lmp = λ .+ DataFrames.DataFrame(lmps)
    return lmp[!, sort(propertynames(lmp))]
end

function check_variable_unbounded(
    model::DecisionModel,
    ::Type{T},
    ::Type{U},
) where {T <: PSI.VariableType, U <: PSY.Component}
    return check_variable_unbounded(model::DecisionModel, PSI.VariableKey(T, U))
end

function check_variable_unbounded(model::DecisionModel, var_key::PSI.VariableKey)
    psi_cont = PSI.get_optimization_container(model)
    variable = PSI.get_variable(psi_cont, var_key)
    for var in variable
        if JuMP.has_lower_bound(var) || JuMP.has_upper_bound(var)
            return false
        end
    end
    return true
end

function check_variable_bounded(
    model::DecisionModel,
    ::Type{T},
    ::Type{U},
) where {T <: PSI.VariableType, U <: PSY.Component}
    return check_variable_bounded(model, PSI.VariableKey(T, U))
end

function check_variable_bounded(model::DecisionModel, var_key::PSI.VariableKey)
    psi_cont = PSI.get_optimization_container(model)
    variable = PSI.get_variable(psi_cont, var_key)
    for var in variable
        if !JuMP.has_lower_bound(var) || !JuMP.has_upper_bound(var)
            return false
        end
    end
    return true
end

function check_flow_variable_values(
    model::DecisionModel,
    ::Type{T},
    ::Type{U},
    device_name::String,
    limit::Float64,
) where {T <: PSI.VariableType, U <: PSY.Component}
    psi_cont = PSI.get_optimization_container(model)
    variable = PSI.get_variable(psi_cont, T(), U)
    for var in variable[device_name, :]
        if !(JuMP.value(var) <= (limit + 1e-2))
            return false
        end
    end
    return true
end

function check_flow_variable_values(
    model::DecisionModel,
    ::Type{T},
    ::Type{U},
    device_name::String,
    limit_min::Float64,
    limit_max::Float64,
) where {T <: PSI.VariableType, U <: PSY.Component}
    psi_cont = PSI.get_optimization_container(model)
    variable = PSI.get_variable(psi_cont, T(), U)
    for var in variable[device_name, :]
        if !(JuMP.value(var) <= (limit_max + 1e-2)) ||
           !(JuMP.value(var) >= (limit_min - 1e-2))
            return false
        end
    end
    return true
end

function check_flow_variable_values(
    model::DecisionModel,
    ::Type{T},
    ::Type{U},
    ::Type{V},
    device_name::String,
    limit_min::Float64,
    limit_max::Float64,
) where {T <: PSI.VariableType, U <: PSI.VariableType, V <: PSY.Component}
    psi_cont = PSI.get_optimization_container(model)
    time_steps = PSI.get_time_steps(psi_cont)
    pvariable = PSI.get_variable(psi_cont, T(), V)
    qvariable = PSI.get_variable(psi_cont, U(), V)
    for t in time_steps
        fp = JuMP.value(pvariable[device_name, t])
        fq = JuMP.value(qvariable[device_name, t])
        flow = sqrt((fp)^2 + (fq)^2)
        if !(flow <= (limit_max + 1e-2)^2) || !(flow >= (limit_min - 1e-2)^2)
            return false
        end
    end
    return true
end

function check_flow_variable_values(
    model::DecisionModel,
    ::Type{T},
    ::Type{U},
    ::Type{V},
    device_name::String,
    limit::Float64,
) where {T <: PSI.VariableType, U <: PSI.VariableType, V <: PSY.Component}
    psi_cont = PSI.get_optimization_container(model)
    time_steps = PSI.get_time_steps(psi_cont)
    pvariable = PSI.get_variable(psi_cont, T(), V)
    qvariable = PSI.get_variable(psi_cont, U(), V)
    for t in time_steps
        fp = JuMP.value(pvariable[device_name, t])
        fq = JuMP.value(qvariable[device_name, t])
        flow = sqrt((fp)^2 + (fq)^2)
        if !(flow <= (limit + 1e-2)^2)
            return false
        end
    end
    return true
end

function PSI.jump_value(int::Int)
    @warn("This is for testing purposes only.")
    return int
end

function _check_constraint_bounds(bounds::PSI.ConstraintBounds, valid_bounds::NamedTuple)
    @test bounds.coefficient.min == valid_bounds.coefficient.min
    @test bounds.coefficient.max == valid_bounds.coefficient.max
    @test bounds.rhs.min == valid_bounds.rhs.min
    @test bounds.rhs.max == valid_bounds.rhs.max
end

function _check_variable_bounds(bounds::PSI.VariableBounds, valid_bounds::NamedTuple)
    @test bounds.bounds.min == valid_bounds.min
    @test bounds.bounds.max == valid_bounds.max
end

function check_duration_on_initial_conditions_values(
    model,
    ::Type{T},
) where {T <: PSY.Component}
    initial_conditions_data =
        PSI.get_initial_conditions_data(PSI.get_optimization_container(model))
    duration_on_data = PSI.get_initial_condition(
        PSI.get_optimization_container(model),
        InitialTimeDurationOn(),
        T,
    )
    for ic in duration_on_data
        name = PSY.get_name(ic.component)
        on_var = PSI.get_initial_condition_value(initial_conditions_data, OnVariable(), T)[
            1,
            name,
        ]
        duration_on = JuMP.value(PSI.get_value(ic))
        if on_var == 1.0 && PSY.get_status(ic.component)
            @test duration_on == PSY.get_time_at_status(ic.component)
        elseif on_var == 1.0 && !PSY.get_status(ic.component)
            @test duration_on == 0.0
        end
    end
end

function check_duration_off_initial_conditions_values(
    model,
    ::Type{T},
) where {T <: PSY.Component}
    initial_conditions_data =
        PSI.get_initial_conditions_data(PSI.get_optimization_container(model))
    duration_off_data = PSI.get_initial_condition(
        PSI.get_optimization_container(model),
        InitialTimeDurationOff(),
        T,
    )
    for ic in duration_off_data
        name = PSY.get_name(ic.component)
        on_var = PSI.get_initial_condition_value(initial_conditions_data, OnVariable(), T)[
            1,
            name,
        ]
        duration_off = JuMP.value(PSI.get_value(ic))
        if on_var == 0.0 && !PSY.get_status(ic.component)
            @test duration_off == PSY.get_time_at_status(ic.component)
        elseif on_var == 0.0 && PSY.get_status(ic.component)
            @test duration_off == 0.0
        end
    end
end

function check_energy_initial_conditions_values(model, ::Type{T}) where {T <: PSY.Component}
    ic_data = PSI.get_initial_condition(
        PSI.get_optimization_container(model),
        InitialEnergyLevel(),
        T,
    )
    for ic in ic_data
        name = PSY.get_name(ic.component)
        e_value = JuMP.value(PSI.get_value(ic))
        @test PSY.get_initial_energy(ic.component) == e_value
    end
end

function check_energy_initial_conditions_values(model, ::Type{T}) where {T <: PSY.HydroGen}
    ic_data = PSI.get_initial_condition(
        PSI.get_optimization_container(model),
        InitialEnergyLevel(),
        T,
    )
    for ic in ic_data
        name = PSY.get_name(ic.component)
        e_value = JuMP.value(PSI.get_value(ic))
        @test PSY.get_initial_storage(ic.component) == e_value
    end
end

function check_status_initial_conditions_values(model, ::Type{T}) where {T <: PSY.Component}
    initial_conditions =
        PSI.get_initial_condition(PSI.get_optimization_container(model), DeviceStatus(), T)
    initial_conditions_data =
        PSI.get_initial_conditions_data(PSI.get_optimization_container(model))
    for ic in initial_conditions
        name = PSY.get_name(ic.component)
        status = PSI.get_initial_condition_value(initial_conditions_data, OnVariable(), T)[
            1,
            name,
        ]
        @test JuMP.value(PSI.get_value(ic)) == status
    end
end

function check_active_power_initial_condition_values(
    model,
    ::Type{T},
) where {T <: PSY.Component}
    initial_conditions =
        PSI.get_initial_condition(PSI.get_optimization_container(model), DevicePower(), T)
    initial_conditions_data =
        PSI.get_initial_conditions_data(PSI.get_optimization_container(model))
    for ic in initial_conditions
        name = PSY.get_name(ic.component)
        power = PSI.get_initial_condition_value(
            initial_conditions_data,
            ActivePowerVariable(),
            T,
        )[
            1,
            name,
        ]
        @test JuMP.value(PSI.get_value(ic)) == power
    end
end

function check_active_power_abovemin_initial_condition_values(
    model,
    ::Type{T},
) where {T <: PSY.Component}
    initial_conditions = PSI.get_initial_condition(
        PSI.get_optimization_container(model),
        PSI.DeviceAboveMinPower(),
        T,
    )
    initial_conditions_data =
        PSI.get_initial_conditions_data(PSI.get_optimization_container(model))
    for ic in initial_conditions
        name = PSY.get_name(ic.component)
        power = PSI.get_initial_condition_value(
            initial_conditions_data,
            PSI.PowerAboveMinimumVariable(),
            T,
        )[
            1,
            name,
        ]
        @test JuMP.value(PSI.get_value(ic)) == power
    end
end

function check_initialization_variable_count(
    model,
    ::S,
    ::Type{T},
) where {S <: PSI.VariableType, T <: PSY.Component}
    container = PSI.get_optimization_container(model)
    initial_conditions_data = PSI.get_initial_conditions_data(container)
    no_component = length(PSY.get_components(T, model.sys, x -> x.available))
    variable = PSI.get_initial_condition_value(initial_conditions_data, S(), T)
    rows, cols = size(variable)
    @test rows * cols == no_component * PSI.INITIALIZATION_PROBLEM_HORIZON
end

function check_variable_count(
    model,
    ::S,
    ::Type{T},
) where {S <: PSI.VariableType, T <: PSY.Component}
    no_component = length(PSY.get_components(T, model.sys, x -> x.available))
    time_steps = PSI.get_time_steps(PSI.get_optimization_container(model))[end]
    variable = PSI.get_variable(PSI.get_optimization_container(model), S(), T)
    @test length(variable) == no_component * time_steps
end

function check_initialization_constraint_count(
    model,
    ::S,
    ::Type{T};
    filter_func=x -> x.available,
    meta=PSI.CONTAINER_KEY_EMPTY_META,
) where {S <: PSI.ConstraintType, T <: PSY.Component}
    container = model.internal.ic_model_container
    no_component = length(PSY.get_components(T, model.sys, filter_func))
    time_steps = PSI.get_time_steps(container)[end]
    constraint = PSI.get_constraint(container, S(), T, meta)
    @test length(constraint) == no_component * time_steps
end

function check_constraint_count(
    model,
    ::S,
    ::Type{T};
    filter_func=x -> x.available,
    meta=PSI.CONTAINER_KEY_EMPTY_META,
) where {S <: PSI.ConstraintType, T <: PSY.Component}
    no_component = length(PSY.get_components(T, model.sys, filter_func))
    time_steps = PSI.get_time_steps(PSI.get_optimization_container(model))[end]
    constraint = PSI.get_constraint(PSI.get_optimization_container(model), S(), T, meta)
    @test length(constraint) == no_component * time_steps
end

function check_constraint_count(
    model,
    ::PSI.RampConstraint,
    ::Type{T},
) where {T <: PSY.Component}
    container = PSI.get_optimization_container(model)
    set_name =
        PSY.get_name.(
            PSI._get_ramp_constraint_devices(
                container,
                get_components(T, model.sys, x -> x.available),
            ),
        )
    check_constraint_count(
        model,
        PSI.RampConstraint(),
        T;
        meta="up",
        filter_func=x -> x.name in set_name,
    )
    check_constraint_count(
        model,
        PSI.RampConstraint(),
        T;
        meta="dn",
        filter_func=x -> x.name in set_name,
    )
    return
end

function check_constraint_count(
    model,
    ::PSI.DurationConstraint,
    ::Type{T},
) where {T <: PSY.Component}
    container = PSI.get_optimization_container(model)
    resolution = PSI.get_resolution(container)
    steps_per_hour = 60 / Dates.value(Dates.Minute(resolution))
    fraction_of_hour = 1 / steps_per_hour
    duration_devices = filter!(
        x -> !(
            PSY.get_time_limits(x).up <= fraction_of_hour &&
            PSY.get_time_limits(x).down <= fraction_of_hour
        ),
        collect(get_components(T, model.sys, x -> x.available)),
    )
    set_name = PSY.get_name.(duration_devices)
    check_constraint_count(
        model,
        PSI.DurationConstraint(),
        T;
        meta="up",
        filter_func=x -> x.name in set_name,
    )
    return check_constraint_count(
        model,
        PSI.DurationConstraint(),
        T;
        meta="dn",
        filter_func=x -> x.name in set_name,
    )
end

function check_constraint_count(
    model,
    ::PSI.MustRunConstraint,
    ::Type{T},
) where {T <: PSY.Component}
    container = PSI.get_optimization_container(model)
    _devices =
        filter!(x -> x.must_run, collect(get_components(T, model.sys, x -> x.available)))
    set_name = PSY.get_name.(_devices)
    return check_constraint_count(
        model,
        PSI.MustRunConstraint(),
        T,
        filter_func=x -> x.name in set_name,
    )
end


text:
using RobustMeans
using Test

@testset "RobustMean.jl" begin
using Distributions
using Statistics, Random
Random.seed!(3)
n = 3 * 7 * 8
M = 10^5
b = 2.0001
dist = Pareto(b)
μ = mean(dist)
σ = std(dist)
x = rand(dist, M, n)

δ = 3exp(-8)
p = 1 # for Minsker Ndaoud
estimator = [EmpiricalMean(), MedianOfMean(1), TrimmedMean(0, 0), Catoni(σ), Huber(σ), LeeValiant(), MinskerNdaoud(p)]

μ̂ = zeros(M, length(estimator))
for (i, es) in enumerate(estimator)
    for m in 1:M
        μ̂[m, i] = mean(x[m, :], δ, es)
    end
end

for (i, es) in enumerate(estimator)
    if typeof(es) <: MinskerNdaoud
        continue # constant in the bound not explciit.
    else
        @test quantile(μ̂[:, i], 1 - δ) < σ * bound(δ, n, es)
    end
end

end


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

module PermutedDimsArrays

import Base: permutedims, permutedims!
export PermutedDimsArray

# Some day we will want storage-order-aware iteration, so put perm in the parameters
struct PermutedDimsArray{T,N,perm,iperm,AA<:AbstractArray} <: AbstractArray{T,N}
    parent::AA

    function PermutedDimsArray{T,N,perm,iperm,AA}(data::AA) where {T,N,perm,iperm,AA<:AbstractArray}
        (isa(perm, NTuple{N,Int}) && isa(iperm, NTuple{N,Int})) || error("perm and iperm must both be NTuple{$N,Int}")
        isperm(perm) || throw(ArgumentError(string(perm, " is not a valid permutation of dimensions 1:", N)))
        all(map(d->iperm[perm[d]]==d, 1:N)) || throw(ArgumentError(string(perm, " and ", iperm, " must be inverses")))
        new(data)
    end
end

"""
    PermutedDimsArray(A, perm) -> B

Given an AbstractArray `A`, create a view `B` such that the
dimensions appear to be permuted. Similar to `permutedims`, except
that no copying occurs (`B` shares storage with `A`).

See also: [`permutedims`](@ref).

# Examples
```jldoctest
julia> A = rand(3,5,4);

julia> B = PermutedDimsArray(A, (3,1,2));

julia> size(B)
(4, 3, 5)

julia> B[3,1,2] == A[1,2,3]
true
```
"""
function PermutedDimsArray(data::AbstractArray{T,N}, perm) where {T,N}
    length(perm) == N || throw(ArgumentError(string(perm, " is not a valid permutation of dimensions 1:", N)))
    iperm = invperm(perm)
    PermutedDimsArray{T,N,(perm...,),(iperm...,),typeof(data)}(data)
end

Base.parent(A::PermutedDimsArray) = A.parent
Base.size(A::PermutedDimsArray{T,N,perm}) where {T,N,perm} = genperm(size(parent(A)), perm)
Base.axes(A::PermutedDimsArray{T,N,perm}) where {T,N,perm} = genperm(axes(parent(A)), perm)

Base.similar(A::PermutedDimsArray, T::Type, dims::Base.Dims) = similar(parent(A), T, dims)

Base.unsafe_convert(::Type{Ptr{T}}, A::PermutedDimsArray{T}) where {T} = Base.unsafe_convert(Ptr{T}, parent(A))

# It's OK to return a pointer to the first element, and indeed quite
# useful for wrapping C routines that require a different storage
# order than used by Julia. But for an array with unconventional
# storage order, a linear offset is ambiguous---is it a memory offset
# or a linear index?
Base.pointer(A::PermutedDimsArray, i::Integer) = throw(ArgumentError("pointer(A, i) is deliberately unsupported for PermutedDimsArray"))

function Base.strides(A::PermutedDimsArray{T,N,perm}) where {T,N,perm}
    s = strides(parent(A))
    ntuple(d->s[perm[d]], Val(N))
end

@inline function Base.getindex(A::PermutedDimsArray{T,N,perm,iperm}, I::Vararg{Int,N}) where {T,N,perm,iperm}
    @boundscheck checkbounds(A, I...)
    @inbounds val = getindex(A.parent, genperm(I, iperm)...)
    val
end
@inline function Base.setindex!(A::PermutedDimsArray{T,N,perm,iperm}, val, I::Vararg{Int,N}) where {T,N,perm,iperm}
    @boundscheck checkbounds(A, I...)
    @inbounds setindex!(A.parent, val, genperm(I, iperm)...)
    val
end

@inline genperm(I::NTuple{N,Any}, perm::Dims{N}) where {N} = ntuple(d -> I[perm[d]], Val(N))
@inline genperm(I, perm::AbstractVector{Int}) = genperm(I, (perm...,))

"""
    permutedims(A::AbstractArray, perm)

Permute the dimensions of array `A`. `perm` is a vector specifying a permutation of length
`ndims(A)`.

See also: [`PermutedDimsArray`](@ref).

# Examples
```jldoctest
julia> A = reshape(Vector(1:8), (2,2,2))
2×2×2 Array{Int64,3}:
[:, :, 1] =
 1  3
 2  4

[:, :, 2] =
 5  7
 6  8

julia> permutedims(A, [3, 2, 1])
2×2×2 Array{Int64,3}:
[:, :, 1] =
 1  3
 5  7

[:, :, 2] =
 2  4
 6  8
```
"""
function permutedims(A::AbstractArray, perm)
    dest = similar(A, genperm(axes(A), perm))
    permutedims!(dest, A, perm)
end

"""
    permutedims(m::AbstractMatrix)

Permute the dimensions of the matrix `m`, by flipping the elements across the diagonal of
the matrix. Differs from `LinearAlgebra`'s [`transpose`](@ref) in that the
operation is not recursive.

# Examples
```jldoctest; setup = :(using LinearAlgebra)
julia> a = [1 2; 3 4];

julia> b = [5 6; 7 8];

julia> c = [9 10; 11 12];

julia> d = [13 14; 15 16];

julia> X = [[a] [b]; [c] [d]]
2×2 Array{Array{Int64,2},2}:
 [1 2; 3 4]     [5 6; 7 8]
 [9 10; 11 12]  [13 14; 15 16]

julia> permutedims(X)
2×2 Array{Array{Int64,2},2}:
 [1 2; 3 4]  [9 10; 11 12]
 [5 6; 7 8]  [13 14; 15 16]

julia> transpose(X)
2×2 Transpose{Transpose{Int64,Array{Int64,2}},Array{Array{Int64,2},2}}:
 [1 3; 2 4]  [9 11; 10 12]
 [5 7; 6 8]  [13 15; 14 16]
```
"""
permutedims(A::AbstractMatrix) = permutedims(A, (2,1))

"""
    permutedims(v::AbstractVector)

Reshape vector `v` into a `1 × length(v)` row matrix.
Differs from `LinearAlgebra`'s [`transpose`](@ref) in that
the operation is not recursive.

# Examples
```jldoctest; setup = :(using LinearAlgebra)
julia> permutedims([1, 2, 3, 4])
1×4 Array{Int64,2}:
 1  2  3  4

julia> V = [[[1 2; 3 4]]; [[5 6; 7 8]]]
2-element Array{Array{Int64,2},1}:
 [1 2; 3 4]
 [5 6; 7 8]

julia> permutedims(V)
1×2 Array{Array{Int64,2},2}:
 [1 2; 3 4]  [5 6; 7 8]

julia> transpose(V)
1×2 Transpose{Transpose{Int64,Array{Int64,2}},Array{Array{Int64,2},1}}:
 [1 3; 2 4]  [5 7; 6 8]
```
"""
permutedims(v::AbstractVector) = reshape(v, (1, length(v)))

"""
    permutedims!(dest, src, perm)

Permute the dimensions of array `src` and store the result in the array `dest`. `perm` is a
vector specifying a permutation of length `ndims(src)`. The preallocated array `dest` should
have `size(dest) == size(src)[perm]` and is completely overwritten. No in-place permutation
is supported and unexpected results will happen if `src` and `dest` have overlapping memory
regions.

See also [`permutedims`](@ref).
"""
function permutedims!(dest, src::AbstractArray, perm)
    Base.checkdims_perm(dest, src, perm)
    P = PermutedDimsArray(dest, invperm(perm))
    _copy!(P, src)
    return dest
end

function Base.copyto!(dest::PermutedDimsArray{T,N}, src::AbstractArray{T,N}) where {T,N}
    checkbounds(dest, axes(src)...)
    _copy!(dest, src)
end
Base.copyto!(dest::PermutedDimsArray, src::AbstractArray) = _copy!(dest, src)

function _copy!(P::PermutedDimsArray{T,N,perm}, src) where {T,N,perm}
    # If dest/src are "close to dense," then it pays to be cache-friendly.
    # Determine the first permuted dimension
    d = 0  # d+1 will hold the first permuted dimension of src
    while d < ndims(src) && perm[d+1] == d+1
        d += 1
    end
    if d == ndims(src)
        copyto!(parent(P), src) # it's not permuted
    else
        R1 = CartesianIndices(axes(src)[1:d])
        d1 = findfirst(isequal(d+1), perm)::Int  # first permuted dim of dest
        R2 = CartesianIndices(axes(src)[d+2:d1-1])
        R3 = CartesianIndices(axes(src)[d1+1:end])
        _permutedims!(P, src, R1, R2, R3, d+1, d1)
    end
    return P
end

@noinline function _permutedims!(P::PermutedDimsArray, src, R1::CartesianIndices{0}, R2, R3, ds, dp)
    ip, is = axes(src, dp), axes(src, ds)
    for jo in first(ip):8:last(ip), io in first(is):8:last(is)
        for I3 in R3, I2 in R2
            for j in jo:min(jo+7, last(ip))
                for i in io:min(io+7, last(is))
                    @inbounds P[i, I2, j, I3] = src[i, I2, j, I3]
                end
            end
        end
    end
    P
end

@noinline function _permutedims!(P::PermutedDimsArray, src, R1, R2, R3, ds, dp)
    ip, is = axes(src, dp), axes(src, ds)
    for jo in first(ip):8:last(ip), io in first(is):8:last(is)
        for I3 in R3, I2 in R2
            for j in jo:min(jo+7, last(ip))
                for i in io:min(io+7, last(is))
                    for I1 in R1
                        @inbounds P[I1, i, I2, j, I3] = src[I1, i, I2, j, I3]
                    end
                end
            end
        end
    end
    P
end

function Base.showarg(io::IO, A::PermutedDimsArray{T,N,perm}, toplevel) where {T,N,perm}
    print(io, "PermutedDimsArray(")
    Base.showarg(io, parent(A), false)
    print(io, ", ", perm, ')')
    toplevel && print(io, " with eltype ", eltype(A))
end

end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
using Documenter
using BPFnative
using Dates

makedocs(
    modules = [BPFnative],
    sitename= "BPFnative.jl",
    authors = "Julian Samaroo",
    format = Documenter.HTML(),
)

deploydocs(
    repo = "github.com/lawless-m/BPFnative.jl.git",
 #   devbranch = "master",
 #   push_preview = true,
)


text:
# Unit tests

push!(LOAD_PATH, "../src")
ENV["PDEPUMIINTERFACE_TESTING"] = true

using Base.Test
using PumiInterface
using ODLCommonTools
import ODLCommonTools.sview
using SummationByParts
using PdePumiInterface
using PumiConfig

import Base.isapprox

# generic fallback
function isapprox(a, b)
  return a == b
end

function not_isapprox(args...; kwargs...)
  return !isapprox(args...; kwargs...)
end

@testset "Testing PUMIInterface.jl" begin

  # do tests here
  @testset "Testing apf.initilization" begin
  # create mesh
  dmg_name = ".null"
  #smb_name = "../../mesh_files/quarter_vortex3l.smb"
  # smb_name = "../../mesh_files/quarter_vortex8l.smb"
  smb_name = "./tri8l.smb"

  order = 1  # linear elements
  m_ptr, dim = apf.loadMesh(dmg_name, smb_name, order)
  mshape_ptr, num_Entities, n_arr = apf.initMesh(m_ptr)
#  num_Entities, m_ptr, mshape_ptr, dim, n_arr = apf.init2(dmg_name, smb_name, order)
  println("length(n_arr) = ", length(n_arr))

  @test ( num_Entities[1] )== 9
  @test ( num_Entities[2] )== 16
  @test ( num_Entities[3] )== 8
  @test ( m_ptr )!=C_NULL
  @test ( mshape_ptr )!=C_NULL

  # check the apf.countJ() function
  @test ( apf.countJ(m_ptr, 0) )== 9
  @test ( apf.countJ(m_ptr, 1) )== 16
  @test ( apf.countJ(m_ptr, 2) )== 8
  @test ( apf.countJ(m_ptr, 3) )== 0

  # check numberings
  vertN_ptr = n_arr[1] #getVertNumbering()
  edgeN_ptr = n_arr[2] #getEdgeNumbering()
  faceN_ptr = n_arr[3] #getFaceNumbering()

  numberings = [vertN_ptr, edgeN_ptr, faceN_ptr]

#  apf.printEdgeVertNumbers(edgeN_ptr, vertN_ptr)

#  resetVertIt()
#  resetEdgeIt()
#  resetFaceIt()

  vertit = apf.MeshIterator(m_ptr, 0)
  edgeit = apf.MeshIterator(m_ptr, 1)
  faceit = apf.MeshIterator(m_ptr, 2)

  vert = apf.deref(vertit)
  edge = apf.deref(edgeit)
  face = apf.deref(faceit)
#  vert = getVert()
#  edge = getEdge()
#  face = getFace()

  vertnum = apf.getNumberJ(vertN_ptr, vert, 0, 0)

  edgenum = apf.getNumberJ(edgeN_ptr, edge, 0, 0)

  facenum = apf.getNumberJ(faceN_ptr, face, 0, 0)

  @test ( apf.getNumberJ(vertN_ptr, vert, 0, 0) )== 0
  @test ( apf.getNumberJ(edgeN_ptr, edge, 0, 0) )== 0
  @test ( apf.getNumberJ(faceN_ptr, face, 0, 0) )== 0

  for i=1:num_Entities[1]  # loop over verts
    entity = apf.iterate(vertit)
    @test ( apf.getNumberJ(vertN_ptr, entity, 0, 0) )== i-1
#    incrementVertIt()
  end
  apf.free(vertit)
#  resetVertIt()

  for i=1:num_Entities[2]  # loop over edges
    entity = apf.iterate(edgeit)
    @test ( apf.getNumberJ(edgeN_ptr, entity, 0, 0) )== i-1
#    incrementEdgeIt()
  end
  apf.free(edgeit)
#  resetEdgeIt()

  for i=1:num_Entities[3]  # loop over faces
    entity = apf.iterate(faceit)
    @test ( apf.getNumberJ(faceN_ptr, entity, 0, 0) )== i-1
#    incrementFaceIt()
  end
  apf.free(faceit)
#  resetFaceIt()

  @test ( apf.getMeshDimension(m_ptr) )== 2
  @test ( apf.getType(m_ptr, vert) )== 0
  @test ( apf.getType(m_ptr, edge) )== 1
  @test ( apf.getType(m_ptr, face) )== 2

  apf.printEdgeVertNumbers(edgeN_ptr, vertN_ptr)
  apf.printFaceVertNumbers( faceN_ptr, vertN_ptr)
  # need to recreate these test post apf.reordering

  # check downward adjacency function
  down_entities, num_down = apf.getDownward(m_ptr, face, 1) # face -> edges
  down_nums = zeros(Int, num_down)
  for i=1:num_down
    down_nums[i] = apf.getNumberJ(edgeN_ptr, down_entities[i], 0, 0)
  end
   @test ( down_nums )== [1, 5, 0]

  down_entities, num_down = apf.getDownward(m_ptr, face, 0) # face -> verts
  down_nums = zeros(Int, num_down)
  for i=1:num_down
    down_nums[i] = apf.getNumberJ(vertN_ptr, down_entities[i], 0, 0)
  end
   @test ( down_nums )== [0, 2, 1]

  down_entities, num_down = apf.getDownward(m_ptr, edge, 0) # edge -> verts
  down_nums = zeros(Int, num_down)
  for i=1:num_down
    down_nums[i] = apf.getNumberJ(vertN_ptr, down_entities[i], 0, 0)
  end
   @test ( down_nums )== [1, 0]

  # check upward adjacency function
  println("testing vert -> edge")
  num_up = apf.countAdjacent(m_ptr, vert, 1)  # vert -> edge
  up_entities = apf.getAdjacent( num_up)
  up_nums = zeros(Int, num_up)
  for i=1:num_up
    up_nums[i] = apf.getNumberJ(edgeN_ptr, up_entities[i], 0, 0)
  end

#  println("up_nums = ", up_nums)
  @test ( up_nums )== [2, 1, 0]

 num_up = apf.countAdjacent(m_ptr, vert, 2)  # vert -> face
  up_entities = apf.getAdjacent( num_up)
  up_nums = zeros(Int, num_up)
  for i=1:num_up
    up_nums[i] = apf.getNumberJ(faceN_ptr, up_entities[i], 0, 0)
  end

#  println("up_nums = ", up_nums)
  @test ( up_nums )== [1, 0]

  num_up = apf.countAdjacent(m_ptr, edge, 2)  # edge -> face
  up_entities = apf.getAdjacent( num_up)
  up_nums = zeros(Int, num_up)
  for i=1:num_up
    up_nums[i] = apf.getNumberJ(faceN_ptr, up_entities[i], 0, 0)
  end

#  println("up_nums = ", up_nums)
  @test ( up_nums )== [0]

  # check second order adjacencies

  # check FieldShape and EntityShape functions
  @test ( apf.hasNodesIn(mshape_ptr, 0) )== true
  @test ( apf.hasNodesIn(mshape_ptr, 1) )== false
  @test ( apf.hasNodesIn(mshape_ptr, 2) )== false

  @test ( apf.countNodesOn(mshape_ptr, 0) )== 1
  @test ( apf.countNodesOn(mshape_ptr, 1) )== 0
  @test ( apf.countNodesOn(mshape_ptr, 2) )== 0

  # check shape functions, jacobians, integration points here?

  # check second order adjacencies

  # check vertices connected through edges
  in_dim = 1
  out_dim = 0
  cnt = apf.countBridgeAdjacent(m_ptr, vert, in_dim, out_dim)
  adj = Array{Ptr{Void}}(cnt)
  adj_num = Array{Int}(cnt)
  apf.getBridgeAdjacent(adj)
  for i=1:cnt
    adj_num[i] = apf.getNumberJ(numberings[out_dim + 1], adj[i], 0, 0)
  end
  sort!(adj_num)
  @test ( adj_num )== [1,2,3]

  # target_dimension = bridge dimension not defined
  in_dim = 1
  out_dim = 2
  cnt = apf.countBridgeAdjacent(m_ptr, vert, in_dim, out_dim)
  adj = Array{Ptr{Void}}(cnt)
  adj_num = Array{Int}(cnt)
  apf.getBridgeAdjacent(adj)
  for i=1:cnt
    adj_num[i] = apf.getNumberJ(numberings[out_dim + 1], adj[i], 0, 0)
  end
  sort!(adj_num)
  @test ( adj_num )== [0,1]

  # edge connected to vertices by edges
  in_dim = 0
  out_dim = 1
  cnt = apf.countBridgeAdjacent(m_ptr, edge, in_dim, out_dim)
  adj = Array{Ptr{Void}}(cnt)
  adj_num = Array{Int}(cnt)
  apf.getBridgeAdjacent(adj)
  for i=1:cnt
    adj_num[i] = apf.getNumberJ(numberings[out_dim + 1], adj[i], 0, 0)
  end
  sort!(adj_num)
  @test ( adj_num )== sort!([1,2,3,4,5])

  # edge connected to edges by faces
  in_dim = 2
  out_dim = 1
  cnt = apf.countBridgeAdjacent(m_ptr, edge, in_dim, out_dim)
  adj = Array{Ptr{Void}}(cnt)
  adj_num = Array{Int}(cnt)
  apf.getBridgeAdjacent(adj)
  for i=1:cnt
    adj_num[i] = apf.getNumberJ(numberings[out_dim + 1], adj[i], 0, 0)
  end
  sort!(adj_num)
  @test ( adj_num )== [1, 5]

  # face connected to edges by vertex
  in_dim = 0
  out_dim = 1
  cnt = apf.countBridgeAdjacent(m_ptr, face, in_dim, out_dim)
  adj = Array{Ptr{Void}}(cnt)
  adj_num = Array{Int}(cnt)
  apf.getBridgeAdjacent(adj)
  for i=1:cnt
    adj_num[i] = apf.getNumberJ(numberings[out_dim + 1], adj[i], 0, 0)
  end
  sort!(adj_num)
  ans = [0,1,2,3,4,5,6,7,8,9]

  @test ( adj_num )== ans

  # face connected to faces by vert
  in_dim = 0
  out_dim = 2
  cnt = apf.countBridgeAdjacent(m_ptr, face, in_dim, out_dim)
  adj = Array{Ptr{Void}}(cnt)
  adj_num = Array{Int}(cnt)
  apf.getBridgeAdjacent(adj)
  for i=1:cnt
    adj_num[i] = apf.getNumberJ(numberings[out_dim + 1], adj[i], 0, 0)
  end
  sort!(adj_num)
  @test ( adj_num )== sort!([1, 2,3,4,5,6])

  # face connected to faces by edge
  in_dim = 1
  out_dim = 2
  cnt = apf.countBridgeAdjacent(m_ptr, face, in_dim, out_dim)
  adj = Array{Ptr{Void}}(cnt)
  adj_num = Array{Int}(cnt)
  apf.getBridgeAdjacent(adj)
  for i=1:cnt
    adj_num[i] = apf.getNumberJ(numberings[out_dim + 1], adj[i], 0, 0)
  end
  sort!(adj_num)
  @test ( adj_num )== [1, 3]

  mshape2 = apf.getFieldShape(1, 4, 2)
  node_entities = apf.getNodeEntities(m_ptr, mshape2, face)
  verts = node_entities[1:3]
  edges = node_entities[4:12]
  faces = node_entities[13:18]
  for i in verts
    @test ( apf.getType(m_ptr, i) )== apf.VERTEX
  end
  for i in edges
    @test ( apf.getType(m_ptr, i) )== apf.EDGE
  end
  for i in faces
    @test ( apf.getType(m_ptr, i) )== apf.TRIANGLE
  end

  @test ( apf.countPeers(m_ptr, apf.VERTEX) )== 0
  @test ( apf.countPeers(m_ptr, apf.EDGE) )== 0

  for i in verts
    @test ( apf.countRemotes(m_ptr, i) )== 0
  end
  for i in edges
    @test ( apf.countRemotes(m_ptr, i) )== 0
  end

  # test MeshEntity iteration
  nvert = 0
  for e in apf.MeshIterator(m_ptr, 0)
    nvert += 1
  end
  @test nvert == num_Entities[1]

  # test FieldEntity iteration
  # count number of vertices in 1st order lagrange field
  fshape = apf.getFieldShape(0, 1, apf.getDimension(m_ptr))
  it = apf.FieldEntityIt(m_ptr, fshape)
  @test Base.iteratorsize(it) == Base.HasLength()
  @test Base.eltype(it) == Tuple{Ptr{Void}, Int}
  state = Base.start(it)
  old_state, new_state = Base.next(it, state)
  @test typeof(old_state) == Base.eltype(it)
  apf.free(it)

  nvert = 0
  for (e, edim) in apf.FieldEntityIt(m_ptr, fshape)
    @test edim == 0
    nvert += 1
  end
  @test nvert == num_Entities[1]
  println("nvert = ", nvert)

  # test enumerate
  i = 1
  for (idx, e) in enumerate(apf.FieldEntityIt(m_ptr, fshape))
    @test idx == i
    i += 1
  end

  # check coordinates
#=
  coords_vert = zeros(3,1)
  apf.getVertCoords(vert, coords_vert, 3, 1)
  coords = coords_vert[:,1]
  @test isapprox( coords, [-1.0, -1.0, 0])
  apf.getVertCoords(coords_vert, 3, 1)
  coords = coords_vert[:,1]
  @test isapprox( coords, [-1.0, -1.0, 0])

  coords_edge = zeros(3,2)
  apf.getEdgeCoords(edge, coords_edge, 3, 2)
  @test isapprox( coords_edge, [-1.0 0; -1 -1; 0 0])
  apf.getEdgeCoords(coords_edge, 3, 2)
  @test isapprox( coords_edge, [-1.0 0; -1 -1; 0 0])

  coords_el = zeros(3,3)
  apf.getFaceCoords(face, coords_el, 3, 3)
  @test isapprox( coords_el, [-1.0 0 0; -1 -1 0; 0 0 0])
  apf.getFaceCoords(coords_el, 3, 3)
 @test isapprox( coords_el, [-1.0 0 0; -1 -1 0; 0 0 0])

=#

 # check creating numberings
 n_ptr = apf.createNumberingJ(m_ptr, "testnumbering", mshape_ptr, 1)

 apf.numberJ(n_ptr, vert, 0, 0, 1)
 @test ( apf.getNumberJ(n_ptr, vert, 0, 0) )== 1
 apf.destroyNumbering(n_ptr)

  # test writeVtk all fields option
  # the mesh has a first order lagrange field, so make a 2nd order numbering
  nname = "testOmegaNumbering"
  nshape = apf.getFieldShape(2, 1, 2)
  n_ptr = apf.createNumberingJ(m_ptr, nname, nshape, 1)

  for (i, el_i) in enumerate(apf.MeshIterator(m_ptr, 2))
    for j=0:2
      apf.numberJ(n_ptr, el_i, j, 0, i)
    end
  end

  apf.writeVtkFiles("test_vtk", m_ptr)
  apf.writeVtkFiles("test_vtk_all", m_ptr, writeall=true)

  found_name = false
  for line in eachline("test_vtk/test_vtk.pvtu")
    if contains(line, nname)
      found_name = true
      break
    end
  end

  found_name_all = false
  for line in eachline("test_vtk_all/test_vtk_all.pvtu")
    if contains(line, "faceNums")
      found_name_all = true
      break
    end
  end

  @test !found_name
  @test found_name_all
  apf.destroyNumbering(n_ptr)

 # test mesh warping functions

 # find vertex at -1, -1
# resetVertIt()
 it = apf.MeshIterator(m_ptr, 0)
 coords = zeros(3)
 entity = C_NULL
 for i=1:num_Entities[1]
   vert_i = apf.iterate(it)
   apf.getPoint(m_ptr, vert_i, 0, coords)
   if coords[1] < -0.5 && coords[2] < -0.5
     entity = vert_i
     break
   end
 end
  apf.free(it)
 @test ( entity )!=C_NULL

 coords[1] *= 2
 coords[2] *= 2
 apf.setPoint(m_ptr, entity, 0, coords)
 apf.acceptChanges(m_ptr)
 apf.Verify(m_ptr)

 coords2 = copy(coords)
 apf.getPoint(m_ptr, entity, 0, coords2)
 for i=1:2
   @test  abs(coords[i] - coords2[i])  < 1e-12
 end

 apf.writeVtkFiles("warp_test", m_ptr)

 # test periodic mesh
  @test ( apf.hasMatching(m_ptr) )== false

  smb_name = "tri3_px.smb"
  m_ptr, dim = apf.loadMesh(dmg_name, smb_name, order)
  mshape_ptr, num_Entities, n_arr = apf.initMesh(m_ptr)
#  num_Entities, m_ptr, mshape_ptr = apf.init2(dmg_name, smb_name, order)
  shr_ptr = apf.getSharing(m_ptr)
  numVert = num_Entities[1]
  numEdge = num_Entities[2]
  numFace = num_Entities[3]
  # verify all meshentities are owned
#  resetAllIts2(m_ptr)
  nowned = 0

  @test ( apf.hasMatching(m_ptr) )== true

  it = apf.MeshIterator(m_ptr, 0)
  for i =1:numVert
    entity = apf.iterate(it)
    if apf.isOwned(shr_ptr, entity)
      nowned += 1
    end

#    incrementVertIt()
  end
  apf.free(it)
  @test ( nowned )== numVert - 4

  nowned = 0
  it = apf.MeshIterator(m_ptr, 1)
  for i=1:numEdge
    entity = apf.iterate(it)
    if apf.isOwned(shr_ptr, entity)
      nowned += 1
    end
#    incrementEdgeIt()
  end
  apf.free(it)
  @test ( nowned )== numEdge - 3

  ncopies = zeros(Int, 2)
  nmatches = zeros(Int, 2)
  part_nums = Array{Cint}(1)
  copies = Array{Ptr{Void}}(1)
  matches = Array{Ptr{Void}}(1)
#  resetAllIts2(m_ptr)
  it = apf.MeshIterator(m_ptr, 1)
  for i=1:numEdge
    entity = apf.iterate(it)
    n = apf.countCopies(shr_ptr, entity)
    n2 = apf.countMatches(m_ptr, entity)
    ncopies[n+1] += 1  # either 0 or 1 copy
    nmatches[n+1] += 1 # either 0 or 1 match
    if n > 0
      apf.getCopies(part_nums, copies)
      @test ( part_nums[1] )== 0
      e_type = apf.getType(m_ptr, copies[1])
      @test ( e_type )== apf.EDGE
    end
    if n2 > 0
      apf.getMatches(part_nums, matches)
      @test ( part_nums[1] )== 0
      e_type = apf.getType(m_ptr, copies[1])
      @test ( e_type )== apf.EDGE
    end
#    incrementEdgeIt()
  end
  apf.free(it)

  @test ( ncopies[2] )== 6  # 3 x 3 element mesh has 6 shared edges
  @test ( nmatches[2] )== 6

  # test that convertSimMesh works correctly
  if HAVE_SIMMETRIX
    startdir = pwd()
    dirname = "test_convertSim"
    srcdir = joinpath(pwd(), "meshes", "UnitCubeCurve")
    if isdir(dirname)
      rm(dirname, recursive=true)
    end
    mkdir(dirname)
    cd(dirname)
    fname = "UnitCubeCurve"
    cp(joinpath(srcdir, "$fname.smd"), "./$fname.smd")
    cp(joinpath(srcdir, "$fname.sms"), "./$fname.sms")
    exepath = joinpath(pwd(), "..", "..", "install", "bin")
    run(`$exepath/convertSimMesh $fname.smd $fname.sms $fname.smb`)

    order = 2  # linear elements
    m_ptr, dim = apf.loadMesh("$fname.smd", "$fname.smb", order)
    # check that the quadratic parameter field is present
    f_ptr = apf.findField(m_ptr, "edge_params")
    @test f_ptr != C_NULL
    cd(startdir)
  end

end  # end context
end

#=
@testset "Testing PdePumiInterface3.jl" begin

  # do tests here
  @testset "Testing apf.initilization" begin
  dmg_name =  ".null"
  smb_name = "cube8l.smb"
  order = 1
  sbp = TetSBP{Float64}(degree=order)

  mesh = PumiMesh3{Float64}(dmg_name, smb_name, order, sbp, dofpernode=4)

  @test ( PdePumiInterface.constructNodemap(mesh, 0) )== [1, 3, 2]
  @test ( PdePumiInterface.constructNodemap(mesh, 1) )== [3, 2, 1]
  @test ( PdePumiInterface.constructNodemap(mesh, 2) )== [2, 1, 3]

  order = 2
  mesh = PumiMesh3{Float64}(dmg_name, smb_name, order, sbp, dofpernode=4)

  @test ( PdePumiInterface.constructNodemap(mesh, 0) )== [1, 3, 2, 6, 5, 4]
  @test ( PdePumiInterface.constructNodemap(mesh, 1) )== [3, 2, 1, 5, 4, 6]
  @test ( PdePumiInterface.constructNodemap(mesh, 2) )== [2, 1, 3, 4, 6, 5]

  order = 3
  mesh = PumiMesh3{Float64}(dmg_name, smb_name, order, sbp, dofpernode=4)

  @test ( PdePumiInterface.constructNodemap(mesh, 0) )== [1, 3, 2, 9, 8, 7, 6, 5, 4, 10]
  @test ( PdePumiInterface.constructNodemap(mesh, 1) )== [3, 2, 1, 7, 6, 5, 4, 9, 8, 10]
  @test ( PdePumiInterface.constructNodemap(mesh, 2) )== [2, 1, 3, 5, 4, 9, 8, 7, 6, 10]

  order = 4
  mesh = PumiMesh3{Float64}(dmg_name, smb_name, order, sbp, dofpernode=4)

  @test ( PdePumiInterface.constructNodemap(mesh, 0) )== [1, 3, 2, 12, 11, 10, 9, 8, 7, 6, 5, 4, 13, 15, 14]
  @test ( PdePumiInterface.constructNodemap(mesh, 1) )== [3, 2, 1, 9, 8, 7, 6, 5, 4, 12, 11, 10, 15, 14, 13]
  @test ( PdePumiInterface.constructNodemap(mesh, 2) )== [2, 1, 3, 6, 5, 4, 12, 11, 10, 9, 8, 7, 14, 13, 15]

  end

end
=#

include("test_gmi.jl")
test_gmi()

include("test_math.jl")
include("test_funcs.jl")  # include functions used by both 2d and 3D
include("common_functions.jl")
include("test_adapt.jl") # mesh adaptation
include("test_parallel_types.jl")
include("pdepumiinterface.jl")
include("pdepumiinterface3.jl")

MPI.Barrier(MPI.COMM_WORLD)
if MPI.Initialized()
  MPI.Finalize()
end


text:
# Autogenerated wrapper script for NASM_jll for x86_64-linux-musl
export nasm, ndisasm

JLLWrappers.@generate_wrapper_header("NASM")
JLLWrappers.@declare_executable_product(nasm)
JLLWrappers.@declare_executable_product(ndisasm)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        nasm,
        "bin/nasm",
    )

    JLLWrappers.@init_executable_product(
        ndisasm,
        "bin/ndisasm",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
## Log-likelihood
# Unique session
function loglikelihood(mdl::Tm, latent::NT, data::StructVector) where {Tm <: AbstractMCMC.AbstractModel, NT <: NamedTuple}
    θ = deepcopy(latent)
    P = [cycle!(θ, mdl, obs) for obs in data]
    L = logpdf.(P, data.a)
    return L
end

# Multiple sessions
function loglikelihood(mdl::Tm, latent::NT, data::Array{Ts}) where {Tm <: AbstractMCMC.AbstractModel, NT <: NamedTuple, Ts <: StructVector}
    nsess = length(data)
    L = Float64[]
    for sess = 1:nsess
        θ = deepcopy(latent) # Re-initialize the latent variables
        P = [cycle!(θ, mdl, obs) for obs in data[sess]]
        L = append!(L, logpdf.(P, data[sess].a))
    end
    return L
end

## DIC
# Number of effective parameters
neff_params(::Val{:pD}, D_samples, D_avg) = mean(D_samples) - D_avg
neff_params(::Val{:pV}, D_samples, _) = 0.5*var(D_samples)

function dic(P::Posterior, pDIC = :pD)
    latent_avg = average(P.latent)
    lp_avg = sum(P.loglikelihood_func(latent_avg))
    lp_samples = sum.(P.loglikelihood_func.(P.latent))
    D_samples = -2 .* lp_samples
    D_avg = -2 * lp_avg
    pD = neff_params(Val(pDIC), D_samples, D_avg)
    return D_avg + 2 * pD, pD
end

## WAIC
function waic(P::Posterior)
    lp_y = hcat([P.loglikelihood_func(P.latent[i]) for i = 1:P.nsamples]...)
    N, K = size(lp_y)
    lppd = sum(logsumexp(lp_y, dims=2)) - N*log(K)
    pW = sum(var(lp_y, dims=2))
    return -2 * lppd + 2 * pW, pW
end

## BIC
function bic(P::Posterior)
    latent_avg = average(P.latent)
    lp_avg = sum(P.loglikelihood_func(latent_avg))
    return - 2 * lp_avg + P.nparams * log(P.ndata)
end

max_stars_repo_name: sqwayer/AnimalBehavior

max_issues_repo_name: sqwayer/AnimalBehavior


text:
include("../imports.jl")

@testset "`SanityCheck`" begin

    @testset "Checks pass" begin
        cb = SanityCheck()
        learner = testlearner(cb)
        @test_nowarn fit!(learner, 1)
    end

    @testset "No training iterator" begin
        cb = SanityCheck([CHECKS[1]], usedefault = false)
        learner = testlearner(cb)
        @test_nowarn fit!(learner, 1)
        cb = SanityCheck([CHECKS[3]], usedefault = false)
        learner = testlearner(cb)
        learner.data.training = nothing
        @test_throws Exception (@suppress fit!(learner, 1))
    end

    @testset "Iterators over tuples" begin
        cb = SanityCheck([CHECKS[3]], usedefault = false)
        learner = testlearner(cb)
        @test_nowarn fit!(learner, 1)

        cb = SanityCheck([CHECKS[3]], usedefault = false)
        learner = testlearner(cb)
        learner.data.training = 1:16
        @test_throws SanityCheckException (@suppress fit!(learner, 1))
    end
@testset "Optimization step" begin
        cb = SanityCheck([CHECKS[4]], usedefault = false)
        learner = testlearner(cb)
        @test_nowarn fit!(learner, 1)

        cb = SanityCheck([CHECKS[4]], usedefault = false)
        learner = testlearner(cb)
        learner.model = x -> "hi"
        @test_throws SanityCheckException (@suppress fit!(learner, 1))

        cb = SanityCheck([CHECKS[4]], usedefault = false)
        learner = testlearner(cb)
        learner.model = x -> "hi"
        @test_throws SanityCheckException (@suppress fit!(learner, 1))

        cb = SanityCheck([CHECKS[4]], usedefault = false)
        learner = testlearner(cb)
        learner.lossfn = x -> true
        @test_throws SanityCheckException (@suppress fit!(learner, 1))
    end
end


text:
# functions to fit models to data

using NLopt, DataFrames, Random, Printf

const tnd_sd_scale = 0.1 # SD(tnd) = tnd_sd_scale * <tnd> (as in Palmer, Huk & Shadlen, 2005)
#const tnd_sd_scale = 0 # SD(tnd) = tnd_sd_scale * <tnd>
const simsrand = 0

###############################################################################
# objective functions for fitting model parameters

# Bernoulli likelihood
plogq(p, q) = 0 .<= p .< 1e-300 ? -500.0p : p .* log.(max.(pmin, q))  # ensure 0 * log 0 = 0
#bernllh(m_pr, d_pr, d_n) = d_n .* (plogq(d_pr, m_pr) + plogq(1 - d_pr, 1 - m_pr))
bernllh(m_pr, d_pr, d_n) = d_n .* (d_pr .* log.(max.(pmin, m_pr)) .+ (1 .- d_pr) .* log.(max.(pmin, 1-m_pr)))
# Gaussian likelihood for data drawn from N(m_tμ, m_tσ^2)
gaussfullllh(m_tμ, m_tσ, d_tμ, d_tσ, d_n) = d_n .* (-0.5logtwoπ .- log.(m_tσ) .-
    (abs2.(d_tσ) .+ abs2.(d_tμ .- m_tμ)) ./ (2abs2.(m_tσ)))
# Gaussian likelihood for data mean drawn d_n times from N(m_tμ, d_tsem^2)
gausssemllh(m_tμ, d_tμ, d_tsem, d_n) = d_n .* (-0.5logtwoπ .- log.(d_tsem) .-
    abs2.(d_tμ .- m_tμ)./(2abs2.(d_tsem)))

# the objective function
abstract type FittingObjective end

###############################################################################
# objective functions that don't model sequential choices

# - fit choices by Bernoulli likelihood
# - fit RTs with Gaussian likelihood for correct / incorrect choices combined
# - only fit a subset or both conditions (identified at initialization)
# - the RT variance is augmented by (tnd_sd_scale * tnd)^2
struct PsychChron1ObjFn <: FittingObjective
    m::ModelBase
    evidence::Matrix{Float64}
    ids::Vector{Int}
    uniqueids::Vector{Int}
    burnin::Int
    d_psych::Vector{DataFrame}
    conditions::Vector{Int}

    PsychChron1ObjFn(m, evidence, ids, burnin, d_ident, d_categ, datasets,
                     conds::AbstractVector{Symbol}) =
        new(m, evidence, ids, sort(unique(ids)), burnin,
            DataFrame[by(df, :stimid, df -> _tcstats(df, datasets))
                      for df in [d_ident, d_categ]],
            map(x -> Dict(:ident => 1, :categ => 2)[x], conds))
    # defaults to fitting both conditions simultaneously
    PsychChron1ObjFn(m, evidence, ids, burnin, d_ident, d_categ, datasets) =
        PsychChron1ObjFn(m, evidence, ids, burnin, d_ident, d_categ, datasets,
                         [:ident, :categ])
end
# function returns avg. llh across modeled trials
function llh(o::PsychChron1ObjFn, ϕ::Vector)
    # simulate model and gather summary stats
    Random.seed!(simsrand)
    ms = simmodel(o.m, ϕ, o.evidence, o.ids)
    tndvar = abs2(tnd_sd_scale * gettnd(o.m, ϕ))
    # assemble likelihood
    llh, n = 0.0, 0.0
    for i in o.conditions
        d = o.d_psych[i]
        m = by(ms[i][(o.burnin+1):end, :], :stimid, _tcstats)
        for id in o.uniqueids
            m_prμ, m_tμ = Vector{Float64}(m[m[:stimid] .== id, [:prμ, :tμ]][1,:])
            d_n, d_prμ, d_tμ, d_tsem =
                Vector{Float64}(d[d[:stimid] .== id, [:n, :prμ, :tμ, :tsem]][1,:])
            llh += bernllh(m_prμ, d_prμ, d_n) +
                # we use d_n = 1, as we are already using the SEM
                gausssemllh(m_tμ, d_tμ, √(abs2(d_tsem) + tndvar), 1)
            n += d_n
        end
    end
    return llh / n
end
# alternatives: fit single condition
PsychChron1IdentObjFn(m, evidence, ids, burnin, d_ident, d_categ, datasets) =
    PsychChron1ObjFn(m, evidence, ids, burnin, d_ident, d_categ, datasets, [:ident])
PsychChron1CategObjFn(m, evidence, ids, burnin, d_ident, d_categ, datasets) =
    PsychChron1ObjFn(m, evidence, ids, burnin, d_ident, d_categ, datasets, [:categ])

# - fit choices by Bernoulli likelihood
# - fit RTs for Gaussian likelihood for correct / incorrect choices separately
# - the RT variance is augmented by (tnd_sd_scale * tnd)^2
struct PsychChron2ObjFn <: FittingObjective
    m::ModelBase
    evidence::Matrix{Float64}
    ids::Vector{Int}
    uniqueids::Vector{Int}
    burnin::Int
    d_psych::Vector{DataFrame}

    PsychChron2ObjFn(m, evidence, ids, burnin, d_ident, d_categ, datasets) =
        new(m, evidence, ids, sort(unique(ids)), burnin,
            DataFrame[by(df, :stimid, df -> _tcstats(df, datasets))
                      for df in [d_ident, d_categ]])
end
# function returns avg. llh across modeled trials
function llh(o::PsychChron2ObjFn, ϕ::Vector)
    # simulate model
    Random.seed!(simsrand)
    ms = simmodel(o.m, ϕ, o.evidence, o.ids)
    tndvar = abs2(tnd_sd_scale * gettnd(o.m, ϕ))
    # assemble likelihood
    llh, n = 0.0, 0.0
    for i in [1, 2]
        d = o.d_psych[i]
        m = by(ms[i][(o.burnin+1):end, :], :stimid, _tcstats)
        for id in o.uniqueids
            m_prμ, m_tcorrμ, m_tincorrμ = Vector{Float64}(
                m[m[:stimid] .== id, [:prμ, :tcorrμ, :tincorrμ]][1,:])
            d_n, d_ncorr, d_prμ, d_tcorrμ, d_tincorrμ, d_tcorrsem, d_tincorrsem =
                Vector{Float64}(
                    d[d[:stimid] .== id,
                      [:n, :ncorr, :prμ, :tcorrμ, :tincorrμ, :tcorrsem, :tincorrsem]][1,:])
            llh += bernllh(m_prμ, d_prμ, d_n) +
                # we use d_n = 1, as we are already using the SEM
                gausssemllh(m_tcorrμ, d_tcorrμ, √(abs2(d_tcorrsem) + tndvar), 1) +
                gausssemllh(m_tincorrμ, d_tincorrμ, √(abs2(d_tincorrsem) + tndvar), 1)
            n += d_n
        end
    end
    return llh / n
end

# - fit choices by Bernoulli likelihood
# - fit RTs with Gaussian SD likelihood for correct / incorrect choices combined
# - only fit a subset or both conditions (identified at initialization)
# (same as PsychChron2ObjFn, only that RTs are modelled by SD, not SEM)
struct PsychChron3ObjFn <: FittingObjective
    m::ModelBase
    evidence::Matrix{Float64}
    ids::Vector{Int}
    uniqueids::Vector{Int}
    burnin::Int
    d_psych::Vector{DataFrame}
    conditions::Vector{Int}

    PsychChron3ObjFn(m, evidence, ids, burnin, d_ident, d_categ, datasets,
                     conds::AbstractVector{Symbol}) =
        new(m, evidence, ids, sort(unique(ids)), burnin,
            DataFrame[by(df, :stimid, df -> _tcstats(df, datasets))
                      for df in [d_ident, d_categ]],
            map(x -> Dict(:ident => 1, :categ => 2)[x], conds))
    # defaults to fitting both conditions simultaneously
    PsychChron3ObjFn(m, evidence, ids, burnin, d_ident, d_categ, datasets) =
        PsychChron3ObjFn(m, evidence, ids, burnin, d_ident, d_categ, datasets,
                         [:ident, :categ])
end
# function returns avg. llh across modeled trials
function llh(o::PsychChron3ObjFn, ϕ::Vector)
    # simulate model and gather summary stats
    Random.seed!(simsrand)
    ms = simmodel(o.m, ϕ, o.evidence, o.ids)
    # assemble likelihood
    llh, n = 0.0, 0.0
    for i in o.conditions
        d = o.d_psych[i]
        m = by(ms[i][(o.burnin+1):end, :], :stimid, _tcstats)
        for id in o.uniqueids
            m_prμ, m_tμ, m_tσ =
                Vector{Float64}(m[m[:stimid] .== id, [:prμ, :tμ, :tσ]][1,:])
            d_n, d_prμ, d_tμ, d_tσ =
                Vector{Float64}(d[d[:stimid] .== id, [:n, :prμ, :tμ, :tσ]][1,:])
            llh += bernllh(m_prμ, d_prμ, d_n) +
                gaussfullllh(m_tμ, m_tσ, d_tμ, d_tσ, d_n)
            n += d_n
        end
    end
    return llh / n
end
# alternatives: fit single condition
PsychChron3IdentObjFn(m, evidence, ids, burnin, d_ident, d_categ, datasets) =
    PsychChron3ObjFn(m, evidence, ids, burnin, d_ident, d_categ, datasets, [:ident])
PsychChron3CategObjFn(m, evidence, ids, burnin, d_ident, d_categ, datasets) =
    PsychChron3ObjFn(m, evidence, ids, burnin, d_ident, d_categ, datasets, [:categ])

###############################################################################
# objective functions that model sequential choices

# - after correct choices, fit choices by Bernoulli likelihood, condition on
#   previous stimulus id
# - after incorrect choices, fit choices by Bernoulli likelihood, independent
#   of pervious stimulus
# - fit RTs with Gaussian likelihood for correct / incorrect choices combined
# - the RT variance is augmented by (tnd_sd_scale * tnd)^2
struct PsychChronSeq1ObjFn <: FittingObjective
    m::ModelBase
    evidence::Matrix{Float64}
    ids::Vector{Int}
    uniqueids::Vector{Int}
    burnin::Int
    d_nonseq::Vector{DataFrame}
    d_seq_corr::Vector{DataFrame}
    d_seq_incorr::Vector{DataFrame}

    PsychChronSeq1ObjFn(m, evidence, ids, burnin, d_ident, d_categ, datasets) =
        new(m, evidence, ids, sort(unique(ids)), burnin,
            DataFrame[by(df, :stimid, df -> _tcstats(df, datasets))
                      for df in [d_ident, d_categ]],
            DataFrame[by(df[df[:prevrew] .== true, :], [:stimid, :previd],
                         df -> _tcstats(df, datasets)) for df in [d_ident, d_categ]],
            DataFrame[by(df[df[:prevrew] .== false, :], :stimid,
                         df -> _tcstats(df, datasets)) for df in [d_ident, d_categ]])
end
function llh(o::PsychChronSeq1ObjFn, ϕ::Vector)
    Random.seed!(simsrand)
    ms = simmodel(o.m, ϕ, o.evidence, o.ids)
    tndvar = abs2(tnd_sd_scale * gettnd(o.m, ϕ))
    # assemble likelihood
    llh, n = 0.0, 0.0
    for i in [1, 2]
        mi = ms[i][(o.burnin+1):end, :]
        m = by(mi, :stimid, _tcstats)
        mcorr = by(mi[mi[:prevrew] .== true, :], [:stimid, :previd], _tcstats)
        if any(mi[:prevrew] .== false)
            mincorr = by(mi[mi[:prevrew] .== false, :], :stimid, _tcstats)
        else
            # create dataframe that violates any(mincorr[:stimid] .== id) below
            mincorr = DataFrame(stimid = Int[0])
        end
        d = o.d_nonseq[i]
        dcorr = o.d_seq_corr[i]
        dincorr = o.d_seq_incorr[i]
        # in all cases, iterate over current stimulus id
        for id in o.uniqueids
            # reaction times
            m_tμ = m[m[:stimid] .== id, :tμ][1]
            d_n, d_tμ, d_tsem =
                Vector{Float64}(d[d[:stimid] .== id, [:n, :tμ, :tsem]][1,:])
            # we use d_n = 1, as we are already using the SEM
            llh += gausssemllh(m_tμ, d_tμ, √(abs2(d_tsem) + tndvar), 1)
            # choices after incorrect trials
            if any(mincorr[:stimid] .== id)
                m_prμ = mincorr[mincorr[:stimid] .== id, :prμ][1]
            else
                m_prμ = id > 4 ? 1.0 : 0.0    # can happen for too high 'gain'
            end
            d_nincorr, d_prμ =
                Vector{Float64}(dincorr[dincorr[:stimid] .== id, [:n, :prμ]][1,:])
            llh += bernllh(m_prμ, d_prμ, d_nincorr)
            # choices after correct trials
            d_ncorr = 0.0
            for previd in o.uniqueids
                m_prμ = mcorr[(mcorr[:previd] .== previd) .&
                              (mcorr[:stimid] .== id), :prμ][1]
                d_ni, d_prμ = Vector{Float64}(
                    dcorr[(dcorr[:previd] .== previd) .&
                          (dcorr[:stimid] .== id), [:n, :prμ]][1,:])
                llh += bernllh(m_prμ, d_prμ, d_ni)
                d_ncorr += d_ni
            end
            # check if n is consistent
            @assert d_n == d_ncorr + d_nincorr
            n += d_n
        end
    end
    return llh / n
end

# - after correct choices, fit choices by Bernoulli likelihood and RTs with
#   Gaussian likelhood (indep. of corr/incorr), conditioned on previous stimid
# - after incorrect choices, do the same without conditioning on prev. stimid
# - the RT variance is augmented by (tnd_sd_scale * tnd)^2
#
# This model provides larger llh's than the above models, due to the scaling of
# the SEM with n. For this model, the SEM's will be unproportionally larger
# (based on less data) than for the above models, causing the llh's to be larger.
# Overall, this causes a stronger focus on p(choice) than RTs.
struct PsychChronSeq2ObjFn <: FittingObjective
    m::ModelBase
    evidence::Matrix{Float64}
    ids::Vector{Int}
    uniqueids::Vector{Int}
    burnin::Int
    d_seq_corr::Vector{DataFrame}
    d_seq_incorr::Vector{DataFrame}

    PsychChronSeq2ObjFn(m, evidence, ids, burnin, d_ident, d_categ, datasets) =
        new(m, evidence, ids, sort(unique(ids)), burnin,
            DataFrame[by(df[df[:prevrew] .== true, :], [:stimid, :previd],
                         df -> _tcstats(df, datasets)) for df in [d_ident, d_categ]],
            DataFrame[by(df[df[:prevrew] .== false, :], :stimid,
                         df -> _tcstats(df, datasets)) for df in [d_ident, d_categ]])
end
function llh(o::PsychChronSeq2ObjFn, ϕ::Vector)
    Random.seed!(simsrand)
    ms = simmodel(o.m, ϕ, o.evidence, o.ids)
    tndvar = abs2(tnd_sd_scale * gettnd(o.m, ϕ))
    # assemble likelihood
    llh, n = 0.0, 0.0
    for i in [1, 2]
        mi = ms[i][(o.burnin+1):end, :]
        mcorr = by(mi[mi[:prevrew] .== true, :], [:stimid, :previd], _tcstats)
        if any(mi[:prevrew] .== false)
            mincorr = by(mi[mi[:prevrew] .== false, :], :stimid, _tcstats)
        else
            # create dataframe that violates any(mincorr[:stimid] .== id) below
            mincorr = DataFrame(stimid = Int[0], tμ = Float64[mean(mcorr[:tμ])])
        end
        # llh for trials after correct choices
        dcorr = o.d_seq_corr[i]
        for previd in o.uniqueids
            for id in o.uniqueids
                m_prμ, m_tμ = Vector{Float64}(
                    mcorr[(mcorr[:previd] .== previd) .&
                          (mcorr[:stimid] .== id), [:prμ, :tμ]][1,:])
                d_n, d_prμ, d_tμ, d_tsem = Vector{Float64}(
                    dcorr[(dcorr[:previd] .== previd) .&
                          (dcorr[:stimid] .== id), [:n, :prμ, :tμ, :tsem]][1,:])
                llh += bernllh(m_prμ, d_prμ, d_n) +
                    # we use d_n = 1, as we are already using the SEM
                    gausssemllh(m_tμ, d_tμ, √(abs2(d_tsem) + tndvar), 1)
                n += d_n
            end
        end
        # llh for trials after incorrect choices
        dincorr = o.d_seq_incorr[i]
        for id in o.uniqueids
            if any(mincorr[:stimid] .== id)
                m_prμ, m_tμ = Vector{Float64}(
                    mincorr[mincorr[:stimid] .== id, [:prμ, :tμ]][1,:])
            else
                # for large p(correct), there might have been too few previous
                # incorrect choices - in this case, take the average tμ, and the
                # correct pr
                m_prμ = id > 4 ? 1.0 : 0.0
                m_tμ = mean(mincorr[:tμ])
            end
            d_n, d_prμ, d_tμ, d_tsem = Vector{Float64}(
                dincorr[dincorr[:stimid] .== id, [:n, :prμ, :tμ, :tsem]][1,:])
            llh += bernllh(m_prμ, d_prμ, d_n) +
                # we use d_n = 1, as we are already using the SEM
                gausssemllh(m_tμ, d_tμ, √(abs2(d_tsem) + tndvar), 1)
            n += d_n
        end
    end
    return llh / n
end

# - after correct choices, fit choices by Bernoulli likelihood, condition on
#   previous stimulus id
# - after incorrect choices, fit choices by Bernoulli likelihood, independent
#   of pervious stimulus
# - fit RTs with Gaussian likelihood for correct / incorrect choices combined
# (same as PsychChronSeq1ObjFn, only that RTs are modelled by SD, not SEM)
struct PsychChronSeq3ObjFn <: FittingObjective
    m::ModelBase
    evidence::Matrix{Float64}
    ids::Vector{Int}
    uniqueids::Vector{Int}
    burnin::Int
    d_nonseq::Vector{DataFrame}
    d_seq_corr::Vector{DataFrame}
    d_seq_incorr::Vector{DataFrame}

    PsychChronSeq3ObjFn(m, evidence, ids, burnin, d_ident, d_categ, datasets) =
        new(m, evidence, ids, sort(unique(ids)), burnin,
            DataFrame[by(df, :stimid, df -> _tcstats(df, datasets))
                      for df in [d_ident, d_categ]],
            DataFrame[by(df[df[:prevrew] .== true, :], [:stimid, :previd],
                         df -> _tcstats(df, datasets)) for df in [d_ident, d_categ]],
            DataFrame[by(df[df[:prevrew] .== false, :], :stimid,
                         df -> _tcstats(df, datasets)) for df in [d_ident, d_categ]])
end
function llh(o::PsychChronSeq3ObjFn, ϕ::Vector)
    Random.seed!(simsrand)
    ms = simmodel(o.m, ϕ, o.evidence, o.ids)
    # assemble likelihood
    llh, n = 0.0, 0.0
    for i in [1, 2]
        mi = ms[i][(o.burnin+1):end, :]
        m = by(mi, :stimid, _tcstats)
        mcorr = by(mi[mi[:prevrew] .== true, :], [:stimid, :previd], _tcstats)
        if any(mi[:prevrew] .== false)
            mincorr = by(mi[mi[:prevrew] .== false, :], :stimid, _tcstats)
        else
            # create dataframe that violates any(mincorr[:stimid] .== id) below
            mincorr = DataFrame(stimid = Int[0])
        end
        d = o.d_nonseq[i]
        dcorr = o.d_seq_corr[i]
        dincorr = o.d_seq_incorr[i]
        # in all cases, iterate over current stimulus id
        for id in o.uniqueids
            # reaction times
            m_tμ, m_tσ = Vector{Float64}(m[m[:stimid] .== id, [:tμ, :tσ]][1,:])
            d_n, d_tμ, d_tσ = Vector{Float64}(d[d[:stimid] .== id, [:n, :tμ, :tσ]][1,:])
            llh += gaussfullllh(m_tμ, m_tσ, d_tμ, d_tσ, d_n)
            # choices after incorrect trials
            if any(mincorr[:stimid] .== id)
                m_prμ = mincorr[mincorr[:stimid] .== id, :prμ][1]
            else
                m_prμ = id > 4 ? 1.0 : 0.0    # can happen for too high 'gain'
            end
            d_nincorr, d_prμ =
                Vector{Float64}(dincorr[dincorr[:stimid] .== id, [:n, :prμ]][1,:])
            llh += bernllh(m_prμ, d_prμ, d_nincorr)
            # choices after correct trials
            d_ncorr = 0.0
            for previd in o.uniqueids
                m_prμ = mcorr[(mcorr[:previd] .== previd) .&
                              (mcorr[:stimid] .== id), :prμ][1]
                d_ni, d_prμ = Vector{Float64}(
                    dcorr[(dcorr[:previd] .== previd) .&
                          (dcorr[:stimid] .== id), [:n, :prμ]][1,:])
                llh += bernllh(m_prμ, d_prμ, d_ni)
                d_ncorr += d_ni
            end
            # check if n is consistent
            @assert d_n == d_ncorr + d_nincorr
            n += d_n
        end
    end
    return llh / n
end

###############################################################################
# objective functions for the interleaved dataset

# - fit choices by Bernoulli likelihood
# - fit RTs with Gaussian likelihood for correct / incorrect choices combined
# - only fit a subset or both conditions (identified at initialization)
# - the RT variance is augmented by (tnd_sd_scale * tnd)^2
struct PsychChron1ItlvdObjFn <: FittingObjective
    m::ModelBase
    evidence::Matrix{Float64}
    ids::Vector{Int}
    burnin::Int
    d_psych::DataFrame
    modeledids::Vector{Int}

    PsychChron1ItlvdObjFn(m, evidence, ids, burnin, d_itlvd, datasets, modeledids) =
        new(m, evidence, ids, burnin,
            by(d_itlvd, :stimid, df -> _tcstats(df, datasets)), modeledids)
    # defaults to fitting all stimulus ids simultaneously
    PsychChron1ItlvdObjFn(m, evidence, ids, burnin, d_itlvd, datasets) =
        PsychChron1ItlvdObjFn(m, evidence, ids, burnin, d_itlvd, datasets,
                               collect(1:32))
end
# function returns avg. llh across modeled trials
function llh(o::PsychChron1ItlvdObjFn, ϕ::Vector)
    # simulate model and gather summary stats
    Random.seed!(simsrand)
    ms = simitlvdmodel(o.m, ϕ, o.evidence, o.ids)
    tndvar = abs2(tnd_sd_scale * gettnd(o.m, ϕ))
    # assemble likelihood
    llh, n = 0.0, 0.0
    d = o.d_psych
    m = by(ms[(o.burnin+1):end, :], :stimid, _tcstats)
    for id in o.modeledids
        m_prμ, m_tμ = Vector{Float64}(m[m[:stimid] .== id, [:prμ, :tμ]][1,:])
        d_n, d_prμ, d_tμ, d_tsem =
            Vector{Float64}(d[d[:stimid] .== id, [:n, :prμ, :tμ, :tsem]][1,:])
        #println("id $id : mprμ = $m_prμ  dprμ = $d_prμ    mtμ = $m_tμ  dtμ = $d_tμ")
        llh += bernllh(m_prμ, d_prμ, d_n) +
            # we use d_n = 1, as we are already using the SEM
            gausssemllh(m_tμ, d_tμ, √(abs2(d_tsem) + tndvar), 1)
        n += d_n
    end
    return llh / n
end

###############################################################################
# model fitting function

# returns all tested and best-fitting model parameters and objective value for given model
function fitmodel(m::ModelBase, o::FittingObjective, ϕini::AbstractVector; verbose=true)
    maxevals = 5000
    ϕmin, ϕmax = minimumϕ(m), maximumϕ(m)
    ϕn = length(ϕini)
    # store tested parameters + LLHs
    ϕs = Array{Float64}(undef, maxevals+1, ϕn+1)

    # formats output vector of floats
    fmtvec(x) = foldl((a, b) -> a * ", " * b, [@sprintf("%.2f", i) for i in x])

    # verbose likelihood function
    evalcount = -1  # start at -1 to not record first llhv(.) test call
    function llhv(ϕ, g)
        boundllh = llh(o, max.(ϕmin, min.(ϕmax, ϕ)))
        if verbose && mod(evalcount, 20) == 0
            @printf("%8d  ϕ = [%s]   llh = %f\n", evalcount, fmtvec(ϕ), boundllh)
        end
        evalcount += 1
        if evalcount >= 1
            ϕs[evalcount,:] = [ϕ; boundllh]
        end
        return boundllh
    end

    # Call the objective function once, to make sure it doesn't bug out.
    # We do this here, as within NLopt.optimize, it doesn't generate an stack
    # trace, which makes debugging impossible.
    !verbose || print("Testing llh function...")
    llhv(ϕini, Nothing)
    !verbose || println(" passed")

    # use subplex method, which is able to cope with noisy objective functions
    opt = Opt(:LN_SBPLX, ϕn)
    max_objective!(opt, llhv)
    lower_bounds!(opt, ϕmin)
    upper_bounds!(opt, ϕmax)
    ftol_rel!(opt, 1e-8)
    xtol_rel!(opt, 1e-32)
    maxtime!(opt, 11 * 60 * 60)  # 11h max
    maxeval!(opt, maxevals)     # max evaluations

    !verbose || println("    iter  result")
    optllh, optϕ, ret = NLopt.optimize(opt, ϕini)
    # add final value as last element
    evalcount += 1
    ϕs[evalcount,:] = [optϕ; optllh]

    # return dataframe
    return DataFrame(Any[ϕs[1:evalcount,i] for i in 1:ϕn+1], [ϕnames(m); :llh]), ret
end
fitmodel(m::ModelBase, o::FittingObjective; verbose=true) = fitmodel(m, o, initialϕ(m); verbose=verbose)


text:
include("BaseAPI.jl")
include("TournamentAPI.jl")
include("SignalsAPI.jl")

# include("Staking.jl")


text:
@testset "Small graphs" begin
    @testset "by string" begin
        g = smallgraph("diamondgraph")
        @test nv(g) == 4 && ne(g) == 5

        g = smallgraph("diamond")
        @test nv(g) == 4 && ne(g) == 5
    end

    smallgraphs = [
        (:diamond        , 4, 5, false),
        (:bull           , 5, 5, false),
        (:chvatal        , 12, 24, false),
        (:cubical        , 8, 12, false),
        (:desargues      , 20, 30, false),
        (:dodecahedral   , 20, 30, false),
        (:frucht         , 12, 18, false),
        (:heawood        , 14, 21, false),
        (:house          , 5,6, false),
        (:housex         , 5, 8, false),
        (:icosahedral    , 12, 30, false),
        (:krackhardtkite , 10, 18, false),
        (:moebiuskantor  , 16, 24, false),
        (:octahedral     , 6, 12, false),
        (:pappus         , 18, 27, false),
        (:petersen       , 10, 15, false),
        (:sedgewickmaze  , 8, 10, false),
        (:tetrahedral    , 4, 6, false),
        (:truncatedcube  , 24, 36, false),
        (:truncatedtetrahedron , 12, 18, false),
        (:truncatedtetrahedron_dir , 12, 18, true),
        (:tutte          , 46, 69, false)
    ]

    @testset "$(s[1])" for s in smallgraphs
        g = smallgraph(s[1])
        (nvg, neg, dir) = s[2:4]
        @test nv(g) == nvg
        @test ne(g) == neg
        @test is_directed(g) == dir
    end

    @testset "karate" begin
        g = smallgraph(:karate)
        degree_sequence = sort([16, 9, 10, 6, 3, 4, 4, 4, 5, 2, 3, 1, 2, 5, 2,  2,  2,
                                 2, 2,  3, 2, 2, 2, 5, 3, 3, 2, 4, 3, 4, 4, 6, 12, 17])
        @test nv(g) == 34 && ne(g) == 78 && sort(degree(g)) == degree_sequence
    end
    @testset "nonexistent graph" begin
        @test_throws ArgumentError g = smallgraph(:nonexistent)
    end
end


text:
function layError(proj::projectile, target::AbstractTarget)
    range = euclidean(proj.position, target.position)
    lay = 300/range+0.05
end

function randomError(σdisp::Error,σaiming::Error)
    σrandom=Error(0.0,0.0)
    σrandom.horizontal = sqrt(σdisp.horizontal^2+σaiming.horizontal^2)

    σrandom.vertical = sqrt(σdisp.vertical^2+σaiming.vertical^2)

    return σrandom

end