text:
# test mercer kernels and RKHS theory.

using LinearAlgebra
import Random

import GSL
import PyPlot
import Interpolations
import Distributions
using SpecialFunctions

import HCubature

import SignalTools

include("../src/misc/declarations.jl")

include("../src/RKHS/RKHS.jl")
include("../src/RKHS/kernel.jl")
include("../src/RKHS/interpolators.jl")
include("../src/RKHS/mercer.jl")

include("../src/warp_map/Rieszanalysis.jl")
include("../src/misc/utilities.jl")

PyPlot.close("all")

Random.seed!(25)

fig_num = 1

##### regression
θ = GaussianKernel1DType(1.1)

σ² = 1e-3
N = 14

x_range = LinRange(-π,2.34*π,N)
x_range = LinRange(-π-4, 2.34*π+4,N)

X = collect( [x_range[n]] for n = 1:N )

### specify oracle function.
c = randn(N)
f_hat = xx->sum( c[j]*evalkernel(xx, X[j], θ)^2 for j = 1:N )

# query locations.
Nq = 1000
xq_range = LinRange(-π-4, 2.34*π+2, Nq)
xq = collect( [xq_range[n]] for n = 1:Nq )

### set up mercer expansion.

# parameters.
M = 100 #  this is the truncation limit for the Mercer expansion.
α = 1.0

# pre-compute.
ϵ, β, δ_sq = setupmercer(θ, α)
λ_array = geteigenvalues(θ, M, α)

# unnormalized eigenfunction.
ϕ_tilde = (ii,xx)->sqrt(λ_array[ii])*evaleigenfunction(θ, α, β, δ_sq, ii, xx[1])

# normalized eigenfunction.
ϕ = (ii,xx)->evaleigenfunction(θ, α, β, δ_sq, ii, xx[1])

### verify Mercer expansion.
#x = randn()
#z = randn()
x = 11.0
z = 14.0

ϕ_x_array = collect( ϕ(m,x) for m = 1:M )
ϕ_z_array = collect( ϕ(m,z) for m = 1:M )

k_x_z = evalkernel([x],[z],θ)
mercer_sum = sum(λ_array[n]*ϕ_x_array[n]*ϕ_z_array[n] for n = 1:M)
println("k(x,z) = ", k_x_z, ", mercer sum = ", mercer_sum)
println("discrepancy = ", abs(k_x_z-mercer_sum) )
println()

### verify kernel to SOS.
Q = getQmat(c, ϕ_tilde, X, M)

f_hat2 = xx->sum( Q[i,j] *evaleigenfunction(θ, α, β, δ_sq, i, xx)
                *evaleigenfunction(θ, α, β, δ_sq, j, xx)
                *sqrt(λ_array[i])*sqrt(λ_array[j]) for i = 1:M for j = 1:M )

#
y_Q_SOS = f_hat2.(xq_range)
y_kernel_SOS = f_hat.(xq)
discrepancy = norm( y_Q_SOS - y_kernel_SOS )
println("discrepancy between kernel and Q SOS expressions = ", discrepancy)

# # Visualize regression result.
PyPlot.figure(fig_num)
fig_num += 1

PyPlot.plot(xq_range, y_Q_SOS, ".", label = "y_Q_SOS")
PyPlot.plot(xq_range, y_kernel_SOS, label = "y_kernel_SOS")

title_string = "SOS RKHS - Gaussian Mercer series"
PyPlot.title(title_string)
PyPlot.legend()

# seems like for x > 9.4, things are so ill-conditioned that the kernel and
#   Q representations no longer agree. We could do log-space operations, but
#   this still mean the SDP problem, which involves a ill-conditioned kernel matrix,
#   is difficult to enforce to be posdef, at least numerically.

# continue investigating, just avoid using Q.

@assert 1==2

#### orthogonality check.
μ = xx->α/sqrt(π)*exp(-α^2*xx^2)
val, err = evalindefiniteintegral(μ)
println("integral of μ is ", val)

i = 70
j = 9
ϵ, β, δ_sq = setupmercer(θ, α)
#integrand_func = xx->μ(xx)*evaleigenfunction(θ, α, β, δ_sq, i, xx)*evaleigenfunction(θ, α, β, δ_sq, j, xx)
#val, err = evalindefiniteintegral(integrand_func)

#N_MC = 2000000
N_MC = 20000
σ = 1/(sqrt(2)*α)
measure_dist = Distributions.Normal(0.0, σ)
X = collect( rand(measure_dist) for n = 1:N_MC )
MC_func = xx->evaleigenfunction(θ, α, β, δ_sq, i, xx)*evaleigenfunction(θ, α, β, δ_sq, j, xx)

# check orthogonality of eigenfunctions wrt to the measure μ.
MC_integral = sum( MC_func(X[n])/N_MC for n = 1:N_MC)

### explore RKHS norm, and the nature of ℋ. Verdict: the norm is a series, not integration.
k_z = xx->evalkernel(xx, [z], θ)
g = xx->sum( η.c[j]*evalkernel(xx, η.X[j], θ) for j = 1:length(η.c) )

N_MC = 20000
MC_norm_g = sum( g( [X[n]] )^2/N_MC for n = 1:N_MC)

K = constructkernelmatrix(η.X, η.θ)
RKHS_norm_g = dot(η.c, K*η.c)
#RKHS_norm_g = sum( dot( η.c[j], g([η.X[j]]) ) for j = 1:length(η.c) )
println("RKHS norm(g) = ", RKHS_norm_g, ", MC_norm_g = ", MC_norm_g)

### explore SOS.
x = randn()
M = 100 #  this is the truncation limit for the Mercer expansion.
V = randn(M,M)
V = V./norm(V)
ϕ_x = evaleigenfunctions(θ, M, α, x)

f = xx->evaleigenfunctions(θ, M, α, xx)'*V'*V*evaleigenfunctions(θ, M, α, xx)
q = xx->evaleigenfunctions(θ, M, α, xx)

y = f.(x_range)
f_xq = f.(xq_range)

# # Visualize regression result.
PyPlot.figure(fig_num)
fig_num += 1

PyPlot.plot(x_range, y, ".", label = "observed")
#PyPlot.plot(xq, yq, label = "query - RKHS")

PyPlot.plot(xq_range, f_xq, label = "true")

title_string = "SOS RKHS - Gaussian Mercer series"
PyPlot.title(title_string)
PyPlot.legend()

max_stars_repo_name: RoyCCWang/PatchMixtureKriging

max_issues_repo_name: RoyCCWang/PatchMixtureKriging


text:
using ProximalOperators

function agd(X, y, x0, f, g, eps, λ, maxiter, stepsize, num_data,prox_operator)

    x_k = x0
    x_old = x0
    oval = Inf

    val_array = zeros(0)
    stop_array = zeros(0)

    gval = g(x_k)
    for i in 1:maxiter

        # Logistic-regression with l2-norm of lambda
        fval = f(x_k)

        y_k = x_k + (i-1) / (i+2) * (x_k-x_old)
        gval = g(y_k)
        z_k = y_k - stepsize * gval

        x_old = x_k
        x_k, f_k = prox(prox_operator, z_k, stepsize)
        oval = fval

        stoppingCriteria = norm(gval,Inf)
        append!( val_array, fval )
        append!( stop_array, stoppingCriteria )
        @printf("%i \t\t %f \t\t %f \t %e \n",i,fval,(oval-fval),stoppingCriteria)

    end
    return (x_k,"iter",val_array,stop_array)
end

ext: jl

max_stars_repo_name: zagazao/num-opt

max_issues_repo_name: zagazao/num-opt

max_forks_repo_name: zagazao/num-opt


text:
## this should be in SymEngine/src/subs.jl

function Base.call(ex::Basic, args...)
    xs = free_symbols(ex)
    subs(ex, collect(zip(xs, args))...)
end

## function (ex::Basic)(args...)
##     xs = SymEngine.free_symbols(ex)
##     subs(ex, collect(zip(xs, args))...)
## end


text:
using PSFModels
using Documenter

setup = quote
    using PSFModels
end

DocMeta.setdocmeta!(PSFModels, :DocTestSetup, setup; recursive = true)

makedocs(;
    modules=[PSFModels],
    authors="Miles Lucas <mdlucas@hawaii.edu> and contributors",
    repo="https://github.com/JuliaAstro/PSFModels.jl/blob/{commit}{path}#L{line}",
    sitename="PSFModels.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://juliaastro.github.io/PSFModels.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Examples" => "examples.md",
        "Benchmarks" => "bench.md",
        "API/Reference" => "api.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaAstro/PSFModels.jl",
    push_preview=true,
    devbranch="main"
)


text:
# All the simulations
state_list=[mono,
            ps,
            sastre,
            bbc,
            sid,
            mono_org,
            ps_org,
            sastre_org,
            bbc_org,
            sid_org
            ];

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
using StochasticDiffEq, Test, Random
using DiffEqGPU
using CUDA

function weak_error(prob,alg,numtraj, batchsize, f_true,trange;abstol=1,reltol=0,ensemblealg=EnsembleCPUArray())
  sol = @time solve(prob,alg,ensemblealg,
    dt=0.01f0,adaptive=true,abstol=abstol,reltol=reltol,
    trajectories=numtraj,batch_size=batchsize,
    saveat = trange
    )
  computed_exp = (sol.u/numtraj)[1,:]
  true_exp = f_true.(trange)
  return sum((computed_exp-true_exp).^2)/length(trange)
end

function prob_func(prob, i, repeat)
    #(i%50000 == 0) && @show i
    remake(prob,seed=seeds[i])
end

function reduction(u,batch,I)
  u.+sum(batch),false
end

function output_func(sol,i)
  #h1(asinh(sol[end][1])),false
  h1.(asinh.(sol)),false
end

# prob 1
u₀ = [0.0f0]
tspan = (0.0f0,2.0f0)
h1(z) = z^3-6*z^2+8*z

function f1!(du,u,p,t)
 @inbounds begin
     du[1] = 1//2*u[1]+sqrt(u[1]^2 +1)
 end
 nothing
end

function g1!(du,u,p,t)
 @inbounds begin
     du[1] = sqrt(u[1]^2 +1)
 end
 nothing
end

f_true1(t) = t^3-3*t^2+2*t

prob1 = SDEProblem(f1!,g1!,u₀,tspan)
ensemble_prob1 = EnsembleProblem(prob1;
        output_func = output_func,
        prob_func = prob_func,
        reduction = reduction,
        u_init=Vector{eltype(prob1.u0)}([0.0])
        )

# prob 2
u₀ = [0.1f0,0.1f0]
function f2!(du,u,p,t)
  @inbounds begin
    du[1] = 3//2*u[1]
    du[2] = 3//2*u[2]
  end
  nothing
end
function g2!(du,u,p,t)
  @inbounds begin
    du[1] = 1//10*u[1]
    du[2] = 1//10*u[2]
  end
  nothing
end

f_true2(t) = 1//10*exp(3//2*t) #1//100*exp(301//100*t)

h2(z) = z #z^2

prob2 = SDEProblem(f2!,g2!,u₀,tspan)
ensemble_prob2 = EnsembleProblem(prob2;
        output_func = (sol,i) -> (h2.(sol),false),
        prob_func = prob_func,
        reduction = reduction,
        u_init=Vector{eltype(prob2.u0)}([0.0, 0.0])
        )

tsave = 0.0f0:0.05f0:2f0

probs = Vector{EnsembleProblem}(undef, 2)
probs[1] = ensemble_prob1
probs[2] = ensemble_prob2

ftrue = Vector{}(undef, 2)
ftrue[1] = f_true1
ftrue[2] = f_true2

numtraj = Int(1e5)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

for i in 1:2
  @show i

  err1 = weak_error(probs[i],DRI1NM(),numtraj,Int(1e4),ftrue[i],tsave,abstol=1f0,reltol=1f0, ensemblealg=EnsembleGPUArray())
  @show err1
  # err2 = weak_error(probs[i],DRI1NM(),numtraj,Int(1e4),ftrue[i],tsave,abstol=0.1f0,reltol=0.1f0, ensemblealg=EnsembleGPUArray())
  # @show err2
  err3 = weak_error(probs[i],DRI1NM(),numtraj,Int(1e4),ftrue[i],tsave,abstol=0.01f0,reltol=0.01f0, ensemblealg=EnsembleGPUArray())
  @show err3
  @test err1 > err3
  println("")
end

#
# sol = @time solve(probs[1],DRI1NM(),EnsembleGPUArray(),
#   dt=0.001f0,adaptive=false,abstol=0.1f0,reltol=0.1f0,
#   trajectories=numtraj,batch_size=Int(1e1),
#   saveat = tsave
#   )
# computed_exp = (sol.u/numtraj)[1,:]
# true_exp = f_true1.(tsave)
#
# @show  sum((computed_exp-true_exp).^2)/length(tsave)
#
# using Plots; plot(tsave, true_exp)
# plot!(tsave, computed_exp)

max_stars_repo_name: JuliaDiffEq/StochasticDiffEq


text:
function ellec_bulirsch(k2::Real)
ca = sqrt(eps(k2))
if k2 == 1.0
  return 1.0
end
m=1.0
a=1.0
b=1.0-k2
kc=sqrt(b)
c=a
a+=b

b=2.0*(c*kc+b)
c=a
m0=m
m+=kc
a+=b/m
while abs(m0-kc) > (ca*m0)
  kc=2.0*sqrt(kc*m0)
  b=2.0*(c*kc+b)
  c=a
  m0=m
  m+=kc
  a+=b/m
#  println(a,' ',b,' ',c,' ',m0,' ',kc)
end
return pi*0.25*a/m
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
using Codecs
using Base.Test
using Compat

include("tests.jl")


