text:
using Inflation
using SymPy
using LinearAlgebra
using Test

@testset "inflation models" begin
    @testset "single-field quadratic" begin
        d = 1
        function G(Phi,p)
            return [1.0]
        end

        m = symbols("m")
        params = [m]
        pvalues = [1.0e-10]
        function V(Phi,p)
            g = G(Phi,p)
            m = p[1]
            v = 0.0
            for i in 1:d
                for j in 1:d
                    v += (Phi[i] * m * Phi[j])/2.0
                end
            end
            return v
        end

        opt = SetupOptions(
        allow_caching = false
        )
        funcs = inflation_setup(d,V,G,params,options=opt)
        Phi0 = [18.0]
        Pi0  = [0.0]

        sol1 = background_evolve(Phi0,Pi0,pvalues,funcs,verbose=true)
        @test sol1.retcode == :Terminated
        @test sol1.t[end] ≈ 81.9 rtol=0.02
        eh = funcs["Eh"]([sol1.u[end];pvalues]...)
        @test eh ≈ 1.0 atol=1e-8
        tsol = transport_perturbations(sol1,pvalues,funcs,verbose=true)
        ns = tsol["ns"]
        as = tsol["as"]
        As = tsol["As"]
        r = tsol["r"]
        # results from the same model in Multi-Mode-Code
        MMC_Ne = 81.607534
        #MMC_al = 14.67
        MMC_As = 5.098295961e-9
        MMC_Aiso = 1.591443943
        MMC_Pent = 9.881312917
        MMC_Pcross = -10
        MMC_r = 0.1440518
        MMC_ns = 0.9633723789
        MMC_alpha_s = -6.854134226e-04
        MMC_f_NL = 2.842
        @test ns ≈ MMC_ns rtol=0.01
        @test as ≈ MMC_alpha_s rtol=0.05
        @test As ≈ MMC_As rtol=0.05
        @test r ≈ MMC_r rtol=0.02
    end
    @testset "single-field α-attractors" begin
        d = 1
        α = symbols("a",real=true)
        m = symbols("m",real=true)
        params = [α,m]
        pvalues = [1/2,1e-5]

        function G(Phi,p)
            return [6*α/(1-Phi[1]^2)^2]
        end
        function V(Phi,p)
            return α/2.0*(m^2*Phi[1]^2)
        end
        Phi0 = [.999]
        Pi0 = [0.0]

        opt = SetupOptions(
        allow_caching = false,
        attempt_HCA = false
        )

        funcs = inflation_setup(d,V,G,params,options=opt)
        sol1 = background_evolve(Phi0,Pi0,pvalues,funcs,verbose=true)

        #analytics
        Nexit_to_end = 55
        anal_ns = 1 - 2/Nexit_to_end
        anal_r  = 12*pvalues[1]/Nexit_to_end^2
        @test sol1.retcode == :Terminated
        @test sol1.t[end] ≈ 376.666 rtol=0.002
        tsol = transport_perturbations(sol1,pvalues,funcs,verbose=true)
        ns = tsol["ns"]
        as = tsol["as"]
        As = tsol["As"]
        r = tsol["r"]
        @test ns ≈ anal_ns rtol=0.02
        @test r ≈ anal_r rtol=0.08
        #@test as ≈ MMC_alpha_s rtol=0.05
        #@test As ≈ MMC_As rtol=0.05
    end
    @testset "double-field α-attractors" begin
        "see fig. 1 of 1803.09841 for this model"
        d = 2
        α,m,R = symbols("a m r",real=true)
        params = [α,m,R]
        pvalues = [1/6,1e-5,10]
        #α = 1/600
        function G(Phi,p)
            α = p[1]
            g = [ 6*α/(1-Phi[1]^2 - Phi[2]^2)^2 0; 0 6*α/(1-Phi[1]^2 - Phi[2]^2)^2]
            return g
        end
        #function Ginv(Phi,p)
        #    α = p[1]
        #    gi = [ :(1/(6*$α)*(1-$(Phi[1])^2 - $(Phi[2])^2)^2) :(0); :(0) :(1/(6*$α)*(1-$(Phi[1])^2 - $(Phi[2])^2)^2)]
        #    return gi
        #end
        function V(Phi,p)
            m = [p[2]^2 0.0; 0.0 p[3]*p[2]^2]
            v = 0
            for i in 1:d
                for j in 1:d
                    v += α/2.0*(m[i,j]*Phi[i]*Phi[j])
                end
            end
            return v
        end
        r0 = 0.999
        #r0 = 0.99
        θ0 = pi/4
        Phi0 = [r0*cos(θ0), r0*sin(θ0)]
        Pi0 = [0.0, 0.0]

        opt = SetupOptions(
        allow_caching = false,
        attempt_HCA = false
        )
        funcs = inflation_setup(d,V,G,params,options=opt)#,Ginv=Ginv)
        sol1 = background_evolve(Phi0,Pi0,pvalues,funcs,verbose=true)

        #analytics
        Nexit_to_end = 55
        Nang = 0
        Neff = Nexit_to_end - Nang
        anal_ns = 1 - 2/Neff
        anal_r  = 12*pvalues[1]/Neff^2
        @test sol1.retcode == :Terminated
        @test sol1.t[end] ≈ 129.567 rtol=0.002
        tsol = transport_perturbations(sol1,pvalues,funcs,verbose=true)
        ns = tsol["ns"]
        as = tsol["as"]
        As = tsol["As"]
        r = tsol["r"]
        @test ns ≈ anal_ns rtol=0.002
        @test r ≈ anal_r rtol=0.12 # this is about 10% higher than analytic estimate
        #@test as ≈ MMC_alpha_s rtol=0.05
        #@test As ≈ MMC_As rtol=0.05
    end
    @testset "2-field hyperinflation" begin
        d = 2
        L = symbols("l")
        params = [L]
        pvalues = [0.01]
        function G(Phi,p)
            #return [1 0;0 L^2*sinh(Phi[1]/L)^2]
            L = p[1]
            sh = sinh((Phi[1])/L)
            return [1 0;0 L^2*sh^2]
        end

        function V(Phi,p)
            L = p[1]
            return L^2 * Phi[1]^2 / 2
        end
        opt = SetupOptions(
        allow_caching = false,
        attempt_HCA = false
        )
        funcs = inflation_setup(d,V,G,params,options=opt)
        Phi0 = [0.6925,1e-24]
        Pi0 = [-0.0297881,4.8921e-29]

        sol2 = background_evolve(Phi0,Pi0,pvalues,funcs,verbose=true)
        @test sol2.retcode == :Terminated
        @test sol2.t[end] ≈ 23.05 rtol=0.02
        eh = funcs["Eh"]([sol2.u[end][1:end];pvalues]...)
        @test eh ≈ 1.0 atol=1e-8
    end
    @testset "2-field quadratic" begin
        d = 2
        params = [symbols(x,real=true,positive=true) for x in ["m","r"]]
        pvalues = [1e-5,10]
        function G(Phi,p)
            return Matrix{Int64}(I,d,d)
        end

        function V(Phi,p)
            m = [p[1]^2 0 ; 0 p[2]*p[1]^2]
            v = SymPy.Sym(0.0)
            for i in 1:d
                for j in 1:d
                    v += (Phi[i] * m[i,j] * Phi[j])/2.0
                end
            end
            return v
        end

        Phi0 = [15.0,10.0]
        Pi0  = [0.0,0.0]
        opt = SetupOptions(
        allow_caching = false
        )
        funcs = inflation_setup(d,V,G,params,options=opt)
        sol2 = background_evolve(Phi0,Pi0,pvalues,funcs,verbose=true)
        @test sol2.retcode == :Terminated
        @test sol2.t[end] ≈ 82.407 rtol=0.02
        eh = funcs["Eh"]([sol2.u[end][1:end];pvalues]...)
        @test eh ≈ 1.0 atol=1e-8
        mT_Ne = 82.407
        mT_As = 9.8619e-9
        mT_ns = 0.919087
        mT_r = 0.144267
        tsol = transport_perturbations(sol2,pvalues,funcs,verbose=true)
        ns = tsol["ns"]
        as = tsol["as"]
        As = tsol["As"]
        r = tsol["r"]
        @test ns ≈ mT_ns rtol=0.02
        #@test as ≈ MMC_alpha_s rtol=0.05
        #MMC_alpha_s = -6.535329799e-04
        #MMC_f_NL = -1.5301e-02
        @test As ≈ mT_As rtol=0.05
        @test mT_r ≈ r rtol=0.02
    end
    @testset "5-field quadratic" begin
        d = 5
        params = [symbols("m$i",real=true,positive=true) for i in 1:d]
        pvalues = [1e-5,1e-5,sqrt(1e-9),1e-5,1e-4]
        function G(Phi,p)
            return Matrix{Int64}(I,d,d)
        end

        function V(Phi,p)
            m = p.^2
            v = 0.0
            for i in 1:d
                v += (Phi[i]^2 * m[i])/2.0
            end
            return v
        end

        Phi0 = [15.0,10.0,-5.0,10.0,3.0]
        Pi0  = [0.0,0.0,0.0,0.0,0.0]

        opt = SetupOptions(
        allow_caching = false
        )
        funcs = inflation_setup(d,V,G,params,options=opt)
        sol5 = background_evolve(Phi0,Pi0,pvalues,funcs,verbose=true)
        @test sol5.retcode == :Terminated
        @test sol5.t[end] ≈ 115.62897 rtol=0.02
        eh = funcs["Eh"]([sol5.u[end][1:end];pvalues]...)
        @test eh ≈ 1.0 atol=1e-8
        MMC_Ne = 1.1562897e+02
        #MMC_al = 14.67
        MMC_As = 5.105241913e-09
        MMC_Aiso = 4.217972911e-13
        MMC_Pent = 0
        MMC_Pcross = -2.026595982e-20
        MMC_r = 0.1440332439
        MMC_ns = 0.9631782590
        MMC_alpha_s = -6.535329799e-04
        MMC_f_NL = -1.5301e-02
        tsol = transport_perturbations(sol5,pvalues,funcs,verbose=true)
        ns = tsol["ns"]
        as = tsol["as"]
        As = tsol["As"]
        r = tsol["r"]
        @test ns ≈ MMC_ns rtol=0.01
        @test as ≈ MMC_alpha_s rtol=0.05
        @test As ≈ MMC_As rtol=0.05
        @test r ≈ MMC_r rtol=0.02
    end
end


text:
@userplot triplot

@recipe function f(h::triplot;
                    NodeName=["A","B","C"],m=9,tlim=[0.1,0.9], Norm=false)
    x,y = h.args

    if Norm==true
       x=x./sum(x,dims=2)

    end

    N=[1 0 0;0 1 0; 0 0 1;1 0 0]
    T=[0 0;0.5 sqrt(3)/2;1 0]  # Transformation matrix

    i=x*T  # Cartesian coordinates of given data
    n=N*T  # Cartesian coordinates of nodes

#Minor grid construction

c=zeros(m,3)
a=zeros(m,3)
b=zeros(m,3)
ba=zeros(m,3)

r=range(tlim[1],stop=tlim[2],length=m)

c[:,1]=r[end:-1:1]
c[:,3]=r

a[:,2]=r
a[:,1]=r[end:-1:1]

b[:,3]=r
b[:,2]=r[end:-1:1]

ba[:,2]=r
ba[:,3]=r[end:-1:1]

A=a*T
B=b*T
C=c*T
BA=ba*T

    legend := false

   grid := false
   showaxis := false
    @series begin
            seriestype := :scatter
           markershape --> :circle
           markersize --> 5
           markercolor --> :red
            i[:,1], i[:,2]
    end

    @series begin
       seriestype := :path
       linecolor := :black
        linewidth := 2
        series_annotation := [text("$(NodeName[1])",:right),text("$(NodeName[2])",:bottom),
                                text("$(NodeName[3])",:top),text("$(NodeName[1])",:right)]
        n[:,1],n[:,2]
    end

    @series begin
        seriestype := :path
        linecolor := :gray
        linewidth := 0.5
        [A[:,1] C[:,1]]',[A[:,2] C[:,2]]'
    end

        @series begin
        seriestype := :path
        linecolor := :gray
        linewidth := 0.5
        [C[:,1] B[:,1]]',[C[:,2] B[:,2]]'
    end

            @series begin
        seriestype := :path
        linecolor := :gray
        linewidth := 0.5
        [BA[:,1] A[:,1]]',[BA[:,2] A[:,2]]'
    end
end

ext: jl


text:
using BayesNets
using DataFrames
using Base.Test
using LightGraphs

srand(0)

"""
A simple variant of isapprox that is true if the isapprox comparison works
elementwise in the vector
"""
function elementwise_isapprox{F<:AbstractFloat}(x::AbstractArray{F},
        y::AbstractArray{F},
        rtol::F=sqrt(eps(F)),
        atol::F=zero(F))

    if length(x) != length(y)
        return false
    end

    for (a,b) in zip(x,y)
        if !isapprox(a,b,rtol=rtol, atol=atol)
            return false
        end
    end

    true
end

