text:
import Base:writemime

function html_body(app)
    """

      <title>webRTC.jl</title>
      <meta name="viewport" content="width=device-width; initial-scale=1;">
      <meta name="description" content="Transfer files(of any size) and messages in the browser using WebRTC.">
      <meta name="author" content="erbbysam">

      <link href="bootstrap/css/bootstrap.css" rel="stylesheet"> </link>

      <link type="text/css" href="style.css" rel="stylesheet"></link>

            <button type="button" class="navbar-toggle" data-toggle="collapse" data-target=".navbar-collapse">

            </button>
            webRTC.jl

To get started create a new room:

              <form id="webrtc_room_form" class="form-signin">
                <input id="existing" type="text" placeholder="room name" value="" />
                <button class="btn btn-lg btn-primary btn-block connect-btn" type="submit">Create</button>

              </form>

              or

Generate a random room name

            <form id="webrtc_input_form" class="form-signin">

              <input id="username" type="text" class="form-control" placeholder="username" value="" autofocus>

              <button class="btn btn-lg btn-primary btn-block connect-btn" type="submit">connect</button>
            </form>

              <button type="button" class="btn btn-primary btn-xs btn_off" data-toggle="offcanvas">Toggle User List</button>

<input type="file" id="select_file" name="select_file" class="small" />

                <input id="chatinput" style="display:none;" type="text" placeholder="Message:"/>

                  Please wait a few seconds to connect.

    """
end

function Base.writemime(io::IO, ::MIME"text/html", app::WebRTCTextApplication)
    print(io, html_body(app))
end


text:
module Merly
import Base.|
using Sockets
using HTTP
using MbedTLS

|(x::String, y::String)= string(x,"|",y)

GET="GET"
POST="POST"
PUT="PUT"
DELETE="DELETE"
HEAD = "HEAD"
OPTIONS = "OPTIONS"
PATCH = "PATCH"
CONNECT	= "CONNECT"
TRACE = "TRACE"

tonumber = Dict{String,Char}(
 "CONNECT" => '1'
,"TRACE"   => '2'
,"HEAD"    => '3'
,"DELETE"  => '4'
,"PUT"     => '5'
,"OPTIONS" => '6'
,"POST"    => '7'
,"GET"     => '8'
,"PATCH"   => '9'
)

# Dict("g" => (g(x::Int) = x + 5)) #mas rapido
# Dict("g" => function g(x::Int); x + 5; end)

myendpoints = Dict{Int64,Array{NamedTuple{(:route, :toexec, :urlparams),Tuple{Union{String,Regex},Function,Union{Nothing,Dict{Int64,String}}}},1}}(

0 => [(route= "", toexec= function nf(request,HTTP); HTTP.Response(404); end , urlparams= nothing)]

)

constantheaders  = Pair{String,String}[]
CORSenabled = Union{Nothing,Function}[nothing]

include("utilities.jl")
include("files.jl")
include("core.jl")
include("mimetypes.jl")
include("routes.jl")
include("allformats.jl")

addnotfound("NotFound",myendpoints)

export webserverfiles, webserverpath, notfound, headersalways, File, @page, @route, GET,POST,PUT,DELETE,HEAD,OPTIONS,PATCH,Get,Post,Put,Delete,Connect,Trace,Head,Patch,Options,start,useCORS,formats

end # module

ext: jl


text:
using Test
import Pluto
import Pluto: update_run!, WorkspaceManager, ClientSession, ServerSession, Notebook, Cell
import Distributed

@testset "Bonds" begin

    🍭 = ServerSession()
    🍭.options.evaluation.workspace_use_distributed = false
    fakeclient = ClientSession(:fake, nothing)
    🍭.connected_clients[fakeclient.id] = fakeclient

    @testset "AbstractPlutoDingetjes.jl" begin
        🍭.options.evaluation.workspace_use_distributed = true
        notebook = Notebook([
                # 1
                Cell("""
                begin
                    import AbstractPlutoDingetjes
                    const APD = AbstractPlutoDingetjes
                    import AbstractPlutoDingetjes.Bonds
                end
                """),
                # 2
                Cell("""
                begin
                    struct HTMLShower
                        f::Function
                    end
                    Base.show(io::IO, m::MIME"text/html", hs::HTMLShower) = hs.f(io)
                end
                """),
                Cell("""
                APD.is_inside_pluto()
                """),
                # 4
                Cell("""
                HTMLShower() do io
                    write(io, string(APD.is_inside_pluto(io)))
                end
                """),
                Cell("""
                HTMLShower() do io
                    write(io, string(APD.is_supported_by_display(io, Bonds.initial_value)))
                end
                """),
                Cell("""
                HTMLShower() do io
                    write(io, string(APD.is_supported_by_display(io, Bonds.transform_value)))
                end
                """),
                Cell("""
                HTMLShower() do io
                    write(io, string(APD.is_supported_by_display(io, Bonds.validate_value)))
                end
                """),
                Cell("""
                HTMLShower() do io
                    write(io, string(APD.is_supported_by_display(io, sqrt)))
                end
                """),
                # 9
                Cell("""
                @bind x_simple html"<input type=range>"
                """),
                Cell("""
                x_simple
                """),
                # 11
                Cell("""
                begin
                    struct OldSlider end

                    Base.show(io::IO, m::MIME"text/html", os::OldSlider) = write(io, "<input type=range value=1>")

                    Base.get(os::OldSlider) = 1
                end
                """),
                Cell("""
                @bind x_old OldSlider()
                """),
                Cell("""
                x_old
                """),
                # 14
                Cell("""
                begin
                    struct NewSlider end

                    Base.show(io::IO, m::MIME"text/html", os::NewSlider) = write(io, "<input type=range value=1>")

                    Bonds.initial_value(os::NewSlider) = 1

                    Bonds.possible_values(s::NewSlider) = [1,2,3]
                end
                """),
                Cell("""
                @bind x_new NewSlider()
                """),
                Cell("""
                x_new
                """),
                # 17
                Cell("""
                begin
                    struct TransformSlider end

                    Base.show(io::IO, m::MIME"text/html", os::TransformSlider) = write(io, "<input type=range value=1>")

                    Bonds.initial_value(os::TransformSlider) = "x"
                    Bonds.possible_values(os::TransformSlider) = 1:10
                    Bonds.transform_value(os::TransformSlider, from_js) = repeat("x", from_js)
                end
                """),
                Cell("""
                @bind x_transform TransformSlider()
                """),
                Cell("""
                x_transform
                """),
                # 20
                Cell("""
                begin
                    struct BadTransformSlider end

                    Base.show(io::IO, m::MIME"text/html", os::BadTransformSlider) = write(io, "<input type=range value=1>")

                    Bonds.initial_value(os::BadTransformSlider) = "x"
                    Bonds.possible_values(os::BadTransformSlider) = 1:10
                    Bonds.transform_value(os::BadTransformSlider, from_js) = error("bad")
                end
                """),
                Cell("""
                @bind x_badtransform BadTransformSlider()
                """),
                Cell("""
                x_badtransform
                """),
                # 23
                Cell("""
                count = Ref(0)
                """),
                Cell("""
                @bind x_counter NewSlider() #or OldSlider(), same idea
                """),
                Cell("""
                let
                    x_counter
                    count[] += 1
                end
                """),
                # 26
                Cell("""
                @assert x_old == 1
                """),
                Cell("""
                @assert x_new == 1
                """),
                Cell("""
                @assert x_transform == "x"
                """),
                # 29
                Cell("""
                begin
                    struct PossibleValuesTest
                        possible_values::Any
                    end

                    Base.show(io::IO, m::MIME"text/html", ::PossibleValuesTest) = write(io, "hello")

                    Bonds.possible_values(pvt::PossibleValuesTest) = pvt.possible_values
                end
                """),
                Cell("@bind pv1 PossibleValuesTest(Bonds.NotGiven())"),
                Cell("@bind pv2 PossibleValuesTest(Bonds.InfinitePossibilities())"),
                Cell("@bind pv3 PossibleValuesTest([1,2,3])"),
                Cell("@bind pv4 PossibleValuesTest((x+1 for x in 1:10))"),
                # 34
                Cell("@bind pv5 PossibleValuesTest(1:10)"),
            ])
        fakeclient.connected_notebook = notebook

        function set_bond_value(name, value, is_first_value=false)
            notebook.bonds[name] = Dict("value" => value)
            Pluto.set_bond_values_reactive(;
                session=🍭,
                notebook=notebook,
                bound_sym_names=[name],
                is_first_values=[is_first_value],
                run_async=false,
            )
        end

        # before loading AbstractPlutoDingetjes, test the default behaviour:
        update_run!(🍭, notebook, notebook.cells[9:10])
        @test noerror(notebook.cells[9])
        @test noerror(notebook.cells[10])

        @test Pluto.possible_bond_values(🍭, notebook, :x_simple) == :NotGiven
        @test notebook.cells[10].output.body == "missing"
        set_bond_value(:x_simple, 1, true)
        @test notebook.cells[10].output.body == "1"

        update_run!(🍭, notebook, notebook.cells)

        @test noerror(notebook.cells[1])
        @test noerror(notebook.cells[2])
        @test noerror(notebook.cells[3])
        @test noerror(notebook.cells[4])
        @test noerror(notebook.cells[5])
        @test noerror(notebook.cells[6])
        @test noerror(notebook.cells[7])
        @test noerror(notebook.cells[8])

        @test notebook.cells[3].output.body == "true"
        @test notebook.cells[4].output.body == "true"
        @test notebook.cells[5].output.body == "true"
        @test notebook.cells[6].output.body == "true"
        @test notebook.cells[7].output.body == "false"
        @test notebook.cells[8].output.body == "false"

        @test noerror(notebook.cells[9])
        @test noerror(notebook.cells[10])
        @test noerror(notebook.cells[11])
        @test noerror(notebook.cells[12])
        @test noerror(notebook.cells[13])
        @test noerror(notebook.cells[14])
        @test noerror(notebook.cells[15])
        @test noerror(notebook.cells[16])
        @test noerror(notebook.cells[17])
        @test noerror(notebook.cells[18])
        @test noerror(notebook.cells[19])
        @test noerror(notebook.cells[20])
        @test noerror(notebook.cells[21])
        @test noerror(notebook.cells[22])
        @test noerror(notebook.cells[23])
        @test noerror(notebook.cells[24])
        @test noerror(notebook.cells[25])
        @test noerror(notebook.cells[26])
        @test noerror(notebook.cells[27])
        @test noerror(notebook.cells[28])
        @test noerror(notebook.cells[29])
        @test noerror(notebook.cells[30])
        @test noerror(notebook.cells[31])
        @test noerror(notebook.cells[32])
        @test noerror(notebook.cells[33])
        @test noerror(notebook.cells[34])
        @test length(notebook.cells) == 34

        @test Pluto.possible_bond_values(🍭, notebook, :x_new) == [1,2,3]
        @test_throws Exception Pluto.possible_bond_values(🍭, notebook, :asdfasdfx_new)
        @test Pluto.possible_bond_values(🍭, notebook, :pv1) == :NotGiven
        @test Pluto.possible_bond_values(🍭, notebook, :pv2) == :InfinitePossibilities
        @test Pluto.possible_bond_values(🍭, notebook, :pv3) == [1,2,3]
        @test Pluto.possible_bond_values(🍭, notebook, :pv4) == 2:11
        @test Pluto.possible_bond_values(🍭, notebook, :pv5) === 1:10
        @test Pluto.possible_bond_values_length(🍭, notebook, :pv1) == :NotGiven
        @test Pluto.possible_bond_values_length(🍭, notebook, :pv2) == :InfinitePossibilities
        @test Pluto.possible_bond_values_length(🍭, notebook, :pv3) == 3
        @test Pluto.possible_bond_values_length(🍭, notebook, :pv4) == 10
        @test Pluto.possible_bond_values_length(🍭, notebook, :pv5) == 10

        @test notebook.cells[10].output.body == "missing"
        set_bond_value(:x_simple, 1, true)
        @test notebook.cells[10].output.body == "1"

        @test notebook.cells[13].output.body == "1"
        set_bond_value(:x_old, 1, true)
        @test notebook.cells[13].output.body == "1"
        set_bond_value(:x_old, 99, false)
        @test notebook.cells[13].output.body == "99"

        @test notebook.cells[16].output.body == "1"
        set_bond_value(:x_new, 1, true)
        @test notebook.cells[16].output.body == "1"
        set_bond_value(:x_new, 99, false)
        @test notebook.cells[16].output.body == "99"

        @test notebook.cells[19].output.body == "\"x\""
        set_bond_value(:x_transform, 1, true)
        @test notebook.cells[19].output.body == "\"x\""
        set_bond_value(:x_transform, 3, false)
        @test notebook.cells[19].output.body == "\"xxx\""

        @test notebook.cells[22].output.body != "missing"
        @info "The following error is expected:"
        set_bond_value(:x_badtransform, 1, true)
        @test notebook.cells[22].output.body != "missing"

        @test notebook.cells[25].output.body == "1"
        set_bond_value(:x_counter, 1, true)
        @test notebook.cells[25].output.body == "1"
        set_bond_value(:x_counter, 7, false)
        @test notebook.cells[25].output.body == "2"

        WorkspaceManager.unmake_workspace((🍭, notebook))
        🍭.options.evaluation.workspace_use_distributed = false

        # test that the notebook file is runnable:

        test_proc = Distributed.addprocs(1)[1]

        Distributed.remotecall_eval(Main, test_proc, quote
            import Pkg
            try
                Pkg.UPDATED_REGISTRY_THIS_SESSION[] = true
            catch; end
            Pkg.activate(mktempdir())
            Pkg.add("AbstractPlutoDingetjes")
        end)
        @test Distributed.remotecall_eval(Main, test_proc, quote
            include($(notebook.path))
            true
        end)
        Distributed.rmprocs(test_proc)
    end
