text:
using Distributed
if haskey(ENV, "BENCHMARK_PROCS")
    const np, nt = parse.(Ref(Int), split(ENV["BENCHMARK_PROCS"], ":"))
    for i in workers()
        addprocs(np; exeflags="-t $nt")
    end
else
    const np = 2
    const nt = 1
end
if haskey(ENV, "BENCHMARK_REMOTES")
    const remotes = split(ENV["BENCHMARK_REMOTES"], ":")
    if !isempty(remotes)
        for i in 1:np
            addprocs(remotes; exeflags="-t $nt")
        end
    end
end
@everywhere using Dagger
import Dagger: Computation, reduceblock
using Dates, Random, Statistics, LinearAlgebra, InteractiveUtils

const output_format = get(ENV, "BENCHMARK_OUTPUT_FORMAT", "jls")
if output_format == "jld"
    using JLD
elseif output_format == "jls"
    using Serialization
else
    error("Unknown output format: $output_format")
end

const render = get(ENV, "BENCHMARK_RENDER", "")
if render == "live"
    const live = true
    using Luxor, ProfileSVG
    using Mux
elseif render == "offline"
    const live = false
    using Luxor, ProfileSVG
    using FFMPEG, FileIO, ImageMagick
end
const RENDERS = Dict{Int,Dict}()
const live_port = parse(Int, get(ENV, "BENCHMARK_LIVE_PORT", "8000"))

const graph = parse(Bool, get(ENV, "BENCHMARK_GRAPH", "0"))
const profile = parse(Bool, get(ENV, "BENCHMARK_PROFILE", "0"))

_benches = get(ENV, "BENCHMARK", "cpu,cpu+dagger")
const benches = []
for bench in split(_benches, ',')
    if endswith(bench, "+dagger")
        accel = split(bench, '+')[1]
        dagger = true
    else
        accel = bench
        dagger = false
    end
    push!(benches, (name=bench, accel=accel, dagger=dagger))
end

if any(x->x.accel == "cuda", benches)
    @everywhere using DaggerGPU, CUDA
elseif any(x->x.accel == "amdgpu", benches)
    @everywhere using DaggerGPU, AMDGPU
end

const scales = eval(Meta.parse(get(ENV, "BENCHMARK_SCALE", "1:5:50")))

using BenchmarkTools

AmulB_compatible(x,y) =
    try; size(x), size(y), x * y; true; catch err false end

AmulB_compatible(x::Computation,y::Computation) = AmulB_compatible(domain(x), domain(y))