testdir = joinpath(dirname(@__DIR__), "test")
include(joinpath(testdir, "test_utils.jl"))
include(joinpath(testdir, "test_cpds.jl"))
include(joinpath(testdir, "test_tables.jl"))
include(joinpath(testdir, "test_factors.jl"))
include(joinpath(testdir, "test_bayesnets.jl"))
include(joinpath(testdir, "test_gibbs.jl"))
include(joinpath(testdir, "test_sampling.jl"))
include(joinpath(testdir, "test_inference.jl"))
include(joinpath(testdir, "test_learning.jl"))
include(joinpath(testdir, "test_io.jl"))
include(joinpath(testdir, "test_ndgrid.jl"))

include(joinpath(testdir, "test_discrete_bayes_nets.jl"))
include(joinpath(testdir, "test_genbn.jl"))

include(joinpath(testdir, "test_docs.jl"))


text:
Base.@kwdef struct LearningMethodRegistry
    methods::Dict{Any,Any} = Dict()
end

"""
    registerlearningmethod!(registry, methodfn, blocktypes)

Register a learning method constructor `methodfn` as compatible with
`blocktypes` in `registry::LearningMethodRegistry`.

`blocks` should be the least specific set of types that a `methodfn`
can handle. `methodfn` needs to have a method that takes concrete block
instances as the only non-keyword argument, i.e. `methodfn(blocks; kwargs...)`.
"""
function registerlearningmethod!(reg::LearningMethodRegistry, f, blocktypes)
	reg.methods[f] = typify(blocktypes)
	return reg
end

"""
	findlearningmethods(blocktypes)
	findlearningmethods(registry, blocktypes)

Find learning methods compatible with block types `TBlocks` in
`registry::LearningMethodRegistry`.

#### Examples

```julia
julia> findlearningmethods((Image, Label))
[ImageClassificationSingle,]

julia> findlearningmethods((Image, Any))
[ImageClassificationSingle, ImageClassificationMulti, ImageSegmentation, ImageKeypointRegression, ...]
```
"""
function findlearningmethods(reg::LearningMethodRegistry, blocktypes=Any)
	return [methodfn for (methodfn, methodblocks) in reg.methods if methodblocks <: typify(blocktypes)]
end

const FASTAI_METHOD_REGISTRY = LearningMethodRegistry()


text:
"""
An empty interactor. This is an empty type that has no functionality apart
from 'existing'. Its use case is testing and debugging population-dynamic
functions; linguistically, empty interactors are inert. They can be rendezvoused,
technically, but nothing happens.
"""
mutable struct EmptyInteractor <: AbstractInteractor
  meta::Any
end

"""
Construct an empty interactor.
"""
function EmptyInteractor()
  EmptyInteractor(nothing)
end

"""
    act!(x::EmptyInteractor, y::AbstractInteractor)

Make an `EmptyInteractor` attempt an action on another interactor. Does nothing
and returns `false`.
"""
function act!(x::EmptyInteractor, y::AbstractInteractor)
  return false
end

"""
    act!(x::AbstractInteractor, y::EmptyInteractor)

Make an interactor attempt an action on an `EmptyInteractor`. Does nothing
and returns `false`.
"""
function act!(x::AbstractInteractor, y::EmptyInteractor)
  return false
end

"""
    act!(x::EmptyInteractor, y::EmptyInteractor)

Make an `EmptyInteractor` attempt an action on another `EmptyInteractor`. Does nothing
and returns `false`.
"""
function act!(x::EmptyInteractor, y::EmptyInteractor)
  return false
end


text:
# indicator of the L0 norm ball with given (integer) radius

export IndBallL0

"""
**Indicator of a ``L_0`` pseudo-norm ball**

    IndBallL0(r=1)

Returns the indicator function of the set
```math
S = \\{ x : \\mathrm{nnz}(x) \\leq r \\}.
```
Parameter `r` must be a positive integer.
"""
struct IndBallL0{I <: Integer} <: ProximableFunction
  r::I
  function IndBallL0{I}(r::I) where {I <: Integer}
    if r <= 0
      error("parameter r must be a positive integer")
    else
      new(r)
    end
  end
end

is_set(f::IndBallL0) = true

IndBallL0(r::I) where {I <: Integer} = IndBallL0{I}(r)

function (f::IndBallL0)(x::AbstractArray{T}) where T <: RealOrComplex
  if count(!isequal(0), x) > f.r
    return +Inf
  end
  return zero(T)
end

function prox!(y::AbstractArray{T}, f::IndBallL0, x::AbstractArray{T}, gamma::Real=1.0) where T <: RealOrComplex
  p = []
  if ndims(x) == 1
    p = partialsortperm(x, 1:f.r, by=abs, rev=true)
  else
    p = partialsortperm(x[:], 1:f.r, by=abs, rev=true)
  end
  sort!(p)
  idx = 1
  for i = 1:length(p)
    y[idx:p[i]-1] .= zero(T)
    y[p[i]] = x[p[i]]
    idx = p[i]+1
  end
  y[idx:end] .= zero(T)
  return zero(T)
end

fun_name(f::IndBallL0) = "indicator of an L0 pseudo-norm ball"
fun_dom(f::IndBallL0) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::IndBallL0) = "x ↦ 0 if countnz(x) ⩽ r, +∞ otherwise"
fun_params(f::IndBallL0) = "r = $(f.r)"

function prox_naive(f::IndBallL0, x::AbstractArray{T}, gamma::Real=1.0) where T <: RealOrComplex
  p = sortperm(abs.(x)[:], rev=true)
  y = similar(x)
  y[p[1:f.r]] .= x[p[1:f.r]]
  y[p[f.r+1:end]] .= 0.0
  return y, 0.0
end


text:
module JVODE

using LinearAlgebra

using Parameters: @unpack
using DiffEqBase: ODEProblem, @..

using Reexport: @reexport
@reexport using DiffEqBase
export Adams, BDF

###
### Types
###
abstract type JVNLSolver <: DiffEqBase.AbstractNLSolverAlgorithm
end
struct JVFunctional <: JVNLSolver
end
struct JVNewton <: JVNLSolver
end

abstract type JVODEAlgorithm <: DiffEqBase.AbstractODEAlgorithm
end
struct Adams{NLAlg} <: JVODEAlgorithm
    nlsolve::NLAlg
end
Adams() = Adams(JVFunctional())
struct BDF{NLAlg} <: JVODEAlgorithm
    nlsolve::NLAlg
end
BDF() = BDF(JVNewton())

# TODO: Optional Inputs and Outputs (in `cvode.h`)

# Basic JVODE constants
const ADAMS_Q_MAX = 12            # max value of q for lmm == Adams
const BDF_Q_MAX   =  5            # max value of q for lmm == BDF
const Q_MAX       =  ADAMS_Q_MAX  # max value of q for either lmm
const L_MAX       =  (Q_MAX+1)    # max value of L for either lmm
const NUM_TESTS   =  5            # number of error test quantities

mutable struct JVOptions{Rtol,Atol}
    reltol::Rtol
    abstol::Atol
end

mutable struct JVIntegrator{Alg,uType,tType,uEltype,solType,Rtol,Atol,F,P} <: DiffEqBase.AbstractODEIntegrator{Alg,true,uType,tType}
    sol::solType
    opts::JVOptions{Atol,Rtol}
    zn::NTuple{L_MAX,uType} # Nordsieck history array

    # vectors with length `length(u0)`
    ewt::uType # error weight vector
    u::uType
    acor::uType
    tempv::uType
    ftemp::uType

    # step data
    q::Int                        # current order
    qprime::Int                   # order to be used on the next step {q-1, q, q+1}
    qwait::Int                    # number of steps to wait before order change
    #L::Int                        # L = q+1
    dt::tType                     # current step size
    dtprime::tType                # next step size
    eta::tType                    # eta = dtprime / dt
    dtscale::tType                # the step size information in `zn`
    t::tType                      # current time
    tau::NTuple{L_MAX,tType}      # tuple of previous `q+1` successful step sizes
    tq::NTuple{NUM_TESTS,tType}   # tuple of test quantities
    coeff::NTuple{L_MAX,uEltype}  # coefficients of l(x)
    rl2::uEltype                  # 1/l[2]
    gamma::uEltype                # gamma = h * rl2
    gammap::uEltype               # `gamma` at the last setup call
    gamrat::uEltype               # gamma/gammap
    crate::uEltype                # estimated corrector convergence rate
    acnrm::uEltype                # | acor | wrms
    mnewt::Int                    # Newton iteration counter

    # Limits
    qmax::Int                     # q <= qmax
    mxstep::Int                   # maximum number of internal steps for one user call
    maxcor::Int                   # maximum number of `nlsolve`
    mxhnil::Int                   # maximum number of warning messages issued to the
                                  # user that `t + h == t` for the next internal step

    dtmin::tType                  # |h| >= hmin
    dtmax_inv::tType              # |h| <= 1/hmax_inv
    etamax::tType                 # eta <= etamax

    # counters
    nst::Int                      # number of internal steps taken
    nfe::Int                      # number of f calls
    ncfn::Int                     # number of corrector convergence failures
    netf::Int                     # number of error test failures
    nni::Int                      # number of Newton iterations performed
    nsetups::Int                  # number of setup calls
    nhnil::Int                    # number of messages issued to the user that
                                  # `t + h == t` for the next iternal step
    lrw::Int                      # number of real words in CVODE work vectors
    liw::Int                      # no. of integer words in CVODE work vectors

    # saved vales
    qu::Int                       # last successful q value used
    nstlp::Int                    # step number of last setup call
    dtu::tType                    # last successful h value used
    saved_tq5::tType              # saved value of tq[5]
    jcur::Bool                    # Is the Jacobian info used by
                                  # linear solver current?
    tolsf::Float64                # tolerance scale factor
    setupNonNull::Bool            # Does setup do something?

    # Arrays for Optional Input and Optional Output

    #long int *cv_iopt::Int  # long int optional input, output */
    #real     *cv_ropt::Int  # real optional input, output     */
    function JVIntegrator(prob::ODEProblem, ::solType, ::Alg, ::Rtol, ::Atol) where {solType,Alg,Rtol,Atol}
        @unpack f, u0, tspan, p = prob
        obj = new{Alg,typeof(u0),eltype(tspan),eltype(u0),solType,Rtol,Atol,typeof(f),typeof(p)}()
        return obj
    end
end

###
### Default Constants
###
const HMIN_DEFAULT     = 0.0
const HMAX_INV_DEFAULT = 0.0
const MXHNIL_DEFAULT   = 10
const MXSTEP_DEFAULT   = 5000

###
### Routine-Specific Constants
###
macro define(name, val)
    esc(:(const $name = $val))
end
begin

    # CVodeDky */

    @define FUZZ_FACTOR 100.0

    # CVHin */

    @define HLB_FACTOR 100.0
    @define HUB_FACTOR 0.1
    @define H_BIAS     0.5
    @define MAX_ITERS  4

    # CVSet */

    @define CORTES 0.1

    # CVStep return values */

    @define SUCCESS_STEP      0
    @define REP_ERR_FAIL     -1
    @define REP_CONV_FAIL    -2
    @define SETUP_FAILED     -3
    @define SOLVE_FAILED     -4

    # CVStep control constants */

    @define PREDICT_AGAIN    -5
    @define DO_ERROR_TEST     1

    # CVStep */

    @define THRESH 1.5
    @define ETAMX1 10000.0
    @define ETAMX2 10.0
    @define ETAMX3 10.0
    @define ETAMXF 0.2
    @define ETAMIN 0.1
    @define ETACF  0.25
    @define ADDON  0.000001
    @define BIAS1  6.0
    @define BIAS2  6.0
    @define BIAS3  10.0
    @define ONEPSM 1.000001

    @define SMALL_NST    10   # nst > SMALL_NST => use ETAMX3          */
    @define MXNCF        10   # max no. of convergence failures during */
    # one step try                           */
    @define MXNEF         7   # max no. of error test failures during  */
    # one step try                           */
    @define MXNEF1        3   # max no. of error test failures before  */
    # forcing a reduction of order           */
    @define SMALL_NEF     2   # if an error failure occurs and         */
    # SMALL_NEF <= nef <= MXNEF1, then       */
    # reset eta =  MIN(eta, ETAMXF)          */
    @define LONG_WAIT    10   # number of steps to wait before         */
    # considering an order change when       */
    # q==1 and MXNEF1 error test failures    */
    # have occurred                          */

    # CVnls return values */

    @define SOLVED            0
    @define CONV_FAIL        -1
    @define SETUP_FAIL_UNREC -2
    @define SOLVE_FAIL_UNREC -3

    # CVnls input flags */

    @define FIRST_CALL      0
    @define PREV_CONV_FAIL -1
    @define PREV_ERR_FAIL  -2

    # CVnls other constants */

    @define FUNC_MAXCOR 3  # maximum no. of corrector iterations   */
    # for iter == FUNCTIONAL                */
    @define NEWT_MAXCOR 3  # maximum no. of corrector iterations   */
    # for iter == NEWTON                    */

    @define CRDOWN 0.3 # constant used in the estimation of the   */
    # convergence rate (crate) of the          */
    # iterates for the nonlinear equation      */
    @define DGMAX  0.3 # iter == NEWTON, |gamma/gammap-1| > DGMAX */
    # => call lsetup                           */

    @define RDIV        2  # declare divergence if ratio del/delp > RDIV  */
    @define MSBP       20  # max no. of steps between lsetup calls        */

    @define TRY_AGAIN  99  # control constant for CVnlsNewton - should be */
    # distinct from CVnls return values            */