end


text:
module ArrayFire

using Libdl,Random,SparseArrays,LinearAlgebra,FFTW,DSP,Statistics

include("common.jl")
include("array.jl")
include("util.jl")
include("wrap.jl")
include("random.jl")
include("graphics.jl")
include("indexing.jl")
include("fft.jl")

end


text:
# This file is a part of BinaryCommitteeMachineFBP.jl. License is MIT: http://github.com/carlobaldassi/BinaryCommitteeMachineFBP.jl/LICENCE.md

using StatsFuns

primitive type MagT64 <: Mag64 64 end

f2mT(a::Float64) = f2m(MagT64, a)

include("AtanhErf.jl")
using .AtanhErf

const mInf = 30.0

magformat(::Type{MagT64}) = :tanh
parseinner(::Type{Val{:tanh}}, s::AbstractString) = mtanh(MagT64, parse(Float64, s))

convert(::Type{MagT64}, y::Float64) = f2mT(clamp(atanh(y), -mInf, mInf))
convert(::Type{Float64}, y::MagT64) = tanh(m2f(y))

MagT64(y::Float64) = convert(MagT64, y)

forcedmag(::Type{MagT64}, y::Float64) = f2mT(atanh(y))

mtanh(::Type{MagT64}, x::Float64) = f2mT(x)
atanh(x::MagT64) = m2f(x)

MagT64(pp::Real, pm::Real) = f2mT(clamp((log(pp) - log(pm)) / 2, -mInf, mInf))

isfinite(a::MagT64) = !isnan(m2f(a))

⊗(a::MagT64, b::MagT64) = f2mT(m2f(a) + m2f(b))
function ⊘(a::MagT64, b::MagT64)
    xa = m2f(a)
    xb = m2f(b)
    return f2mT(ifelse(xa == xb, 0.0, xa - xb)) # NOTE: the ifelse is for the Inf case
end

reinforce(m0::MagT64, γ::Float64) = f2mT(m2f(m0) * γ)

damp(newx::MagT64, oldx::MagT64, λ::Float64) = f2mT(m2f(newx) * (1 - λ) + m2f(oldx) * λ)

lr(x::Float64) = log1p(exp(-2abs(x)))
log2cosh(x::Float64) = abs(x) + lr(x)

function (*)(x::MagT64, y::MagT64)
    ax = m2f(x)
    ay = m2f(y)

    if ax ≥ ay && ax ≥ -ay
        t1 = 2ay
    elseif ax ≥ ay && ax < -ay
        t1 = -2ax
    elseif ax < ay && ax ≥ -ay
        t1 = 2ax
    else # ax < ay && ax < -ay
        t1 = -2ay
    end

    t2 = isinf(ax) || isinf(ay) ?
         0.0 : lr(ax + ay) - lr(ax - ay)

    return f2mT((t1 + t2) / 2)
end

merf(::Type{MagT64}, x::Float64) = f2mT(atanherf(x))

function auxmix(H::MagT64, a₊::Float64, a₋::Float64)
    aH = m2f(H)

    aH == 0.0 && return f2mT(0.0)

    xH₊ = aH + a₊
    xH₋ = aH + a₋

    # we need to compute
    #   t1 = abs(xH₊) - abs(a₊) - abs(xH₋) + abs(a₋)
    #   t2 = lr(xH₊) - lr(a₊) - lr(xH₋) + lr(a₋)
    # but we also need to take into account infinities
    if isinf(aH)
        if !isinf(a₊) && !isinf(a₋)
            t1 = sign(aH) * (a₊ - a₋) - abs(a₊) + abs(a₋)
            t2 = -lr(a₊) + lr(a₋)
        elseif isinf(a₊) && !isinf(a₋)
            if sign(a₊) == sign(aH)
                t1 = -sign(aH) * (a₋) + abs(a₋)
                t2 = lr(a₋)
            else
                t1 = -2mInf
                t2 = 0.0
            end
        elseif !isinf(a₊) && isinf(a₋)
            if sign(a₋) == sign(aH)
                t1 = sign(aH) * (a₊) - abs(a₊)
                t2 = -lr(a₊)
            else
                t1 = 2mInf
                t2 = 0.0
            end
        else # isinf(a₊) && isinf(a₋)
            if (sign(a₊) == sign(aH) && sign(a₋) == sign(aH)) || (sign(a₊) ≠ sign(aH) && sign(a₋) ≠ sign(aH))
                t1 = 0.0
                t2 = 0.0
            elseif sign(a₊) == sign(aH) # && sign(a₋) ≠ sign(aH)
                t1 = 2mInf
                t2 = 0.0
            else # sign(a₋) == sign(aH) && sign(a₊) ≠ sign(aH)
                t1 = -2mInf
                t2 = 0.0
            end
        end
    else # !isinf(aH)
        t1 = 0.0
        t1 += isinf(a₊) ? 0.0 : abs(xH₊) - abs(a₊)
        t1 -= isinf(a₋) ? 0.0 : abs(xH₋) - abs(a₋)
        t2 = lr(xH₊) - lr(a₊) - lr(xH₋) + lr(a₋)
    end

    return f2mT((t1 + t2) / 2)
end

exactmix(H::MagT64, p₊::MagT64, p₋::MagT64) = auxmix(H, m2f(p₊), m2f(p₋))

function erfmix(H::MagT64, m₊::Float64, m₋::Float64)
    aerf₊ = atanherf(m₊)
    aerf₋ = atanherf(m₋)
    return auxmix(H, aerf₊, aerf₋)
end

# log((1 + x * y) / 2)
function log1pxy(x::MagT64, y::MagT64)
    ax = m2f(x)
    ay = m2f(y)

    return !isinf(ax) && !isinf(ay) ? abs(ax + ay) - abs(ax) - abs(ay) + lr(ax + ay) - lr(ax) - lr(ay) :
            isinf(ax) && !isinf(ay) ? sign(ax) * ay - abs(ay) - lr(ay) :
           !isinf(ax) &&  isinf(ay) ? sign(ay) * ax - abs(ax) - lr(ax) :
           sign(ax) == sign(ay)     ? 0.0 : -Inf # isinf(ax) && isinf(ay)
end

# cross entropy with magnetizations:
#
# -(1 + x) / 2 * log((1 + y) / 2) - (1 - x) / 2 * log((1 - y) / 2)
#
# == -x * atanh(y) - log(1 - y^2) / 2 + log(2)
#
# with atanh's:
#
# == -ay * tanh(ax) + log(2cosh(ay))
function mcrossentropy(x::MagT64, y::MagT64)
    tx = tanh(m2f(x))
    ay = m2f(y)
    return !isinf(ay)          ? -abs(ay) * (sign0(ay) * tx - 1) + lr(ay) :
           sign(tx) ≠ sign(ay) ? Inf : 0.0
end

function logZ(u0::MagT64, u::Vector{MagT64})
    a0 = m2f(u0)
    if !isinf(a0)
        s1 = a0
        s2 = abs(a0)
        s3 = lr(a0)
        hasinf = 0
    else
        s1 = s2 = s3 = 0.0
        hasinf = sign(a0)
    end
    for ui in u
        ai = m2f(ui)
        if !isinf(ai)
            s1 += ai
            s2 += abs(ai)
            s3 += lr(ai)
        elseif hasinf == 0
            hasinf = sign(ai)
        elseif hasinf ≠ sign(ai)
            return -Inf
        end
    end
    return abs(s1) - s2 + lr(s1) - s3
end


text:
using JLD
using DiffPrecTest
using DataFrames, CSV
using Random, LinearAlgebra, Distributions

pArr = [100, 200]
elemArr = [(5,5), (8, 7), (50, 25), (22, 20), (32, 30)]
methodArr = ["Sym-N", "Asym-N", "YinXia", "Sym-B", "Asym-B", "O-Sym-N", "O-Asym-N", "O-Sym-B", "O-Asym-B"]

df = DataFrame(p = Int[], row = Int[], col = Int[], trueValue=Float64[], method=String[], bias=Float64[], coverage=Float64[], lenCI=Float64[])

for ip=1:2
  for iElem=1:5
    p = pArr[ip]
    Random.seed!(7689)

    # generate model
    Ωx = Matrix{Float64}(I, p, p)
    for l=1:p-1
      Ωx[l  , l+1] = 0.3
      Ωx[l+1, l  ] = 0.3
    end
    for l=1:p-2
      Ωx[l  , l+2] = 0.2
      Ωx[l+2, l  ] = 0.2
    end
    Ωy = Matrix{Float64}(I, p, p) * 1.1
    for l=1:p-1
      Ωy[l  , l+1] = 0.3
      Ωy[l+1, l  ] = 0.3
    end
    for l=1:p-2
      Ωy[l  , l+2] = -0.1
      Ωy[l+2, l  ] = -0.1
    end
    Δ = Ωx - Ωy
    # generate Delta
    d = Vector{Float64}(undef, p)
    rand!(Uniform(0.5, 2.5), d)
    d .= sqrt.(d)
    D = Diagonal(d)
    Σx = inv(Symmetric(D * Ωx * D))
    Σy = inv(Symmetric(D * Ωy * D))

    tΔ = D * Δ * D

    res = load("../sim7_res_$(ip)_$(iElem).jld", "results")

    ri, ci = elemArr[iElem]
    indE = (ci - 1) * p + ri

    for j in 1:9
      sr = computeSimulationResult([res[j, i] for i=1:1000], tΔ[indE])
      push!(df, [pArr[ip], elemArr[iElem][1], elemArr[iElem][2], tΔ[indE], methodArr[j], sr.bias, sr.coverage, sr.lenCoverage])
    end
  end
end

CSV.write("sim7.csv", df)


text:
import DifferentiablePhysics2D
import Test

Test.@testset "DifferentiablePhysics2D.jl" begin
end


text:
using MySQLClient
using Test

@testset "MySQLClient.jl" begin
    # Write your tests here.
end


text:
using ImageTracking
using Base.Test

# write your own tests here
# @test 1 == 2

include("optical_flow.jl")
include("haar.jl")


text:
using RunTests
using Base.Test

@testmodule SixthSetOfTests begin
  function test_hello()
    @test true
  end
end


text:
using Misfits
using Test
using ForwardDiff
using BenchmarkTools
using LinearAlgebra
using Calculus
using LinearMaps

# =================================================
# generalized least-squares
# =================================================
x=abs2.(randn(10,10))
Q=x*x'
pa=Misfits.P_gls(LinearMap(Q))
x=randn(10);
y=randn(10);
dfdx1=similar(x);
@btime Misfits.func_grad!(dfdx1,x,y, pa)
@btime Misfits.func_grad!(nothing,x,y, pa)
f(x)=Misfits.func_grad!(nothing,x,y, pa)
dfdx2=Calculus.gradient(f,x)

@test dfdx1 ≈ dfdx2

# =================================================
# squared_euclidean!
# =================================================
x=randn(1000,10);
y=randn(1000,10);
w=randn(1000,10);
dfdx1=similar(x);
dfdx2=similar(x);

@btime Misfits.error_squared_euclidean!(dfdx1,x,y,w)
@btime Misfits.error_squared_euclidean!(nothing,x,y,w)
f(x)=Misfits.error_squared_euclidean!(nothing,x,y,w)
ForwardDiff.gradient!(dfdx2,f, x);

@test dfdx1 ≈ dfdx2

# =================================================
# error_after_scaling
# =================================================
# real 1D
x = randn(10,1); α = randn(1); y = α[1] .* x;
J, α1 = Misfits.error_after_scaling(x,y)
@test isapprox(α1.*x, y)
#@test_approx_eq J 0.0

# real 2D
x = randn(10,10); α = randn(1); y = α[1] .* x;
J, α1 = Misfits.error_after_scaling(x,y)
@test isapprox(α1.*x, y)

# complex 2D
x = complex.(randn(10,10), randn(10,10));
α = complex.(randn(1), randn(1)); y = α[1] .* x;
J, α1 = Misfits.error_after_scaling(x,y)
@test isapprox(α1.*x, y)

# =================================================
# weighted norm
# =================================================
x=randn(100,10);
w=randn(100,10);
dfdx1=similar(x);
@btime Misfits.error_weighted_norm!(dfdx1,x,w)
xvec=vec(x)
dfdx2=similar(xvec);
#ForwardDiff.gradient!(dfdx2,x -> Misfits.error_weighted_norm!(nothing, reshape(x,100,10), w), xvec);
dfdx2=Calculus.gradient(x -> Misfits.error_weighted_norm!(nothing, reshape(x,100,10), w), xvec)

