text:
@testset "Rename" begin
    @testset "Rename named tuples" begin
        @test rename((a = 1,), :a, :b) === (b = 1,)
        @test rename((a = 1,), :a, :a) === (a = 1,)
        @test rename((a = 1, c = 2.0), :a, :b) === (b = 1, c = 2.0)
        @test rename((a = 1, c = 2.0), (:a, :c), (:aa, :cc)) === (aa = 1, cc = 2.0)
        @test rename((a = 1, c = 2.0), (:a, :c), (:c, :a)) === (c = 1, a = 2.0)

        # Not sure if this one is correct or should throw an error...
        @test rename((a = 1,), :b, :a) === (a = 1,)

        @test_throws ErrorException rename((a = 1, c = 2.0), :a, :c)
    end

    @testset "Rename tables" begin
        va = [2,   3,   1]
        vb = [4.0, 6.0, 2.0]
        t = Table{(:a, :b), NamedTuple{(:a,:b), Tuple{Int, Float64}}, Tuple{Vector{Int}, Vector{Float64}}, Tuple{}}((va, vb), ())
        r = Rename(:a, :c)

        @test r(accelerate(t, UniqueIndex{(:a,)})) == Table(c = va, b = vb)
        @test r(accelerate(t, UniqueIndex{(:b,)})) == Table(c = va, b = vb)
        @test r(accelerate(t, UniqueIndex{(:a,:b)})) == Table(c = va, b = vb)
        @test r(accelerate(t, UniqueIndex{(:b,:a)})) == Table(c = va, b = vb)

        @test r(accelerate(t, SortIndex{(:a,)})) == Table(c = va, b = vb)
        @test r(accelerate(t, SortIndex{(:b,)})) == Table(c = va, b = vb)
        @test r(accelerate(t, SortIndex{(:a,:b)})) == Table(c = va, b = vb)
        @test r(accelerate(t, SortIndex{(:b,:a)})) == Table(c = va, b = vb)

        @test r(accelerate(t, UniqueSortIndex{(:a,)})) == Table(c = va, b = vb)
        @test r(accelerate(t, UniqueSortIndex{(:b,)})) == Table(c = va, b = vb)
        @test r(accelerate(t, UniqueSortIndex{(:a,:b)})) == Table(c = va, b = vb)
        @test r(accelerate(t, UniqueSortIndex{(:b,:a)})) == Table(c = va, b = vb)

        @test r(accelerate(t, HashIndex{(:a,)})) == Table(c = va, b = vb)
        @test r(accelerate(t, HashIndex{(:b,)})) == Table(c = va, b = vb)
        @test r(accelerate(t, HashIndex{(:a,:b)})) == Table(c = va, b = vb)
        @test r(accelerate(t, HashIndex{(:b,:a)})) == Table(c = va, b = vb)

        @test r(accelerate(t, UniqueHashIndex{(:a,)})) == Table(c = va, b = vb)
        @test r(accelerate(t, UniqueHashIndex{(:b,)})) == Table(c = va, b = vb)
        @test r(accelerate(t, UniqueHashIndex{(:a,:b)})) == Table(c = va, b = vb)
        @test r(accelerate(t, UniqueHashIndex{(:b,:a)})) == Table(c = va, b = vb)
    end
end


text:
abstract type Sim end

struct POMDPSim <: Sim
    simulator::Simulator
    pomdp::POMDP
    policy::Policy
    updater::Updater
    initial_belief::Any
    initial_state::Any
    metadata::Dict{Symbol}
end

struct MDPSim <: Sim
    simulator::Simulator
    mdp::MDP
    policy::Policy
    initial_state::Any
    metadata::Dict{Symbol}
end

"""
    Sim(p::POMDP, policy::Policy, metadata=Dict(:note=>"a note"))
    Sim(p::POMDP, policy::Policy[, updater[, initial_belief[, initial_state]]]; kwargs...)

Create a `Sim` object that represents a POMDP simulation.
"""
function Sim(pomdp::POMDP,
                    policy::Policy,
                    up=updater(policy),
                    initial_belief=initial_state_distribution(pomdp),
                    initial_state=nothing;
                    rng::AbstractRNG=Base.GLOBAL_RNG,
                    max_steps::Int=typemax(Int),
                    simulator::Simulator=HistoryRecorder(rng=rng, max_steps=max_steps),
                    metadata::Dict{Symbol}=Dict{Symbol, Any}()
                   )

    if initial_state == nothing && state_type(pomdp) != Void
        is = rand(rng, initial_belief)
    else
        is = initial_state
    end
    return POMDPSim(simulator, pomdp, policy, up, initial_belief, is, metadata)
end

"""
    Sim(p::MDP, policy::Policy, metadata=Dict(:note=>"a note"))
    Sim(p::MDP, policy::Policy[, initial_state]; kwargs...)

Create a `Sim` object that represents a MDP simulation.

A vector of `Sim` objects can be executed with `run` or `run_parallel`.

## Keyword Arguments
- `rng::AbstractRNG=Base.GLOBAL_RNG`
- `max_steps::Int=typemax(Int)`
- `simulator::Simulator=HistoryRecorder(rng=rng, max_steps=max_steps)`
- `metadata::Dict{Symbol}=Dict{Symbol, Any}()` a dictionary of metadata for the sim that will be recorded, e.g. `Dict(:solver_iterations=>500)`.
"""
function Sim(mdp::MDP,
             policy::Policy,
             initial_state=nothing;
             rng::AbstractRNG=Base.GLOBAL_RNG,
             max_steps::Int=typemax(Int),
             simulator::Simulator=HistoryRecorder(rng=rng, max_steps=max_steps),
             metadata::Dict{Symbol}=Dict{Symbol, Any}()
            )

    if initial_state == nothing && state_type(pomdp) != Void
        is = initial_state(mdp, rng)
    else
        is = initial_state
    end
    return MDPSim(simulator, mdp, policy, is, metadata)
end

POMDPs.simulate(s::POMDPSim) = simulate(s.simulator, s.pomdp, s.policy, s.updater, s.initial_belief, s.initial_state)
POMDPs.simulate(s::MDPSim) = simulate(s.simulator, s.mdp, s.policy, s.initial_state)

function default_process(s::Sim, r::Float64)
    stuff = metadata_as_pairs(s)
    return push!(stuff, :reward=>r)
end
default_process(s::Sim, hist::SimHistory) = default_process(s, discounted_reward(hist))

metadata_as_pairs(s::Sim) = convert(Array{Any}, collect(s.metadata))

run_parallel(queue::AbstractVector) = run_parallel(default_process, queue)

function run_parallel(process::Function, queue::AbstractVector;
                      progress=Progress(length(queue), desc="Simulating..."))

    #=
    frame_lines = pmap(progress, queue) do sim
        result = simulate(sim)
        return process(sim, result)
    end
    =#

    # based on the simple implementation of pmap here: https://docs.julialang.org/en/latest/manual/parallel-computing
    np = nprocs()
    n = length(queue)
    i = 1
    prog = 0
    frame_lines = Vector{Any}(n)
    nextidx() = (idx=i; i+=1; idx)
    prog_lock = ReentrantLock()
    @sync begin
        for p in 1:np
            if np == 1 || p != myid()
                @async begin
                    while true
                        idx = nextidx()
                        if idx > n
                            break
                        end
                        frame_lines[idx] = remotecall_fetch(p, queue[idx]) do sim
                            result = simulate(sim)
                            return process(sim, result)
                        end
                        if progress isa Progress
                            lock(prog_lock)
                            update!(progress, prog+=1)
                            unlock(prog_lock)
                        end
                    end
                end
            end
        end
    end

    return create_dataframe(frame_lines)
end

Base.run(queue::AbstractVector) = run(default_process, queue)

function Base.run(process::Function, queue::AbstractVector; show_progress=true)
    lines = []
    if show_progress
        @showprogress for sim in queue
            result = simulate(sim)
            push!(lines, process(sim, result))
        end
    else
        for sim in queue
            result = simulate(sim)
            push!(lines, process(sim, result))
        end
    end
    return create_dataframe(lines)
end

function create_dataframe(lines::Vector)
    master = Dict{Symbol, DataArray}()
    for line in lines
        push_line!(master, line)
    end
    return DataFrame(master)
end

function _push_line!(d::Dict{Symbol, DataArray}, line)
    if isempty(d)
        len = 0
    else
        len = length(first(values(d)))
    end
    for (key, val) in line
        if !haskey(d, key)
            d[key] = DataArray(typeof(val), len)
        end
        data = d[key]
        if !(typeof(val) <: eltype(data))
            d[key] = convert(DataArray{Any,1}, data)
        end
        push!(d[key], val)
    end
    for da in values(d)
        if length(da) < len + 1
            push!(da, NA)
        end
    end
    return d
end
push_line!(d::Dict{Symbol, DataArray}, line::Dict) = _push_line!(d, line)
push_line!(d::Dict{Symbol, DataArray}, line::DataFrame) = _push_line!(d, n=>first(line[n]) for n in names(line))
push_line!(d::Dict{Symbol, DataArray}, line::Vector{Any}) = _push_line!(d, line)
push_line!(d::Dict{Symbol, DataArray}, line::Tuple) = _push_line!(d, line)

ext: jl


text:
# Heisenberg-1d.jl
# The MIT License (MIT)
# Copyright (c) 2020 Sota Yoshida

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.

using LinearAlgebra
using ArgParse
using SIMD
#using Arpack
#using SparseArrays
#using KrylovKit
#using BenchmarkTools
#using Profile

function prodS(vals::Array{Int64,1},numS::Int64,N::Int64)
    tlist = copy(vals)
    for i = 1:N-1
        nlist = []
        for tmp in tlist
            for val in vals
                if typeof(tmp) == Int
                    ctmp = [tmp]
                else
                    ctmp = copy(tmp)
                end
                push!(nlist,append!(ctmp,[val]))
            end
        end
        tlist = nlist
    end
    Dim = numS^N
    if length(tlist) != Dim
        println("These must be the same: length(tlist)=",length(tlist), " NumS^N=", Dim)
    end
    return tlist
end

function makeV(S::Float64,N::Int64)
    vals=collect(-Int64(2*S):2:Int64(2*S))
    numS = length(vals)
    prodS(vals,numS,N)
end

function convcheck(k::Int64,numstate::Int64,tDim::Int64,alphas::FA,betas::FA,en::FA,tol=1.e-14) where{FA<:Array{Float64,1}}
    if k < numstate
        return en,false
    else
        T = zeros(Float64,k,k)
        for i = 1:k
            T[i,i] =alphas[i]
        end
        for j = 2:k
            i = j-1
            T[i,j] = betas[i];T[j,i] = betas[i]
        end
        iend = minimum([tDim,numstate])

        nen = eigvals(T)[1:iend]
        if k == numstate; return nen,false; end
        dif = maximum( abs.(en[1:iend] .- nen[1:iend]))
        if dif < tol
            return nen,true
        else
            return nen,false
        end
    end
end

function operateH(v::Array{Float64,1},t::Float64,U::Float64,
                  non0s::Array{Array{Int64,1},1},facs::Array{Array{Float64,1},1}
                  ,Ufacs::Array{Float64,1},tDim::Int64)
    ### diagonal part
    w = Ufacs .* v
    ### non-diagonal part
    tvec = 0.0 * v
    for ith = 1:tDim-1
        for (ind, jth) in enumerate(non0s[ith])
            tvec[jth] += 0.5* t * facs[ith][ind] * v[ith]
            tvec[ith] += 0.5* t * facs[ith][ind] * v[jth]
        end
    end
    w += tvec
    return w
end

function Lanczos(t::Float64,U::Float64,non0s::AAI,facs::AAF,Ufacs::FA,tDim::Int64,
                 numstate::Int64,maxiter=1000) where{FA<:Array{Float64,1},
                                                     AAI<:Array{Array{Int64,1},1},
                                                     AAF<:Array{Array{Float64,1},1}}
    vkm1 = zeros(Float64,tDim)
    vk   = zeros(Float64,tDim); vk[1]=1.0
    betakm1 = 0.0; betak=0.0; alphas = Float64[]; betas = Float64[]
    en = [1.0]; TF = false

    for k = 1:maxiter
        w = operateH(vk,t,U,non0s,facs,Ufacs,tDim)  - betak .* vkm1
        alphak = dot(vk,w)
        push!(alphas,alphak)
        en,TF = convcheck(k,numstate,tDim,alphas,betas,en)
        if TF;println(" En=",en[1:minimum([numstate,tDim])]," Lanczos it=$k"); break; end
        if k==maxiter
            try;println("En=",en[1:minimum([tDim,numstate])]);catch;println("En=$en");end
            #println("vec=$vk")
        end
        w = w - alphak .* vk
        betakm1 = betak; betak = sqrt(dot(w,w)); push!(betas,betak)
        if betak == 0.0
            break
        end
        vkm1 = vk; vk = w/betak
    end
end

function calcfac(S::F,Sz_cre::F,Sz_ani::F) where{F<:Float64}
    sqrt( (S*(S+1.0) - Sz_cre*(Sz_cre+1.0)) * (S*(S+1.0) - Sz_ani*(Sz_ani-1.0)))
end

function neighbor(ret::IA,dif::IA,N::Int64) where{IA<:Array{Int64,1}}
    ij = Int64[ ]
    for i=1:length(dif)
        if dif[i] != 0
            push!(ij,i)
        end
    end
    if abs(ij[1]-ij[2])== 1
        if dif[ij[1]] == 1
            ret[1]= 1; ret[2]=ij[1]; ret[3]=ij[2]
        else
            ret[1]= 1; ret[3]=ij[1]; ret[2]=ij[2]
        end
    else
        ret[1]=0 ; ret[2]=0; ret[3]=0
    end
    nothing
end

function neighborPBC(ret::IA,dif::IA,N::Int64) where{IA<:Array{Int64,1}}
    ij = Int64[ ]
    for i=1:length(dif)
        if dif[i] != 0
            push!(ij,i)
        end
    end
    if abs(ij[1]-ij[2])== 1 || abs(ij[1]-ij[2])==N-1
        if dif[ij[1]] == 1
            ret[1]= 1; ret[2]=ij[1]; ret[3]=ij[2]
        else
            ret[1]= 1; ret[3]=ij[1]; ret[2]=ij[2]
        end
    else
        ret[1]=0 ; ret[2]=0; ret[3]=0
    end
    nothing