end

###
### CVODE Implementation
###
function DiffEqBase.__init(prob::ODEProblem, alg::JVODEAlgorithm; reltol=1e-3, abstol=1e-6)
    # analogous to `CVodeMalloc`
    @unpack f, u0, tspan, p = prob
    # copy input paramters
    uType = typeof(u0)
    uEltype = eltype(u0)
    tType = eltype(tspan)
    timeseries = uType[]
    ts = tType[]

    destats = DiffEqBase.DEStats(0)
    # TODO: dense output
    dense = true
    #ks = Vector{uType}(undef, 0)
    #id = InterpolationData(f,timeseries,ts,ks,dense,cache)

    sol = DiffEqBase.build_solution(prob,alg,ts,timeseries,
                                    dense=dense, #k=ks, interp=id, TODO
                                    calculate_error = false, destats=destats)
    integrator = JVIntegrator(prob, sol, alg, reltol, abstol)

    integrator.sol = sol
    integrator.t = prob.tspan |> first
    integrator.opts = JVOptions(reltol, abstol)
    integrator.dtmin = HMIN_DEFAULT
    integrator.dtmax_inv = HMAX_INV_DEFAULT
    integrator.mxhnil = MXHNIL_DEFAULT
    integrator.mxstep = MXSTEP_DEFAULT
    integrator.maxcor = alg.nlsolve isa JVNewton ? NEWT_MAXCOR : FUNC_MAXCOR

    maxorder = alg isa Adams ? ADAMS_Q_MAX : BDF_Q_MAX

    # allocate the vectors
    integrator.u = copy(u0)
    integrator.zn = ntuple(i->i<=maxorder+1 ? similar(integrator.u) : similar(integrator.u, 0), Val(L_MAX))
    integrator.ewt = similar(integrator.u)
    integrator.acor = similar(integrator.u)
    integrator.tempv = similar(integrator.u)
    integrator.ftemp = similar(integrator.u)

    # set the `ewt` vector
    setewt!(integrator, integrator.u)

    # set step paramters
    integrator.q = 1
    integrator.qwait = integrator.q + 1
    integrator.etamax = maxorder

    # TODO: set the linear solver

    # initialize `zn[1]` in the history array
    copyto!(integrator.zn[1], u0)

    # initialize all counters
    integrator.nst = integrator.nfe = integrator.ncfn = integrator.netf = integrator.nni = integrator.nsetups = integrator.nhnil = integrator.lrw = integrator.liw = 0

    # initialize misc
    integrator.qu = 0
    integrator.dtu = zero(integrator.dtu)
    integrator.tolsf = 1
    return integrator
end

DiffEqBase.has_reinit(integrator::JVIntegrator) = true
function DiffEqBase.reinit!(integrator::JVIntegrator, u0=integrator.uprev;
                            alg=integrator.alg, tspan=integrator.tspan,
                            reltol=integrator.reltol, abstol=integrator.abstol,
                            p=integrator.p) # TODO: opts
    # TODO: `CVReInit`
end

function DiffEqBase.step!(integrator::JVIntegrator, dt=nothing)#, stop_at_tdt=false)
    @unpack f, p, tspan = integrator.sol.prob
    @unpack zn, t, dtmax_inv, dtmin, nst, mxstep, mxhnil = integrator
    tout = dt === nothing ? tspan[end] : t + dt
    if nst === 0 # first step
        f(zn[2], zn[1], p, t)
        integrator.nfe = 1
        integrator.dt = zero(integrator.dt)
        # TODO: user initdt
        if iszero(integrator.dt)
            dt_ok = initdt!(integrator, tout)
            dt_ok || error("tout=$tout too close to t0=$t to start integration")
        end
        #rh = abs(integrator.dt) * dtmax_inv
        #rh > one(rh) && (integrator.dt /= rh)
        abs(integrator.dt) < dtmin && (integrator.dt *= dtmin / abs(integrator.dt))
        integrator.dtscale = integrator.dt
        lmul!(integrator.dt, zn[2])
    end
    # If not the first call, check if tout already reached
    if dt !== nothing && nst > 0 && (integrator.t - tout)*integrator.dt >= zero(integrator.dt)
        # integrator.t = tout
        # interpolate???
        return integrator
    end

    nstloc = 0
    while true
        nextdt = integrator.dt
        nextq  = integrator.q

        # reset and check ewt
        nst > 0 && setewt!(integrator, zn[1])

        # check for too many steps
        if nstloc >= mxstep
            @warn "At t=$t, mxstep=$mxstep steps taken on this call before reaching tout=$tout"
            copyto!(integrator.u, integrator.zn[1])
            break
        end

        # TODO: check for too much accuracy requested

        # check for `dt` below roundoff
        t = integrator.t
        if t + integrator.dt == t
            integrator.nhnil += 1
            integrator.nhnil <= mxhnil && @warn "internal t=$t and step size dt=$(integrator.dt) are such that t + dt = t on the next step. The solver will continue anyway."
            integrator.nhnil == mxhnil && @warn "The above warning has been issued $mxhnil times and will not be issued again for this problem."
        end

        kflag = jvstep(integrator)
        t = integrator.t
        if kflag !== SUCCESS_STEP
            handlefailure(integrator, kflag)
            copyto!(integrator.u, zn[1])
            break
        end

        nstloc += 1

        if dt === nothing
            copyto!(integrator.u, zn[1])
            nextq  = integrator.qprime
            nextdt = integrator.dtprime
        end

        if (t-tout)*integrator.dt >= zero(integrator.t)
            copyto!(integrator.u, zn[1])
            nextq  = integrator.qprime
            nextdt = integrator.dtprime
            break
        end
    end
    return integrator
end

function jvstep(integrator::JVIntegrator)
    # TODO
    return SUCCESS_STEP
end

function initdt!(integrator::JVIntegrator, tout)::Bool
    @unpack zn, t = integrator
    # test for tout too close to t
    tdiff = tout - t
    iszero(tdiff) && return false
    tdir = sign(tdiff)
    tdist = abs(tdiff)
    tround = eps() * max(abs(t), abs(tout))
    tdist < 2tround && return false

    # Set lower and upper bounds on h0, and take geometric mean. Exit with this
    # value if the bounds cross each other.
    hlb = HLB_FACTOR * tround
    hub = let integrator=integrator, tdist=tdist
        temp1 = integrator.tempv
        temp2 = integrator.acor
        map!(abs, temp1, zn[1])
        map!(abs, temp2, zn[2])
        @.. temp1 = temp2 / muladd(HUB_FACTOR, temp1, integrator.opts.abstol)
        hub_inv = norm(temp1, Inf)
        hub = HUB_FACTOR * tdist
        hub*hub_inv > one(hub) && (hub = inv(hub_inv))
        hub
    end
    hg = sqrt(hlb*hub)
    if hub < hlb
        tdir == -1 && (hg = -hg)
        integrator.dt = hg
        return true
    end

    # Loop up to MAX_ITERS times to find h0.
    # Stop if new and previous values differ by a factor < 2.
    # Stop if hnew/hg > 2 after one iteration, as this probably means
    # that the ydd value is bad because of cancellation error.
    count = 0
    yddnorm = (integrator, hg) -> begin
        # This routine computes an estimate of the second derivative of y
        # using a difference quotient, and returns its WRMS norm.
        @unpack zn, tempv, ewt, u, t = integrator
        @unpack f, p = integrator.sol.prob
        @.. u = muladd(hg, zn[2], zn[1])
        f(tempv, u, p, t+hg)
        integrator.nfe += 1
        @. tempv = tempv - zn[2]
        lmul!(inv(hg), tempv)
        yddnrm = wrmsnorm(tempv, ewt)
    end
    local hnew
    while true
        hgs = hg*tdir
        yddnrm = yddnorm(integrator, hgs)
        hnew = (yddnrm*hub*hub > 2) ? sqrt(2/yddnrm) : sqrt(hg*hub)
        count += 1
        count >= MAX_ITERS && break
        hrat = hnew/hg
        (hrat > 0.5) && (hrat < 2) && break
        if (count >= 2) && (hrat > 2)
            hnew = hg
            break
        end
        hg = hnew
    end

    # Apply bounds, bias factor, and attach sign */

    h0 = H_BIAS*hnew
    h0 < hlb && (h0 = hlb)
    h0 > hub && (h0 = hub)
    tdir == -1 && (h0 = -h0)
    integrator.dt = h0
    return true
end

setewt!(integrator::JVIntegrator, u) = (@.. integrator.ewt = inv(integrator.opts.reltol * abs(u) + integrator.opts.abstol); return)
wrmsnorm(x, w) = sum(((x,w),)->abs2(x*w), zip(x, w))/length(x) |> sqrt

"""
    jvode(f, u0, tspan; rtol=1e-3, atol=1e-6)

DVODE: Variable-coefficient Ordinary Differential Equation solver,
with fixed-leading-coefficient implementation.
This version is in double precision.

DVODE solves the initial value problem for stiff or nonstiff
systems of first order ODEs,
    dy/dt = f(t,y) ,  or, in component form,
    dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(NEQ)) (i = 1,...,NEQ).
DVODE is a package based on the EPISODE and EPISODEB packages, and
on the ODEPACK user interface standard, with minor modifications.
----------------------------------------------------------------------
Authors:
              Peter N. Brown and Alan C. Hindmarsh
              Center for Applied Scientific Computing, L-561
              Lawrence Livermore National Laboratory
              Livermore, CA 94551
and
              George D. Byrne
              Illinois Institute of Technology
              Chicago, IL 60616
----------------------------------------------------------------------
References:

1. P. N. Brown, G. D. Byrne, and A. C. Hindmarsh, "VODE: A Variable
   Coefficient ODE Solver," SIAM J. Sci. Stat. Comput., 10 (1989),
   pp. 1038-1051.  Also, LLNL Report UCRL-98412, June 1988.
2. G. D. Byrne and A. C. Hindmarsh, "A Polyalgorithm for the
   Numerical Solution of Ordinary Differential Equations,"
   ACM Trans. Math. Software, 1 (1975), pp. 71-96.
3. A. C. Hindmarsh and G. D. Byrne, "EPISODE: An Effective Package
   for the Integration of Systems of Ordinary Differential
   Equations," LLNL Report UCID-30112, Rev. 1, April 1977.
4. G. D. Byrne and A. C. Hindmarsh, "EPISODEB: An Experimental
   Package for the Integration of Systems of Ordinary Differential
   Equations with Banded Jacobians," LLNL Report UCID-30132, April
   1976.
5. A. C. Hindmarsh, "ODEPACK, a Systematized Collection of ODE
   Solvers," in Scientific Computing, R. S. Stepleman et al., eds.,
   North-Holland, Amsterdam, 1983, pp. 55-64.
6. K. R. Jackson and R. Sacks-Davis, "An Alternative Implementation
   of Variable Step-Size Multistep Formulas for Stiff ODEs," ACM
   Trans. Math. Software, 6 (1980), pp. 295-318.
"""
function jvode(f, u0, tspan; rtol=1e-3, atol=1e-6)
    maxorder = (12, 5)
    maxstep0 = 500
    maxwarns = 10
end

end # module


text:
"""
    Metrics(metrics...)

Callback that tracks metrics during training.

`metrics` can be both [`AbstractMetric`](#)s or functions
like `f(ŷs, ys)` which will be converted to [`Metric`](#)s.

A metric tracking `lossfn` is included by default.

## Examples

- `metrics = Metrics(accuracy)`
- `metrics = Metrics(Metric(Flux.mse, device = gpu), Metric(Flux.mae, device = gpu))`

"""
struct Metrics <: Callback
    metrics::Tuple
    function Metrics(metrics...)
        return new(Tuple(m isa AbstractMetric ? m : Metric(m) for m in (Loss(), metrics...)))
    end
end

runafter(::Metrics) = (Recorder,)
stateaccess(::Metrics) = (
    cbstate = (metricsstep = Write(), metricsepoch = Write(), history = Read()),
    batch = Read()
)

Base.show(io::IO, metrics::Metrics) = print(io, "Metrics(", join(string.(metrics.metrics), ", "), ")")

# store metrics in `cbstate` so other callbacks can access them
function on(::Init, ::Phase, metrics::Metrics, learner)
    learner.cbstate.metricsstep = MVHistory()
    learner.cbstate.metricsepoch = DefaultDict{Phase, MVHistory}(() -> MVHistory())
end

on(::EpochBegin, ::Phase, metrics::Metrics, learner) = foreach(reset!, metrics.metrics)

function on(::BatchEnd, phase, metrics::Metrics, learner)
    metricsstep = learner.cbstate.metricsstep
    step = learner.cbstate.history.steps
    for metric in metrics.metrics
        step!(metric, learner)
        if phase isa AbstractTrainingPhase
            push!(metricsstep, Symbol(metricname(metric)), step, stepvalue(metric))
        end
    end
