text:
module StreamReader

export PartsIterator,
    ReaderIterator,
    IOReaderIterator,
    IOStringReaderIterator,
    ReaderListIterator,
    DEFAULT_PART_SIZE,
    loadedpct,
    calculate

const DEFAULT_PART_SIZE = 1024
const N = Void
const N_TYPE = Type{N}
const STR_TYPE = isdefined(:AbstractString)? AbstractString : String

"""
Calculate sizes.
Return tuple of:
    1: Number of parts
    2: The parts size
    3: The last part size
"""
function calculate(total_size::Integer, part_size::Integer=0)
    part_size = (part_size == 0 ? DEFAULT_PART_SIZE : part_size)

    if part_size > total_size
        part_size = total_size
    end

    parts = convert(Integer, ceil(total_size / part_size))
    last_part_size = total_size % part_size
    if last_part_size == 0
        last_part_size = part_size
    end
    (parts, part_size, last_part_size)
end

"""
Iterator over parts
"""
type PartsIterator
    size::Integer
    part_size::Integer
    loaded::Integer
    length::Integer
    part::Integer
    current_size::Integer
    left::Integer
    last_part_size::Integer

    PartsIterator(s::Integer, ps::Integer = 0) = begin
        l, ps, lps = calculate(s, ps)
        new(s, ps, 0, l, 0, 0, 0, lps)
    end
end

Base.start(p::PartsIterator) = begin
    s = p.size < p.part_size ? p.size : p.part_size
    p.loaded = 0
    p.part = 0
    p.left = p.size
    s
end

Base.done(p::PartsIterator, state) = p.loaded == p.size

Base.length(p::PartsIterator) = p.length

Base.next(p::PartsIterator, state) = begin
    p.current_size = state
    p.left -= state
    p.part += 1
    p.loaded += state
    next_state = p.left < state ? p.left : state
    state, next_state
end

function loadedpct(pi::PartsIterator)
    pi.loaded * 100 / pi.size
end

type ReaderIterator
    parts::PartsIterator
    prestart::Union(N_TYPE,Function)
    postdone::Union(N_TYPE,Function)
    read::Union(N_TYPE,Function)
    iterating::Bool
    num::Integer

    ReaderIterator(bp::PartsIterator; prestart::Union(N_TYPE,Function) = N,
        postdone::Union(N_TYPE,Function) = N) = new(bp, prestart,
        postdone, N, false, 0)
    ReaderIterator(ps::Function, bp::PartsIterator; kwargs...) = ReaderIterator(bp;
        prestart=ps, kwargs...)
end

Base.start(sr::ReaderIterator) = begin
    state = start(sr.parts)

    if sr.prestart != N
        sr.prestart(sr)
    end

    sr.iterating = true
    sr.num = 0

    return state
end

Base.done(ri::ReaderIterator, state) = begin
    d = done(ri.parts, state)

    if d && ri.iterating
        ri.iterating = false

        if ri.postdone != N
            ri.postdone(ri)
        end
    end

    return d
end

Base.next(ri::ReaderIterator, state) = begin
    ri.num += 1
    ri.read(state), next(ri.parts, state)[2]
end

loaded(ri::ReaderIterator) = loaded(ri.parts)
loadedpct(ri::ReaderIterator) = loadedpct(ri.parts)

type ReaderListIterator
    makereader::Function
    count::Integer
    num::Integer
    reader::Union(N_TYPE, ReaderIterator)

    ReaderListIterator(mr::Function, c::Integer) = new(mr, c, 0, N)
end

Base.start(ri::ReaderListIterator) = begin
    ri.num = 0
    ri.reader = N
    return 1
end

Base.done(ri::ReaderListIterator, state) = begin
    return state > ri.count
end

Base.next(ri::ReaderListIterator, state) = begin
    ri.num = state
    ri.reader = ri.makereader(ri)
    (state, ri.reader), state+1
end

function setdatatype{T}(ri::ReaderIterator, datatype::Type{T})
    oldread = ri.read
    ri.read = (size) -> datatype(oldread(size))
end

function IOReaderIterator{T,D}(io::IO, bp::PartsIterator; read_type::Type{T} = Uint8,
    datatype::Type{D} = N_TYPE, prestart::Union(N_TYPE,Function) = N, kwargs...)
    ReaderIterator(bp; kwargs...) do ri
        ri.read = (size) -> Base.read(io, read_type, size)

        if datatype != N_TYPE
            setdatatype(ri, datatype)
        end

        if prestart != N
            prestart(ri)
        end
    end
end

IOReaderIterator(io::IO, len::Integer, bsize::Integer=0; kwargs...) = IOReaderIterator(io,
    PartsIterator(len, bsize); kwargs...)

function IOStringReaderIterator{D<:STR_TYPE}(io::IO, pi::PartsIterator;
    datatype::Type{D} = UTF8String, kwargs...)
    IOReaderIterator(io, pi; datatype=datatype, kwargs...)
end

IOStringReaderIterator(io::IO, len::Integer, bsize::Integer=0; kwargs...) = IOStringReaderIterator(io,
    PartsIterator(len, bsize); kwargs...)

end


text:
using OrdinaryDiffEq, DiffEqProblemLibrary, Base.Test

prob = prob_ode_linear