end

function connectable(tf::IT,dif::IT) where{IT<:Array{Int64,1}}
    tf[1] = ifelse( sum( abs.(dif)) != 4 || sum(dif) != 0, 0,1)
    nothing
end

function main(t::F,U::F,S::F,N::I,numstate::I,PBC::I) where{F<:Float64,I<:Int64}
    basisvec = makeV(S,N); println("Full Dim.=",length(basisvec))

    Slist = [ i for i=-S*N:1.0:S*N]
    vecsSz = [ [] for i=1:length(Slist)]
    for tvec in basisvec
        push!(vecsSz[Int64(0.5 * sum(tvec)+S*N)+1],tvec)
    end

    for i = 1:length(Slist)
        Sz = Slist[i]
        #if Sz != 0.0; continue ;end
        tvec = vecsSz[i]
        tDim = length(tvec)

        non0s=[Int64[] for ith=1:tDim]
        facs=[Float64[] for ith=1:tDim]
        #println("typeof(non0s) ",typeof(non0s), " typeof(facs) ", typeof(facs))
        Ufacs=zeros(Float64,tDim)
        tfac = 0.0; ret=zeros(Int64,3)
        print("Sz = $Sz tDim=$tDim\t")
        Fneighbor = ifelse(PBC==1, neighborPBC, neighbor)
        r = 0.0
        @inbounds for ith=1:tDim
            for j = 1:N-1
                Ufacs[ith] += 0.25*(tvec[ith][j]*tvec[ith][j+1])
            end
            Ufacs[ith] += ifelse(PBC==1, 0.25*tvec[ith][N]*tvec[ith][1],0)
        end
        @inbounds for jth=2:tDim
            @simd for ith = 1:jth
                connectable(ret,tvec[ith]-tvec[jth])
                if ret[1] == 1
                    Fneighbor(ret,tvec[ith]-tvec[jth],N)
                    if ret[1]==1
                        push!(non0s[ith],jth)
                        push!(facs[ith],calcfac(S,0.5*tvec[jth][ret[2]],0.5*tvec[jth][ret[3]]))
                        ###  #Sz_cre = 0.5 * ket[cre]; Sz_ani = 0.5 * ket[ani]
                    end ## cre = ret[2]; ani=ret[3]
                end
            end
        end
        if tDim == 1; println("En=$Ufacs");end
        Lanczos(t,U,non0s,facs,Ufacs,tDim,numstate)
    end
    nothing
end

function run()
    for (arg,val) in parsed_args
        print("$arg  => $val \t")
    end
    print("\n")
    t = parsed_args["t"]; S = parsed_args["S"];N = parsed_args["N"]
    nums = parsed_args["nums"]; PBC = parsed_args["PBC"]
    U = t  # (for debug) This must be identical in the 1D Heisenberg Hamiltonian"
    @time main(t,U,S,N,nums,PBC)
end

s = ArgParseSettings()
@add_arg_table! s begin
    "--t"
        help = "coefficient in the 1D Heisenberg Hamiltonian"
        arg_type = Float64
        default = 1.0
    "--S"
        help = "Spin (in FLoat) 0.5(S=1/2), 1.0 (S=1), etc."
        arg_type = Float64
        default = 0.5
    "--N"
        help = "particle number N"
        arg_type = Int64
        default = 10
    "--nums"
        help = "How many states do you need?"
        arg_type = Int64
        default = 1
    "--PBC"
        help = "Periodic Boundary Condition 1:true 0:false"
        arg_type = Int64
        default = 1
end

parsed_args = parse_args(s)
run()


text:
using NESSie: eye!, pluseye!, isdegenerate, seek, reverseindex, unpack, vertexnormals, cos, cathetus, sign, distance

@testset "eye! and pluseye!" begin
    for T in (Int, testtypes...)
        m = -ones(T, 3, 3)
        pluseye!(m)
        @test typeof(m) == Array{T, 2}
        @test size(m) == (3, 3)
        @test m == [0 -1 -1; -1 0 -1; -1 -1 0]
        pluseye!(m, 2)
        @test typeof(m) == Array{T, 2}
        @test size(m) == (3, 3)
        @test m == [2 -1 -1; -1 2 -1; -1 -1 2]
        eye!(m)
        @test typeof(m) == Array{T, 2}
        @test size(m) == (3, 3)
        @test m == [1 0 0; 0 1 0; 0 0 1]
        eye!(m, 2)
        @test typeof(m) == Array{T, 2}
        @test size(m) == (3, 3)
        @test m == [2 0 0; 0 2 0; 0 0 2]
        m = zeros(T, 2, 3)
        eye!(m)
        @test typeof(m) == Array{T, 2}
        @test size(m) == (2, 3)
        @test m == [1 0 0; 0 1 0]
        eye!(m, 2)
        @test typeof(m) == Array{T, 2}
        @test size(m) == (2, 3)
        @test m == [2 0 0; 0 2 0]
        m = zeros(T, 3, 2)
        eye!(m)
        @test typeof(m) == Array{T, 2}
        @test size(m) == (3, 2)
        @test m == [1 0; 0 1; 0 0]
        eye!(m, 2)
        @test typeof(m) == Array{T, 2}
        @test size(m) == (3, 2)
        @test m == [2 0; 0 2; 0 0]
    end
end

@testset "props and isdegenerate" begin
    @test_throws AssertionError isdegenerate(Triangle([1., 1.], [1., 1., 1.], [1., 1., 1.]))
    @test_throws AssertionError isdegenerate(Triangle([1., 1., 1.], [1., 1.], [1., 1., 1.]))
    @test_throws AssertionError isdegenerate(Triangle([1., 1., 1.], [1., 1., 1.], [1., 1.]))
    for T in testtypes
        # degenerate triangles
        @test_throws AssertionError Triangle(T[0, 0, 0], T[0, 0, 0], T[1, 0, 0])
        @test_throws AssertionError Triangle(T[0, 0, 0], T[0, 0, 1], T[0, 0, 0])
        @test_throws AssertionError Triangle(T[1, 0, 0], T[0, 0, 0], T[0, 0, 0])
        @test_throws AssertionError Triangle(T[0, 1, 0], T[0, 2, 0], T[0, 3, 0])
        @test_throws AssertionError Triangle(T[0, 1, 0], T[1e-11, 1, 0], T[0, 3, 0])
        # simple 2D triangle
        elem = Triangle(T[0, 0, 0], T[0, 0, 3], T[0, 3, 0])
        @test typeof(elem.center) == Vector{T}
        @test length(elem.center) == 3
        @test typeof(elem.normal) == Vector{T}
        @test length(elem.normal) == 3
        @test typeof(elem.distorig) == T
        @test typeof(elem.area) == T
        @test elem.center ≈ [0, 1, 1]
        @test elem.normal ≈ [-1, 0, 0]
        @test elem.distorig ≈ 0
        @test elem.area ≈ 4.5
        # simple 3D triangle
        elem = Triangle(T[3, 0, 0], T[0, 4, 0], T[0, 0, 5])
        @test typeof(elem.center) == Vector{T}
        @test length(elem.center) == 3
        @test typeof(elem.normal) == Vector{T}
        @test length(elem.normal) == 3
        @test typeof(elem.distorig) == T
        @test typeof(elem.area) == T
        @test elem.center ≈ map(T, ([1, 4/3, 5/3]))
        @test elem.normal ≈ map(T, (√769 \ [20, 15, 12]))
        @test elem.distorig ≈ map(T, 60/√769)
        @test elem.area ≈ T(√769/2)
    end
end

@testset "seek" begin
    fname, fh = mktemp()
    write(fh, "Lorem\nipsum\ndolor\nsit\namet.\n")
    seekstart(fh)
    seek(fh, "Foo")
    @test eof(fh)
    seekstart(fh)
    seek(fh, "Foo", false)
    @test eof(fh)
    seekstart(fh)
    seek(fh, "dolor")
    @test readline(fh) == "sit"
    seekstart(fh)
    seek(fh, "dolor", false)
    @test readline(fh) == "dolor"
    close(fh)
    rm(fname)
end

@testset "reverseindex" begin
    for T in testtypes
        v1 = T[1, 2, 3]
        v2 = T[4, 5, 6]
        v3 = T[7, 8, 9]
        d = reverseindex(Vector{T}[])
        @test typeof(d) == Dict{UInt, UInt}
        @test d == Dict()
        d = reverseindex(Vector{T}[v1, v2, v3])
        @test typeof(d) == Dict{UInt, UInt}
        @test length(d) == 3
        @test d[objectid(v1)] == 1
        @test d[objectid(v2)] == 2
        @test d[objectid(v3)] == 3
        d = reverseindex(Vector{T}[v1, v1, v2])
        @test typeof(d) == Dict{UInt, UInt}
        @test length(d) == 2
        @test d[objectid(v1)] == 2
        @test d[objectid(v2)] == 3
        @test_throws KeyError d[objectid(v3)]
    end
end

@testset "unpack" begin
    for T in testtypes
        d = unpack(Vector{T}[])
        @test typeof(d) == Vector{T}
        @test d == []
        d = unpack(Vector{T}[T[1], T[2], T[3]])
        @test typeof(d) == Vector{T}
        @test d == [1, 2, 3]
        d = unpack(Vector{T}[T[1, 2], T[3, 4]])
        @test typeof(d) == Vector{T}
        @test d == [1, 2, 3, 4]
        d = unpack(Vector{T}[T[1, 2, 3], T[4, 5, 6]])
        @test typeof(d) == Vector{T}
        @test d == [1, 2, 3, 4, 5, 6]
        d = unpack(Vector{T}[T[1, 2, 3, 4, 5, 6]])
        @test typeof(d) == Vector{T}
        @test d == [1, 2, 3, 4, 5, 6]
    end
end

@testset "vertexnormals" begin
    for T in testtypes
        nodes    = Vector{T}[T[0, 0, 0], T[0, 0, 3], T[0, 3, 0], T[1, -3, 3]]
        elements = [Triangle(nodes[1], nodes[2], nodes[3]),
                    Triangle(nodes[1], nodes[4], nodes[2])]
        d = vertexnormals(Model{T, Triangle{T}}())
        @test typeof(d) == Vector{Vector{T}}
        @test d == []
        d = vertexnormals(Model(
            Vector{T}[nodes[1], nodes[2], nodes[3]],
            Triangle{T}[elements[1]]
        ))
        @test typeof(d) == Vector{Vector{T}}
        @test length(d) == 3
        @test d[1] ≈ [-1, 0, 0]
        @test d[2] ≈ [-1, 0, 0]
        @test d[3] ≈ [-1, 0, 0]
        d = vertexnormals(Model(
            Vector{T}[nodes[1], nodes[2], nodes[4]],
            Triangle{T}[elements[2]]
        ))
        @test typeof(d) == Vector{Vector{T}}
        @test length(d) == 3
        @test d[1] ≈ map(T, √90 \ [-9, -3, 0])
        @test d[2] ≈ map(T, √90 \ [-9, -3, 0])
        @test d[3] ≈ map(T, √90 \ [-9, -3, 0])
        d = vertexnormals(Model(nodes, elements))
        @test typeof(d) == Vector{Vector{T}}
        @test length(d) == 4
        @test d[1] ≈ map(T, √360 \ [-9 - √90, -3, 0])
        @test d[2] ≈ map(T, √360 \ [-9 - √90, -3, 0])
        @test d[3] ≈ [-1, 0, 0]
        @test d[4] ≈ map(T, √90 \ [-9, -3, 0])
    end
end

@testset "cos" begin
    for T in testtypes
        u = T[1, 0, 0]
        @test typeof(cos(u, T[1, 0, 0])) == T
        @test cos(u, T[1, 0, 0]) ≈ one(T)
        @test cos(u, T[0, 1, 0]) ≈ zero(T)
        @test cos(u, T[-1, 0, 0]) ≈ -one(T)
        @test cos(u, T[0, -1, 0]) ≈ zero(T)
        @test cos(u, T[1, 1, 0]) ≈ T(cos(π/4))
        @test cos(u, T[-1, 1, 0]) ≈ T(cos(3π/4))
        @test cos(u, T[1, 0, 1]) ≈ T(cos(π/4))
    end
end

@testset "meshunion" begin
    for T in testtypes
        # empty system
        empty = Model{T, Tetrahedron{T}}()
        res = meshunion(empty, empty)
        @test typeof(res.nodes) == Vector{Vector{T}}
        @test typeof(res.elements) == Vector{Tetrahedron{T}}
        @test length(res.nodes) == 0
        @test length(res.elements) == 0

        # small system
        nodesΩ = Vector{T}[[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]]
        modelΩ = Model(
            nodesΩ,
            [Tetrahedron(nodesΩ...), Tetrahedron(nodesΩ...)]
        )
        nodesΣ = Vector{T}[[0, 0, 0], [-1, 0, 0], [0, -1, 0], [0, 0, -1]]
        modelΣ = Model(
            nodesΣ,
            [Tetrahedron(nodesΣ...)]
        )
        res = meshunion(modelΩ, modelΣ)
        oids = Set([objectid(e) for e in res.nodes])
        @test typeof(res) == Model{T, Tetrahedron{T}}
        @test length(res.nodes) == 7
        @test length(res.elements) == 3
        for node in nodesΩ ∪ nodesΣ
            @test node ∈ res.nodes
        end
        @test res.elements[1] == modelΩ.elements[1]
        @test res.elements[2] == modelΩ.elements[2]
        @test res.elements[3].v1 === nodesΩ[1]
        @test res.elements[3].v2 === nodesΣ[2]
        @test res.elements[3].v3 === nodesΣ[3]
        @test res.elements[3].v4 === nodesΣ[4]
        for elem in res.elements
            for v in (elem.v1, elem.v2, elem.v3, elem.v4)
                v ∈ oids
            end
        end
    end
end