text:
const x = 1
x = π # ERROR: invalid ridefinition of constant x

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Variational                       #
# --------------------------------- #
@rule MvNormalMeanPrecision(:Λ, Marginalisation) (q_out::Any, q_μ::Any) = begin
    m_out, v_out   = mean_cov(q_out)
    m_mean, v_mean = mean_cov(q_μ)

    df = ndims(q_μ) + 2
    S  = Matrix(Hermitian(cholinv(v_mean + v_out + (m_mean - m_out) * (m_mean - m_out)')))

    return Wishart(df, S)
end


text:
# problem set 8 solutions

using Random
using LinearAlgebra
using Statistics
using Distributions
using Optim
using DataFrames
using DataFramesMeta
using CSV
using HTTP
using GLM
using FreqTables
using ForwardDiff
using LineSearches
using SMM
using MultivariateStats

include("lgwt.jl")

function wrapper()
    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    # question 1 - read in data and simple regression
    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS8-factor/nlsy.csv"
    df = CSV.read(HTTP.get(url).body)
    est1 = lm(@formula(logwage ~ black + hispanic + female + schoolt + gradHS + grad4yr), df)
    println(est1)

    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    # question 2 - compute correlation of  ASVAB scores
    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    asvabs = convert(Array,df[:,end-5:end])
    cors = cor(asvabs)
    cordf = DataFrame(cor1 = cors[:,1], cor2 = cors[:,2], cor3 = cors[:,3], cor4 = cors[:,4], cor5 = cors[:,5], cor6 = cors[:,6])
    display(cordf)
    # these are all pretty highly correlated; some more than others

    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    # question 3 - add ASVAB scores to regression
    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    est2 = lm(@formula(logwage ~ black + hispanic + female + schoolt + gradHS + grad4yr + asvabAR + asvabCS + asvabMK + asvabNO + asvabPC + asvabWK), df)
    println(est2)

    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    # question 4 - PCA of ASVAB scores
    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    asvabMat = convert(Array,df[:,end-5:end])'
    M = fit(PCA, asvabMat; maxoutdim=1)
    asvabPCA = MultivariateStats.transform(M, asvabMat)
    df = @transform(df, asvabPCA = asvabPCA[:])
    est3 = lm(@formula(logwage ~ black + hispanic + female + schoolt + gradHS + grad4yr + asvabPCA), df)
    println(est3)

    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    # question 5 - Factor Analysis of ASVAB scores
    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    asvabMat = convert(Array,df[:,end-5:end])'
    M = fit(FactorAnalysis, asvabMat; maxoutdim=1)
    asvabFac = MultivariateStats.transform(M, asvabMat)
    df = @transform(df, asvabFac = asvabFac[:])
    est4 = lm(@formula(logwage ~ black + hispanic + female + schoolt + gradHS + grad4yr + asvabFac), df)
    println(est4)

    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    # question 6 - Estimate full factor model by MLE
    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    X = [df.black df.hispanic df.female df.schoolt df.gradHS df.grad4yr ones(size(df,1),1)]
    y = df.logwage
    Xfac = [df.black df.hispanic df.female ones(size(df,1),1)]
    asvabs = [df.asvabAR df.asvabCS df.asvabMK df.asvabNO df.asvabPC df.asvabWK]

    # likelihood obj function
    function factor_model(θ, X, Xfac, Meas, y, R)
        K = size(X,2)    # no. covariates in outcome equation
        L = size(Xfac,2) # no. covariates in measurement eqns
        J = size(Meas,2) # no. measurements
        N = length(y)    # no. observations

        # organize the parameter vector
        # first J*L elements are the coefficients of the measurement system
        γ = reshape(θ[1:J*L],L,J)
        # next K elements are the coefficients of the outcome equation
        β = θ[J*L+1:J*L+K]
        # next J+1 elements are the factor loadings
        α = θ[J*L+K+1:J*L+K+J+1]
        # last J+1 elements are the variances (stdevs)
        σ = θ[end-J:end]

        # get grid nodes and weights for integral
        ξ,ω = lgwt(R,-5,5)

        # now build the likelihood
        T = promote_type(eltype(X),eltype(θ))
        like  = zeros(T,N)
        for r=1:R
            Mlike = zeros(T,N,J)
            for j=1:J
                Mres        = Meas[:,j] .- (Xfac*γ[:,j] .+ α[j]*ξ[r])
                sdj         = sqrt(σ[j]^2)
                Mlike[:,j] .= (1 ./ sdj).*pdf.(Normal(0,1),Mres./sdj)
            end
            Yres  = y .- (X*β .+ α[end]*ξ[r])
            sdy   = sqrt(σ[end]^2)
            Ylike = (1 ./ sdy).*pdf.(Normal(0,1),Yres./sdy)

            like += ω[r].*prod(Mlike; dims=2).*Ylike.*pdf(Normal(0,1),ξ[r])
        end
        loglike = -sum(log.(like))
        return loglike
    end

    # SMM objective function
    function factor_model_smm(θ, X, Xfac, Meas, y, D)
        K = size(X,2)    # no. covariates in outcome equation
        L = size(Xfac,2) # no. covariates in measurement eqns
        J = size(Meas,2) # no. measurements
        N = length(y)    # no. observations

        # organize the parameter vector
        # first J*L elements are the coefficients of the measurement system
        γ = reshape(θ[1:J*L],L,J)
        # next K elements are the coefficients of the outcome equation
        β = θ[J*L+1:J*L+K]
        # next J+1 elements are the factor loadings
        α = θ[J*L+K+1:J*L+K+J+1]
        # last J+1 elements are the variances (stdevs)
        σ = θ[end-J:end]

        # simulate the model given the parameter values
        Random.seed!(1234)
        T = promote_type(eltype(X),eltype(θ))
        M̃ = zeros(T,N,J)
        ỹ = zeros(T,N)
        for d=1:D
            ξ = randn(N)
            for j=1:J
                M̃[:,j] .+= (Xfac*γ[:,j] .+ α[j].*ξ .+ sqrt(σ[j].^2).*randn(N)).*(1/D)
            end
            ỹ .+= (X*β .+ α[end].*ξ .+ sqrt(σ[end].^2).*randn(N)).*(1/D)
        end
        g = vcat(Meas[:],y) .- vcat(M̃[:],ỹ)
        J = g'*I*g
        return J
    end
    svals = vcat(Xfac\asvabs[:,1],
                 Xfac\asvabs[:,2],
                 Xfac\asvabs[:,3],
                 Xfac\asvabs[:,4],
                 Xfac\asvabs[:,5],
                 Xfac\asvabs[:,6],
                 X\y,
                 rand(7),
                 .5*ones(7)
                )
    stata = vcat(-.9343684, -.6783938, -.0167224, .362379, -.4309737, -.4427523, .399414, -.003017, -.8084006, -.6365506, .1950629, .2234858, -.3418862, -.4675659, .2531868, .0513019, -.7595222, -.5958526, .2988241, .1545899, -.9353962, -.750025, .0496955, .3432385, -.1870249, -.0805532, -.1349807, -.0024857, .1064984, .2389176, 2.154472, .7704785, .5507262, .816202, .6431685, .7398663, .6622146, .1051398, sqrt(.2306496),sqrt(.6136927),sqrt(.1903943),sqrt(.5309717),sqrt(.3110141),sqrt(.3743002),sqrt(.2188733))
    # run the optimizer for MLE
    td = TwiceDifferentiable(th -> factor_model(th, X, Xfac, asvabs, y, 9), svals; autodiff = :forward)
    println("optimizing quadrature full factor model")
    θ̂_optim_ad = optimize(td, svals, Newton(linesearch = BackTracking()), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=1))
    θ̂_mle_optim_ad = θ̂_optim_ad.minimizer
    loglikeval = θ̂_optim_ad.minimum
    # evaluate the Hessian at the estimates
    H  = Optim.hessian!(td, θ̂_mle_optim_ad)
    θ̂_mle_optim_ad_se = sqrt.(diag(inv(H)))

    # run the optimizer for MLE from Stata starting values
    td = TwiceDifferentiable(th -> factor_model(th, X, Xfac, asvabs, y, 9), stata; autodiff = :forward)
    println("optimizing quadrature full factor model")
    θ̂_optim_ad_fr_stata = optimize(td, stata, Newton(linesearch = BackTracking()), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=1))
    θ̂_mle_optim_ad_fr_stata = θ̂_optim_ad_fr_stata.minimizer
    loglikeval_fr_stata = θ̂_optim_ad_fr_stata.minimum
    # evaluate the Hessian at the estimates
    H  = Optim.hessian!(td, θ̂_mle_optim_ad)
    θ̂_mle_optim_ad_se_fr_stata = sqrt.(diag(inv(H)))

    # run the optimizer for SMM from Stata starting values
    td = TwiceDifferentiable(th -> factor_model_smm(th, X, Xfac, asvabs, y, 3_000), stata; autodiff = :forward)
    println("optimizing SMM full factor model")
    θ̂_sm_optim_ad = optimize(td, stata, Newton(linesearch = BackTracking()), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=1))
    θ̂_smm_optim_ad = θ̂_sm_optim_ad.minimizer
    errval = θ̂_sm_optim_ad.minimum
    # evaluate the Hessian at the estimates
    H  = Optim.hessian!(td, θ̂_smm_optim_ad)
    θ̂_smm_optim_ad_se = sqrt.(diag(inv(H)))

    # print output
    eqs = vcat(fill("asvabAR",size(Xfac,2),1),
               fill("asvabCS",size(Xfac,2),1),
               fill("asvabMK",size(Xfac,2),1),
               fill("asvabNO",size(Xfac,2),1),
               fill("asvabPC",size(Xfac,2),1),
               fill("asvabWK",size(Xfac,2),1),
               fill("wage",size(X,2),1),
               repeat(["asvabAR", "asvabCS", "asvabMK", "asvabNO", "asvabPC", "asvabWK", "wage"],2,1),
              "overall")
    namer = vcat(repeat(["black","hispanic","female","intercept"],size(asvabs,2),1),
                 ["black","hispanic","female","yrs school","gradHS","grad4yr","intercept"],
                 fill("loading",size(asvabs,2)+1,1),
                 fill("stdev",size(asvabs,2)+1,1),
                "loglike")
    results = DataFrame(equation = vec(eqs), variable = vec(namer), coef_mle = vcat(vec(θ̂_mle_optim_ad),loglikeval), se_mle = vcat(vec(θ̂_mle_optim_ad_se),missing), coef_mle_fr_stata = vcat(vec(θ̂_mle_optim_ad_fr_stata),loglikeval_fr_stata), se_mle_fr_stata = vcat(vec(θ̂_mle_optim_ad_se_fr_stata),missing), coef_smm = vcat(vec(θ̂_smm_optim_ad),missing), se_smm = vcat(vec(θ̂_smm_optim_ad_se),missing),coef_stata = vcat(vec(stata),17524.022))
    println(results)

    return nothing
end
wrapper()


text:
############################################################################################
"""
$(TYPEDEF)
NUTS - Container used throughout sampling process.

# Fields
$(TYPEDFIELDS)
"""
mutable struct NUTS{R<:ℓObjectiveResult,D<:AbstractDifferentiableTune,C<:KineticEnergy} <:
               MCMCKernel
    "Target result stored for caching."
    result::R
    "Differentiation tuning struct."
    diff::D
    "Energy used for Hamiltonian."
    energy::C
    function NUTS(
        result::R, diff::D, energy::C
    ) where {R<:ℓObjectiveResult,D<:AbstractDifferentiableTune,C<:KineticEnergy}
        return new{R,D,C}(result, diff, energy)
    end
end
function NUTS(sym::S; kwargs...) where {S<:Union{Symbol,NTuple{k,Symbol} where k}}
    return MCMC(NUTS, sym; kwargs...)
end

function update!(kernel::NUTS, objective::Objective, up::BaytesCore.UpdateTrue)
    ## Update log-target result with current (latent) data
    kernel.diff = update(kernel.diff, objective)
    kernel.result = ModelWrappers.log_density_and_gradient(objective, kernel.diff)
    return nothing
end
function update!(kernel::NUTS, objective::Objective, up::BaytesCore.UpdateFalse)
    return nothing
end

############################################################################################
"""
$(SIGNATURES)
No-U-turn Hamiltonian Monte Carlo transition, using Hamiltonian `H`, starting at evaluated
log density position `Q`, using stepsize `ϵ`. Parameters of `kernel` govern the details
of tree construction. Return two values, the new evaluated log density position, and tree statistics.

# Examples
```julia
```

"""
function sample_tree(
    _rng::Random.AbstractRNG,
    trajectory::TrajectoryNUTS,
    phasepoint::PhasePoint;
    directions=rand(_rng, Directions),
)
    phasepointᵖ, v, termination, depth = sample_trajectory(
        _rng, trajectory, phasepoint, directions
    )
    tree_statistics = DiagnosticsNUTS(
        ℓdensity(trajectory.H, phasepointᵖ),
        depth,
        termination,
        acceptance_rate(v),
        trajectory.ϵ,
        v.steps,
        directions,
    )
    return phasepointᵖ, tree_statistics
end

############################################################################################
"""
$(SIGNATURES)
Propagate forward one proposal step.

# Examples
```julia
```

"""
function propagate(
    _rng::Random.AbstractRNG, kernel::NUTS, tune::MCMCTune, objective::Objective
)
    @unpack ϵ = tune.stepsize
    ## Set new DiffObjective
    diff = DiffObjective(objective, kernel.diff)
    ## If needed, update kinetic energy
    update!(kernel.energy, tune.proposal)
    H = Hamiltonian(kernel.energy, diff)
    ## Create new trajectory and phasepoint, and evaluate Hamiltonian at phasepoint
    phasepoint = PhasePoint(kernel.result, rand_ρ(_rng, H.K))
    trajectory = TrajectoryNUTS(H, ℓdensity(H, phasepoint), ϵ)
    ## Sample from NUTS trajectory
    phasepointᵖ, sampler_statistic = sample_tree(_rng, trajectory, phasepoint)
    ## Pack container and return output
    return phasepointᵖ.result,
    is_divergent(sampler_statistic.termination),
    BaytesCore.AcceptStatistic(sampler_statistic.acceptance_rate, true),
    sampler_statistic
end