sol =solve(prob,DP5(),dt=1//2^(2),save_timeseries=false,dense=false)
sol2=solve(prob,DP5(),dt=1//2^(2),save_timeseries=false,dense=false,saveat=[1/2])

@test symdiff(sol.t,sol2.t) == [1/2]

sol3=solve(prob,DP5(),dt=1//2^(2),save_timeseries=false,dense=false,saveat=[1/2],tstops=[1/2])

@test sol3.t == [0.0,0.5,1.00]

sol =solve(prob,DP5(),dt=1//2^(2),save_timeseries=true,dense=true)
sol2=solve(prob,DP5(),dt=1//2^(2),save_timeseries=true,dense=true,saveat=[1/2])

sol2(.49)

interpd = sol2(0:1//2^(4):1)

#plot(0:1//2^(4):1,interpd)

@test symdiff(sol.t,sol2.t) == [1/2]

sol =solve(prob,RK4(),dt=1/2^(2),save_timeseries=true,dense=false)
sol2=solve(prob,RK4(),dt=1/2^(2),save_timeseries=true,dense=false,saveat=[.125,.6,.61,.8])

@test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]

sol =solve(prob,Rosenbrock32(),dt=1/2^(2),save_timeseries=true,dense=false)
sol2=solve(prob,Rosenbrock32(),dt=1/2^(2),save_timeseries=true,dense=false,saveat=[.125,.6,.61,.8])

@test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]

sol =solve(prob,Trapezoid(),dt=1/2^(2),save_timeseries=true,dense=false)
sol2=solve(prob,Trapezoid(),dt=1/2^(2),save_timeseries=true,dense=false,saveat=[.125,.6,.61,.8])

@test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]

prob = prob_ode_2Dlinear

sol =solve(prob,DP5(),dt=1//2^(2),save_timeseries=false,dense=true)
sol2=solve(prob,DP5(),dt=1//2^(2),save_timeseries=false,dense=true,saveat=[1/2])

@test symdiff(sol.t,sol2.t) == [1/2]

sol =solve(prob,DP5(),dt=1//2^(2),save_timeseries=true,dense=false)
sol2=solve(prob,DP5(),dt=1//2^(2),save_timeseries=true,dense=false,saveat=[1/2])

@test symdiff(sol.t,sol2.t) == [1/2]

sol =solve(prob,RK4(),dt=1/2^(2),save_timeseries=true,dense=false)
sol2=solve(prob,RK4(),dt=1/2^(2),save_timeseries=true,dense=false,saveat=[.125,.6,.61,.8])

@test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]

sol =solve(prob,Rosenbrock32(),dt=1/2^(2),save_timeseries=true,dense=false)
sol2=solve(prob,Rosenbrock32(),dt=1/2^(2),save_timeseries=true,dense=false,saveat=[.125,.6,.61,.8])

@test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]

sol =solve(prob,Trapezoid(),dt=1/2^(2),save_timeseries=true,dense=false)
sol2=solve(prob,Trapezoid(),dt=1/2^(2),save_timeseries=true,dense=false,saveat=[.125,.6,.61,.8])

@test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]

sol=solve(prob,Trapezoid(),dt=1/2^(2),save_timeseries=true,dense=false,saveat=[0,.125,.6,.61,.8])

@test !(sol.t[2] ≈ 0)

# Test Iterators

sol2=solve(prob,DP5(),dt=1//2^(2),save_timeseries=false,dense=false,saveat=0:1//100:1)

@test sol2.t ≈ collect(0:1//100:1)

sol2=solve(prob,DP5(),dt=1//2^(2),save_timeseries=false,dense=false,saveat=linspace(0,1,100))

@test sol2.t ≈ linspace(0,1,100)

ext: jl


text:
module BCTRNN

using Distributions
using OrdinaryDiffEq
using StochasticDiffEq
using DiffEqSensitivity
import DiffEqFlux: initial_params, paramlength, FastChain, FastDense, sciml_train, BFGS
using GalacticOptim
using Zygote
using Flux
import Flux.Data: DataLoader
using NNlib: sigmoid
using Random
using ModelingToolkit
using DataInterpolations
using Measurements
using FastBroadcast
using ComponentArrays

using Plots

#using BlackBoxOptim

import IterTools: ncycle

struct VariableLowerBound end
struct VariableUpperBound end
Symbolics.option_to_metadata_type(::Val{:lower}) = VariableLowerBound
Symbolics.option_to_metadata_type(::Val{:upper}) = VariableUpperBound

function get_bounds(sys, s_in; T=Float32, default_lb=T(-Inf), default_ub=T(Inf))
  cell_lb = T[]
  cell_ub = T[]

  params = collect(parameters(sys))[s_in+1:end]
  states = collect(ModelingToolkit.states(sys))
  for v in vcat(params,states)
    contains(string(v), "OutPin") && continue
    lower = hasmetadata(v, VariableLowerBound) ? getmetadata(v, VariableLowerBound) : default_lb
    upper = hasmetadata(v, VariableUpperBound) ? getmetadata(v, VariableUpperBound) : default_ub
    push!(cell_lb, lower)
    push!(cell_ub, upper)
  end
  return cell_lb, cell_ub
end

include("layers.jl")
include("systems/ltc.jl")
include("systems/ltc_orig.jl")
include("systems/ltc_mtk.jl")
include("systems/ltc_gj.jl")
include("systems/ltc_syn_state.jl")
include("systems/ltc_syn_state_mtk.jl")
include("systems/stateful_mtk.jl")
include("systems/dLeak.jl")
include("wiring/basic.jl")
include("model.jl")
include("optimize.jl")
include("losses.jl")
include("callback.jl")
include("utils.jl")

export DataLoader

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "SeqGen"
version = v"1.3.4"

# Collection of sources required to build SFML
sources = [
    GitSource("https://github.com/rambaut/Seq-Gen.git", "bc9d8070b2cd1f1352f74282e5b209302eae38a1")
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/Seq-Gen/source
head -39 seq-gen.c > license
install_license ${WORKSPACE}/srcdir/Seq-Gen/source/license

# compile seq-gen
make -j${nproc}

# move to bindir
install -Dvm 755 "seq-gen" "${bindir}/seq-gen${exeext}"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("seq-gen", :seqgen),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: cmhamel/Yggdrasil

max_issues_repo_name: cmhamel/Yggdrasil

max_forks_repo_name: cmhamel/Yggdrasil


text:
using Statistics
using Plots
using FFTW
using Statistics
using Optim
using Images, FileIO, ImageIO
using Printf
using Revise
using Profile
using LinearAlgebra
using JLD2
using Random
using Distributions
using FITSIO
using LineSearches
using SparseArrays

push!(LOAD_PATH, pwd()*"/../../../main")
using DHC_2DUtils
push!(LOAD_PATH, pwd()*"/../../../scratch_NM")
using Deriv_Utils_New
using Data_Utils
using Visualization
using ReconFuncs

numfile = 1000

logbool = false
apdbool = true
isobool = false
invcovtype = "Full+Eps"

direc = "scratch_NM/StandardizedExp/Nx64/" #Change
datfile = direc * "Data_" * string(numfile) * ".jld2" #Replace w SLURM array
loaddf = load(datfile)
true_img = loaddf["true_img"]
init = loaddf["init"]

fname_save = direc * "SFDTargSFDCov/" * string(logbool) * "_" * string(apdbool) * "_" * string(isobool) * "/LambdaVary/" * string(numfile) * "rand_init"  #Change
Nx=size(true_img)[1]
filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(S1iso, Nf) = size(filter_hash["S1_iso_mat"])
(S2iso, Nfsq) = size(filter_hash["S2_iso_mat"])
dhc_args = Dict(:doS2=>false, :doS20=>true, :apodize=>apdbool, :iso=>isobool) #Iso #CHANGE: Change sig for the sfd data since the noise model is super high and the tiny values make sense
if dhc_args[:iso]
    coeff_mask = falses(2+S1iso+S2iso)
    coeff_mask[S1iso+3:end] .= true
else #Not iso
    coeff_mask = falses(2+Nf+Nf^2)
    coeff_mask[Nf+3:end] .= triu(trues(Nf, Nf))[:]
end

white_noise_args = Dict(:loc=>0.0, :sig=>std(true_img), :Nsam=>1000, :norm=>false, :smooth=>false, :smoothval=>0.8) #Iso #only if you're using noise based covar
optim_settings = Dict([("iterations", 1000), ("norm", false), ("minmethod", ConjugateGradient(linesearch=LineSearches.BackTracking()))])
recon_settings = Dict([("target_type", "ground_truth"), ("covar_type", "sfd_dbn"), ("log", logbool), ("GaussianLoss", false), ("iteration_init", randn(Nx, Nx)), ("TransformedGaussianLoss", false), ("RandomInit", true), ("Invcov_matrix", invcovtype), ("epsvalue", 1e-8),
  ("optim_settings", optim_settings), ("white_noise_args", white_noise_args)]) #Add constraints

recon_settings["datafile"] = datfile

regs_true = DHC_compute_wrapper(log.(true_img), filter_hash; norm=false, dhc_args...)[coeff_mask]
s2mean, s2icov = meancov_generator(true_img, filter_hash, dhc_args, coeff_mask, recon_settings, safety=nothing)
#Check
#s2mean - s_true

if recon_settings["log"] & dhc_args[:apodize]
    lval = std(apodizer(log.(true_img))).^(-2)
elseif !recon_settings["log"] & dhc_args[:apodize]
    s_init = DHC_compute_wrapper(apodizer(init), filter_hash, norm=false; dhc_args...)[coeff_mask]
    l1init = ( 0.5 .* (s_init - s2mean)' * s2icov * (s_init - s2mean))[1]
    println("L1init", l1init)
    est1 = std(apodizer(true_img)).^(-2)
    wn_exp = mean(wind_2d(Nx).^2)*Nx*Nx*(loaddf["std"])^2
    lval = minimum([est1, 0.01*(l1init - 1.0)/(wn_exp)])
elseif recon_settings["log"] & !dhc_args[:apodize]
    lval = std(log.(true_img)).^(-2)
else
    s_init = DHC_compute_wrapper(init, filter_hash, norm=false; dhc_args...)[coeff_mask]
    l1init = ( 0.5 .* (s_init - s2mean)' * s2icov * (s_init - s2mean))[1]
    println("L1init", l1init)
    est1 = std(true_img).^(-2)
    wn_exp = Nx*Nx*(loaddf["std"])^2
    lval = minimum([est1, 0.01*(l1init - 1.0)/(wn_exp)])
end

recon_settings["lambda"] = lval
println("Regularizer Lambda=", round(lval, sigdigits=3))

recon_settings["s_targ_mean"] = s2mean
recon_settings["s_invcov"] = s2icov
recon_settings["safemode"] = false
recon_settings["fname_save"] = fname_save * ".jld2"
recon_settings["optim_settings"] = optim_settings
resobj, recon_img = reconstruction_wrapper(true_img, imfilter(init, Kernel.gaussian(0.8)), filter_hash, dhc_args, coeff_mask, recon_settings)

plot(heatmap(true_img), heatmap(recon_settings["iteration_init"]), heatmap(recon_img))

DHC_compute_wrapper(recon_img, filter_hash, norm=false; dhc_args...) .- DHC_compute_wrapper(true_img, filter_hash, norm=false; dhc_args...)

p= plot([t.value for t in Optim.trace(resobj)])
plot!(title="Loss: S20 | Targ = S(SFD) | Cov = SFD Covariance", xlabel = "No. Iterations", ylabel = "Loss")
savefig(p, fname_save * "_trace.png")

Visualization.plot_diffscales([apodizer(true_img), apodizer(init), apodizer(recon_img), apodizer(recon_img) - apodizer(true_img)], ["GT", "Init", "Reconstruction", "Residual"], fname=fname_save*".png")
Visualization.plot_synth_QA(apodizer(true_img), apodizer(init), apodizer(recon_img), fname=fname_save*"_6p.png")

println("Input Data", datfile)
println("Output File", fname_save)
println("Mean Abs Res: Init-True = ", mean(abs.(init - true_img)), " Recon-True = ", mean(abs.(recon_img - true_img)))
println("Mean Abs Frac Res", mean(abs.((init - true_img)./true_img)), " Recon-True=", mean(abs.((recon_img - true_img)./true_img)))
println("Mean Abs Res: Init-True = ", mean((init - true_img).^2).^0.5, " Recon-True = ", mean((recon_img - true_img).^2).^0.5)

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
@testset "1636.sort-array-by-increasing-frequency.jl" begin
    @test frequency_sort([-1,1,-6,4,5,-6,1,4,1]) == [5,-1,4,4,-6,-6,1,1,1]
    @test frequency_sort([2,3,1,3,2]) == [1,3,3,2,2]
    @test frequency_sort([1,1,2,2,2,3]) == [3,1,1,2,2,2]
end


text:
using OrdinaryDiffEq, Test

function test_ode(u, p, t)
    [p[1] - (1 - p[1])*u[1]]
end

function test_ode(du, u, p, t)
    du[1] = p[1] - (1 - p[1])*u[1]
    return nothing
end

test_solution(t) = t <= 5 ? t : 5. * exp(-(t-5))

tspan = (0.,10.)
testtimes = range(tspan[1], stop=tspan[2], length=1001)
pullback_condition(u, t, i) = t - 5
pullback_affect!(i) = i.p[1] = abs(1 - i.p[1])
cb = ContinuousCallback(pullback_condition, pullback_affect!)

# DPRKN6 is left out because second order ODE
algs = [Tsit5, Rosenbrock23, DP5, DP8, SSPRK43, SSPRK432, OwrenZen3, SSPRK22, SSPRK33,
        OwrenZen4, OwrenZen5, Rosenbrock32, Rodas5, Rodas4, Rodas42] ## Works for these
bad_algs = []
lazy_alg = [BS5, Vern6, Vern7, Vern8, Vern9]

nonstandard_interp_algs = union(algs,bad_algs,lazy_alg)

passed = fill(false, 2length(nonstandard_interp_algs))

cur_itr = 0
for inplace in [false,true], alg in algs
    global cur_itr
    prob = ODEProblem{inplace}(test_ode, [0.], tspan, [1.])
    sol = solve(prob, alg(); callback=cb,dt=0.0013)
    pass = all(isapprox(sol(t)[1], test_solution(t); atol=0.05) for t in testtimes)
    x = maximum(sol(t)[1] - test_solution(t) for t in testtimes)
    @show x, x > 0.05, pass
    cur_itr += 1
    passed[cur_itr] = pass
end

for inplace in [false,true], alg in lazy_alg
    global cur_itr
    prob = ODEProblem{inplace}(test_ode, [0.], tspan, [1.])
    sol = solve(prob, alg(); callback=cb,dt=0.0013)
    fail = all(isapprox(sol(t)[1], test_solution(t); atol=0.05) for t in testtimes)

    prob = ODEProblem{inplace}(test_ode, [0.], tspan, [1.])
    sol = solve(prob, alg(lazy=false); callback=cb,dt=0.0013)
    pass = all(isapprox(sol(t)[1], test_solution(t); atol=0.05) for t in testtimes)

    cur_itr += 1
    passed[cur_itr] = pass && !fail
end

any(.!(passed)) && @warn("The following algorithms failed the continuous callback test: $(vcat(algs,algs,lazy_alg,lazy_alg)[.!(passed)])")
@test all(passed)


text:
@testset "report invalid builtin call" begin
    analyzer, = report_call((Int, Type{Int}, Any)) do a, b, c
        isa(a, b, c)
    end
    @test length(get_reports(analyzer)) === 1
    report = first(get_reports(analyzer))
    @test report isa UnimplementedBuiltinCallErrorReport &&
        widenconst.(report.argtypes) == [Int, Type{Int}, Any]

    @testset "constant propagation" begin
        m = gen_virtual_module()
        Core.eval(m, quote
            struct T
                v
            end
            access_field(t, sym) = getfield(t, sym)
        end)

        analyzer, = Core.eval(m, quote
            $report_call(t->access_field(t,:v), (T,))
        end)
        @test isempty(get_reports(analyzer))

        analyzer, = Core.eval(m, quote
            $report_call(t->access_field(t,:w), (T,))
        end)
        @test length(get_reports(analyzer)) === 1
        er = first(get_reports(analyzer))
        @test er isa NoFieldErrorReport
        @test er.typ === m.T
        @test er.name === :w

        analyzer, = Core.eval(m, quote
            $report_call(t->access_field(t,:v), (T,))
        end)
        @test isempty(get_reports(analyzer))
    end
end

@testset "malformed getfield" begin
    let
        # shouldn't error
        analyzer, = report_call((Any,)) do a
            getfield(a)
        end
        @test length(get_reports(analyzer)) == 1
        @test first(get_reports(analyzer)) isa UnimplementedBuiltinCallErrorReport
    end
end

@testset "getfield with abstract global variable" begin
    # nested module access will be resolved as a direct call of `getfield`
    let
        res = @analyze_toplevel begin
            module foo

            const bar = sum

            module baz

            using ..foo

            foo.bar("julia") # -> NoMethodErrorReports

            end # module bar

            end # module foo
        end

        @test isempty(res.toplevel_error_reports)
        test_sum_over_string(res)
    end

    # this should work even if the accessed variable is not constant
    let
        res = @analyze_toplevel begin
            module foo

            bar = sum

            module baz

            using ..foo

            foo.bar("julia") # -> NoMethodErrorReports

            end # module bar

            end # module foo
        end

        @test isempty(res.toplevel_error_reports)
        test_sum_over_string(res)
    end
end

@testset "special case `return_type`" begin
    # don't report invalid method calls simulated in `return_type_tfunc`
    let
        analyzer, = report_call(()->CC.return_type(sum, Tuple{String}))
        @test isempty(get_reports(analyzer))
    end

    # report invalid call of `return_type` itself
    let
        analyzer, = report_call(()->CC.return_type(sum))
        @test length(get_reports(analyzer)) == 1
        @test isa(first(get_reports(analyzer)), InvalidReturnTypeCall)
    end

    # end to end
    let
        # this shouldn't report "no matching method found for call signature: Base.iterate(itr::DataType)"
        # , which otherwise will be caught in `abstract_cal` in `return_type_tfunc`
        analyzer, = report_call(()->Dict('a' => 1,
                                              :b => 2)
                                     )
        @test isempty(get_reports(analyzer))
    end
end


text:
##
# Christoph Aurnhammer <aurnhammer@coli.uni-saarland.de>
#
# Functions implementing lmerERP analysis
##

## FUNCTION DEFINITION
function read_data(path)
    # Load data. Transform grouping variables into pooledarrrays.
    data = transform!(DataFrame(CSV.File(path)),
        :Subject => PooledArray => :Subject,
        :ItemNum => PooledArray => :ItemNum);

    # !!!!!!!!!!!!!!!!!!!!
    # data = data[(data.Condition .== "C") .| (data.Condition .== "D"),:]
    data = data[(data.Condition .== "A"),:]

    # Computer interaction, then zscore
    data[:,:Surprisal] = log.(data.Cloze .+ 0.01)
    data[:,:inter] = zscore(data[:,:rcnoun] .* data[:,:Surprisal])

    # data[:,:logarc] = log.(normalize(data.Association_RC).+0.01)
    data[:,:logarc] = log.(data.Association_RC)
    data[:,:sqrtarc] = zscore(sqrt.(data.Association_RC))
    data[:,:squaredarc] = zscore(data.Association_RC.^2)

    data[:,:interaction] = zscore(data[:,:Association_RC] .* (data[:,:Cloze] .+ 1));
    data[:,:interaction_log] = zscore(data[:,:logarc] .* (data[:,:Surprisal]));

    data[:,:Cloze] = zscore((data[:,:Cloze] .+ 1));
    data[:,:Association_RC] = zscore(data[:,:Association_RC]);
    data[:,:Surprisal] = zscore(data[:,:Surprisal]);
    data[:,:logarc] = zscore(data[:,:logarc])

    data[:,:logrcnoun] = zscore(log.(data[:,:rcnoun]))
    data[:,:rcnoun] = zscore(data[:,:rcnoun])
    data[:,:rcverb] = zscore(data[:,:rcverb])

    # Melt data to long format. Select electrodes on the fly.
    m_vars = [:Condition, :ItemNum, :Subject, :Timestamp, :TrialNum, :Cloze, :Cloze_Balanced, :Association, :Association_RC, :Association_MC, :Association_weighted, :RT, :ACC, :interaction, :Surprisal, :rcverb, :rcnoun, :inter, :logrcnoun];
    id_vars = [:Fz, :Cz, :Pz];
    data = stack(data, id_vars, m_vars, variable_name=:Electrode, value_name=:EEG);

    # rename predictors
    rename!(data, :ItemNum => :Item, :Cloze => :cl, :Association_RC => :arc, :Association_MC => :amc, :Association => :a, :Association_weighted => :aw, :Surprisal => :srp);

    data
end

function normalize(vector)
    min_v = minimum(vector)
    max_v = maximum(vector)

    [(x - min_v)/(max_v - min_v) for x in vector]
end

function item_coef(b_array,index,n_sub,n_it,n_pred)
    collect(Iterators.flatten(fill(reshape(b_array,(n_it*n_pred))[index:n_pred:(n_it*n_pred)],n_sub)))
end

function subj_coef(b_array,index,n_sub,n_it,n_pred)
    repeat(reshape(b_array,(n_sub*n_pred))[index:n_pred:(n_sub*n_pred)], inner=(n_it))
end

function fit_models(data, form, f_contr ; parallel)
    ## Allocate output df, run models, enter numsub * numitems rows into output in one go.
    # output alloction
    num_groupings = length(f_contr);             # number of grouping factors (in ranefs)
    numpred = length(form.rhs) - num_groupings;  # number of predictors (including zval_intercept)
    pred_names = string.(form.rhs[1:numpred]);   # get predictor names
    tl = length(unique(data[:Timestamp]));       # num of unique timestamps
    el = length(unique(data[:Electrode]));       # num of unique electrodes
    su = length(unique(data[:Subject]));         # num of unique subjects
    it = length(unique(data[:Item]));            # num of unique items
    n = tl * el * su * it;                       # nrows of final data frame
    r = su * it;                                 # nrows that are written
    subj = unique(data[:Subject]);               # array of unique subjects
    items = unique(data[:Item]);                 # array of unique items
    # Preallocate output dataframe
    output_names = vcat([:Timestamp, :Electrode, :Subject, :Item], [Symbol.(string("zval_",string(x))) for x in Symbol.(pred_names)],
                        [Symbol.(string("pval_",string(x))) for x in Symbol.(pred_names)],
                        [Symbol.(string("SE_",string(x))) for x in Symbol.(pred_names)],
                        [x for x in Symbol.(pred_names)],
                        [Symbol.(string("s_",string(x))) for x in Symbol.(pred_names)],
                        [Symbol.(string("i_",string(x))) for x in Symbol.(pred_names)], :AIC, :BIC);

    ncols = (4 + numpred * 4 + numpred * num_groupings + 2);
    output = DataFrame(zeros(n, ncols));
    names!(output, output_names);
    output[:Electrode] = string.(output[:Electrode]);

    # model fitting
    counter = 0;
    startind = 1;
    for (i, t) in enumerate(unique(data[:Timestamp]))
        if parallel == true
            output = parallel_lmer(data, form, f_contr, counter, startind, tl, el, t, r, subj, items, su, it, numpred, output)
            counter = counter + 3
            startind = startind + r * 3
        else
            output = serial_lmer(data, form, f_contr, counter, startind, tl, el, t, r, subj, items, su, it, numpred, output)
            counter = counter + 3
            startind = startind + r * 3
        end
    end

    output
end

function parallel_lmer(data, form, f_contr, counter, startind, tl, el, t, r, subj, items, su, it, numpred, output)
    Threads.@threads for e in unique(data[:Electrode])
                # subset data for timestamp and electrode
                te_dt = data[(data.Timestamp .== t) .& (data.Electrode .== e),:];
                print("> ", round((counter*100) / (tl*el), digits=1), "% done. Model = ", counter, " / ", tl*el, "; Timestamp = ", t, ".", "\r")
                m = fit(MixedModel, form, te_dt, contrasts=f_contr);

                zvals = m.beta ./ m.stderror
                aic = m.objective + 2 * dof(m)
                bic = m.objective + dof(m) * log(nobs(m))

                # determine input indices
                if e == "Cz"
                    sind = startind+r
                elseif e == "Pz"
                    sind = startind+(r)*2
                else
                    sind = startind
                end

                output[sind:(sind+r-1),:] = DataFrame(vcat([fill(t,r), fill(e,r), repeat(subj,inner=it), repeat(items,outer=su)], [fill(y,r) for y in zvals], [fill(y,r) for y in m.pvalues], [fill(y,r) for y in m.stderror],
                                                    [fill(y,r) for y in m.beta], [subj_coef(m.b[2],x,su,it,numpred) for x in 1:numpred], [item_coef(m.b[1],x,su,it,numpred) for x in 1:numpred], [fill(aic,r), fill(bic, r)]))
    end

    output
end

function serial_lmer(data, form, f_contr, counter, startind, tl, el, t, r, subj, items, su, it, numpred, output)
    for (j, e) in enumerate(unique(data[:Electrode]))
            counter = counter + 1
            # subset data for timestamp and electrode
            te_dt = data[(data.Timestamp .== t) .& (data.Electrode .== e),:];

            # fit model
            print("> ", round((counter*100) / (tl*el), digits=1), "% done. Model = ", counter, " / ", tl*el, "; Timestamp = ", t, ". Electrode = ", e, ".\r")
            m = fit(MixedModel, form, te_dt, contrasts=f_contr);
            zvals = m.beta ./ m.stderror
            aic = m.objective + 2 * dof(m)
            bic = m.objective + dof(m) * log(nobs(m))

            # TO DO : use ranef names. Currently only supports fixef = ranef models.
            output[startind:(startind+r-1),:] = DataFrame(vcat([fill(t,r), fill(e,r), repeat(subj,inner=it), repeat(items,outer=su)], [fill(y,r) for y in zvals], [fill(y,r) for y in m.pvalues], [fill(y,r) for y in m.stderror],
                                                [fill(y,r) for y in m.beta], [subj_coef(m.b[2],x,su,it,numpred) for x in 1:numpred], [item_coef(m.b[1],x,su,it,numpred) for x in 1:numpred], [fill(aic,r), fill(bic, r)]))

            startind = startind + r
    end

    output
end

function combine_datasets(eeg_data, lmer_data, preds)
    # collapse eeg data to item * subject subset of non-excluded trials. Rename predictor cols to z*. (zscored predictors)
    stim = combine(groupby(eeg_data, [:Subject, :Item, :Condition, :TrialNum]), [x => mean => string("z_", x) for x in preds]);
    lmer_data = innerjoin(lmer_data, stim, on = [:Subject, :Item]);
    lmer_data = innerjoin(lmer_data, eeg_data[:,[:Timestamp, :Electrode, :Subject, :Item, :EEG]], on = [:Timestamp, :Electrode, :Subject, :Item]);

    lmer_data
end

function compute_waveforms(lmer_data, preds)
    # Exclude interecept in combinations
    if preds[1] == Symbol("1")
        combi = [[x] for x in preds[2:length(preds)]];
    else
        combi = [[x] for x in preds];
    end

    # TODO: improve to include more than two combis
    # Genereate predictor combinations
    combi_cp = combi
    global excludelist = [:plc]
    global combi
    for x in combi_cp
        for y in combi_cp
            if (x != y) & (sum([[x, y] == z for z in excludelist]) == 0)
                # TODO: avoid global
                global combi = vcat(combi, [[x, y]]);
                global excludelist = vcat(excludelist, [[y, x]]);
            end
        end
    end

    # vcat intercept with combination
    if preds[1] == Symbol("1")
        combi = vcat([[preds[1]]], combi);
        combi =  vcat(combi, [[[x] for x in preds[2:length(preds)]]]);
    else
        combi = vcat(combi, [[[x] for x in preds[2:length(preds)]]])
    end

    for x in combi
        if x[1] == Symbol("1") # intercept
            lmer_data[:est_1] = lmer_data[Symbol("1")] .+ lmer_data[:s_1]  .+ lmer_data[:i_1]
        elseif typeof(x[1]) == Symbol # intercept + single pred
            lmer_data[Symbol(string("est_"), x[1])] = lmer_data[:est_1] + (lmer_data[x[1]] .+ lmer_data[Symbol("s_", x[1])] .+ lmer_data[Symbol("i_", x[1])]) .* lmer_data[Symbol("z_", x[1])];
        elseif typeof(x[1]) == Array{Symbol,1} # intercept + combination of preds
            # start with intercept
            temp_data = lmer_data[:est_1]
            colname = "est_"
            # for each element within combination
            for y in x
                colname = string(colname, y[1])
                temp_data = temp_data .+ (lmer_data[Symbol(y[1])] .+ lmer_data[Symbol("s_", y[1])] .+ lmer_data[Symbol("i_", y[1])]) .* lmer_data[Symbol("z_", y[1])]
            end
            lmer_data[Symbol(colname)] = temp_data
        end
    end

    lmer_data
end

function compute_residuals(lmer_data, preds)
    # Residuals: observed - estimated
    for x in preds
        lmer_data[Symbol(string("res_", string(x[5:lastindex(x)])))] = lmer_data[:EEG] .- lmer_data[Symbol(x)]
    end

    lmer_data
end

function compute_coefs(lmer_data, preds)
    # Compute sum of intercept and coefs (fixef + subject ranef + item ranef)
    for x in preds
        if x == Symbol("1")
            lmer_data[:coef_1] = lmer_data[Symbol("1")] .+ lmer_data[:s_1] .+ lmer_data[:i_1]
        else
            lmer_data[Symbol(string("coef_intercept_", string(x)))] = lmer_data[:coef_1] .+ lmer_data[x] .+ lmer_data[Symbol(string("s_", string(x)))] .+ lmer_data[Symbol(string("i_", string(x)))]
        end
    end

    lmer_data
end


text:
using NumericalContinuation
using NumericalContinuationShootingProblems
using StaticArrays
using Test

NC = NumericalContinuation

@testset "Shooting problems" begin
    hopf = (out, u, p, t) -> (out[1] = p[1]*u[1] - u[2] + p[2]*u[1]*(u[1]^2 + u[2]^2);
                              out[2] = u[1] + p[1]*u[2] + p[2]*u[2]*(u[1]^2 + u[2]^2))
    prob = NC.ProblemStructure()
    @test_throws ArgumentError add_shootingproblem!(prob, "hopf", hopf, [1.0, 0.0], [1.0, -1.0], [0, 2π], abstol=1e-8, pnames=["α", "β", "γ"])
    add_shootingproblem!(prob, "hopf", hopf, [1.0, 0.0], [1.0, -1.0], [0, 2π], abstol=1e-8)
    NC.initialize!(prob)
    out = ones(6)
    u0 = NC.get_u0(Float64, NC.get_vars(prob))
    d0 = NC.get_data(NC.get_data(prob))
    get_funcs(prob)[:embedded](out, u0, data=d0, prob=prob)
    @test isapprox(out, zero(out), atol=1e-5)

    hopf2 = (u, p, t) -> SVector(p[1]*u[1] - u[2] + p[2]*u[1]*(u[1]^2 + u[2]^2),
                                 u[1] + p[1]*u[2] + p[2]*u[2]*(u[1]^2 + u[2]^2))
    prob = NC.ProblemStructure()
    add_shootingproblem!(prob, "hopf", hopf2, SVector(1.0, 0.0), SVector(1.0, -1.0), [0, 2π], abstol=1e-8)
    NC.initialize!(prob)
    out = ones(6)
    u0 = NC.get_u0(Float64, NC.get_vars(prob))
    d0 = NC.get_data(NC.get_data(prob))
    get_funcs(prob)[:embedded](out, u0, data=d0, prob=prob)
    @test isapprox(out, zero(out), atol=1e-5)
end


text:
using RoME
using Test
# using RoMEPlotting, Distributions

import IncrementalInference: getSample

mutable struct NorthSouthPartial{T} <: FunctorSingleton
  Z::T
  partial::Tuple{Int}
  NorthSouthPartial{D}() where D = new{D}()
  NorthSouthPartial{D}(Z::D) where {D <: IIF.SamplableBelief} = new{D}(Z, (2,))
end
NorthSouthPartial(Z::D) where {D <: IIF.SamplableBelief} = NorthSouthPartial{D}(Z)

getSample(ns::NorthSouthPartial, N=1) = (reshape(rand(ns.Z, N),1,N),)

@testset "test standard multimodal range bearing factor setup..." begin

# Start with an empty graph
global N = 100
global fg = initfg(sessionname="MULTIMODAL_2D_TUTORIAL")

# Add landmarks with Bearing range measurements
addVariable!(fg, :l1, Point2, labels=[:LANDMARK;])
addVariable!(fg, :l2, Point2, labels=[:LANDMARK;])

addFactor!(fg, [:l1], Prior(MvNormal([10.0;0.0], Matrix(Diagonal([1.0;1.0].^2)))) )
addFactor!(fg, [:l2], Prior(MvNormal([30.0;0.0], Matrix(Diagonal([1.0;1.0].^2)))) )

val, = predictbelief(fg, :l1, [:l1f1;])
setVal!(fg, :l1, val)

val, = predictbelief(fg, :l2, [:l2f1;])
setVal!(fg, :l2, val)

addVariable!(fg, :x0, Pose2)
# addFactor!(fg, [:x0], Prior(MvNormal([0.0;0.0;0], Matrix(Diagonal([1.0;1.0;0.01].^2)))) )
addFactor!(fg, [:x0;], NorthSouthPartial(Normal(0,1.0)))

global p2br = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))
addFactor!(fg, [:x0; :l1; :l2], p2br, multihypo=[1.0; 0.5; 0.5])

predictbelief(fg, :x0, ls(fg, :x0))

# writeGraphPdf(fg)
end

@testset "test multimodal bearing range factors calculate pose position properly..." begin

ensureAllInitialized!(fg)

tree, smt, hist = solveTree!(fg)
# global tree = batchSolve!(fg, N=N, recursive=true)
# global tree = wipeBuildNewTree!(fg)
# global infv = [inferOverTreeR!(fg, tree, N=N) for i in 1:4]
# @test length(infv) == 4

# X0 = getVertKDE(fg, :x0)
# X0pts = getPoints(X0)
global X0pts = getVal(fg, :x0)

@test size(X0pts,1) == 3
@test size(X0pts,2) == N

# find at least two of the four possible modes
@test sum([0 < sum(-20 .< X0pts[1,:] .< 0);
           0 < sum(  0 .< X0pts[1,:] .< 20);
           0 < sum( 20 .< X0pts[1,:] .< 40);
           0 < sum( 40 .< X0pts[1,:] .< 60)]) >= 2

end

@testset "test multimodal landmark locations are computed correclty..." begin

# Start with an empty graph
global fg = initfg(sessionname="MULTIMODAL_2D_TUTORIAL")

addVariable!(fg, :x0, Pose2)
addFactor!(fg, [:x0], PriorPose2(MvNormal([0.0;0.0;0], Matrix(Diagonal([1.0;1.0;0.01].^2)))) ) # TODO IIF.Prior with IIF 0.3.9

# Add landmarks with Bearing range measurements
addVariable!(fg, :l1, Point2, labels=[:LANDMARK;])
addFactor!(fg, [:l1], PriorPose2(MvNormal([40.0;0.0], Matrix(Diagonal([1.0;1.0].^2)))) ) # TODO IIF.Prior with IIF 0.3.9

addVariable!(fg, :l2, Point2, labels=[:LANDMARK;])
addFactor!(fg, [:l2;], NorthSouthPartial(Normal(0,1.0)))
# addFactor!(fg, [:l2], PriorPose2(MvNormal([30.0;0.0], Matrix(Diagonal([1.0;1.0].^2)))) ) # TODO IIF.Prior with IIF 0.3.9

global p2br = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))
addFactor!(fg, [:x0; :l1; :l2], p2br, multihypo=[1.0; 0.5; 0.5])

# @test true

tree, smt, hist = solveTree!(fg)
# setVal!(getVert(fg, :l2), zeros(2,N))
# writeGraphPdf(fg)
# global tree = wipeBuildNewTree!(fg)
# global infv = [inferOverTreeR!(fg, tree, N=N) for i in 1:4]
# @test length(infv) == 4

# L1 = getVertKDE(fg, :l1)
global L2 = getVertKDE(fg, :l2)
# X0 = getVertKDE(fg, :x0)

global L2pts = getPoints(L2)

@test size(L2pts,1) == 2
@test size(L2pts,2) == N

# plot(L2, dims=[1;2])

global numM1 = sum(10 .< L2pts[1,:] .< 30)
global numM2 = sum(30 .< L2pts[1,:] .< 50)

@test 30 < numM1 < 70
@test 30 < numM2 < 70
@test 0.97*N <= numM1 + numM2

end


text:
include("../src/SFFM.jl")
using LinearAlgebra, Plots, JLD2, SparseArrays

include("../examples/paperNumerics/exampleModelDef.jl")

@load pwd()*"/examples/paperNumerics/dump/sims.jld2" sims

Δ = 1.6
Nodes = collect(approxBounds[1, 1]:Δ:approxBounds[1, 2]) # collect(approxBounds[1, 1]:Δ:approxBounds[1, 2])
for NBases in 1:5
# NBases = 2
Basis = "lagrange"
mesh = SFFM.MakeMesh(
    model = approxModel,
    Nodes = Nodes,
    NBases = NBases,
    Basis = Basis,
)

simprobs = SFFM.Sims2Dist(model=simModel,mesh=mesh,sims=sims,type="probability")
whichphase = 4
p2 = plot!(
    mesh.CellNodes[[1;end],:],# [1;1]*simprobs.x',
    [1;1]*simprobs.distribution[:,:,whichphase],
    label = :false,
    color = :red,
    xlims = (0,2),
    legend = :bottomright,
    title = "Phase 10",
    seriestype = :line,
    linestyle = :dot,
    markershape = :x,
    markersize = 4,
    xlabel = "x",
    ylabel = "Cumulative probability",
    windowsize = (600,400),
    grid = false,
    tickfontsize = 10,
    guidefontsize = 12,
    titlefontsize = 18,
    legendfontsize = 10,
)

All = SFFM.MakeAll(model = approxModel, mesh = mesh, approxType = "projection")
Matrices = SFFM.MakeMatrices(model=approxModel,mesh=mesh,probTransform=false)
MatricesR = SFFM.MakeMatricesR(model=approxModel,mesh=mesh)
B = SFFM.MakeB(model=approxModel,mesh=mesh,Matrices=Matrices,probTransform=false)
Dr = SFFM.MakeDR(
    Matrices=Matrices,
    MatricesR=MatricesR,
    model=approxModel,
    mesh=mesh,
    B=B,
)

Ψ = SFFM.PsiFun(D=Dr.DDict)
Ψbase = SFFM.PsiFun(D=All.D)

# construct initial condition
theNodes = mesh.CellNodes[:,convert(Int,ceil(5/Δ))]
basisValues = zeros(length(theNodes))
for n in 1:length(theNodes)
    basisValues[n] = prod(5.0.-theNodes[[1:n-1;n+1:end]])./prod(theNodes[n].-theNodes[[1:n-1;n+1:end]])
end
V = SFFM.vandermonde(NBases=NBases)
if mesh.Basis=="lagrange"
    basisValues = basisValues'*V.V*V.V'*2/Δ
else
    basisValues = basisValues./V.w
    basisValues = basisValues'*V.inv'*2/Δ
end
initpm = [
    zeros(sum(approxModel.C.<=0)) # LHS point mass
    zeros(sum(approxModel.C.>=0)) # RHS point mass
]
initprobs = zeros(Float64,mesh.NBases,mesh.NIntervals,approxModel.NPhases)
initprobs[:,convert(Int,ceil(5/Δ)),3] = basisValues #*All.Matrices.Local.V.V*All.Matrices.Local.V.V'.*2/Δ
initprobsBase = zeros(Float64,mesh.NBases,mesh.NIntervals,approxModel.NPhases)
initprobsBase[:,convert(Int,ceil(5/Δ)),3] = basisValues.*V.w'/2*Δ
# initprobs[1,convert(Int,ceil(5/Δ)),3] = sqrt(2)
initdist = (
    pm = initpm,
    distribution = initprobs,
    x = mesh.CellNodes,
    type = "cumulative"
) # convert to a distribution object so we can apply Dist2Coeffs
# convert to Coeffs α in the DG context
x0 = [0;0;initdist.distribution[:];0;0]' # SFFM.Dist2Coeffs(model = approxModel, mesh = mesh, Distn = initdist) # #
x0Base = [0;0;initprobsBase[:];0;0]' # SFFM.Dist2Coeffs(model = approxModel, mesh = mesh, Distn = initdist) # #
# the initial condition on Ψ is restricted to + states so find the + states
plusIdx = [
    mesh.Fil["p+"];
    repeat(mesh.Fil["+"]', mesh.NBases, 1)[:];
    mesh.Fil["q+"];
]
# get the elements of x0 in + states only
x0 = x0[plusIdx]'
x0Base = x0Base[plusIdx]'
# check that it is equal to 1 (or at least close)
println(sum(x0))

# compute x-distribution at the time when Y returns to 0
w = x0*Ψ

wbase = x0Base*Ψbase

# this can occur in - states only, so find the - states
minusIdx = [
    mesh.Fil["p-"];
    repeat(mesh.Fil["-"]', mesh.NBases, 1)[:];
    mesh.Fil["q-"];
]
# then map to the whole state space for plotting
z = zeros(
    Float64,
    mesh.NBases*mesh.NIntervals * approxModel.NPhases +
        sum(approxModel.C.<=0) + sum(approxModel.C.>=0)
)
zbase = copy(z)
z[minusIdx] = w
zbase[minusIdx] = wbase
# check that it is equal to 1
println(sum(z))
println(sum(zbase))

# convert to a distribution object for plotting
DGProbs = SFFM.Coeffs2Dist(
    model = approxModel,
    mesh = mesh,
    Coeffs = z,
    type="probability",
)
DGProbsBase = SFFM.Coeffs2Dist(
    model = approxModel,
    mesh = mesh,
    Coeffs = zbase,
    type="probability",
)
display(DGProbsBase.pm)
display(DGProbs.pm)
display(simprobs.pm)
plot!(p2, DGProbsBase.x[:], DGProbsBase.distribution[:,:,whichphase][:],label = :false,#"DG: N_k = "*string(NBases),
    color = :blue,
    xlims = (0,2),
    seriestype = :scatter
    )
1
# plot them
if Basis=="lagrange"
    if mesh.NBases==1
        p2 = plot!(p2,
            DGProbs.x,
            [1;1]*reshape(z[3:end-2], mesh.NBases, mesh.NIntervals, approxModel.NPhases)[:,:,whichphase], #DGProbs.distribution[:,:,whichphase],
            label = :false,#"DG: N_k = "*string(NBases),
            color = :blue,
            xlims = (0,2),
            seriestype = :line,
            linestyle = :dash,
        )
        p2 = plot!(p2,
            DGProbsBase.x,
            [1;1]*reshape(zbase[3:end-2], mesh.NBases, mesh.NIntervals, approxModel.NPhases)[:,:,whichphase], #DGProbsBase.distribution[:,:,whichphase],
            label = :false,#"DG: N_k = "*string(NBases),
            color = :black,
            xlims = (0,2),
            seriestype = :line,
            linestyle = :dot,
        )
    else
        p2 = plot!(p2,
            mesh.CellNodes[[1;end],:],
            [1;1]*0.5*Δ*V.w'*reshape(z[3:end-2], mesh.NBases, mesh.NIntervals, approxModel.NPhases)[:,:,whichphase], #DGProbs.distribution[:,:,whichphase],
            label = :false,#"DG: N_k = "*string(NBases),
            color = :blue,
            xlims = (0,2),
            seriestype = :line,
            linestyle = :dash,
        )
        p2 = plot!(p2,
            mesh.CellNodes[[1;end],:],
            [1;1]*ones(1,NBases)*reshape(zbase[3:end-2], mesh.NBases, mesh.NIntervals, approxModel.NPhases)[:,:,whichphase], #DGProbs.distribution[:,:,whichphase],
            label = :false,#"DG: N_k = "*string(NBases),
            color = :black,
            xlims = (0,2),
            seriestype = :line,
            linestyle = :dot,
        )
    end
else
    p2 = plot!(p2,
        [mesh.Nodes[1:end-1]';mesh.Nodes[2:end]'],
        [1;1]*DGProbs.distribution[:,:,whichphase],
        label = :false,#"DG: N_k = "*string(NBases),
        color = :blue,
        xlims = (0,2),
        seriestype = :line,
        linestyle = :dash,
    )
    p2 = plot!(p2,
        [mesh.Nodes[1:end-1]';mesh.Nodes[2:end]'],
        [1;1]*DGProbsBase.distribution[:,:,whichphase],
        label = :false,#"DG: N_k = "*string(NBases),
        color = :black,
        xlims = (0,2),
        seriestype = :line,
        linestyle = :dot,
    )
end
d = zeros(1, mesh.NIntervals, approxModel.NPhases)
for i in 1:approxModel.NPhases
    d[:,:,i] = 0.5*Δ*V.w'*reshape(z[3:end-2], mesh.NBases, mesh.NIntervals, approxModel.NPhases)[:,:,i]
end
DGProbs = (pm = DGProbs.pm,
    distribution = d,
    x = DGProbs.x,
    type = "probability")
eR = SFFM.starSeminorm(d1 = DGProbs, d2 = simprobs)
e = SFFM.starSeminorm(d1 = DGProbsBase, d2 = simprobs)
ediff = SFFM.starSeminorm(d1 = DGProbsBase, d2 = DGProbs)
println("DR: ",eR)
println("D: ",e)
println("diff: ",ediff)
end
1
# p = SFFM.PlotSFM(model=approxModel,mesh=mesh,Dist=DGProbs)
# SFFM.PlotSFM!(p;model=approxModel,mesh=mesh,Dist=simprobs,color=2)
# for sp in 1:4
#     plot!(subplot=sp, xlims=(0,2))
# end
# display(plot!())
# plot(z)

## simpler set up to start
## Model 1
T = [-1.0 1.0; 1.0 -1.0]
C = [1.0; -1.2]
# r = (
#     r = function (x)
#         [abs.((x .<= 1)-(x .> 1)) abs.(((x .> 1))-2*(x .<= 1).*(x .> 0)-(x.==0))]#[ones(size(x)) ones(size(x))]#
#     end,
#     R = function (x)
#         [x.*abs.(((x .<= 1)-(x .> 1))) x.*abs.((((x .> 1))-2*(x .<= 1).*(x .> 0)-(x.==0)))]#[x x]
#     end,
# )

# r = (
#     r = function (x)
#         [1.5*(x.>-1) -2*(x.>-1)]#[ones(size(x)) ones(size(x))]#
#     end,
#     R = function (x)
#         [1.5*(x.>-1).*x -2*(x.>-1).*x]
#     end,
# )

## Model 3
# T = [-1.0 1.0; 1.0 -1.0]*4
# C = [1.0; -2.0]
# r = (
#     r = function (x)
#         [-2*(cos.(x).+1.05) 1*(x.>-1)]
#     end,
#     R = function (x)
#         [-2*(sin.(x).+1.05.*x) 1*(x.>-1).*x]
#     end,
# )
let
T = [-1.0 1.0; 1.0 -1.0]*0.5
C = [1.0; -2]
r = (
    r = function (x)
        [(cos.(x).+1.05)./(2.05) -3*(x.>-1)]
    end,
    R = function (x)
        [(sin.(x).+1.05.*x)./(2.05) -3*(x.>-1).*x]
    end,
)

Bounds = [0 10; -Inf Inf]
model = SFFM.Model(T = T, C = C, r = r, Bounds = Bounds)

## Define mesh

Δ = 1
Nodes = collect(Bounds[1,1]:Δ:Bounds[1,2])
# NBases = 3
for NBases = 1:5
Basis = "lagrange"
mesh = SFFM.MakeMesh(model = model, Nodes = Nodes, NBases = NBases, Basis=Basis)

## matrices
All = SFFM.MakeAll(model = model, mesh = mesh, approxType = "projection")
Matrices = All.Matrices
MatricesR = SFFM.MakeMatricesR(model=model,mesh=mesh)
Matrices2 = SFFM.MakeMatrices2(model=model,mesh=mesh)#,probTransform=false)
B = SFFM.MakeB(model=model,mesh=mesh,Matrices=Matrices2,probTransform=false)
Dr = SFFM.MakeDR(
    Matrices=Matrices2,
    MatricesR=MatricesR,
    model=model,
    mesh=mesh,
    B=B,
)
## sims for gound truth
x₀ = 1.21
NSim = 50000
IC = (φ = ones(Int, NSim), X = x₀ .* ones(NSim), Y = zeros(NSim))
y = 10
Ψr = SFFM.PsiFun(D=Dr.DDict,s=0)
Ψ = SFFM.PsiFun(D=All.D,s=0)
simBounds = [0 Inf; -Inf Inf]
simModel = SFFM.Model(T = T, C = C, r = r, Bounds = simBounds)

@time sims =
    SFFM.SimSFFM(model = simModel, StoppingTime = SFFM.FirstExitY(u = 0, v = Inf), InitCondition = IC)

simdensity = SFFM.Sims2Dist(model=simModel,mesh=mesh,sims=sims,type="density")
# @time sims =
    # SFFM.SimSFFM(model = model, StoppingTime = SFFM.FirstExitY(u = -Inf, v = y), InitCondition = IC)

## DG stationary dist
# construct initial condition
theNodes = mesh.CellNodes[:,convert(Int,ceil(x₀/Δ))]
basisValues = zeros(length(theNodes))
for n in 1:length(theNodes)
    basisValues[n] = prod(x₀.-theNodes[[1:n-1;n+1:end]])./prod(theNodes[n].-theNodes[[1:n-1;n+1:end]])
end
V = SFFM.vandermonde(NBases=NBases)
if mesh.Basis=="lagrange"
    basisValues = basisValues'*V.V*V.V'*2/Δ
else
    basisValues = basisValues./V.w
    basisValues = basisValues'*V.inv'*2/Δ
end
initpm = [
    zeros(sum(model.C.<=0)) # LHS point mass
    zeros(sum(model.C.>=0)) # RHS point mass
]
initprobs = zeros(Float64,mesh.NBases,mesh.NIntervals,model.NPhases)
initprobs[:,convert(Int,max(1,ceil(x₀/Δ))),1] = basisValues #*All.Matrices.Local.V.V*All.Matrices.Local.V.V'.*2/Δ
initprobsBase = zeros(Float64,mesh.NBases,mesh.NIntervals,model.NPhases)
initprobsBase[:,convert(Int,max(1,ceil(x₀/Δ))),1] = basisValues.*V.w'/2*Δ
# initprobs[1,convert(Int,ceil(5/Δ)),3] = sqrt(2)
initdist = (
    pm = initpm,
    distribution = initprobs,
    x = mesh.CellNodes,
    type = "cumulative"
) # convert to a distribution object so we can apply Dist2Coeffs
# convert to Coeffs α in the DG context
x0 = [0;initdist.distribution[:];0]' # SFFM.Dist2Coeffs(model = approxModel, mesh = mesh, Distn = initdist) # #
x0Base = [0;initprobsBase[:];0]' # SFFM.Dist2Coeffs(model = approxModel, mesh = mesh, Distn = initdist) # #
# the initial condition on Ψ is restricted to + states so find the + states
plusIdx = [
    mesh.Fil["p+"];
    repeat(mesh.Fil["+"]', mesh.NBases, 1)[:];
    mesh.Fil["q+"];
]
# get the elements of x0 in + states only
x0 = x0[plusIdx]'
x0Base = x0Base[plusIdx]'
# check that it is equal to 1 (or at least close)
println(sum(x0))

# compute x-distribution at the time when Y returns to 0
w = x0*Ψr

wbase = x0Base*Ψ

# this can occur in - states only, so find the - states
minusIdx = [
    mesh.Fil["p-"];
    repeat(mesh.Fil["-"]', mesh.NBases, 1)[:];
    mesh.Fil["q-"];
]
# then map to the whole state space for plotting
z = zeros(
    Float64,
    mesh.NBases*mesh.NIntervals * model.NPhases +
        sum(model.C.<=0) + sum(model.C.>=0)
)
zbase = copy(z)
z[minusIdx] = w
zbase[minusIdx] = wbase
# check that it is equal to 1
println(sum(z))
println(sum(zbase))

DRdensity = SFFM.Coeffs2Dist(model=model,mesh=mesh,Coeffs=z,type="density")
Ddensity = SFFM.Coeffs2Dist(model=model,mesh=mesh,Coeffs=zbase,type="density")

# p = SFFM.PlotSFM(model=model,mesh=mesh,Dist=simdensity)
# p = SFFM.PlotSFM!(p;model=model,mesh=mesh,Dist=DRdensity,color=:red)
# p = SFFM.PlotSFM!(p;model=model,mesh=mesh,Dist=Ddensity,color=:black)

simprobs = SFFM.Sims2Dist(model=model,mesh=mesh,sims=sims,type="probability")
DRprobs = SFFM.Coeffs2Dist(model=model,mesh=mesh,Coeffs=z,type="probability")
Dprobs = SFFM.Coeffs2Dist(model=model,mesh=mesh,Coeffs=zbase,type="probability")

display(simprobs.pm)
display(DRprobs.pm)
display(Dprobs.pm)
d = zeros(1, mesh.NIntervals, model.NPhases)
for i in 1:model.NPhases
    d[:,:,i] = 0.5*Δ*V.w'*reshape(z[2:end-1], mesh.NBases, mesh.NIntervals, model.NPhases)[:,:,i]
end
DGProbs = (pm = DRprobs.pm,
    distribution = d,
    x = DRprobs.x,
    type = "probability")
eR = SFFM.starSeminorm(d1 = DGProbs, d2 = simprobs)
e = SFFM.starSeminorm(d1 = Dprobs, d2 = simprobs)
ediff = SFFM.starSeminorm(d1 = DGProbs, d2 = Dprobs)
# plot(DRprobs.x[:], DRprobs.distribution[:,:,2][:],label="DR")
xvals = [mesh.Nodes[1:end-1]'; mesh.Nodes[2:end]']
plot(xvals, [1;1]*Dprobs.distribution[:,:,2],label=false, color=:black,
    seriestype = :line,
)
plot!(xvals, [1;1]*simprobs.distribution[:,:,2],label=false, color=:red)

if Basis=="lagrange"
    if mesh.NBases==1
        p2 = plot!(
            xvals,
            [1;1]*0.5*Δ*2*reshape(z[2:end-1], mesh.NBases, mesh.NIntervals, model.NPhases)[:,:,2], #DGProbs.distribution[:,:,2],
            label = :false,#"DG: N_k = "*string(NBases),
            color = :blue,
            xlims = (0,10),
            seriestype = :line,
            linestyle = :dash,
        )
    else
        p2 = plot!(
            xvals,
            [1;1]*0.5*Δ*V.w'*reshape(z[2:end-1], mesh.NBases, mesh.NIntervals, model.NPhases)[:,:,2], #DGProbs.distribution[:,:,2],
            label = :false,#"DG: N_k = "*string(NBases),
            color = :blue,
            xlims = (0,10),
            seriestype = :line,
            linestyle = :dash,
        )
    end
else
    plot!(xvals, [1;1]*DRprobs.distribution[:,:,2],label=false, color=:blue)
    # p2 = plot!(
    #     DGProbs.x,
    #     DGProbs.distribution[:,:,2],
    #     label = :false,#"DG: N_k = "*string(NBases),
    #     color = :blue,
    #     xlims = (0,10),
    #     seriestype = :line,
    #     linestyle = :dash,
    # )
end

println("DR: ", eR)
println("D: ", e)
println("diff: ", ediff)

plot!([0],[0], color = :blue, label = "DR")
plot!([0],[0], color = :black,label = "D")
plot!([0],[0], color = :red,label = "SIM")
end

end

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
function extract_kwargs(kwargs, keys)
    return (k => v for (k, v) in pairs(kwargs) if k in keys)
end

extract_fields(x) = Dict(String(name) => getfield(x, name) for name in fieldnames(typeof(x)))

# Trailing directory separator matters for `basename`, so let's strip it to
# avoid confusion
strip_backslash(path::AbstractString) = isdirpath(path) ? dirname(path) : path


text:
using ESeriesRounding
using Documenter

DocMeta.setdocmeta!(ESeriesRounding, :DocTestSetup, :(using ESeriesRounding); recursive=true)

makedocs(;
    modules=[ESeriesRounding],
    authors="KronosTheLate",
    repo="https://github.com/KronosTheLate/ESeriesRounding.jl/blob/{commit}{path}#{line}",
    sitename="ESeriesRounding.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://KronosTheLate.github.io/ESeriesRounding.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/KronosTheLate/ESeriesRounding.jl",
)


text:
#####
##### Viscous fluxes
#####

@inline viscous_flux_ux(i, j, k, grid, ν_ccc, u) = ν_ccc * ℑxᶜᵃᵃ(i, j, k, grid, Axᵃᵃᶜ) * ∂xᶜᵃᵃ(i, j, k, grid, u)
@inline viscous_flux_uy(i, j, k, grid, ν_ffc, u) = ν_ffc * ℑyᵃᶠᵃ(i, j, k, grid, Ayᵃᵃᶜ) * ∂yᵃᶠᵃ(i, j, k, grid, u)
@inline viscous_flux_uz(i, j, k, grid, ν_fcf, u) = ν_fcf * ℑzᵃᵃᶠ(i, j, k, grid, Azᵃᵃᵃ) * ∂zᵃᵃᶠ(i, j, k, grid, u)

@inline viscous_flux_vx(i, j, k, grid, ν_ffc, v) = ν_ffc * ℑxᶠᵃᵃ(i, j, k, grid, Axᵃᵃᶜ) * ∂xᶠᵃᵃ(i, j, k, grid, v)
@inline viscous_flux_vy(i, j, k, grid, ν_ccc, v) = ν_ccc * ℑyᵃᶜᵃ(i, j, k, grid, Ayᵃᵃᶜ) * ∂yᵃᶜᵃ(i, j, k, grid, v)
@inline viscous_flux_vz(i, j, k, grid, ν_cff, v) = ν_cff * ℑzᵃᵃᶠ(i, j, k, grid, Azᵃᵃᵃ) * ∂zᵃᵃᶠ(i, j, k, grid, v)

@inline viscous_flux_wx(i, j, k, grid, ν_fcf, w) = ν_fcf * ℑxᶠᵃᵃ(i, j, k, grid, Axᵃᵃᶜ) * ∂xᶠᵃᵃ(i, j, k, grid, w)
@inline viscous_flux_wy(i, j, k, grid, ν_cff, w) = ν_cff * ℑyᵃᶠᵃ(i, j, k, grid, Ayᵃᵃᶜ) * ∂yᵃᶠᵃ(i, j, k, grid, w)
@inline viscous_flux_wz(i, j, k, grid, ν_ccc, w) = ν_ccc * ℑzᵃᵃᶜ(i, j, k, grid, Azᵃᵃᵃ) * ∂zᵃᵃᶜ(i, j, k, grid, w)

#####
##### Viscous dissipation operators
#####

"""
    ∂ⱼνᵢⱼ∂ᵢu(i, j, k, grid, νˣ, νʸ, νᶻ, u)

Calculates viscous dissipation for the u-velocity via

    1/V * [δxᶠᵃᵃ(νˣ * ℑxᶜᵃᵃ(Ax) * δxᶜᵃᵃ(u)) + δyᵃᶜᵃ(νʸ * ℑyᵃᶠᵃ(Ay) * δyᵃᶠᵃ(u)) + δzᵃᵃᶜ(νᶻ * ℑzᵃᵃᶠ(Az) * δzᵃᵃᶠ(u))]

which will end up at the location `fcc`.
"""
@inline function ∂ⱼνᵢⱼ∂ᵢu(i, j, k, grid, νˣ, νʸ, νᶻ, u)
    return 1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, viscous_flux_ux, νˣ, u) +
                                    δyᵃᶜᵃ(i, j, k, grid, viscous_flux_uy, νʸ, u) +
                                    δzᵃᵃᶜ(i, j, k, grid, viscous_flux_uz, νᶻ, u))
end

"""
    ∂ⱼνᵢⱼ∂ᵢu(i, j, k, grid, νˣ, νʸ, νᶻ, v)

Calculates viscous dissipation for the v-velocity via

    1/Vᵛ * [δxᶜᵃᵃ(νˣ * ℑxᶠᵃᵃ(Ax) * δxᶠᵃᵃ(v)) + δyᵃᶠᵃ(νʸ * ℑyᵃᶜᵃ(Ay) * δyᵃᶜᵃ(v)) + δzᵃᵃᶜ(νᶻ * ℑzᵃᵃᶠ(Az) * δzᵃᵃᶠ(v))]

which will end up at the location `cfc`.
"""
@inline function ∂ⱼνᵢⱼ∂ᵢv(i, j, k, grid, νˣ, νʸ, νᶻ, v)
    return 1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, viscous_flux_vx, νˣ, v) +
                                    δyᵃᶠᵃ(i, j, k, grid, viscous_flux_vy, νʸ, v) +
                                    δzᵃᵃᶜ(i, j, k, grid, viscous_flux_vz, νᶻ, v))
end

"""
    ∂ⱼνᵢⱼ∂ᵢu(i, j, k, grid::Grid, νˣ, νʸ, νᶻ, w)

Calculates viscous dissipation for the w-velocity via

    1/Vʷ * [δxᶜᵃᵃ(νˣ * ℑxᶠᵃᵃ(Ax) * δxᶠᵃᵃ(w)) + δyᵃᶜᵃ(νʸ * ℑyᵃᶠᵃ(Ay) * δyᵃᶠᵃ(w)) + δzᵃᵃᶠ(νᶻ * ℑzᵃᵃᶜ(Az) * δzᵃᵃᶜ(w))]

which will end up at the location `ccf`.
"""
@inline function ∂ⱼνᵢⱼ∂ᵢw(i, j, k, grid, νˣ, νʸ, νᶻ, w)
    return 1/Vᵃᵃᶠ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, viscous_flux_wx, νˣ, w) +
                                    δyᵃᶜᵃ(i, j, k, grid, viscous_flux_wy, νʸ, w) +
                                    δzᵃᵃᶠ(i, j, k, grid, viscous_flux_wz, νᶻ, w))
end

#####
##### Viscous dissipation for isotropic viscosity
#####

@inline ∂ⱼνᵢⱼ∂ᵢu(i, j, k, grid, ν, u) = ∂ⱼνᵢⱼ∂ᵢu(i, j, k, grid, ν, ν, ν, u)
@inline ∂ⱼνᵢⱼ∂ᵢv(i, j, k, grid, ν, v) = ∂ⱼνᵢⱼ∂ᵢv(i, j, k, grid, ν, ν, ν, v)
@inline ∂ⱼνᵢⱼ∂ᵢw(i, j, k, grid, ν, w) = ∂ⱼνᵢⱼ∂ᵢw(i, j, k, grid, ν, ν, ν, w)


text:
function play_slot_machine(p::Player)
    bet = 0
    print_dict(p, "slot_machine_0")

    while true
        s = readline(p.socket)
        if s == ""
            p.status = reception
            return
        elseif s == "50" && p.balance >= 50
            bet = 50
            break
        elseif s == "10" && p.balance >= 10
            bet = 10
            break
        elseif s == "5"
            if p.balance < 5
                print_dict(p, "slot_machine_1")
                bet = 1
            else
                bet = 5
            end
            break
        else
            print_dict(p, "slot_machine_2")
        end
    end

    print_dict(p, "slot_machine_3")

    if rand(1:10) == 10
        print_dict(p, "slot_machine_4")
        p.balance += bet
    else
        print_dict(p, "slot_machine_5")
        p.balance -= bet
    end

    if p.balance < 0
        print_dict(p, "debt_0")
        p.status = reception
    end
end


text:
include("material.jl")

struct Particle
    A::Int  # Mass Number
    Z::Int  # Element number
    material::Material  # Target material
end

function stoppingpower(p::Particle, E)
    Erel = E + p.mev
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
"""
    PostOperator

Abstract supertype of all post operators.

### Notes

All post operators should provide the following methods:
```julia
init(op::PostOperator, system, options)

post(op::PostOperator, system, ...)
```
"""
abstract type PostOperator end


text:
module TestVariableRSOCtoSOC

using Test

using VecMathOptInterface
const MOI = VecMathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

function test_rotatedsoc4()
    mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    bridged_mock = MOI.Bridges.Variable.RSOCtoSOC{Float64}(mock)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [√2, 0.0, 1.0, 1.0],
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1.0],
        )
    MOI.Test.test_conic_RotatedSecondOrderCone_out_of_order(
        bridged_mock,
        MOI.Test.Config(),
    )

    ceqs = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64},
            MOI.LessThan{Float64},
        }(),
    )
    @test length(ceqs) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), ceqs[1], "c")

    var_names = ["a", "b", "c", "d"]
    MOI.set(
        mock,
        MOI.VariableName(),
        MOI.get(mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    socs = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{MOI.VectorOfVariables,MOI.SecondOrderCone}(),
    )
    @test length(socs) == 1
    MOI.set(mock, MOI.ConstraintName(), socs[1], "soc")

    s2 = √2
    s = """
    variables: a, b, c, d
    soc: [a, b, c, d] in VecMathOptInterface.SecondOrderCone(4)
    c: $s2*a <= 2.0
    maxobjective: c + d
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(mock, model, var_names, ["soc", "c"])
    var_names = ["t", "u", "x", "y"]
    MOI.set(
        bridged_mock,
        MOI.VariableName(),
        MOI.get(bridged_mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    rsocs = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorOfVariables,
            MOI.RotatedSecondOrderCone,
        }(),
    )
    @test length(rsocs) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), rsocs[1], "rsoc")

    s = """
    variables: t, u, x, y
    rsoc: [t, u, x, y] in VecMathOptInterface.RotatedSecondOrderCone(4)
    c: t + u <= 2.0
    maxobjective: x + y
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        bridged_mock,
        model,
        var_names,
        ["rsoc", "c"],
    )

    tuxy = MOI.get(bridged_mock, MOI.ListOfVariableIndices())

    message = string(
        "Cannot delete variable as it is constrained with other",
        " variables in a `MOI.VectorOfVariables`.",
    )
    for i in eachindex(tuxy)
        err = MOI.DeleteNotAllowed(tuxy[i], message)
        @test_throws err MOI.delete(bridged_mock, tuxy[i])
    end

    _test_delete_bridged_variables(
        bridged_mock,
        tuxy,
        MOI.RotatedSecondOrderCone,
        4,
        ((MOI.VectorOfVariables, MOI.SecondOrderCone, 0),),
    )
    return