@testset "obspoints_line" begin
    for T in testtypes
        Ξ = [ξ for ξ in obspoints_line(T[1, 1, 1], T[1, 1, 1], 3)]
        @test typeof(Ξ) == Vector{Vector{T}}
        @test length(Ξ) == 3
        for ξ in Ξ
            @test ξ == T[1, 1, 1]
        end
        Ξ = [ξ for ξ in obspoints_line(T[0, 0, 0], T[2, 2, 2], 3)]
        @test typeof(Ξ) == Vector{Vector{T}}
        @test length(Ξ) == 3
        @test Ξ[1] == T[0, 0, 0]
        @test Ξ[2] == T[1, 1, 1]
        @test Ξ[3] == T[2, 2, 2]
        Ξ = [ξ for ξ in obspoints_line(T[2, 2, 2], T[0, 0, 0], 3)]
        @test typeof(Ξ) == Vector{Vector{T}}
        @test length(Ξ) == 3
        @test Ξ[3] == T[0, 0, 0]
        @test Ξ[2] == T[1, 1, 1]
        @test Ξ[1] == T[2, 2, 2]
    end
end

@testset "obspoints_plane" begin
    for T in testtypes
        Ξ = [collect(Ξ) for Ξ in obspoints_plane(T[1, 1, 1], T[1, 1, 1], T[1, 1, 1], 3, 3)]
        @test typeof(Ξ) == Vector{Vector{Vector{T}}}
        @test length(Ξ) == 3
        for row in Ξ
            @test row == [T[1, 1, 1], T[1, 1, 1], T[1, 1, 1]]
        end
        Ξ = [collect(Ξ) for Ξ in obspoints_plane(T[1, 0, 0], T[0, 0, -1], T[0, 1, 0], 3, 3)]
        @test typeof(Ξ) == Vector{Vector{Vector{T}}}
        @test length(Ξ) == 3
        @test Ξ[1] == [T[0, 0, -1], T[0, 0.5, -0.5], T[0, 1, 0]]
        @test Ξ[2] == [T[0.5, 0, -0.5], T[0.5, 0.5, 0], T[0.5, 1, 0.5]]
        @test Ξ[3] == [T[1, 0, 0], T[1, 0.5, 0.5], T[1, 1, 1]]
        Ξ = [collect(Ξ) for Ξ in obspoints_plane(T[1, 0, 0], T[1, 1, 1], T[0, 1, 0], 3, 3)]
        @test typeof(Ξ) == Vector{Vector{Vector{T}}}
        @test length(Ξ) == 3
        @test Ξ[1] == [T[1, 1, 1], T[0.5, 1, 0.5], T[0, 1, 0]]
        @test Ξ[2] == [T[1, 0.5, 0.5], T[0.5, 0.5, 0], T[0, 0.5, -0.5]]
        @test Ξ[3] == [T[1, 0, 0], T[0.5, 0, -0.5], T[0, 0, -1]]
    end
end

@test_skip cathetus
@test_skip sign
@test_skip distance
@test_skip ddot


text:
struct ThreadSafeVector{T} <: AbstractVector{T}
    V::Matrix{T}
    function ThreadSafeVector{T}(V::Matrix{T}) where T
        @assert size(V, 2) == Threads.nthreads()
        new{T}(V)
    end
end

ThreadSafeVector(V::Matrix) = ThreadSafeVector{eltype(V)}(V)
ThreadSafeVector(v::Vector) = ThreadSafeVector(repmat(v, 1, Threads.nthreads()))

@inline size(v::ThreadSafeVector) = (size(v.V, 1), )
@inline getindex(v::ThreadSafeVector{T}, i::Integer) where T = v.V[i, Threads.threadid()]
@inline setindex!(v::ThreadSafeVector{T}, x, i::Integer) where T = v.V[i, Threads.threadid()] = x

threadsafezeros(::Type{T}, n::Integer) where T = ThreadSafeVector{T}(zeros(T, n, Threads.nthreads()))
threadsafeones(::Type{T}, n::Integer) where T = ThreadSafeVector{T}(ones(T, n, Threads.nthreads()))

abstract type AbstractLowRankMatrix{T} <: AbstractMatrix{T} end

"""
Store the singular value decomposition of a matrix:

    A = UΣV'

"""
struct LowRankMatrix{T} <: AbstractLowRankMatrix{T}
    U::Matrix{T}
    Σ::Diagonal{T, Vector{T}}
    V::Matrix{T}
    temp::ThreadSafeVector{T}
end

LowRankMatrix(U::Matrix{T}, Σ::Diagonal{T}, V::Matrix{T}) where T = LowRankMatrix(U, Σ, V, threadsafezeros(T, length(Σ.diag)))

size(L::LowRankMatrix) = size(L.U, 1), size(L.V, 1)
rank(L::LowRankMatrix) = length(L.Σ.diag)
norm(L::LowRankMatrix) = first(L.Σ)
cond(L::LowRankMatrix) = ifelse(rank(L) < minimum(size(L)), Inf, first(L.Σ)/last(L.Σ))

istriu(L::LowRankMatrix) = false
istril(L::LowRankMatrix) = false
issymmetric(L::LowRankMatrix) = false
ishermitian(L::LowRankMatrix) = false

function getindex(L::LowRankMatrix{T},i::Integer,j::Integer) where T
    ret = zero(T)
    U, Σ, V, r = L.U, L.Σ, L.V, rank(L)
    for k = r:-1:1
        ret += U[i,k]*Σ[k,k]*V[j,k]
    end

    ret
end

function convert(::Type{LowRankMatrix{T}},L::LowRankMatrix) where T
    LowRankMatrix(convert(Matrix{T}, L.U), convert(Diagonal{T}, L.Σ), convert(Matrix{T}, L.V))
end

convert(::Type{LowRankMatrix},A::AbstractMatrix{T})  where T = svdtrunc(A)

function getrank(σ::Vector{T})  where T<:Real
    r = length(σ)
    r == 0 && (return 0)
    tol = r*eps(first(σ))
    while r ≥ 1
        if σ[r] > tol return r end
        r -= 1
    end

    r
end

function svdtrunc(A::AbstractMatrix)
    SVD = svdfact(A)
    r = getrank(SVD[:S])
    LowRankMatrix(SVD[:U][:,1:r], Diagonal(SVD[:S][1:r]), SVD[:V][:,1:r])
end

function lrzeros(::Type{T}, m::Int, n::Int) where T
    U = zeros(T, m, 0)
    Σ = Diagonal(zeros(T, 0))
    V = zeros(T, n, 0)
    LowRankMatrix(U, Σ, V)
end