############################################################################################
"""
$(SIGNATURES)
Tune initial step size.

# Examples
```julia
```

"""
function get_acceptrate(
    _rng::Random.AbstractRNG, kernel::NUTS, objective::Objective, Σ::M
) where {M<:AbstractMatrix} # LinearAlgebra.Diagonal(ones(length(objective.tagged)))
    ## Update Hamiltonian with current Σ and log-posterior adjusted for (latent/new) data
    H = Hamiltonian(kernel.energy, DiffObjective(objective, kernel.diff))
    ## Create new trajectory and phasepoint, and evaluate Hamiltonian at phasepoint
    phasepoint = PhasePoint(kernel.result, rand_ρ(_rng, H.K))
    H₀ = ℓdensity(H, phasepoint)
    ## Function of stepsize
    return function acceptrate(ϵ::T) where {T<:Real}
        # Make leapfrog step
        phasepointᵖ = leapfrog(H, phasepoint, ϵ)
        # Calculate acceptance rate
        return exp(ℓdensity(H, phasepointᵖ) - H₀)
    end
end

############################################################################################
# Export
export NUTS, sample_tree


text:
using Revise; using BenchmarkTools; using FunctionWrappers:FunctionWrapper;

struct TT
   Ga::Function
end

struct TTFW{T<:FunctionWrapper}
   Ga::T
end

function FWftest(G::Expr)
   p=(:a1, :a2)
   unpack_expr = Meta.parse(string(join(string.(p), ","), ", = a"))

   @eval return(TTFW(((ta, a) -> ($unpack_expr; [$G for t in ta])) |> FunctionWrapper{Vector{Float64}, Tuple{Vector{Float64}, Vector{Float64}}}))
end

macro FUmtest(z)
   p=(:a1, :a2)
   unpack_expr = Meta.parse(string(join(string.(p), ","), ", = a"))

   return quote TT((ta, a) -> ($unpack_expr; [$z for t in ta])) end
end

function localscopetesting()
   testvector = Vector(0.0:1.0:10.0)
   params = [1.0, 1.0];

   # FW1 = FWftest(:(a1*t + a2))
   # FW2 = FWftest(:(a1 + a2))
   # println(FW1.Ga(testvector, params))
   # println(FW2.Ga(testvector, params))

   # println("\n")

   FUM1 = @FUmtest a1*t+a2
   # FUM2 = @FUmtest a1+a2
   println(FUM1.Ga(testvector, params))
   # println(FUM2.Ga(testvector, params))

   #=
   FW1 (FunctionWrappers) works without declaring a global (FunctionWrapper compiled within current world age)
   FUM Works fine because it is literally a macro expansion that happens at parse time, no need to declare global
   =#
end

# localscopetesting();
# @macroexpand @FUmtest

# Global scope testing
# testvector = Vector(0.0:0.01:100.0)
# params = [1.0, 1.0];

# FW1 = FWftest()
# FUM = @FUmtest

# @btime FW1.Ga($testvector, $params) # 9.040 μs (3 allocations: 78.30 KiB)
# @btime FUM.Ga($testvector, $params) # 9.040 μs (3 allocations: 78.30 KiB)

# println("\nDone")

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
using LinearAlgebra

# Base transition matrix
T = [-1 1 0 0; 0 -1 1 0; 0 0 -1 1; 1 0 0 -1]