end

end  # module

TestVariableRSOCtoSOC.runtests()


text:
function eigmaxreal(A::Matrix)
  if size(A,1)<= 2
    return maximum(real(eigvals(A)))
  else
    return eigs(A,nev=1,which = :LR)[1][1]
  end
end
function eigminreal(A::Matrix)
  if size(A,1)<= 2
    return minimum(real(eigvals(A)))
  else
    return eigs(A,nev=1,which = :SR)[1][1]
  end
end


text:
# Dummy Module for testing macro expansion.
module TestModule
import CachedArrays
end

# The general idea here is that we manually perform macroexpansion on expressions passed
# to `CachedArrays.@annotate`. This will do all the the things like symbol and function
# resolution.
#
# We can then use the `@expected` macro to manually specify what we expect the macro
# to expand to for testing.
#
# This is a little brittle and requires a bit of care for the following reasons:
#
# 1. Gensyms need to be replaced in a deterministic manner. See `alias_gensyms` below.
macro annotatetest(expr)
    return :(
        alias_gensyms(MacroTools.prettify(
            macroexpand(TestModule, $(QuoteNode(expr)); recursive = true); alias = false
        ))
    )
end

macro expected(expr)
    expanded = QuoteNode(macroexpand(TestModule, expr; recursive = true))
    return quote
        alias_gensyms(MacroTools.prettify($expanded; alias = false))
    end