end

function on(::EpochEnd, phase, metrics::Metrics, learner)
    metricsepoch = learner.cbstate.metricsepoch
    epoch = learner.cbstate.history.epochs
    for metric in metrics.metrics
        push!(metricsepoch[phase], Symbol(metricname(metric)), epoch, epochvalue(metric))
    end
end

# AbstractMetric interface

"""
    abstract type AbstractMetric

Abstract type for metrics passed to [`Metrics`](#).

For most use cases, you should use [`Metric`](#), the standard
implementation.

## Interface

If [`Metric`](#) doesn't fit your use case, you can create
a new subtype of `AbstractMetric` and implement the following
methods to make it compatible with [`Metrics`](#):

- [`reset!`](#)`(metric)`
- [`step!`](#)`(metric, learner)`
- [`stepvalue`](#)`(metric)`
- [`epochvalue`](#)`(metric)`
- [`metricname`](#)`(metric)`
"""
abstract type AbstractMetric end

mutable struct Metric{T}
    metricfn
    statistic::OnlineStat{T}
    _statistic
    name
    device
    last::Union{Nothing, T}
end

Base.show(io::IO, metric::Metric{T}) where T = print(io, "Metric(", metric.name, ")")

"""
    Metric(metricfn[; statistic, device, name])

Implementation of [`AbstractMetric`](#) that can be used with the
[`Metrics`](#) callback.

## Arguments

- `metricfn(ŷs, ys)` should return a number.
- `statistic` is a `OnlineStats.Statistic` that is updated after every step.
    The default is `OnlineStats.Mean()`
- `name` is used for printing.
- `device = cpu` is a function applied to `ys` and `ŷs` before calling
    `metricfn`. If `metricfn` works on the GPU, you may want to pass
    `Flux.gpu` here for better performance.

## Examples

- `Metric(accuracy)`
- `Metric(Flux.mse, device = gpu)`
- `Metric(Flux.mae, device = gpu)`

"""
function Metric(
        metricfn;
        name = uppercasefirst(string(metricfn)),
        statistic = Mean(),
        device = cpu)

    return Metric(
        metricfn,
        deepcopy(statistic),
        statistic,
        name,
        device,
        nothing,
    )
end

function reset!(metric::Metric)
    metric.statistic = deepcopy(metric._statistic)
end

function step!(metric::Metric, learner)
    ŷs, ys = learner.batch.ŷs, learner.batch.ys
    metric.last = metric.metricfn(ŷs, ys)
    OnlineStats.fit!(metric.statistic, metric.last)
end

stepvalue(metric::Metric) = metric.last
epochvalue(metric::Metric) = OnlineStats.value(metric.statistic)
metricname(metric::Metric) = metric.name

# Loss Metric

mutable struct Loss <: AbstractMetric
    sum
    last
    count
end

Loss() = Loss(0, Inf, 0)
Base.show(io::IO, loss::Loss) = print(io, "Loss()")

function reset!(metric::Loss)
    metric.last = nothing
    metric.sum = 0.
    metric.count = 0
end

function step!(metric::Loss, learner)
    metric.last = learner.batch.loss
    metric.sum += learner.batch.loss
    metric.count += 1
end

stepvalue(metric::Loss) = metric.last
epochvalue(metric::Loss) = metric.sum / metric.count
metricname(metric::Loss) = "Loss"


text:
@testset "137.single-number-ii.jl" begin
    @test single_number([2, 2, 3, 2]) == 3
    @test single_number([0, 1, 0, 1, 0, 1, 99]) == 99
end


text:
@testset "default_material" begin
    m = RayTracer.material()
    @test m.color ≈ RayTracer.color(1,1,1)
    @test m.ambient ≈ 0.1
    @test m.diffuse ≈ 0.9
    @test m.specular ≈ 0.9
    @test m.shininess ≈ 200
end

@testset "lighning_tests" begin
    m = RayTracer.material()
    position = RayTracer.point(0,0,0)

    @testset "lighning_with_eye_between_lignt_and_surface" begin
        eyev = RayTracer.vector(0,0,-1)
        normalv = RayTracer.vector(0,0,-1)
        light = RayTracer.point_light(RayTracer.point(0,0,-10),RayTracer.color(1,1,1))
        result = RayTracer.lighting(m,light,position,eyev,normalv)
        @test result ≈ RayTracer.color(1.9,1.9,1.9)
    end

    @testset "lighning_with_eye_between_lignt_and_surface_offset_45_degrees" begin
        eyev = RayTracer.vector(0,√2/2,-√2/2)
        normalv = RayTracer.vector(0,0,-1)
        light = RayTracer.point_light(RayTracer.point(0,0,-10),RayTracer.color(1,1,1))
        result = RayTracer.lighting(m,light,position,eyev,normalv)
        @test result ≈ RayTracer.color(1,1,1)
    end

    @testset "lighning_with_eye_opposite_surface_light_offset_45_degrees" begin
        eyev = RayTracer.vector(0,0,-1)
        normalv = RayTracer.vector(0,0,-1)
        light = RayTracer.point_light(RayTracer.point(0,10,-10),RayTracer.color(1,1,1))
        result = RayTracer.lighting(m,light,position,eyev,normalv)
        @test result ≈ RayTracer.color(0.1+0.9*√2/2,0.1+0.9*√2/2,0.1+0.9*√2/2)
    end

    @testset "lighning_with_eye_path_reflection_vector" begin
        eyev = RayTracer.vector(0,-√2/2,-√2/2)
        normalv = RayTracer.vector(0,0,-1)
        light = RayTracer.point_light(RayTracer.point(0,10,-10),RayTracer.color(1,1,1))
        result = RayTracer.lighting(m,light,position,eyev,normalv)
        @test result ≈ RayTracer.color(0.1+0.9*√2/2+0.9,0.1+0.9*√2/2+0.9,0.1+0.9*√2/2+0.9)
    end

    @testset "lighning_with_light_behind_surface" begin
        eyev = RayTracer.vector(0,0,-1)
        normalv = RayTracer.vector(0,0,-1)
        light = RayTracer.point_light(RayTracer.point(0,0,10),RayTracer.color(1,1,1))
        result = RayTracer.lighting(m,light,position,eyev,normalv)
        @test result ≈ RayTracer.color(0.1,0.1,0.1)
    end

end


text:
using Distributed
using ClusterManagers
using JLD2 # extra dependency
using Distributions

addprocs(32, exeflags = "--project")

@everywhere using SurvivalSignature, LinearAlgebra

@everywhere joinpath(pwd(), "demo", "models", "reduced_gb.jl") |> include

@everywhere function floyd_warshall!(D)

    n = size(D, 1)

    @inbounds for k = 1:n, i = 1:n, j = 1:n
        sum = D[i, k] + D[k, j]
        (sum < D[i, j]) && (D[i, j] = sum)
    end

    return D
end

@everywhere function efficiency(D)
    n = size(D, 1)

    floyd_warshall!(D)
    D[D.==0] .= Inf

    D = D .^ -1

    D[diagind(D)] == 0

    1 / (n * (n - 1)) * sum(D)
end

@everywhere E = efficiency(copy(adj))

@everywhere function φ(system::Array{Float64,2}, x::Vector)

    A = copy(system)

    failed = setdiff(nodes, x)

    @inbounds A[failed, :] .= Inf
    @inbounds A[:, failed] .= Inf

    efficiency(A) / E > 0.5
end

@load "demo/data/reduced_gb_exact.jld2"

samples = [1e1, 1e2, 1e3, 1e4, 1e5, 1e6]
n = 1000
absolute_error = zeros(length(samples), n)
relative_error = zeros(length(samples), n)
mse = zeros(length(samples), n)

for (i, s) ∈ enumerate(samples)
    for j = 1:n
        signature, _ =
            survivalsignature(adj, types, φ, Int(s), 0.001, percolation_preprocessor)

        absolute_error[i, j] = norm(Φ - signature)
        relative_error[i, j] = norm(Φ - signature) / norm(Φ)
        mse[i, j] = mean((Φ - signature) .^ 2)
    end
end

@save "reduced_gb_errors.jld2" samples absolute_error relative_error mse


text:
using Amincoacid_to_RNA_Translation
using Test


text:
"""
`optimalPLD_solveSAO_constraints( n, revenues, numPermutations, lambda, orderings, MIP_timeLimit, isRelaxation, x_initial, A_feas, b_feas )`

Formulate and solve the first-choice PLD formulation from Bertsimas and Mišić (see Section 3.2 of that paper).

`n` = number of candidate products.

`revenues` = array of n numbers representing marginal revenue/profit of each of the n products. (These are the pi_i's in the Bertsimas and Mišić paper.)

`numPermutations` = number of rankings / customer types.

`lambda` = array specifying the probability of each ranking / customer type; should be nonnegative and sum to 1.

`orderings` = 2D array of integers encoding each ranking. orderings[k,j] contains the index of ranking k's jth top product; (n+1) encodes the no-purchase option. For example, if n = 5 and

	orderings[k,:] = [3,5,2,6,1,4]

then the customer prefers 3 < 5 < 2 < 0 (no-purchase option).

`MIP_timeLimit` = maximum time to run Gurobi for. If set to 0, then no time limit is set.

`isRelaxation` = boolean to indicate whether we are solving LP relaxation (true) or the integer problem (false).

`x_initial` = array with with n values (either 0 or 1) indicating the initial values of the x_i variables to use to warm-start Gurobi. If empty, no warm-start solution will be given to Gurobi. Only allowed for integer problem (isRelaxation = false).

`A_feas`, `b_feas` = 2D array and 1D array specifying the system of linear inequalities Ax <= b for the PLD problem. This can be used to impose side constraints (this is the Cx <= d constraint system in the formulation in Bertsimas and Mišić 2018). If b_feas is empty, then no constraints will be added (i.e., it will solve the unconstrained problem).

Returns:

`x_val`, `y_val` = arrays of optimal x_i and y^k_i variables for either the integer problem or the relaxation.

`expRevenue` = expected per-customer revenue/profit of the best integer solution; this is the optimal objective if the problem is solved to full optimality.

`upperbound` = best upper bound on the optimal expected per customer revenue/profit found at termination; equal to the optimal objective (modulo Gurobi's default tolerance) if the problem is solved to full optimality.

`solutiontime` = time required to solve the problem.

`MIPgap` = optimality gap upon termination.

"""
function optimalPLD_solveSAO_constraints( n, revenues, numPermutations, lambda, orderings, MIP_timeLimit, isRelaxation, x_initial, A_feas, b_feas )

	if (MIP_timeLimit > 0)
		m = Model( solver = GurobiSolver(TimeLimit = MIP_timeLimit))
	else
		m = Model( solver = GurobiSolver())
	end

	#weights = rand(Gamma())

	if (!isRelaxation)
		@variable(m, x[1:n], Bin)
	else
		@variable(m, 0 <= x[1:n] <= 1)
	end

	if (!isRelaxation)
		if (!isempty(x_initial))
			for i=1:n
				setvalue(x[i], x_initial[i]);
			end
		end
	end

	@variable(m, y[1:(n+1), 1:numPermutations] >= 0 )
	#@variable(m, y[1:(n+1), 1:numSamples], Bin )

	@objective(m, Max, sum( lambda[k] * revenues[i] * y[i,k] for i=1:n, k=1:numPermutations)  )

	# CONSTRAINTS
	for t=1:length(b_feas)
		@constraint(m, sum(A_feas[t,i]*x[i] for i=1:n) <= b_feas[t])
	end

	for k=1:numPermutations
		for i=1:n
			@constraint(m, y[i,k] <= x[i])
		end
		for i2=1:n
			if (orderings[k,i2] <= n)
				@constraint(m, sum( y[orderings[k,i3],k] for i3=(i2+1):(n+1)) <= 1 - x[orderings[k,i2]] )
			else
				# VVM: this should be correct, but let's wait to see...
				@constraint(m, sum( y[orderings[k,i3],k] for i3=(i2+1):(n+1)) <= 0)
			end
		end
		@constraint(m, sum( y[i,k] for i=1:(n+1)) == 1)
	end

	solve(m)

	println(getvalue(x))

	x_val = getvalue(x)[:]
	y_val = getvalue(y)[:,:]

	solutiontime = Gurobi.get_runtime(m.internalModel.inner);
	if (!isRelaxation)
		upperBound = Gurobi.get_dblattr(m.internalModel.inner, "ObjBound");
		MIPgap = Gurobi.get_dblattr(m.internalModel.inner, "MIPGap");
	else
		upperBound = -1;
		MIPgap = -1;
	end

	expRevenue = getobjectivevalue(m)

	return x_val, y_val, expRevenue, upperBound, solutiontime, MIPgap
end

max_stars_repo_name: vvmisic/EFCPLO

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: vvmisic/EFCPLO

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: vvmisic/EFCPLO

max_forks_repo_licenses: MIT Unlicense


text:
"""
    MultiplicationOperation <: AbstractGroupOperation

Group operation that consists of multiplication.
"""
struct MultiplicationOperation <: AbstractGroupOperation end

const MultiplicationGroupTrait = TraitList{<:IsGroupManifold{<:MultiplicationOperation}}