@test vec(dfdx1) ≈ vec(dfdx2)

# =================================================
# test derivative_vector_magnitude
# =================================================

# some func
function f(x, z)
    y=x./norm(x)
    J=sum((y-z).^2)
    return J
end

# test derivative_vector_magnitude
function g!(g, x, z)
    xn=norm(x)
    rmul!(x, inv(xn))
    g1=similar(g)
    for i in eachindex(g1)
        g1[i]=2. * (x[i]-z[i])
    end
         rmul!(x, xn)
         nx=length(x)
         X=zeros(nx,nx)
         @time Misfits.derivative_vector_magnitude!(g,g1,x,X)
    return g
end

x=randn(10)
z=randn(10)
g1=zero(x)
f1(x)=f(x,z)
ForwardDiff.gradient!(g1,f1, x)
g2=zero(x)
@time g!(g2,x,z)
@test g1 ≈ g2

# error invariant of translation or global phase
#x=randn(100); y=randn() .* circshift(x,20)
#J, α = Misfits.error_after_autocorr_scaling(x,y)
#@test J < 1e-15

#
#x=randn(100,10);
#dfdx1=similar(x);
#@time Misfits.error_pairwise_corr_dist(dfdx1,x)
#xvec=vec(x)
#dfdx2=similar(xvec);
#Inversion.finite_difference!(x -> Misfits.error_pairwise_corr_dist(nothing, reshape(x,100,10)), xvec, dfdx2, :central)
#
#@test dfdx1 ≈ reshape(dfdx2,100,10)
#
#
#x=randn(10,10);
#dfdx1=similar(x);
#Misfits.error_autocorr_pairwise_corr_dist(dfdx1,x)
#xvec=vec(x)
#dfdx2=similar(xvec);
#Inversion.finite_difference!(x -> Misfits.error_autocorr_pairwise_corr_dist(nothing, reshape(x,10,10)), xvec, dfdx2, :central)
#
#@test dfdx1 ≈ reshape(dfdx2,10,10)


text:
mutable struct PushRelabelHeap{V<:Integer,T<:Integer}
    data::MutableBinaryHeap{Pair{V,T},LessThan2}
    handles::Vector{V}
end

function PushRelabelHeap{V,T}(n::Integer) where {V<:Integer,T<:Integer}
    handles = zeros(V, n)
    data = MutableBinaryHeap{Pair{V,T},LessThan2}()
    return PushRelabelHeap{V,T}(data, handles)
end

function push!(h::PushRelabelHeap{V,T}, k::V, v::T) where {V<:Integer,T<:Integer}
    a = push!(h.data, k=>v)
    h.handles[k] = a
end
push!(h::PushRelabelHeap{V,T}, k::Integer, v::Integer) where {V,T} = push!(h, V(k), T(v))

function update!(h::PushRelabelHeap{V,T}, k::V, v::T) where {V<:Integer,T<:Integer}
    a = h.handles[k]
    a > 0 && update!(h.data, a, k=>v)
end
update!(h::PushRelabelHeap{V,T}, k::Integer, v::Integer) where {V<:Integer,T<:Integer} = update!(h, V(k), T(v))

function pop!(h::PushRelabelHeap)
    k = pop!(h.data)[1]
    h.handles[k] = 0
    return k
end

length(h::PushRelabelHeap) = length(h.data)

"""
Implementation of the push relabel algorithm with gap and highest excess heuristics.
Takes O(V²√E) time.

Maintains the following auxillary arrays:
- height -> Stores the labels of all vertices
- count  -> Stores the number of vertices at each height
- excess -> Stores the difference between incoming and outgoing flow for all vertices
- active -> Stores the status of all vertices. (e(v)>0 => active[v] = true)
- Q      -> The heap that stores active vertices waiting to be discharged.

Requires arguments:

- g::ADiGraph                # the input graph
- source                            # the source vertex
- target                            # the target vertex
- capacity_matrix::AbstractMatrix{T}    # edge flow capacities
"""
function push_relabel_impl(
        g::ADiGraph,               # the input graph
        source,                           # the source vertex
        target,                           # the target vertex
        capacity_matrix::Vector{Vector{T}},   # edge flow capacities
        pos
    ) where T<:Number

    n = nv(g)
    # flow_matrix = zeros(T, n, n)
    flow_matrix = deepcopy(capacity_matrix) # initialize flow matrix
    for i=1:n
        fill!(flow_matrix[i], 0)
    end
    height = zeros(Int, n)
    height[source] = n

    count = zeros(Int, 2*n+1)
    count[0+1] = n-1
    count[n+1] = 1

    excess = zeros(T, n)
    excess[source] = typemax(T)

    active = falses(n)
    active[source] = true
    active[target] = true

    Q = PushRelabelHeap{Int,Int}(n)

    for (k, u) in enumerate(neighbors(g, source))
        push_flow!(g, source, u, k, capacity_matrix, flow_matrix, excess, height, active, Q, pos)
    end

    while length(Q) > 0
        v = pop!(Q)
        active[v] = false
        discharge!(g, v, capacity_matrix, flow_matrix, excess, height, active, count, Q, pos)
    end

    flow = zero(T)
    for (k, j) in enumerate(neighbors(g, target))
        k2 = pos[target][k]
        flow += flow_matrix[j][k2]
    end
    # return sum([flow_matrix[v,target] for v in badj(g, target) ]), flow_matrix
    return flow, flow_matrix
end

"""
Pushes inactive nodes into the queue and activates them.

Requires arguments:

- Q::PushRelabelHeap
- v
- active::AbstractVector{Bool}
- excess::AbstractVector{T}
"""
function enqueue_vertex!(
        Q::PushRelabelHeap,
        v,                                # input vertex
        active::AbstractVector{Bool},
        excess::Vector{T},
        height::Vector{Int}) where T<:Number

    @inbounds if !active[v] && excess[v] > 0
        active[v] = true
        push!(Q, v, height[v])
    end
    return nothing
end

"""
Pushes as much flow as possible through the given edge.

Requires arguements:

- g::ADiGraph              # the input graph
- u                               # input from-vertex
- v                               # input to-vetex
- capacity_matrix::AbstractMatrix{T}
- flow_matrix::AbstractMatrix{T}
- excess::AbstractVector{T}
- height::AbstractVector{Int}
- active::AbstractVector{Bool}
- Q::PushRelabelHeap
"""
function push_flow!(
        g::ADiGraph,             # the input graph
        u,                              # input from-vertex
        v,                                 # input to-vetex
        k,                                 #index of v as neig of u
        capacity_matrix,
        flow_matrix,
        excess::Vector{T},
        height::Vector{Int},
        active::AbstractVector{Bool},
        Q::PushRelabelHeap,
        pos
    ) where T<:Number

    height[u] <= height[v] && return
    flow = min(excess[u], capacity_matrix[u][k] - flow_matrix[u][k])
    flow == 0 && return

    flow_matrix[u][k] += flow
    k2 = pos[u][k]
    flow_matrix[v][k2] -= flow

    excess[u] -= flow
    excess[v] += flow

    #no need to update u since it is not in heap
    # update!(Q, v, excess[v])

    enqueue_vertex!(Q, v, active, excess, height)
end

"""
Implements the gap heuristic. Relabels all vertices above a cutoff height.
Reduces the number of relabels required.

Requires arguments:

- g::ADiGraph                # the input graph
- h                                 # cutoff height
- excess::AbstractVector{T}
- height::AbstractVector{Int}
- active::AbstractVector{Bool}
- count::AbstractVector{Int}
- Q::PushRelabelHeap
"""
function gap!(
        g::ADiGraph,               # the input graph
        h,                                # cutoff height
        excess::AbstractVector{T},
        height::AbstractVector{Int},
        active::AbstractVector{Bool},
        count::AbstractVector{Int},
        Q::PushRelabelHeap
    ) where T<:Number

    n = nv(g)
    @inbounds for v in vertices(g)
        hv = height[v]
        hv < h && continue
        count[hv+1] -= 1
        hv = max(hv, n + 1)
        count[hv+1] += 1
        height[v] = hv
        update!(Q, v, height[v])
        enqueue_vertex!(Q, v, active, excess, height)
    end
end

"""
Relabels a vertex with respect to its neighbors, to produce an admissable
edge.

Requires arguments:

- g::ADiGraph                 # the input graph
- v                                  # input vertex to be relabeled
- capacity_matrix::AbstractMatrix{T}
- flow_matrix::AbstractMatrix{T}
- excess::AbstractVector{T}
- height::AbstractVector{Int}
- active::AbstractVector{Bool}
- count::AbstractVector{Int}
- Q::AbstractVector{Int}
"""
function relabel!(
        g::ADiGraph,                # the input graph
        v,                                 # input vertex to be relabeled
        capacity_matrix,
        flow_matrix,
        excess::Vector{T},
        height::Vector{Int},
        active::AbstractVector{Bool},
        count::Vector{Int},
        Q::PushRelabelHeap
    ) where T<:Number

    n = nv(g)
    count[height[v]+1] -= 1
    height[v] = 2*n
    @inbounds for (k,to) in enumerate(neighbors(g, v))
        if capacity_matrix[v][k] > flow_matrix[v][k]
            height[v] = min(height[v], height[to]+1)
        end
    end
    count[height[v]+1] += 1
    update!(Q, v, height[v]) #TODO remove since v is not in Q
    enqueue_vertex!(Q, v, active, excess, height)
end

"""
Drains the excess flow out of a vertex. Runs the gap heuristic or relabels the
vertex if the excess remains non-zero.

Requires arguments:

- g::ADiGraph                 # the input graph
- v                                  # vertex to be discharged
- capacity_matrix::AbstractMatrix{T}
- flow_matrix::AbstractMatrix{T}
- excess::AbstractVector{T}
- height::AbstractVector{Int}
- active::AbstractVector{Bool}
- count::AbstractVector{Int}
- Q::PushRelabelHeap
"""
function discharge!(
        g::ADiGraph,                # the input graph
        v,                                 # vertex to be discharged
        capacity_matrix,
        flow_matrix,
        excess::Vector{T},
        height::Vector{Int},
        active::AbstractVector{Bool},
        count::Vector{Int},
        Q::PushRelabelHeap,
        pos
    ) where T<:Number

    for (k, u) in enumerate(neighbors(g, v))
        excess[v] == 0 && break
        push_flow!(g, v, u, k, capacity_matrix, flow_matrix, excess, height, active, Q, pos)
    end

    if excess[v] > 0
        if count[height[v]+1] == 1
            gap!(g, height[v], excess, height, active, count, Q)
        else
            relabel!(g, v, capacity_matrix, flow_matrix, excess, height, active, count, Q)
        end
    end
end


text:
using Cards
using StatsBase

mutable struct Deck
    deck::Hand
end

suits = [♣, ♢, ♡, ♠]
deck = [Card(r,s) for s in suits for r = 2:14]

function suit(c::Card)
    Suit((0x30 & c.value) >>> 4)
end

function rank(c::Card)
    (c.value & 0x0f) % Int8
end

function blackjackSum(hand::Array{Card})
    ranks = [rank(card) for card ∈ hand]
    numAces = length(findall(x -> x ∈ [Card(14,i) for i in 1:4], hand))
    # Returns sum with high aces and low aces
    out = sum(x -> min(10, x), ranks) + numAces, sum(x -> min(10, x), ranks) - 9*numAces
    return out
end

function shuffleDeck(nDecks::Integer)
    return [Card(r,s) for s in suits for r = 2:14 for _ in 1:nDecks]
end

function Deal(deck::Array{Card,1}, cardsPerHand::Integer, numHands::Integer)
    out = Array{Card}(undef, numHands, cardsPerHand)
    for i ∈ 1:numHands
        for j ∈ 1:cardsPerHand
            out[i,j] = rand(deck)
            deleteat!(deck, findall(x->x == out[i,j], deck))
        end
    end
    println(length(deck))
    return out
end

deck = shuffleDeck(1)
hand = Deal(deck, 2, 1)
blackjackSum(hand)

max_stars_repo_name: blazickjp/ReinforcementLearning

max_issues_repo_name: blazickjp/ReinforcementLearning

max_forks_repo_name: blazickjp/ReinforcementLearning


text:
"""
    LaTeXEngine

Subtypes of this type represent LaTeX engines and should be singletons.
"""
abstract type LaTeXEngine end

Base.show(io::IO, ::T) where T <: LaTeXEngine = print(io, nameof(T), " engine")

const _CURRENT_ENGINE = Ref{LaTeXEngine}()

"""
    texengine()

Returns the current LaTeX engine.
"""
texengine() = _CURRENT_ENGINE[]

"""
    texengine!(eng)

Sets the LaTeX engine for this session. `eng` can be [`PDFLaTeX`](@ref), [`XeLaTeX`](@ref), or [`LuaLaTeX`](@ref), e.g.
```julia
texengine!(PDFLaTeX)
```
"""
function texengine!(eng::LaTeXEngine)
    _CURRENT_ENGINE[] = eng
    return nothing