end

# Deterministically replace gensyms with things that are easier to read.
# Good for testing, not for general usage.
const gensyms = [:animal, :butterfly, :caracal, :dingo]

function alias_gensyms(expr)
    map = Dict{Symbol,Symbol}()
    return MacroTools.postwalk(expr) do ex
        if MacroTools.isgensym(ex)
            return get!(map, ex, gensyms[length(map) + 1])
        end
        return ex
    end
end

@testset "Testing Library Utilities" begin
    @testset "Testing @wrapper" begin
        manager = CachedArrays.CacheManager(
            CachedArrays.AlignedAllocator(),
            CachedArrays.AlignedAllocator();
            localsize = 1_000_000_000,
            remotesize = 1_000_000_000,
            minallocation = 12,
        )

        x = CachedArray{Float32}(undef, manager, 10, 20)
        @test size(x) == (10, 20)
        rx = CachedArrays.writable(x)
        rx .= zero(Float32)

        y = transpose(x)
        @test size(y) == (20, 10)

        # Cannot index because the wrapped array is locked.
        @test_throws ErrorException y[1]
        @test_throws ErrorException y[1] = 2

        # Convert to a readable form.
        # `getindex` works but `setindex!` should not.
        readable_y = CachedArrays.readable(y)
        @test all(iszero, readable_y)
        @test_throws ErrorException readable_y[1] = 10

        # Convert to a writable form.
        # Changes to the writable form should be visible in the original array.
        writable_y = CachedArrays.writable(y)
        @test all(iszero, writable_y)
        writable_y .+= 1
        @test all(isequal(1), writable_y)
        @test all(isequal(1), rx)

        # Make sure we can release wrapper types as well.
        released_y = CachedArrays.release(y)
        @test_throws ErrorException released_y[1]
        @test_throws ErrorException released_y[1] = 10

        # Now - does this work on nested structs.
        z = reshape(transpose(x), 5, :)
        @test_throws ErrorException z[1]
        @test_throws ErrorException z[1] = 10

        readable_z = CachedArrays.readable(z)
        @test all(isequal(1), readable_z)
        @test_throws ErrorException readable_z[1] = 10

        writable_z = CachedArrays.writable(z)
        @test all(isequal(1), writable_z)
        writable_z .+= 1
        @test all(isequal(2), writable_z)
        @test all(isequal(2), rx)

        released_z = CachedArrays.release(z)
        @test_throws ErrorException released_z[1]
        @test_throws ErrorException released_z[1] = 10

        # Nest in the opposite order
        w = transpose(reshape(x, 5, :))
        @test_throws ErrorException w[1]
        @test_throws ErrorException w[1] = 10

        readable_w = CachedArrays.readable(w)
        @test all(isequal(2), readable_w)
        @test_throws ErrorException readable_w[1] = 10

        writable_w = CachedArrays.writable(w)
        @test all(isequal(2), writable_w)
        writable_w .+= 1
        @test all(isequal(3), writable_w)
        @test all(isequal(3), rx)

        released_w = CachedArrays.release(w)
        @test_throws ErrorException released_w[1]
        @test_throws ErrorException released_w[1] = 10
    end

    @testset "Testing @annotate" begin
        # Test "maybe_process_call" works on keywords and ignores nonkeywords.
        for keyword in CachedArrays.KEYWORDS
            sym = Symbol("__$(keyword)__")
            @test CachedArrays.maybe_process_call(sym) == :(CachedArrays.$(Symbol(keyword)))
        end

        @test CachedArrays.maybe_process_call(:__invoke__) == :__invoke__
        @test CachedArrays.maybe_process_call(:__recurse__) == :__recurse__
        @test CachedArrays.maybe_process_call(:__notakeyword__) == :__notakeyword__
        @test CachedArrays.maybe_process_call(:release) == :release

        #####
        ##### Test Cases
        #####

        # # Detect if MacroTools gets updated and add the `alias` keyword.
        # @test_throws MethodError MacroTools.prettify(:(myexpr); alias = false)

        ### Test 1
        fn = @annotatetest CachedArrays.@annotate function Mod.fn(
            x::MaybeTranspose{<:UnreadableCachedArray}, desc
        )
            return __recurse__(__readable__(x), desc)
        end
        ex = @expected function Mod.fn(x::MaybeTranspose{<:UnreadableCachedArray}, desc;)
            return Mod.fn(CachedArrays.readable(x), desc)
        end
        @test fn == ex

        ### Test 2
        fn = @annotatetest CachedArrays.@annotate function Mod.another_fn(
            f, data::CachedArray
        )
            return __invoke__(f, __writable__(data))
        end

        # Test replacement order
        @test gensyms[1] == :animal
        ex = @expected function Mod.another_fn(f, data::CachedArray;)
            return begin
                animal = (f, CachedArrays.writable(data))
                invoke(
                    Mod.another_fn,
                    CachedArrays.Tuple{CachedArrays.maybesuper.(animal)...},
                    animal...;
                )
            end
        end

        ### Test 3
        fn = @annotatetest CachedArrays.@annotate function (dot::SomeModule.Callable)(
            x::UnreadableCachedArray, ys::ReadableCachedArray; kw...
        )
            return dot(__readable__(x), ys; kw...)
        end
        ex = @expected function (dot::SomeModule.Callable)(
            x::UnreadableCachedArray, ys::ReadableCachedArray; kw...
        )
            return dot(CachedArrays.readable(x), ys; kw...)
        end
        @test fn == ex

        ### Test 4
        # Ensure that functions scoped outside the body of the macro are called correctly.
        fn = @annotatetest CachedArrays.@annotate function (dot::SomeModule.Callable)(
            x::UnreadableCachedArray
        )
            y = CachedArrays.@noescape __manager__(x) __recurse__(__readable__(x))
            CachedArrays.onobjects(dot) do o
                somefunction!(o)
            end
            return __release__(y)
        end

        ex = @expected function (dot::SomeModule.Callable)(x::UnreadableCachedArray;)
            y = CachedArrays.@noescape CachedArrays.manager(x) (dot::SomeModule.Callable)(
                CachedArrays.readable(x)
            )
            CachedArrays.onobjects(dot) do o
                somefunction!(o)
            end
            return CachedArrays.release(y)
        end
        @test fn == ex

        ### Add corner cases as they are encountered.
    end