Base.:*(e::Identity{MultiplicationOperation}) = e
Base.:*(::Identity{MultiplicationOperation}, p) = p
Base.:*(p, ::Identity{MultiplicationOperation}) = p
Base.:*(e::Identity{MultiplicationOperation}, ::Identity{MultiplicationOperation}) = e
Base.:*(::Identity{MultiplicationOperation}, e::Identity{AdditionOperation}) = e
Base.:*(e::Identity{AdditionOperation}, ::Identity{MultiplicationOperation}) = e

Base.:/(p, ::Identity{MultiplicationOperation}) = p
Base.:/(::Identity{MultiplicationOperation}, p) = inv(p)
Base.:/(e::Identity{MultiplicationOperation}, ::Identity{MultiplicationOperation}) = e

Base.:\(p, ::Identity{MultiplicationOperation}) = inv(p)
Base.:\(::Identity{MultiplicationOperation}, p) = p
Base.:\(e::Identity{MultiplicationOperation}, ::Identity{MultiplicationOperation}) = e

LinearAlgebra.det(::Identity{MultiplicationOperation}) = true
LinearAlgebra.adjoint(e::Identity{MultiplicationOperation}) = e

function identity_element!(
    ::MultiplicationGroupTrait,
    G::AbstractDecoratorManifold,
    p::AbstractMatrix,
)
    return copyto!(p, I)
end

function identity_element!(
    ::MultiplicationGroupTrait,
    G::AbstractDecoratorManifold,
    p::AbstractArray,
)
    if length(p) == 1
        fill!(p, one(eltype(p)))
    else
        throw(DimensionMismatch("Array $p cannot be set to identity element of group $G"))
    end
    return p
end

function is_identity(
    ::MultiplicationGroupTrait,
    G::AbstractDecoratorManifold,
    q::Number;
    kwargs...,
)
    return isapprox(G, q, one(q); kwargs...)
end
function is_identity(
    ::MultiplicationGroupTrait,
    G::AbstractDecoratorManifold,
    q::AbstractArray{<:Any,0};
    kwargs...,
)
    return isapprox(G, q[], one(q[]); kwargs...)
end
function is_identity(
    ::MultiplicationGroupTrait,
    G::AbstractDecoratorManifold,
    q::AbstractMatrix;
    kwargs...,
)
    return isapprox(G, q, I; kwargs...)
end

LinearAlgebra.mul!(q, ::Identity{MultiplicationOperation}, p) = copyto!(q, p)
LinearAlgebra.mul!(q, p, ::Identity{MultiplicationOperation}) = copyto!(q, p)
function LinearAlgebra.mul!(
    q::AbstractMatrix,
    ::Identity{MultiplicationOperation},
    ::Identity{MultiplicationOperation},
)
    return copyto!(q, I)
end
function LinearAlgebra.mul!(
    q,
    ::Identity{MultiplicationOperation},
    ::Identity{MultiplicationOperation},
)
    return copyto!(q, one(q))
end
function LinearAlgebra.mul!(
    q::Identity{MultiplicationOperation},
    ::Identity{MultiplicationOperation},
    ::Identity{MultiplicationOperation},
)
    return q
end
Base.one(e::Identity{MultiplicationOperation}) = e

Base.inv(::MultiplicationGroupTrait, G::AbstractDecoratorManifold, p) = inv(p)
function Base.inv(
    ::MultiplicationGroupTrait,
    G::AbstractDecoratorManifold,
    e::Identity{MultiplicationOperation},
)
    return e
end

inv!(::MultiplicationGroupTrait, G::AbstractDecoratorManifold, q, p) = copyto!(q, inv(G, p))
function inv!(
    ::MultiplicationGroupTrait,
    G::AbstractDecoratorManifold,
    q,
    ::Identity{MultiplicationOperation},
)
    return identity_element!(G, q)
end
function inv!(
    ::MultiplicationGroupTrait,
    G::AbstractDecoratorManifold,
    q::Identity{MultiplicationOperation},
    e::Identity{MultiplicationOperation},
)
    return q
end

compose(::MultiplicationGroupTrait, G::AbstractDecoratorManifold, p, q) = p * q

function compose!(::MultiplicationGroupTrait, G::AbstractDecoratorManifold, x, p, q)
    return mul!_safe(x, p, q)
end

function inverse_translate(
    ::MultiplicationGroupTrait,
    G::AbstractDecoratorManifold,
    p,
    q,
    ::LeftAction,
)
    return p \ q
end
function inverse_translate(
    ::MultiplicationGroupTrait,
    G::AbstractDecoratorManifold,
    p,
    q,
    ::RightAction,
)
    return q / p
end

function inverse_translate!(
    ::MultiplicationGroupTrait,
    G::AbstractDecoratorManifold,
    x,
    p,
    q,
    conv::ActionDirection,
)
    return copyto!(x, inverse_translate(G, p, q, conv))
end

function exp_lie!(
    ::MultiplicationGroupTrait,
    G::AbstractDecoratorManifold,
    q,
    X::Union{Number,AbstractMatrix},
)
    copyto!(q, exp(X))
    return q
end

function log_lie!(
    ::MultiplicationGroupTrait,
    G::AbstractDecoratorManifold,
    X::AbstractMatrix,
    q::AbstractMatrix,
)
    return log_safe!(X, q)
end
function log_lie!(
    ::MultiplicationGroupTrait,
    G::AbstractDecoratorManifold,
    X,
    ::Identity{MultiplicationOperation},
)
    return zero_vector!(G, X, identity_element(G))
end

function lie_bracket(::MultiplicationGroupTrait, G::AbstractDecoratorManifold, X, Y)
    return mul!(X * Y, Y, X, -1, true)
end

function lie_bracket!(::MultiplicationGroupTrait, G::AbstractDecoratorManifold, Z, X, Y)
    mul!(Z, X, Y)
    mul!(Z, Y, X, -1, true)
    return Z
end


text:
# Ordinary Maximum Likelihood Estimation

### Type
#= struct SemML{
        I <: AbstractArray,
        T <: AbstractArray,
        U, V,
        W <: Union{Nothing, AbstractArray}} <: LossFunction
    imp_inv::I
    mult::T # is this type known?
    objective::U
    grad::V
    meandiff::W
end =#

struct SemML{
        INV <: AbstractArray,
        C <: Union{AbstractArray, Nothing},
        L <: Union{AbstractArray, Nothing},
        M <: Union{AbstractArray, Nothing},
        M2 <: Union{AbstractArray, Nothing},
        U,
        V} <: LossFunction
    inverses::INV #preallocated inverses of imp_cov
    choleskys::C #preallocated choleskys
    mult::M
    logdets::L #logdets of implied covmats
    meandiff::M2
    #data_rowwise::DP  -> add to semobs
    objective::U
    grad::V
end

### Constructor
function SemML(observed::T, objective, grad) where {T <: SemObs}
    isnothing(observed.obs_mean) ?
        meandiff = nothing :
        meandiff = copy(observed.obs_mean)
    return SemML(
        copy(observed.obs_cov),
        nothing,
        copy(observed.obs_cov),
        nothing,
        meandiff,
        copy(objective),
        copy(grad),
        )
end

### Loss
# generic (not optimized)
function (semml::SemML)(par, model)
    B = copy(model.imply.imp_cov)
    C = similar(model.imply.imp_cov)

    a = cholesky!(Hermitian(B); check = false)
    if !isposdef(a)
        F = Inf
    else
        a_inv = inv(a)
        mul!(C, a_inv, model.observed.obs_cov)
        F = logdet(a) +
             tr(C)
        if !isnothing(model.imply.imp_mean)
        meandiff = model.observed.obs_mean - model.imply.imp_mean
        F_mean = meandiff'*a_inv*meandiff
        F += F_mean
    end
    end
    return F
end

# maximum speed for finite differences
function (semml::SemML)(par, model::Sem{O, I, L, D}) where
            {O <: SemObs, L <: Loss, I <: Imply, D <: SemFiniteDiff}
    semml.inverses .= model.imply.imp_cov
    a = cholesky!(Symmetric(semml.inverses); check = false)
    if !isposdef(a)
        F = Inf
    else
        ld = logdet(a)
        semml.inverses .= LinearAlgebra.inv!(a)
        #inv_cov = inv(a)
        mul!(semml.mult, semml.inverses, model.observed.obs_cov)
        F = ld + tr(semml.mult)

        if !isnothing(model.imply.imp_mean)
            @. semml.meandiff = model.observed.obs_mean - model.imply.imp_mean
            F_mean = semml.meandiff'*semml.inverses*semml.meandiff
            F += F_mean
        end
    end
    return F
end

# maximum speed for analytic diff
function (semml::SemML)(par, model::Sem{O, I, L, D}) where
            {O <: SemObs, L <: Loss, I <: Imply, D <: SemAnalyticDiff}
    semml.inverses .= model.imply.imp_cov
    a = cholesky!(Symmetric(semml.inverses); check = false)
    if !isposdef(a)
        F = Inf
    else
        ld = logdet(a)
        semml.inverses .= LinearAlgebra.inv!(a)
        #inv_cov = inv(a)
        mul!(semml.mult, semml.inverses, model.observed.obs_cov)
        F = ld + tr(semml.mult)

        if !isnothing(model.imply.imp_mean)
            @. semml.meandiff = model.observed.obs_mean - model.imply.imp_mean
            F_mean = semml.meandiff'*semml.inverses*semml.meandiff
            F += F_mean
        end
    end
    return F
end