function AcmulB_compatible(x,y)
    AmulB_compatible(domain(x)', domain(y))
end

function array_suite_inner(ctx, X, Y, f=x->x)
    suite = BenchmarkGroup(["array"])

    suite["alloc"] = @benchmarkable (compute($f($X)); GC.gc())

    Y = compute(ctx, f(X))
    suite["X.+1.0"]  = @benchmarkable (compute($ctx, $f($Y.+1)); GC.gc())
    suite["X+X"]     = @benchmarkable (compute($ctx, $f($Y+$Y)); GC.gc())
    suite["X'"]      = @benchmarkable (compute($ctx, $f($Y')); GC.gc())
    suite["sin.(X)"] = @benchmarkable (compute($ctx, $f(sin.($Y))); GC.gc())

    if ndims(Y) == 2
        suite["X+X'"] = @benchmarkable (compute($ctx, $f($Y+$Y')); GC.gc())
        suite["X*X"]  = @benchmarkable (compute($ctx, $f($Y*$Y)); GC.gc())
        suite["X'*X"] = @benchmarkable (compute($ctx, $f($Y'*$Y)); GC.gc())
    end

    suite
end

function nnmf(X, W, H)
    # H update
    H = (H .* (W' * (X ./ (W * H)))
         ./ (sum(W; dims=1))')
    # W update
    W = (W .* ((X ./ (W * H)) * (H'))
         ./ (sum(H; dims=2)'))
    # error estimate
    X - W * H
end

theory_flops(nrow, ncol, nfeatures) = 11 * ncol * nrow * nfeatures + 2 * (ncol + nrow) * nfeatures

function nmf_suite(ctx; dagger, accel)
    suite = BenchmarkGroup()

    #= TODO: Re-enable
    pairs = Set{Tuple}()
    for i in (1, 2, 5, 10, 64, 100, 1028)
        for bscale in (1, 0.5, 0.3, 0.25, 0.1)
            b = ceil(Int, i*bscale)
            if !((b,i) in pairs)
                push!(pairs, (b,i))
                X = zeros(Blocks(b), i)
                Y = compute(ctx, f(X))
                suite["vector $i/$b"] = array_suite_inner(ctx, X, Y, f)
            end
            if !(((b,b),(i,i)) in pairs)
                push!(pairs, ((b,b),(i,i)))
                X = zeros(Blocks(b,b), i, i)
                Y = compute(ctx, f(X))
                suite["matrix ($i,$i)/$b"] = array_suite_inner(ctx, X, Y, f)
            end
        end
    end
    =#

    X = Ref{Any}()
    W = Ref{Any}()
    H = Ref{Any}()

    for scale in scales
        ncol = 2001 * scale
        nrow = 1002
        nfeatures = 12

        if !dagger
            suite["NNMF scaled by: $scale"] = @benchmarkable begin
                nnmf($X[], $W[], $H[])
            end setup=begin
                _scale = $scale
                @info "Starting non-Dagger NNMF (scale by $_scale)"
                if $accel == "cuda"
                    $X[] = CUDA.rand(Float32, $nrow, $ncol)
                    $W[] = CUDA.rand(Float32, $nrow, $nfeatures)
                    $H[] = CUDA.rand(Float32, $nfeatures, $ncol)
                elseif $accel == "amdgpu"
                    $X[] = ROCArray(rand(Float32, $nrow, $ncol))
                    $W[] = ROCArray(rand(Float32, $nrow, $nfeatures))
                    $H[] = ROCArray(rand(Float32, $nfeatures, $ncol))
                elseif $accel == "cpu"
                    $X[] = rand(Float32, $nrow, $ncol)
                    $W[] = rand(Float32, $nrow, $nfeatures)
                    $H[] = rand(Float32, $nfeatures, $ncol)
                end
            end teardown=begin
                $X[] = nothing
                $W[] = nothing
                $H[] = nothing
                @everywhere GC.gc()
            end
        else
            RENDERS[scale] = Dict{Int,Vector}()
            nw = length(workers())
            nsuite = BenchmarkGroup()
            while nw > 0
                opts = if accel == "cuda"
                    Dagger.Sch.SchedulerOptions(;proctypes=[
                        DaggerGPU.CuArrayDeviceProc
                    ])
                elseif accel == "amdgpu"
                    Dagger.Sch.SchedulerOptions(;proctypes=[
                        DaggerGPU.ROCArrayProc
                    ])
                elseif accel == "cpu"
                    Dagger.Sch.SchedulerOptions()
                else
                    error("Unknown accelerator $accel")
                end
                p = sum([length(Dagger.get_processors(OSProc(id))) for id in 2:(nw+1)])
                #bsz = ncol ÷ length(workers())
                bsz = ncol ÷ 64
                nsuite["Workers: $nw"] = @benchmarkable begin
                    _ctx = Context($ctx, workers()[1:$nw])
                    compute(_ctx, nnmf($X[], $W[], $H[]); options=$opts)
                end setup=begin
                    _nw, _scale = $nw, $scale
                    @info "Starting $_nw worker Dagger NNMF (scale by $_scale)"
                    if $accel == "cuda"
                        # FIXME: Allocate with CUDA.rand if possible
                        $X[] = Dagger.mapchunks(CUDA.cu, compute(rand(Blocks($nrow, $bsz), Float32, $nrow, $ncol); options=$opts))
                        $W[] = Dagger.mapchunks(CUDA.cu, compute(rand(Blocks($nrow, $bsz), Float32, $nrow, $nfeatures); options=$opts))
                        $H[] = Dagger.mapchunks(CUDA.cu, compute(rand(Blocks($nrow, $bsz), Float32, $nfeatures, $ncol); options=$opts))
                    elseif $accel == "amdgpu"
                        $X[] = Dagger.mapchunks(ROCArray, compute(rand(Blocks($nrow, $ncol÷$p), Float32, $nrow, $ncol); options=$opts))
                        $W[] = Dagger.mapchunks(ROCArray, compute(rand(Blocks($nrow, $ncol÷$p), Float32, $nrow, $nfeatures); options=$opts))
                        $H[] = Dagger.mapchunks(ROCArray, compute(rand(Blocks($nrow, $ncol÷$p), Float32, $nfeatures, $ncol); options=$opts))
                    elseif $accel == "cpu"
                        $X[] = compute(rand(Blocks($nrow, $bsz), Float32, $nrow, $ncol); options=$opts)
                        $W[] = compute(rand(Blocks($nrow, $bsz), Float32, $nrow, $nfeatures); options=$opts)
                        $H[] = compute(rand(Blocks($nrow, $bsz), Float32, $nfeatures, $ncol); options=$opts)
                    end
                end teardown=begin
                    if render != "" && !live
                        Dagger.continue_rendering[] = false
                        for i in 1:5
                            isready(Dagger.render_results) && break
                            sleep(1)
                        end
                        if isready(Dagger.render_results)
                            video_paths = take!(Dagger.render_results)
                            try
                                video_data = Dict(key=>read(video_paths[key]) for key in keys(video_paths))
                                push!(get!(()->[], RENDERS[$scale], $nw), video_data)
                            catch err
                                @error "Failed to process render results" exception=(err,catch_backtrace())
                            end
                        else
                            @warn "Failed to fetch render results"
                        end
                    end
                    $X[] = nothing
                    $W[] = nothing
                    $H[] = nothing
                    @everywhere GC.gc()
                end
                nw ÷= 2
            end
            suite["NNMF scaled by: $scale"] = nsuite
        end
    end

    suite
end

function main()
    nw = length(workers())
    output_prefix = "result-$(np)workers-$(nt)threads-$(Dates.now())"

    suites = Dict()
    graph_opts = if graph && render != ""
        (log_sink=Dagger.LocalEventLog(), log_file=output_prefix*".dot")
    elseif render != ""
        (log_sink=Dagger.LocalEventLog(),)
    else
        NamedTuple()
    end
    ctx = Context(collect((1:nw) .+ 1); profile=profile, graph_opts...)
    for bench in benches
        name = bench.name
        println("creating $name benchmarks")
        suites[name] = nmf_suite(ctx; dagger=bench.dagger, accel=bench.accel)
    end
    if render != ""
        Dagger.show_gantt(ctx; width=1800, window_length=5, delay=2, port=live_port, live=live)
        if live
            # Make sure server code is compiled
            sleep(1)
            run(pipeline(`curl -s localhost:$live_port/`; stdout=devnull))
            run(pipeline(`curl -s localhost:$live_port/profile`; stdout=devnull))
            @info "Rendering started on port $live_port"
        end
    end
    res = Dict()
    for bench in benches
        name = bench.name
        println("running $name benchmarks")
        res[name] = try
            run(suites[name]; samples=3, seconds=10*60, gcsample=true)
        catch err
            @error "Error running $name benchmarks" exception=(err,catch_backtrace())
            nothing
        end
    end
    for bench in benches
        println("benchmark results for $(bench.name): $(minimum(res[bench.name]))")
    end

    println("saving results in $output_prefix.$output_format")
    if output_format == "jld"
        JLD.save(output, "results", res, "peakflops", peakflops(), "renders", RENDERS)
    elseif output_format == "jls"
        outdict = Dict("results"=>res, "peakflops"=>peakflops(), "renders"=>RENDERS)
        open(output_prefix*".jls", "w") do io
            serialize(io, outdict)
        end
    end

    if parse(Bool, get(ENV, "BENCHMARK_VISUALIZE", "0"))
        run(`$(Base.julia_cmd()) $(joinpath(pwd(), "visualize.jl")) -- $(output_prefix*"."*output_format)`)
    end

    println("Done.")

    # TODO: Compare with multiple results
    if length(ARGS) == 1
        compare_file = ARGS[1]
        println("Comparing with results in $compare_file")
        reference_results = JLD.load(compare_file, "result")
        @show judge(mean(reference_results), mean(res))
    end
end

main()


text:
using Documenter
using MPI
using MPIPreferences

# generate example markdown
EXAMPLES = [
    "Hello world" => "examples/01-hello.md",
    "Broadcast" => "examples/02-broadcast.md",
    "Reduce" => "examples/03-reduce.md",
    "Send/receive" => "examples/04-sendrecv.md",
    "Job Scheduling" => "examples/05-job_schedule.md",
    "Scatterv and Gatherv" => "examples/06-scatterv.md",
    "Active RMA" => "examples/07-rma_active.md",
    "Passive RMA" => "examples/08-rma_passive.md",
    "Graph Communication" => "examples/09-graph_communication.md",
]

examples_md_dir = joinpath(@__DIR__,"src/examples")
isdir(examples_md_dir) || mkdir(examples_md_dir)

for (example_title, example_md) in EXAMPLES
    example_jl = example_md[1:end-2]*"jl"
    @info "Building $example_md"
    open(joinpath(@__DIR__, "src", example_md), "w") do mdfile
        println(mdfile, """
            ```@meta
            EditURL = "https://github.com/JuliaParallel/MPI.jl/blob/master/docs/$(example_jl)"
            ```
            """
        )
        println(mdfile, "# $example_title")
        println(mdfile)
        println(mdfile, "```julia")
        println(mdfile, "# $example_jl")
        println(mdfile, readchomp(joinpath(@__DIR__,example_jl)))
        println(mdfile, "```")
        println(mdfile)

        println(mdfile, "```")
        println(mdfile, "> mpiexecjl -n 4 julia $example_jl")
        cd(@__DIR__) do
            write(mdfile, mpiexec(cmd -> read(`$cmd -n 4 $(Base.julia_cmd()) --project $example_jl`)))
        end
        println(mdfile, "```")
    end
end

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

makedocs(
    sitename = "MPI.jl",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    modules = [MPI, MPIPreferences],
    pages = Any[
        "index.md",
        "configuration.md",
        "usage.md",
        "knownissues.md",
        "Examples" => EXAMPLES,
        "Reference" => [
            "reference/mpipreferences.md",
            "reference/library.md",
            "reference/environment.md",
            "reference/comm.md",
            "reference/buffers.md",
            "reference/pointtopoint.md",
            "reference/collective.md",
            "reference/onesided.md",
            "reference/topology.md",
            "reference/io.md",
            "reference/advanced.md",
        ],
        "refindex.md",
    ],
    strict = true,
)

deploydocs(
    repo = "github.com/JuliaParallel/MPI.jl.git",
    push_preview = true,
)

ext: jl


text:
module SLMC_Flux
    using Flux
end


text:
#Equilibrium Moisture Content
#The National Fire Danger Rating Sysytem: Basic Equations
#Jack D Cohen, John E. Deeming
#GTR PSW-82

function emc(relative_humidity::Float64,temp::Float64)
    if relative_humidity < 10.0
        emc = 0.03229 + 0.281073 * relative_humidity - 0.000578 * temp * relative_humidity
    elseif relative_humidity >= 10.0 && relative_humidity < 50.0
        emc = 2.22749+0.160107 * relative_humidity - 0.014784 * temp
    else
        emc = 21.0606 + 0.005565 * relative_humidity ^ 2 - 0.00035 * relative_humidity * temp - 0.483199 * relative_humidity
    end
end


text:
# To test a simplified quadrotor.
# Modeling the quadrotor frame as a cube and the props
# as another rigid body with a revolute joint between the two.

include("../src/plotSol.jl")
include("../src/simulateDiff.jl")
include("../src/OrientationConversion.jl")
##
clearconsole()
j1 = Nothing
j2 = Nothing
##
mb = 0.155; mp = 0.02;

Ib =  [ 0.0009    0.0000   -0.0000;
        0.0000    0.0053    0.0000;
       -0.0000    0.0000    0.0054]#*(10^-3)

Ip =    [ 0.0029    0.0000    0.0000;
          0.0000    0.5824    0.0000;
          0.0000    0.0000    0.5805]*(10^-4)

InFrame = InertialFrameAsRB()
QuadCube = RigidBody(mb,Ib,2)
Props = RigidBody(mp,Ip,3)

x0Cube = [zeros(3);[1;zeros(3)];zeros(3);zeros(4)]
initialiseRigidBody!(QuadCube,x0Cube)
x0Props = [[0.0;0.0;0.5];[1;zeros(3)];zeros(3);zeros(4)]
initialiseRigidBody!(Props,x0Props)

axis = [0.0 0.0 1.0][:]
rjCube = [0.0 0.0 0.5][:]
rjProp = [0.0 0.0 0.0][:]

j1 = Joint(InFrame,QuadCube,zeros(3),zeros(3))
j2 = Joint(QuadCube,Props,rjCube,rjProp,type = "Revolute",axis = axis);#,jointTorque = [0.0;0.0;0.0001])

j = [j1,j2];

## Linearization
# x0, u0 = getXU_0(j);
# A,B = linearizeDiff(x0, u0, j);
# println("A = ",A)
# println("B = ",B)
# sleep(1200)
## Simulation
tEnd = 0.01
tSpan = 0.01
g = [0.0,0.0,0.0]; # Gravity Vector.
@btime simulateDiff(tEnd, tInt, j, g=g);
sleep(1000);
tSim, solFinal = simulate(tEnd,tSpan,j,g=g)

solQuad = solFinal[1]
solProp = solFinal[2]
##
ωCube, ωProp, ωRel, jointLoc = checkRevJoint(solQuad, solProp, rjCube, rjProp);

##
plotErrNorm(tSim,solQuad.β)
plotErrNorm(tSim,solProp.β)
# Check if joint location has moved
# jointLoc = Matrix{Float64}(undef,length(tSim),3)
# ωCube = Matrix{Float64}(undef,length(tSim),3)
# ωProp = Matrix{Float64}(undef,length(tSim),3)
# ωPropInCube = Matrix{Float64}(undef,length(tSim),3)
# ωRel = Matrix{Float64}(undef,length(tSim),3)
# for i=1:length(tSim)
#     QuadCube.dcm = quat2dcm(solQuad.β[i,:])
#     ωCube[i,:] = angVel(solQuad.β[i,:],solQuad.βdot[i,:])
#     ωProp[i,:] = angVel(solProp.β[i,:],solProp.βdot[i,:])
#     ωPropInCube[i,:] = quat2dcm(solQuad.β[i,:])*transpose(quat2dcm(solProp.β[i,:]))*ωProp[i,:]
#     ωRel[i,:] = ωPropInCube[i,:] - ωCube[i,:]
#     jointLoc[i,:] = solQuad.r[i,:] + transpose(QuadCube.dcm)*rjCube - solProp.r[i,:]
# end
##
close("all");
plotPos(tSim,jointLoc)
plotPos(tSim, solQuad.r)
# plotVel(tSim,solQuad.v - solProp.v)
# plotAngVel(tSim,ωCube)
plotAngVel(tSim,ωProp)
# plotAngVel(tSim, ωRel)
# plotPos(tSim, solQuad.r)
# plotPos(tSim, solProp.r)
# plotQuat(tSim, solProp.β)
# plotPos(tSim, solQuad.r + solProp.r)
##

# plotAngVel(tSim,ωRel)
# plotAngVel(tSim,ωCube - ωProp)

# # Check revJoint axis
# axisSol = Matrix{Float64}(undef,length(tSim),3)
# errAxis = Matrix{Float64}(undef,length(tSim),3)
# errAxisNorm = Vector{Float64}(undef,length(tSim))
# for i=1:length(tSim)
#     axisSol[i,:] = solProp.β[i,2:4]/norm(solProp.β[i,2:4],2)
#
#     # β1inv = [solQuad.β[i,1];solQuad.β[i,2:4]]
#     # βRev = quaternionProduct(solProp.β[i,:],β1inv)
#     # axisSol[i,:] = βRev[2:4]/norm(βRev[2:4])
#
#     # axisSol[i,:] = (solQuad.β[i,:])*solProp.β[i,2:4]/norm(solProp.β[i,2:4],2)
#
#     # axisSol[i,:] = axixRev(solQuad.β[i,:],solProp.β[i,:])
#     errAxis[i,:] = axisSol[i,:] - axis
#     errAxisNorm[i] = norm(axisSol[i,:],2) - 1.0
# end
# plot(tSim[2:end],errAxisNorm[2:end])
# plotPos(tSim,errAxis)

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
using OrdinaryDiffEq, BenchmarkTools, StaticArrays
abstract type ContinuousDynamicalSystem end

function lorenz1(u0=[0.0, 10.0, 0.0]; σ = 10.0, ρ = 28.0, β = 8//3)
    @inline @inbounds function eom_lorenz!(du, u)
        du[1] = σ*(u[2]-u[1])
        du[2] = u[1]*(ρ-u[3]) - u[2]
        du[3] = u[1]*u[2] - β*u[3]
    end
    @inline @inbounds function jacob_lorenz(u)
        i = one(eltype(u))
        o = zero(eltype(u))
        @SMatrix [-σ*i           σ*i    zero(i);
                  (ρ*i - u[3])   (-i)   (-u[1]);
                  u[2]           u[1]   (-β*i) ]
    end# should give exponents [0.9056, 0, -14.5723]
    return ContinuousDS1(u0, eom_lorenz!, jacob_lorenz)
end
mutable struct ContinuousDS1{T<:AbstractVector, F, J} <: ContinuousDynamicalSystem
    state::T
    eom!::F
    jacob::J
end

function lorenz2(u0=[0.0, 10.0, 0.0]; σ = 10.0, ρ = 28.0, β = 8//3)
    @inline @inbounds function eom_lorenz!(du, u)
        du[1] = σ*(u[2]-u[1])
        du[2] = u[1]*(ρ-u[3]) - u[2]
        du[3] = u[1]*u[2] - β*u[3]
    end
    i = one(eltype(u0))
    o = zero(eltype(u0))
    J = zeros(eltype(u0), 3, 3)
    J[1,:] .= [-σ*i,        σ*i,    o]
    J[2,:] .= [ρ*i - u0[3],   -i,   -u0[1]]
    J[3,:] .= [u0[2],        u0[1],   (-β*i)]

    @inline @inbounds function jacob_lorenz!(J, u)
        J[2,1] = ρ - u[3]
        J[2,3] = -u[1]
        J[3,1] = u[2]; J[3,2] = u[1]
    end
    return ContinuousDS2(u0, eom_lorenz!, jacob_lorenz!, J)
end
mutable struct ContinuousDS2{T<:Real, F, J} <: ContinuousDynamicalSystem
    state::Vector{T}
    eom!::F
    jacob::J
    J::Matrix{T}
end

### Test jacobian function speed
function bjac()
    ds1 = lorenz1()
    b1 = @benchmark $ds1.jacob($ds1.state)
    println("Jacobian of 1: ")
    Juno.render(ds1.jacob(ds1.state))
    Juno.render(b1)

    ds2 = lorenz2()
    b2 = @benchmark $ds2.jacob($ds2.J, $ds2.state)
    println("Jacobian of 2: ")
    Juno.render(ds2.J)
    Juno.render(b2)
end
bjac()

### TANGENT BUNDLE!!!
function tangentbundle_setup_integrator(ds::ContinuousDS1, t_final;
  diff_eq_kwargs=Dict())

    D = dimension(ds)
    f! = ds.eom!
    jac = ds.jacob

    tbeom! = (t, u, du) -> begin
        f!(view(du, :, 1), u)
        A_mul_B!(
            view(du, :, 2:D+1),
            jac(view(u, :, 1)),
            view(u, :, 2:D+1)
        )
    end

    S = [ds.state eye(eltype(ds.state), D)]
    tbprob = ODEProblem(tbeom!, S, (zero(t_final), t_final))
    tb_integ = init(tbprob, Tsit5(); diff_eq_kwargs..., save_everystep=false)
    return tb_integ
end

function tangentbundle_setup_integrator(ds::ContinuousDS2, t_final;
  diff_eq_kwargs=Dict())

    D = dimension(ds)
    f! = ds.eom!

    tbeom! = (t, u, du) -> begin
        us = view(u, :, 1)
        f!(view(du, :, 1), us)
        ds.jacob(ds.J, us)
        A_mul_B!(
            view(du, :, 2:D+1),
            ds.J,
            view(u, :, 2:D+1)
        )
    end

    S = [ds.state eye(eltype(ds.state), D)]
    tbprob = ODEProblem(tbeom!, S, (zero(t_final), t_final))
    tb_integ = init(tbprob, Tsit5(); diff_eq_kwargs..., save_everystep=false)
    return tb_integ
end

# Benchmark it

function lyapunovspectrum1(ds::ContinuousDynamicalSystem, N::Real=1000;
    Ttr::Real = 0.0, diff_eq_kwargs::Dict = Dict(), dt::Real = 0.1)
    # Initialize
    tstops = dt:dt:N*dt
    D = dimension(ds)
    λ = zeros(eltype(ds.state), D)
    Q = eye(eltype(ds.state), D)
    # Transient evolution:
    Ttr != 0 && evolve!(ds, Ttr; diff_eq_kwargs = diff_eq_kwargs)
    # Create integrator for dynamics and tangent space:
    integ = tangentbundle_setup_integrator(
    ds, tstops[end]; diff_eq_kwargs = diff_eq_kwargs)

    # Main algorithm
    for τ in tstops
        integ.u[:, 2:end] .= Q # update tangent dynamics state (super important!)
        u_modified!(integ, true)
        # Integrate
        while integ.t < τ
            step!(integ)
        end
        # Perform QR (on the tangent flow):
        Q, R = qr(view(integ.u, :, 2:D+1))
        # Add correct (positive) numbers to Lyapunov spectrum
        for j in 1:D
            λ[j] += log(abs(R[j,j]))
        end
    end
    λ./(integ.t) #return spectrum
end

function lyapunovspectrum1(ds::ContinuousDynamicalSystem, N::Real=1000;
    Ttr::Real = 0.0, diff_eq_kwargs::Dict = Dict(), dt::Real = 0.1)
    # Initialize
    tstops = dt:dt:N*dt
    D = dimension(ds)
    λ = zeros(eltype(ds.state), D)
    Q = eye(eltype(ds.state), D)
    # Transient evolution:
    Ttr != 0 && evolve!(ds, Ttr; diff_eq_kwargs = diff_eq_kwargs)
    # Create integrator for dynamics and tangent space:
    integ = tangentbundle_setup_integrator(
    ds, tstops[end]; diff_eq_kwargs = diff_eq_kwargs)

    # Main algorithm
    for τ in tstops
        integ.u[:, 2:end] .= Q # update tangent dynamics state (super important!)
        u_modified!(integ, true)
        # Integrate
        while integ.t < τ
            step!(integ)
        end
        # Perform QR (on the tangent flow):
        Q, R = qr(view(integ.u, :, 2:D+1))
        # Add correct (positive) numbers to Lyapunov spectrum
        for j in 1:D
            λ[j] += log(abs(R[j,j]))
        end
    end
    λ./(integ.t) #return spectrum
end
dimension(ds) = 3

function binteg()
    A = Dict(:abstol=> 1e-6, :reltol=>1e-6)
    ds1 = lorenz1()
    ls1 = lyapunovspectrum1(ds1, 5000.0, diff_eq_kwargs=A)
    b1 = @benchmark lyapunovspectrum1($ds1, 500.0, diff_eq_kwargs=$A)
    println("Lyapunovs 1: $(ls1)")
    Juno.render(b1)

    ds2 = lorenz2()
    ls2 = lyapunovspectrum1(ds2, 5000.0, diff_eq_kwargs=A)
    b2 = @benchmark lyapunovspectrum1($ds2, 500.0, diff_eq_kwargs=$A)
    println("Lyapunovs 2: $(ls2)")

    Juno.render(b2)

end
binteg()

#
# b1 = @benchmark ds1.jacob(ds1.state)
# println("Jacobian of 1: ")
# Juno.render(b1)
#
#
# ds2 = lorenz2()
# b2 = @benchmark ds2.jacob(ds2.J, ds2.state)
# println("Jacobian of 2: ")
# Juno.render(b2)


text:
function test_dynamic_table()
    # Section 2.3.2: The dynamic table can contain duplicate entries (i.e., entries
    # with the same name and same value). Therefore, duplicate entries MUST NOT be
    # treated as an error by a decoder.
    table = HPack.DynamicTable()
    HPack.add_header!(table, (":method", "GET"))
    HPack.add_header!(table, (":method", "GET"))

    # Section 2.3.3: Indices strictly greater than the sum of the lengths of both
    # tables MUST be treated as a decoding error.

    table = HPack.DynamicTable()
    HPack.set_max_table_size!(table, 1)
    @test_throws DecodeError HPack.get_header(table, 100)
end

test_dynamic_table()


text:
module IdentityRanges

using OffsetArrays
using Base.Broadcast: DefaultArrayStyle

export IdentityRange

"""
    IdentityRange(start, stop) -> idr
    IdentityRange(r::AbstractUnitRange) -> idr

Defines an `AbstractUnitRange` where `idr[i] == i` for any valid `i`,
or equivalently `axes(idr, 1)` returns a range with the same values
as present in `idr`.

These are particularly useful for creating `view`s of arrays that
preserve the supplied axes:
```jldoctest
julia> a = rand(8);

julia> v1 = view(a, 3:5);

julia> axes(v1, 1)
Base.OneTo(3)

julia> idr = IdentityRange(3:5)
IdentityRange(3:5)

julia> v2 = view(a, idr);

julia> axes(v2, 1)
3:5
```
"""
struct IdentityRange{T<:Integer} <: AbstractUnitRange{T}
    start::T
    stop::T
    IdentityRange{T}(start, stop) where {T} = new{T}(start, stop)
end
IdentityRange(start::T, stop::T) where {T<:Integer} = IdentityRange{T}(start, stop)

Base.axes(r::IdentityRange) = (r.start:r.stop,)
Base.unsafe_indices(r::IdentityRange) = (r.start:r.stop,)

_length(r::IdentityRange{T}) where {T} = max(zero(T), convert(T, r.stop-r.start+1))
Base.length(r::IdentityRange) = _length(r)
Base.length(r::IdentityRange{T}) where {T<:Union{Int,Int64}} = _length(r)

let smallint = (Int === Int64 ?
                Union{Int8,UInt8,Int16,UInt16,Int32,UInt32} :
                Union{Int8,UInt8,Int16,UInt16})
    Base.length(r::IdentityRange{T}) where {T<:smallint} = Int(_length(r))
end

Base.first(r::IdentityRange) = r.start
Base.last(r::IdentityRange) = r.stop

Base.eltype(r::IdentityRange{T}) where {T} = typeof(one(T) + one(T))
function Base.iterate(r::IdentityRange{T}) where {T}
    last(r) < first(r) && return nothing
    x = oftype(one(T) + one(T), first(r))
    return (x,x)
end
function Base.iterate(r::IdentityRange{T}, i) where T
    x = i + one(T)
    return (x > oftype(x,last(r)) ? nothing : (x,x))
end

@inline function Base.getindex(v::IdentityRange{T}, i::Integer) where T
    @boundscheck ((i >= first(v)) & (i <= last(v))) || Base.throw_boundserror(v, i)
    convert(T, i)
end

@inline function Base.getindex(r::IdentityRange{R}, s::AbstractUnitRange{S}) where {R,S<:Integer}
    @boundscheck checkbounds(r, s)
    IdentityRange{R}(first(s), last(s))
end
Base.firstindex(r::IdentityRange) = r.start
Base.lastindex(r::IdentityRange) = r.stop

Base.intersect(r::IdentityRange, s::IdentityRange) = IdentityRange(max(first(r), first(s)),
                                                                   min(last(r), last(s)))

Base.:(==)(r::IdentityRange, s::IdentityRange) = (first(r) == first(s)) & (step(r) == step(s)) & (last(r) == last(s))
Base.:(==)(r::IdentityRange, s::OrdinalRange) = (first(r) == first(s) == first(axes(s)[1])) & (step(r) == step(s)) & (last(r) == last(s))
Base.:(==)(s::OrdinalRange, r::IdentityRange) = r == s

function Base.:+(r::IdentityRange, s::IdentityRange)
    indsr = axes(r, 1)
    indsr == axes(s, 1) || throw(DimensionMismatch("axes $indsr and $(axes(s, 1)) do not match"))
    OffsetArray(convert(UnitRange, r)+convert(UnitRange, s), indsr)
end
function Base.:-(r::IdentityRange, s::IdentityRange)
    indsr = axes(r, 1)
    indsr == axes(s, 1) || throw(DimensionMismatch("axes $indsr and $(axes(s, 1)) do not match"))
    OffsetArray(fill(first(r)-first(s), length(r)), indsr)
end

# Unary operation

function Base.:-(r::IdentityRange)
    indsr = axes(r, 1)
    OffsetArray(-indsr, indsr)
end

# Binary operations

for T in (Real, Number)
    @eval begin
        function Broadcast.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::IdentityRange, x::$T)
            indsr = axes(r, 1)
            OffsetArray(indsr.+x, indsr)
        end
        Broadcast.broadcasted(::DefaultArrayStyle{1}, ::typeof(+),  x::$T, r::IdentityRange) = r .+ x
        function Broadcast.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::IdentityRange, x::$T)
            indsr = axes(r, 1)
            OffsetArray(indsr.-x, indsr)
        end
        function Broadcast.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), x::$T, r::IdentityRange)
            indsr = axes(r, 1)
            OffsetArray(x.-indsr, indsr)
        end
    end
end
function Base.:*(r::IdentityRange, x::Number)
    indsr = axes(r, 1)
    OffsetArray(indsr.*x, indsr)
end
Base.:*(x::Number, r::IdentityRange) = r*x
function Base.:/(r::IdentityRange, x::Number)
    indsr = axes(r, 1)
    OffsetArray(indsr./x, indsr)
end
function Base.:\(x::Number, r::IdentityRange)
    indsr = axes(r, 1)
    OffsetArray(x .\ indsr, indsr)
end
Broadcast.broadcasted(::DefaultArrayStyle{1}, ::typeof(*), r::IdentityRange, x::Number) = r * x
Broadcast.broadcasted(::DefaultArrayStyle{1}, ::typeof(*), x::Number, r::IdentityRange) = r * x
Broadcast.broadcasted(::DefaultArrayStyle{1}, ::typeof(/), r::IdentityRange, x::Number) = r / x
Broadcast.broadcasted(::DefaultArrayStyle{1}, ::typeof(\), x::Number, r::IdentityRange) = x \ r

Base.collect(r::IdentityRange) = convert(Vector, first(r):last(r))
Base.sortperm(r::IdentityRange) = r
function Base.reverse(r::IdentityRange)
    indsr = axes(r, 1)
    OffsetArray(reverse(indsr), indsr)
end

Base.promote_rule(::Type{IdentityRange{T1}},::Type{IdentityRange{T2}}) where {T1,T2} =
    IdentityRange{promote_type(T1,T2)}
Base.convert(::Type{IdentityRange{T}}, r::IdentityRange{T}) where {T<:Integer} = r
Base.convert(::Type{IdentityRange{T}}, r::AbstractUnitRange) where {T<:Integer} =
    IdentityRange{T}(first(r), last(r))
Base.convert(::Type{IdentityRange}, r::AbstractUnitRange{T}) where {T<:Integer} =
    convert(IdentityRange{T}, r)

Base.show(io::IO, r::IdentityRange) = print(io, "IdentityRange(", first(r), ":", last(r), ")")

IdentityRange{R}(r::AbstractUnitRange{T}) where {R,T} = convert(IdentityRange{R}, r)
IdentityRange(r::AbstractUnitRange{T}) where {T} = convert(IdentityRange, r)

# Deprecations
import Base: +, -
@deprecate +(r::IdentityRange, x::Real) broadcast(+, r, x)
@deprecate +(x::Real, r::IdentityRange) broadcast(+, x, r)
@deprecate +(r::IdentityRange, x::Number) broadcast(+, r, x)
@deprecate +(x::Number, r::IdentityRange) broadcast(+, x, r)
@deprecate -(r::IdentityRange, x::Real) broadcast(-, r, x)
@deprecate -(x::Real, r::IdentityRange) broadcast(-, x, r)
@deprecate -(r::IdentityRange, x::Number) broadcast(-, r, x)
@deprecate -(x::Number, r::IdentityRange) broadcast(-, x, r)

end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(x^2*cosh(a+b*x+c*x^2), x) == :((1/2)*x*c^-1*sinh(a+b*x+c*x^2)+-1//4*b*c^-2*sinh(a+b*x+c*x^2)+-1//8*pi^(1/2)*c^-3//2*erfi((1/2)*c^-1//2*(b+2*c*x))*exp(a+-1//4*b^2*c^-1)+1//8*pi^(1/2)*c^-3//2*erf((1/2)*c^-1//2*(b+2*c*x))*exp(-1a+1//4*b^2*c^-1)+1//16*pi^(1/2)*b^2*c^-5//2*erf((1/2)*c^-1//2*(b+2*c*x))*exp(-1a+1//4*b^2*c^-1)+1//16*pi^(1/2)*b^2*c^-5//2*erfi((1/2)*c^-1//2*(b+2*c*x))*exp(a+-1//4*b^2*c^-1))
@test integrate(x*cosh(a+b*x+c*x^2), x) == :((1/2)*c^-1*sinh(a+b*x+c*x^2)+-1//8*b*pi^(1/2)*c^-3//2*erf((1/2)*c^-1//2*(b+2*c*x))*exp(-1a+1//4*b^2*c^-1)+-1//8*b*pi^(1/2)*c^-3//2*erfi((1/2)*c^-1//2*(b+2*c*x))*exp(a+-1//4*b^2*c^-1))
@test integrate(cosh(a+b*x+c*x^2), x) == :(1//4*pi^(1/2)*c^-1//2*erf((1/2)*c^-1//2*(b+2*c*x))*exp(-1a+1//4*b^2*c^-1)+1//4*pi^(1/2)*c^-1//2*erfi((1/2)*c^-1//2*(b+2*c*x))*exp(a+-1//4*b^2*c^-1))
@test integrate(x^-2*cosh(a+b*x+c*x^2)+-1*b*x^-1*sinh(a+b*x+c*x^2), x) == :(-1*x^-1*cosh(a+b*x+c*x^2)+(1/2)*pi^(1/2)*c^(1/2)*erfi((1/2)*c^-1//2*(b+2*c*x))*exp(a+-1//4*b^2*c^-1)+-1//2*pi^(1/2)*c^(1/2)*erf((1/2)*c^-1//2*(b+2*c*x))*exp(-1a+1//4*b^2*c^-1))
@test integrate(x^2*cosh(a+b*x+-1*c*x^2), x) == :(-1//2*x*c^-1*sinh(a+b*x+-1*c*x^2)+-1//4*b*c^-2*sinh(a+b*x+-1*c*x^2)+-1//8*pi^(1/2)*c^-3//2*erf((1/2)*c^-1//2*(b+-2*c*x))*exp(a+1//4*b^2*c^-1)+1//8*pi^(1/2)*c^-3//2*erfi((1/2)*c^-1//2*(b+-2*c*x))*exp(-1a+-1//4*b^2*c^-1)+-1//16*pi^(1/2)*b^2*c^-5//2*erf((1/2)*c^-1//2*(b+-2*c*x))*exp(a+1//4*b^2*c^-1)+-1//16*pi^(1/2)*b^2*c^-5//2*erfi((1/2)*c^-1//2*(b+-2*c*x))*exp(-1a+-1//4*b^2*c^-1))
@test integrate(x*cosh(a+b*x+-1*c*x^2), x) == :(-1//2*c^-1*sinh(a+b*x+-1*c*x^2)+-1//8*b*pi^(1/2)*c^-3//2*erf((1/2)*c^-1//2*(b+-2*c*x))*exp(a+1//4*b^2*c^-1)+-1//8*b*pi^(1/2)*c^-3//2*erfi((1/2)*c^-1//2*(b+-2*c*x))*exp(-1a+-1//4*b^2*c^-1))
@test integrate(cosh(a+b*x+-1*c*x^2), x) == :(-1//4*pi^(1/2)*c^-1//2*erf((1/2)*c^-1//2*(b+-2*c*x))*exp(a+1//4*b^2*c^-1)+-1//4*pi^(1/2)*c^-1//2*erfi((1/2)*c^-1//2*(b+-2*c*x))*exp(-1a+-1//4*b^2*c^-1))
@test integrate(x^-2*cosh(a+b*x+-1*c*x^2)+-1*b*x^-1*sinh(a+b*x+-1*c*x^2), x) == :(-1*x^-1*cosh(a+b*x+-1*c*x^2)+(1/2)*pi^(1/2)*c^(1/2)*erf((1/2)*c^-1//2*(b+-2*c*x))*exp(a+1//4*b^2*c^-1)+-1//2*pi^(1/2)*c^(1/2)*erfi((1/2)*c^-1//2*(b+-2*c*x))*exp(-1a+-1//4*b^2*c^-1))
@test integrate(x^2*cosh(1//4+x+x^2), x) == :(-1//4*sinh(1//4+x+x^2)+(1/2)*x*sinh(1//4+x+x^2)+-1//16*pi^(1/2)*erfi(1/2+x)+3//16*pi^(1/2)*erf(1/2+x))
@test integrate(x*cosh(1//4+x+x^2), x) == :((1/2)*sinh(1//4+x+x^2)+-1//8*pi^(1/2)*erf(1/2+x)+-1//8*pi^(1/2)*erfi(1/2+x))
@test integrate(cosh(1//4+x+x^2), x) == :(1//4*pi^(1/2)*erf(1/2+x)+1//4*pi^(1/2)*erfi(1/2+x))
@test integrate(x^2*cosh(a+b*x+c*x^2)^2, x) == :(1//6*x^3+-1//16*b*c^-2*sinh(2a+2*b*x+2*c*x^2)+1//8*x*c^-1*sinh(2a+2*b*x+2*c*x^2)+-1//64*2^(1/2)*pi^(1/2)*c^-3//2*erfi((1/2)*2^(1/2)*c^-1//2*(b+2*c*x))*exp(2a+-1//2*b^2*c^-1)+1//64*2^(1/2)*pi^(1/2)*c^-3//2*erf((1/2)*2^(1/2)*c^-1//2*(b+2*c*x))*exp(-2a+(1/2)*b^2*c^-1)+1//64*2^(1/2)*pi^(1/2)*b^2*c^-5//2*erf((1/2)*2^(1/2)*c^-1//2*(b+2*c*x))*exp(-2a+(1/2)*b^2*c^-1)+1//64*2^(1/2)*pi^(1/2)*b^2*c^-5//2*erfi((1/2)*2^(1/2)*c^-1//2*(b+2*c*x))*exp(2a+-1//2*b^2*c^-1))
@test integrate(x*cosh(a+b*x+c*x^2)^2, x) == :(1//4*x^2+1//8*c^-1*sinh(2a+2*b*x+2*c*x^2)+-1//32*b*2^(1/2)*pi^(1/2)*c^-3//2*erf((1/2)*2^(1/2)*c^-1//2*(b+2*c*x))*exp(-2a+(1/2)*b^2*c^-1)+-1//32*b*2^(1/2)*pi^(1/2)*c^-3//2*erfi((1/2)*2^(1/2)*c^-1//2*(b+2*c*x))*exp(2a+-1//2*b^2*c^-1))
@test integrate(cosh(a+b*x+c*x^2)^2, x) == :((1/2)*x+1//16*2^(1/2)*pi^(1/2)*c^-1//2*erf((1/2)*2^(1/2)*c^-1//2*(b+2*c*x))*exp(-2a+(1/2)*b^2*c^-1)+1//16*2^(1/2)*pi^(1/2)*c^-1//2*erfi((1/2)*2^(1/2)*c^-1//2*(b+2*c*x))*exp(2a+-1//2*b^2*c^-1))
@test integrate(x^2*cosh(a+b*x+-1*c*x^2)^2, x) == :(1//6*x^3+-1//8*x*c^-1*sinh(2a+-2*c*x^2+2*b*x)+-1//16*b*c^-2*sinh(2a+-2*c*x^2+2*b*x)+-1//64*2^(1/2)*pi^(1/2)*c^-3//2*erf((1/2)*2^(1/2)*c^-1//2*(b+-2*c*x))*exp(2a+(1/2)*b^2*c^-1)+1//64*2^(1/2)*pi^(1/2)*c^-3//2*erfi((1/2)*2^(1/2)*c^-1//2*(b+-2*c*x))*exp(-2a+-1//2*b^2*c^-1)+-1//64*2^(1/2)*pi^(1/2)*b^2*c^-5//2*erf((1/2)*2^(1/2)*c^-1//2*(b+-2*c*x))*exp(2a+(1/2)*b^2*c^-1)+-1//64*2^(1/2)*pi^(1/2)*b^2*c^-5//2*erfi((1/2)*2^(1/2)*c^-1//2*(b+-2*c*x))*exp(-2a+-1//2*b^2*c^-1))
@test integrate(x*cosh(a+b*x+-1*c*x^2)^2, x) == :(1//4*x^2+-1//8*c^-1*sinh(2a+-2*c*x^2+2*b*x)+-1//32*b*2^(1/2)*pi^(1/2)*c^-3//2*erf((1/2)*2^(1/2)*c^-1//2*(b+-2*c*x))*exp(2a+(1/2)*b^2*c^-1)+-1//32*b*2^(1/2)*pi^(1/2)*c^-3//2*erfi((1/2)*2^(1/2)*c^-1//2*(b+-2*c*x))*exp(-2a+-1//2*b^2*c^-1))
@test integrate(cosh(a+b*x+-1*c*x^2)^2, x) == :((1/2)*x+-1//16*2^(1/2)*pi^(1/2)*c^-1//2*erf((1/2)*2^(1/2)*c^-1//2*(b+-2*c*x))*exp(2a+(1/2)*b^2*c^-1)+-1//16*2^(1/2)*pi^(1/2)*c^-1//2*erfi((1/2)*2^(1/2)*c^-1//2*(b+-2*c*x))*exp(-2a+-1//2*b^2*c^-1))
@test integrate(x^2*cosh(1//4+x+x^2)^2, x) == :(-1//16*sinh(1/2+2x+2*x^2)+1//6*x^3+1//8*x*sinh(1/2+2x+2*x^2)+1//32*2^(1/2)*pi^(1/2)*erf((1/2)*2^(1/2)*(1+2x)))
@test integrate(x*cosh(1//4+x+x^2)^2, x) == :(1//4*x^2+1//8*sinh(1/2+2x+2*x^2)+-1//32*2^(1/2)*pi^(1/2)*erf((1/2)*2^(1/2)*(1+2x))+-1//32*2^(1/2)*pi^(1/2)*erfi((1/2)*2^(1/2)*(1+2x)))
@test integrate(cosh(1//4+x+x^2)^2, x) == :((1/2)*x+1//16*2^(1/2)*pi^(1/2)*erf((1/2)*2^(1/2)*(1+2x))+1//16*2^(1/2)*pi^(1/2)*erfi((1/2)*2^(1/2)*(1+2x)))
@test integrate((d+e*x)^2*cosh(a+b*x+c*x^2), x) == :((1/2)*e*c^-1*(d+e*x)*sinh(a+b*x+c*x^2)+1//4*e*c^-2*(-1*b*e+2*c*d)*sinh(a+b*x+c*x^2)+-1//8*pi^(1/2)*c^-3//2*e^2*erfi((1/2)*c^-1//2*(b+2*c*x))*exp(a+-1//4*b^2*c^-1)+1//8*pi^(1/2)*c^-3//2*e^2*erf((1/2)*c^-1//2*(b+2*c*x))*exp(-1a+1//4*b^2*c^-1)+1//16*pi^(1/2)*c^-5//2*(-1*b*e+2*c*d)^2*erf((1/2)*c^-1//2*(b+2*c*x))*exp(-1a+1//4*b^2*c^-1)+1//16*pi^(1/2)*c^-5//2*(-1*b*e+2*c*d)^2*erfi((1/2)*c^-1//2*(b+2*c*x))*exp(a+-1//4*b^2*c^-1))
@test integrate((d+e*x)*cosh(a+b*x+c*x^2), x) == :((1/2)*e*c^-1*sinh(a+b*x+c*x^2)+1//8*pi^(1/2)*c^-3//2*(-1*b*e+2*c*d)*erf((1/2)*c^-1//2*(b+2*c*x))*exp(-1a+1//4*b^2*c^-1)+1//8*pi^(1/2)*c^-3//2*(-1*b*e+2*c*d)*erfi((1/2)*c^-1//2*(b+2*c*x))*exp(a+-1//4*b^2*c^-1))
@test integrate((d+e*x)^2*cosh(a+b*x+c*x^2)^2, x) == :(1//6*e^-1*(d+e*x)^3+1//8*e*c^-1*(d+e*x)*sinh(2a+2*b*x+2*c*x^2)+1//16*e*c^-2*(-1*b*e+2*c*d)*sinh(2a+2*b*x+2*c*x^2)+-1//64*2^(1/2)*pi^(1/2)*c^-3//2*e^2*erfi((1/2)*2^(1/2)*c^-1//2*(b+2*c*x))*exp(2a+-1//2*b^2*c^-1)+1//64*2^(1/2)*pi^(1/2)*c^-5//2*(-1*b*e+2*c*d)^2*erf((1/2)*2^(1/2)*c^-1//2*(b+2*c*x))*exp(-2a+(1/2)*b^2*c^-1)+1//64*2^(1/2)*pi^(1/2)*c^-5//2*(-1*b*e+2*c*d)^2*erfi((1/2)*2^(1/2)*c^-1//2*(b+2*c*x))*exp(2a+-1//2*b^2*c^-1)+1//64*2^(1/2)*pi^(1/2)*c^-3//2*e^2*erf((1/2)*2^(1/2)*c^-1//2*(b+2*c*x))*exp(-2a+(1/2)*b^2*c^-1))
@test integrate(cosh(a+b*x+c*x^2)^2*(d+e*x), x) == :(1//4*e^-1*(d+e*x)^2+1//8*e*c^-1*sinh(2a+2*b*x+2*c*x^2)+1//32*2^(1/2)*pi^(1/2)*c^-3//2*(-1*b*e+2*c*d)*erf((1/2)*2^(1/2)*c^-1//2*(b+2*c*x))*exp(-2a+(1/2)*b^2*c^-1)+1//32*2^(1/2)*pi^(1/2)*c^-3//2*(-1*b*e+2*c*d)*erfi((1/2)*2^(1/2)*c^-1//2*(b+2*c*x))*exp(2a+-1//2*b^2*c^-1))

max_stars_repo_path: test/6 Hyperbolic functions/6.2 Hyperbolic cosine/6.2.4 (d+e x)^m cosh(a+b x+c x^2)^n.jl


text:
@testset "debiasedwhittle" begin
    @test_throws ArgumentError DebiasedWhittleLikelihood(OU,ones(10),-1)
    @test_throws ArgumentError DebiasedWhittleLikelihood(OU,ones(10,2),1)
    for (model,par,ts) in zip((OU,TwoOU,CorrelatedOU,TwoCorrelatedOU),(ones(2),ones(4),[1,1,0.5],[1,1,0.5,1,1,0.5]),(ones(10),ones(10),ones(10,2),ones(10,2)))
        w = DebiasedWhittleLikelihood(model,ts,1)
        @test_throws ArgumentError w(nothing,ones(length(par)+1),nothing,par)
        @test_throws ArgumentError w(nothing,nothing,ones(length(par),length(par)+1),par)
        @test w(par) isa Float64
        @test w(1,ones(length(par)),ones(length(par),length(par)),par) isa Float64
        @test w(nothing,ones(length(par)),ones(length(par),length(par)),par) == nothing
        @test grad(central_fdm(5, 1), w, par)[1] ≈ gradgen(w, par)
    end
end


text:
using WaveletMatrices
using IndexableBitVectors
using Base.Test
using FactCheck

srand(12345)

function seq_select(a, vec, j)
    if j ≤ 0
        return 0
    end
    for i in 1:length(vec)
        if vec[i] == a
            j -= 1
        end
        if j == 0
            return i
        end
    end
    return 0
end

facts("WaveletMatrix") do
    context("ordered") do
        wm = WaveletMatrix([0x00, 0x01, 0x02, 0x03])
        @fact length(wm) --> 4
        # getindex
        @fact wm[1] --> 0x00
        @fact wm[2] --> 0x01
        @fact wm[3] --> 0x02
        @fact wm[4] --> 0x03
        @fact typeof(wm[4]) --> UInt8
        # 0x00
        @fact rank(0x00, wm, 0) --> 0
        @fact rank(0x00, wm, 1) --> 1
        @fact rank(0x00, wm, 2) --> 1
        @fact rank(0x00, wm, 3) --> 1
        @fact rank(0x00, wm, 4) --> 1
        @fact rank(0x00, wm, 5) --> 1
        # 0x01
        @fact rank(0x01, wm, 0) --> 0
        @fact rank(0x01, wm, 1) --> 0
        @fact rank(0x01, wm, 2) --> 1
        @fact rank(0x01, wm, 3) --> 1
        @fact rank(0x01, wm, 4) --> 1
        @fact rank(0x01, wm, 5) --> 1
        # 0x02
        @fact rank(0x02, wm, 0) --> 0
        @fact rank(0x02, wm, 1) --> 0
        @fact rank(0x02, wm, 2) --> 0
        @fact rank(0x02, wm, 3) --> 1
        @fact rank(0x02, wm, 4) --> 1
        @fact rank(0x02, wm, 5) --> 1
        # 0x03
        @fact rank(0x03, wm, 0) --> 0
        @fact rank(0x03, wm, 1) --> 0
        @fact rank(0x03, wm, 2) --> 0
        @fact rank(0x03, wm, 3) --> 0
        @fact rank(0x03, wm, 4) --> 1
        @fact rank(0x03, wm, 5) --> 1

        @fact select(0x00, wm, 0) --> 0
        @fact select(0x01, wm, 0) --> 0
        @fact select(0x02, wm, 0) --> 0
        @fact select(0x03, wm, 0) --> 0

        @fact select(0x00, wm, 1) --> 1
        @fact select(0x01, wm, 1) --> 2
        @fact select(0x02, wm, 1) --> 3
        @fact select(0x03, wm, 1) --> 4

        @fact select(0x00, wm, 2) --> 0
        @fact select(0x01, wm, 2) --> 0
        @fact select(0x02, wm, 2) --> 0
        @fact select(0x03, wm, 2) --> 0

        @fact freq(0x00, wm, 1, 4) --> 1
        @fact freq(0x00, wm, 2, 4) --> 0
        @fact freq(0x01, wm, 1, 2) --> 1
        @fact freq(0x01, wm, 1, 1) --> 0
        @fact freq(0x02, wm, 4, 2) --> 0
        @fact freq(0x02, wm, 3, 3) --> 1
        @fact freq(0x03, wm, 4, 1) --> 0
        @fact freq(0x03, wm, 4, 4) --> 1
    end

    context("unordered") do
        wm = WaveletMatrix([0x01, 0x03, 0x02, 0x00])
        @fact length(wm) --> 4
        # getindex
        @fact wm[1] --> 0x01
        @fact wm[2] --> 0x03
        @fact wm[3] --> 0x02
        @fact wm[4] --> 0x00
        # 0x00
        @fact rank(0x00, wm, 0) --> 0
        @fact rank(0x00, wm, 1) --> 0
        @fact rank(0x00, wm, 2) --> 0
        @fact rank(0x00, wm, 3) --> 0
        @fact rank(0x00, wm, 4) --> 1
        @fact rank(0x00, wm, 5) --> 1
        # 0x01
        @fact rank(0x01, wm, 0) --> 0
        @fact rank(0x01, wm, 1) --> 1
        @fact rank(0x01, wm, 2) --> 1
        @fact rank(0x01, wm, 3) --> 1
        @fact rank(0x01, wm, 4) --> 1
        @fact rank(0x01, wm, 5) --> 1
        # 0x02
        @fact rank(0x02, wm, 0) --> 0
        @fact rank(0x02, wm, 1) --> 0
        @fact rank(0x02, wm, 2) --> 0
        @fact rank(0x02, wm, 3) --> 1
        @fact rank(0x02, wm, 4) --> 1
        @fact rank(0x02, wm, 5) --> 1
        # 0x03
        @fact rank(0x03, wm, 0) --> 0
        @fact rank(0x03, wm, 1) --> 0
        @fact rank(0x03, wm, 2) --> 1
        @fact rank(0x03, wm, 3) --> 1
        @fact rank(0x03, wm, 4) --> 1
        @fact rank(0x03, wm, 5) --> 1

        @fact select(0x00, wm, 1) --> 4
        @fact select(0x01, wm, 1) --> 1
        @fact select(0x02, wm, 1) --> 3
        @fact select(0x03, wm, 1) --> 2

        @fact select(0x00, wm, 2) --> 0
        @fact select(0x01, wm, 2) --> 0
        @fact select(0x02, wm, 2) --> 0
        @fact select(0x03, wm, 2) --> 0
    end

    context("homogeneous") do
        wm = WaveletMatrix([0x01, 0x01, 0x01, 0x01])
        @fact length(wm) --> 4
        # getindex
        for i in 1:4
            @fact wm[i] --> 0x01
        end
        @fact rank(0x01, wm, 0) --> 0
        @fact rank(0x01, wm, 1) --> 1
        @fact rank(0x01, wm, 2) --> 2
        @fact rank(0x01, wm, 3) --> 3
        @fact rank(0x01, wm, 4) --> 4
        @fact rank(0x01, wm, 5) --> 4

        @fact select(0x01, wm, 1) --> 1
        @fact select(0x01, wm, 2) --> 2
        @fact select(0x01, wm, 3) --> 3
        @fact select(0x01, wm, 4) --> 4
        @fact select(0x01, wm, 5) --> 0

        @fact select(0x00, wm, 1) --> 0
        @fact select(0x02, wm, 1) --> 0

        for i in 1:4, j in 1:4
            @fact freq(0x00, wm, i, j) --> 0
            @fact freq(0x01, wm, i, j) --> max(0, j - i + 1)
        end
    end

    context("Vector{Uint64}") do
        x = [
            0x0000000000000001,
            0xffffffffffffffff,
            0x0000000000000000,
            0x0000000000000001,
            0xffffffffffffffff,
            0xffffff171f410d54,
            0x0000000000000009,
            0x0000000000000000,
            0x1fffff171f410d54,
            0x0000000000000008,
        ]
        wm = WaveletMatrix(x)
        @fact length(wm) --> 10
        for i in 1:endof(x)
            @fact wm[i] --> x[i]
        end
        @fact typeof(wm[1]) --> UInt64
        for v in x, i in 1:length(x)
            @fact rank(v, wm, i) --> count(v′ -> v′ == v, x[1:i])
        end
        for v in x, j in 1:length(x)+1
            @fact select(v, wm, j) --> seq_select(v, x, j)
        end
    end

    context("2-bit encoding") do
        x = rand(0x00:0x03, 100)
        wm = WaveletMatrix{2}(x)
        @fact length(wm) --> 100
        for i in 1:endof(x)
            @fact wm[i] --> x[i]
        end
        for a in 0x00:0x03, i in 1:100
            @fact rank(a, wm, i) --> count(a′ -> a′ == a, x[1:i])
        end
        for a in 0x00:0x03, j in 1:100+1
            @fact select(a, wm, j) --> seq_select(a, x, j)
        end
    end

    context("17-bit encoding") do
        x = rand(0x00000000:0x00000011, 500)
        wm = WaveletMatrix{17}(x)
        @fact length(wm) --> 500
        for i in 1:endof(x)
            @fact wm[i] --> x[i]
        end
        for a in 0x00000000:0x00000011, i in 1:500
            @fact rank(a, wm, i) --> count(a′ -> a′ == a, x[1:i])
        end
        for a in 0x00000000:0x00000011, j in 1:500+1
            @fact select(a, wm, j) --> seq_select(a, x, j)
        end
    end

    context("destructive") do
        x = rand(0x00:0x03, 1000)
        wm = WaveletMatrix{2}(x, destructive=false)
        @fact length(x) --> 1000
        wm = WaveletMatrix{2}(x, destructive=true)
        @fact length(x) --> 0
    end

    context("random") do
        # there is at least one duplicated byte (pigeonhole principle)
        len = typemax(UInt8) + 1
        bytes = rand(UInt8, len)
        wm = WaveletMatrix(copy(bytes))
        for i in 1:len
            @fact wm[i] --> bytes[i]
        end
        for byte in 0x00:0xff, i in 1:len
            @fact rank(byte, wm, i) --> count(b -> b == byte, bytes[1:i])
        end
        for byte in 0x00:0xff, j in 0:len+1
            @fact select(byte, wm, j) --> seq_select(byte, bytes, j)
        end
        for byte in 0x00:0xff
            @fact freq(byte, wm, len, 1) --> 0
            @fact freq(byte, wm, 1, len) --> count(b -> b == byte, bytes)
            i, j = rand(1:len), rand(1:len)
            @fact freq(byte, wm, i, j) --> count(b -> b == byte, bytes[i:j])
        end
    end
end


text:
module runExamplesContactForceLaw

import Modia3D

# Test Modia3D - Examples

# examples/dynamics
collContactForceLawPath = joinpath(Modia3D.path, "test", "collision", "contactForceLaw")
include(joinpath(collContactForceLawPath, "contactForceLaw_2_boxes.jl"))
include(joinpath(collContactForceLawPath, "contactForceLaw_2_boxes_Prismatic.jl"))
include(joinpath(collContactForceLawPath, "contactForceLaw_2_boxes2.jl"))
include(joinpath(collContactForceLawPath, "contactForceLaw_Ball.jl"))
include(joinpath(collContactForceLawPath, "contactForceLaw_ballWithBall.jl"))
include(joinpath(collContactForceLawPath, "contactForceLaw_ballWithBox_45Deg.jl"))
include(joinpath(collContactForceLawPath, "contactForceLaw_ballWithBox_Prismatic.jl"))
include(joinpath(collContactForceLawPath, "contactForceLaw_Billard.jl"))
include(joinpath(collContactForceLawPath, "contactForceLaw_Billard_1Ball.jl"))
include(joinpath(collContactForceLawPath, "contactForceLaw_newtons_cradle.jl"))
include(joinpath(collContactForceLawPath, "contactForceLaw_rollingBall.jl"))
include(joinpath(collContactForceLawPath, "contactForceLaw_rollingBall_2balls.jl"))
include(joinpath(collContactForceLawPath, "contactForceLaw_rollingBall_2ballsb.jl"))

println("\n... success of runExamplesContactForceLaw.jl")
end


text:
using SnoopCompile
using SnoopCompile: countchildren

function hastv(typ)
    isa(typ, UnionAll) && return true
    if isa(typ, DataType)
        for p in typ.parameters
            hastv(p) && return true
        end
    end
    return false
end

trees = invalidation_trees(@snoopr using Revise)

function summary(trees)
    npartial = ngreater = nlesser = nambig = nequal = 0
    for methinvs in trees
        method = methinvs.method
        for fn in (:mt_backedges, :backedges)
            list = getfield(methinvs, fn)
            for item in list
                sig = nothing
                if isa(item, Pair)
                    sig = item.first
                    root = item.second
                else
                    sig = item.mi.def.sig
                    root = item
                end
                # if hastv(sig)
                #     npartial += countchildren(invtree)
                # else
                    ms1, ms2 = method.sig <: sig, sig <: method.sig
                    if ms1 && !ms2
                        ngreater += countchildren(root)
                    elseif ms2 && !ms1
                        nlesser += countchildren(root)
                    elseif ms1 && ms2
                        nequal += countchildren(root)
                    else
                        # if hastv(sig)
                        #     npartial += countchildren(root)
                        # else
                            nambig += countchildren(root)
                        # end
                    end
                # end
            end
        end
    end
    @assert nequal == 0
    println("$ngreater | $nlesser | $nambig |") # $npartial |")
end

summary(trees)


text:
## ---
using CSV: read
using DataFrames
using Plots
using StatsBase
using DataStructures: SortedSet

## ---

df=read("data/Womens Clothing E-Commerce Reviews.csv", DataFrame)

#get 70% of the rows for training
trainingdataPercent=0.70
trainingRowIds=rand(1:nrow(df),convert(Int, round(nrow(df)/trainingdataPercent)))
validationRowIds=setdiff(1:nrow(df),trainingRowIds)

@assert isempty(intersect(trainingRowIds, validationRowIds))

trainingData=df[trainingRowIds,:]
validationData=df[validationRowIds, :]

## ---
names(trainginData)

#make sure both sets has all the dept names
traingDeptName=unique(trainginData."Department Name")
validDeptName=unique(validationData."Department Name")
@assert isempty(setdiff(traingDeptName, validDeptName))

#get an idea of the frequency of the review messages for the dept
countmap(trainginData."Department Name")

reviews=hcat(trainingData.Title, trainingData."Review Text")

labels=SortedSet()
for (t, r) in (trainingData.Title, trainingData."Review Text")
    push!(labels, (t*r)...)
end
labels

max_stars_repo_name: daryoush/CudaAndFluxExperiments

max_issues_repo_name: daryoush/CudaAndFluxExperiments

max_forks_repo_name: daryoush/CudaAndFluxExperiments


text:
function BracketFinder()
    return BracketFinder(())
end

function BracketFinder(arg0::jdouble, arg1::jint)
    return BracketFinder((jdouble, jint), arg0, arg1)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_evaluations(obj::BracketFinder)
    return jcall(obj, "getEvaluations", jint, ())
end

function get_f_hi(obj::BracketFinder)
    return jcall(obj, "getFHi", jdouble, ())
end

function get_f_lo(obj::BracketFinder)
    return jcall(obj, "getFLo", jdouble, ())
end

function get_f_mid(obj::BracketFinder)
    return jcall(obj, "getFMid", jdouble, ())
end

function get_hi(obj::BracketFinder)
    return jcall(obj, "getHi", jdouble, ())
end

function get_lo(obj::BracketFinder)
    return jcall(obj, "getLo", jdouble, ())
end

function get_max_evaluations(obj::BracketFinder)
    return jcall(obj, "getMaxEvaluations", jint, ())
end

function get_mid(obj::BracketFinder)
    return jcall(obj, "getMid", jdouble, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function search(obj::BracketFinder, arg0::UnivariateFunction, arg1::GoalType, arg2::jdouble, arg3::jdouble)
    return jcall(obj, "search", void, (UnivariateFunction, GoalType, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
import Pluto: Pluto, without_pluto_file_extension, generate_html, @asynclog
using Base64
using SHA
using FromFile

@from "./MoreAnalysis.jl" import bound_variable_connections_graph
@from "./Export.jl" import try_get_exact_pluto_version, try_fromcache, try_tocache
@from "./Types.jl" import NotebookSession, RunningNotebook, FinishedNotebook
@from "./Configuration.jl" import PlutoDeploySettings
@from "./FileHelpers.jl" import find_notebook_files_recursive
myhash = base64encode ∘ sha256
function path_hash(path)
    myhash(read(path))
end

showall(xs) = Text(join(string.(xs), "\n"))

###
# Shutdown
function process(
    s::NotebookSession{String,Nothing,<:Any};
    server_session::Pluto.ServerSession,
    settings::PlutoDeploySettings,
    output_dir::AbstractString,
    start_dir::AbstractString,
    progress,
)::NotebookSession

    if s.run isa RunningNotebook
        Pluto.SessionActions.shutdown(server_session, s.run.notebook)
    end

    try
        remove_static_export(s.path; settings, output_dir)
    catch e
        @warn "Failed to remove static export files" s.path exception =
            (e, catch_backtrace())
    end

    @info "### ✓ $(progress) Shutdown complete" s.path

    NotebookSession(; path=s.path, current_hash=nothing, desired_hash=nothing, run=nothing)
end

###
# Launch
function process(
    s::NotebookSession{Nothing,String,<:Any};
    server_session::Pluto.ServerSession,
    settings::PlutoDeploySettings,
    output_dir::AbstractString,
    start_dir::AbstractString,
    progress,
)::NotebookSession

    path = s.path
    abs_path = joinpath(start_dir, path)

    @info "###### ◐ $(progress) Launching..." s.path

    jl_contents = read(abs_path, String)
    new_hash = myhash(jl_contents)
    if new_hash != s.desired_hash
        @warn "Notebook file does not have desired hash. This probably means that the file changed too quickly. Continuing and hoping for the best!" s.path new_hash s.desired_hash
    end

    keep_running = settings.SliderServer.enabled
    skip_cache = keep_running || path ∈ settings.Export.ignore_cache

    cached_state = skip_cache ? nothing : try_fromcache(settings.Export.cache_dir, new_hash)

    run = if cached_state !== nothing
        @info "Loaded from cache, skipping notebook run" s.path new_hash
        original_state = cached_state
        FinishedNotebook(; path, original_state)
    else
        try
            # open and run the notebook
            notebook = Pluto.SessionActions.open(server_session, abs_path; run_async=false)
            # get the state object
            original_state = Pluto.notebook_to_js(notebook)
            # shut down the notebook
            if !keep_running
                @info "Shutting down notebook process" s.path
                Pluto.SessionActions.shutdown(server_session, notebook)
            end
            try_tocache(settings.Export.cache_dir, new_hash, original_state)
            if keep_running
                bond_connections = bound_variable_connections_graph(notebook)
                @info "Bond connections" s.path showall(collect(bond_connections))

                RunningNotebook(; path, notebook, original_state, bond_connections)
            else
                FinishedNotebook(; path, original_state)
            end
        catch e
            (e isa InterruptException) || rethrow(e)
            @error "$(progress) Failed to run notebook!" path exception =
                (e, catch_backtrace())
            # continue
            nothing
        end
    end

    generate_static_export(
        path,
        run.original_state,
        jl_contents;
        settings,
        start_dir,
        output_dir,
    )

    @info "### ✓ $(progress) Ready" s.path new_hash

    NotebookSession(;
        path=s.path,
        current_hash=new_hash,
        desired_hash=s.desired_hash,
        run=run,
    )
end

###
# Update if needed
function process(
    s::NotebookSession{String,String,<:Any};
    server_session::Pluto.ServerSession,
    settings::PlutoDeploySettings,
    output_dir::AbstractString,
    start_dir::AbstractString,
    progress,
)::NotebookSession

    if s.current_hash != s.desired_hash
        @info "Updating notebook... will shut down and relaunch" s.path

        # Simple way to update: shut down notebook and start new one
        if s.run isa RunningNotebook
            Pluto.SessionActions.shutdown(server_session, s.run.notebook)
        end

        @info "Shutdown complete" s.path

        result = process(
            NotebookSession(;
                path=s.path,
                current_hash=nothing,
                desired_hash=s.desired_hash,
                run=nothing,
            );
            server_session,
            settings,
            output_dir,
            start_dir,
            progress,
        )

        result
    else
        s
    end
end

###
# Leave it shut down
process(s::NotebookSession{Nothing,Nothing,<:Any}; kwargs...)::NotebookSession = s

should_shutdown(::NotebookSession{String,Nothing,<:Any}) = true
should_shutdown(::NotebookSession) = false
should_update(s::NotebookSession{String,String,<:Any}) = s.current_hash != s.desired_hash
should_update(::NotebookSession) = false
should_launch(::NotebookSession{Nothing,String,<:Any}) = true
should_launch(::NotebookSession) = false

will_process(s) = should_update(s) || should_launch(s) || should_shutdown(s)

"""
Core Action: Generate static export for a Pluto Notebook

# Arguments:
1. slider_server_url: URL of the slider server. This will be the URL of your server, if you deploy
2. offer_binder: Flag to enable the Binder button
3. binder_url: URL of the binder link that will be invoked. Use a compatible pluto-enabled binder
4. baked_state: Whether to export pluto state within the html or in a separate file.
5. pluto_cdn_root: URL where pluto will go to find the static frontend assets
"""
function generate_static_export(
    path,
    original_state,
    jl_contents;
    settings,
    output_dir,
    start_dir,
)
    pluto_version = try_get_exact_pluto_version()
    export_jl_path = let
        relative_to_notebooks_dir = path
        joinpath(output_dir, relative_to_notebooks_dir)
    end
    export_html_path = let
        relative_to_notebooks_dir = without_pluto_file_extension(path) * ".html"
        joinpath(output_dir, relative_to_notebooks_dir)
    end
    export_statefile_path = let
        relative_to_notebooks_dir = without_pluto_file_extension(path) * ".plutostate"
        joinpath(output_dir, relative_to_notebooks_dir)
    end

    mkpath(dirname(export_jl_path))
    mkpath(dirname(export_html_path))
    mkpath(dirname(export_statefile_path))

    slider_server_running_somewhere =
        settings.Export.slider_server_url !== nothing ||
        (settings.SliderServer.serve_static_export_folder && settings.SliderServer.enabled)

    notebookfile_js = if settings.Export.offer_binder || slider_server_running_somewhere
        if settings.Export.baked_notebookfile
            "\"data:text/julia;charset=utf-8;base64,$(base64encode(jl_contents))\""
        else
            repr(basename(export_jl_path))
        end
    else
        "undefined"
    end
    slider_server_url_js = if slider_server_running_somewhere
        abs_path = joinpath(start_dir, path)
        url_of_root = relpath(start_dir, dirname(abs_path)) # e.g. "." or "../../.."
        repr(something(settings.Export.slider_server_url, url_of_root))
    else
        "undefined"
    end
    binder_url_js = if settings.Export.offer_binder
        repr(something(settings.Export.binder_url, Pluto.default_binder_url))
    else
        "undefined"
    end
    statefile_js = if !settings.Export.baked_state
        open(export_statefile_path, "w") do io
            Pluto.pack(io, original_state)
        end
        repr(basename(export_statefile_path))
    else
        statefile64 = base64encode() do io
            Pluto.pack(io, original_state)
        end

        "\"data:;base64,$(statefile64)\""
    end

    html_contents = generate_html(;
        pluto_cdn_root=settings.Export.pluto_cdn_root,
        version=pluto_version,
        notebookfile_js,
        statefile_js,
        slider_server_url_js,
        binder_url_js,
        disable_ui=settings.Export.disable_ui,
    )
    write(export_html_path, html_contents)

    if (settings.Export.offer_binder || settings.Export.slider_server_url !== nothing) &&
       !settings.Export.baked_notebookfile
        write(export_jl_path, jl_contents)
    end

    @debug "Written to $(export_html_path)"
end

tryrm(x) = isfile(x) && rm(x)

function remove_static_export(path; settings, output_dir)
    export_jl_path = let
        relative_to_notebooks_dir = path
        joinpath(output_dir, relative_to_notebooks_dir)
    end
    export_html_path = let
        relative_to_notebooks_dir = without_pluto_file_extension(path) * ".html"
        joinpath(output_dir, relative_to_notebooks_dir)
    end
    export_statefile_path = let
        relative_to_notebooks_dir = without_pluto_file_extension(path) * ".plutostate"
        joinpath(output_dir, relative_to_notebooks_dir)
    end

    if !settings.Export.baked_state
        tryrm(export_statefile_path)
    end
    tryrm(export_html_path)
    if (settings.Export.offer_binder || settings.Export.slider_server_url !== nothing) &&
       !settings.Export.baked_notebookfile
        tryrm(export_jl_path)
    end
end


text:
module CAOS

using BioTools
using BioTools.BLAST
using BioSequences
using JSON
using FastaIO
using Tokenize

export generate_caos_rules,
       classify_new_sequence,
       parse_tree,
       load_tree,
       get_nodes,
       remove_from_tree!,
       remove_blanks,
       get_max_depth,
       find_sequence,
       get_sPu_and_sPr,
       get_cPu_and_cPr,
       get_group_taxa_at_node,
       get_group_combos,
       classification,
       get_neighbors,
       get_first_taxa_from_tree,
       get_descendents,
       downsample_taxa,
       get_adjusted_start,
       get_all_neighbors,
       add_nodes!,
       add_blanks,
       convert_to_struct,
       classify_sequence

include("caos_functions.jl")
include("tree_functions.jl")
include("utils.jl")
include("classification.jl")
include("gap_imputation.jl")
include("user_functions.jl")

end # module


text:
# Copyright (c) 2020 Can Aknesil

# 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.

function switch_block_constraint(i, X, Y)
    """
    set_property BEL D6LUT [get_cells {apuf_driver/APUF_TEST_UNIT/SWITCH_BLOCKS[$i].SWITCH_BLOCK/EQUAL_PATHS[0].EQUAL_PATH_MUX/Q_INST_0}]
    set_property LOC SLICE_X$(X)Y$(Y) [get_cells {apuf_driver/APUF_TEST_UNIT/SWITCH_BLOCKS[$i].SWITCH_BLOCK/EQUAL_PATHS[0].EQUAL_PATH_MUX/Q_INST_0}]
    set_property BEL C6LUT [get_cells {apuf_driver/APUF_TEST_UNIT/SWITCH_BLOCKS[$i].SWITCH_BLOCK/EQUAL_PATHS[1].EQUAL_PATH_MUX/Q_INST_0}]
    set_property LOC SLICE_X$(X)Y$(Y) [get_cells {apuf_driver/APUF_TEST_UNIT/SWITCH_BLOCKS[$i].SWITCH_BLOCK/EQUAL_PATHS[1].EQUAL_PATH_MUX/Q_INST_0}]
    set_property BEL B6LUT [get_cells {apuf_driver/APUF_TEST_UNIT/SWITCH_BLOCKS[$i].SWITCH_BLOCK/EQUAL_PATHS[2].EQUAL_PATH_MUX/Q_INST_0}]
    set_property LOC SLICE_X$(X)Y$(Y) [get_cells {apuf_driver/APUF_TEST_UNIT/SWITCH_BLOCKS[$i].SWITCH_BLOCK/EQUAL_PATHS[2].EQUAL_PATH_MUX/Q_INST_0}]
    set_property BEL A6LUT [get_cells {apuf_driver/APUF_TEST_UNIT/SWITCH_BLOCKS[$i].SWITCH_BLOCK/EQUAL_PATHS[3].EQUAL_PATH_MUX/Q_INST_0}]
    set_property LOC SLICE_X$(X)Y$(Y) [get_cells {apuf_driver/APUF_TEST_UNIT/SWITCH_BLOCKS[$i].SWITCH_BLOCK/EQUAL_PATHS[3].EQUAL_PATH_MUX/Q_INST_0}]
    """
end

function first_switch_block_constraint(X, Y)
    """
    set_property BEL D6LUT [get_cells {apuf_driver/APUF_TEST_UNIT/FIRST_SWITCH_BLOCK/EQUAL_PATHS[0].EQUAL_PATH_MUX/Q_INST_0}]
    set_property LOC SLICE_X$(X)Y$(Y) [get_cells {apuf_driver/APUF_TEST_UNIT/FIRST_SWITCH_BLOCK/EQUAL_PATHS[0].EQUAL_PATH_MUX/Q_INST_0}]
    set_property BEL C6LUT [get_cells {apuf_driver/APUF_TEST_UNIT/FIRST_SWITCH_BLOCK/EQUAL_PATHS[1].EQUAL_PATH_MUX/Q_INST_0}]
    set_property LOC SLICE_X$(X)Y$(Y) [get_cells {apuf_driver/APUF_TEST_UNIT/FIRST_SWITCH_BLOCK/EQUAL_PATHS[1].EQUAL_PATH_MUX/Q_INST_0}]
    set_property BEL B6LUT [get_cells {apuf_driver/APUF_TEST_UNIT/FIRST_SWITCH_BLOCK/EQUAL_PATHS[2].EQUAL_PATH_MUX/Q_INST_0}]
    set_property LOC SLICE_X$(X)Y$(Y) [get_cells {apuf_driver/APUF_TEST_UNIT/FIRST_SWITCH_BLOCK/EQUAL_PATHS[2].EQUAL_PATH_MUX/Q_INST_0}]
    set_property BEL A6LUT [get_cells {apuf_driver/APUF_TEST_UNIT/FIRST_SWITCH_BLOCK/EQUAL_PATHS[3].EQUAL_PATH_MUX/Q_INST_0}]
    set_property LOC SLICE_X$(X)Y$(Y) [get_cells {apuf_driver/APUF_TEST_UNIT/FIRST_SWITCH_BLOCK/EQUAL_PATHS[3].EQUAL_PATH_MUX/Q_INST_0}]
    """
end

Y = map(n->string(n), 23:-1:1)
X = "66"
I = map(n->string(n), 1:23)

println(first_switch_block_constraint(X, "24"))

for i in 1:23
    println(switch_block_constraint(I[i], X, Y[i]))
end


text:
module Waveforms

export
    squarewave,
    squarewave1,
    sawtoothwave,
    sawtoothwave1,
    trianglewave,
    trianglewave1

@doc raw"""
    squarewave(x)

Compute ``2\pi``-periodic square wave of `x` with a peak amplitude ``1``.
"""
squarewave(x::Real) = ifelse(mod2pi(x) < π, 1.0, -1.0)

@doc raw"""
    squarewave(x, θ)

Compute ``2\pi``-periodic square wave of `x` with a duty cycle `θ`
and a peak amplitude ``1``.
"""
function squarewave(x::Real, θ::Real)
    0 ≤ θ ≤ 1 || throw(DomainError(θ, "squwarewave(x, θ) is only defined for 0 ≤ θ ≤ 1."))
    ifelse(mod2pi(x) < 2π * θ, 1.0, -1.0)
end

@doc raw"""
    squarewave1(x)

Compute ``1``-periodic square wave of `x` with a peak amplitude ``1``.
"""
squarewave1(x::Real) = ifelse(mod(x, 1) < 1/2, 1.0, -1.0)

@doc raw"""
    squarewave1(x, θ)

Compute ``1``-periodic square wave of `x` with a duty cycle `θ`
and a peak amplitude ``1``.
"""
function squarewave1(x::Real, θ::Real)
    0 ≤ θ ≤ 1 || throw(DomainError(θ, "squwarewave1(x, θ) is only defined for 0 ≤ θ ≤ 1."))
    ifelse(mod(x, 1) < θ, 1.0, -1.0)
end

@doc raw"""
    trianglewave(x)

Compute ``2\pi``-periodic triangle wave of `x` with a peak amplitude ``1``.
"""
function trianglewave(x::Real)
    modx = mod2pi(x + π/2)
    ifelse(modx < π, 2modx/π - 1, -2modx/π + 3)
end

@doc raw"""
    trianglewave1(x)

Compute ``1``-periodic triangle wave of `x` with a peak amplitude ``1``.
"""
function trianglewave1(x::Real)
    modx = mod(x + 1/4, 1.0)
    ifelse(modx < 1/2, 4modx - 1, -4modx + 3)
end

@doc raw"""
    sawtoothwave(x)

Compute ``2\pi``-periodic sawtooth wave of `x` with a peak amplitude ``1``.
"""
sawtoothwave(x::Real) = rem2pi(x, RoundNearest) / π

@doc raw"""
    sawtoothwave1(x)

Compute ``1``-periodic sawtooth wave of `x` with a peak amplitude ``1``.
"""
sawtoothwave1(x::Real) = rem(x, 1.0, RoundNearest) * 2

end # module


text:
# A version of Stiefel’s Conjugate Residual method for linesearch.

"""A version of Stiefel’s Conjugate Residual method for linesearch.
CR(A, b, ϵa, ϵr, itmax; quad) solves the symmetric linear system 'A * x = b'
or the least-squares problem : 'min ‖b - A * x‖²'
If quad = true, the values of the quadratic model are computed
A can be positive definite or not.
"""
function CRlin(A, b, ϵa::Float64=1e-8, ϵr::Float64=1e-6, itmax::Int=0, ε::Float64=1e-6; quad::Bool=false)
    n = size(b, 1) # size of the problem
    (size(A, 1) == n & size(A, 2) == n) || error("Inconsistent problem size")
    @info(loggerCRlin, @sprintf("CRlin: system of %d equations in %d variables", n, n))

    x = zeros(n) # initial estimation x = 0
    xNorm = 0.0
    xNorms = [xNorm] # Values of ‖x‖
    r = b # initial residual r = b - Ax = b
    rNorm = norm(r, 2) # ‖r‖
    rNorm² = rNorm * rNorm
    pr = rNorm²
    s = A * r
    ρ = dot(r, s)
    p = r
    pNorm² = rNorm²
    q = s
    ϵ = ϵa + ϵr * rNorm
    pAp = ρ # = dot(p, q) = dot(r, s)

    iter = 0
    itmax == 0 && (itmax = 2 * n)

    if quad
        m = 0.0
        mvalues = [m] # values of the quadratic model
        @info(loggerCRlin, @sprintf("%5s %7s %8s", "Iter", "‖r‖", "q"))
        @info(loggerCRlin, @sprintf("%5d %7.1e %8.1e", iter, rNorm, m))
    end

    solved = rNorm ≤ ϵ
    tired = iter ≥ itmax

    while ! (solved || tired)
        iter += 1

        if (pAp ≤ ε * pNorm²) || (ρ ≤ ε * rNorm²)
            @debug(loggerCRlin, @sprintf("nonpositive curvature detected: pAp = %8.1e and rAr = %8.1e", pAp, ρ))
            iter == 1 && return b
            return x
        end

        α = ρ / dot(q, q) # step
        x = x + α * p
        xNorm = norm(x, 2)
        push!(xNorms, xNorm)
        r = r - α * q  # residual
        rNorm² = abs(rNorm² - α * ρ)
        rNorm = sqrt(rNorm²)

        if quad
            m = -dot(b, x) + 0.5 * dot(x, A * x)
            push!(mvalues, m)
            @info(loggerCRlin, @sprintf("%5d %7.1e %8.1e", iter, rNorm, m))
        end

        solved = rNorm <= ϵ
        tired = iter >= itmax
        (solved || tired) && continue

        s = A * r
        ρbar = ρ
        ρ = dot(r, s)
        β = ρ / ρbar # step for the direction calculus
        p = r + β * p # descent direction
        pNorm² = rNorm² + 2 * β * pr - 2 * β * α * pAp + β^2 * pNorm²
        pr = rNorm² + β * pr - β * α * pAp # pᵀr
        q = s + β * q
        pAp = ρ + β^2 * pAp # dot(p, q)

    end

    return x
end


text:
"""
```
function solve_adaptive_ϕ(cloud::ParticleCloud, proposed_fixed_schedule::Vector{Float64},
                          i::Int64, j::Int64, ϕ_prop::Float64, ϕ_n1::Float64,
                          tempering_target::Float64, resampled_last_period::Bool)
```

Solves for next Φ. Returns ϕ_n, resampled_last_period, j, ϕ_prop.
"""
function solve_adaptive_ϕ(cloud::ParticleCloud, proposed_fixed_schedule::Vector{Float64},
                          i::Int64, j::Int64, ϕ_prop::Float64, ϕ_n1::Float64,
                          tempering_target::Float64, resampled_last_period::Bool)
    n_Φ = length(proposed_fixed_schedule)

    if resampled_last_period
        # The ESS_bar is reset to target an evenly weighted particle population
        ESS_bar = tempering_target*length(cloud)
        resampled_last_period = false
    else
        ESS_bar = tempering_target*cloud.ESS[i-1]
    end

    # Setting up the optimal ϕ solving function for endogenizing the tempering schedule
    optimal_ϕ_function(ϕ) = compute_ESS(get_loglh(cloud), get_weights(cloud), ϕ, ϕ_n1,
                                        old_loglh = get_old_loglh(cloud)) - ESS_bar

    # Find a ϕ_prop such that the optimal ϕ_n will lie somewhere between ϕ_n1 and ϕ_prop
    # Do so by iterating through a proposed_fixed_schedule and finding the first
    # ϕ_prop such that the ESS would fall by more than the targeted amount ESS_bar
    while optimal_ϕ_function(ϕ_prop) >= 0 && j <= n_Φ
        ϕ_prop = proposed_fixed_schedule[j]
        j += 1
    end

    # Note: optimal_ϕ_function(ϕ_n1) > 0 because ESS_{t-1} is always positive
    # When ϕ_prop != 1.0, there are still ϕ increments strictly below 1 that
    # give the optimal ϕ step, ϕ_n.
    # When ϕ_prop == 1.0 but optimal_ϕ_function(ϕ_prop) < 0, there still exists
    # an optimal ϕ step, ϕ_n, that does not equal 1.
    # Thus the interval [optimal_ϕ_function(ϕ_n1), optimal_ϕ_function(ϕ_prop)] always
    # contains a 0 by construction.

    # Modification makes it such that ϕ_n is the minimum of ϕ_prop (the fixed schedule)
    # at a given stage or the root-solved ϕ such that the ESS drops by the target amount.
    # Thus the ϕ_schedule should be strictly bounded above by the fixed schedule
    # i.e. the adaptive ϕ schedule should not outpace the fixed schedule at the end
    # (when the fixed schedule tends to drop by less than 5% per iteration)

    if ϕ_prop != 1. || optimal_ϕ_function(ϕ_prop) < 0
        ϕ_n = fzero(optimal_ϕ_function, [ϕ_n1, ϕ_prop], xtol = 0.)
        push!(cloud.tempering_schedule, ϕ_n)
    else
        ϕ_n = 1.
        push!(cloud.tempering_schedule, ϕ_n)
    end

    return ϕ_n, resampled_last_period, j, ϕ_prop
end

"""
```
function solve_adaptive_ϕ(cloud::Cloud, proposed_fixed_schedule::Vector{Float64},
                          i::Int64, j::Int64, ϕ_prop::Float64, ϕ_n1::Float64,
                          tempering_target::Float64, resampled_last_period::Bool)
```

Solves for next Φ. Returns ϕ_n, resampled_last_period, j, ϕ_prop.
"""
function solve_adaptive_ϕ(cloud::Cloud, proposed_fixed_schedule::Vector{Float64},
                          i::Int64, j::Int64, ϕ_prop::Float64, ϕ_n1::Float64,
                          tempering_target::Float64, resampled_last_period::Bool)
    n_Φ = length(proposed_fixed_schedule)

    if resampled_last_period
        # The ESS_bar is reset to target an evenly weighted particle population
        ESS_bar = tempering_target * length(cloud)
        resampled_last_period = false
    else
        ESS_bar = tempering_target*cloud.ESS[i-1]
    end

    # Setting up the optimal ϕ solving function for endogenizing the tempering schedule
    optimal_ϕ_function(ϕ) = compute_ESS(get_loglh(cloud), get_weights(cloud), ϕ, ϕ_n1,
                                        old_loglh = get_old_loglh(cloud)) - ESS_bar

    # Find a ϕ_prop such that the optimal ϕ_n will lie somewhere between ϕ_n1 and ϕ_prop
    # Do so by iterating through a proposed_fixed_schedule and finding the first
    # ϕ_prop such that the ESS would fall by more than the targeted amount ESS_bar
    while optimal_ϕ_function(ϕ_prop) >= 0 && j <= n_Φ
        ϕ_prop = proposed_fixed_schedule[j]
        j += 1
    end

    # Note: optimal_ϕ_function(ϕ_n1) > 0 because ESS_{t-1} is always positive
    # When ϕ_prop != 1.0, there are still ϕ increments strictly below 1 that
    # give the optimal ϕ step, ϕ_n.
    # When ϕ_prop == 1.0 but optimal_ϕ_function(ϕ_prop) < 0, there still exists
    # an optimal ϕ step, ϕ_n, that does not equal 1.
    # Thus the interval [optimal_ϕ_function(ϕ_n1), optimal_ϕ_function(ϕ_prop)] always
    # contains a 0 by construction.

    # Modification makes it such that ϕ_n is the minimum of ϕ_prop (the fixed schedule)
    # at a given stage or the root-solved ϕ such that the ESS drops by the target amount.
    # Thus the ϕ_schedule should be strictly bounded above by the fixed schedule
    # i.e. the adaptive ϕ schedule should not outpace the fixed schedule at the end
    # (when the fixed schedule tends to drop by less than 5% per iteration)

    if ϕ_prop != 1. || optimal_ϕ_function(ϕ_prop) < 0
        ϕ_n = fzero(optimal_ϕ_function, [ϕ_n1, ϕ_prop], xtol = 0.)
        push!(cloud.tempering_schedule, ϕ_n)
    else
        ϕ_n = 1.
        push!(cloud.tempering_schedule, ϕ_n)
    end

    return ϕ_n, resampled_last_period, j, ϕ_prop
end

"""
```
mvnormal_mixture_draw(θ_old, σ; cc, α, θ_prop) where {T<:AbstractFloat}
```

Create a `DegenerateMvNormal` distribution object, `d`, from a parameter vector, `p`, and a
standard deviation matrix (obtained from SVD), `σ`.

Generate a draw from the mixture distribution of:
1. A `DegenerateMvNormal` centered at θ_old with the standard deviation matrix `σ`, scaled by `cc^2` and with mixture proportion `α`.
2. A `DegenerateMvNormal` centered at the same mean, but with a standard deviation matrix of the diagonal entries of `σ` scaled by `cc^2` with mixture proportion `(1 - α)/2`.
3. A `DegenerateMvNormal`  with the same standard deviation matrix `σ` but centered at the new proposed mean, `θ_prop`, scaled by `cc^2`, and with mixture proportion `(1 - α)/2`.

If no `θ_prop` is given, but an `α` is specified, then the mixture will consist of `α` of
the standard distribution and `(1 - α)` of the diagonalized distribution.

### Arguments
- `θ_old::Vector{T}`: The mean of the desired distribution
- `σ::Matrix{T}`: The standard deviation matrix of the desired distribution

### Keyword Arguments
- `cc::T`: The standard deviation matrix scaling factor
- `α::T`: The mixing proportion
- `θ_prop::Vector{T}`: The proposed parameter vector to be used as part of the mixture distribution, set by default to be the weighted mean of the particles, prior to mutation.

### Outputs
- `θ_new::Vector{T}`: The draw from the mixture distribution to be used as the MH proposed step
- `new mixture_density::T`: The mixture density conditional on θ_old evaluated at `θ_new` to be used in calculating the MH move probability
- `old mixture_density::T`: The mixture density conditional on θ_new evaluated at `θ_old` to be used in calculating the MH move probability

"""
function mvnormal_mixture_draw(θ_old::Vector{T}, d_prop::Distribution;
                               c::T = 1.0, α::T = 1.0) where T<:AbstractFloat
    @assert 0 <= α <= 1
    d_bar = MvNormal(d_prop.μ, c^2 * d_prop.Σ)

    # Create mixture distribution conditional on the previous parameter value, θ_old
    d_old      = MvNormal(θ_old, c^2 * d_prop.Σ)
    d_diag_old = MvNormal(θ_old, diagm(0 => diag(c^2 * d_prop.Σ))) # FORTRAN: multiplies by 3
    d_mix_old  = MixtureModel(MvNormal[d_old, d_diag_old, d_bar], [α, (1 - α)/2, (1 - α)/2])

    θ_new = rand(d_mix_old)

    # Create mixture distribution conditional on the new parameter value, θ_new
    d_new      = MvNormal(θ_new, c^2 * d_prop.Σ)
    d_diag_new = MvNormal(θ_new, diagm(0 => diag(c^2 * d_prop.Σ)))
    d_mix_new  = MixtureModel(MvNormal[d_new, d_diag_new, d_bar], [α, (1 - α)/2, (1 - α)/2])

    # The density of θ_new | θ_old and the density of θ_old | θ_new
    # taken with respect to their respective mixture densities
    new_mix_density = logpdf(d_mix_old, θ_new)
    old_mix_density = logpdf(d_mix_new, θ_old)

    return θ_new, new_mix_density, old_mix_density
end

"""
```
function compute_ESS(loglh::Vector{T}, current_weights::Vector{T}, ϕ_n::T, ϕ_n1::T;
                     old_loglh::Vector{T} = zeros(length(loglh))) where {T<:AbstractFloat}
```
Compute ESS given log likelihood, current weights, ϕ_n, ϕ_{n-1}, and old log likelihood.
"""
function compute_ESS(loglh::Vector{T}, current_weights::Vector{T}, ϕ_n::T, ϕ_n1::T;
                     old_loglh::Vector{T} = zeros(length(loglh))) where T<:AbstractFloat
    inc_weights  = exp.((ϕ_n1 - ϕ_n) * old_loglh + (ϕ_n - ϕ_n1) * loglh)
    new_weights  = current_weights .* inc_weights
    norm_weights = new_weights / sum(new_weights)
    ESS          = 1 / sum(norm_weights .^ 2)
    return ESS
end

"""
```
function vector_particles_to_cloud(m::AbstractDSGEModel, particles::Vector{Particles})
```
Converts a ParticleCloud (old) to Cloud (new).
"""
function vector_particles_to_cloud(m::AbstractDSGEModel, particles::Vector{Particle})
    cloud = Cloud(m, length(particles))
    N = size(cloud.particles, 2)
    for i in 1:size(cloud.particles,1)
        cloud.particles[i, 1:length(m.parameters)] = particles[i].value
        cloud.particles[i, ind_loglh(N)] = particles[i].loglh
        cloud.particles[i, ind_logpost(N)] = particles[i].logpost
        update!(m, particles[i].value)
        cloud.particles[i, ind_logprior(N)] = prior(m)
        cloud.particles[i, ind_old_loglh(N)] = particles[i].old_loglh
        cloud.particles[i, ind_accept(N)] = particles[i].accept
        cloud.particles[i, ind_weight(N)] = particles[i].weight
    end
    return cloud
end


text:
export ImageView

mutable struct ImageView
    mImage::Image
    mCreateInfo::vk.VkImageViewCreateInfo
    mHandle::vk.VkImageView

    function ImageView(ofImg, type::vk.VkImageViewType)
        imgInfo = ofImg.mCreateInfo;
        range = vk.VkImageSubresourceRange(
            aspectsOf(imgInfo.format), #aspectMask::VkImageAspectFlags
            0, #baseMipLevel::UInt32
            imgInfo.mipLevels, #levelCount::UInt32
            0, #baseArrayLayer::UInt32
            imgInfo.arrayLayers #layerCount::UInt32
        )
        return ImageView(ofImg, type, range)
    end

    function ImageView(ofImg, type::vk.VkImageViewType, range::vk.VkImageSubresourceRange)
        this = new()
        this.mImage = ofImg
        imgInfo = ofImg.mCreateInfo
        components = vk.VkComponentMapping(
            vk.VK_COMPONENT_SWIZZLE_IDENTITY, #r::VkComponentSwizzle
            vk.VK_COMPONENT_SWIZZLE_IDENTITY, #g::VkComponentSwizzle
            vk.VK_COMPONENT_SWIZZLE_IDENTITY, #b::VkComponentSwizzle
            vk.VK_COMPONENT_SWIZZLE_IDENTITY #a::VkComponentSwizzle
        )
        this.mCreateInfo = vk.VkImageViewCreateInfo(
            vk.VK_STRUCTURE_TYPE_IMAGE_VIEW_CREATE_INFO, #sType::VkStructureType
            C_NULL, #pNext::Ptr{Cvoid}
            0, #flags::VkImageViewCreateFlags
            handle(ofImg), #image::VkImage
            type, #viewType::VkImageViewType
            imgInfo.format, #format::VkFormat
            components, #components::VkComponentMapping
            range, #subresourceRange::VkImageSubresourceRange
        )

        newImageView = Ref{vk.VkImageView}()
        infoRef = Ref(this.mCreateInfo)
        GC.@preserve infoRef begin
            if (vk.vkCreateImageView(getLogicalDeviceOf(this.mImage), infoRef, C_NULL, newImageView) != vk.VK_SUCCESS)
                error("Failed to create image view!")
            end
        end
        this.mHandle = newImageView[]
        return this
    end
end

@class ImageView [:handle, :destroy]

function destroy(this::ImageView)
    vk.vkDestroyImageView(getLogicalDeviceOf(this.mImage), this.mHandle, C_NULL)
end

function createView(image::Image, type::vk.VkImageViewType, range::vk.VkImageSubresourceRange = vk.VkImageSubresourceRange(0, 0, 0, 0, 0))
    if range == vk.VkImageSubresourceRange(0, 0, 0, 0, 0)
        return ImageView(image, type)
    end
    return ImageView(image, type, range)
end

function handle(this::ImageView)::vk.VkImageView
    return this.mHandle
end

function getImage(this::ImageView)::Image
    return this.mImage
end

function getLayers(this::ImageView)::UInt32
    return this.mCreateInfo.subresourceRange.layerCount
end