end


text:
import CUDAnative, CUDAdrv
import CUDAnative: cufunction, DevicePtr
import CUDAdrv: CuEvent, CuStream, CuDefaultStream

const FREE_STREAMS = CuStream[]
const STREAMS = CuStream[]
const STREAM_GC_THRESHOLD = Ref{Int}(16)

# This code is loaded after an `@init` step
if haskey(ENV, "KERNELABSTRACTIONS_STREAMS_GC_THRESHOLD")
    global STREAM_GC_THRESHOLD[] = parse(Int, ENV["KERNELABSTRACTIONS_STREAMS_GC_THRESHOLD"])
end

## Stream GC
# Simplistic stream gc design in which when we have a total number
# of streams bigger than a threshold, we start scanning the streams
# and add them back to the freelist if all work on them has completed.
# Alternative designs:
# - Enqueue a host function on the stream that adds the stream back to the freelist
# - Attach a finalizer to events that adds the stream back to the freelist
# Possible improvements
# - Add a background task that occasionally scans all streams
# - Add a hysterisis by checking a "since last scanned" timestamp
# - Add locking
function next_stream()
    if !isempty(FREE_STREAMS)
        return pop!(FREE_STREAMS)
    end

    if length(STREAMS) > STREAM_GC_THRESHOLD[]
        for stream in STREAMS
            if CUDAdrv.query(stream)
                push!(FREE_STREAMS, stream)
            end
        end
    end

    if !isempty(FREE_STREAMS)
        return pop!(FREE_STREAMS)
    end

    stream = CUDAdrv.CuStream(CUDAdrv.STREAM_NON_BLOCKING)
    push!(STREAMS, stream)
    return stream