#= # analytic differentiation
function (semml::SemML)(par, model::Sem{O, I, L, D}, E, G) where
            {O <: SemObs, L <: Loss, I <: Imply, D <: SemAnalyticDiff}

    a = cholesky(Hermitian(model.imply.imp_cov); check = false)
    if !isposdef(a)
        if E != nothing
            semml.objective .= Inf
        end
        if G != nothing
            semml.grad .= 0.0
        end
    else
        ld = logdet(a)
        inv_cov = inv(a)
        if E != nothing
            mul!(semml.mult, inv_cov, model.observed.obs_cov)
            semml.objective .= ld + tr(semml.mult)
        end
        if G != nothing
            model.diff.B!(model.diff.B, par) # B = inv(I-A)
            model.diff.E!(model.diff.E, par) # E = B*S*B'
            let B = model.diff.B, E = model.diff.E,
                Σ_inv = inv_cov, F = model.diff.F,
                D = model.observed.obs_cov
                for i = 1:size(par, 1)
                    S_der = sparse(model.diff.S_ind_vec[i]..., model.diff.matsize...)
                    A_der = sparse(model.diff.A_ind_vec[i]..., model.diff.matsize...)

                    term = F*B*A_der*E*F'
                    Σ_der = Array(F*B*S_der*B'F' + term + term')

                    semml.grad[i] = tr(Σ_inv*Σ_der) + tr((-Σ_inv)*Σ_der*Σ_inv*D)
                end
            end
        end
    end
end
 =#


text:
"""
REINFORCE learner (Williams 1992) with baseline
(state values V)

See also REINFORCElearner.jl

Specifically for 2 actions and reward only at the end of episode!

Note
----
The policy parameters θ as well as the feature vector ϕ are
column vectors. Here we use a (states x actions) matrix for simpler access.
also, the feature vector ϕ is a one-hot vector in our case. therefore we
do not need it to keep a vector but access the entries by the
state/action pair [s_t, a_t]
"""

mutable struct REINFORCEbaseline
    nr_states::Int
    nr_actions::Int
    α::Float64
    α_eff::Float64
    γ::Float64
    γ_k::Float64
    action_selection_prob_in_episode::Array{Any,1}
    states_in_episode::Array{Any,1}
    actions_in_episode::Array{Any,1}
    θ::Array{Float64, 2}
    target::Array{Float64, 1}
    modulatingfactor_in_episode::Array{Any,1}
    V_baseline::Array{Float64, 1}
    α_baseline::Float64
    do_debug_print::Bool
end
function REINFORCEbaseline(; nr_states = 10, nr_actions = 4, α = 0.1, γ = .9, # initvalue = 0.,
    α_baseline = α, do_debug_print = false)
    θ = zeros(Float64, nr_states, nr_actions) #.+ initvalue,
    V_baseline = zeros(Float64, nr_states)

    REINFORCEbaseline(nr_states, nr_actions, α, α, γ, 1., [], [], [], θ, [0.], [],
                    V_baseline, α_baseline, do_debug_print)
end
export REINFORCEbaseline
function gettarget!(learner::REINFORCEbaseline, state, reward)
    V_s = reward * learner.γ_k  #  note: assuming a single reward, set during the forward pass above
    TD_error = V_s - learner.V_baseline[state]
    learner.V_baseline[state] += learner.α_baseline * TD_error
    learner.target[1] = copy(TD_error)
end


text:
import LaTeXStrings
import PGFPlotsX

"""
"""
function plotsinglelabelregressiontrueversuspredicted(
        estimator::AbstractFittable,
        features_df::DataFrames.AbstractDataFrame,
        labels_df::DataFrames.AbstractDataFrame,
        single_label_name::Symbol;
        includeorigin::Bool = false,
        legend_pos::AbstractString = "outer north east",
        )::PGFPlotsXPlot
    ytrue = singlelabelregressionytrue(
        labels_df[single_label_name],
        )
    predictionsalllabels = predict(estimator, features_df)
    ypred = singlelabelregressionypred(
        predictionsalllabels[single_label_name],
        )
    truevalueversuspredictedvalue_linearplotobject = PGFPlotsX.@pgf(
        PGFPlotsX.Plot(
            {
                only_marks,
                style = "black, fill = black",
                },
            PGFPlotsX.Coordinates(
                ypred,
                ytrue,
                ),
            )
        )
    if includeorigin
        perfectlinevalues = sort(
            unique(
                vcat(
                    0,
                    ytrue,
                    ),
                );
            rev = false,
            )
    else
        perfectlinevalues = sort(
            unique(
                ytrue,
                );
            rev = false,
            )
    end
    perfectline_linearplotobject = PGFPlotsX.@pgf(
        PGFPlotsX.Plot(
            {
                no_marks,
                style = "dotted, red, color=red, fill=red",
                },
            PGFPlotsX.Coordinates(
                perfectlinevalues,
                perfectlinevalues,
                ),
            ),
    )
    estimated_intercept,
        estimated_x_coefficient = simple_linear_regression(
            Float64.(ypred), # X
            Float64.(ytrue), # Y
            )
    bestfitline_linearplotobject = PGFPlotsX.@pgf(
        PGFPlotsX.Plot(
            {
                no_marks,
                style = "dashed, blue, color=blue, fill=blue",
                },
            PGFPlotsX.Coordinates(
                perfectlinevalues,
                estimated_intercept .+ estimated_x_coefficient*perfectlinevalues,
                ),
            ),
        )
    p = PGFPlotsX.@pgf(
        PGFPlotsX.Axis(
            {
                xlabel = LaTeXStrings.LaTeXString(
                    "Predicted value"
                    ),
                ylabel = LaTeXStrings.LaTeXString(
                    "True value"
                    ),
                legend_pos = legend_pos,
                },
            truevalueversuspredictedvalue_linearplotobject,
            perfectline_linearplotobject,
            bestfitline_linearplotobject,
            ),
        )
    wrapper = PGFPlotsXPlot(p)
    return wrapper
end


text:
module TestSlidingWindows

if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

using Bio.Windows,
    Bio.Seq,
    TestFunctions

function generate_testcase_data()
    n = rand(1:1000)
    arr = collect(1:n)
    str = randstring(n)
    nuc = DNASequence(random_dna(n))
    rib = RNASequence(random_rna(n))
    aas = AminoAcidSequence(random_aa(n))
    return n, arr, str, nuc, rib, aas
end

function generate_test_window(n::Int)
    return rand(1:n), rand(1:n)
end

function generate_test_itrs(width, step, arr, str, nuc, rib, aas)
    arritr = eachwindow(arr, width, step)
    stritr = eachwindow(str, width, step)
    nucitr = eachwindow(nuc, width, step)
    ribitr = eachwindow(rib, width, step)
    aasitr = eachwindow(aas, width, step)
    return arritr, stritr, nucitr, ribitr, aasitr
end

@testset "Sliding-Windows" begin
    @testset "Iterator Creation" begin
        @testset "Sensible values" begin
            for _ in 1:10
                n, arr, str, nuc, rib, aas = generate_testcase_data()
                winwidth, winstep = generate_test_window(n)
                arritr, stritr, nucitr, ribitr, aasitr =
                    generate_test_itrs(winwidth, winstep, arr, str, nuc, rib, aas)
                @test arritr == EachWindowIterator{Array{Int, 1}}(arr, winwidth, winstep)
                @test stritr == EachWindowIterator{ASCIIString}(str, winwidth, winstep)
                @test nucitr == EachWindowIterator{DNASequence}(nuc, winwidth, winstep)
                @test ribitr == EachWindowIterator{RNASequence}(rib, winwidth, winstep)
                @test aasitr == EachWindowIterator{AminoAcidSequence}(aas, winwidth, winstep)
                @test arritr.width == winwidth
                @test stritr.width == winwidth
                @test nucitr.width == winwidth
                @test ribitr.width == winwidth
                @test aasitr.width == winwidth
                @test arritr.step == winstep
                @test stritr.step == winstep
                @test nucitr.step == winstep
                @test ribitr.step == winstep
                @test aasitr.step == winstep
                @test arritr.data == arr
                @test stritr.data == str
                @test nucitr.data == nuc
                @test ribitr.data == rib
                @test aasitr.data == aas
            end
        end

        @testset "Bad values" begin
            for _ in 1:10
                n, arr, str, nuc = generate_testcase_data()
                winwidth = rand(n + 1:n + rand(1:1000))
                winstep = rand(n + 1:n + rand(1:1000))
                @test_throws Exception eachwindow(arr, winwidth, winstep)
                @test_throws Exception eachwindow(str, winwidth, winstep)
                @test_throws Exception eachwindow(nuc, winwidth, winstep)
                @test_throws Exception eachwindow(rib, winwidth, winstep)
                @test_throws Exception eachwindow(aas, winwidth, winstep)
            end
        end
    end

    @testset "Iteration" begin
        @testset "N & size of Windows, N of missed elements" begin
            for _ in 1:10
                n, arr, str, nuc, rib, aas = generate_testcase_data()
                winwidth, winstep = generate_test_window(n)
                arritr, stritr, nucitr, ribitr, aasitr =
                    generate_test_itrs(winwidth, winstep, arr, str, nuc, rib, aas)
                arrres = collect(arritr)
                strres = collect(stritr)
                nucres = collect(nucitr)
                ribres = collect(ribitr)
                aasres = collect(aasitr)
                expectedMissed = n - StepRange(winwidth, winstep, n).stop
                @test length(arrres) == size(arritr)
                @test length(strres) == size(stritr)
                @test length(nucres) == size(nucitr)
                @test length(ribres) == size(ribitr)
                @test length(aasres) == size(aasitr)
                @test missed(arritr) == expectedMissed
                @test missed(stritr) == expectedMissed
                @test missed(nucitr) == expectedMissed
                @test missed(ribitr) == expectedMissed
                @test missed(aasitr) == expectedMissed
                for win in arrres
                    @test length(win[2]) == winwidth
                end
                for win in strres
                    @test length(win[2]) == winwidth
                end
                for win in nucres
                    @test length(win[2]) == winwidth
                end
                for win in ribres
                    @test length(win[2]) == winwidth
                end
                for win in aasres
                    @test length(win[2]) == winwidth
                end
            end
        end

        @testset "Content of windows generated" begin
            @testset "DNA" begin
                nuc = DNASequence(random_dna(10))
                for width = 1:5, step = 1:5
                    nucitr = eachwindow(nuc, width, step)
                    i, j = 1, width
                    for win in nucitr
                        @test win[1] == i
                        @test win[2] == sub(nuc, i:j)
                        @test win[2] == nuc[i:j]
                        i += step
                        j += step
                    end
                end
            end

            @testset "RNA" begin
                rib = RNASequence(random_rna(10))
                for width = 1:5, step = 1:5
                    ribitr = eachwindow(rib, width, step)
                    i, j = 1, width
                    for win in ribitr
                        @test win[1] == i
                        @test win[2] == sub(rib, i:j)
                        @test win[2] == rib[i:j]
                        i += step
                        j += step
                    end
                end
            end

            @testset "Amino Acids" begin
                aas = AminoAcidSequence(random_aa(10))
                for width = 1:5, step = 1:5
                    aasitr = eachwindow(aas, width, step)
                    i, j = 1, width
                    for win in aasitr
                        @test win[1] == i
                        @test win[2] == sub(aas, i:j)
                        @test win[2] == aas[i:j]
                        i += step
                        j += step
                    end
                end
            end

            @testset "Arrays" begin
                arr = collect(1:10)
                for width = 1:5, step = 1:5
                    arritr = eachwindow(arr, width, step)
                    i, j = 1, width
                    for win in arritr
                        @test win[1] == i
                        @test win[2] == sub(arr, i:j)
                        i += step
                        j += step
                    end
                end
            end

            @testset "Strings" begin
                str = randstring(10)
                for width = 1:5, step = 1:5
                    stritr = eachwindow(str, width, step)
                    i, j = 1, width
                    for win in stritr
                        @test win[1] == i
                        @test win[2] == str[i:j]
                        i += step
                        j += step
                    end
                end
            end

        end
    end
end

end


text:
module VortexHelperDashBubble
using ..Ahorn, Maple

@mapdef Entity "VortexHelper/DashBubble" DashBubble(x::Integer, y::Integer,
    spiked::Bool = false, singleUse::Bool = false, wobble::Bool = true)

const placements = Ahorn.PlacementDict(
    "Dash Bubble (Vortex Helper)" => Ahorn.EntityPlacement(
        DashBubble,
		"point"
    ),
    "Dash Bubble (Spiked) (Vortex Helper)" => Ahorn.EntityPlacement(
        DashBubble,
		"point",
		Dict{String, Any}(
			"spiked" => true
        )
    )
)

function getSprite(entity::DashBubble)
    spiked = get(entity.data, "spiked", false)
    if spiked
        return "objects/VortexHelper/dashBubble/spiked00"
    end
    return "objects/VortexHelper/dashBubble/idle00"
end

function Ahorn.selection(entity::DashBubble)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::DashBubble, room::Maple.Room)
    Ahorn.drawSprite(ctx, getSprite(entity), 0, 0)
end

end

max_stars_repo_name: catapillie/VortexHelper

max_issues_repo_name: catapillie/VortexHelper

max_forks_repo_name: catapillie/VortexHelper


text:
using Distributions, PDMats, MCMC

# Using Optim, M would be ~ 106, acceptance rate ~ 12%
#using Optim

using Base.Test
#srand(1)

m = [0.0, 2.0]
sig = 0.4
mu = [1.0, 1.0]
sigma = 4

# Target distribution shape, unscaled
d = DiagNormal(m, PDiagMat(sig^2 * ones(2)))
g(x) = exp(sum(logpdf(d, x)))
log_g(x) = sum(logpdf(d, x))

# Candidate (starting) Distribution, use MvNormal([1, 1], [1 0; 0 1])
d0 = DiagNormal(mu, PDiagMat(sigma^2 * ones(2)))
g0(x) = exp(sum(logpdf(d0, x)))
log_g0(x) = sum(logpdf(d0, x))

# Compute M such that candidate dominates target
h(x) = -g(x)/g0(x)
#res = optimize(h, ones(2))
#M = -res.f_minimum
M = 120.0
log_M = log(M)

mcmodel = model(log_g, init=ones(2))

mcchain = run(mcmodel, ARS(log_g0, log_M), SerialMC(1000:1:10000))

println()
println("M: ", M)
println()

println("Acceptance rate g0: ", acceptance(mcchain))
println()
describe(mcchain)


text:
include("include.jl")
include("read_lp.jl")

####################################
### PLOT A TOY PROBLEM
####################################
using JuMP, Ipopt, PyPlot
m = Model()

A = [[1.0 1.0 0.0]; [1.0 0.0 -1.0]]
b, c, lbounds, ubounds = [0.0 0.0]', [0.0 0.0 0.0]', [0.0 0.0 0.0]', [1e30 1e30 1e30]'

ls = Dict(1e-3 => "-", 1e-6 => "-.", 1e-9 => "--", 0.0 => ":")
delta_set = [1e-3, 1e-6, 1e-9, 0.0]

ylims = [1e-2, 1e10]

subplot(121)
for delta = delta_set
  max_iter = 100
  hist = run_IPOPT(A, b, c, lbounds, ubounds, max_iter, delta)
  dual_hist = compute_dual_hist(hist)

  semilogy(1:length(dual_hist), dual_hist, color="black", linestyle=ls[delta], label="δ = $delta", basey=10)
end
title("IPOPT")
ax = gca()
ax[:set_ylim](ylims)
xlabel("iteration")
ylabel("maximum dual variable")

legend()

function run_one_phase(jump_model)
    nlp_raw = MathProgNLPModel(jump_model)
    nlp = Class_CUTEst(nlp_raw)
    timer = class_advanced_timer()
    start_advanced_timer(timer)
    start_advanced_timer(timer, "INIT")
    intial_it = init(nlp, my_par, timer)
    pause_advanced_timer(timer, "INIT")
    @assert(is_feasible(intial_it, my_par.ls.comp_feas))
    iter, status, hist, t, err = one_phase_IPM(intial_it, my_par, timer);
    pause_advanced_timer(timer)

    return hist
end

subplot(122)
for delta = delta_set
  jump_model = read_lp_into_JuMP(A, b, c, lbounds, ubounds, delta, true)
  hist = run_one_phase(jump_model)
  dual_hist = compute_dual_hist(hist)

  semilogy(1:length(dual_hist), dual_hist, color="black", linestyle=ls[delta], label="δ = $delta", basey=10)
end
title("Well behaved IPM")
xlabel("iteration")
ax = gca()
ax[:set_ylim](ylims)

#ax[:legend]()[:draggable]()