end

texengine!(::Type{T}) where T <: LaTeXEngine = texengine!(T())

"""
    runlatex(<:LaTeXEngine, input_file::AbstractString, output_path::AbstractString; extra_args=String[])

Compiles `input_file` into a dvi document of the same file name in `output_path`. The first argument needs to be an instance of a subtype of `LaTeXEngine` and determines the associated LaTeX command to run. The keyword argument `extra_args` is a `Vector` of `String`s that contains additional arguments/flags to pass to the LaTeX engine, and defaults to being empty. Make sure that strings in `extra_args` contain no spaces, or are otherwise valid shell arguments.
"""
function runlatex(
    engine::LaTeXEngine,
    input_file::AbstractString,
    output_path::AbstractString;
    extra_args::Vector{String}=String[],
)
    error("You need to overload the `runlatex` function for `$(nameof(typeof(engine)))`.")
end

"""
    dviext(::LaTeXEngine)

Returns the file extension of the dvi output of the [`LaTeXEngine`](@ref).
"""
function dviext(engine::LaTeXEngine)
    error("You need to overload the `dviext` function for `$(nameof(typeof(engine)))`.")
end

function _tex2dvi(filename::AbstractString, engine::LaTeXEngine; extra_args::Vector{String}=String[])
    output_path = joinpath(splitpath(filename)[begin:end - 1]...)
    runlatex(engine, filename * ".tex", output_path; extra_args=extra_args)
    return filename * "." * dviext(engine)
end


text:
export BCQLearner

mutable struct BCQLearner{
    BA1<:NeuralNetworkApproximator,
    BA2<:NeuralNetworkApproximator,
    BC1<:NeuralNetworkApproximator,
    BC2<:NeuralNetworkApproximator,
    V<:NeuralNetworkApproximator,
    R<:AbstractRNG,
} <: Any
    policy::BA1
    target_policy::BA2
    qnetwork1::BC1
    qnetwork2::BC2
    target_qnetwork1::BC1
    target_qnetwork2::BC2
    vae::V
    γ::Float32
    τ::Float32
    λ::Float32
    p::Int
    batch_size::Int
    start_step::Int
    update_freq::Int
    update_step::Int
    rng::R
    # Logging
    actor_loss::Float32
    critic_loss::Float32
end

"""
    BCQLearner(;kwargs...)

See [Off-Policy Deep Reinforcement Learning without Exploration](https://arxiv.org/abs/1812.02900)

# Keyword arguments
- `policy`, used to get action with perturbation. This can be implemented using a `PerturbationNetwork` in a `NeuralNetworkApproximator`.
- `target_policy`, similar to `policy`, but used to estimate the target. This can be implemented using a `PerturbationNetwork` in a `NeuralNetworkApproximator`.
- `qnetwork1`, used to get Q-values.
- `qnetwork2`, used to get Q-values.
- `target_qnetwork1`, used to estimate the target Q-values.
- `target_qnetwork2`, used to estimate the target Q-values.
- `vae`, used for sampling action. This
can be implemented using a `VAE` in a `NeuralNetworkApproximator`.
- `γ::Float32 = 0.99f0`, reward discount rate.
- `τ::Float32 = 0.005f0`, the speed at which the target network is updated.
- `λ::Float32 = 0.75f0`, used for Clipped Double Q-learning.
- `p::Int = 10`, the number of state-action pairs used when calculating the Q value.
- `batch_size::Int = 32`
- `start_step::Int = 1000`
- `update_freq::Int = 50`, the frequency of updating the `approximator`.
- `update_step::Int = 0`
- `rng = Random.GLOBAL_RNG`
"""
function BCQLearner(;
    policy,
    target_policy,
    qnetwork1,
    qnetwork2,
    target_qnetwork1,
    target_qnetwork2,
    vae,
    γ = 0.99f0,
    τ = 0.005f0,
    λ = 0.75f0,
    p = 10,
    batch_size = 32,
    start_step = 1000,
    update_freq = 50,
    update_step = 0,
    rng = Random.GLOBAL_RNG,
)
    copyto!(policy, target_policy)  # force sync
    copyto!(qnetwork1, target_qnetwork1)  # force sync
    copyto!(qnetwork2, target_qnetwork2)  # force sync
    BCQLearner(
        policy,
        target_policy,
        qnetwork1,
        qnetwork2,
        target_qnetwork1,
        target_qnetwork2,
        vae,
        γ,
        τ,
        λ,
        p,
        batch_size,
        start_step,
        update_freq,
        update_step,
        rng,
        0.0f0,
        0.0f0,
    )
end

function (l::BCQLearner)(env)
    s = send_to_device(device(l.policy), state(env))
    s = Flux.unsqueeze(s, ndims(s) + 1)
    s = repeat(s, outer = (1, 1, l.p))
    action = l.policy(s, decode(l.vae.model, s))
    q_value = l.qnetwork1(vcat(s, action))
    idx = argmax(q_value)
    action[idx]
end

function RLBase.update!(l::BCQLearner, batch::NamedTuple{SARTS})
    update_vae!(l, batch)
    if l.update_step >= l.start_step
        update_learner!(l, batch)
    end
end

function update_vae!(l::BCQLearner, batch::NamedTuple{SARTS})
    D = device(l.vae)
    s, a, r, t, s′ = (send_to_device(D, batch[x]) for x in SARTS)
    a = reshape(a, :, l.batch_size)
    vae_grad = gradient(Flux.params(l.vae)) do
        recon_loss, kl_loss = vae_loss(l.vae.model, s, a)
        0.5f0 * kl_loss + recon_loss
    end
    update!(l.vae, vae_grad)
end

function update_learner!(l::BCQLearner, batch::NamedTuple{SARTS})
    D = device(l.qnetwork1)
    s, a, r, t, s′ = (send_to_device(D, batch[x]) for x in SARTS)

    γ, τ, λ = l.γ, l.τ, l.λ

    repeat_s′ = repeat(s′, outer = (1, 1, l.p))
    repeat_a′ = l.target_policy(repeat_s′, decode(l.vae.model, repeat_s′))

    q′_input = vcat(repeat_s′, repeat_a′)
    q′ = maximum(
        λ .* min.(l.target_qnetwork1(q′_input), l.target_qnetwork2(q′_input)) +
        (1 - λ) .* max.(l.target_qnetwork1(q′_input), l.target_qnetwork2(q′_input)),
        dims = 3,
    )

    y = r .+ γ .* (1 .- t) .* vec(q′)

    # Train Q Networks
    a = reshape(a, :, l.batch_size)
    q_input = vcat(s, a)

    q_grad_1 = gradient(Flux.params(l.qnetwork1)) do
        q1 = l.qnetwork1(q_input) |> vec
        loss = mse(q1, y)
        ignore() do
            l.critic_loss = loss
        end
        loss
    end
    update!(l.qnetwork1, q_grad_1)

    q_grad_2 = gradient(Flux.params(l.qnetwork2)) do
        q2 = l.qnetwork2(q_input) |> vec
        loss = mse(q2, y)
        ignore() do
            l.critic_loss += loss
        end
        loss
    end
    update!(l.qnetwork2, q_grad_2)

    # Train Policy
    p_grad = gradient(Flux.params(l.policy)) do
        sampled_action = decode(l.vae.model, s)
        perturbed_action = l.policy(s, sampled_action)
        actor_loss = -mean(l.qnetwork1(vcat(s, perturbed_action)))
        ignore() do
            l.actor_loss = actor_loss
        end
        actor_loss
    end
    update!(l.policy, p_grad)

    # polyak averaging
    for (dest, src) in zip(
        Flux.params([l.target_policy, l.target_qnetwork1, l.target_qnetwork2]),
        Flux.params([l.policy, l.qnetwork1, l.qnetwork2]),
    )
        dest .= (1 - τ) .* dest .+ τ .* src
    end
end


text:
# Necessary to have initialize set u_modified to false if all don't do anything
# otherwise unnecessary save
INITIALIZE_DEFAULT(cb,u,t,integrator) = u_modified!(integrator, false)
FINALIZE_DEFAULT(cb,u,t,integrator) = nothing

@enum RootfindOpt::Int8 begin
  NoRootFind    = 0
  LeftRootFind  = 1
  RightRootFind = 2
end

function Base.convert(::Type{RootfindOpt}, b::Bool)
  return b ? LeftRootFind : NoRootFind
end

"""
```julia
ContinuousCallback(condition,affect!,affect_neg!;
                   initialize = INITIALIZE_DEFAULT,
                   finalize = FINALIZE_DEFAULT,
                   idxs = nothing,
                   rootfind=LeftRootFind,
                   save_positions=(true,true),
                   interp_points=10,
                   abstol=10eps(),reltol=0,repeat_nudge=1//100)
```

```julia
function ContinuousCallback(condition,affect!;
                   initialize = INITIALIZE_DEFAULT,
                   finalize = FINALIZE_DEFAULT,
                   idxs = nothing,
                   rootfind=LeftRootFind,
                   save_positions=(true,true),
                   affect_neg! = affect!,
                   interp_points=10,
                   abstol=10eps(),reltol=0,repeat_nudge=1//100)
```

Contains a single callback whose `condition` is a continuous function. The callback is triggered when this function evaluates to 0.

# Arguments
- `condition`: This is a function `condition(u,t,integrator)` for declaring when
  the callback should be used. A callback is initiated if the condition hits
  `0` within the time interval. See the [Integrator Interface](@ref integrator) documentation for information about `integrator`.
- `affect!`: This is the function `affect!(integrator)` where one is allowed to
  modify the current state of the integrator. If you do not pass an `affect_neg!`
  function, it is called when `condition` is found to be `0` (at a root) and
  the cross is either an upcrossing (from negative to positive) or a downcrossing
  (from positive to negative). You need to explicitly pass `nothing` as the
  `affect_neg!` argument if it should only be called at upcrossings, e.g.
  `ContinuousCallback(condition, affect!, nothing)`. For more information on what can
  be done, see the [Integrator Interface](@ref integrator) manual page. Modifications to
  `u` are safe in this function.
- `affect_neg!=affect!`: This is the function `affect_neg!(integrator)` where one is allowed to
  modify the current state of the integrator. This is called when `condition` is
  found to be `0` (at a root) and the cross is an downcrossing (from positive to
  negative). For more information on what can
  be done, see the [Integrator Interface](@ref integrator) manual page. Modifications to
  `u` are safe in this function.
- `rootfind=LeftRootFind`: This is a flag to specify the type of rootfinding to do for finding
  event location. If this is set to `LeftRootfind`, the solution will be backtracked to the point where
  `condition==0` and if the solution isn't exact, the left limit of root is used. If set to
  `RightRootFind`, the solution would be set to the right limit of the root. Otherwise the systems and
  the `affect!` will occur at `t+dt`.
- `interp_points=10`: The number of interpolated points to check the condition. The
  condition is found by checking whether any interpolation point / endpoint has
  a different sign. If `interp_points=0`, then conditions will only be noticed if
  the sign of `condition` is different at `t` than at `t+dt`. This behavior is not
  robust when the solution is oscillatory, and thus it's recommended that one use
  some interpolation points (they're cheap to compute!).
  `0` within the time interval.
- `save_positions=(true,true)`: Boolean tuple for whether to save before and after the `affect!`.
  This saving will occur just before and after the event, only at event times, and
  does not depend on options like `saveat`, `save_everystep`, etc. (i.e. if
  `saveat=[1.0,2.0,3.0]`, this can still add a save point at `2.1` if true).
  For discontinuous changes like a modification to `u` to be
  handled correctly (without error), one should set `save_positions=(true,true)`.
- `idxs=nothing`: The components which will be interpolated into the condition. Defaults
  to `nothing` which means `u` will be all components.
- `initialize`: This is a function `(c,u,t,integrator)` which can be used to initialize
  the state of the callback `c`. It should modify the argument `c` and the return is
  ignored.
- `finalize`: This is a function `(c,u,t,integrator)` which can be used to finalize
  the state of the callback `c`. It can modify the argument `c` and the return is ignored.
- `abstol=1e-14` & `reltol=0`: These are used to specify a tolerance from zero for the rootfinder:
  if the starting condition is less than the tolerance from zero, then no root will be detected.
  This is to stop repeat events happening just after a previously rootfound event.
- `repeat_nudge = 1//100`: This is used to set the next testing point after a
  previously found zero. Defaults to 1//100, which means after a callback the next
  sign check will take place at t + dt*1//100 instead of at t to avoid repeats.
"""
struct ContinuousCallback{F1,F2,F3,F4,F5,T,T2,T3,I,R} <: AbstractContinuousCallback
  condition::F1
  affect!::F2
  affect_neg!::F3
  initialize::F4
  finalize::F5
  idxs::I
  rootfind::RootfindOpt
  interp_points::Int
  save_positions::BitArray{1}
  dtrelax::R
  abstol::T
  reltol::T2
  repeat_nudge::T3
  ContinuousCallback(condition::F1,affect!::F2,affect_neg!::F3,
                     initialize::F4,finalize::F5,idxs::I,rootfind,
                     interp_points,save_positions,dtrelax::R,abstol::T,reltol::T2,
                     repeat_nudge::T3) where {F1,F2,F3,F4,F5,T,T2,T3,I,R} =
                       new{F1,F2,F3,F4,F5,T,T2,T3,I,R}(condition,
                                               affect!,affect_neg!,
                                               initialize,finalize,idxs,rootfind,interp_points,
                                               BitArray(collect(save_positions)),
                                               dtrelax,abstol,reltol,repeat_nudge)