end

struct CudaEvent <: Event
    event::CuEvent
end

function Event(::CUDA)
    stream = CUDAdrv.CuDefaultStream()
    event = CuEvent(CUDAdrv.EVENT_DISABLE_TIMING)
    CUDAdrv.record(event, stream)
    CudaEvent(event)
end

wait(ev::CudaEvent, progress=nothing) = wait(CPU(), ev, progress)
function wait(::CPU, ev::CudaEvent, progress=nothing)
    if progress === nothing
        CUDAdrv.synchronize(ev.event)
    else
        while !CUDAdrv.query(ev.event)
            progress()
            # do we need to `yield` here?
        end
    end
end

# Use this to synchronize between computation using the CuDefaultStream
wait(::CUDA, ev::CudaEvent, progress=nothing) = __enqueue_wait(ev, CUDAdrv.CuDefaultStream())

# There is no efficient wait for CPU->GPU synchronization, so instead we
# do a CPU wait, and therefore block anyone from submitting more work.
# We maybe could do a spinning wait on the GPU and atomic flag to signal from the CPU,
# but which stream would we target?
wait(::CUDA, ev::CPUEvent,  progress=nothing) = wait(CPU(), ev, progress)

function __enqueue_wait(ev::CudaEvent, stream::CuStream)
    CUDAdrv.wait(ev.event, stream)