function (+)(L1::LowRankMatrix{T}, L2::LowRankMatrix{T}) where T
    QRU = qrfact!(hcat(L1.U, L2.U))
    QRV = qrfact!(hcat(L1.V, L2.V))
    SVD = svdfact!(QRU[:R]*Diagonal(vcat(L1.Σ.diag, L2.Σ.diag))*QRV[:R]')
    r = getrank(SVD[:S])

    LowRankMatrix((QRU[:Q]*SVD[:U])[:,1:r], Diagonal(SVD[:S][1:r]), (QRV[:Q]*SVD[:V])[:,1:r])
end

function (-)(L1::LowRankMatrix{T}, L2::LowRankMatrix{T}) where T
    QRU = qrfact!(hcat(L1.U, L2.U))
    QRV = qrfact!(hcat(L1.V, L2.V))
    SVD = svdfact!(QRU[:R]*Diagonal(vcat(L1.Σ.diag, -L2.Σ.diag))*QRV[:R]')
    r = getrank(SVD[:S])

    LowRankMatrix((QRU[:Q]*SVD[:U])[:,1:r], Diagonal(SVD[:S][1:r]), (QRV[:Q]*SVD[:V])[:,1:r])
end

function (*)(L1::LowRankMatrix{T}, L2::LowRankMatrix{T}) where T
    SVD = svdfact!(L1.Σ*L1.V'*L2.U*L2.Σ)
    r = getrank(SVD[:S])

    LowRankMatrix((L1.U*SVD[:U])[:,1:r], Diagonal(SVD[:S][1:r]), (L2.V*SVD[:V])[:,1:r])
end

function (*)(A::AbstractMatrix{T}, L::LowRankMatrix{T}) where T
    QRU = qrfact!(A*L.U)
    SVD = svdfact!(QRU[:R]*L.Σ)
    r = getrank(SVD[:S])

    LowRankMatrix((QRU[:Q]*SVD[:U])[:,1:r], Diagonal(SVD[:S][1:r]), (L.V*SVD[:V])[:,1:r])
end

function (*)(L::LowRankMatrix{T}, A::AbstractMatrix{T}) where T
    QRV = qrfact!(A'*L.V)
    SVD = svdfact!(L.Σ*QRV[:R]')
    r = getrank(SVD[:S])

    LowRankMatrix((L.U*SVD[:U])[:,1:r], Diagonal(SVD[:S][1:r]), (QRV[:Q]*SVD[:V])[:,1:r])
end

function (*)(A::Diagonal{T}, L::LowRankMatrix{T}) where T
    QRU = qrfact!(A*L.U)
    SVD = svdfact!(QRU[:R]*L.Σ)
    r = getrank(SVD[:S])

    LowRankMatrix((QRU[:Q]*SVD[:U])[:,1:r], Diagonal(SVD[:S][1:r]), (L.V*SVD[:V])[:,1:r])
end

function (*)(L::LowRankMatrix{T}, A::Diagonal{T}) where T
    QRV = qrfact!(A'*L.V)
    SVD = svdfact!(L.Σ*QRV[:R]')
    r = getrank(SVD[:S])

    LowRankMatrix((L.U*SVD[:U])[:,1:r], Diagonal(SVD[:S][1:r]), (QRV[:Q]*SVD[:V])[:,1:r])
end

(*)(B::T, L::LowRankMatrix{T}) where {T<:Number} = LowRankMatrix(L.U, B*L.Σ, L.V)
(*)(L::LowRankMatrix{T}, B::T) where {T<:Number} = LowRankMatrix(L.U, L.Σ*B, L.V)
(/)(L::LowRankMatrix{T}, B::T) where {T<:Number} = LowRankMatrix(L.U, L.Σ/B, L.V)
function broadcast(::typeof(*), B::T, L::LowRankMatrix{T}) where T
    LowRankMatrix(L.U, B.*L.Σ, L.V)
end
function broadcast(::typeof(*), L::LowRankMatrix{T}, B::T) where T
    LowRankMatrix(L.U, L.Σ.*B, L.V)
end
function broadcast(::typeof(/), L::LowRankMatrix{T}, B::T) where T
    LowRankMatrix(L.U, L.Σ./B, L.V)
end


text:
struct NormalProblem{T}
    func::Function
    domain::T
end

struct DifficultProblem{T}
    func::Function
    domain::T
end

prblm_1 = NormalProblem{IntervalBox}((x,y,z)->x^3+y^2+z, IntervalBox(2..3, 5..6, 4..8))
prblm_2 = NormalProblem{IntervalBox}((x,y,z)->((x^3+y^2)^4+z^2, IntervalBox(2..3, 5..6, 4..8))
prblm_3 = NormalProblem{IntervalBox}((x,y,z)->((x^3+y^2)^4+z^2, IntervalBox(2..3, 5..6, 4..8))


text:
using OrdinaryDiffEq, RecursiveArrayTools, Test, StaticArrays, DiffEqCallbacks

f = function (u,p,t)
  - u + sin(-t)
end

prob = ODEProblem(f,1.0,(0.0,-10.0))

condition= function (u,t,integrator) # Event when event_f(u,t,k) == 0
  - t - 2.95
end

affect! = function (integrator)
  integrator.u = integrator.u + 2
end

callback = ContinuousCallback(condition,affect!)

sol = solve(prob,Tsit5(),callback=callback)
@test length(sol) < 20

condition= function (out, u,t,integrator) # Event when event_f(u,t,k) == 0
  out[1] = - t - 2.95
end

affect! = function (integrator, idx)
  if idx == 1
    integrator.u = integrator.u + 2
  end
end

callback = VectorContinuousCallback(condition,affect!,1)

sol = solve(prob,Tsit5(),callback=callback)

f = function (du,u,p,t)
  du[1] = - u[1] + sin(t)
end

prob = ODEProblem(f,[1.0],(0.0,10.0))

condition= function (u,t,integrator) # Event when event_f(u,t,k) == 0
  t - 2.95
end

affect! = function (integrator)
  integrator.u = integrator.u .+ 2
end

callback = ContinuousCallback(condition,affect!)

sol = solve(prob,Tsit5(),callback=callback,abstol=1e-8,reltol=1e-6)

condition= function (out, u,t,integrator) # Event when event_f(u,t,k) == 0
  out[1] = t - 2.95
end

affect! = function (integrator, idx)
  if idx == 1
    integrator.u = integrator.u .+ 2
  end
end

callback = VectorContinuousCallback(condition,affect!,1)

sol = solve(prob,Tsit5(),callback=callback,abstol=1e-8,reltol=1e-6)

f = function (du,u,p,t)
  du[1] = u[2]
  du[2] = -9.81
end

condition= function (u,t,integrator) # Event when event_f(u,t,k) == 0
  u[1]
end

affect! = nothing
affect_neg! = function (integrator)
  integrator.u[2] = -integrator.u[2]
end

callback = ContinuousCallback(condition,affect!,affect_neg!,interp_points=100)

u0 = [50.0,0.0]
tspan = (0.0,15.0)
prob = ODEProblem(f,u0,tspan)

sol = solve(prob,Tsit5(),callback=callback,adaptive=false,dt=1/4)

condition= function (out, u,t,integrator) # Event when event_f(u,t,k) == 0
  out[1] = u[1]
end

affect! = nothing
affect_neg! = function (integrator, idx)
  if idx == 1
    integrator.u[2] = -integrator.u[2]
  end
end

vcb = VectorContinuousCallback(condition,affect!,1,interp_points=100)

sol = solve(prob,Tsit5(),callback=vcb,adaptive=false,dt=1/4)

condition_single = function (u,t,integrator) # Event when event_f(u,t,k) == 0
  u
end

affect! = nothing
affect_neg! = function (integrator)
  integrator.u[2] = -integrator.u[2]
end

callback_single = ContinuousCallback(condition_single,affect!,affect_neg!,interp_points=100,idxs=1)

u0 = [50.0,0.0]
tspan = (0.0,15.0)
prob = ODEProblem(f,u0,tspan)

sol = solve(prob,Tsit5(),callback=callback_single,adaptive=false,dt=1/4)
sol = solve(prob,Tsit5(),callback=callback_single,save_everystep=false)
t = sol.t[end÷2] # this is the callback time point
sol = solve(prob,Tsit5(),callback=callback_single,saveat=t)
@test count(x->x==t, sol.t) == 2
sol = solve(prob,Tsit5(),callback=callback_single,saveat=t-eps(t))
@test count(x->x==t, sol.t) == 2
# check interpolation @ discontinuity
@test sol(t,continuity=:right)[2] > 0
@test sol(t,continuity=:left)[2] < 0

#plot(sol,denseplot=true)

condition_single = function (out,u,t,integrator) # Event when event_f(u,t,k) == 0
  out[1] = u[1]
end

affect! = nothing
affect_neg! = function (integrator,idx)
  if idx == 1
    integrator.u[2] = -integrator.u[2]
  end
end

callback_single = VectorContinuousCallback(condition_single,affect!,affect_neg!,1,interp_points=100)

sol = solve(prob,Tsit5(),callback=callback_single,adaptive=false,dt=1/4)
sol = solve(prob,Tsit5(),callback=callback_single,save_everystep=false)
t = sol.t[end÷2] # this is the callback time point
sol = solve(prob,Tsit5(),callback=callback_single,saveat=t)
@test count(x->x==t, sol.t) == 2
sol = solve(prob,Tsit5(),callback=callback_single,saveat=t-eps(t))
@test count(x->x==t, sol.t) == 2
# check interpolation @ discontinuity
@test sol(t,continuity=:right)[2] > 0
@test sol(t,continuity=:left)[2] < 0

sol = solve(prob,Vern6(),callback=callback)
sol = solve(prob,Vern6(),callback=vcb)
#plot(sol,denseplot=true)
sol = solve(prob,BS3(),callback=vcb)
sol = solve(prob,BS3(),callback=callback)

sol33 = solve(prob,Vern7(),callback=callback)
sol33 = solve(prob,Vern7(),callback=vcb)

bounced = ODEProblem(f,sol[8],(0.0,1.0))
sol_bounced = solve(bounced,Vern6(),callback=callback,dt=sol.t[9]-sol.t[8])
#plot(sol_bounced,denseplot=true)
sol_bounced(0.04) # Complete density
@test maximum(maximum.(map((i)->sol.k[9][i]-sol_bounced.k[2][i],1:length(sol.k[9])))) == 0

sol2= solve(prob,Vern6(),callback=callback,adaptive=false,dt=1/2^4)
#plot(sol2)

sol2= solve(prob,Vern6())

sol3= solve(prob,Vern6(),saveat=[.5])

## Saving callback

condition = function (u,t,integrator)
  true
end
affect! = function (integrator) end

save_positions = (true,false)
saving_callback = DiscreteCallback(condition,affect!,save_positions=save_positions)

sol4 = solve(prob,Tsit5(),callback=saving_callback)

@test sol2(3) ≈ sol(3)

affect! = function (integrator)
  u_modified!(integrator,false)
end
saving_callback2 = DiscreteCallback(condition,affect!,save_positions=save_positions)
sol4 = solve(prob,Tsit5(),callback=saving_callback2)

cbs = CallbackSet(saving_callback,saving_callback2)
sol4_extra = solve(prob,Tsit5(),callback=cbs)

@test length(sol4_extra) == 2length(sol4) - 1

condition= function (u,t,integrator)
  u[1]
end

vcondition= function (out,u,t,integrator)
  out[1] = u[1]
end

affect! = function (integrator, retcode = nothing)
  if retcode === nothing
    terminate!(integrator)
  else
    terminate!(integrator, retcode)
  end
end

vaffect! = function (integrator, idx, retcode = nothing)
  if idx == 1
    if retcode === nothing
      terminate!(integrator)
    else
      terminate!(integrator, retcode)
    end
  end
end

terminate_callback = ContinuousCallback(condition,affect!)
custom_retcode_callback = ContinuousCallback(condition,x->affect!(x,:Custom))
vterminate_callback = VectorContinuousCallback(vcondition,vaffect!,1)
vcustom_retcode_callback = VectorContinuousCallback(vcondition,(x,idx)->vaffect!(x,idx,:Custom),1)

tspan2 = (0.0,Inf)
prob2 = ODEProblem(f,u0,tspan2)

sol5 = solve(prob2,Tsit5(),callback=terminate_callback)
sol5_1 = solve(prob2,Tsit5(),callback=custom_retcode_callback)

@test sol5.retcode == :Terminated
@test sol5_1.retcode == :Custom
@test sol5[end][1] < 3e-12
@test sol5.t[end] ≈ sqrt(50*2/9.81)

sol5 = solve(prob2,Tsit5(),callback=vterminate_callback)
sol5_1 = solve(prob2,Tsit5(),callback=vcustom_retcode_callback)

@test sol5.retcode == :Terminated
@test sol5_1.retcode == :Custom
@test sol5[end][1] < 3e-12
@test sol5.t[end] ≈ sqrt(50*2/9.81)

affect2! = function (integrator)
  if integrator.t >= 3.5
    terminate!(integrator)
  else
    integrator.u[2] = -integrator.u[2]
  end
end

vaffect2! = function (integrator,idx)
  if idx == 1
    if integrator.t >= 3.5
      terminate!(integrator)
    else
      integrator.u[2] = -integrator.u[2]
    end
  end
end

terminate_callback2 = ContinuousCallback(condition,nothing,affect2!,interp_points=100)
vterminate_callback2 = VectorContinuousCallback(vcondition,nothing,vaffect2!,1,interp_points=100)

sol5 = solve(prob2,Vern7(),callback=terminate_callback2)

@test sol5[end][1] < 1.3e-10
@test sol5.t[end] ≈ 3*sqrt(50*2/9.81)

sol5 = solve(prob2,Vern7(),callback=vterminate_callback2)

@test sol5[end][1] < 1.3e-10
@test sol5.t[end] ≈ 3*sqrt(50*2/9.81)

condition= function (u,t,integrator) # Event when event_f(u,t,k) == 0
  t-4
end

affect! = function (integrator)
  terminate!(integrator)
end

terminate_callback3 = ContinuousCallback(condition,affect!,interp_points=1000)

bounce_then_exit = CallbackSet(callback,terminate_callback3)

sol6 = solve(prob2,Vern7(),callback=bounce_then_exit)

@test sol6[end][1] > 0
@test sol6[end][1] < 100
@test sol6.t[end] ≈ 4

# Test ContinuousCallback hits values on the steps
t_event = 100.0
f_simple(u,p,t) = 1.00001*u
event_triggered = false
condition_simple(u,t,integrator) = t_event-t
function affect_simple!(integrator)
  global event_triggered
  event_triggered = true
end
cb = ContinuousCallback(condition_simple, nothing, affect_simple!)
prob = ODEProblem(f_simple, [1.0], (0.0, 2.0*t_event))
sol = solve(prob,Tsit5(),callback=cb, adaptive = false, dt = 10.0)
@test event_triggered

# https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/issues/328
ode = ODEProblem((du, u, p, t) -> (@. du .= -u), ones(5), (0.0, 100.0))
sol = solve(ode, AutoTsit5(Rosenbrock23()), callback=TerminateSteadyState())
sol1 = solve(ode, Tsit5(), callback=TerminateSteadyState())
@test sol.u == sol1.u

# DiscreteCallback
f = function (du,u,p,t)
  du[1] = -0.5*u[1] + 10
  du[2] = -0.5*u[2]
end

u0 = [10,10.]
tstop = [5.;8.]
prob = ODEProblem(f, u0, (0, 10.))
condition = (u,t,integrator) -> t in tstop
affect! = (integrator) -> integrator.u .= 1.0
save_positions = (true,true)
cb = DiscreteCallback(condition, affect!, save_positions=save_positions)
sol1 = solve(prob,Tsit5(),callback = cb,tstops=tstop,saveat=tstop)
@test count(x->x==tstop[1], sol1.t) == 2
@test count(x->x==tstop[2], sol1.t) == 2
sol2 = solve(prob,Tsit5(),callback = cb,tstops=tstop,saveat=prevfloat.(tstop))
@test count(x->x==tstop[1], sol2.t) == 2
@test count(x->x==tstop[2], sol2.t) == 2


text:
# small scrip to compute the solution of Lippmann-Schwinger equation
# We test the types introduced in FastConvolution.jl
# we test that the application is fast and that the construction
# is performed fast.

using PyPlot
using IterativeSolvers
using SpecialFunctions

include("../src/FastConvolution.jl")
include("../src/Preconditioner.jl")

# setting the number of threads for the FFT and BLAS
# libraries (please set them to match the number of
# physical cores in your system)
FFTW.set_num_threads(4);
BLAS.set_num_threads(4);

#Defining Omega
h = 0.005
k = 1/(h)

# size of box
a = 1
x = collect(-a/2:h:a/2)
y = collect(-a/2:h:a/2)
(n,m) = length(x), length(y)
N = n*m
X = repmat(x, 1, m)[:]
Y = repmat(y', n,1)[:]
# we solve \triangle u + k^2(1 + nu(x))u = 0

# We use the modified quadrature in Ruan and Rohklin
(ppw,D) = referenceValsTrapRule();
D0 = D[1];

# Defining the smooth perturbation of the slowness
nu(x,y) = 0.3*exp.(-40*(x.^2 + y.^2)).*(abs.(x).<0.48).*(abs.(y).<0.48);

## You can choose between Duan Rohklin trapezoidal quadrature
#fastconv = buildFastConvolution(x,y,h,k,nu)

# or Greengard Vico Quadrature (this is not optimized and is 2-3 times slower)
fastconv = buildFastConvolution(x,y,h,k,nu, quadRule = "Greengard_Vico");

# assembling the sparsifiying preconditioner
@time As = buildSparseA(k,X,Y,D0, n ,m);

# assembling As*( I + k^2G*nu)
@time Mapproxsp = As + k^2*(buildSparseAG(k,X,Y,D0, n ,m)*spdiagm(nu(X,Y)));

# defining the preconditioner
precond = SparsifyingPreconditioner(Mapproxsp, As)#, solverType="MKLPardiso")

# building the RHS from the incident field
u_inc = exp.(k*im*X);
rhs = -k^2*FFTconvolution(fastconv, nu(X,Y).*u_inc) ;

#rhs = u_inc;

# allocating the solution
u = zeros(Complex128,N);

# solving the system using GMRES
@time info =  gmres!(u, fastconv, rhs, Pl=precond, log = true)
#println(info[2].residuals[:])

# plotting the solution
figure(1)
clf()
imshow(real(reshape(u+u_inc,n,m)))

max_stars_repo_name: Forgotten/Lippmann-Schingwer

max_issues_repo_name: Forgotten/Lippmann-Schingwer

max_forks_repo_name: Forgotten/Lippmann-Schingwer


text:
## Documentation and Examples

# Documentation is entered in files like this:
# @sjdoc SomeHead "SomeHead does this."

# @sjseealso(Sym, Sym1, Sym2, etc.)
# Sets Seealso's for Sym to Sym1, ...

# @sjseealso_group(Sym1,Sym2,...) set's each symbols see alsos to
# all of the others.

# Examples are entered like this. All these pairs together constitute and example
# @sjexamp(SomeHead,
#    "first input", "first output",
#    "second ....
# Examples are printed along with the document.

# Repeating @sjexamp pushes more examples onto an array for SomeHead

# do_example evaluates the examples (only input strings) and
# prints input and output, and pushes input strings onto the history
# so the user can play with them.

# Single doc string for a symbol. Retrieved by:
# symata>? SomeHead

const SJEXAMPLES = Dict{Symbol,Array{Any,1}}()
const SJSEEALSO = Dict{Symbol,Array{Any,1}}()

## Following moved to mxpr.jl
#const SJDOCS = Dict{Symbol,Any}()
# function sjdocfun(sym,str)
#     SJDOCS[sym] = Markdown.parse(str)
#     nothing
# end
# macro sjdoc(sym,str)
#     SJDOCS[sym] = Markdown.parse(str)
#     nothing
# end

# Called from symataevaluate, which parses cli input to see if we have a help query
# or an expression. There is no io stream specified when printing here
function check_doc_query(ex)
    if isa(ex,Expr) && ex.head == :tuple && ex.args[1] == MAGIC_HELP_QUERY_SYMBOL
        if length(ex.args) == 1
            println("Documentation is available for these symbols.")
            print(wrapout(list_documented_symbols()))
            return true
        end
        q = ex.args[2]
        print_doc(q)
        return true
    else
        return false
    end
end

print_doc(s::AbstractString) = print_doc(Symbol(s))

function print_doc(qs...)
    length(qs) == 0 && (println("Try Help(sym) for these symbols."); return list_documented_symbols())
    for q in qs
        check_autoload(q)
        if haskey(SJDOCS,q)
            display(SJDOCS[q])
            println()
            format_see_alsos(q)
            format_sjexamples(q)
        elseif ! write_definition_no_attributes(q)
            println("No documentation for '", string(q),"'.")
        end
        as = get_attributes(q)
        if length(as) > 0
            symprintln("\n Attributes(", string(q), ") = ", mxpr(:List,as...))
        end
        if getkerneloptions(:show_sympy_docs) print_sympy_doc(q) end
    end
    return Null
end

# FIXME. We need to have all doc printing functions return strings.
# Then, we either print them or pass them to a front end.
# Commented out parts remain TODO
function retrieve_doc(qs...)
#    length(qs) == 0 && (println("Try Help(sym) for these symbols."); return list_documented_symbols())
    for iq in qs
        q = Symbol(iq)
        check_autoload(q)
        if haskey(SJDOCS,q)
            return(SJDOCS[q])
#            format_see_alsos(q)
#            format_sjexamples(q)
        else
            return("No documentation for '" * string(q) * "'.")
        end
    end
end

function print_sympy_doc(sjsymin::SymString)
    sjsym = Symbol(sjsymin)
    if have_pyfunc_symbol(sjsym)
        println("\nSymPy documentation")
        pysym = lookup_pyfunc_symbol(sjsym)
        spysym = string(pysym)
        printcom = "println(sympy[:$spysym][:__doc__])"
        try eval(Meta.parse(printcom))
        catch
            Null
        end
    end
    return Null
end

function print_matching_topics(r::Regex)
    i = 0
    lastt = :none
    lastdoc = :none
    for (t, doc) in SJDOCS
        if occursin(r, string(t)) || occursin(r, string(doc))
            i += 1
            println(t)
            lastt = t
            lastdoc = doc
        end
    end
    for sym in get_builtins()  # check builtin symbols with no doc string.
        if (! haskey(SJDOCS, sym)) && occursin(r, string(sym))
            println(string(sym))
        end
    end
    if i == 1
        if occursin(r, string(lastt)) # only print whole doc if we match name
            println(lastdoc)
        end
    end
    return nothing
end

function print_all_docs()
    syms = sort!(collect(keys(SJDOCS)))
    for x in syms
        println("##### ", x)
        print(SJDOCS[x])
        println()
        format_see_alsos(x)
        println()
        if haskey(SJEXAMPLES,x)
            exs = SJEXAMPLES[x]
            println("Examples")
            for strs in exs
                println()
                println("```")
                format_sjexample(strs)
                println("```")
            end
            println()
        end
        println()
    end
    return Null
end

function list_documented_symbols()
    syms = sort!(collect(keys(SJDOCS)))
    len = length(syms)
    nargs = newargs()
    @inbounds for i in 1:len
        if syms[i] == :ans continue end
        push!(nargs,syms[i])
    end
    return MListA(nargs)
end

#### See Alsos

# Set seealsos. clobbers existing.
macro sjseealso(source_sym, targ_syms...)
    SJSEEALSO[source_sym] = Any[targ_syms...]
end

# each sym seealso's all others
# This does not clobber existing seealsos
# FIXME: Do not add symbols more than once.
macro sjseealso_group(syms...)
    local na
    for src in syms
        if ! haskey(SJSEEALSO,src)
            na = Array{Any}(undef, 0)
            SJSEEALSO[src] = na
        else
            na = SJSEEALSO[src]
        end
        for targ in syms
            if src == targ continue end
            if ! in(targ,na) push!(na,targ) end
        end
    end
end

## Print see alsos for sym
function format_see_alsos(sym)
    ! haskey(SJSEEALSO,sym) && return
    sas = copy(SJSEEALSO[sym])
    length(sas) == 0 && return
    sort!(sas)
    print("See also ")
    len = length(sas)
    if len == 1
        println(string(sas[1]),".")
        return nothing
    elseif len == 2
        println(string(sas[1]), " and ", string(sas[2]),".")
        return nothing
    else
        for i in 1:len-1
            print(string(sas[i]),", ")
        end
        println("and ", string(sas[end]),".")
    end
    return nothing
end

#### Examples

"""
    macro sjexamp(sym,strs...)

define examples for a symbol (representing a Head). This is used
in source files to define examples for a symbol (representing a Head).
Maybe this should be a function
"""
macro sjexamp(sym,strs...)
    if ! haskey(SJEXAMPLES,sym)
        SJEXAMPLES[sym] = Array{Any}(undef, 0)
    end
    exs = SJEXAMPLES[sym]
    ar = Array{Any}(undef, 0)
    for i in 1:length(strs)
        push!(ar,eval(strs[i]))
    end
    push!(exs,ar)
end

"""
    format_sjexamples(sym)

Print Symata examples associated with a symbol.
"""
function format_sjexamples(sym)
    if haskey(SJEXAMPLES,sym)
        exs = SJEXAMPLES[sym]
        println("Examples")
        for strs in exs
            println()
            format_sjexample(strs)
        end
    end
    return nothing
end

"""
    format_sjexample(lines)

Format just one example for printing with docs. Do not evaluate.
Each item in example is either an explanatory string, or a 2-tuple of input and output.
"""
function format_sjexample(lines)
    local ins, outs, expl
    for i in 1:length(lines)
        ln = lines[i]
        if isa(ln,AbstractString)
            println(ln)
        else
            expl = ""
            if length(ln) == 2
                (ins,outs) = ln
            else
                (ins,outs,expl) = ln
            end
            if length(expl) > 0
                println("# ",expl)
            end
            println("symata> ", ins)
            if length(outs) > 0 println(outs) end
            if i < length(lines)
                println()
            end
        end
    end
    return nothing
end

# evaluate example code, and push corresponding input strings onto
# the history, so the user can modify the example.
function do_example(lines)
    len = length(lines)
    local ins, outs, expl
    for i in 1:len
        ln = lines[i]
        if isa(ln,AbstractString)
            println(ln)
        else
            expl = ""
            if length(ln) == 2
                (ins,outs) = ln
            else
                (ins,outs,expl) = ln
            end
            if length(expl) > 0
                println("# ",expl)
            end
            println("symata> ", ins)
            ex = MacroTools.rmlines(Meta.parse(ins))
            res = (eval(Expr(:macrocall, Symbol("@symfull"), ex)))
            if res != nothing println(res) end
            sj_add_history(ins)
            if i < len
                println()
            end
        end
    end
    return nothing
end

## evaluate the nth example for symbol sym
function do_example_n(sym,n)
    if haskey(SJEXAMPLES,sym)
        exs = SJEXAMPLES[sym]
        if length(exs) < n
            println("There is no example $n for ", string(sym), ".")
            return nothing
        else
            do_example(exs[n])
        end
    else
        println("There are no examples for ", string(sym), ".")
    end
    return nothing
end

## evaluate all examples for symbol sym
function do_examples(sym)
    if haskey(SJEXAMPLES,sym)
        exs = SJEXAMPLES[sym]
        for strs in exs
            do_example(strs)
            println()
            println("--------------------------")
            println()
        end
    end
    return nothing
end

## Information on hierarchy of fields:

# julia> Base.active_repl.interface.modes[2]
# "Prompt(\"symata> \",...)"

# julia> typeof(Base.active_repl.interface.modes[2])
# Base.LineEdit.Prompt

# julia> typeof(Base.active_repl.interface.modes[2].hist)
# Base.REPL.REPLHistoryProvider

# Push a string onto history in symata (symjulia) mode
# modified from add_history in REPL.jl
function sj_add_history(s::AbstractString)
    sj_add_history(symata_repl_history(), s)
end

function sj_add_history(hist::REPL.REPLHistoryProvider, s::AbstractString)
    str = rstrip(String(s))
    isempty(strip(str)) && return ## TODO: efficiency
    mode = :symata
    length(hist.history) > 0 &&  # we could be more clever and not push the same example sequence twice
    mode == hist.modes[end] && str == hist.history[end] && return
    push!(hist.modes, mode)
    push!(hist.history, str)
    hist.history_file == nothing && return
    entry = """
    # time: $(Libc.strftime("%Y-%m-%d %H:%M:%S %Z", time()))
    # mode: $mode
    $(Base.replace(str, r"^"ms => "\t"))
    """
    # TODO: write-lock history file
    seekend(hist.history_file)
    print(hist.history_file, entry)
    flush(hist.history_file)
    return nothing
end

### Help

@sjdoc Help """
    Help(sym), Help(\"sym\"), or ? topic

print documentation for the symbol `sym`. For example `Help(Expand)`.
Due to parsing restrictions at the repl, for some topics, the input must be a string.

    h"topic" gives a case-insensitive regular expression search.

In the REPL, or Jupyter notebook, hit TAB to see all the available completions.

    Help(All => True)

print all of the documentation.

    Help(regex)

prints a list of topics whose documentation text matches the
regular expression regex. For example `Help(r\"Set\"i)` lists all topics that
match "Set" case-independently.
"""

@mkapprule Help  nodefault => true

@doap Help() = print_doc("Help")

function do_Help(mx::Mxpr{:Help}, r::Mxpr{:Rule})
    if r[1] == :All && r[2] == true
        print_all_docs()
    end
    return Null
end

do_Help(mx::Mxpr{:Help},args...) =  print_doc(args...)

function do_Help(mx::Mxpr{:Help}, r::T) where T<:Regex
    print_matching_topics(r)
    return nothing
end

### Example

@sjdoc Example """
    Example(s)

run (evaluates) all examples for the symbol `s`, typically a function or variable.
Input, output, and comments are displayed. Input strings
for the example are pushed onto the terminal history so they can be retrieved and
edited and re-evaluated.

    Example(s,n)

run the `n`th example for symbol `s`. When viewing documentation strings via `?`
the examples are printed along with the documentation string, but are not evaluated.

    Example()

Returns a list of all example topics.
"""

@mkapprule Example  :nargs => 0:2

@doap Example() = mxprcf(:List,Any[sort(collect(keys(SJEXAMPLES)))...])
@doap Example(topic) = do_examples(mx[1])  ## <-- mx[1] is topic ?
@doap Example(topic, n::Int) = do_example_n(mx[1],n)


text:
using PUDLData
using Documenter

makedocs(;
    modules=[PUDLData],
    repo="https://github.com/nickrobinson251/PUDLData.jl/blob/{commit}{path}#{line}",
    sitename="PUDLData.jl",
    format=Documenter.HTML(; prettyurls=false),
    pages=[
        "Home" => "index.md",
    ],
    strict=true,
    checkdocs=:exports,
)

deploydocs(;
    repo="github.com/nickrobinson251/PUDLData.jl",
)


text:
export MetropolisSampler

abstract type MCMCRule end
mutable struct MetropolisSampler{T<:MCMCRule} <: MonteCarloSampler
    chain_length::Int
    burn_length::Int
    passes::Int
    seed::UInt
    rule::T
end

"""
    MetropolisSampler(rule, chain_length, passes; burn=0, seed=rand)

Constructs a Metropolis-Hastings sampler which samples Markov chains of length
`chain_length + burn`, ignoring the first `burn` samples. Transition rules are
specified by the rule `rule`. To reduce auto-correlation, `passes` number of
metropolis-hastings steps are performed for each sample returned (minimum 1).

Effectively, this means that the chain is actually `passes * (chain_lenght + burn)`
long, but only 1 every passes elements are stored and used to compute expectation
values.

Initial seed can be set bu specifying `seed`.
"""
function MetropolisSampler(rule, chain_length, passes; burn=0, seed=rand(UInt))
    @assert passes > 0
    @assert chain_length > 0

    if iseven(passes)
        @warn """
            LocalRule sampling with an even number of switches per step
            is known to be non-ergodic, because it can potentially never
            change the configuration.
            Setting the number of switches to $(passes+1).
            """
        passes += 1
    end
    MetropolisSampler{typeof(rule)}(chain_length, burn, passes, seed, rule)
end

## Show method
function Base.show(io::IO, s::MetropolisSampler)
    burn_time = @sprintf "%2.2f" (s.burn_length/(s.burn_length+s.chain_length))*100

    println(io, "Metropolis Sampler:")
    println(io, "\t- chain_length     : $(s.chain_length)")
    println(io, "\t- burn_length      : $(s.burn_length) ($burn_time% warmup time)")
    println(io, "\t- passes           : $(s.passes)")
    print(io,   "\t- seed             : $(s.seed)")
end

## Cache
mutable struct MetropolisSamplerCache{A<:AbstractRNG,RC,H,S,RV,V,M} <: SamplerCache{MetropolisSampler}
    rng::A
    loc_chain_length::Int
    steps_done::Int
    passes_done::Int
    passes_accepted::Int
    rule_cache::RC

    hilb::H
    σ_old::S

    logψ_σ::RV
    logψ_σp::RV
    log_prob_bias::RV

    ψtmp::V
    prob::RV
    mask::M
end

function _sampler_cache(s::MetropolisSampler, v, hilb, net, par_cache)
    loc_chain_length = Int(ceil(s.chain_length/num_workers(par_cache)))
    loc_seed         = worker_local_seed(s.seed, par_cache)

    logψ_σ   = real(out_similar(net))
    logψ_σp  = real(out_similar(net))
    log_bias = real(out_similar(net)) .= 0.0

    out      = out_similar(net)
    prob     = logψ_σp - logψ_σ
    mask     = similar(logψ_σ, Bool)

    rng      = build_rng_generator_T(prob, loc_seed)

    MetropolisSamplerCache(rng, loc_chain_length, 0, 0, 0,
                      RuleSamplerCache(s.rule, s, v, net, par_cache),
                      hilb, deepcopy(v),
                      logψ_σ, logψ_σp, log_bias,
                      out, prob, mask)
end

# Construct the cache of a specific rule.
# By default rules have no cache.
RuleSamplerCache(r::MCMCRule, s, v, net, par_cache) = nothing
rule(s::MetropolisSampler) = s.rule
rule_cache(sc::MetropolisSamplerCache) = sc.rule_cache

function init_sampler!(s::MetropolisSampler, net, σ,
                       c::MetropolisSamplerCache)
    c.steps_done = - s.burn_length + 1
    c.passes_done = 0
    c.passes_accepted = 0
    rand!(c.rng, σ, c.hilb)
    init_sampler_rule_cache!(rule_cache(c), s, net, σ, c)
    c.log_prob_bias .= 0

    while c.steps_done <= 0
        samplenext!(σ, σ, s, net, c)
    end

    return c
end

init_sampler_rule_cache!(rc, s, net, σ, c) =
    nothing

chain_length(s::MetropolisSampler, c::MetropolisSamplerCache) = c.loc_chain_length

done(s::MetropolisSampler, σ, c) = c.steps_done >= chain_length(s, c)

function samplenext!(σ_out::T, σ_in::T, s::MetropolisSampler,
                     net::NeuralNetwork, c) where {T<:Union{AStateBatch, ADoubleStateBatch}}
    # Check termination condition, and return if verified
    done(s, σ, c) && return false

    # Caches for batched state output
    logψ_σ  = c.logψ_σ
    logψ_σp = c.logψ_σp

    # Copy the old state in the output and tmp cache
    # Only do this when when input is not same as ouput
    # to avoid aliasing issues (especially on GPU)
    σ_out === σ_in || state_copy!(σ_out, σ_in)

    # Compute the old value
    ψtmp     = log_prob_ψ!(logψ_σ, c.ψtmp, net, σ_in)

    for i=1:s.passes
        state_copy!(c.σ_old, σ_out)
        # Apply the transition rule
        propose_step!(σ_out, s, net, c, rule_cache(c))

        ψtmp     = log_prob_ψ!(logψ_σp, c.ψtmp, net, σ_out)

        rand!(c.rng, c.prob)
        c.prob .-= exp.(logψ_σp .- logψ_σ .+ c.log_prob_bias)

        # the mask encodes the non-switched states
        # when rand >= prob
        # ex : prob = 0.1 => 0.1..1.0 >= 0.1
        c.mask .= c.prob .>= 0

        # Copy the old configurations for chains that did not swithc
        state_copy!(σ_out, c.σ_old, c.mask)
        # Copy the new log_σ values for states that switched
        state_copy_invertmask!(logψ_σ, logψ_σp, c.mask)

        # Store number of accepted moves
        c.passes_accepted += length(c.prob) - sum(c.mask)
        c.passes_done     += length(c.prob)
    end
    c.steps_done += 1
    return true
end

# General implementation of Rule-step proposal on a batch, serially applying
# the rule on all batches.
@inline function propose_step!(
                       σp::Union{AStateBatch,ADoubleStateBatch},
                       s::MetropolisSampler,
                       net::NeuralNetwork, c, rc)
    for σ=states(σp)
        propose_step!(σ, s, net, c, rc)
    end
end

## Show method
function Base.show(io::IO, c::MetropolisSamplerCache)
    println(io, "MetropolisSamplerCache{...}:")
    println(io, "\t- Space                 : $(c.hilb)")
    println(io, "\t- Stored chain length   : $(c.loc_chain_length)")
    println(io, "\t- MCMC steps done       : $(c.steps_done)")
    println(io, "\t- Passes done           : $(c.passes_done)")

    acceptance_ratio = @sprintf "%2.2f" (c.passes_accepted/c.passes_done)*100
    print(io, "\t- Passes accepted       : $(c.passes_accepted) ($acceptance_ratio %)")
end


text:
using Dash, DashHtmlComponents, DashCoreComponents, DashBootstrapComponents

include("learn_PlotlyJS.jl")

function show_plot(app, plt)
    if isnothing(app)
        app = dash(external_stylesheets=[dbc_themes.SPACELAB])
    end

    @show typeof(plt)

    app.layout = html_div() do #dbc_container() do
        dcc_graph(
            id = "plotly",
            figure = plt,
        )
    end

    return app
end

app = dash(external_stylesheets=[dbc_themes.SPACELAB])

#app = show_plot(app, clustering_alpha_shapes())
app = show_plot(app, linescatter5())
#app = show_plot(app, scatter_3d())

run_server(
    app,
    "0.0.0.0",
    8050,
    debug=true, # enables hot reload and more
)


text:
using CompScienceMeshes
using BEAST

Γ1 = readmesh(joinpath(dirname(pathof(BEAST)),"../examples/sphere2.in"))
Γ2 = CompScienceMeshes.translate(Γ1, [10,0,0])

SL = Maxwell3D.singlelayer(wavenumber=1.0)
X1 = raviartthomas(Γ1)
X2 = raviartthomas(Γ2)

x1 = refspace(X1)
x2 = refspace(X2)

function BEAST.quaddata(op::typeof(SL), tref::typeof(x1), bref::typeof(x2),
    tels, bels, qs::BEAST.DoubleNumQStrat)

    qs = BEAST.DoubleNumWiltonSauterQStrat(qs.outer_rule, qs.inner_rule, 1, 1, 1, 1, 1, 1)
    BEAST.quaddata(op, tref, bref, tels, bels, qs)
end

function BEAST.quadrule(op::typeof(SL), tref::typeof(x1), bref::typeof(x2),
    i ,τ, j, σ, qd, qs::BEAST.DoubleNumQStrat)

    return BEAST.DoubleQuadRule(
        qd.tpoints[1,i],
        qd.bpoints[1,j])
end

BEAST.quadinfo(SL,X1,X2)
BEAST.quadinfo(SL,X1,X2, quadstrat=BEAST.DoubleNumQStrat(2,3))

@time Z1 = assemble(SL,X1,X2);
@time Z2 = assemble(SL,X1,X2,quadstrat=BEAST.DoubleNumQStrat(2,3));

ba1 = blockassembler(SL,X1,X2)
ba2 = blockassembler(SL,X1,X2,quadstrat=BEAST.DoubleNumQStrat(2,3))

T = scalartype(SL,X1,X2)
n1 = numfunctions(X1)
n2 = numfunctions(X2)

W1 = zeros(T,n1,n2);
W2 = zeros(T,n1,n2);

idcs1 = collect(1:n1)
idcs2 = collect(1:n2)
@time ba1(idcs1,idcs2, (v,m,n)->(W1[m,n]+=v))
@time ba2(idcs1,idcs2, (v,m,n)->(W2[m,n]+=v));


text:
struct IPSReactionStruct
  pairwise :: Bool
  input1   :: Int # center
  input2   :: Int # adjacent
  output1  :: Int # center
  output2  :: Int # adjacent
  class    :: Int
  rate     :: Float64
  # formula    :: Expr
  # parameter  :: Symbol
  # paramidx   :: Int

  function IPSReactionStruct(pairwise, i1, i2, o1, o2, class, rate)
    if has_invalid_type(pairwise, i1, i2, o1, o2)
      throw(ErrorException("""
      Particle types must be positive integers. Non-pairwise reactions must have the adjacent types be zero.

      pairwise? $pairwise
      center:   $i1 -> $o1
      adjacent: $i2 -> $o2
      """))
    end

    if has_same_types(pairwise, i1, i2, o1, o2)
      throw(ErrorException("""
      At least one particle type must be different.

      pairwise? $pairwise
      center:   $i1 -> $o1
      adjacent: $i2 -> $o2
      """))
    end

    if has_invalid_class(class)
      throw(ErrorException("Assigned class must be a positive integer; class = $class"))
    end

    if has_invalid_rate(rate)
      throw(ErrorException("Reaction rate must be a non-negative real number; rate = $rate"))
    end

    return new(pairwise, i1, i2, o1, o2, class, rate)
  end
end

function Base.show(io::IO, x::IPSReactionStruct)
  i1, i2 = x.input1, x.input2
  o1, o2 = x.output1, x.output2

  if x.pairwise
    str = "$(i1) + $(i2) --> $(o1) + $(o2)"
  else
    str = "$(i1) --> $(o1)"
  end

  print(io, str, "  class = ", x.class)
end

function Base.show(io::IO, m::MIME"text/plain", x::IPSReactionStruct)
  show(io, x)
end

@inline function has_invalid_type(pairwise, i1, i2, o1, o2)
  if pairwise
    # all types must be positive integers
    is_invalid = (i1 < 1) || (i2 < 1) || (o1 < 1) || (o2 < 1)
  else
    # center types must be positve
    # adjacent type must be zero
    is_invalid = (i1 < 1) || (o1 < 1) || (i2 != 0) || (o2 != 0)
  end

  return is_invalid
end

@inline function has_same_types(pairwise, i1, i2, o1, o2)
  if pairwise
    is_invalid = (i1 == i2) && (i2 == o1) && (o1 == o2)
  else
    is_invalid = (i1 == o1)
  end

  return is_invalid
end

@inline function has_invalid_class(class)
  # class must be a positive integer
  class < 1
end

@inline function has_invalid_rate(rate)
  # rate must be a non-negative real number
  rate < 0
end

ispairwise(x::IPSReactionStruct) = x.pairwise

@inline get_reactant_types(x::IPSReactionStruct) = (x.input1, x.input2)
@inline get_product_types(x::IPSReactionStruct) = (x.output1, x.output2)

##### pretty printing #####

function _print_formula(x::IPSReactionStruct)
  ex = formula(x)
  return string(ex.args[1], " ", "→", " ", ex.args[2])
end

# struct InteractingParticleSystem{DG,enumType}
struct InteractingParticleSystem{enumType}
  reactions::Vector{IPSReactionStruct}
  # rxn_rates::Vector{T}
  # dep_graph::DG
  # spc_graph::DG
  # rxn_graph::DG
  isactive::Vector{Bool}
  enumeration::enumType
end

function InteractingParticleSystem(reactions::Vector{IPSReactionStruct}, isactive, enumeration)
  num_reactions = length(reactions)

  # dep_graph = rxnrxn_depgraph(DGLazy(), reactions, isactive, enumeration)

  # return InteractingParticleSystem(reactions, rxn_rates, dep_graph, spc_graph, rxn_graph)
  # return InteractingParticleSystem(reactions, dep_graph, isactive, enumeration)
  return InteractingParticleSystem(reactions, isactive, enumeration)
end

function Base.summary(io::IO, x::InteractingParticleSystem)
  types = length(x.isactive) - 1
  str = length(types) > 1 ? "types" : "type"
  print(io, "Interacting Particle System w/ $(types) ", str)
end

function Base.show(io::IO, x::InteractingParticleSystem)
  summary(io, x)
  println()
  show(io, x.reactions)
end

function Base.show(io::IO, m::MIME"text/plain", x::InteractingParticleSystem)
  summary(io, x)
  println()
  show(io, m, x.reactions)
end

##### execute_jump! #####

execute_jump!(state, model::InteractingParticleSystem, j) = execute_jump!(state, model.reactions[j], model.enumeration)

function execute_jump!(lattice, reaction, enumeration)
  # xclass = lattice.xclass
  # xdiff  = lattice.xdiff

  # save counts in each class
  # for s in eachindex(xclass)
  #   xdiff[s] = false
  #   xclass[s] = length(lattice.classes[s].members)
  # end

  if ispairwise(reaction)
    # println("pairwise")
    execute_pairwise!(lattice, reaction, enumeration)
  else
    # println("onsite")
    execute_onsite!(lattice, reaction, enumeration)
  end

  # check which counts have changed
  # for s in eachindex(xclass)
  #   if xclass[s] != length(lattice.classes[s].members)
  #     xdiff[s] = true
  #   end
  # end

  return nothing
end

function execute_pairwise!(lattice, reaction, enumeration)
  # unpack information
  i1, i2 = get_reactant_types(reaction)
  j1, j2 = get_product_types(reaction)

  s = reaction.class
  # N = lattice.N

  # sample a particle of type i1 with nbclass k; i.e. sample class s
  x = sample_from_class(lattice, enumeration, s)
  #
  # if get_ptype(x) == 1
  #   println("execute reaction:")
  #   @show reaction
  #   println("center particle:")
  #   display(x)
  #   println("class members:")
  #   for z in lattice.site[enumeration.class[s]]
  #     print("  ")
  #     display(z)
  #   end
  #   error("sampled an open site")
  # end

  # sample a particle of type i2 from x's neighborhood
  y = sample_neighbor(x, lattice, enumeration, i2)
  # we use up the particles (k1, i1) and (k2, i2)

  # change the types for x and y
  i1 != j1 && transform!(x, j1)
  i2 != j2 && transform!(y, j2)

  # update neighborhood and sample classes for x and y
  update_classes_particle!(enumeration, x, i1, j1, i2, j2)
  update_classes_particle!(enumeration, y, i2, j2, i1, j1)

  # update the remaining particles surrounding x and y
  i1 != j1 && update_classes_neighbors!(lattice, enumeration, x, y, i1, j1)
  i2 != j2 && update_classes_neighbors!(lattice, enumeration, y, x, i2, j2)

  # we produced new particles so update the counts
  # N[get_neighbor_class(x), get_ptype(x)] += 1
  # N[get_neighbor_class(y), get_ptype(y)] += 1

  return nothing
end

function execute_onsite!(lattice, reaction, enumeration)
  # unpack reaction information
  i, _ = get_reactant_types(reaction)
  j, _ = get_product_types(reaction)

  s = reaction.class
  # N = lattice.N

  # sample a particle of type i with nbclass k; i.e. sample class s
  x = sample_from_class(lattice, enumeration, s)

  # we use up the x particle
  # N[get_neighbor_class(x), get_ptype(x)] -= 1

  # change the type for x
  transform!(x, j)

  # update neighborhood and sample classes for x and y
  update_classes_particle!(enumeration, x, i, j, 0, 0)

  # update the remaining particles surrounding x and y
  update_classes_neighbors!(lattice, enumeration, x, x, i, j)

  # update the counts
  # N[get_neighbor_class(x), get_ptype(x)] += 1

  return nothing
end
function rate(model::InteractingParticleSystem, lattice, j)
  rate(model.reactions[j], lattice, model.enumeration)
end

function rate(reaction::IPSReactionStruct, lattice, enumeration)
  # grab reactant indices
  class = enumeration.class
  s = reaction.class

  # grab per particle rate
  per_particle_rate = reaction.rate

  # grab the count for the center particle
  N = length(class[s]) # need to fix this

  return N * per_particle_rate
end


text:
export Strang
"""
`Strang` matrix

A special `SymTridiagonal` matrix named after Gilbert Strang

```julia
julia> Strang(6)
6x6 Strang{Float64}:
  2.0  -1.0   0.0   0.0   0.0   0.0
 -1.0   2.0  -1.0   0.0   0.0   0.0
  0.0  -1.0   2.0  -1.0   0.0   0.0
  0.0   0.0  -1.0   2.0  -1.0   0.0
  0.0   0.0   0.0  -1.0   2.0  -1.0
  0.0   0.0   0.0   0.0  -1.0   2.0
```
"""
struct Strang{T} <: AbstractArray{T, 2}
    n :: Int
end
Strang(n::Int) = Strang{Float64}(n)
strang(T, n)= n >1 ? SymTridiagonal(2ones(T, n),-ones(T, n-1)) :
              n==1 ? Diagonal([2one(T)]) : error("Invalid dimension ", n)
function getindex(S::Strang{T}, i, j) where T
    i == j && return 2
    abs(i - j) == 1 && return -1
    0
end
getindex(S::Strang{T}, I...) where T = getindex(S,ind2sub(size(S),I...)...)
size(S::Strang, r::Int) = r==1 || r==2 ? S.n : throw(ArgumentError("Invalid dimension $r"))
size(S::Strang) = S.n, S.n
Matrix(S::Strang{T}) where T = Matrix(strang(T, S.n))
*(A::VecOrMat{T}, S::Strang{T}) where T = A*Matrix(S)
*(S::Strang{T}, A::VecOrMat{T}) where T = Matrix(S)*A


text:
# Test print_tree()

struct Num{I} end
Num(I::Int) = Num{I}()
Base.show(io::IO, ::Num{I}) where {I} = print(io, I)
AbstractTrees.children(x::Num{I}) where {I} = (Num(I+1), Num(I+1))

struct SingleChildInfiniteDepth end
AbstractTrees.children(::SingleChildInfiniteDepth) = (SingleChildInfiniteDepth(),)

@testset "Truncation" begin

    # test that `print_tree(headnode, maxdepth)` truncates the output at right depth
    # julia > print_tree(Num(0), 3)
    # 0
    # ├─ 1
    # │  ├─ 2
    # │  │  ├─ 3
    # │  │  └─ 3
    # │  └─ 2
    # │     ├─ 3
    # │     └─ 3
    # └─ 1
    #    ├─ 2
    #    │  ├─ 3
    #    │  └─ 3
    #    └─ 2
    #       ├─ 3
    #       └─ 3
    #

    for maxdepth in [3,5,8]
        ptxt = repr_tree(Num(0), maxdepth)

        n1  = sum([1 for c in ptxt if c=="$(maxdepth-1)"[1]])
        n2  = sum([1 for c in ptxt if c=="$maxdepth"[1]])
        n3  = sum([1 for c in ptxt if c=="$(maxdepth+1)"[1]])

        @test n1==2^(maxdepth-1)
        @test n2==2^maxdepth
        @test n3==0
    end

    # test that `print_tree(headnode)` prints truncation characters under each
    # node at the default maxdepth level = 5
    truncation_char = AbstractTrees.TreeCharSet().trunc
    ptxt = repr_tree(Num(0))

    n1  = sum([1 for c in ptxt if c=='5'])
    n2  = sum([1 for c in ptxt if c=='6'])
    n3  = sum([1 for c in ptxt if c==truncation_char])
    @test n1==2^5
    @test n2==0
    @test n3==2^5

    lines = split(ptxt, '\n')
    for i in 1:length(lines)
        if ~isempty(lines[i]) && lines[i][end] == '5'
            @test lines[i+1][end] == truncation_char
        end
    end

    # test correct number of lines printed 1
    ptxt = repr_tree(SingleChildInfiniteDepth())
    numlines = sum([1 for c in split(ptxt, '\n') if ~isempty(strip(c))])
    @test numlines == 7 # 1 (head node) + 5 (default depth) + 1 (truncation char)

    # test correct number of lines printed 2
    ptxt = repr_tree(SingleChildInfiniteDepth(), 3)
    numlines = sum([1 for c in split(ptxt, '\n') if ~isempty(strip(c))])
    @test numlines == 5 # 1 (head node) + 3 (depth) + 1 (truncation char)

    # test correct number of lines printed 3
    ptxt = repr_tree(SingleChildInfiniteDepth(), 3, indicate_truncation=false)
    numlines = sum([1 for c in split(ptxt, '\n') if ~isempty(strip(c))])
    @test numlines == 4 # 1 (head node) + 3 (depth)
end


text:
module MimiDICE2016

using Mimi
using XLSX: readxlsx

include("helpers.jl")
include("parameters.jl")

include("marginaldamage.jl")

include("components/totalfactorproductivity_component.jl")
include("components/grosseconomy_component.jl")
include("components/emissions_component.jl")
include("components/co2cycle_component.jl")
include("components/radiativeforcing_component.jl")
include("components/climatedynamics_component.jl")
include("components/damages_component.jl")
include("components/neteconomy_component.jl")
include("components/welfare_component.jl")

export constructdice, getdiceexcel

const model_years = 2015:5:2510

function constructdice(params)

    m = Model()
    set_dimension!(m, :time, model_years)

    #--------------------------------------------------------------------------
    # Add components in order
    #--------------------------------------------------------------------------

    add_comp!(m, totalfactorproductivity, :totalfactorproductivity)
    add_comp!(m, grosseconomy, :grosseconomy)
    add_comp!(m, emissions, :emissions)
    add_comp!(m, co2cycle, :co2cycle)
    add_comp!(m, radiativeforcing, :radiativeforcing)
    add_comp!(m, climatedynamics, :climatedynamics)
    add_comp!(m, damages, :damages)
    add_comp!(m, neteconomy, :neteconomy)
    add_comp!(m, welfare, :welfare)

    #--------------------------------------------------------------------------
    # Make internal parameter connections
    #--------------------------------------------------------------------------

    # Socioeconomics
    connect_param!(m, :grosseconomy, :AL, :totalfactorproductivity, :AL)
    connect_param!(m, :grosseconomy, :I, :neteconomy, :I)
    connect_param!(m, :emissions, :YGROSS, :grosseconomy, :YGROSS)

    # Climate
    connect_param!(m, :co2cycle, :E, :emissions, :E)
    connect_param!(m, :radiativeforcing, :MAT, :co2cycle, :MAT)
    connect_param!(m, :climatedynamics, :FORC, :radiativeforcing, :FORC)

    # Damages
    connect_param!(m, :damages, :TATM, :climatedynamics, :TATM)
    connect_param!(m, :damages, :YGROSS, :grosseconomy, :YGROSS)
    connect_param!(m, :neteconomy, :YGROSS, :grosseconomy, :YGROSS)
    connect_param!(m, :neteconomy, :DAMAGES, :damages, :DAMAGES)
	connect_param!(m, :neteconomy, :SIGMA, :emissions, :SIGMA)
    connect_param!(m, :welfare, :CPC, :neteconomy, :CPC)

    #--------------------------------------------------------------------------
    # Set external parameter values
    #--------------------------------------------------------------------------
    for (name, value) in params
        set_param!(m, name, value)
    end

    return m

end

function getdiceexcel(;datafile = joinpath(dirname(@__FILE__), "..", "data", "DICE2016R-090916ap-v2.xlsm"))
    params = getdice2016excelparameters(datafile)

    m = constructdice(params)

    return m
end

# get_model function for standard Mimi API: use the Excel version
get_model = getdiceexcel

end # module


text:
module Ipopt

using Libdl
using LinearAlgebra

if VERSION < v"1.3" || (
    haskey(ENV, "JULIA_IPOPT_LIBRARY_PATH") &&
    haskey(ENV, "JULIA_IPOPT_EXECUTABLE_PATH")
)
    const _DEPS_FILE = joinpath(dirname(@__DIR__), "deps", "deps.jl")
    if !isfile(_DEPS_FILE)
        error(
            "Ipopt not properly installed. Please run import Pkg; Pkg.build(\"Ipopt\")",
        )
    end
    include(_DEPS_FILE)
else
    import Ipopt_jll: libipopt
end

function amplexefun(::String)
    return error(
        """
  amplexefun has been removed because it caused problems with other packages.

   - On Julia 1.3 or later, use `Ipopt_jll.amplexe` instead.
   - On Julia 1.0, install Ipopt.jl version 0.6.x via `] add Ipopt@0.6`. Make
     sure to restart Julia for the changes to take effect.
  """,
    )
end

export createProblem,
    addOption,
    openOutputFile,
    setProblemScaling,
    setIntermediateCallback,
    solveProblem,
    IpoptProblem

mutable struct IpoptProblem
    ref::Ptr{Cvoid}             # Reference to the internal data structure
    n::Int                      # Num vars
    m::Int                      # Num cons
    x::Vector{Float64}          # Starting and final solution
    g::Vector{Float64}          # Final constraint values
    mult_g::Vector{Float64}     # lagrange multipliers on constraints
    mult_x_L::Vector{Float64}   # lagrange multipliers on lower bounds
    mult_x_U::Vector{Float64}   # lagrange multipliers on upper bounds
    obj_val::Float64            # Final objective
    status::Int                 # Final status

    # Callbacks
    eval_f::Function
    eval_g::Function
    eval_grad_f::Function
    eval_jac_g::Function
    eval_h::Union{Function,Nothing}
    intermediate::Union{Function,Nothing}

    sense::Symbol               # For MathProgBase

    function IpoptProblem(
        ref::Ptr{Cvoid},
        n,
        m,
        eval_f,
        eval_g,
        eval_grad_f,
        eval_jac_g,
        eval_h,
    )
        prob = new(
            ref,
            n,
            m,
            zeros(Float64, n),
            zeros(Float64, m),
            zeros(Float64, m),
            zeros(Float64, n),
            zeros(Float64, n),
            0.0,
            0,
            eval_f,
            eval_g,
            eval_grad_f,
            eval_jac_g,
            eval_h,
            nothing,
            :Min,
        )
        finalizer(freeProblem, prob)
        return prob
    end
end

# From Ipopt/src/Interfaces/IpReturnCodes_inc.h
const ApplicationReturnStatus = Dict(
    0 => :Solve_Succeeded,
    1 => :Solved_To_Acceptable_Level,
    2 => :Infeasible_Problem_Detected,
    3 => :Search_Direction_Becomes_Too_Small,
    4 => :Diverging_Iterates,
    5 => :User_Requested_Stop,
    6 => :Feasible_Point_Found,
    -1 => :Maximum_Iterations_Exceeded,
    -2 => :Restoration_Failed,
    -3 => :Error_In_Step_Computation,
    -4 => :Maximum_CpuTime_Exceeded,
    -5 => :Maximum_WallTime_Exceeded,
    -10 => :Not_Enough_Degrees_Of_Freedom,
    -11 => :Invalid_Problem_Definition,
    -12 => :Invalid_Option,
    -13 => :Invalid_Number_Detected,
    -100 => :Unrecoverable_Exception,
    -101 => :NonIpopt_Exception_Thrown,
    -102 => :Insufficient_Memory,
    -199 => :Internal_Error,
)

###########################################################################
# Callback wrappers
###########################################################################

function eval_f_wrapper(
    n::Cint,
    x_ptr::Ptr{Float64},
    new_x::Cint,
    obj_ptr::Ptr{Float64},
    user_data::Ptr{Cvoid},
)
    prob = unsafe_pointer_to_objref(user_data)::IpoptProblem
    x = unsafe_wrap(Array, x_ptr, Int(n))
    new_obj = convert(Float64, prob.eval_f(x))::Float64
    unsafe_store!(obj_ptr, new_obj)
    return Cint(1)  # Return TRUE for success.
end

function eval_g_wrapper(
    n::Cint,
    x_ptr::Ptr{Float64},
    new_x::Cint,
    m::Cint,
    g_ptr::Ptr{Float64},
    user_data::Ptr{Cvoid},
)
    prob = unsafe_pointer_to_objref(user_data)::IpoptProblem
    new_g = unsafe_wrap(Array, g_ptr, Int(m))
    x = unsafe_wrap(Array, x_ptr, Int(n))
    prob.eval_g(x, new_g)
    return Cint(1)  # Return TRUE for success.
end

function eval_grad_f_wrapper(
    n::Cint,
    x_ptr::Ptr{Float64},
    new_x::Cint,
    grad_f_ptr::Ptr{Float64},
    user_data::Ptr{Cvoid},
)
    prob = unsafe_pointer_to_objref(user_data)::IpoptProblem
    new_grad_f = unsafe_wrap(Array, grad_f_ptr, Int(n))
    x = unsafe_wrap(Array, x_ptr, Int(n))
    prob.eval_grad_f(x, new_grad_f)
    return Cint(1)  # Return TRUE for success.
end

function eval_jac_g_wrapper(
    n::Cint,
    x_ptr::Ptr{Float64},
    new_x::Cint,
    m::Cint,
    nele_jac::Cint,
    iRow::Ptr{Cint},
    jCol::Ptr{Cint},
    values_ptr::Ptr{Float64},
    user_data::Ptr{Cvoid},
)
    prob = unsafe_pointer_to_objref(user_data)::IpoptProblem
    # Determine mode
    mode = (values_ptr == C_NULL) ? (:Structure) : (:Values)
    x = unsafe_wrap(Array, x_ptr, Int(n))
    rows = unsafe_wrap(Array, iRow, Int(nele_jac))
    cols = unsafe_wrap(Array, jCol, Int(nele_jac))
    values = unsafe_wrap(Array, values_ptr, Int(nele_jac))
    prob.eval_jac_g(x, mode, rows, cols, values)
    return Cint(1)  # Return TRUE for success.
end

function eval_h_wrapper(
    n::Cint,
    x_ptr::Ptr{Float64},
    new_x::Cint,
    obj_factor::Float64,
    m::Cint,
    lambda_ptr::Ptr{Float64},
    new_lambda::Cint,
    nele_hess::Cint,
    iRow::Ptr{Cint},
    jCol::Ptr{Cint},
    values_ptr::Ptr{Float64},
    user_data::Ptr{Cvoid},
)
    prob = unsafe_pointer_to_objref(user_data)::IpoptProblem
    if prob.eval_h === nothing  # Did the user specify a Hessian?
        return Cint(0)  # Return FALSE for failure.
    end
    # Determine mode
    mode = (values_ptr == C_NULL) ? (:Structure) : (:Values)
    x = unsafe_wrap(Array, x_ptr, Int(n))
    lambda = unsafe_wrap(Array, lambda_ptr, Int(m))
    rows = unsafe_wrap(Array, iRow, Int(nele_hess))
    cols = unsafe_wrap(Array, jCol, Int(nele_hess))
    values = unsafe_wrap(Array, values_ptr, Int(nele_hess))
    prob.eval_h(x, mode, rows, cols, obj_factor, lambda, values)
    return Cint(1)  # Return TRUE for success.
end

function intermediate_wrapper(
    alg_mod::Cint,
    iter_count::Cint,
    obj_value::Float64,
    inf_pr::Float64,
    inf_du::Float64,
    mu::Float64,
    d_norm::Float64,
    regularization_size::Float64,
    alpha_du::Float64,
    alpha_pr::Float64,
    ls_trials::Cint,
    user_data::Ptr{Cvoid},
)
    prob = unsafe_pointer_to_objref(user_data)::IpoptProblem
    ret = prob.intermediate(
        Int(alg_mod),
        Int(iter_count),
        obj_value,
        inf_pr,
        inf_du,
        mu,
        d_norm,
        regularization_size,
        alpha_du,
        alpha_pr,
        Int(ls_trials),
    )
    # Return TRUE if the optimization should continue.
    return ret ? Cint(1) : Cint(0)
end

###########################################################################
# C function wrappers
###########################################################################

function createProblem(
    n::Int,
    x_L::Vector{Float64},
    x_U::Vector{Float64},
    m::Int,
    g_L::Vector{Float64},
    g_U::Vector{Float64},
    nele_jac::Int,
    nele_hess::Int,
    eval_f,
    eval_g,
    eval_grad_f,
    eval_jac_g,
    eval_h = nothing,
)
    @assert n == length(x_L) == length(x_U)
    @assert m == length(g_L) == length(g_U)
    # Wrap callbacks
    eval_f_cb = @cfunction(
        eval_f_wrapper,
        Cint,
        (Cint, Ptr{Float64}, Cint, Ptr{Float64}, Ptr{Cvoid}),
    )
    eval_g_cb = @cfunction(
        eval_g_wrapper,
        Cint,
        (Cint, Ptr{Float64}, Cint, Cint, Ptr{Float64}, Ptr{Cvoid}),
    )
    eval_grad_f_cb = @cfunction(
        eval_grad_f_wrapper,
        Cint,
        (Cint, Ptr{Float64}, Cint, Ptr{Float64}, Ptr{Cvoid}),
    )
    eval_jac_g_cb = @cfunction(
        eval_jac_g_wrapper,
        Cint,
        (
            Cint,
            Ptr{Float64},
            Cint,
            Cint,
            Cint,
            Ptr{Cint},
            Ptr{Cint},
            Ptr{Float64},
            Ptr{Cvoid},
        ),
    )
    eval_h_cb = @cfunction(
        eval_h_wrapper,
        Cint,
        (
            Cint,
            Ptr{Float64},
            Cint,
            Float64,
            Cint,
            Ptr{Float64},
            Cint,
            Cint,
            Ptr{Cint},
            Ptr{Cint},
            Ptr{Float64},
            Ptr{Cvoid},
        ),
    )
    problem_p = ccall(
        (:CreateIpoptProblem, libipopt),
        Ptr{Cvoid},
        (
            Cint,
            Ptr{Float64},
            Ptr{Float64},    # Num vars, var lower and upper bounds
            Cint,
            Ptr{Float64},
            Ptr{Float64},    # Num constraints, con lower and upper bounds
            Cint,
            Cint,                          # Num nnz in constraint Jacobian and in Hessian
            Cint,                                # 0 for C, 1 for Fortran
            Ptr{Cvoid},
            Ptr{Cvoid},              # Callbacks for eval_f, eval_g
            Ptr{Cvoid},
            Ptr{Cvoid},
            Ptr{Cvoid},  # Callbacks for eval_grad_f, eval_jac_g, eval_h
        ),
        n,
        x_L,
        x_U,
        m,
        g_L,
        g_U,
        nele_jac,
        nele_hess,
        1,
        eval_f_cb,
        eval_g_cb,
        eval_grad_f_cb,
        eval_jac_g_cb,
        eval_h_cb,
    )

    if problem_p == C_NULL
        if n == 0
            error(
                "IPOPT: Failed to construct problem because there are 0 " *
                "variables. If you intended to construct an empty problem, " *
                "one work-around is to add a variable fixed to 0.",
            )
        else
            error("IPOPT: Failed to construct problem for some unknown reason.")
        end
    end
    return IpoptProblem(
        problem_p,
        n,
        m,
        eval_f,
        eval_g,
        eval_grad_f,
        eval_jac_g,
        eval_h,
    )
end

# Note: you should have no need to call this directly. `IpoptProblem` has a
# finalizer.
function freeProblem(prob::IpoptProblem)
    if prob.ref == C_NULL
        return
    end
    ccall((:FreeIpoptProblem, libipopt), Cvoid, (Ptr{Cvoid},), prob.ref)
    prob.ref = C_NULL
    return
end

function addOption(prob::IpoptProblem, keyword::String, value::String)
    #/** Function for adding a string option.  Returns FALSE the option
    # *  could not be set (e.g., if keyword is unknown) */
    if !(isascii(keyword) && isascii(value))
        error("IPOPT: Non ASCII parameters not supported")
    end
    ret = ccall(
        (:AddIpoptStrOption, libipopt),
        Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Ptr{UInt8}),
        prob.ref,
        keyword,
        value,
    )
    if ret == 0
        error("IPOPT: Couldn't set option '$keyword' to value '$value'.")
    end
    return
end

function addOption(prob::IpoptProblem, keyword::String, value::Float64)
    #/** Function for adding a Number option.  Returns FALSE the option
    # *  could not be set (e.g., if keyword is unknown) */
    if !isascii(keyword)
        error("IPOPT: Non ASCII parameters not supported")
    end
    ret = ccall(
        (:AddIpoptNumOption, libipopt),
        Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Float64),
        prob.ref,
        keyword,
        value,
    )
    if ret == 0
        error("IPOPT: Couldn't set option '$keyword' to value '$value'.")
    end
    return
end

function addOption(prob::IpoptProblem, keyword::String, value::Integer)
    #/** Function for adding an Int option.  Returns FALSE the option
    # *  could not be set (e.g., if keyword is unknown) */
    if !isascii(keyword)
        error("IPOPT: Non ASCII parameters not supported")
    end
    ret = ccall(
        (:AddIpoptIntOption, libipopt),
        Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint),
        prob.ref,
        keyword,
        value,
    )
    if ret == 0
        error("IPOPT: Couldn't set option '$keyword' to value '$value'.")
    end
    return
end

function openOutputFile(prob::IpoptProblem, file_name::String, print_level::Int)
    #/** Function for opening an output file for a given name with given
    # *  printlevel.  Returns false, if there was a problem opening the
    # *  file. */
    if !isascii(file_name)
        error("IPOPT: Non ASCII parameters not supported")
    end
    ret = ccall(
        (:OpenIpoptOutputFile, libipopt),
        Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint),
        prob.ref,
        file_name,
        print_level,
    )
    if ret == 0
        error("IPOPT: Couldn't open output file.")
    end
    return
end

#/** Optional function for setting scaling parameter for the NLP.
# *  This corresponds to the get_scaling_parameters method in TNLP.
# *  If the pointers x_scaling or g_scaling are NULL, then no scaling
# *  for x resp. g is done. */
function setProblemScaling(
    prob::IpoptProblem,
    obj_scaling::Float64,
    x_scaling::Union{Nothing,Vector{Float64}} = nothing,
    g_scaling::Union{Nothing,Vector{Float64}} = nothing,
)
    x_scale_arg = (x_scaling === nothing) ? C_NULL : x_scaling
    g_scale_arg = (g_scaling === nothing) ? C_NULL : g_scaling
    ret = ccall(
        (:SetIpoptProblemScaling, libipopt),
        Cint,
        (Ptr{Cvoid}, Float64, Ptr{Float64}, Ptr{Float64}),
        prob.ref,
        obj_scaling,
        x_scale_arg,
        g_scale_arg,
    )
    if ret == 0
        error("IPOPT: Error setting problem scaling.")
    end
    return
end

function setIntermediateCallback(prob::IpoptProblem, intermediate::Function)
    intermediate_cb = @cfunction(
        intermediate_wrapper,
        Cint,
        (
            Cint,
            Cint,
            Float64,
            Float64,
            Float64,
            Float64,
            Float64,
            Float64,
            Float64,
            Float64,
            Cint,
            Ptr{Cvoid},
        ),
    )
    ret = ccall(
        (:SetIntermediateCallback, libipopt),
        Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}),
        prob.ref,
        intermediate_cb,
    )
    prob.intermediate = intermediate
    if ret == 0
        error("IPOPT: Something went wrong setting the intermediate callback.")
    end
    return
end

function solveProblem(prob::IpoptProblem)
    final_objval = Ref{Cdouble}(0.0)
    ret = ccall(
        (:IpoptSolve, libipopt),
        Cint,
        (
            Ptr{Cvoid},
            Ptr{Float64},
            Ptr{Float64},
            Ptr{Float64},
            Ptr{Float64},
            Ptr{Float64},
            Ptr{Float64},
            Any,
        ),
        prob.ref,
        prob.x,
        prob.g,
        final_objval,
        prob.mult_g,
        prob.mult_x_L,
        prob.mult_x_U,
        prob,
    )
    prob.obj_val = final_objval[]
    prob.status = Int(ret)
    return prob.status
end

include("MPB_wrapper.jl")
include("MOI_wrapper.jl")

end


text:
using Rimu.Hamiltonians: n_to_k, correlation_factor

struct WFunction{M}
    values::NTuple{M,Float64}
end
function WFunction(M, cutoff)
    return WFunction(Tuple(Rimu.Hamiltonians.w_function.(0:M-1, cutoff)))
end
(w::WFunction)(n) = w.values[abs(n) + 1]

struct CorrelationFactor{M}
    values::NTuple{M,Float64}
end
function CorrelationFactor(M, cutoff; length=M)
    CorrelationFactor(Tuple(correlation_factor.(1:length, cutoff, M)))
end
function (correlation_factor::CorrelationFactor{M})(n) where {M}
    absn = abs(n)
    if 0 < absn
        return sign(n) * correlation_factor.values[absn]
    else
        return 0.0
    end
end

struct TFunction{M}
    cutoff::Int
    t::Float64
    u::Float64
    correlation_factor::CorrelationFactor{M}
    w::WFunction{M}
end
function TFunction(M, cutoff, t, u)
    w = WFunction(M, cutoff)
    correlation_factor = CorrelationFactor(M, cutoff)
    return TFunction(cutoff, float(t), float(u), correlation_factor, w)
end
@fastmath function (t_fun::TFunction{M})(_, _, p, q, k) where {M}
    t, u = t_fun.t, t_fun.u
    k_pi = n_to_k(k, M)
    pmq_pi = n_to_k(p - q, M)
    cor_k = t_fun.correlation_factor(k)
    return u/M + 2u/M * (cor_k * k_pi - cor_k * pmq_pi) + 2u^2/t * t_fun.w(k)
end

struct QFunction{M}
    cutoff::Int
    t::Float64
    u::Float64
    correlation_factor::CorrelationFactor{M}
end
function QFunction(M, cutoff, t, u)
    return QFunction{M}(cutoff, float(t), float(u), CorrelationFactor(M, cutoff))
end
function (q_fun::QFunction{M})(_, _, k, l) where {M}
    t, u = q_fun.t, q_fun.u
    cor_k = q_fun.correlation_factor(k)
    cor_l = q_fun.correlation_factor(l)

    return -u^2/(t * M^2) * cor_k * cor_l
end

struct Transcorrelated{A,O} <: Hamiltonian{A,Float64}
    address::A
    t::Float64
    u::Float64
    cutoff::Int
    three_body_term::Bool
    terms::O
end

function Transcorrelated(address; t::Real=1, u::Real=1, cutoff::Int=1, three_body_term=true)
    M = num_modes(address)
    C = num_components(address)

    t_fun = TFunction(num_modes(address), cutoff, t, u)

    ts = @SVector fill(float(t), C)
    kinetic_term = ParticleCountTerm(
        KineticEnergyFunction(address, ts, continuum_dispersion)
    )
    interaction_term = MomentumTwoBodyTerm(t_fun; fold=false)
    if three_body_term
        interaction_term += MomentumThreeBodyTerm(QFunction(M, cutoff, t, u))
    end
    terms = kinetic_term + interaction_term
    return Transcorrelated(address, Float64(t), Float64(u), cutoff, three_body_term, terms)
end

function Base.show(io::IO, h::Transcorrelated)
    print(io, "Transcorrelated(")
    print(IOContext(io, :compact => true), h.address)
    print(io, "; t=$(h.t), u=$(h.u), cutoff=$(h.cutoff), three_body_term=$(h.three_body_term))")
end

terms(h::Transcorrelated) = h.terms

starting_address(h::Transcorrelated) = h.address
LOStructure(::Transcorrelated) = AdjointKnown()
function Base.adjoint(h::Transcorrelated)
    return Transcorrelated(
        h.address, h.t, h.u, h.cutoff, h.three_body_term, adjoint(h.terms)
    )
end

basis(::Transcorrelated) = MomentumSpace()


text:
using Accessors
using Accessors: IndexLens, PropertyLens, ComposedOptic

struct Lens!!{L}
    pure::L
end

function Base.show(io::IO, l::Lens!!)
    print(io, "Lens!!", l.pure)
end

(l::Lens!!)(o) = l.pure(o)

@inline function Accessors.set(o, l::Lens!!{<: ComposedOptic}, val)
    set(o, Lens!!(l.pure.outer) ∘ Lens!!(l.pure.inner), val)
end

@inline function Accessors.set(o, l::Lens!!{PropertyLens{prop}}, val) where {prop}
    if ismutable(o)
        setproperty!(o, prop, val)
    else
        setproperties(o, NamedTuple{(prop,)}((val,)))
    end
end

@inline function Accessors.set(o, l::Lens!!{typeof(identity)}, val)
    if ismutable(o)
        o .= val
    else
        val
    end
end

@inline function Accessors.set(o, l::Lens!!{<:IndexLens}, val)
    if ismutable(o)
        setindex!(o, val, l.pure.indices...)
    else
        Base.setindex(o, val, l.pure.indices...)
    end
end

@inline function Accessors.modify(f, o, l::Lens!!)
    set(o, l, f(l(o)))
end

@inline function Accessors.modify(f, o, l::Lens!!{<:ComposedOptic})
    o_inner = l.pure.inner(o)
    modify(f, o_inner, Lens!!(l.pure.outer))
end

using Accessors: setmacro, opticmacro, modifymacro

macro set!!(ex)
    setmacro(Lens!!, ex; overwrite=true)
end

macro optic!!(ex)
    opticmacro(Lens!!, ex)
end

macro modify!!(f, ex)
    modifymacro(Lens!!, f, ex)
end

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

function unescape(ast)
    leaf(x) = x
    @inline function branch(f, head, args)
        default() = Expr(head, map(f, args)...)

        head == :escape ? f(args[1]) : default()
    end
    foldast(leaf, branch)(ast)
end

function opticize(ast)
    leaf(x) = x
    @inline function branch(f, head, args)
        default() = Expr(head, map(f, args)...)
        head == :call || return default()
        first(args) == :~ || return default()
        length(args) == 3 || return default()

        # If we get here, we know we're working with something like `lhs ~ rhs`
        lhs = args[2]
        rhs = args[3]

        lhs′ = @match lhs begin
            :(($(x::Symbol), $o)) => :(($x, $o))
            :(($(x::Var), $o)) => :(($x, $o))
            _ => begin
                (x, o) = unescape.(Accessors.parse_obj_optic(lhs))
                :(($x, $o))
            end
        end
        rhs′ = f(rhs)
        :($lhs′ ~ $rhs′)
    end

    foldast(leaf, branch)(ast)
end