Λ, V = eigen(T)
Λ, W = eigen(T')

# Right eigenvectors
V⁴ = round.(V[:, 1] ./ V[1, 1])
V³ = round.(V[:, 2] ./ V[1, 2])
V² = round.(V[:, 3] ./ V[1, 3])
V¹ = round.(V[:, 4] ./ V[1, 4])

# Left eigenvectors
W⁴ = round.(W[:, 1] ./ W[1, 1])
W³ = round.(W[:, 2] ./ W[1, 2])
W² = round.(W[:, 3] ./ W[1, 3])
W¹ = round.(W[:, 4] ./ W[1, 4])

# Flip Order (take into account later)
W[:, 1] .= W¹
W[:, 2] .= W²
W[:, 3] .= W³
W[:, 4] .= W⁴
inv(W) * 4

# Quick Check (right eigenvector)
T * V - V * Diagonal(Λ)
maximum(abs.(T * V - V * Diagonal(Λ))) ≤ eps(10.0)
# Quick Check (left eigenvector)
T' * W - W * Diagonal(reverse(Λ))    # since we flipped the order so that index 1 corresponds to the lowest eigenvalue
W' * T - Diagonal(reverse(Λ))' * W'  # same thing, recall that ' means adjoint, thus takes the conjucate
maximum(abs.(T' * W - W * Diagonal(reverse(Λ)))) ≤ eps(10.0)

# The following transformation is what is used
S = real.(V)
S[:, 1] .= real.(W¹)
S[:, 2] .= real.(0.5 * (W² + W³))
S[:, 3] .= real.(0.5 * im * (W² - W³))
S[:, 4] .= real.(W⁴)

display(S)
display(inv(S))
display(inv(S) * T * S)

# The following partitions the magnitude of
# the anti-symmetric component and the symmetric component
# of the transition matrix
sT = (T + T') / 2
aT = (T - T') / 2
λ = 1.0
ω = 1.0
inv(S) * (λ * sT + ω * aT) * S

## check kronecker products corresponding to random phase
kroneckered_phase = kron(sT, I + 0 * sT) + kron(I + 0 * sT, sT)
Λ, V = eigen(kroneckered_phase)

## check out higher dimensional one
idmat = I + 0 * sT # the identity matrix
kroneckered_phase = kron(sT, idmat, idmat) + kron(idmat, sT, idmat) + kron(idmat, idmat, sT)
Λ, V = eigen(kroneckered_phase)

max_stars_repo_name: sandreza/StatisticalNonlocality

max_issues_repo_name: sandreza/StatisticalNonlocality

max_forks_repo_name: sandreza/StatisticalNonlocality


text:
module ThinWallBCTest

using Gridap
using Test
using LineSearches: BackTracking

using GridapMHD: conductive_thin_wall
using GridapMHD: analytical_shercliff_u
using GridapMHD: analytical_shercliff_j
using GridapMHD: compute_u_j_errors

# Problem setting
ρ = 1.0
ν = 1.0
σ = 1.0
wall_thickness = 0.10
wall_conductivity = 1.0
c_w = wall_conductivity * wall_thickness/σ

U0 = 10.0
B0 = 10.0
L = 1.0

Re = U0 * L / ν
Ha = B0 * L * sqrt(σ/(ρ*ν))
K = Ha / (1-0.825*Ha^(-1/2)-Ha^(-1))
∂p∂z = -Re * K / L^3

f_u = VectorValue(0.0,0.0, -∂p∂z) * L/U0^2

xh, trian, dΩ = conductive_thin_wall(nx=5,ny=5,c_w=c_w,Re=Re,Ha=Ha,α=10.0,f_u=f_u)

uh, ph, jh, φh = xh
divj = (∇⋅jh)

# Scale unknowns
uh = uh * U0
jh = jh * σ * B0 * U0

# Analyical solutions
side_wall_semilength = 1.0
hartmann_wall_semilength = 1.0
number_fourier_sumands = 20
K = Ha / (1-0.825*Ha^(-1/2)-Ha^(-1))
∂p∂z = -Re * K / L^3
u0(x) = analytical_shercliff_u(side_wall_semilength, hartmann_wall_semilength,
  wall_thickness, wall_conductivity, σ, ν*ρ, ∂p∂z, Ha, number_fourier_sumands, x)

j0(x) = analytical_shercliff_j(side_wall_semilength, hartmann_wall_semilength,
  wall_thickness, wall_conductivity, σ, ν*ρ, ∂p∂z, Ha, number_fourier_sumands, x)

writevtk(trian, "results.vtu",
  cellfields=["uh"=>uh, "ph"=>ph, "jh"=>jh, "φh"=>φh, "u"=>u0, "j"=>j0])

eu_l2, ej_l2 = compute_u_j_errors(uh, jh, u0, j0, dΩ)
e_divj = sqrt(sum(∫(divj*divj)*dΩ))

@test eu_l2 < 0.3
@test ej_l2 < 2
@test e_divj < 1e-12

end #module


text:
export Expectation, Moment

const ExpFloat = Float32

ExpCacheDict = Dict{Pair{ProbCircuit, LogisticCircuit}, Array{ExpFloat, 2}}
MomentCacheDict = Dict{Tuple{ProbCircuit, LogisticCircuit, Int64}, Array{ExpFloat, 2}}

struct ExpectationCache
    f::ExpCacheDict
    fg::ExpCacheDict
end

ExpectationCache() = ExpectationCache(ExpCacheDict(), ExpCacheDict())

struct MomentCache
    f::ExpCacheDict
    fg::MomentCacheDict
end
MomentCache() = MomentCache( ExpCacheDict(),  MomentCacheDict())

# Find a better way to cache n_choose_k values
max_k = 31
choose_cache = [ 1.0 * binomial(i,j) for i=0:max_k+1, j=0:max_k+1 ]
@inline function choose(n::Int, m::Int)
    return choose_cache[n+1, m+1]
end

"""
    Expectation(pc::ProbCircuit, lc::LogisticCircuit, data)

Compute Expected Prediction of a Logistic/Regression Circuit w.r.t to a ProbabilistcCircuit

Missing values should be denoted by missing
See: On Tractable Computation of Expected Predictions [arxiv.org/abs/1910.02182](https://arxiv.org/abs/1910.02182)
"""
function Expectation(pc::ProbCircuit, lc::LogisticCircuit, data)
    # 1. Get probability of each observation
    log_likelihoods = marginal(pc, data)
    p_observed = exp.( log_likelihoods )

    # 2. Expectation w.r.t. P(x_m, x_o)
    cache = ExpectationCache()
    results_unnormalized = exp_g(pc, children(lc)[1], data, cache) # skipping the bias node of lc

    # 3. Expectation w.r.t P(x_m | x_o)
    results = transpose(results_unnormalized) ./ p_observed

    # 4. Add Bias terms
    biases = lc.thetas
    results .+= biases

    results, cache
end

"""
    Moment(pc::ProbCircuit, lc::LogisticCircuit, data, moment::Int)

Compute higher moments of Expected Prediction for the pair of Logistic/Regression Circuit, ProbabilistcCircuit
"""
function Moment(pc::ProbCircuit, lc::LogisticCircuit, data, moment::Int)
    # 1. Get probability of each observation
    log_likelihoods = marginal(pc, data)
    p_observed = exp.( log_likelihoods )

    # 2. Moment w.r.t. P(x_m, x_o)
    cache = MomentCache()
    biases = lc.thetas
    results_unnormalized = zeros(ExpFloat, num_examples(data), num_classes(lc))

    for z = 0:moment-1
        results_unnormalized .+= choose(moment, z) .* (biases .^ (z)) .* transpose(moment_g(pc, children(lc)[1], data, moment - z, cache))
    end

    # 3. Moment w.r.t P(x_m | x_o)
    results = results_unnormalized ./ p_observed

    # 4. Add Bias^moment terms
    results .+= biases .^ (moment)

    results, cache
end

# exp_f (pr-constraint) is originally from:
#   Arthur Choi, Guy Van den Broeck, and Adnan Darwiche. Tractable learning for structured probability spaces: A case study in learning preference distributions. In Proceedings of IJCAI, 2015.

function exp_f(n::Union{PlainSumNode, StructSumNode}, m::Logistic⋁Node, data, cache::Union{ExpectationCache, MomentCache})
    @inbounds get!(cache.f, Pair(n, m)) do
        value = zeros(ExpFloat, 1 , num_examples(data) )
        @inline pthetas(i) = ExpFloat.(exp(n.log_probs[i])) #[exp(n.log_probs[i]) for i in 1:num_children(n)]
        @fastmath @simd for i in 1:num_children(n)
            @simd for j in 1:num_children(m)
                value .+= (pthetas(i) .* exp_f(children(n)[i], children(m)[j], data, cache))
            end
        end
        return value
    end
end

function exp_f(n::Union{PlainMulNode, StructMulNode}, m::Logistic⋀Node, data, cache::Union{ExpectationCache, MomentCache})
    @inbounds get!(cache.f, Pair(n, m)) do
        # value = ones(ExpFloat, 1 , num_examples(data) )
        # @fastmath for (i,j) in zip(children(n), children(m))
        #     value .*= exp_f(i, j, data, cache)
        # end
        # return value
        @fastmath @inbounds exp_f(children(n)[1], children(m)[1], data, cache) .* exp_f(children(n)[2], children(m)[2], data, cache)
    end
    # @fastmath @inbounds exp_f(children(n)[1], children(m)[1], data, cache) .* exp_f(children(n)[2], children(m)[2], data, cache)
end

@inline function exp_f(n::Union{PlainProbLiteralNode, StructProbLiteralNode}, m::LogisticLiteralNode, data, cache::Union{ExpectationCache, MomentCache})
    @inbounds get!(cache.f, Pair(n, m)) do
        value = zeros(ExpFloat, 1 , num_examples(data) )
        var = lit2var(literal(m))
        X = data
        if ispositive(n) && ispositive(m)
            # value[1, X[:, var] .== -1 ] .= 1.0  # missing observation always agrees
            # value[1, X[:, var] .== 1 ] .= 1.0 # positive observations
            value[1, .!isequal.(X[:, var], 0)] .= ExpFloat(1.0) # positive or missing observations
        elseif isnegative(n) && isnegative(m)
            # value[1, X[:, var] .== -1 ] .= 1.0  # missing observation always agrees
            # value[1, X[:, var] .== 0 ] .= 1.0 # negative observations
            value[1, .!isequal.(X[:, var], 1)] .= ExpFloat(1.0) # negative or missing observations
        end
        return value
    end
end

"""
Has to be a Logistic⋁Node with only one child, which is a leaf node
"""
@inline function exp_f(n::Union{PlainProbLiteralNode, StructProbLiteralNode}, m::Logistic⋁Node, data, cache::Union{ExpectationCache, MomentCache})
    @inbounds get!(cache.f, Pair(n, m)) do
        exp_f(n, children(m)[1], data, cache)
    end
end

#######################################################################
######## exp_g, exp_fg
########################################################################

@inline function exp_g(n::Union{PlainSumNode, StructSumNode}, m::Logistic⋁Node, data, cache::ExpectationCache)
    exp_fg(n, m, data, cache) # exp_fg and exp_g are the same for OR nodes
end

# function exp_g(n::Prob⋀, m::Logistic⋀Node, data, cache::ExpectationCache)
#     value = zeros(ExpFloat, classes(m) , num_examples(data))
#     @fastmath for (i,j) in zip(children(n), children(m))
#         value .+= exp_fg(i, j, data, cache)
#     end
#     return value
#     # exp_fg(children(n)[1], children(m)[1], data, cache) .+ exp_fg(children(n)[2], children(m)[2], data, cache)
# end

function exp_fg(n::Union{PlainSumNode, StructSumNode}, m::Logistic⋁Node, data, cache::ExpectationCache)
    @inbounds get!(cache.fg, Pair(n, m)) do
        value = zeros(ExpFloat, num_classes(m) , num_examples(data) )
        @inline pthetas(i) = ExpFloat.(exp(n.log_probs[i])) #pthetas = [exp(n.log_probs[i]) for i in 1:num_children(n)]
        @fastmath @simd for i in 1:num_children(n)
            for j in 1:num_children(m)
                value .+= (pthetas(i) .* ExpFloat.(m.thetas[j,:])) .* exp_f(children(n)[i], children(m)[j], data, cache)
                value .+= pthetas(i) .* exp_fg(children(n)[i], children(m)[j], data, cache)
            end
        end
        return value
    end
end

function exp_fg(n::Union{PlainMulNode, StructMulNode}, m::Logistic⋀Node, data, cache::ExpectationCache)
    @inbounds get!(cache.fg, Pair(n, m)) do
        # Assuming 2 children
        value = exp_f(children(n)[1], children(m)[1], data, cache) .* exp_fg(children(n)[2], children(m)[2], data, cache)
        value .+= exp_f(children(n)[2], children(m)[2], data, cache) .* exp_fg(children(n)[1], children(m)[1], data, cache)
        return value
    end
    # @fastmath @inbounds exp_f(children(n)[1], children(m)[1], data, cache) .* exp_fg(children(n)[2], children(m)[2], data, cache)
    #     .+ exp_f(children(n)[2], children(m)[2], data, cache) .* exp_fg(children(n)[1], children(m)[1], data, cache)
end

"""
Has to be a Logistic⋁Node with only one child, which is a leaf node
"""
@inline function exp_fg(n::Union{PlainProbLiteralNode, StructProbLiteralNode}, m::Logistic⋁Node, data, cache::ExpectationCache)
    @inbounds get!(cache.fg, Pair(n, m)) do
        ExpFloat.(m.thetas[1,:]) .* exp_f(n, m, data, cache)
    end
end

@inline function exp_fg(n::Union{PlainProbLiteralNode, StructProbLiteralNode}, m::LogisticLiteralNode, data, cache::ExpectationCache)
    #dont know how many classes, boradcasting does the job
    zeros(ExpFloat, 1 , num_examples(data))
end

#######################################################################
######## moment_g, moment_fg
########################################################################

@inline function moment_g(n::Union{PlainSumNode, StructSumNode}, m::Logistic⋁Node, data, moment::Int, cache::MomentCache)
    get!(cache.fg, (n, m, moment)) do
        moment_fg(n, m, data, moment, cache)
    end
end

"""
Calculating  E[g^k * f]
"""
function moment_fg(n::Union{PlainSumNode, StructSumNode}, m::Logistic⋁Node, data, moment::Int, cache::MomentCache)
    if moment == 0
        return exp_f(n, m, data, cache)
    end

    get!(cache.fg, (n, m, moment)) do
        value = zeros(ExpFloat, num_classes(m) , num_examples(data) )
        pthetas = [exp(n.log_probs[i]) for i in 1:num_children(n)]
        @fastmath @simd for i in 1:num_children(n)
            for j in 1:num_children(m)
                for z in 0:moment
                    value .+= pthetas[i] .* choose(moment, z) .* m.thetas[j,:].^(moment - z) .* moment_fg(children(n)[i], children(m)[j], data, z, cache)
                end
            end
        end
        return value
    end
end

@inline function moment_fg(n::Union{PlainProbLiteralNode, StructProbLiteralNode}, m::Logistic⋁Node, data, moment::Int, cache::MomentCache)
    get!(cache.fg, (n, m, moment)) do
        m.thetas[1,:].^(moment) .* exp_f(n, m, data, cache)
    end
end

@inline function moment_fg(n::Union{PlainProbLiteralNode, StructProbLiteralNode}, m::LogisticLiteralNode, data, moment::Int, cache::MomentCache)
    #dont know how many classes, boradcasting does the job
    if moment == 0
        exp_f(n, m, data, cache)
    else
        zeros(ExpFloat, 1, num_examples(data))
    end
end

function moment_fg(n::Union{PlainMulNode, StructMulNode}, m::Logistic⋀Node, data, moment::Int, cache::MomentCache)
    if moment == 0
        return exp_f(n, m, data, cache)
    end
    get!(cache.fg, (n, m, moment)) do
        value = moment_fg(children(n)[1], children(m)[1], data, 0, cache) .* moment_fg(children(n)[2], children(m)[2], data, moment, cache)

        for z in 1:moment
            value .+= choose(moment, z) .* moment_fg(children(n)[1], children(m)[1], data, z, cache) .* moment_fg(children(n)[2], children(m)[2], data, moment - z, cache)
        end
        return value
    end
end

max_stars_repo_name: Juice-jl/DiscriminativeCircuits

max_issues_repo_name: Juice-jl/DiscriminativeCircuits

max_forks_repo_name: Juice-jl/DiscriminativeCircuits


text:
mutable struct SpacePartitioningTree{
    T<:AbstractFloat,
    I<:Integer}

    cuts::I
    leafsize::I

    dimensionlist::Vector{I}
    recursiondepth::I
    cutlist::Vector{T}
end
SpacePartitioningTree(T::DataType, I::DataType) = SpacePartitioningTree{T, I}(zero(I), zero(I), zeros(I, 0), zero(I), zeros(I, 0))
isinitialized(x::SpacePartitioningTree) = !(iszero(x.cuts) && iszero(x.leafsize) && isempty(x.dimensionlist) && iszero(x.recursiondepth) && isempty(x.cutlist))

"""
    DataSet{T<:AbstractFloat, I<:Integer}

*AHMI-internal, not part of stable public API.*

Holds the MCMC output. For construction use constructor: function DataSet{T<:Real}(data::Matrix{T}, logprob::Vector{T}, weights::Vector{T})
# Variables
- 'data' : An P x N array with N data points with P parameters.
- 'logprob' : The logarithmic probability for each samples stored in an array
- 'weights' : How often each sample occurred. Set to an array of ones if working directly on MCMC output
- 'ids' : Array which is used to assign each sample to a batch, required for the cov. weighed uncertainty estimation
- .sortids : an array of indices which stores the original ordering of the samples (the space partitioning tree reorders the samples), required to calculate an effective sample size.
- 'N' : number of samples
- 'P' : number of parameters
- 'nsubsets' : the number of batches
- 'iswhitened' : a boolean value which indicates whether the data set is iswhitened
- 'isnew' : a boolean value which indicates whether the data set was swapped out with a different one (it is possible to redo the integration with a different sample set using previously generated hyper-rectangles)
- 'partitioningtree' : The space partitioning tree, used to efficiently identify samples in a point cloud
- 'startingIDs' : The Hyper-Rectangle Seed Samples are stored in this array
- 'tolerance' : A threshold required for the hyper-rectangle creation process.
"""
mutable struct DataSet{T<:AbstractFloat, I<:Integer}
    data::Array{T, 2}
    logprob::Array{T, 1}
    weights::Array{T, 1}
    ids::Array{I, 1}    #used to divide the dataset into sub-sets
    sortids::Array{I, 1}#used to calculate the ess on the unsorted dataset
    N::I
    P::I
    nsubsets::I    #number of sub-sets
    subsetsize::T
    iswhitened::Bool
    isnew::Bool
    partitioningtree::SpacePartitioningTree
    startingIDs::Array{I, 1}
    tolerance::T
end

function DataSet(
    data:: Array{T, 2},
    logprob::Array{T, 1},
    weights::Array{I, 1},
    nsubsets::Integer = 0,
    subsetsize::T = zero(T)
    )::DataSet{T, Int} where {T<:AbstractFloat, I<:Integer}

    DataSet(data, logprob, convert(Array{T, 1}, weights), nsubsets, subsetsize)
end

function DataSet(
    data:: Array{T, 2},
    logprob::Array{T, 1},
    weights::Array{T, 1},
    nsubsets::Integer = 0,
    subsetsize::T = zero(T)
    )::DataSet{T, Int} where {T<:AbstractFloat}

    P, N = size(data)

    if iszero(nsubsets)
        nsubsets = 10
    end

    maxbatchsize = sum(weights) / 10 / nsubsets
    if iszero(subsetsize)
        subsetsize = 100.0
    end
    subsetsize = minimum(typeof(maxbatchsize), (maxbatchsize, subsetsize))

    ids = zeros(Int, N)
    cnt = 1

    batch_currentsize = 0.0

    for i = 1:N
        ids[i] = cnt
        batch_currentsize += weights[i]

        if batch_currentsize >= subsetsize
            cnt += 1
            if cnt > nsubsets
                cnt = 1
            end
            batch_currentsize = 0.0
        end
    end
    DataSet(data, logprob, weights, ids, [i for i=1:N], N, P, nsubsets, subsetsize, false, true, SpacePartitioningTree(T, Int), zeros(Int, 0), T(0))
end

Base.show(io::IO, data::DataSet) = print(io, "DataSet: $(data.N) samples, $(data.P) parameters")
Base.eltype(data::DataSet{T, I}) where {T<:AbstractFloat, I<:Integer} = (T, I)

"""
    HMISettings

*AHMI-internal, not part of stable public API.*

holds the settings for the hm_integrate function. There are several default constructors available:
HMIFastSettings()
HMIStandardSettings()
HMIPrecisionSettings()

#Variables
- 'whitening_method::Symbol' : which whitening method to use
- 'max_startingIDs::Integer' : influences how many starting ids are allowed to be generated
- 'max_startingIDs_fraction::AbstractFloat' : how many points are considered as possible starting points as a fraction of total points available
- 'rect_increase::AbstractFloat' : describes the procentual rectangle volume increase/decrease during hyperrectangle creation. Low values can increase the precision if enough points are available but can cause systematically wrong results if not enough points are available.
- 'use_all_rects::Bool' : All rectangles are used for the integration process no matter how big their overlap is. If enabled the rectangles are weighted by their overlap.
- 'useMultiThreading' : activate multithreading support.
- 'warning_minstartingids' : the required minimum amount of starting samples
- 'dotrimming' : determines whether the integral estimates are trimmed (1σ trim) before combining them into a final result (more robust)
- 'uncertainty_estimators' : A dictionary of different uncertainty estimator functions. Currently three functions are available: hm_combineresults_legacy! (outdated, overestimates uncertainty significantly in higher dimensions), hm_combineresults_covweighted! (very fast) and hm_combineresults_analyticestimation! (recommended)
end
"""
mutable struct HMISettings
    whitening_function!::Function
    max_startingIDs::Integer
    max_startingIDs_fraction::AbstractFloat
    rect_increase::AbstractFloat
    useMultiThreading::Bool
    warning_minstartingids::Integer
    dotrimming::Bool
    uncertainty_estimators::Dict{String, Function}
end
HMIFastSettings() =      return HMISettings(cholesky_whitening!, 100,   0.1, 0.1, true, 16, true, Dict("cov. weighted result" => hm_combineresults_covweighted!))
HMIStandardSettings() =  return HMISettings(cholesky_whitening!, 1000,  0.5, 0.1, true, 16, true, Dict("cov. weighted result" => hm_combineresults_covweighted!, "analytic result" => hm_combineresults_analyticestimation!))
HMIPrecisionSettings() = return HMISettings(cholesky_partial_whitening!, 10000, 2.5, 0.1, true, 16, true, Dict("cov. weighted result" => hm_combineresults_covweighted!))
# HMIPrecisionSettings() = return HMISettings(cholesky_whitening!, 10000, 2.5, 0.1, true, 16, true, Dict("cov. weighted result" => hm_combineresults_covweighted!, "analytic result" => hm_combineresults_analyticestimation!))

"""
    WhiteningResult{T<:AbstractFloat}

*AHMI-internal, not part of stable public API.*

Stores the information obtained during the Whitening Process
# Variables
- 'determinant' : The determinant of the whitening matrix
- 'targetprobfactor' : The suggested target probability factor
- 'whiteningmatrix' : The whitening matrix
- 'meanvalue' : the mean vector of the input data
"""
struct WhiteningResult{T<:AbstractFloat}
    determinant::T
    targetprobfactor::T
    whiteningmatrix::Matrix{T}
    meanvalue::Vector{T}
end
WhiteningResult(T::DataType) = WhiteningResult(zero(T), zero(T), zeros(T, 0, 0), zeros(T, 0))
Base.show(io::IO, wres::WhiteningResult) = print(io, "Whitening Result: Determinant: $(wres.determinant), Target Prob. Factor: $(wres.targetprobfactor)")
isinitialized(x::WhiteningResult) = !(iszero(x.determinant) && iszero(x.targetprobfactor) && isempty(x.whiteningmatrix) && isempty(x.meanvalue))

"""
    SearchResult{T<:AbstractFloat, I<:Integer}

*AHMI-internal, not part of stable public API.*

Stores the results of the space partitioning tree's search function

# Variables
- 'pointIDs' : the IDs of samples found, might be empty because it is optional
- 'points' : The number of points found.
- 'maxLogProb' : the maximum log. probability of the points found.
- 'minLogProb' : the minimum log. probability of the points found.
- 'maxWeightProb' : the weighted minimum log. probability found.
- 'minWeightProb' : the weighted maximum log. probfactor found.
"""
mutable struct SearchResult{T<:AbstractFloat, I<:Integer}
    pointIDs::Vector{I}
    points::I
    maxLogProb::T
    minLogProb::T
    maxWeightProb::T
    minWeightProb::T
end

function SearchResult(T::DataType, I::DataType)
    @assert T<:AbstractFloat
    @assert I<:Integer
    return SearchResult{T, I}(zeros(I, 0), I(0), T(0), T(0), T(0), T(0))
end
Base.show(io::IO, sres::SearchResult) = print(io, "Search Result: Points: $(sres.points), Max. Log. Prob.: $(sres.maxLogProb), Min. Log. Prob.: $(sres.minLogProb)")

"""
    PointCloud{T<:AbstractFloat, I<:Integer}

*AHMI-internal, not part of stable public API.*

Stores the information of the points of an e.g. HyperRectVolume
# Variables
- 'maxLogProb' : The maximum log. probability of one of the points inside the hyper-rectangle
- 'minLogProb' : The minimum log. probability of one of the points inside the hyper-rectangle
- 'maxWeightProb' : the weighted max. log. probability
- 'minWeightProb' : the weighted min. log. probability
- 'probfactor' : The probability factor of the hyper-rectangle
- 'probweightfactor' : The weighted probability factor
- 'points' : The number of points inside the hyper-rectangle
- 'pointIDs' : the IDs of the points inside the hyper-rectangle, might be empty because it is optional and costs performance
- 'searchres' : used to boost performance
"""
mutable struct PointCloud{T<:AbstractFloat, I<:Integer}
    maxLogProb::T
    minLogProb::T
    maxWeightProb::T
    minWeightProb::T
    probfactor::T
    probweightfactor::T
    points::I
    pointIDs::Vector{I}
    searchres::SearchResult{T, I}
end
PointCloud(T::DataType, I::DataType) = PointCloud(zero(T), zero(T), zero(T), zero(T), zero(T), zero(T), zero(I), zeros(I, 0), SearchResult(T, I))
Base.show(io::IO, cloud::PointCloud) = print(io, "Point Cloud with $(cloud.points) points, probability factor: $(cloud.probfactor)")

"""
    IntegrationVolume{T<:AbstractFloat, I<:Integer}

*AHMI-internal, not part of stable public API.*

# Variables
- 'pointcloud' : holds the point cloud of the integration volume
- 'spatialvolume' : the boundaries of the integration volume
- 'volume' : the volume

Hold the point cloud and the spatial volume for integration.
"""
mutable struct IntegrationVolume{T<:AbstractFloat, I<:Integer, V<:SpatialVolume}
    pointcloud::PointCloud{T, I}
    spatialvolume::V
    volume::T
end
Base.show(io::IO, vol::IntegrationVolume) = print(io, "Hyperrectangle: $(vol.pointcloud.points) points, $(vol.volume) Volume")

mutable struct IntermediateResults{T<:AbstractFloat}
    integrals::Array{T, 1}
    volumeID::Array{Int, 1}
    Y::Array{T, 2}
end
IntermediateResults(T::DataType, n::Integer) = IntermediateResults(zeros(T, n), [Int(i) for i=1:n], zeros(T, 0, 0))
Base.length(x::IntermediateResults) = length(x.integrals)

mutable struct HMIEstimate{T<:AbstractFloat}
    estimate::T
    uncertainty::T
    weights::Array{T, 1}
end
HMIEstimate(T::DataType) = HMIEstimate(zero(T), zero(T), zeros(T, 0))
function HMIEstimate(a::HMIEstimate{T}, b::HMIEstimate{T})::HMIEstimate{T} where {T<:AbstractFloat}
    val = mean([a.estimate, b.estimate], AnalyticWeights([1 / a.uncertainty^2, 1 / b.uncertainty^2]))
    unc = 1 / sqrt(1 / a.uncertainty^2 + 1 / b.uncertainty^2)
    HMIEstimate(val, unc, [a.weights..., b.weights...])
end
Base.show(io::IO, ires::HMIEstimate) = print(io, "$(round(ires.estimate, sigdigits=6))  +-  $(round(ires.uncertainty, sigdigits=6))")

mutable struct HMIResult{T<:AbstractFloat}
    result1::HMIEstimate{T}
    result2::HMIEstimate{T}
    final::HMIEstimate{T}
    dat1::Dict{String, Any}
    dat2::Dict{String, Any}
end
HMIResult(T::DataType) = HMIResult(HMIEstimate(T), HMIEstimate(T), HMIEstimate(T), Dict{String, Any}(), Dict{String, Any}())

"""
    HMIData{T<:AbstractFloat, I<:Integer}

*AHMI-internal, not part of stable public API.*

Includes all the informations of the integration process, including a list of hyper-rectangles, the results of the whitening transformation,
the starting ids, and the average number of points and volume of the created hyper-rectangles.

# Variables
- 'dataset1' : Data Set 1
- 'dataset2' : Data Set 2
- 'whiteningresult' : contains the whitening matrix and its determinant, required to scale the final integral estimate
- 'volumelist1' : An array of integration volumes created using dataset1, but filled with samples from dataset2
- 'volumelist2' : An array of integration volumes created using dataset2, but filled with samples from dataset1
- 'cubelist1' : An array of small hyper-cubes created around seeding samples of dataset 1
- 'cubelist2' : An array of small hyper-cubes created around seeding samples of dataset 2
- 'iterations1' : The number of volume adapting iterations for the creating volumelist1
- 'iterations2' : The number of volume adapting iterations for the creating volumelist2
- 'rejectedrects1' : An array of ids, indicating which hyper-rectangles of volumelist1 were rejected due to trimming
- 'rejectedrects2' : An array of ids, indicating which hyper-rectangles of volumelist2 were rejected due to trimming
- 'integralestimates' : A dictionary containing the final integral estimates with uncertainty estimation using different uncertainty estimators. Also includes all intermediate results required for the integral estimate combination
"""
mutable struct HMIData{T<:AbstractFloat, I<:Integer, V<:SpatialVolume}
    dataset1::DataSet{T, I}
    dataset2::DataSet{T, I}
    whiteningresult::WhiteningResult{T}
    volumelist1::Vector{IntegrationVolume{T, I, V}}
    volumelist2::Vector{IntegrationVolume{T, I, V}}
    cubelist1::Vector{V}
    cubelist2::Vector{V}
    iterations1::Vector{I}
    iterations2::Vector{I}
    rejectedrects1::Vector{I}
    rejectedrects2::Vector{I}
    integrals1::IntermediateResults{T}
    integrals2::IntermediateResults{T}
    integralestimates::Dict{String, HMIResult}
end

function HMIData(
    dataset1::DataSet{T, I},
    dataset2::DataSet{T, I},
    dataType::DataType = HyperRectVolume{T})::HMIData where {T<:AbstractFloat, I<:Integer}

    HMIData(
        dataset1,
        dataset2,
        WhiteningResult(T),
        Vector{IntegrationVolume{T, I, dataType}}(undef, 0),
        Vector{IntegrationVolume{T, I, dataType}}(undef, 0),
        Vector{dataType}(undef, 0),
        Vector{dataType}(undef, 0),
        zeros(I, 0),
        zeros(I, 0),
        zeros(I, 0),
        zeros(I, 0),
        IntermediateResults(T, 0),
        IntermediateResults(T, 0),
        Dict{String, HMIResult}()
    )
end

function HMIData(dataset::DataSet{T, I})::HMIData{T, I} where {T<:AbstractFloat, I<:Integer}
    HMIData(split_dataset(dataset)...)
end

@deprecate HMIData(bat_samples::Tuple{DensitySampleVector, MCMCBasicStats, AbstractVector{<:MCMCIterator}}) HMIData(bat_samples[1])

function HMIData(samples::DensitySampleVector)
    logprob = samples.logd
    weights = samples.weight
    samples = convert(Array{eltype(logprob), 2}, flatview(samples.v))
    ds = DataSet(samples, logprob, weights)
    HMIData(ds)
end

function Base.show(io::IO, ires::HMIData)
    output = "Parameters: $(ires.dataset1.P)\tTotal Samples: $(ires.dataset1.N + ires.dataset2.N)"
    output *= "\nData Set 1: $(length(ires.volumelist1)) Volumes\nData Set 2: $(length(ires.volumelist2)) Volumes"

    for pair in ires.integralestimates
        output *= "\n\nIntegral Estimate ($(pair[1])):\n\t $(pair[2].final)"
    end

    println(io, output)
end


text:
############################################################################################
# Basic Functionality
_modelProb = ModelWrapper(ProbModel(), val_dist)
@testset "Models - basic functionality" begin
    ## Type Check 1 - Constrain/Unconstrain
    theta_unconstrained_vec = randn(length(_modelProb))
    theta_unconstrained = unflatten(_modelProb, theta_unconstrained_vec)
    @test typeof(theta_unconstrained) == typeof(_modelProb.val)
    theta_constrained = constrain(_modelProb.info.transform, theta_unconstrained)
    theta_constrained2 = unflatten_constrain(_modelProb, theta_unconstrained_vec)
    @test typeof(theta_constrained) == typeof(_modelProb.val)
    @test typeof(theta_constrained2) == typeof(_modelProb.val)
    ## Type Check 2 - Flatten/Unflatten
    _θ1, _ = flatten(_modelProb.info.reconstruct, theta_constrained)
    _θ2, _ = flatten(_modelProb.info.reconstruct, theta_constrained2)
    @test sum(abs.(_θ1 - _θ2)) ≈ 0 atol = _TOL
    ## Check if densities match
    @test log_prior(_modelProb) + log_abs_det_jac(_modelProb) ≈
          log_prior_with_transform(_modelProb)
    ## Utility functions
    unconstrain(_modelProb)
    flatten(_modelProb)
    unconstrain_flatten(_modelProb)
    simulate(_modelProb)
    fill(_modelProb, _modelProb.val)
    fill!(_modelProb, _modelProb.val)
    subset(_modelProb, keys(_modelProb.val))
    unflatten!(_modelProb, flatten(_modelProb))
    unflatten_constrain!(_modelProb, unconstrain_flatten(_modelProb))
    sample(_RNG, _modelProb)
    sample(_modelProb)
    sample!(_RNG, _modelProb)
    sample!(_modelProb)
end

############################################################################################
## Model with transforms in lower dimensions

_modelExample = ModelWrapper(ExampleModel(), _val_examplemodel)
_tagged = Tagged(_modelExample)
@testset "Models - Model with transforms in lower dimensions" begin
    ## Model Length accounting discrete parameter
    unconstrain(_modelExample)
    flatten(_modelExample)
    unconstrain_flatten(_modelExample)
    ## Type Check 1 - Constrain/Unconstrain
    theta_unconstrained_vec = randn(length(_modelExample))
    theta_unconstrained = unflatten(_modelExample, theta_unconstrained_vec)
    @test typeof(theta_unconstrained) == typeof(_modelExample.val)
    theta_constrained = constrain(_modelExample.info.transform, theta_unconstrained)
    theta_constrained2 = unflatten_constrain(_modelExample, theta_unconstrained_vec)
    @test typeof(theta_constrained) == typeof(_modelExample.val)
    @test typeof(theta_constrained2) == typeof(_modelExample.val)
    ## Type Check 2 - Flatten/Unflatten
    _θ1, _ = flatten(_modelExample.info.reconstruct, theta_constrained)
    _θ2, _ = flatten(_modelExample.info.reconstruct, theta_constrained2)
    @test sum(abs.(_θ1 - _θ2)) ≈ 0 atol = _TOL
    ## Check if densities match
    @test log_prior(_modelExample) + log_abs_det_jac(_modelExample) ≈
          log_prior_with_transform(_modelExample)
    ## Check utility functions
    @test length(_modelExample) == 23
    @test ModelWrappers.paramnames(_modelExample) == keys(_val_examplemodel)
    fill(_modelExample, _tagged, _modelExample.val)
    fill!(_modelExample, _tagged, _modelExample.val)
end

############################################################################################


text:
# This file was generated, do not modify it. # hide
pprint(sort!(model_info,[:CrossEntropyLoss]));

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
#addprocs(2)
#@show nprocs()

using RoME
using Test
using TensorCast
import Manifolds
using Manifolds: ProductManifold, SpecialEuclidean, ProductRepr, SpecialOrthogonal, TranslationGroup, identity_element
using DistributedFactorGraphs
using Statistics

@error("must restore testG2oParser.jl")
@error("must restore testParametric.jl")
@error("add test for generateCanonicalFG_Beehive, norm( simulated - default ) < tol")

testfiles = [
  # any wip fail-fast dev testing
  "testPartialPose3.jl";

  # important tests are that are broken and must be restored.
  "testBearingRange2D.jl";
  "testBearing2D.jl";
  "testMultimodalRangeBearing.jl"; # restore after Bearing factors are fixed

  # "testParametric.jl"; # deferred to v0.16.x
  # "testG2oParser.jl";  # deferred to v0.16.x

  # tests most likely to fail on numerics
  "testScalarFields.jl";
  "testPoint2Point2Init.jl";
  "threeDimLinearProductTest.jl";
  "testBeehiveGrow.jl"; # also starts multiprocess

  # regular tests expected to pass
  "testpackingconverters.jl";
  "testInflation380.jl";
  "testPoint2Point2.jl";
  "testParametricSimulated.jl";
  "testBasicPose2Conv.jl";
  "testGraphGenerators.jl";
  "testTreeInitCommonMsg_IIF913.jl";
  "testPose3Pose3NH.jl";
  "testHexagonal2D_CliqByCliq.jl";      # special case debugging
  "testhigherdimroots.jl";
  "testDidsonFunctions.jl";
  "testBasicPose2Stationary.jl";
  "TestPoseAndPoint2Constraints.jl";
  "testPartialRangeCrossCorrelations.jl";
  "testDynPoint2D.jl";
  "testDeltaOdo.jl";
  "testFixedLagFG.jl";
  "testDynPose2D.jl";
  "testPartialPriorYawPose2.jl";
  "TestDefaultFGInitialization.jl";
  "testAccumulateFactors.jl";
  "testDeadReckoningTether.jl";
  "testFluxModelsPose2.jl";
  "testGenerateHelix.jl";
]

## Tests not ready yet
# "HexagonalLightGraphs.jl"
# "testCameraFunctions.jl"
# "testmultiplefeatures.jl"

for testf in testfiles
  println("[TEST] $testf")
  include(testf)
  println("[SUCCESS] $testf ==========================================================")
  println()
  println()
  println()
end

#


text:
using Flux
using Flux: onehotbatch, onecold, flatten
using Flux.Losses: logitcrossentropy
using Flux.Data: DataLoader
using Parameters: @with_kw
using Statistics: mean
using CUDA
using ProgressMeter: @showprogress
using Formatting
using MLDatasets: CIFAR10, SVHN2

using DeepUncertainty

if CUDA.has_cuda()
    @info "CUDA is on"
    CUDA.allowscalar(false)
end

@with_kw mutable struct Args
    batchsize::Int = 64
    lr::Float64 = 3e-4
    epochs::Int = 50
    valsplit::Float64 = 0.1
    sample_size = 10
    complexity_constant = 1e-8
end

function accuracy(preds, labels)
    acc = sum(onecold(preds |> cpu) .== onecold(labels |> cpu))
    return acc
end

loss(ŷ, y) = logitcrossentropy(ŷ, y)

num_params(model) = sum(length, Flux.params(model))

round4(x) = round(x, digits = 4)

function get_data(args)
    x, y = CIFAR10.traindata()
    train_x = float(x)
    train_y = onehotbatch(y, 0:9)
    train_data = (train_x, train_y)

    test_x, test_y = CIFAR10.testdata()
    test_x = float(test_x)
    test_y = onehotbatch(test_y, 0:9)
    test_data = (test_x, test_y)

    # # OOD dataset
    # ood_test_x, ood_test_y = SVHN2.testdata()
    # ood_test_x = float(ood_test_x)
    # ood_test_y = onehotbatch(ood_test_y, 1:10)

    train_loader =
        DataLoader(train_data, batchsize = args.batchsize, shuffle = true, partial = false)
    test_loader = DataLoader(test_data, batchsize = args.batchsize)
    #ood_test_loader = DataLoader(ood_test_data, batchsize=args.batchsize)

    return train_loader, test_loader
end

# VGG16 and VGG19 models
function vgg16()
    dropout = 0.2
    Chain(
        MCConv((3, 3), 3 => 64, dropout, relu, pad = (1, 1), stride = (1, 1)),
        BatchNorm(64),
        MCConv((3, 3), 64 => 64, dropout, relu, pad = (1, 1), stride = (1, 1)),
        BatchNorm(64),
        MaxPool((2, 2)),
        MCConv((3, 3), 64 => 128, dropout, relu, pad = (1, 1), stride = (1, 1)),
        BatchNorm(128),
        MCConv((3, 3), 128 => 128, dropout, relu, pad = (1, 1), stride = (1, 1)),
        BatchNorm(128),
        MaxPool((2, 2)),
        MCConv((3, 3), 128 => 256, dropout, relu, pad = (1, 1), stride = (1, 1)),
        BatchNorm(256),
        MCConv((3, 3), 256 => 256, dropout, relu, pad = (1, 1), stride = (1, 1)),
        BatchNorm(256),
        MCConv((3, 3), 256 => 256, dropout, relu, pad = (1, 1), stride = (1, 1)),
        BatchNorm(256),
        MaxPool((2, 2)),
        MCConv((3, 3), 256 => 512, dropout, relu, pad = (1, 1), stride = (1, 1)),
        BatchNorm(512),
        MCConv((3, 3), 512 => 512, dropout, relu, pad = (1, 1), stride = (1, 1)),
        BatchNorm(512),
        MCConv((3, 3), 512 => 512, dropout, relu, pad = (1, 1), stride = (1, 1)),
        BatchNorm(512),
        MaxPool((2, 2)),
        MCConv((3, 3), 512 => 512, dropout, relu, pad = (1, 1), stride = (1, 1)),
        BatchNorm(512),
        MCConv((3, 3), 512 => 512, dropout, relu, pad = (1, 1), stride = (1, 1)),
        BatchNorm(512),
        MCConv((3, 3), 512 => 512, dropout, relu, pad = (1, 1), stride = (1, 1)),
        BatchNorm(512),
        MaxPool((2, 2)),
        flatten,
        Dense(512, 4096, relu),
        Dropout(0.5),
        Dense(4096, 4096, relu),
        Dropout(0.5),
        Dense(4096, 10),
    )
end

function test(args, loader, model)
    loss = 0
    acc = 0
    ece = 0
    entropy = 0
    ntot = 0
    for (x, y) in loader
        predictions = []

        x, y = x |> gpu, y |> gpu
        # Loop through each model's predictions
        for ensemble = 1:args.sample_size
            logits = model(x)
            push!(predictions, logits)
        end

        # Get the mean predictions
        mean_preds = Flux.batch(predictions)
        mean_predictions = mean(mean_preds, dims = ndims(mean_preds))
        logits = dropdims(mean_predictions, dims = ndims(mean_predictions))

        n = size(logits)[end]
        loss += logitcrossentropy(logits, y) * n
        acc += accuracy(logits, y)
        ece += expected_calibration_error(logits, onecold(y)) * n
        entropy += mean(calculate_entropy(logits)) * n
        ntot += n
    end

    mean_loss = loss / ntot |> round4
    mean_acc = acc / ntot * 100 |> round4
    mean_ece = ece / ntot |> round4
    mean_entropy = entropy / ntot |> round4

    @info (format(
        "Mean Loss: {} Mean Accuracy: {} Mean ECE: {} Mean Entropy: {}",
        mean_loss,
        mean_acc,
        mean_ece,
        mean_entropy,
    ))
    @info "==========================================================="
    return nothing
end

function train(; kws...)
    # Initialize the hyperparameters
    args = Args(; kws...)

    # Load the train, validation data
    train_loader, test_loader = get_data(args)

    @info("Constructing Model")
    m = vgg16() |> gpu

    ## Training
    # Defining the optimizer
    opt = ADAM(args.lr)
    ps = Flux.params(m)

    test(args, test_loader, m)

    @info("Training....")
    # Starting to train models
    for epoch = 1:args.epochs
        @info "Epoch $epoch"

        loss_fn(x, y) = logitcrossentropy(m(x), y)

        @showprogress for (x, y) in train_loader
            x, y = x |> gpu, y |> gpu
            gs = Flux.gradient(ps) do
                loss_fn(x, y)
            end
            Flux.update!(opt, ps, gs)
        end
        test(args, test_loader, m)
    end
    test(args, test_loader, m)

    return m
end

m = train()


text:
struct SingularException <: Exception end

struct DggesException <: Exception
    error_nbr::Int64
end

Base.showerror(io::IO, e::DggesException) = print(io, "dgges error ", e.error_nbr)


text:
abstract type AbstractOptimizationSolution{T, N} <: AbstractNoTimeSolution{T, N} end

struct OptimizationSolution{T, N, uType, P, A, Tf, O} <: AbstractOptimizationSolution{T, N}
    u::uType # minimizer
    prob::P # optimization problem
    alg::A # algorithm
    minimum::Tf
    retcode::Symbol
    original::O # original output of the optimizer
end

function build_solution(prob::AbstractOptimizationProblem,
                        alg, u, minimum;
                        retcode = :Default,
                        original = nothing,
                        kwargs...)

    T = eltype(eltype(u))
    N = ndims(u)

    OptimizationSolution{T, N, typeof(u), typeof(prob), typeof(alg),
                         typeof(minimum), typeof(original)}(
                         u, prob, alg, minimum, retcode, original)
end

function Base.show(io::IO, A::AbstractOptimizationSolution)
    println(io,string("retcode: ",A.retcode))
    print(io,"u: ")
    show(io, A.u)
    println(io)
    print(io,"Final objective value:     $(A.minimum)\n")
    return
end

Base.@propagate_inbounds function Base.getproperty(x::AbstractOptimizationSolution,s::Symbol)
    if s === :minimizer
        return getfield(x,:u)
    end
    return getfield(x,s)
end

function Base.summary(io::IO, A::AbstractOptimizationSolution)
    type_color, no_color = get_colorizers(io)
    print(io,
        type_color, nameof(typeof(A)),
        no_color, " with uType ",
        type_color, eltype(A.u),
        no_color)
end


text:
"""
Source: https://github.com/FermiQC/Fermi.jl/discussions/117

Credit to: Gustavo Aroeira (https://github.com/gustavojra)
"""
function molecular_orbital_hamiltonian_coefficients(;
  molecule,
  basis="sto-3g",
  diis=diis(Molecule(molecule)),
  oda=oda(Molecule(molecule))
)

  @suppress begin
    Fermi.Options.set("molstring", ITensorChemistry.molecule(Molecule(molecule)))
    Fermi.Options.set("basis", basis)
    Fermi.Options.set("diis", diis)
    Fermi.Options.set("oda", oda)
  end
  hf_wfn = @suppress Fermi.HartreeFock.RHF()
  aoints = @suppress Fermi.Integrals.IntegralHelper(orbitals=hf_wfn.orbitals)

  hα = @suppress aoints["T"] + aoints["V"]
  gα = @suppress aoints["ERI"]

  # Convert `g` from Chemist convention to Physicist convention
  gα = 0.5 * permutedims(gα, (3, 2, 1, 4))

  nαocc = hf_wfn.molecule.Nα
  hartree_fock_energy = hf_wfn.energy
  nuclear_energy = aoints.molecule.Vnuc

  return (; hα, gα, nαocc, hartree_fock_energy, nuclear_energy)
end

function _molecular_orbital_hamiltonian(hα, gα, nuclear_energy)
  # Representation of the second quantized quantum chemistry Hamiltonian.
  hamiltonian = OpSum()
  hamiltonian += nuclear_energy, "Id", 1
  nα = size(hα, 1)
  for i in 1:nα, j in 1:nα
    if norm(hα[i, j]) > 1e-15
      hamiltonian .+= hα[i, j], "c†↑", i, "c↑", j
      hamiltonian .+= hα[i, j], "c†↓", i, "c↓", j
    end
  end
  for i in 1:nα, j in 1:nα, k in 1:nα, l in 1:nα
    if norm(gα[i, j, k, l]) > 1e-15
      if (i ≠ j) && (k ≠ l) # Otherwise the terms are exactly zero
        hamiltonian .+= gα[i, j, k, l], "c†↑", i, "c†↑", j, "c↑", k, "c↑", l
        hamiltonian .+= gα[i, j, k, l], "c†↓", i, "c†↓", j, "c↓", k, "c↓", l
      end
      hamiltonian .+= gα[i, j, k, l], "c†↑", i, "c†↓", j, "c↓", k, "c↑", l
      hamiltonian .+= gα[i, j, k, l], "c†↓", i, "c†↑", j, "c↑", k, "c↓", l
    end
  end
  return hamiltonian
end

function proc(p, ptot, i)
  return p == mod1(i, ptot)
end

# i ≥ j
f(i, j) = (i - 1) * i ÷ 2 + j
function proc(p, ptot, i, j)
  if i ≤ j
    return p == mod1(f(j, i), ptot)
  end
  return p == mod1(f(i, j), ptot)
end

function proc_sorted(p, ptot, i, j, k, l)
  if j == k
    return proc(p, ptot, j)
  end
  return proc(p, ptot, i, j)
end

function proc(p, ptot, i, j, k, l)
  return proc_sorted(p, ptot, sort((i, j, k, l))...)
end

# From: https://arxiv.org/abs/2103.09976
function _molecular_orbital_hamiltonian(hα, gα, nuclear_energy, nsub_hamiltonians)
  # Representation of the second quantized quantum chemistry Hamiltonian.
  # Split over `nsub_hamiltonians` to be parallelized over.
  hamiltonian = [OpSum() for _ in 1:nsub_hamiltonians]
  hamiltonian[1] += nuclear_energy, "Id", 1
  nα = size(hα, 1)
  for p in 1:nsub_hamiltonians
    for i in 1:nα, j in 1:nα
      pij = 1 / 2 * (proc(p, nsub_hamiltonians, i) + proc(p, nsub_hamiltonians, j))
      if pij ≠ 0 && norm(hα[i, j]) > 1e-10
        hamiltonian[p] .+= pij * hα[i, j], "c†↑", i, "c↑", j
        hamiltonian[p] .+= pij * hα[i, j], "c†↓", i, "c↓", j
      end
    end
    for i in 1:nα, j in 1:nα, k in 1:nα, l in 1:nα
      pijkl = proc(p, nsub_hamiltonians, i, j, k, l)
      if pijkl ≠ 0 && norm(gα[i, j, k, l]) > 1e-10
        if (i ≠ j) && (k ≠ l) # Otherwise the terms are exactly zero
          hamiltonian[p] .+= gα[i, j, k, l], "c†↑", i, "c†↑", j, "c↑", k, "c↑", l
          hamiltonian[p] .+= gα[i, j, k, l], "c†↓", i, "c†↓", j, "c↓", k, "c↓", l
        end
        hamiltonian[p] .+= gα[i, j, k, l], "c†↑", i, "c†↓", j, "c↓", k, "c↑", l
        hamiltonian[p] .+= gα[i, j, k, l], "c†↓", i, "c†↑", j, "c↑", k, "c↓", l
      end
    end
  end
  return hamiltonian
end

function molecular_orbital_hamiltonian(nsub_hamiltonians=nothing; kwargs...)
  (; hα, gα, nαocc, hartree_fock_energy, nuclear_energy) = molecular_orbital_hamiltonian_coefficients(; kwargs...)

  if isnothing(nsub_hamiltonians)
    hamiltonian = _molecular_orbital_hamiltonian(hα, gα, nuclear_energy)
  else
    hamiltonian = _molecular_orbital_hamiltonian(hα, gα, nuclear_energy, nsub_hamiltonians)
  end

  nα = size(hα, 1)
  αocc_state = [i in 1:nαocc ? 1 : 0 for i in 1:nα]
  occ_to_state = Dict([0 => 1, 1 => 4])
  state = [occ_to_state[n] for n in αocc_state]
  return (; hamiltonian, state, hartree_fock_energy)
end


text:
"""
    SplitBinion{T <: Real} = Hyperbolic{T}

A split-binion is a hyperbolic Cayley-Dickson construct with real numbers.
"""
const SplitBinion{T <: Real} = Hyperbolic{T}

function SplitBinion(a::T, b::T) where T <: Real
    SplitBinion{T}(a, b)
end

function SplitBinion(a::T) where T <: Real
    SplitBinion{T}(a, zero(T))
end

function SplitBinion(a::Real, b::Real)
    SplitBinion(promote(a, b)...)
end

"""
    unreal(z::SplitBinion)

The unreal part of a split-binion is a 1-vector.
"""
function unreal(z::SplitBinion{T}) where T <: Real
    T[z.r]
end

function show(io::IO, z::SplitBinion)
    print(io, "[1: ")
    print(io, z.l)
    print(io, ", s: ")
    print(io, z.r)
    print(io, "]")
end

function random(::Type{SplitBinion{T}}) where T <: Real
    SplitBinion{T}(rand(T), rand(T))
end

function random(::Type{SplitBinion{BigInt}})
    SplitBinion{BigInt}(randomBigInt(), randomBigInt())
end

function random(::Type{SplitBinion{BigFloat}})
    SplitBinion{BigFloat}(randomBigFloat(), randomBigFloat())
end

function random(::Type{SplitBinion{Rational{BigInt}}})
    SplitBinion{Rational{BigInt}}(randomBigRational(), randomBigRational())
end

function convert(::Type{SplitBinion{T}}, a::Real) where T <: Real
    SplitBinion{T}(a, zero(T))
end

function convert(::Type{SplitBinion{T}}, z::SplitBinion) where T <: Real
    SplitBinion{T}(z.l, z.r)
end

function promote_rule(::Type{SplitBinion{T}}, ::Type{S}) where {T <: Real, S <: Real}
    SplitBinion{promote_type(T, S)}
end

function promote_rule(::Type{SplitBinion{T}}, ::Type{SplitBinion{S}}) where {T <: Real, S <: Real}
    SplitBinion{promote_type(T, S)}
end

function selfstar(z::SplitBinion)
    (z + star(z)) / 2
end

function anti_selfstar(z::SplitBinion)
    (z - star(z)) / 2
end


text:
mutable struct ScaledOrthogonalMatrix <: MyXMLElement
    el::XMLOrNothing
    scale::Parameter
    U::MatrixParameter
    id::String

    function ScaledOrthogonalMatrix(scale::Parameter, U::MatrixParameter,
                                    id::String;
                                    tol::Float64=1e-10)

        scale_data = scale.val
        U_data = U.mat
        k = length(scale_data)
        k2, p = size(U_data)

        @assert k == k2

        # check U is orthonormal
        for i = 1:k
            ui = @view U_data[i, :]

            if abs(dot(ui, ui) - 1.0) > tol
                error("matrix is not orthonormal")
            end

            for j = (i + 1):k
                uj = @view U_data[j, :]
                if abs(dot(ui, uj)) > tol
                    error("matrix is not orthonormal")
                end
            end
        end

        return new(nothing, scale, U, id)
    end

end

function ScaledOrthogonalMatrix(L::Matrix{Float64}, id, scale_id, U_id)
    Lsvd = svd(L)
    return ScaledOrthogonalMatrix(Parameter(Lsvd.S, scale_id),
                                  MatrixParameter(Lsvd.Vt, U_id),
                                  id)
end

function make_xml(som::ScaledOrthogonalMatrix)
    make_xml(som.scale)
    make_xml(som.U)
    el = new_element(bn.SCALED_MATRIX)
    set_id!(el, som.id)
    scale_el = new_child(el, bn.SCALE)
    add_ref_el(scale_el, som.scale.el)
    mat_el = new_child(el, bn.MATRIX)
    add_ref_el(mat_el, som.U.el)

    som.el = el
    return el
end

function get_loggables(som::ScaledOrthogonalMatrix)
    make_xml(som)
    return [get_loggables(som.scale); som.el]
end

function set_mat!(som::ScaledOrthogonalMatrix, mat::AbstractMatrix{Float64})
    s = svd(mat)
    set_mat!(som.U, s.Vt)
    set_value(som.scale, s.S)
end

function set_scale!(som::ScaledOrthogonalMatrix, scale::Vector{Float64})
    set_value(som.scale, scale)
end

function size(som::ScaledOrthogonalMatrix, args...)
    return size(som.U, args...)
end

function name(::ScaledOrthogonalMatrix)
    return bn.SCALED_MATRIX
end


text:
## Figures for the Quadratic Equation slides
using Plots
using LaTeXStrings
pyplot() # nicer output
theme(:dark)
# Need for 3d arrow
using PyPlot

##
A = [3 -1.5; -1.5 2]
b = [-1; 1]
quadratic_ls(A,b) = x-> 0.5*x'*A*x - x'*b
ezsurf(x, y, f) = begin
    X = repeat(x', length(y), 1)
    Y = repeat(y, 1, length(x))
    # Evaluate each f(x, y)
    Z = map((x,y) -> f([x,y]), X, Y)
    plot(x, y, Z, st=:surface)
end
ezsurf(-3:0.1:3,-3:0.1:3,
  quadratic_ls(A,b))
plot!(colorbar=false,size=(500,450))
# Show the gradient at x = [1,1]
x = [-2,2]
g = A*x - b
g = 0.1.*g

plot!([x[1],x[1]+g[1]],[x[2],x[2]+g[2]],[quadratic_ls(A,b)(x), quadratic_ls(A,b)(x+g)],
  line=(arrow=:arrow),label="",color=2)
  plot!([x[1],x[1]-g[1]],[x[2],x[2]-g[2]],[quadratic_ls(A,b)(x), quadratic_ls(A,b)(x-g)],
    line=(arrow=:arrow),label="",color=3)
##
anim = @animate for i=1:360
  plot!(camera=(i,30))
end
gif(anim, "quadratic_ls1.gif")
##
using Printf
using LinearAlgebra
""" Run gradient descent to solve a linear system of equations.

This function assumes that A is symmetric positive definite. """
function gradient_descent_linsys(A,b;step=0.1, maxiter=10*size(A,1),tol=1e-6)
  x = zeros(size(A,1)) # start at 0
  for i=1:maxiter
    g = A*x - b
    @printf("iteration %4i  residual = %.3e\n", i, norm(g))
    if norm(g) <= tol
      break
    end
    x -= step*g
  end
  return x
end
gradient_descent_linsys(A,b)
## Animate the algorithm
ezsurf(-0.5:0.025:0.5,-0.5:0.025:0.5,
  quadratic_ls(A,b))
plot!(colorbar=false,size=(500,450))
plot!(camera=(170,30))
x = zeros(size(A,1)) # initialize
p = plot3d!(1, label="", marker = 2)
step = 0.1
anim = @animate for i=1:20
  global x
  g = A*x - b
  x -= step*g
  push!(p[1][2], x[1], x[2], quadratic_ls(A,b)(x))
  title!(@sprintf("Iteration %3i  Residual = %.3e", i, norm(g)))
end
gif(anim, "gradient_descent_ls1.gif")


text:
export Alcohol, Side, DrinkComponent, Recipe, DrinkBook, BarCabinet, Systembolag, Supermarket

abstract type DataLayer end    # For dispatch
abstract type Presenter end     # For dispatch
abstract type Ingredient end

struct Alcohol <: Ingredient
    name::String
end

struct Side <: Ingredient
    name::String
end

struct DrinkComponent
    ingredient::Ingredient
    mandatory::Bool
    quantity
end

struct Recipe
    name::String
    ingredients::Vector{DrinkComponent}
    description::String
end

struct DrinkBook
    recipes::Vector{Recipe}
end

struct BarCabinet
    alcohols::Vector{Alcohol}
    sides::Vector{Side}
end

function BarCabinet(ing::Vector{Ingredient})
    BarCabinet(alcohols(ing), sides(ing))
end

struct Systembolag
    alcohols::Vector{Alcohol}
end
Systembolag(ing::Vector{Ingredient}) = Systembolag(alcohols(ing))

struct Supermarket
    sides::Vector{Side}
end
Supermarket(ing::Vector{Ingredient}) = Supermarket(sides(ing))

const ContainerType = Union{DrinkBook, BarCabinet, Systembolag, Supermarket}
# All types in ContainerType implements push!, merge!

const ElementType = Union{Ingredient, Alcohol, Side, DrinkComponent, Recipe}
# All types in ElementType implements name(e), similarname(e)


text:
using KernelFunctions
using AxisArrays
using Distances
using Documenter
using Functors: functor
using Kronecker: Kronecker
using LinearAlgebra
using LogExpFunctions
using PDMats
using Random
using SpecialFunctions
using Test
using Zygote: Zygote
using ForwardDiff: ForwardDiff
using ReverseDiff: ReverseDiff
using FiniteDifferences: FiniteDifferences

using KernelFunctions: SimpleKernel, metric, kappa, ColVecs, RowVecs, TestUtils

using KernelFunctions.TestUtils: test_interface

# The GROUP is used to run different sets of tests in parallel on the GitHub Actions CI.
# If you want to introduce a new group, ensure you also add it to .github/workflows/ci.yml
const GROUP = get(ENV, "GROUP", "")

# Writing tests:
# 1. The file structure of the test should match precisely the file structure of src.
#   Amongst other things, this means that there should be exactly 1 test file per src file.
#   This makes it trivially easy for someone to find the tests associated with a particular
#   src file.
# 2. A consequence of 1 is that there should be exactly 1 test file per src file.
# 3. A test file called foo.jl should have the structure:
#   @testset "foo" begin
#       code
#   end
#
#   Note that the testset is called `foo`, not `foo.jl`. Use whatever testset structure
#   seems appropriate within a given file. eg. if multiple types / functions are defined in
#   a particular source file, you might want multiple testsets in the test file.
# 4. Each directory should have its own testset, in which each test file is `include`d.
# 5. Each test file should create its own state, and shouldn't rely on state defined in
#   other test files. e.g. don't define a matrix used by all of the files in kernels. If
#   two test files are similar enough to share state, perhaps the corresponding source code
#   should be in the same file.
# 6. If you write a src file without any tests, create a corresponding test file with the
#   usual structure, but without any tests.
# 7. Explicitly create a new random number generate for _at_ _least_ each new test file, and
#   use it whenever generating randomness. This ensures complete control over random number
#   generation and makes it clear what randomness depends on other randomness.
# 8. All `using` statements should appear in runtests.jl.
# 9. List out all test files explicitly (eg. don't loop over them). This makes it easy to
#   disable tests by simply commenting them out, and makes it very clear which tests are not
#   currently being run.
# 10. If utility functionality is required, it should be placed in `src/test_utils.jl` so
#   that other packages can benefit from it when implementing new kernels.
@info "Packages Loaded"

include("test_utils.jl")

@testset "KernelFunctions" begin
    if GROUP == "" || GROUP == "Transform"
        @testset "transform" begin
            include("transform/transform.jl")
            print(" ")
            include("transform/scaletransform.jl")
            print(" ")
            include("transform/ardtransform.jl")
            print(" ")
            include("transform/lineartransform.jl")
            print(" ")
            include("transform/functiontransform.jl")
            print(" ")
            include("transform/selecttransform.jl")
            print(" ")
            include("transform/chaintransform.jl")
            print(" ")
            include("transform/periodic_transform.jl")
            print(" ")
            include("transform/with_lengthscale.jl")
            print(" ")
        end
        @info "Ran tests on Transform"
    end

    if GROUP == "" || GROUP == "BaseKernels"
        @testset "basekernels" begin
            include("basekernels/constant.jl")
            print(" ")
            include("basekernels/cosine.jl")
            print(" ")
            include("basekernels/exponential.jl")
            print(" ")
            include("basekernels/exponentiated.jl")
            print(" ")
            include("basekernels/fbm.jl")
            print(" ")
            include("basekernels/gabor.jl")
            print(" ")
            include("basekernels/matern.jl")
            print(" ")
            include("basekernels/nn.jl")
            print(" ")
            include("basekernels/periodic.jl")
            print(" ")
            include("basekernels/piecewisepolynomial.jl")
            print(" ")
            include("basekernels/polynomial.jl")
            print(" ")
            include("basekernels/rational.jl")
            print(" ")
            include("basekernels/sm.jl")
            print(" ")
            include("basekernels/wiener.jl")
            print(" ")
        end
        @info "Ran tests on BaseKernel"
    end

    if GROUP == "" || GROUP == "Kernels"
        @testset "kernels" begin
            include("kernels/kernelproduct.jl")
            include("kernels/kernelsum.jl")
            include("kernels/kerneltensorproduct.jl")
            include("kernels/overloads.jl")
            include("kernels/scaledkernel.jl")
            include("kernels/transformedkernel.jl")
            include("kernels/normalizedkernel.jl")
            include("kernels/neuralkernelnetwork.jl")
            include("kernels/gibbskernel.jl")
        end
        @info "Ran tests on Kernel"
    end

    if GROUP == "" || GROUP == "MultiOutput"
        @testset "multi_output" begin
            include("mokernels/moinput.jl")
            include("mokernels/independent.jl")
            include("mokernels/slfm.jl")
            include("mokernels/intrinsiccoregion.jl")
            include("mokernels/lmm.jl")
        end
        @info "Ran tests on Multi-Output Kernels"
    end

    if GROUP == "" || GROUP == "Others"
        include("utils.jl")

        @testset "distances" begin
            include("distances/pairwise.jl")
            include("distances/dotproduct.jl")
            include("distances/delta.jl")
            include("distances/sinus.jl")
        end
        @info "Ran tests on Distances"

        @testset "matrix" begin
            include("matrix/kernelmatrix.jl")
            include("matrix/kernelkroneckermat.jl")
            include("matrix/kernelpdmat.jl")
        end
        @info "Ran tests on matrix"

        @testset "approximations" begin
            include("approximations/nystrom.jl")
        end

        include("generic.jl")
        include("chainrules.jl")
        include("zygoterules.jl")

        @testset "doctests" begin
            DocMeta.setdocmeta!(
                KernelFunctions,
                :DocTestSetup,
                quote
                    using KernelFunctions
                    using LinearAlgebra
                    using Random
                    using PDMats: PDMats
                end;
                recursive=true,
            )
            doctest(
                KernelFunctions;
                doctestfilters=[
                    r"{([a-zA-Z0-9]+,\s?)+[a-zA-Z0-9]+}",
                    r"(\s?Array{[a-zA-Z0-9]+,\s?1}|\s?Vector{[a-zA-Z0-9]+})",
                    r"(\s?Array{[a-zA-Z0-9]+,\s?2}|\s?Matrix{[a-zA-Z0-9]+})",
                ],
            )
        end
    end
end


text:
using Plots, Random

function angle2point(θ)
    N = length(θ)
    x = zeros(N,2)
    for i in 1:N
        x[i,:] = [cos(θ[i]), sin(θ[i])]
    end
    return x
end

function circle_plot(x; color=:red)
    scatter(x[:,1], x[:,2]; c=color, aspect_ratio=1, legend=false, grid=false, mswidth=0, xlims=[-1.1, 1.1], ylims=[-1.1, 1.1])
end

function state_vector(θx)
    sort!(θx)
    γ = zeros(N)
    γ[1] = θx[1]-θx[end]+2π
    for i in 2:N
        γ[i] = θx[i]-θx[i-1]
    end
    𝝓 = [γ[i]-2π/N for i in 1:N]
    return 𝝓
end

# 0. Initialize points on [0, 2π)
Random.seed!(2020)
N = 32; iter = 5000
θx = 2π*rand(N)
circle_plot(angle2point(θx))

# 1. Learning steps
α = 0.1
𝝓_norm = []
anim = @animate for i in 1:iter
    push!(𝝓_norm, norm(state_vector(θx)))
    θ = 2π*rand()
    θy = (θx .+ θ) .% 2π
    nnx = Dict()
    for j in 1:N
        nnx[j] = []
    end
    for i in 1:N
        arr = min.(abs.(θx .- θy[i]), 2π .- abs.(θx .- θy[i]))
        min_val = minimum(arr)
        min_ind = findall(arr .== min_val)
        for j in min_ind
            nnx[j] = push!(nnx[j], i)
        end
    end
    for j in 1:N
        if nnx[j] != []
            θx[j] += α * sum([θy[i]-θx[j]+(abs(θy[i]-θx[j])<π ? 0 : -sign(θy[i]-θx[j])*2π)  for i in nnx[j]])
        end
    end
    if i % 50 == 0
        circle_plot(angle2point(θx))
        title!("iter=$(i)")
    end
end when i % 50 == 0

gif(anim, "selfSimNN_1D.gif", fps = 5)

plot(𝝓_norm, legend=false)

# circle_plot(angle2point(θx))


text:
# This file was generated, do not modify it. # hide
@load RandomForestRegressor pkg=ScikitLearn

rf_mdl = RandomForestRegressor()
rf = machine(rf_mdl, X, y)
fit!(rf, rows=train)

ypred = predict(rf, rows=test)
round(rms(ypred, y[test]), sigdigits=3)


text:
using PyCall

let

println("Before we start, here's a word from our engine, pygame:")
# Use pycall to import pygame
pygame = pyimport("pygame")

pygame.init()

# define basic colors
white = (255, 255, 255)
black = (0, 0, 0)
red = (255, 0, 0)
green = (76, 175, 80)

# define screen_size, block_size, speed, FPS etc.
screen_size = (400, 400)
block_size = 20

max_speed = 20

FPS = 5

# Import all different sprites
img = pygame.image.load("./sprites/head.png")
head_1 = pygame.image.load("./sprites/head_1.png")

apple_img = pygame.image.load("./sprites/apple.png")

# Setup game window
game_display = pygame.display.set_mode(screen_size)
pygame.display.set_caption("Snake")
clock = pygame.time.Clock()

game_exit = false

# Start snake at middle of screen
head_x = (screen_size[1] / 2) - block_size
head_y = (screen_size[2] / 2) - block_size
x_speed = 0
y_speed = 0

# Define function to draw snake
function snake(snake_list, block_size)
    # Load appropriate image for head according to size of snake and orientation of snake
    head = img

    if length(snake_list) == 1
        head = head_1
        if x_speed > 0
            head = pygame.transform.rotate(head_1, 270)
        elseif x_speed < 0
            head = pygame.transform.rotate(head_1, 90)
        elseif y_speed > 0
            head = pygame.transform.rotate(head_1, 180)
        elseif y_speed < 0
            head = head_1
        end
    else
        head = img
    # Rotate snake head according to where snake is heading
        if x_speed > 0
            head = pygame.transform.rotate(img, 270)
        elseif x_speed < 0
            head = pygame.transform.rotate(img, 90)
        elseif y_speed > 0
            head = pygame.transform.rotate(img, 180)
        elseif y_speed < 0
            head = img
        end
    end
    # Draw image for head, and render the rest of the snake as rects
    game_display.blit(head, snake_list[end])
    for snake_block in snake_list[1:end - 1]
        game_display.fill(green, rect = [snake_block[1], snake_block[2], block_size, block_size])
    end
end

# Create function to get user input and recursively validate the user input
function get_input(text_to_print, in_same_line, allowed_options, required_params)

    # If we want to print on the same line, use `print()`. Else, `use println()`
    in_same_line ? print(text_to_print) : println(text_to_print)

    # Get input and split it into an array
    txtin = split(lowercase(readline()));

    # If the required amount of parameters is not met or if the command is not in the list of allowed options,
    # ask for input again
    if length(txtin) < required_params || !(txtin[1] in allowed_options)
        println("Please enter a valid command (see README)")
        get_input(text_to_print, in_same_line, allowed_options, required_params)
    else
        # If input is valid, return the input as a string.
        return join(txtin, " ")
    end
end

# Ask if user wants to play_again
function play_again()

    ans = get_input("Play again? (y/n): ", true, ["y", "n"], 1)

    # If user wants to play again, run the game again. Else, exit()
    if ans == "y"
        run(`julia $PROGRAM_FILE`)
        pygame.quit()
        exit()
    else
        exit()
    end
end

# Start game
function game_loop(color_mode)
    # Decide random position of apple
    apple_pos = (rand(0 : (screen_size[1] / block_size) - 1) * block_size, rand(0 : (screen_size[1] / block_size) - 1) * block_size)

    # init snake
    snake_list = []

    # Size of snake
    score = 1

    # While user is not exiting, play game
    while !game_exit
        # Check every event pyagem gets
        for event in pygame.event.get()
            # If user wants to quit, quit
            if event.type == pygame.QUIT
                return score
            end

            # According to key pressed, move snake
            if event.type == pygame.KEYDOWN
                if event.key == pygame.K_LEFT
                    x_speed = -max_speed
                    y_speed = 0
                elseif event.key == pygame.K_RIGHT
                    x_speed = max_speed
                    y_speed = 0
                elseif event.key == pygame.K_DOWN
                    x_speed = 0
                    y_speed = max_speed
                elseif event.key == pygame.K_UP
                    x_speed = 0
                    y_speed = -max_speed
                end
            end

        end

        # Move snake's head position
        head_x += x_speed
        head_y += y_speed

        # Draw the background
        if color_mode == "1"
            game_display.fill(white)
        else
            game_display.fill(black)
        end

        # Draw apple sprite at random location
        game_display.blit(apple_img, [apple_pos[1], apple_pos[2], block_size, block_size])

        # Put snake's head in snake_list
        snake_head = [head_x, head_y]
        push!(snake_list, snake_head)

        # If the snake positions that we store are more than the length, remove them (the positions are no more snake_segments)
        if length(snake_list) > score
            deleteat!(snake_list, 1)
        end

        # Loop snake back into game
        for snake_segment in snake_list
            snake_segment[1] >= 0 ? snake_segment[1] %= screen_size[1] : snake_segment[1] = screen_size[1] - snake_segment[1]
            snake_segment[2] >= 0 ? snake_segment[2] %= screen_size[2] : snake_segment[2] = screen_size[2] - snake_segment[2]
        end

        # Loop head back into game
        head_x > 0 ? head_x %= screen_size[1] : head_x = screen_size[1] - head_x
        head_y > 0 ? head_y %= screen_size[2] : head_y = screen_size[2] - head_y

        # If any snake_segment except head is touching head, end game
        for snake_segment in snake_list[1: end - 1]
            if snake_segment == snake_list[end]
                return score
            end
        end

        # Render snake
        snake(snake_list, block_size)

        # Update display
        pygame.display.update()

        # If snake eats apple, increment score by 1, and draw new apple
        for snake_segment in snake_list
            if snake_segment[1] == apple_pos[1] && snake_segment[2] == apple_pos[2]
                apple_pos = (rand(0 : (screen_size[1] / block_size) - 1) * block_size, rand(0 : (screen_size[1] / block_size) - 1) * block_size)
                score += 1
            end
        end

        # Wait for 1/FPS seconds
        clock.tick(FPS)
    end
end

a = """
Hi! Welcome to Snake. Before we start, here's some instructions:
- Don't run into yourself
- Running into the walls loops you back into the screen
- Eat as many apples as you can
- Press any arrow key to start playing

"""
println(a)

mode = get_input("Would you like to play in light mode (1) or dark mode (2)?", false, ["1", "2"], 1)

# Run game
user_score = game_loop(mode)

println("Your score was: $user_score")

play_again()
end

max_stars_repo_name: soumitradev/Snake

max_issues_repo_name: soumitradev/Snake

max_forks_repo_name: soumitradev/Snake


text:
function interact(M, N)
    d = 1.0 + (M * N)^0.75 * 2.01E-5 + (M * N)^1.52 * M * 5.31E-15
    return 0.47 * log(d)
end

@defcomp climateforcing begin
    # Atmospheric CO2 concentration
    acco2 = Parameter(index=[time])

    # Pre-industrial atmospheric CO2 concentration
    co2pre = Parameter(default = 275.0)

    # Atmospheric CH4 concentration
    acch4 = Parameter(index=[time])

    # Pre-industrial atmospheric CH4 concentration
    ch4pre = Parameter(default = 790.0)

    # Indirect radiative forcing increase for CH4
    ch4ind = Parameter(default = 0.4)

    # Atmospheric N2O concentration
    acn2o = Parameter(index=[time])

    # Pre-industrial atmospheric N2O concentration
    n2opre = Parameter(default = 285.0)

    # Atmospheric SF6 concentrations
    acsf6 = Parameter(index=[time])

    # Pre-industrial atmospheric SF6 concentration
    sf6pre = Parameter(default = 0.04)

    # Radiative forcing from CO2
    rfco2 = Variable(index=[time])

    # Radiative forcing from CH4
    rfch4 = Variable(index=[time])

    # Radiative forcing from N2O
    rfn2o = Variable(index=[time])

    # Radiative forcing from N2O
    rfsf6 = Variable(index=[time])

    # Radiative forcing from SO2
    rfso2 = Parameter(index=[time])

    # Radiative forcing
    radforc = Variable(index=[time])

    # EMF22 radiative forcing
    rfemf22 = Variable(index=[time])

    function run_timestep(p, v, d, t)

        if !is_first(t)
            ch4n2o = interact(p.ch4pre, p.n2opre)

            v.rfco2[t] = 5.35 * log(p.acco2[t] / p.co2pre)

            v.rfch4[t] = 0.036 * (1.0 + p.ch4ind) * (sqrt(p.acch4[t]) - sqrt(p.ch4pre)) - interact(p.acch4[t], p.n2opre) + ch4n2o

            v.rfn2o[t] = 0.12 * (sqrt(p.acn2o[t]) - sqrt(p.n2opre)) - interact(p.ch4pre, p.acn2o[t]) + ch4n2o

            v.rfsf6[t] = 0.00052 * (p.acsf6[t] - p.sf6pre)

            v.radforc[t] = v.rfco2[t] + v.rfch4[t] + v.rfn2o[t] + v.rfsf6[t] + p.rfso2[t]

            v.rfemf22[t] = v.rfco2[t] + v.rfch4[t] + v.rfn2o[t]
        end
    end
end


text:
function IterationManager(arg0::jint)
    return IterationManager((jint,), arg0)
end

function IterationManager(arg0::jint, arg1::Incrementor_MaxCountExceededCallback)
    return IterationManager((jint, Incrementor_MaxCountExceededCallback), arg0, arg1)
end

function add_iteration_listener(obj::IterationManager, arg0::IterationListener)
    return jcall(obj, "addIterationListener", void, (IterationListener,), arg0)
end

function fire_initialization_event(obj::IterationManager, arg0::IterationEvent)
    return jcall(obj, "fireInitializationEvent", void, (IterationEvent,), arg0)
end

function fire_iteration_performed_event(obj::IterationManager, arg0::IterationEvent)
    return jcall(obj, "fireIterationPerformedEvent", void, (IterationEvent,), arg0)
end

function fire_iteration_started_event(obj::IterationManager, arg0::IterationEvent)
    return jcall(obj, "fireIterationStartedEvent", void, (IterationEvent,), arg0)
end

function fire_termination_event(obj::IterationManager, arg0::IterationEvent)
    return jcall(obj, "fireTerminationEvent", void, (IterationEvent,), arg0)
end

function get_iterations(obj::IterationManager)
    return jcall(obj, "getIterations", jint, ())
end

function get_max_iterations(obj::IterationManager)
    return jcall(obj, "getMaxIterations", jint, ())
end

function increment_iteration_count(obj::IterationManager)
    return jcall(obj, "incrementIterationCount", void, ())
end

function remove_iteration_listener(obj::IterationManager, arg0::IterationListener)
    return jcall(obj, "removeIterationListener", void, (IterationListener,), arg0)
end

function reset_iteration_count(obj::IterationManager)
    return jcall(obj, "resetIterationCount", void, ())
end