text:
@testset "LinRange" begin
    for R in (LinMRange,LinSRange)
        @testset "$R" begin
            r = R(1, 4, 4)
            b = LinRange(1, 4, 4)
            @test reverse(r) == reverse(b)
            @test R(r) == r
            @test R(1:4) == r
            @test -(r) == -(b)
            @test -(r, R(2, 5, 4)) == -(b, LinRange(2, 5, 4))
            @test +(r, R(2, 5, 4)) == +(b, LinRange(2, 5, 4))

            @test R(1,1,1) == LinRange(1, 1, 1)

            @test intersect(r, r[2]) == intersect(b, b[2])
            @test intersect(r, r[2]) == intersect(b, b[2])

            @test R{Float64}(r) == LinRange{Float64}(r)

            @test_throws ErrorException r.notfield
            # issue #20380
            let r = R(1,4,4)
                @test isa(r[UnitSRange(1, 4)], StaticRanges.AbstractLinRange)
            end

            if R == LinMRange
               setproperty!(r, :start, 2)
                @test r == LinMRange(2, 4, 4)

                setproperty!(r, :stop, 10)
                @test r == LinMRange(2, 10, 4)

                setproperty!(r, :len, 9)
                @test r == LinMRange(2, 10, 9)

                setproperty!(r, :lendiv, 8)
                @test r == LinMRange(2, 10, 9)

                @test_throws ErrorException setproperty!(r, :anything, 3)
            end
        end
    end
end


text:
function make_circle()
    img = zeros(Gray{N0f8},50, 50)
    draw!(img, CirclePointRadius(Point(25,25), 20))
    return img
end


text:
function get_return_type(test::TransferEntropyCausalityTest)
    T = typeof(1.0)

    if test.ηs isa Number
        return T
    else test.ηs isa Vector
        return Vector{T}
    end
end

function get_return_type(test::CrossMappingTest)
    T = typeof(1.0)

    return Vector{T}
end

function get_return_type(test::ConvergentCrossMappingTest)
    T = typeof(1.0)

    return Vector{Vector{T}}
end

function get_return_type(test::JointDistanceDistributionTest)
    T = typeof(1.0)

    return Vector{T}
end

function get_return_type(test::JointDistanceDistributionTTest)
    return OneSampleTTest
end

function get_return_type(test::SMeasureTest)
    T = typeof(1.0)
    return Vector{T}
end

function get_return_type(test::AbstractPredictiveAsymmetryTest)
    T = typeof(1.0)
    return Vector{T}
end

export get_return_type


text:
using Plots

gr()

xs=-pi:0.01:pi

anim = @animate for i = 1:100
    plot(xs,sin.(0.05i*xs))
end

Plots.gif(anim,"anim_fps15.gif",fps=15)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using DashBootstrapComponents, DashHtmlComponents

toast = html_div([
    dbc_button(
        "Open toast",
        id="positioned-toast-toggle",
        color="primary",
        n_clicks=0,
    ),
    dbc_toast(
        "This toast is placed in the top right",
        id="positioned-toast",
        header="Positioned toast",
        is_open=false,
        dismissable=true,
        icon="danger",
        # top: 66 positions the toast below the navbar
        style=Dict("position" => "fixed", "top" => 66, "right" => 10, "width" => 350),
    ),
]);

callback!(
    app,
    Output("positioned-toast", "is_open"),
    Input("positioned-toast-toggle", "n_clicks"),
) do n
    return n > 0
end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
module PlutoBinder

f() = 3

end

max_stars_repo_name: rikhuijzer/PlutoBinder

max_issues_repo_name: rikhuijzer/PlutoBinder

max_forks_repo_name: rikhuijzer/PlutoBinder


text:
using Base: Integer, AbstractVecOrTuple
using MPI, HDF5, DataFrames

####################
# Helper functions #
####################

isa_df(obj) = (@isdefined(DataFrames) && @isdefined(AbstractDataFrame)) && obj isa DataFrames.AbstractDataFrame
isa_gdf(obj) = (@isdefined(DataFrames) && @isdefined(GroupedDataFrame)) && obj isa DataFrames.GroupedDataFrame
isa_array(obj) = obj isa AbstractArray || obj isa HDF5.Dataset

get_worker_idx(comm::MPI.Comm) = MPI.Comm_rank(comm) + 1
get_nworkers(comm::MPI.Comm) = MPI.Comm_size(comm)

get_partition_idx(batch_idx, nbatches, comm::MPI.Comm) =
    get_partition_idx(batch_idx, nbatches, get_worker_idx(comm))

get_partition_idx(batch_idx, nbatches, worker_idx) =
    (worker_idx - 1) * nbatches + batch_idx

get_npartitions(nbatches, comm::MPI.Comm) = nbatches * get_nworkers(comm)

split_len(src_len::Integer, idx::Integer, npartitions::Integer) =
    if npartitions > 1
        # dst_len = Int64(cld(src_len, npartitions))
        dst_len = cld(src_len, npartitions)
        dst_start = min((idx - 1) * dst_len + 1, src_len + 1)
        dst_end = min(idx * dst_len, src_len)
        dst_start:dst_end
    else
        1:src_len
    end

split_len(src_len, batch_idx::Integer, nbatches::Integer, comm::MPI.Comm) = split_len(
    src_len,
    get_partition_idx(batch_idx, nbatches, comm),
    get_npartitions(nbatches, comm),
)

split_on_executor(src::AbstractArray, d::Integer, i) = selectdim(src, d, i)
split_on_executor(src::AbstractDataFrame, d::Integer, i) = @view src[i, :]
split_on_executor(src::GroupedDataFrame, d::Integer, i) = nothing

# In case we are trying to `Distribute` a grouped data frame,
# we can't do that so we will simply return nothing so that the groupby
# partitioned computation will redo the groupby.

split_on_executor(
    src::Nothing,
    dim::Integer,
    batch_idx::Integer,
    nbatches::Integer,
    comm::MPI.Comm,
) = nothing

split_on_executor(
    src::GroupedDataFrame,
    dim::Integer,
    batch_idx::Integer,
    nbatches::Integer,
    comm::MPI.Comm,
) = nothing

split_on_executor(
    src::T,
    dim::Integer,
    batch_idx::Integer,
    nbatches::Integer,
    comm::MPI.Comm,
) where {T} = begin
    npartitions = get_npartitions(nbatches, comm)
    if npartitions > 1
        split_on_executor(
            src,
            dim,
            split_len(
                size(src, dim),
                get_partition_idx(batch_idx, nbatches, comm),
                npartitions,
            ),
        )
    else
        src
    end
end

function merge_on_executor(obj::Vararg{AbstractArray{T,N},M}; key = nothing) where {T,N,M}
    if length(obj) == 1
        obj[1]
    else
        cat(obj...; dims = key)
    end
end

# If this is a dataframe then we ignore the grouping key
function merge_on_executor(obj::Vararg{DataFrame,M}; key = nothing) where {M}
    if length(obj) == 1
        obj[1]
    else
        vcat(obj...)
    end
end

merge_on_executor(obj::Vararg{GroupedDataFrame,M}; key = nothing) where {M} = nothing
merge_on_executor(obj::Vararg{T,M}; key = nothing) where {T,M} = first(obj)

# TODO: Make `merge_on_executor` and `tobuf` and `frombuf`
# dispatch based on the `kind` so we only have to precompile Arrow if we are
# working with dataframes.
function merge_on_executor(kind::Symbol, vbuf::MPI.VBuffer, nchunks::Integer; key)
    chunks = [
        begin
            chunk = view(vbuf.data, (vbuf.displs[i]+1):(vbuf.displs[i]+vbuf.counts[i]))
            if kind == :df
                DataFrames.DataFrame(Arrow.Table(IOBuffer(chunk)))
            elseif kind == :bits
                chunk
            else
                deserialize(IOBuffer(chunk))
            end
        end for i = 1:nchunks
    ]
    merge_on_executor(chunks...; key = key)
end

function get_partition_idx_from_divisions(
    val,
    divisions;
    boundedlower = false,
    boundedupper = false,
)
    # The first and last partitions (used if this lacks a lower or upper bound)
    # must have actual division(s) associated with them. If there is no
    # partition that has divisions, then they will all be skipped and -1 will
    # be returned. So these indices are only used if there are nonempty
    # divisions.
    firstdivisionidx = findfirst(x->!isempty(x), divisions)
    lastdivisionidx = findlast(x->!isempty(x), divisions)

    # The given divisions may be returned from `get_divisions`
    oh = orderinghash(val)
    for (i, div) in enumerate(divisions)
        # Now _this_ is a plausible cause. `get_divisions` can return a bunch
        # of empty arrays and in that case we should just skip that division.
        if isempty(div)
            continue
        end

        isfirstdivision = i == firstdivisionidx
        islastdivision = i == lastdivisionidx
        if ((!boundedlower && isfirstdivision) || oh >= first(div)[1]) &&
           ((!boundedupper && islastdivision) || oh < last(div)[2])
            return i
        end
    end

    # We return -1 since this value doesn't belong to any of the partitions
    # represented by `divisions`.
    -1
end

isoverlapping(a::AbstractRange, b::AbstractRange) = a.start ≤ b.stop && b.start ≤ a.stop

######################################################
# Helper functions for serialization/deserialization #
######################################################

to_jl_value(jl) = Dict("is_banyan_value" => true, "contents" => to_jl_value_contents(jl))

# NOTE: This function is shared between the client library and the PT library
to_jl_value_contents(jl) = begin
    # Handle functions defined in a module
    # TODO: Document this special case
    # if jl isa Function && !(isdefined(Base, jl) || isdefined(Core, jl) || isdefined(Main, jl))
    if jl isa Expr && eval(jl) isa Function
        jl = Dict("is_banyan_udf" => true, "code" => jl)
    end

    # Convert Julia object to string
    io = IOBuffer()
    iob64_encode = Base64EncodePipe(io)
    serialize(iob64_encode, jl)
    close(iob64_encode)
    String(take!(io))
end

# NOTE: This function is shared between the client library and the PT library
from_jl_value_contents(jl_value_contents) = begin
    # Converty string to Julia object
    io = IOBuffer()
    iob64_decode = Base64DecodePipe(io)
    write(io, jl_value_contents)
    seekstart(io)
    res = deserialize(iob64_decode)

    # Handle functions defined in a module
    if res isa Dict && haskey(res, "is_banyan_udf") && res["is_banyan_udf"]
        eval(res["code"])
    else
        res
    end
end

##########################################
# Ordering hash for computing  divisions #
##########################################

# NOTE: `orderinghash` must either return a number or a vector of
# equally-sized numbers
# NOTE: This is an "order-preserving hash function" (google that for more info)
orderinghash(x::Any) = x # This lets us handle numbers and dates
orderinghash(s::AbstractString) = Integer.(codeunits(first(s, 32) * repeat(" ", 32-length(s))))
orderinghash(A::AbstractArray) = orderinghash(first(A))

to_vector(v::Vector) = v
to_vector(v) = [v]