end

function (obj::Kernel{CUDA})(args...; ndrange=nothing, dependencies=nothing, workgroupsize=nothing)
    if ndrange isa Integer
        ndrange = (ndrange,)
    end
    if workgroupsize isa Integer
        workgroupsize = (workgroupsize, )
    end
    if dependencies isa Event
        dependencies = (dependencies,)
    end

    stream = next_stream()
    if dependencies !== nothing
        for event in dependencies
            if event isa CudaEvent
                __enqueue_wait(event, stream)
            end
        end
        for event in dependencies
            if !(event isa CudaEvent)
                wait(CUDA(), event, ()->yield())
            end
        end
    end

    if KernelAbstractions.workgroupsize(obj) <: DynamicSize && workgroupsize === nothing
        # TODO: allow for NDRange{1, DynamicSize, DynamicSize}(nothing, nothing)
        #       and actually use CUDAnative autotuning
        workgroupsize = (256,)
    end
    # If the kernel is statically sized we can tell the compiler about that
    if KernelAbstractions.workgroupsize(obj) <: StaticSize
        maxthreads = prod(get(KernelAbstractions.workgroupsize(obj)))
    else
        maxthreads = nothing
    end

    iterspace, dynamic = partition(obj, ndrange, workgroupsize)

    nblocks = length(blocks(iterspace))
    threads = length(workitems(iterspace))

    ctx = mkcontext(obj, ndrange, iterspace)
    # Launch kernel
    event = CuEvent(CUDAdrv.EVENT_DISABLE_TIMING)
    CUDAnative.@cuda(threads=threads, blocks=nblocks, stream=stream,
                     name=String(nameof(obj.f)), maxthreads=maxthreads,
                     Cassette.overdub(ctx, obj.f, args...))

    CUDAdrv.record(event, stream)
    return CudaEvent(event)
end

Cassette.@context CUDACtx

function mkcontext(kernel::Kernel{CUDA}, _ndrange, iterspace)
    metadata = CompilerMetadata{ndrange(kernel), true}(_ndrange, iterspace)
    Cassette.disablehooks(CUDACtx(pass = CompilerPass, metadata=metadata))
end

@inline function Cassette.overdub(ctx::CUDACtx, ::typeof(__index_Local_Linear))
    return CUDAnative.threadIdx().x
end

@inline function Cassette.overdub(ctx::CUDACtx, ::typeof(__index_Group_Linear))
    return CUDAnative.blockIdx().x
end

@inline function Cassette.overdub(ctx::CUDACtx, ::typeof(__index_Global_Linear))
    I =  @inbounds expand(__iterspace(ctx.metadata), CUDAnative.blockIdx().x, CUDAnative.threadIdx().x)
    # TODO: This is unfortunate, can we get the linear index cheaper
    @inbounds LinearIndices(__ndrange(ctx.metadata))[I]
end

@inline function Cassette.overdub(ctx::CUDACtx, ::typeof(__index_Local_Cartesian))
    @inbounds workitems(__iterspace(ctx.metadata))[CUDAnative.threadIdx().x]
end

@inline function Cassette.overdub(ctx::CUDACtx, ::typeof(__index_Group_Cartesian))
    @inbounds blocks(__iterspace(ctx.metadata))[CUDAnative.blockIdx().x]
end

@inline function Cassette.overdub(ctx::CUDACtx, ::typeof(__index_Global_Cartesian))
    return @inbounds expand(__iterspace(ctx.metadata), CUDAnative.blockIdx().x, CUDAnative.threadIdx().x)
end

@inline function Cassette.overdub(ctx::CUDACtx, ::typeof(__validindex))
    if __dynamic_checkbounds(ctx.metadata)
        I = @inbounds expand(__iterspace(ctx.metadata), CUDAnative.blockIdx().x, CUDAnative.threadIdx().x)
        return I in __ndrange(ctx.metadata)
    else
        return true
    end
end

generate_overdubs(CUDACtx)

###
# CUDA specific method rewrites
###

@inline Cassette.overdub(::CUDACtx, ::typeof(^), x::Float64, y::Float64) = CUDAnative.pow(x, y)
@inline Cassette.overdub(::CUDACtx, ::typeof(^), x::Float32, y::Float32) = CUDAnative.pow(x, y)
@inline Cassette.overdub(::CUDACtx, ::typeof(^), x::Float64, y::Int32)   = CUDAnative.pow(x, y)
@inline Cassette.overdub(::CUDACtx, ::typeof(^), x::Float32, y::Int32)   = CUDAnative.pow(x, y)
@inline Cassette.overdub(::CUDACtx, ::typeof(^), x::Union{Float32, Float64}, y::Int64) = CUDAnative.pow(x, y)

# libdevice.jl
const cudafuns = (:cos, :cospi, :sin, :sinpi, :tan,
          :acos, :asin, :atan,
          :cosh, :sinh, :tanh,
          :acosh, :asinh, :atanh,
          :log, :log10, :log1p, :log2,
          :exp, :exp2, :exp10, :expm1, :ldexp,
          # :isfinite, :isinf, :isnan, :signbit,
          :abs,
          :sqrt, :cbrt,
          :ceil, :floor,)
for f in cudafuns
    @eval function Cassette.overdub(ctx::CUDACtx, ::typeof(Base.$f), x::Union{Float32, Float64})
        @Base._inline_meta
        return CUDAnative.$f(x)
    end
end

###
# GPU implementation of shared memory
###
@inline function Cassette.overdub(ctx::CUDACtx, ::typeof(SharedMemory), ::Type{T}, ::Val{Dims}, ::Val{Id}) where {T, Dims, Id}
    ptr = CUDAnative._shmem(Val(Id), T, Val(prod(Dims)))
    CUDAnative.CuDeviceArray(Dims, ptr)
end

###
# GPU implementation of scratch memory
# - private memory for each workitem
###

@inline function Cassette.overdub(ctx::CUDACtx, ::typeof(Scratchpad), ::Type{T}, ::Val{Dims}) where {T, Dims}
    MArray{__size(Dims), T}(undef)
end

@inline function Cassette.overdub(ctx::CUDACtx, ::typeof(__synchronize))
    CUDAnative.sync_threads()
end

###
# GPU implementation of `@Const`
###
struct ConstCuDeviceArray{T,N,A} <: AbstractArray{T,N}
    shape::Dims{N}
    ptr::DevicePtr{T,A}

    # inner constructors, fully parameterized, exact types (ie. Int not <:Integer)
    ConstCuDeviceArray{T,N,A}(shape::Dims{N}, ptr::DevicePtr{T,A}) where {T,A,N} = new(shape,ptr)
end

Adapt.adapt_storage(to::ConstAdaptor, a::CUDAnative.CuDeviceArray{T,N,A}) where {T,N,A} = ConstCuDeviceArray{T, N, A}(a.shape, a.ptr)

Base.pointer(a::ConstCuDeviceArray) = a.ptr
Base.pointer(a::ConstCuDeviceArray, i::Integer) =
    pointer(a) + (i - 1) * Base.elsize(a)

Base.elsize(::Type{<:ConstCuDeviceArray{T}}) where {T} = sizeof(T)
Base.size(g::ConstCuDeviceArray) = g.shape
Base.length(g::ConstCuDeviceArray) = prod(g.shape)
Base.IndexStyle(::Type{<:ConstCuDeviceArray}) = Base.IndexLinear()

Base.unsafe_convert(::Type{DevicePtr{T,A}}, a::ConstCuDeviceArray{T,N,A}) where {T,A,N} = pointer(a)

@inline function Base.getindex(A::ConstCuDeviceArray{T}, index::Integer) where {T}
    @boundscheck checkbounds(A, index)
    align = Base.datatype_alignment(T)
    CUDAnative.unsafe_cached_load(pointer(A), index, Val(align))::T
end

@inline function Base.unsafe_view(arr::ConstCuDeviceArray{T, 1, A}, I::Vararg{Base.ViewIndex,1}) where {T, A}
    ptr = pointer(arr) + (I[1].start-1)*sizeof(T)
    len = I[1].stop - I[1].start + 1
    return ConstCuDeviceArray{T,1,A}(len, ptr)
end


text:
function heki(y::Union{Array{R,1}, Array{Array{R,1},1}},
             σ::R,
             η::R,
             J::Integer,
             N::Integer,
             priorτ::Function,
             priorθ::Function,
             gmap::Function,
             tmap::Function;
             ρ::R = convert(R,0.5),
             ζ::R = convert(R,2.0),
             γ0::R = convert(R,1e0),
             batched::Bool = false,
             batches::Int = 1,
             parallel::Bool = false,
             verbosity::Int=0,
             rerandomize::Bool = false,
             rerandom_coeff::R = convert(R,0.25)) where R<:Real
    if batched
        @assert typeof(y) == Array{Array{R,1}, 1} "Must supply a list of data subsets for batched HEKI, current type is $(typeof(y))"
        heki_batched(y, σ, η, J, N, priorτ, priorθ, gmap, tmap; ρ = ρ, ζ = ζ, γ0 = γ0, parallel = parallel, verbosity=verbosity, rerandomize=rerandomize, rerandom_coeff=rerandom_coeff, batches=batches, batch_off=batch_off)
    else
        @assert typeof(y) == Array{R,1} "Must supply array of flat data for unbatched HEKI, current type is $(typeof(y))"
        heki_nobatch(y, σ, η, J, N, priorτ, priorθ, gmap, tmap; ρ = ρ, ζ = ζ, γ0 = γ0, parallel = parallel, verbosity=verbosity, rerandomize=rerandomize, rerandom_coeff=rerandom_coeff)
    end
end