end

ContinuousCallback(condition,affect!,affect_neg!;
                   initialize = INITIALIZE_DEFAULT,
                   finalize = FINALIZE_DEFAULT,
                   idxs = nothing,
                   rootfind=LeftRootFind,
                   save_positions=(true,true),
                   interp_points=10,
                   dtrelax=1,
                   abstol=10eps(),reltol=0,
                   repeat_nudge = 1//100) = ContinuousCallback(
                              condition,affect!,affect_neg!,initialize,finalize,
                              idxs,
                              rootfind,interp_points,
                              save_positions,
                              dtrelax,abstol,reltol,repeat_nudge)

function ContinuousCallback(condition,affect!;
                   initialize = INITIALIZE_DEFAULT,
                   finalize = FINALIZE_DEFAULT,
                   idxs = nothing,
                   rootfind=LeftRootFind,
                   save_positions=(true,true),
                   affect_neg! = affect!,
                   interp_points=10,
                   dtrelax=1,
                   abstol=10eps(),reltol=0,repeat_nudge=1//100)

 ContinuousCallback(
            condition,affect!,affect_neg!,initialize,finalize,idxs,
            rootfind,interp_points,
            collect(save_positions),
            dtrelax,abstol,reltol,repeat_nudge)

end

"""
```julia
VectorContinuousCallback(condition,affect!,affect_neg!,len;
                         initialize = INITIALIZE_DEFAULT,
                         finalize = FINALIZE_DEFAULT,
                         idxs = nothing,
                         rootfind=LeftRootFind,
                         save_positions=(true,true),
                         interp_points=10,
                         abstol=10eps(),reltol=0,repeat_nudge = 1//100)
```

```julia
VectorContinuousCallback(condition,affect!,len;
                   initialize = INITIALIZE_DEFAULT,
                   finalize = FINALIZE_DEFAULT,
                   idxs = nothing,
                   rootfind=LeftRootFind,
                   save_positions=(true,true),
                   affect_neg! = affect!,
                   interp_points=10,
                   abstol=10eps(),reltol=0,repeat_nudge=1//100)
```

This is also a subtype of `AbstractContinuousCallback`. `CallbackSet` is not feasible when you have a large number of callbacks,
as it doesn't scale well. For this reason, we have `VectorContinuousCallback` - it allows you to have a single callback for
multiple events.

# Arguments

- `condition`: This is a function `condition(out, u, t, integrator)` which should save the condition value in the array `out`
   at the right index. Maximum index of `out` should be specified in the `len` property of callback. So this way you can have
   a chain of `len` events, which would cause the `i`th event to trigger when `out[i] = 0`.
- `affect!`: This is a function `affect!(integrator, event_index)` which lets you modify `integrator` and it tells you about
   which event occured using `event_idx` i.e. gives you index `i` for which `out[i]` came out to be zero.
- `len`: Number of callbacks chained. This is compulsory to be specified.

Rest of the arguments have the same meaning as in [`ContinuousCallback`](@ref).
"""
struct VectorContinuousCallback{F1,F2,F3,F4,F5,T,T2,T3,I,R} <: AbstractContinuousCallback
  condition::F1
  affect!::F2
  affect_neg!::F3
  len::Int
  initialize::F4
  finalize::F5
  idxs::I
  rootfind::RootfindOpt
  interp_points::Int
  save_positions::BitArray{1}
  dtrelax::R
  abstol::T
  reltol::T2
  repeat_nudge::T3
  VectorContinuousCallback(condition::F1,affect!::F2,affect_neg!::F3,len::Int,
                           initialize::F4,finalize::F5,idxs::I,rootfind,
                           interp_points,save_positions,dtrelax::R,
                           abstol::T,reltol::T2,repeat_nudge::T3) where {F1,F2,F3,F4,F5,T,T2,T3,I,R} =
                       new{F1,F2,F3,F4,F5,T,T2,T3,I,R}(condition,
                                               affect!,affect_neg!,len,
                                               initialize,finalize,idxs,rootfind,interp_points,
                                               BitArray(collect(save_positions)),
                                               dtrelax,abstol,reltol,repeat_nudge)
end

VectorContinuousCallback(condition,affect!,affect_neg!,len;
                         initialize = INITIALIZE_DEFAULT,
                         finalize = FINALIZE_DEFAULT,
                         idxs = nothing,
                         rootfind=LeftRootFind,
                         save_positions=(true,true),
                         interp_points=10,
                         dtrelax=1,
                         abstol=10eps(),reltol=0,repeat_nudge=1//100) = VectorContinuousCallback(
                              condition,affect!,affect_neg!,len,
                              initialize,finalize,
                              idxs,
                              rootfind,interp_points,
                              save_positions,dtrelax,
                              abstol,reltol,repeat_nudge)

function VectorContinuousCallback(condition,affect!,len;
                   initialize = INITIALIZE_DEFAULT,
                   finalize = FINALIZE_DEFAULT,
                   idxs = nothing,
                   rootfind=LeftRootFind,
                   save_positions=(true,true),
                   affect_neg! = affect!,
                   interp_points=10,
                   dtrelax=1,
                   abstol=10eps(),reltol=0,repeat_nudge=1//100)

 VectorContinuousCallback(
            condition,affect!,affect_neg!,len,initialize,finalize,idxs,
            rootfind,interp_points,
            collect(save_positions),
            dtrelax,abstol,reltol,repeat_nudge)

end

"""
```julia
DiscreteCallback(condition,affect!;
                 initialize = INITIALIZE_DEFAULT,
                 finalize = FINALIZE_DEFAULT,
                 save_positions=(true,true))
```

# Arguments

- `condition`: This is a function `condition(u,t,integrator)` for declaring when
  the callback should be used. A callback is initiated if the condition evaluates
  to `true`. See the [Integrator Interface](@ref integrator) documentation for information about `integrator`.
    - `affect!`: This is the function `affect!(integrator)` where one is allowed to
  modify the current state of the integrator. For more information on what can
  be done, see the [Integrator Interface](@ref integrator) manual page.
- `save_positions`: Boolean tuple for whether to save before and after the `affect!`.
  This saving will occur just before and after the event, only at event times, and
  does not depend on options like `saveat`, `save_everystep`, etc. (i.e. if
  `saveat=[1.0,2.0,3.0]`, this can still add a save point at `2.1` if true).
  For discontinuous changes like a modification to `u` to be
  handled correctly (without error), one should set `save_positions=(true,true)`.
- `initialize`: This is a function `(c,u,t,integrator)` which can be used to initialize
  the state of the callback `c`. It should modify the argument `c` and the return is
  ignored.
- `finalize`: This is a function `(c,u,t,integrator)` which can be used to finalize
  the state of the callback `c`. It should can the argument `c` and the return is
  ignored.
"""
struct DiscreteCallback{F1,F2,F3,F4} <: AbstractDiscreteCallback
  condition::F1
  affect!::F2
  initialize::F3
  finalize::F4
  save_positions::BitArray{1}
  DiscreteCallback(condition::F1,affect!::F2,
                   initialize::F3,finalize::F4,save_positions) where {F1,F2,F3,F4} = new{F1,F2,F3,F4}(condition,
                                                                                affect!,initialize,finalize,
                                                                                BitArray(collect(save_positions)))
end
DiscreteCallback(condition,affect!;
                 initialize = INITIALIZE_DEFAULT, finalize = FINALIZE_DEFAULT,
                 save_positions=(true,true)) = DiscreteCallback(condition,affect!,initialize,finalize,save_positions)

"""
$(TYPEDEF)

Multiple callbacks can be chained together to form a `CallbackSet`. A `CallbackSet`
is constructed by passing the constructor `ContinuousCallback`, `DiscreteCallback`,
`VectorContinuousCallback` or other `CallbackSet` instances:

    CallbackSet(cb1,cb2,cb3)

You can pass as many callbacks as you like. When the solvers encounter multiple
callbacks, the following rules apply:

* `ContinuousCallback`s and `VectorContinuousCallback`s are applied before `DiscreteCallback`s. (This is because
  they often implement event-finding that will backtrack the timestep to smaller
  than `dt`).
* For `ContinuousCallback`s and `VectorContinuousCallback`s, the event times are found by rootfinding and only
  the first `ContinuousCallback` or `VectorContinuousCallback` affect is applied.
* The `DiscreteCallback`s are then applied in order. Note that the ordering only
  matters for the conditions: if a previous callback modifies `u` in such a way
  that the next callback no longer evaluates condition to `true`, its `affect`
  will not be applied.

"""
struct CallbackSet{T1<:Tuple,T2<:Tuple} <: DECallback
  continuous_callbacks::T1
  discrete_callbacks::T2
end

CallbackSet(callback::AbstractDiscreteCallback) = CallbackSet((),(callback,))
CallbackSet(callback::AbstractContinuousCallback) = CallbackSet((callback,),())
CallbackSet() = CallbackSet((),())
CallbackSet(cb::Nothing) = CallbackSet()

# For Varargs, use recursion to make it type-stable
CallbackSet(callbacks::Union{DECallback,Nothing}...) = CallbackSet(split_callbacks((), (), callbacks...)...)

"""
    split_callbacks(cs, ds, args...)

Split comma seperated callbacks into sets of continous and discrete callbacks.
"""
@inline split_callbacks(cs, ds) = cs, ds
@inline split_callbacks(cs, ds, c::Nothing, args...) = split_callbacks(cs, ds, args...)
@inline split_callbacks(cs, ds, c::AbstractContinuousCallback, args...) = split_callbacks((cs..., c), ds, args...)
@inline split_callbacks(cs, ds, d::AbstractDiscreteCallback, args...) = split_callbacks(cs, (ds..., d), args...)
@inline function split_callbacks(cs, ds, d::CallbackSet, args...)
  split_callbacks((cs...,d.continuous_callbacks...), (ds..., d.discrete_callbacks...), args...)
end

"""
    initialize!(cb::CallbackSet,u,t,integrator::DEIntegrator)

Recursively apply `initialize!` and return whether any modified u
"""
function initialize!(cb::CallbackSet,u,t,integrator::DEIntegrator)
  initialize!(u,t,integrator,false,cb.continuous_callbacks...,cb.discrete_callbacks...)
end
initialize!(cb::CallbackSet{Tuple{},Tuple{}},u,t,integrator::DEIntegrator) = false
function initialize!(u,t,integrator::DEIntegrator,any_modified::Bool,
                     c::DECallback,cs::DECallback...)
  c.initialize(c,u,t,integrator)
  initialize!(u,t,integrator,any_modified || integrator.u_modified,cs...)
end
function initialize!(u,t,integrator::DEIntegrator,any_modified::Bool,
                     c::DECallback)
  c.initialize(c,u,t,integrator)
  any_modified || integrator.u_modified
end

"""
    finalize!(cb::CallbackSet,u,t,integrator::DEIntegrator)

Recursively apply `finalize!` and return whether any modified u
"""
function finalize!(cb::CallbackSet,u,t,integrator::DEIntegrator)
  finalize!(u,t,integrator,false,cb.continuous_callbacks...,cb.discrete_callbacks...)
end
finalize!(cb::CallbackSet{Tuple{},Tuple{}},u,t,integrator::DEIntegrator) = false
function finalize!(u,t,integrator::DEIntegrator,any_modified::Bool,
                     c::DECallback,cs::DECallback...)
  c.finalize(c,u,t,integrator)
  finalize!(u,t,integrator,any_modified || integrator.u_modified,cs...)
end
function finalize!(u,t,integrator::DEIntegrator,any_modified::Bool,
                     c::DECallback)
  c.finalize(c,u,t,integrator)
  any_modified || integrator.u_modified
end

# Helpers
Base.isempty(cb::CallbackSet) = isempty(cb.continuous_callbacks) && isempty(cb.discrete_callbacks)
Base.isempty(cb::AbstractContinuousCallback) = false
Base.isempty(cb::AbstractDiscreteCallback) = false

has_continuous_callback(cb::DiscreteCallback) = false
has_continuous_callback(cb::ContinuousCallback) = true
has_continuous_callback(cb::VectorContinuousCallback) = true
has_continuous_callback(cb::CallbackSet) = !isempty(cb.continuous_callbacks)
has_continuous_callback(cb::Nothing) = false

#======================================================#
# Callback handling
#======================================================#

function get_tmp(integrator::DEIntegrator, callback)
  _tmp = get_tmp_cache(integrator)
  _tmp === nothing && return nothing
  _cache = first(_tmp)
  if callback.idxs === nothing
    tmp = _cache
  elseif !(callback.idxs isa Number)
    tmp = @view _cache[callback.idxs]
  else
    tmp = nothing
  end
  return tmp
end

function get_condition(integrator::DEIntegrator, callback, abst)
  tmp = get_tmp(integrator, callback)
  ismutable = !(tmp === nothing)
  if abst == integrator.t
    if callback.idxs === nothing
      tmp = integrator.u
    elseif callback.idxs isa Number
      tmp = integrator.u[callback.idxs]
    else
      tmp = @view integrator.u[callback.idxs]
    end
  elseif abst == integrator.tprev
    if callback.idxs === nothing
      tmp = integrator.uprev
    elseif callback.idxs isa Number
      tmp = integrator.uprev[callback.idxs]
    else
      tmp = @view integrator.uprev[callback.idxs]
    end
  else
    if ismutable
      if callback.idxs === nothing
        integrator(tmp,abst,Val{0})
      else
        integrator(tmp,abst,Val{0},idxs=callback.idxs)
      end
    else
      if callback.idxs === nothing
        tmp = integrator(abst,Val{0})
      else
        tmp = integrator(abst,Val{0},idxs=callback.idxs)
      end
    end
    # ismutable && !(callback.idxs isa Number) ? integrator(tmp,abst,Val{0},idxs=callback.idxs) :
    #                                                 tmp = integrator(abst,Val{0},idxs=callback.idxs)
  end
  integrator.sol.destats.ncondition += 1
  if callback isa VectorContinuousCallback
    callback.condition(@view(integrator.callback_cache.tmp_condition[1:callback.len]),tmp,abst,integrator)
    return @view(integrator.callback_cache.tmp_condition[1:callback.len])
  else
    return callback.condition(tmp,abst,integrator)
  end
end

# Use Recursion to find the first callback for type-stability

# Base Case: Only one callback
function find_first_continuous_callback(integrator, callback::AbstractContinuousCallback)
  (find_callback_time(integrator,callback,1)...,1,1)
end

# Starting Case: Compute on the first callback
function find_first_continuous_callback(integrator, callback::AbstractContinuousCallback, args...)
  find_first_continuous_callback(integrator,find_callback_time(integrator,callback,1)...,1,1,args...)
end

function find_first_continuous_callback(integrator,tmin::Number,upcrossing::Number,
                                        event_occured::Bool,event_idx::Int,idx::Int,counter::Int,
                                        callback2)
  counter += 1 # counter is idx for callback2.
  tmin2,upcrossing2,event_occurred2,event_idx2 = find_callback_time(integrator,callback2,counter)

  if event_occurred2 && (tmin2 < tmin || !event_occured)
    return tmin2,upcrossing2,true,event_idx2,counter,counter
  else
    return tmin,upcrossing,event_occured,event_idx,idx,counter
  end
end

function find_first_continuous_callback(integrator,tmin::Number,upcrossing::Number,event_occured::Bool,event_idx::Int,idx::Int,counter::Int,callback2,args...)
  find_first_continuous_callback(integrator,find_first_continuous_callback(integrator,tmin,upcrossing,event_occured,event_idx,idx,counter,callback2)...,args...)
end

@inline function determine_event_occurance(integrator,callback::VectorContinuousCallback,counter)
  event_occurred = false
  if callback.interp_points!=0
    addsteps!(integrator)
  end
  dt = (integrator.t - integrator.tprev) / callback.interp_points
  ts = integrator.tprev:dt:integrator.t
  interp_index = 0
  # Check if the event occured
  previous_condition = @views(integrator.callback_cache.previous_condition[1:callback.len])

  if callback.idxs === nothing
    callback.condition(previous_condition,integrator.uprev,integrator.tprev,integrator)
  else
    callback.condition(previous_condition,integrator.uprev[callback.idxs],integrator.tprev,integrator)
  end
  integrator.sol.destats.ncondition += 1

  ivec = integrator.vector_event_last_time
  prev_sign = @view(integrator.callback_cache.prev_sign[1:callback.len])
  next_sign = @view(integrator.callback_cache.next_sign[1:callback.len])

  if integrator.event_last_time == counter && minimum(ODE_DEFAULT_NORM(ArrayInterface.allowed_getindex(previous_condition,ivec),integrator.t)) <= 100ODE_DEFAULT_NORM(integrator.last_event_error,integrator.t)

    # If there was a previous event, utilize the derivative at the start to
    # chose the previous sign. If the derivative is positive at tprev, then
    # we treat `prev_sign` as negetive, and if the derivative is negative then we
    # treat `prev_sign` as positive, regardless of the postiivity/negativity
    # of the true value due to it being =0 sans floating point issues.

    # Only due this if the discontinuity did not move it far away from an event
    # Since near even we use direction instead of location to reset

    if callback.interp_points==0
      addsteps!(integrator)
    end

    # Evaluate condition slightly in future
    abst = integrator.tprev+integrator.dt * callback.repeat_nudge
    tmp_condition = get_condition(integrator, callback, abst)
    @. prev_sign = sign(previous_condition)
    prev_sign[ivec] = sign(tmp_condition[ivec])
  else
      @. prev_sign = sign(previous_condition)
  end

  prev_sign_index = 1
  abst = integrator.t
  next_condition = get_condition(integrator, callback, abst)
  @. next_sign = sign(next_condition)

  event_idx = findall_events!(next_sign,callback.affect!,callback.affect_neg!,prev_sign)
  if sum(event_idx) != 0
    event_occurred = true
    interp_index = callback.interp_points
  end
  if callback.interp_points!=0 && !isdiscrete(integrator.alg) && sum(event_idx) != length(event_idx) # Use the interpolants for safety checking
    fallback = true
    for i in 2:length(ts)
      abst = ts[i]
      copyto!(next_sign,get_condition(integrator, callback, abst))
      _event_idx = findall_events!(next_sign,callback.affect!,callback.affect_neg!,prev_sign)
      if sum(_event_idx) != 0
        event_occurred = true
        event_idx = _event_idx
        interp_index = i
        fallback = false
        break
      else
        prev_sign_index = i
      end
    end

    if fallback
      # If you get here, then you need to reset the event_idx to the
      # non-interpolated version

      abst = integrator.t
      next_condition = get_condition(integrator, callback, abst)
      @. next_sign = sign(next_condition)
      event_idx = findall_events!(next_sign,callback.affect!,callback.affect_neg!,prev_sign)
      interp_index = callback.interp_points
    end
  end

  event_occurred,interp_index,ts,prev_sign,prev_sign_index,event_idx
end

@inline function determine_event_occurance(integrator,callback::ContinuousCallback,counter)
  event_occurred = false
  if callback.interp_points!=0
    addsteps!(integrator)
  end
  ts = range(integrator.tprev, stop=integrator.t, length=callback.interp_points)
  interp_index = 0
  # Check if the event occured
  if callback.idxs === nothing
    previous_condition = callback.condition(integrator.uprev,integrator.tprev,integrator)
  else
    @views previous_condition = callback.condition(integrator.uprev[callback.idxs],integrator.tprev,integrator)
  end
  integrator.sol.destats.ncondition += 1

  prev_sign = 0.0
  next_sign = 0.0

  if integrator.event_last_time == counter && minimum(ODE_DEFAULT_NORM(previous_condition,integrator.t)) <= 100ODE_DEFAULT_NORM(integrator.last_event_error,integrator.t)

    # If there was a previous event, utilize the derivative at the start to
    # chose the previous sign. If the derivative is positive at tprev, then
    # we treat `prev_sign` as negetive, and if the derivative is negative then we
    # treat `prev_sign` as positive, regardless of the postiivity/negativity
    # of the true value due to it being =0 sans floating point issues.

    # Only due this if the discontinuity did not move it far away from an event
    # Since near even we use direction instead of location to reset

    if callback.interp_points==0
      addsteps!(integrator)
    end

    # Evaluate condition slightly in future
    abst = integrator.tprev+integrator.dt * callback.repeat_nudge
    tmp_condition = get_condition(integrator, callback, abst)
    prev_sign = sign(tmp_condition)
  else
    prev_sign = sign(previous_condition)
  end

  prev_sign_index = 1
  abst = integrator.t
  next_condition = get_condition(integrator, callback, abst)
  next_sign = sign(next_condition)

  if ((prev_sign < 0 && callback.affect! !== nothing) || (prev_sign > 0 && callback.affect_neg! !== nothing)) && prev_sign*next_sign<=0
    event_occurred = true
    interp_index = callback.interp_points
  elseif callback.interp_points!=0 && !isdiscrete(integrator.alg) # Use the interpolants for safety checking
    for i in 2:length(ts)
      abst = ts[i]
      new_sign = get_condition(integrator, callback, abst)
      if ((prev_sign < 0 && callback.affect! !== nothing) || (prev_sign > 0 && callback.affect_neg! !== nothing)) && prev_sign*new_sign<0
        event_occurred = true
        interp_index = i
        break
      else
        prev_sign_index = i
      end
    end
  end
  event_idx = 1

  event_occurred,interp_index,ts,prev_sign,prev_sign_index,event_idx
end

# rough implementation, needs multiple type handling
# always ensures that if r = bisection(f, (x0, x1))
# then either f(nextfloat(r)) == 0 or f(nextfloat(r)) * f(r) < 0
function bisection(f, tup, t_forward::Bool, rootfind::RootfindOpt, abstol, reltol; maxiters=1000)
  if rootfind == LeftRootFind
    NonlinearSolve.solve(NonlinearSolve.NonlinearProblem{false}(f, tup), NonlinearSolve.Falsi(), abstol=abstol, reltol=reltol).left
  else
    NonlinearSolve.solve(NonlinearSolve.NonlinearProblem{false}(f, tup), NonlinearSolve.Falsi(), abstol=abstol, reltol=reltol).right
  end
end

"""
findall_events!(next_sign,affect!,affect_neg!,prev_sign)

Modifies `next_sign` to be an array of booleans for if there is a sign change
in the interval between prev_sign and next_sign
"""
function findall_events!(next_sign::Union{Array,SubArray},affect!::F1,affect_neg!::F2,prev_sign::Union{Array,SubArray}) where {F1,F2}
  @inbounds for i in 1:length(prev_sign)
    next_sign[i] = ((prev_sign[i] < 0 && affect! !== nothing) || (prev_sign[i] > 0 && affect_neg! !== nothing)) && prev_sign[i]*next_sign[i]<=0
  end
  next_sign
end

function findall_events!(next_sign,affect!::F1,affect_neg!::F2,prev_sign) where {F1,F2}
  hasaffect::Bool = affect! !== nothing
  hasaffectneg::Bool = affect_neg! !== nothing
  f = (n,p)-> ((p < 0 && hasaffect) || (p > 0 && hasaffectneg)) && p*n<=0
  A = map!(f,next_sign,next_sign,prev_sign)
  next_sign
end

function find_callback_time(integrator,callback::ContinuousCallback,counter)
  event_occurred,interp_index,ts,prev_sign,prev_sign_index,event_idx = determine_event_occurance(integrator,callback,counter)
  if event_occurred
    if callback.condition === nothing
      new_t = zero(typeof(integrator.t))
    else
      if callback.interp_points!=0
        top_t = ts[interp_index] # Top at the smallest
        bottom_t = ts[prev_sign_index]
      else
        top_t = integrator.t
        bottom_t = integrator.tprev
      end
      if callback.rootfind != NoRootFind && !isdiscrete(integrator.alg)
        zero_func(abst, p=nothing) = get_condition(integrator, callback, abst)
        if zero_func(top_t) == 0
          Θ = top_t
        else
          if integrator.event_last_time == counter &&
            abs(zero_func(bottom_t)) <= 100abs(integrator.last_event_error) &&
            prev_sign_index == 1

            # Determined that there is an event by derivative
            # But floating point error may make the end point negative

            bottom_t += integrator.dt * callback.repeat_nudge
            sign_top = sign(zero_func(top_t))
            sign(zero_func(bottom_t)) * sign_top >= zero(sign_top) && error("Double callback crossing floating pointer reducer errored. Report this issue.")
          end
          Θ = bisection(zero_func, (bottom_t, top_t), isone(integrator.tdir), callback.rootfind, callback.abstol, callback.reltol)
          integrator.last_event_error = ODE_DEFAULT_NORM(zero_func(Θ), Θ)
        end
        #Θ = prevfloat(...)
        # prevfloat guerentees that the new time is either 1 floating point
        # numbers just before the event or directly at zero, but not after.
        # If there's a barrier which is never supposed to be crossed,
        # then this will ensure that
        # The item never leaves the domain. Otherwise Roots.jl can return
        # a float which is slightly after, making it out of the domain, causing
        # havoc.
        new_t = Θ - integrator.tprev
      elseif interp_index != callback.interp_points && !isdiscrete(integrator.alg)
        new_t = ts[interp_index] - integrator.tprev
      else
        # If no solve and no interpolants, just use endpoint
        new_t = integrator.dt
      end
    end
  else
    new_t = zero(typeof(integrator.t))
  end

  new_t,prev_sign,event_occurred,event_idx
end

function find_callback_time(integrator,callback::VectorContinuousCallback,counter)
  event_occurred,interp_index,ts,prev_sign,prev_sign_index,event_idx = determine_event_occurance(integrator,callback,counter)
  if event_occurred
    if callback.condition === nothing
      new_t = zero(typeof(integrator.t))
      min_event_idx = findfirst(isequal(1),event_idx)
    else
      if callback.interp_points!=0
        top_t = ts[interp_index] # Top at the smallest
        bottom_t = ts[prev_sign_index]
      else
        top_t = integrator.t
        bottom_t = integrator.tprev
      end
      if callback.rootfind != NoRootFind && !isdiscrete(integrator.alg)
        min_t = nextfloat(top_t)
        min_event_idx = -1
        for idx in 1:length(event_idx)
          if event_idx[idx] != 0
            zero_func(abst, p=nothing) = ArrayInterface.allowed_getindex(get_condition(integrator, callback, abst),idx)
            if zero_func(top_t) == 0
              Θ = top_t
            else
              if integrator.event_last_time == counter &&
                integrator.vector_event_last_time == idx &&
                abs(zero_func(bottom_t)) <= 100abs(integrator.last_event_error) &&
                prev_sign_index == 1

                # Determined that there is an event by derivative
                # But floating point error may make the end point negative

                bottom_t += integrator.dt * callback.repeat_nudge
                sign_top = sign(zero_func(top_t))
                sign(zero_func(bottom_t)) * sign_top >= zero(sign_top) && error("Double callback crossing floating pointer reducer errored. Report this issue.")
              end
              Θ = bisection(zero_func, (bottom_t, top_t), isone(integrator.tdir), callback.rootfind, callback.abstol, callback.reltol)
              if integrator.tdir * Θ < integrator.tdir * min_t
                integrator.last_event_error = ODE_DEFAULT_NORM(zero_func(Θ), Θ)
              end
            end
            if integrator.tdir * Θ < integrator.tdir * min_t
              min_event_idx = idx
              min_t = Θ
            end
          end
        end
        #Θ = prevfloat(...)
        # prevfloat guerentees that the new time is either 1 floating point
        # numbers just before the event or directly at zero, but not after.
        # If there's a barrier which is never supposed to be crossed,
        # then this will ensure that
        # The item never leaves the domain. Otherwise Roots.jl can return
        # a float which is slightly after, making it out of the domain, causing
        # havoc.
        new_t = min_t -integrator.tprev
      elseif interp_index != callback.interp_points && !isdiscrete(integrator.alg)
        new_t = ts[interp_index] - integrator.tprev
        min_event_idx = findfirst(isequal(1),event_idx)
      else
        # If no solve and no interpolants, just use endpoint
        new_t = integrator.dt
        min_event_idx = findfirst(isequal(1),event_idx)
      end
    end
  else
    new_t = zero(typeof(integrator.t))
    min_event_idx = 1
  end

  if event_occurred && min_event_idx < 0
    error("Callback handling failed. Please file an issue with code to reproduce.")
  end

  new_t,ArrayInterface.allowed_getindex(prev_sign,min_event_idx),event_occurred::Bool,min_event_idx::Int
end

function apply_callback!(integrator,callback::Union{ContinuousCallback,VectorContinuousCallback},cb_time,prev_sign,event_idx)

  if isadaptive(integrator)
    set_proposed_dt!(integrator, integrator.tdir * max(nextfloat(integrator.opts.dtmin), integrator.tdir * callback.dtrelax * integrator.dt))
  end

  change_t_via_interpolation!(integrator,integrator.tprev+cb_time)

  # handle saveat
  _, savedexactly = savevalues!(integrator)
  saved_in_cb = true

  @inbounds if callback.save_positions[1]
    # if already saved then skip saving
    savedexactly || savevalues!(integrator,true)
  end

  integrator.u_modified = true

  if prev_sign < 0
    if callback.affect! === nothing
      integrator.u_modified = false
    else
      callback isa VectorContinuousCallback ? callback.affect!(integrator,event_idx) : callback.affect!(integrator)
    end
  elseif prev_sign > 0
    if callback.affect_neg! === nothing
      integrator.u_modified = false
    else
      callback isa VectorContinuousCallback ? callback.affect_neg!(integrator,event_idx) : callback.affect_neg!(integrator)
    end
  end

  if integrator.u_modified
    reeval_internals_due_to_modification!(integrator)
    @inbounds if callback.save_positions[2]
      savevalues!(integrator,true)
      saved_in_cb = true
    end
    return true,saved_in_cb
  end
  false,saved_in_cb
end

#Base Case: Just one
@inline function apply_discrete_callback!(integrator,callback::DiscreteCallback)
  saved_in_cb = false
  if callback.condition(integrator.u,integrator.t,integrator)
    # handle saveat
    _, savedexactly = savevalues!(integrator)
    saved_in_cb = true
    @inbounds if callback.save_positions[1]
      # if already saved then skip saving
      savedexactly || savevalues!(integrator,true)
    end
    integrator.u_modified = true
    callback.affect!(integrator)
    @inbounds if callback.save_positions[2]
      savevalues!(integrator,true)
      saved_in_cb = true
    end
  end
  integrator.sol.destats.ncondition += 1
  integrator.u_modified,saved_in_cb
end

#Starting: Get bool from first and do next
@inline function apply_discrete_callback!(integrator,callback::DiscreteCallback,args...)
  apply_discrete_callback!(integrator,apply_discrete_callback!(integrator,callback)...,args...)
end

@inline function apply_discrete_callback!(integrator,discrete_modified::Bool,saved_in_cb::Bool,callback::DiscreteCallback,args...)
  bool,saved_in_cb2 = apply_discrete_callback!(integrator,apply_discrete_callback!(integrator,callback)...,args...)
  discrete_modified || bool, saved_in_cb || saved_in_cb2
end

@inline function apply_discrete_callback!(integrator,discrete_modified::Bool,saved_in_cb::Bool,callback::DiscreteCallback)
  bool,saved_in_cb2 = apply_discrete_callback!(integrator,callback)
  discrete_modified || bool, saved_in_cb || saved_in_cb2
end

max_vector_callback_length_int(cs::CallbackSet) = max_vector_callback_length_int(cs.continuous_callbacks...)
max_vector_callback_length_int() = nothing
function max_vector_callback_length_int(continuous_callbacks...)
  all(cb->cb isa ContinuousCallback,continuous_callbacks) && return nothing
  maxlen = -1
  for cb in continuous_callbacks
    if cb isa VectorContinuousCallback && cb.len > maxlen
      maxlen = cb.len
    end
  end
  maxlen
end

function max_vector_callback_length(cs::CallbackSet)
  continuous_callbacks = cs.continuous_callbacks
  maxlen_cb = nothing
  maxlen = -1
  for cb in continuous_callbacks
    if cb isa VectorContinuousCallback && cb.len > maxlen
      maxlen = cb.len
      maxlen_cb = cb
    end
  end
  maxlen_cb
end

"""
$(TYPEDEF)
"""
mutable struct CallbackCache{conditionType,signType}
  tmp_condition::conditionType
  previous_condition::conditionType
  next_sign::signType
  prev_sign::signType
end

function CallbackCache(u,max_len,::Type{conditionType},::Type{signType}) where {conditionType,signType}
    tmp_condition = similar(u, conditionType, max_len)
    previous_condition = similar(u, conditionType, max_len)
    next_sign = similar(u, signType, max_len)
    prev_sign = similar(u, signType, max_len)
    CallbackCache(tmp_condition,previous_condition,next_sign,prev_sign)
end

function CallbackCache(max_len,::Type{conditionType},::Type{signType}) where {conditionType,signType}
    tmp_condition = zeros(conditionType, max_len)
    previous_condition = zeros(conditionType, max_len)
    next_sign = zeros(signType, max_len)
    prev_sign = zeros(signType, max_len)
    CallbackCache(tmp_condition,previous_condition,next_sign,prev_sign)
end


text:
# The functions with underscore in the name are manual implementations, the others are cudnn calls.

using CUDA, NNlib, BenchmarkTools, Knet, Test
using NNlib: softmax, softmax!, ∇softmax, ∇softmax!, logsoftmax, logsoftmax!, ∇logsoftmax, ∇logsoftmax!
using CUDA.CUDNN: _softmax!, _∇softmax!, _logsoftmax!, _∇logsoftmax!

x,y,dx,dy = (CUDA.randn(1024, 256) for i in 1:4)

macro ctime(ex)
    quote
        println($(sprint(Base.show_unquoted,ex)))
        @btime (CUDA.@sync $ex; nothing)
        CUDA.@time $ex
    end |> esc
end

@test softmax!(similar(y),x;dims=1) ≈ _softmax!(similar(y),x;dims=1)
@test softmax!(similar(y),x;dims=2) ≈ _softmax!(similar(y),x;dims=2)
@test softmax!(similar(y),x;dims=:) ≈ _softmax!(similar(y),x;dims=:)
@test logsoftmax!(similar(y),x;dims=1) ≈ _logsoftmax!(similar(y),x;dims=1)
@test logsoftmax!(similar(y),x;dims=2) ≈ _logsoftmax!(similar(y),x;dims=2)
@test logsoftmax!(similar(y),x;dims=:) ≈ _logsoftmax!(similar(y),x;dims=:)
@test ∇softmax!(similar(dx),dy,x,softmax(x;dims=1);dims=1) ≈ _∇softmax!(similar(dx),dy,x,softmax(x;dims=1);dims=1)
@test ∇softmax!(similar(dx),dy,x,softmax(x;dims=2);dims=2) ≈ _∇softmax!(similar(dx),dy,x,softmax(x;dims=2);dims=2)
@test ∇softmax!(similar(dx),dy,x,softmax(x;dims=:);dims=:) ≈ _∇softmax!(similar(dx),dy,x,softmax(x;dims=:);dims=:)
@test ∇logsoftmax!(similar(dx),dy,x,logsoftmax(x;dims=1);dims=1) ≈ _∇logsoftmax!(similar(dx),dy,x,logsoftmax(x;dims=1);dims=1)
@test ∇logsoftmax!(similar(dx),dy,x,logsoftmax(x;dims=2);dims=2) ≈ _∇logsoftmax!(similar(dx),dy,x,logsoftmax(x;dims=2);dims=2)
@test ∇logsoftmax!(similar(dx),dy,x,logsoftmax(x;dims=:);dims=:) ≈ _∇logsoftmax!(similar(dx),dy,x,logsoftmax(x;dims=:);dims=:)

GC.gc(true)
if true
@ctime softmax!(y,x;dims=1)
@ctime _softmax!(y,x;dims=1)
@ctime softmax!(y,x;dims=2)
@ctime _softmax!(y,x;dims=2)
@ctime softmax!(y,x;dims=:)
@ctime _softmax!(y,x;dims=:)
end

#=
softmax!(y, x; dims = 1)
  29.702 μs (24 allocations: 624 bytes)
  0.015599 seconds (7.92 k CPU allocations: 438.803 KiB)
_softmax!(y, x; dims = 1)
  43.810 μs (51 allocations: 1.41 KiB)
  0.000085 seconds (60 CPU allocations: 1.750 KiB) (2 GPU allocations: 2.000 KiB, 2.58% gc time)
softmax!(y, x; dims = 2)
  158.163 μs (24 allocations: 624 bytes)
  0.000186 seconds (33 CPU allocations: 976 bytes)
_softmax!(y, x; dims = 2)
  69.622 μs (59 allocations: 1.56 KiB)
  0.000104 seconds (68 CPU allocations: 1.906 KiB) (2 GPU allocations: 8.000 KiB, 2.52% gc time)
softmax!(y, x; dims = (:))
  399.958 μs (24 allocations: 624 bytes)
  0.000498 seconds (31 CPU allocations: 944 bytes)
_softmax!(y, x; dims = (:))
  83.624 μs (37 allocations: 1.00 KiB)
  0.000127 seconds (44 CPU allocations: 1.312 KiB)
=#

GC.gc(true)
if true
@ctime logsoftmax!(y,x;dims=1)
@ctime _logsoftmax!(y,x;dims=1)
@ctime logsoftmax!(y,x;dims=2)
@ctime _logsoftmax!(y,x;dims=2)
@ctime logsoftmax!(y,x;dims=:)
@ctime _logsoftmax!(y,x;dims=:)
end

#=
logsoftmax!(y, x; dims = 1)
  30.771 μs (24 allocations: 624 bytes)
  0.000073 seconds (33 CPU allocations: 976 bytes)
_logsoftmax!(y, x; dims = 1)
  60.873 μs (71 allocations: 1.97 KiB)
  0.000126 seconds (80 CPU allocations: 2.312 KiB) (3 GPU allocations: 1.002 MiB, 3.46% gc time)
logsoftmax!(y, x; dims = 2)
  130.960 μs (24 allocations: 624 bytes)
  0.000202 seconds (33 CPU allocations: 976 bytes)
_logsoftmax!(y, x; dims = 2)
  78.876 μs (79 allocations: 2.12 KiB)
  0.000107 seconds (88 CPU allocations: 2.469 KiB) (3 GPU allocations: 1.008 MiB, 3.27% gc time)
logsoftmax!(y, x; dims = (:))
  349.391 μs (24 allocations: 624 bytes)
  0.000366 seconds (31 CPU allocations: 944 bytes)
_logsoftmax!(y, x; dims = (:))
  136.404 μs (57 allocations: 1.56 KiB)
  0.000153 seconds (64 CPU allocations: 1.875 KiB) (1 GPU allocation: 1024.000 KiB, 1.47% gc time)
=#

GC.gc(true)
if true
@ctime ∇softmax!(dx,dy,x,y;dims=1)
@ctime _∇softmax!(dx,dy,x,y;dims=1)
@ctime ∇softmax!(dx,dy,x,y;dims=2)
@ctime _∇softmax!(dx,dy,x,y;dims=2)
@ctime ∇softmax!(dx,dy,x,y;dims=:)
@ctime _∇softmax!(dx,dy,x,y;dims=:)
end

#=
∇softmax!(dx, dy, x, y; dims = 1)
  39.541 μs (25 allocations: 736 bytes)
  0.000077 seconds (34 CPU allocations: 1.062 KiB)
_∇softmax!(dx, dy, x, y; dims = 1)
  103.781 μs (50 allocations: 1.44 KiB)
  0.000090 seconds (59 CPU allocations: 1.781 KiB) (2 GPU allocations: 1.001 MiB, 2.14% gc time)
∇softmax!(dx, dy, x, y; dims = 2)
  121.838 μs (25 allocations: 736 bytes)
  0.000141 seconds (34 CPU allocations: 1.062 KiB)
_∇softmax!(dx, dy, x, y; dims = 2)
  52.789 μs (54 allocations: 1.52 KiB)
  0.000092 seconds (63 CPU allocations: 1.859 KiB) (2 GPU allocations: 1.004 MiB, 2.39% gc time)
∇softmax!(dx, dy, x, y; dims = (:))
  293.816 μs (25 allocations: 736 bytes)
  0.000323 seconds (32 CPU allocations: 1.031 KiB)
_∇softmax!(dx, dy, x, y; dims = (:))
  61.927 μs (43 allocations: 1.23 KiB)
  0.000096 seconds (50 CPU allocations: 1.547 KiB) (1 GPU allocation: 1024.000 KiB, 1.38% gc time)
=#

GC.gc(true)
if true
@ctime ∇logsoftmax!(dx,dy,x,y;dims=1)
@ctime _∇logsoftmax!(dx,dy,x,y;dims=1)
@ctime ∇logsoftmax!(dx,dy,x,y;dims=2)
@ctime _∇logsoftmax!(dx,dy,x,y;dims=2)
@ctime ∇logsoftmax!(dx,dy,x,y;dims=:)
@ctime _∇logsoftmax!(dx,dy,x,y;dims=:)
end

#=
∇logsoftmax!(dx, dy, x, y; dims = 1)
  39.390 μs (25 allocations: 736 bytes)
  0.000082 seconds (34 CPU allocations: 1.062 KiB)
_∇logsoftmax!(dx, dy, x, y; dims = 1)
  65.409 μs (66 allocations: 4.25 KiB)
  0.000156 seconds (75 CPU allocations: 4.594 KiB) (1 GPU allocation: 1024 bytes, 1.59% gc time)
∇logsoftmax!(dx, dy, x, y; dims = 2)
  92.615 μs (25 allocations: 736 bytes)
  0.000164 seconds (34 CPU allocations: 1.062 KiB)
_∇logsoftmax!(dx, dy, x, y; dims = 2)
  40.642 μs (70 allocations: 4.33 KiB)
  0.000091 seconds (79 CPU allocations: 4.672 KiB) (1 GPU allocation: 4.000 KiB, 2.18% gc time)
∇logsoftmax!(dx, dy, x, y; dims = (:))
  311.676 μs (25 allocations: 736 bytes)
  0.000396 seconds (32 CPU allocations: 1.031 KiB)
_∇logsoftmax!(dx, dy, x, y; dims = (:))
  85.144 μs (59 allocations: 3.36 KiB)
  0.000136 seconds (66 CPU allocations: 3.672 KiB)
=#

nothing


text:
# UE computation
function updatechoices!()#(ECᵢ, DTC, SR)
    global ECᵢ#, DTC, SR
    newDTC = zero(DTC)
    #newSR = SR
    newSR = Dict(d => Dict(i => zeros(T, nsinks, nclasses) for i in 1:2) for d in divs);
    λ = 1e-4
    for (srcid,src) in enumerate(srcs)
        i = outlinkids(net, src)[1]
        for (snkid,snk) in enumerate(snks)
            for clsid in 1:nclasses
                c = ECᵢ[i,1:Tm,snkid,clsid]
                dtc = DTC[1:Tm, srcid,snkid,clsid]
                newdtc = proportionalize(solveqp(c, dtc, λ))[:,1]
                newDTC[1:Tm, srcid,snkid,clsid] .= newdtc
            end
        end
    end

    for div in divs
        oli = outlinkids(net, div)
        for t in 1:T
            for (snkid,snk) in enumerate(snks)
                for clsid in 1:nclasses
                    sr = [SR[div][i][t, snkid,clsid] for i in 1:2]
                    c = ECᵢ[oli,t,snkid,clsid]
                    newsr = proportionalize(solveqp(c, sr, λ))[:,1]
                    for i in 1:2
                        newSR[div][i][t,snkid,clsid] = newsr[i]
                    end
                end
            end
        end
    end

    return newDTC, newSR
end

function relgap(ECᵢ, DTC)
    _num, _den = 0., 0.
    for (srcid,src) in enumerate(srcs)
        i = outlinkids(net, src)[1]
        for (snkid,snk) in enumerate(snks)
            for clsid in 1:nclasses
                m = minimum(ECᵢ[i,1:Tm,snkid,clsid])
                a = sum(ECᵢ[i,1:Tm,snkid,clsid] .* DTC[1:Tm, srcid,snkid,clsid])
                _num += (a - m)
                _den += m
            end
        end
    end
    #@show (_num/_den)
    return _num/_den
end

# equilibrium computation
function optfn!(x, grad, p)
    if length(grad) > 0
        grad .= (x .- p)
    end
    return 0.5*sum((x .- p).^2)
end

function eqcon!(x, grad)
    if length(grad) > 0
        grad .= 1
    end
    return sum(x) - 1
end

function solveqp(c, dtc, λ)
    p = dtc .- λ * c
    """x = Variable(length(p))
    problem = minimize(0.5*sumsquares(x - p), sum(x) == 1, x >= 0.)
    solve!(problem, () -> SCS.Optimizer(verbose=false))
    if (problem.status != MOI.OPTIMAL)
        @show c, dtc, λ
    end
    return x.value"""
    return projpx(p)
end

function proportionalize(x; digits=6)
    x[x .< 0.] .= 0
    x = round.(x, digits=digits)
    x ./= sum(x)
    return x
end


text:
# BGZF Decompressor
# =================

mutable struct BGZFDecompressor <: TranscodingStreams.Codec
    zstream::CodecZlib.ZStream
    windowbits::Int
    blockstart::Bool
end

function BGZFDecompressor(;windowbits::Integer=CodecZlib.Z_DEFAULT_WINDOWBITS)
    return BGZFDecompressor(CodecZlib.ZStream(), windowbits, true)
end

function Base.show(io::IO, codec::BGZFDecompressor)
    print(io, summary(codec), "(windowbits=$(windowbits))")
end

const BGZFDecompressorStream{S} = TranscodingStream{BGZFDecompressor,S} where S<:IO

function BGZFDecompressorStream(stream::IO)
    return TranscodingStream(BGZFDecompressor(), stream; bufsize=BGZF_MAX_BLOCK_SIZE)
end

"""
    seek(stream::BGZFDecompressorStream, voffset::UInt64)

Seek `stream` to the virtual file offset `voffset`.

A virtual file offset is an unsigned 64-bit integer, which is defined as
`coffset<<16|uoffset`, where `coffset` is an offset into the BGZF file to the
beginning of a block, and `uoffset` is an offset into the uncompressed data.

Seeking a stream to a position from which a BGZF block does not start will
result in an error.

This method is not available in the write mode.
"""
function Base.seek(stream::BGZFDecompressorStream, voffset::UInt64)
    TranscodingStreams.changemode!(stream, :read)
    seek(stream.stream, voffset >> 16)
    TranscodingStreams.initbuffer!(stream.state.buffer1)
    TranscodingStreams.initbuffer!(stream.state.buffer2)
    skip(stream, voffset & 0xffff)
    return
end

# Methods
# -------

function TranscodingStreams.initialize(codec::BGZFDecompressor)
    code = CodecZlib.inflate_init!(codec.zstream, codec.windowbits+16)
    if code != CodecZlib.Z_OK
        CodecZlib.zerror(codec.zstream, code)
    end
    return
end

function TranscodingStreams.finalize(codec::BGZFDecompressor)
    zstream = codec.zstream
    if zstream.state != C_NULL
        code = CodecZlib.inflate_end!(zstream)
        if code != CodecZlib.Z_OK
            CodecZlib.zerror(zstream, code)
        end
    end
    return
end

function TranscodingStreams.minoutsize(::BGZFDecompressor, ::Memory)
    return Int(BGZF_MAX_BLOCK_SIZE)
end

function TranscodingStreams.process(codec::BGZFDecompressor, input::Memory, output::Memory, error::Error)
    zstream = codec.zstream
    if codec.blockstart
        code = CodecZlib.inflate_reset!(zstream)
        if code != CodecZlib.Z_OK
            error[] = ErrorException(CodecZlib.zlib_error_message(zstream, code))
            return 0, 0, :error
        end
    end
    zstream.next_in = input.ptr
    zstream.avail_in = input.size
    zstream.next_out = output.ptr
    zstream.avail_out = output.size
    code = CodecZlib.inflate!(zstream, CodecZlib.Z_NO_FLUSH)
    Δin = Int(input.size - zstream.avail_in)
    Δout = Int(output.size - zstream.avail_out)
    if code == CodecZlib.Z_OK
        codec.blockstart = false
        return Δin, Δout, :ok
    elseif code == CodecZlib.Z_STREAM_END
        codec.blockstart = true
        return Δin, Δout, :end
    else
        error[] = ErrorException(CodecZlib.zlib_error_message(zstream, code))
        return Δin, Δout, :error
    end
end


text:
using QuantumCumulants
using SymbolicUtils
using Test

@testset "average" begin

hf = FockSpace(:cavity)
ha = NLevelSpace(:atom,(:g,:e))
h = hf⊗ha

a = Destroy(h,:a)
σ = Transition(h,:σ,:g,:e)

@test isequal(average(2a), 2*average(a))

@test isequal(QuantumCumulants._conj(average(2im*a'*σ)), (-2im)*average(a*σ'))
@test isequal(average(2*(a+a)), 2*(average(a) + average(a)))
@test isequal(average(a^2), average(a*a))
@test isequal(average(a^2,1), average(a)^2)
@test isequal(average((a'*a)^2), average(a'*a*a'*a))
@test isequal(average(a'*a^2,2), average(a')*average(a^2) + -2*average(a')*average(a)^2 + 2*average(a)*average(a'*a))

@test isequal(simplify(average(σ)+average(σ)), average(2σ))
@test iszero(simplify(average(σ)-average(σ)))

ωc, ωa = cnumbers("ω_c ω_a")
@test isequal(average(ωc),ωc)
@test isequal(average(ωc*a),ωc*average(a))
@test isequal(average(ωc*(a+a')) , ωc*average(a) + ωc*average(a'))

n = average(a'*a)
@test isequal(cumulant_expansion(n,2), n)
@test isequal(cumulant_expansion(n,1), average(a')*average(a))

@test iszero(average(a*σ)*average(a) - average(a)*average(a*σ))

# Test cumulants
hs = FockSpace[]
for i=1:4
    push!(hs, FockSpace(Symbol(:fock, i)))
end
h = ⊗(hs...)

a = Destroy(h,:a,1)
b = Destroy(h,:b,2)
c = Destroy(h,:c,3)
d = Destroy(h,:d,4)

@test isequal(cumulant(a*b),simplify(average(a*b)+ -1*average(a)*average(b)))
@test isequal(cumulant(a*b,1), average(a*b))
@test iszero(simplify(cumulant(a*b*c) - (average(a*b*c) +
            2*average(a)*average(b)*average(c) - average(a)*average(b*c) -
            average(b)*average(a*c) - average(c)*average(a*b))))
@test isequal(simplify(average(a*b*c*d) - cumulant(a*b*c*d)), cumulant_expansion(average(a*b*c*d),3))

end # testset


text:
using Test
using MPI

MPI.Init()
comm = MPI.COMM_WORLD
nnodes = MPI.Comm_size(comm)
ndims = 3

cdims = Cint[0,0,0]
MPI.Dims_create!(nnodes, ndims, cdims)
@test prod(cdims) == nnodes

cdims = Cint[1,0,1]
MPI.Dims_create!(nnodes, ndims, cdims)
@test cdims == Cint[1,nnodes,1]

for dims in ([0,0,0], [0 0 0])
    MPI.Dims_create!(nnodes, dims)
    @test prod(dims) == nnodes
end

MPI.Finalize()