function get_divisions(divisions, npartitions)
    # This function accepts a list of divisions where each division is a tuple
    # of ordering hashes (values returned by `orderinghash` which are either
    # numbers or vectors of numbers). It also accepts a number of partitions to
    # produce divisions for. The result is a list of length `npartitions`
    # containing lists of divisions for each partition. A partition may contain
    # multiple divisions.

    ndivisions = length(divisions)
    if ndivisions == 0
        # If there are no divisions (maybe this dataset or this partition of a
        # dataset is empty), we simply return empty set.
        [[] for _ in 1:npartitions]
    elseif ndivisions >= npartitions
        # If there are more divisions than partitions, we can distribute them
        # easily. Each partition gets 0 or more divisions.
        # TODO: Ensure usage of div here and in sampling (in PT
        # library (here), annotation, and in locations) doesn't result in 0 or
        # instead we use ceiling division
        # ndivisions_per_partition = div(ndivisions, npartitions)
        [
            # This could be an empty array.
            begin
                divisions[split_len(ndivisions, partition_idx, npartitions)]
            end for partition_idx = 1:npartitions
        ]
    else
        # Otherwise, each division must be shared among 1 or more partitions
        allsplitdivisions = []
        # npartitions_per_division = div(npartitions, ndivisions)

        # Iterate through the divisions and split each of them and find the
        # one that contains a split that this partition must own and use as
        # its `partition_divisions`
        for (division_idx, division) in enumerate(divisions)
            # Determine the range (from `firstpartitioni` to `lastpartitioni`) of
            # partitions that own this division
            # islastdivision = division_idx == ndivisions
            # firstpartitioni = ((division_idx-1) * npartitions_per_division) + 1
            # lastpartitioni = islastdivision ? npartitions : division_idx * npartitions_per_division
            # partitionsrange = firstpartitioni:lastpartitioni
            partitionsrange = split_len(npartitions, division_idx, ndivisions)

            # # If the current partition is in that division, compute the
            # # subdivision it should use for its partition
            # if partition_idx in partitionsrange

            # We need to split the division among all the partitions in
            # its range
            ndivisionsplits = length(partitionsrange)

            # Get the `Vector{Number}`s to interpolate between
            divisionbegin = to_vector(first(division))
            divisionend = to_vector(last(division))

            # @show divisionbegin
            # @show divisionend

            # Initialize divisions for each split
            splitdivisions =
                [[copy(divisionbegin), copy(divisionend)] for _ = 1:ndivisionsplits]

            # Adjust the divisions for each split to interpolate. The result
            # of an `orderinghash` call can be an array (in the case of
            # strings), so we must iterate through that array in order to
            # interpolate at the first element in that array where there is a
            # difference.
            for (i, (dbegin, dend)) in enumerate(zip(divisionbegin, divisionend))
                # Find the first index in the `Vector{Number}` where
                # there is a difference that we can interpolate between
                if dbegin != dend
                    # dpersplit = div(dend-dbegin, ndivisionsplits)
                    # Iterate through each split
                    # @show dpersplit
                    # @show dbegin
                    # @show dend
                    start = copy(dbegin)
                    for j = 1:ndivisionsplits
                        # Update the start and end of the division
                        # islastsplit = j == ndivisionsplits
                        splitdivisions[j][1][i] = j == 1 ? dbegin : copy(start)
                        start += cld(dend - dbegin, ndivisionsplits)
                        start = min(start, dend)
                        splitdivisions[j][2][i] = j == ndivisionsplits ? dend : copy(start)
                        # splitdivisions[j][1][i] = dbegin + (dpersplit * (j-1))
                        # splitdivisions[j][2][i] = islastsplit ? dend : dbegin + dpersplit * j

                        # Ensure that the remaining indices are matching between the start and end.
                        if j < ndivisionsplits
                            splitdivisions[j][2][i+1:end] = splitdivisions[j][1][i+1:end]
                        end
                    end

                    # Stop if we have found a difference we can
                    # interpolate between
                    # TODO: If the difference is not that much,
                    # interpolate between multiple consecutive
                    # differeing characters together
                    break
                end
            end

            # Convert back to `Number` if the divisions were originally
            # `Number`s. We support either numbers or lists of numbers for the
            # ordering hashes that we use for the min-max bounds.
            if !(first(division) isa Vector)
                splitdivisions = [
                    # NOTE: When porting this stuff to Python, be sure
                    # to take into account the fact that Julia treats
                    # many values as arrays
                    (first(splitdivisionbegin), first(splitdivisionend)) for
                    (splitdivisionbegin, splitdivisionend) in splitdivisions
                ]
            end

            # # Get the split of the division that this partition should own
            # splitdivision = splitdivisions[1+partition_idx-first(partitionsrange)]

            # # Stop because we have found a division that this partition
            # # is supposed to own a split from
            # break

            # Each partition must have a _list_ of divisions so we must have a list
            # for each partition. So `allsplitdivisions` is an array where
            # each element is either a 1-element array containing a single
            # division or its empty.
            for splitdivision in splitdivisions
                # Check if we have already added this split division before.
                # The last split division may have been empty but we can
                # still check whether there is a last one and if what we're
                # adding is the same or also empty. If it is the same or also
                # empty, then we just add an empty divisions list. Otherwsie,
                # we add in our novel split division.
                if !isempty(allsplitdivisions) && last(allsplitdivisions) == splitdivision
                    push!(allsplitdivisions, [])
                else
                    push!(allsplitdivisions, [splitdivision])
                end
            end

            # end
        end

        allsplitdivisions
    end
end

########################
# Helper MPI functions #
########################

# TODO: Fix below function for reducing values of non-equal sizes
# TODO: Make Allreducev version of below function

# function Reducev(value, op, comm::MPI.Comm)
#     # Reduces values on all processes to a single value on rank 0.
#     #
#     # This function does the same thing as the function MPI_Reduce using
#     # only MPI_Send and MPI_Recv. As shown, it operates with additions on
#     # integers, so you could trivially use MPI_Reduce, but for operations
#     # on variable size structs for which you cannot define an MPI_Datatype,
#     # you can still use this method, by modifying it to use your op
#     # and your datastructure.

#     # TODO: Actually determine buffer
#     tag = 0
#     size = get_nworkers(comm)
#     rank = get_worker_idx(comm)-1
#     lastpower = 1 << log2(size)

#     # each of the ranks greater than the last power of 2 less than size
#     # need to downshift their data, since the binary tree reduction below
#     # only works when N is a power of two.
#     for i in lastpower:(size-1)
#         if rank == i
#             MPI.send(value, i-lastpower, tag, comm)
#         end
#     for i in 0:(size-lastpower-1)
#         if rank == i
#             MPI.Recv!(recvbuffer, i+lastpower, tag, comm)
#             value = op(value, recvbuffer)
#         end
#     end

#     for d in 0:(fastlog2(lastpower)-1)
#         k = 0
#         while k < lastpower
#             k += 1 << (d + 1)
#         end
#         receiver = k
#         sender = k + (1 << d)
#         if rank == receiver
#             MPI.Recv!(recvbuffer, 1, sender, tag)
#             value = op(value, recvbuffer)
#         elseif rank == sender
#             MPI.Send(value, 1, receiver, tag)
#         end
#     end
#     value
# end

# function fastlog2(v::UInt32)
#     multiply_de_bruijn_bit_position::Vector{Int32} = [
#         0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30,
#         8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31
#     ]

#     v |= v >> 1 # first round down to one less than a power of 2
#     v |= v >> 2
#     v |= v >> 4
#     v |= v >> 8
#     v |= v >> 16

#     # TODO: Fix this
#     multiply_de_bruijn_bit_position[(UInt32(v * 0x07C4ACDDU) >> 27) + 1]
# end

# function tobuf(obj)::Tuple{Symbol, MPI.Buffer}
function tobuf(obj)
    # We pass around Julia objects between MPI processes in different ways
    # depending on the data type. For simple isbitstype data we keep it as-is
    # and use the simple C-like data layout for fast transfer. For dataframes,
    # we use Arrow data layout for zero-copy deserialization. For everything
    # else including variably-sized arrays and arbitrary Julia objects, we
    # simply serialize and deserialize using the Serialization module in Julia
    # standard library.

    if isbits(obj)
        (:bits, MPI.Buffer(Ref(obj)))
        # (:bits, MPI.Buffer(obj))
        # (:bits, MPI.Buffer(Ref(obj)))
    elseif isa_array(obj) && isbitstype(first(typeof(obj).parameters)) && ndims(obj) == 1
        # (:bits, MPI.Buffer(obj))
        (:bits, MPI.Buffer(obj))
    elseif isa_df(obj)
        io = IOBuffer()
        Arrow.write(io, obj)
        # (:df, MPI.Buffer(view(io.data, 1:position(io))))
        (:df, MPI.Buffer(view(io.data, 1:io.size)))
    else
        io = IOBuffer()
        serialize(io, obj)
        (:unknown, MPI.Buffer(view(io.data, 1:io.size)))
        # (:unknown, io)
    end
end

function buftovbuf(buf::MPI.Buffer, comm::MPI.Comm)::MPI.VBuffer
    # This function expects that the given buf has buf.data being an array.
    # Basically what it does is it takes the result of a call to tobuf above
    # on each process and constructs a VBuffer with the sum of the sizes of the
    # buffers on different processes.
    sizes = MPI.Allgather(buf.count, comm)
    # NOTE: This function should only be used for variably-sized buffers for
    # receiving data because the returned buffer contains zeroed-out memory.
    VBuffer(similar(buf.data, sum(sizes)), sizes)
end

function bufstosendvbuf(bufs::Vector{MPI.Buffer}, comm::MPI.Comm)::MPI.VBuffer
    sizes = [length(buf.data) for buf in bufs]
    VBuffer(vcat(map(buf -> buf.data, bufs)), sizes)
end

function bufstorecvvbuf(bufs::Vector{MPI.Buffer}, comm::MPI.Comm)::MPI.VBuffer
    # This function expects that each given buf has buf.data being an array and
    # that the number of bufs in bufs is equal to the size of the communicator.
    # sizes = MPI.Allgather(length(buf.data), comm)
    sizes = MPI.Alltoall([length(buf.data) for buf in bufs])
    # NOTE: Ensure that the data fields of the bufs are initialized to have the
    # right data type (e.g., Vector{UInt8} or Vector{Int64})
    # We use `similar` here because we want zeroed out memory to receive data.
    VBuffer(similar(first(bufs).data, sum(sizes)), sizes)
end

function frombuf(kind, obj)
    if kind == :bits && obj isa Ref
        # TODO: Ensure that the "dereferece" here is necessary
        obj[]
    elseif kind == :bits
        obj
    elseif kind == :df
        DataFrames.DataFrame(Arrow.Table(obj), copycols = false)
    else
        deserialize(obj)
    end
end

function getpath(path)
    if startswith(path, "http://") || startswith(path, "https://")
        # TODO: First check for size of file and only download to
        # disk if it doesn't fit in free memory
        # TODO: Add option for Internet locations as to whether or not to
        # cache on disk
        hashed_path = string(hash(path))
        joined_path = "efs/banyan_dataset_" * hashed_path
        # @info "Downloading $path to $joined_path"
        if !isfile(joined_path)
            # NOTE: Even though we are storing in /tmp, this is
            # effectively caching the download. If this is undesirable
            # to a user, a short-term solution is to use a different
            # URL each time (e.g., add a dummy query to the end of the
            # URL)
            Downloads.download(path, joined_path)
        end
        # @show isfile(joined_path)
        joined_path
    elseif startswith(path, "s3://")
        replace(path, "s3://" => "/home/ec2-user/s3fs/")
        # NOTE: We expect that the ParallelCluster instance was set up
        # to have the S3 filesystem mounted at ~/s3fs/<bucket name>
    else
        # Case of local paths to things stored on disk
        "efs/"*path
    end
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
"""Solve a problem that has an exact solution and compare the solutions."""
module ExampleExact

using Fronts
using Plots

# Reference: Philip (1960) Table 1, No. 13
# https://doi.org/10.1071/PH600001
# Exact solution: θ(ϕ) = exp(-ϕ)
prob = DirichletProblem(θ -> 0.5*(1 - log(θ)), i=0, b=1)

θ = solve(prob)

ϕ = range(0, 20, length=200)

plt = plot(θ, label="Fronts", xguide="ϕ", yguide="θ", legend=true)
plot!(ϕ, exp.(-ϕ), label="Exact")
display(plt)

end


text:
using Test: AbstractTestSet, @testset, @test
import Test

mutable struct SingleParameterTestSet <: AbstractTestSet
    parameters::Dict{Symbol,Any}
    results::Vector
    n_passed::Int

    SingleParameterTestSet(::String) = new(Dict(), [], 0)
end

mutable struct ParameterTestSet <: AbstractTestSet
    results::Vector
    ParameterTestSet(::String) = new([])
end

function quietrecord(ts::Test.DefaultTestSet, t::Union{Test.Fail, Test.Error}, pts::SingleParameterTestSet)
    if Test.TESTSET_PRINT_ENABLE[]
        params = join(["$k=$v" for (k,v) in pts.parameters], ",")
        Test.printstyled("$(ts.description): [$params]\n", color=:white)
        # don't print for interrupted tests
        if !(t isa Test.Error) || t.test_type !== :test_interrupted
            print(t)
            println()
        end
    end
    push!(ts.results, t)
    return t
end

function Test.finish(ts::ParameterTestSet)
    # We'll capture them all for the moment, but we probably need to wrap this a little bit with
    # the relevant parameters.
    parent_ts = Test.get_testset()
    parent_ts.n_passed = foldl(+, map(x -> x.n_passed, filter(x -> isa(x, SingleParameterTestSet), ts.results)); init=0)

    actualresults = filter(x -> !isa(x, SingleParameterTestSet), ts.results)
    if length(actualresults) > 0
        for result in actualresults
            Test.record(parent_ts, result)
        end
    else
        failedsets = filter(x -> length(x.results) > 0, ts.results)
        # We just detail errors from the first case.
        if length(failedsets) > 0
            for result in failedsets[1].results
                quietrecord(parent_ts, result, failedsets[1])
            end
            for set in failedsets[2:end]
                for result in set.results
                    push!(parent_ts.results, result)
                end
            end
        end
    end

    return ts
end

Test.record(ts::ParameterTestSet, t::SingleParameterTestSet) = (push!(ts.results, t); ts)
Test.record(ts::ParameterTestSet, t::Union{Test.Broken, Test.Fail, Test.Error}) = (push!(ts.results, t); ts)

Test.record(ts::SingleParameterTestSet, t::Union{Test.Broken, Test.Fail, Test.Error}) = (push!(ts.results, t); t)
Test.record(ts::SingleParameterTestSet, t::Test.Pass) = (ts.n_passed += 1; t)
function Test.finish(ts::SingleParameterTestSet)
    parent_ts = Test.get_testset()
    Test.record(parent_ts, ts)
    return ts
end

struct ParamTestSettings
    samples::Int

    ParamTestSettings(;samples=100) = new(samples)
end

PARAM_TEST_SETTINGS = nothing

function withsettings(f::Base.Callable, settings::ParamTestSettings)
    global PARAM_TEST_SETTINGS
    current = PARAM_TEST_SETTINGS
    try
        PARAM_TEST_SETTINGS = settings
        f()
    finally
        PARAM_TEST_SETTINGS = current
    end
end