function heki_ensemble_update!(τ::Union{Array{R,1}, Array{Array{R,1},1}},
                               θ::Union{Array{R,1}, Array{Array{R,1},1}},
                               w::Union{Array{R,1}, Array{Array{R,1},1}},
                               wm::Array{R,1},
                               y::Array{R,1},
                               γ0::R,
                               σ::R,
                               ρ::R,
                               T::UniformScaling{R},
                               verbosity::Int=0) where R<:Real
        #Analysis
        Cτw = CrossCovarianceOperator(hcat(τ...)', hcat(w...)')
        Cθw = CrossCovarianceOperator(hcat(θ...)', hcat(w...)')
        Cwwf = pheigfact(CovarianceOperator(hcat(w...)')) #this is a low-rank approx to Cww
        #set regularization
        γ = setγ(γ0, ρ, T, Cwwf, y, wm)
        if verbosity >= 2
            println("γ = $γ")
        end
        for j = 1:size(τ)[1]
            yj = y.+σ.*randn(R, length(y))
            τ[j] = τ[j].+Cτw*(wbinv(T/γ,
            Cwwf[:vectors],
            Diagonal(convert(R,1.0)./Cwwf[:values]),
            (yj-w[j])))

            θ[j] = θ[j].+Cθw*(wbinv(T/γ,
            Cwwf[:vectors],
            Diagonal(convert(R,1.0)./Cwwf[:values]),
            (yj-w[j])))
        end
end

function heki_nobatch(y::Array{R,1},
                      σ::R,
                      η::R,
                      J::Integer,
                      N::Integer,
                      priorτ::Function,
                      priorθ::Function,
                      gmap::Function,
                      tmap::Function;
                      ρ::R = convert(R,0.5),
                      ζ::R = convert(R,2.0),
                      γ0::R = convert(R,1.0),
                      parallel::Bool = false,
                      verbosity::Int=0,
                      rerandomize::Bool = false,
                      rerandom_coeff::R = convert(R,0.25)) where {R<:Real}
    #Initialization
    Γ = σ^2 * I #this should use the efficient uniform scaling operator - can use both left & right division also
    T = σ^(-2) * I # Precision Matrix
    τ = [priorτ() for j = 1:J]
    θ = [priorθ() for j = 1:J]
    tmap1arg(x) = tmap(x[1], x[2])
    #Main Optimization loop
    if verbosity > 0
        println("Starting up to $N iterations with $J ensemble members")
    end
    for i = 1:N
        if verbosity > 0
            println("Iteration # $i. Starting Forward Map")
        end
        #Prediction
        #project hierarchical parameters to parameters of interest
        if parallel
            u = pmap(tmap1arg, zip(τ, θ))
        else
            u = map(tmap1arg, zip(τ, θ))
        end
        #predict data
        if parallel
            w = pmap(gmap, u)
        else
            w = map(gmap, u)
        end
        wm = mean(w)
        #Discrepancy principle
        convg = sqrt((y-wm)'*T*(y-wm))
        if verbosity > 0
            println("Iteration # $i. Discrepancy Check; Weighted Norm: $convg, Noise level: $(ζ*η)")
        end
        if convg <= ζ*η
            return (mean(τ), mean(θ))
        end
        heki_ensemble_update!(τ,θ,w,wm,y,γ0,σ,ρ,T,verbosity)
    end
    (mean(τ), mean(θ))
end

function heki_batched(y::Array{Array{R,1},1},
                    σ::R,
                    η::R,
                    J::Integer,
                    N::Integer,
                    priorτ::Function,
                    priorθ::Function,
                    gmap::Function,
                    tmap::Function;
                    ρ::R = convert(R,0.5),
                    ζ::R = convert(R,2.0),
                    γ0::R = convert(R,1.0),
                    parallel::Bool = false,
                    verbosity::Int=0,
                    rerandomize::Bool = false,
                    rerandom_coeff::R = convert(R,0.25),
                    batches::Int=1,
                    batch_off::Int=20) where {R<:Real}
    #Initialization
    Γ = σ^2 * I #this should use the efficient uniform scaling operator - can use both left & right division also
    T = σ^(-2) * I # Precision Matrix
    τ = [priorτ() for j = 1:J]
    θ = [priorθ() for j = 1:J]
    tmap1arg(x) = tmap(x[1], x[2])
    ya = vcat(y...)
    #Main Optimization loop
    if verbosity > 0
        println("Starting up to $N iterations with $J ensemble members")
    end
    for i = 1:N
        if verbosity > 0
            println("Iteration # $i. Starting Forward Map")
        end
        #project hierarchical parameters to parameters of interest
        if parallel
            u = pmap(tmap1arg, zip(τ, θ))
        else
            u = map(tmap1arg, zip(τ, θ))
        end
        #Prediction
        if parallel
            w = pmap(x->vcat([gmap(x,b) for b in 1:length(y)]...), u)
        else
            w = map(x->vcat([gmap(x,b) for b in 1:length(y)]...), u)
        end
        #Discrepancy principle
        wm = mean(w)
        convg = sqrt((ya-wm)'*T*(ya-wm))
        if verbosity > 0
            println("Iteration # $i. Discrepancy Check; Weighted Norm: $convg, Noise level: $(ζ*η)")
        end
        if convg <= ζ*η
            return (mean(τ), mean(θ))
        end
        #set up batch ordering
        if i < batch_off
            parts = kfoldperm(length(y), batches) #randomize batch groupings each iterations
        else
            parts = kfoldperm(length(y), 1) #turn off minibatching after you have approached the minimum
        end
        for (i, p) in enumerate(parts)
            if parallel
                w = pmap(x->vcat([gmap(x,b) for b in p]...), u)
            else
                w = map(x->vcat([gmap(x,b) for b in p]...), u)
            end
            wm = mean(w)
            ensemble_update!(τ,θ,w,wm,vcat(y[p]...),γ0,σ,ρ,T,verbosity)
        end
    end
    (mean(τ), mean(θ))
end


text:
using Test
using GEMPIC

include("test_maxwell_2d_fem.jl")
include("test_vm_1d1v_spin.jl")
include("test_vm_1d2v.jl")
include("test_mesh.jl")
include("test_spline_pp.jl")
include("test_maxwell_1d_fem.jl")
include("test_sampling.jl")
include("test_sampling_spin.jl")
include("test_particle_mesh_coupling_spline_1d.jl")
include("test_particle_mesh_coupling_spline_2d.jl")

if Threads.nthreads() <= 2
    include("test_hamiltonian_splitting.jl")
end

include("test_hamiltonian_splitting_spin.jl")
include("test_hamiltonian_splitting_boris.jl")
include("test_save_particles.jl")


text:
function exact_expectation_inside_log_of_lower_bound_grad(; y=y, idx=idx, S=S, μα=μα, Σα=Σα, v²ₐ=v²ₐ, μβ=μβ, Σβ=Σβ, v²ᵦ=v²ᵦ, μf=μf, Σf=Σf, K=K)

    gradα, gradΣα, gradβ, gradΣβ = expectationterm1_grad(;y=y, μα=μα, Σα=Σα, μβ=μβ, Σβ=Σβ, μf=μf, Σf=Σf, S=S, idx=idx)

    auxα, auxΣα  = expectationterm3_grad(; μα=μα, Σα=Σα, v²ₐ=v²ₐ)

    gradα  += auxα

    gradΣα += auxΣα

    auxβ, auxΣβ = expectationterm4_grad(; μβ=μβ, Σβ=Σβ, v²ᵦ=v²ᵦ)

    gradβ  += auxβ

    gradΣβ += auxΣβ

    return gradα, gradΣα, gradβ, gradΣβ

end

##########
# TERM 1 #
##########

function expectationterm1_grad(;y=y, μα=μα, Σα=Σα, μβ=μβ, Σβ=Σβ, μf=μf, Σf=Σf, S=S, idx=idx)

    N = length(y)

    ## look out for the diagonal and the diff between Σα and Σa
    Ā = Diagonal(μα[idx]) ; Σa = Diagonal(diag(Σα)[idx])

    ## look out for the diagonal and the diff between Σβ and Σb
    b̄ = μβ[idx] ; Σb = Diagonal(diag(Σβ)[idx])

    # subterm0 = -0.5*logdet(S) - 0.5*N*log(2π)
    #
    # subterm1 = -0.5*dot(y, S\y)
    #
    # subterm2 = -0.5*dot(y, -S\Ā*μf) * 2
    #
    # subterm3 = -0.5*dot(y, -S\b̄) * 2
    #
    # # subterm4 = -0.5*sum((μf.^2 .+ diag(Σf)) .* (Ā.diag.^2 .+ Σa.diag) ./ S.diag)
    # subterm4 = -0.5*(μf'*(S\(Ā.^2 + Σa))*μf)  -0.5*tr(Σf*(S\(Ā.^2 + Σa)))
    #
    # subterm5 = -0.5*dot(b̄, S\Ā*μf) * 2
    #
    # subterm6 = -0.5*dot(b̄, S\b̄) -0.5*tr(S\Σb)
    #
    # (subterm0 + subterm1 + subterm2 + subterm3 + subterm4 + subterm5 + subterm6)

    #########
    # gradα #
    #########

    gradα = zeros(length(μα))

    for n = 1:N

        # from subterm (2)
        gradα[idx[n]] +=   y[n] * μf[n] / S[n,n]

        # from subterm (4)
        gradα[idx[n]] += - (μf[n]^2 +  Σf[n,n]) / S[n,n] * Ā[n,n]

        # from subterm (5)
        gradα[idx[n]] += - b̄[n] * μf[n] / S[n,n]

    end

    ##########
    # gradΣα #
    ##########

    gradΣα = zeros(length(μα)); @assert(length(μα) == size(Σα,1))

    for n = 1:N

        # subterm (4)
        # remember we parametrise the n-th entry of the diagonal covariance
        # Σa with elements exp(pₙ), i.e. gradient calculated below is dependent
        # on this parametrisation (see "unpack" function)
        gradΣα[idx[n]] += - 0.5*(μf[n]^2 +  Σf[n,n]) / S[n,n] *  Σa[n,n]

    end

    #########
    # gradβ #
    #########

    gradβ = zeros(length(μβ))

    for n = 1:N

        # from subterm (3)
        gradβ[idx[n]] +=   y[n] / S[n,n]

        # from subterm (5)
        gradβ[idx[n]] += - Ā[n,n] * μf[n] / S[n,n]

        # from subterm (6)
        gradβ[idx[n]] += - b̄[n] / S[n,n]

    end

    ##########
    # gradΣβ #
    ##########

    gradΣβ = zeros(length(μβ)); @assert(length(μβ) == size(Σβ,1))

    for n = 1:N

        # from subterm 6
        gradΣβ[idx[n]] += - 0.5 * Σb[n, n] / S[n,n]

    end

    return gradα, gradΣα, gradβ, gradΣβ

end

##########
# TERM 2 #
##########

# function expectationterm2_grad(;μf=μf, Σf=Σf, K=K)
#     N = length(μf)
#     -0.5*N*log(2π) -0.5*logdet(K) - 0.5*μf'*(K\μf) - 0.5*tr(K\Σf)
# end

##########
# TERM 3 #
##########

function expectationterm3_grad(;μα=μα, Σα=Σα, v²ₐ=v²ₐ)
    # F = length(μα)
    # - 0.5*F*log(2π) - 0.5*F*log(v²ₐ) - 0.5*dot(μα, μα)/v²ₐ - 0.5*tr(Σα)/v²ₐ

    gradα = - μα / v²ₐ

    gradΣα = -0.5*diag(Σα)/v²ₐ

    return gradα, gradΣα

end

##########
# TERM 4 #
##########

function expectationterm4_grad(;μβ=μβ, Σβ=Σβ, v²ᵦ=v²ᵦ)
    # F = length(μβ)
    # - 0.5*F*log(2π) - 0.5*F*log(v²ᵦ) - 0.5*dot(μβ,μβ)/v²ᵦ - 0.5*tr(Σβ)/v²ᵦ

    gradβ = - μβ / v²ᵦ

    gradΣβ = -0.5*diag(Σβ)/v²ᵦ

    return gradβ, gradΣβ

end


text:
# Direct Search with Probabilistic Descent as described in Gratton2014:
#
#  S. Gratton, C. W. Royer, L. N. Vicente, and Z. Zhang, Direct search
#  based on probabilistic descent, preprint 14-11, Dept. Mathematics, Univ. Coimbra.
#  http://www.mat.uc.pt/~lnv/papers/ds-random.pdf
#

"""
Generate `num` random vectors on the `n`-dimensional, unit (hyper)sphere.
This is the Muller-Marsaglia method as described [here](http://mathworld.wolfram.com/HyperspherePointPicking.html).
"""
function sample_unit_hypersphere(n, num = 1)
    X = randn(n, num)
    sqrootsums = 1 ./ sqrt.(sum(abs2, X, dims=1))
    sqrootsums .* X
end

struct RandomDirectionGen <: DirectionGenerator
    numDimensions::Int
    numDirections::Int
end

directions_for_k(rdg::RandomDirectionGen, k) =
    sample_unit_hypersphere(rdg.numDimensions, rdg.numDirections)

"""
Generate half of the directions randomly and then
mirrors by negating them.
"""
struct MirroredRandomDirectionGen <: DirectionGenerator
    numDimensions::Int
    numDirections::Int

    function MirroredRandomDirectionGen(numDims, numDirections)
        iseven(numDirections) ||
            throw(ArgumentError("the number of directions must be even"))
        new(numDims, numDirections)
    end
end

function directions_for_k(rdg::MirroredRandomDirectionGen, k)
    r = sample_unit_hypersphere(rdg.numDimensions, rdg.numDirections÷2)
    [r -r]
end

const DirectSearchProbabilisticDescentDefaultParameters = ParamsDict(
    :NumDirections => 2, # This should be a function of Gamma and Phi for the GSS but 2 is often enough
)

function direct_search_probabilistic_descent(problem::OptimizationProblem, parameters::Parameters)
    params = chain(DirectSearchProbabilisticDescentDefaultParameters, parameters)
    params[:DirectionGenerator] = MirroredRandomDirectionGen(numdims(problem), params[:NumDirections])
    GeneratingSetSearcher(problem, params)
end