text:
# # [Implicit-Explicit (IMEX) Additively-Partitioned Runge-Kutta Timestepping](@id Single-rate-IMEXARK-Timestepping)

# In this tutorial, we shall explore the use of IMplicit-EXplicit (IMEX) methods
# for the solution of nonautonomous (or non time-invariant) equations.
# For our model problem, we shall reuse the acoustic wave test in the GCM
# configuration. See its [code](@ref Acoustic-Wave-Configuration)
# for details on the model and parameters. For the purposes of this tutorial,
# we will only run the experiment for a total of 3600 simulation seconds.
# Details on this test case can be found in Sec. 4.3 of [Giraldo2013](@cite).

using ClimateMachine
const clima_dir = dirname(dirname(pathof(ClimateMachine)));
include(joinpath(
    clima_dir,
    "tutorials",
    "Numerics",
    "TimeStepping",
    "tutorial_acousticwave_config.jl",
));

# The acoustic wave test case used in this tutorial represents a global-scale
# problem with inertia-gravity waves traveling around the entire planet.
# It has a hydrostatically balanced initial state that is given a pressure
# perturbation.
# This initial pressure perturbation causes an acoustic wave to travel to
# the antipode, coalesce, and return to the initial position. The exact solution
# of this test case is simple in that the (linear) acoustic theory allows one
# to verify the analytic speed of sound based on the thermodynamics variables.
# The initial condition is defined as a hydrostatically balanced atmosphere
# with background (reference) potential temperature.

# To fully demonstrate the advantages of using an IMEX scheme over fully explicit
# schemes, we start here by going over a simple, fully explicit scheme. The
# reader can refer to the [Single-rate Explicit Timestepping tutorial](@ref Single-rate-Explicit-Timestepping)
# for detailes on such schemes. Here we use the the 14-stage LSRK method
# [`LSRK144NiegemannDiehlBusch`]((@ref ClimateMachine.ODESolvers.LSRK144NiegemannDiehlBusch)), which contains the largest stability region of
# the low-storage methods available in `ClimateMachine.jl`.

FT = Float64
timeend = FT(100)

ode_solver = ClimateMachine.ExplicitSolverType(
    solver_method = LSRK144NiegemannDiehlBusch,
);

# In the following example, the timestep calculation is based on the CFL condition
# for horizontally-propogating acoustic waves. We use a Courant number ``C = 0.002``
# (denoted by `CFL` in the code bellow) in the horizontal, which corresponds
# to a timestep size of approximately ``1`` second.

CFL = FT(0.002)
cfl_direction = HorizontalDirection()
run_acousticwave(ode_solver, CFL, cfl_direction, timeend);

# However, as it is imaginable, for real-world climate processes a time step
# of 1 second would lead to extemely long time-to-solution simulations.
# How can we do better? To be able to take larger time step, we can treat the
# most restrictive wave speeds (vertical acoustic) implicitly rather than
# explicitly. This motivates the use of an IMplicit-EXplicit (IMEX) methods.

# In general, a single step of an ``s``-stage, ``N``-part additive RK method
# (`ARK_N`) is defined by its generalized Butcher tableau:

# ```math
# \begin{align}
#     \begin{array}{c|c|c|c}
#     \boldsymbol{c} &\boldsymbol{A}_{1} & \cdots & \boldsymbol{A}_{N}\\
#     \hline
#     & \boldsymbol{b}_1^T & \cdots & \boldsymbol{b}_N^T\\
#     \hline
#     & \widehat{\boldsymbol{b}}_1^T & \cdots & \widehat{\boldsymbol{b}}_N^T
#     \end{array} =
#     \begin{array}{c|c c c | c | c c c }
#     c_1 & a^{[ 1 ]}_{1,1} & \cdots & a^{[ 1 ]}_{1,s} & \cdots
#     & a^{[ \nu ]}_{1,1} & \cdots & a^{[ \nu ]}_{1,s}\\
#     \vdots & \vdots & \ddots & \vdots & \cdots
#     & \vdots & \ddots & \vdots \\
#     c_s & a^{[ 1 ]}_{s,1} & \cdots & a^{[ 1 ]}_{s,s} & \cdots
#     & a^{[ \nu ]}_{s,1} & \cdots & a^{[ \nu ]}_{s,s}\\
#     \hline
#     & b^{[ 1 ]}_1 & \cdots & b^{[ 1 ]}_s & \cdots
#     & b^{[ \nu ]}_1 & \cdots & b^{[ \nu ]}_s\\
#     \hline
#     & \widehat{b}^{[ 1 ]}_1 & \cdots & \widehat{b}^{[ 1 ]}_s &
#     & \widehat{b}^{[ \nu ]}_1 & \cdots & \widehat{b}^{[ \nu ]}_s
#     \end{array}
# \end{align}
# ```

# and is given by

# ``
# 	\boldsymbol{q}^{n+1} = \boldsymbol{q}^n + \Delta t \left( \underbrace{\sum_{i=1}^{s}}_{\textrm{Stages}} \underbrace{\sum_{\nu=1}^{N}}_{\textrm{Components}} b_i^{[ \nu ]} {\mathcal{T}}^{[ \nu ]}(\boldsymbol{Q}^i)) \right)
# ``

# where ``s`` denotes the stages and ``N`` the components, and where the stage values are given by:
#
# ``
# 	\boldsymbol{Q}^i = \boldsymbol{q}^n + \Delta t \sum_{j=1}^{s} \sum_{\nu = 1}^{N} a_{i,j}^{[ \nu ]}
# 	{\mathcal{T}}^{[ \nu]}(\boldsymbol{Q}^j).
# ``
#
# Similar to standard RK methods, the stage vectors are approximations to the state at each stage
# of the ARK method. Moreover, the temporal coefficients ``c_i`` satisfy a similar
# row-sum condition, holding for all ``\nu = 1, \cdots, N``:

# ``
# 	c_i = \sum_{j=1}^{s} a_{i, j}^{[ \nu ]}, \quad \forall \nu = 1, \cdots, N.
# ``
#
# The Butcher coefficients ``\boldsymbol{c}``, ``\boldsymbol{b}_{\nu}``, ``\boldsymbol{A}_{\nu}``, and ``\widehat{\boldsymbol{b}}_{\nu}``
# are constrained by certain accuracy and stability requirements, which are summarized in
# [Kennedy2001](@cite).

# A common setting is the case ``N = 2``. This gives the typical context for
# Implicit-Explicit (IMEX) splitting methods, where the tendency ``{\mathcal{T}}``
# is assumed to have the decomposition:
#
# ``
# 	\dot{\boldsymbol{q}} = \mathcal{T}(\boldsymbol{q}) \equiv
# 	{\mathcal{T}}_{s}(\boldsymbol{q}) + {\mathcal{T}}_{ns}(\boldsymbol{q}),
# ``
# where the right-hand side has been split into a "stiff" component ``{\mathcal{T}}_{s}``,
# to be treated implicitly, and a non-stiff part ``{\mathcal{T}}_{ns}`` to be treated explicitly.

# Referencing the canonical form introduced in [Time integration](@ref
# Time-integration) we have that in this particular forumlation
# ``\mathcal{T}_{ns}(t, \boldsymbol{q}) \equiv \mathcal{G}(t, \boldsymbol{q})`` and
# ``\mathcal{T}_{s}(t, \boldsymbol{q}) \equiv \mathcal{F}(t, \boldsymbol{q})``.

# Two different RK methods are applied to ``{\mathcal{T}}_{s}`` and ``{\mathcal{T}}_{ns}``
# separately, which have been specifically designed and coupled. Examples can be found in
# [Giraldo2013](@cite). The Butcher Tableau for an `ARK_2` method will have the
# form
#
# ```math
# \begin{align}
#     \begin{array}{c|c|c}
#     \boldsymbol{c} &\boldsymbol{A}_E &\boldsymbol{A}_I\\
#     \hline
#     & \boldsymbol{b}_E^T & \boldsymbol{b}_I^T \\
#     \hline
#     & \widehat{\boldsymbol{b}}_E^T & \widehat{\boldsymbol{b}}_I^T
#     \end{array},
# \end{align}
# ```
#
# with
#
# ``
# 	\boldsymbol{A}_O = \left\lbrace a_{i, j}^O \right\rbrace, \quad
# 	\boldsymbol{b}_O = \left\lbrace b_{i}^O \right\rbrace, \quad
# 	\widehat{\boldsymbol{b}}_O = \left\lbrace \widehat{b}_{i}^O \right\rbrace,
# ``
#
# where ``O`` denotes the label (either ``E`` for explicit or ``I`` for implicit).

# For the acoustic wave example used here, we use 4th order polynomials in
# our discontinuous Galerkin approximation, with 6 elements in each horizontal
# direction and 4 elements in the vertical direction, on the cubed-sphere.
# This gives an effective minimal node-distance (distance between LGL nodes)
# of roughly 203000 m.
# As in the [previous tutorial](@ref Single-rate-Explicit-Timestepping),
# we can determine our ``\Delta t`` by specifying our desired horizontal
# Courant number ``C`` (the timestep calculation is based on the CFL condition
# for horizontally-propogating acoustic waves). In this very simple test case,
# we can use a value of 0.5, which corresponds to a time-step size of
# around 257 seconds. But for this particular example, even higher values
# might work.

timeend = FT(3600)
ode_solver = ClimateMachine.IMEXSolverType(
    solver_method = ARK2GiraldoKellyConstantinescu,
)
CFL = FT(0.5)
cfl_direction = HorizontalDirection()
run_acousticwave(ode_solver, CFL, cfl_direction, timeend);

# ## References
# - [Giraldo2013](@cite)
# - [Kennedy2001](@cite)


text:
## functions to fit the 2/3D spatial structure of OCs
##

## extract parameters
##

function spatialParameter(ocfile, ntest=10 , nbin=20, niter=10000 ,verbose=true, nburnin= 5000)
    let
        solC= []; sols= []; solm= []
        ηBest= 1e9 ; θBest= 0
        fitfound= false

        oc= CSV.read(ocfile, delim= ";")
        if verbose println("## file $ocfile read...") end

        ## binned density
        r2d,ρ2d,err2d= density2D(oc.Y, oc.Z, nbin)    ## binned radius
        # r3d,ρ3d,err3d= density3D(oc.X, oc.Y, oc.Z, nbin)

        ## local density
        ## not used yet...
        nloc= 3
        rad_2d, loc_2d= locdens2d(oc.Y, oc.Z, nloc)    ## density estimated for each point
        # rad_3d, loc_3d= locdens3d(oc.X,oc.Y, oc.Z, nloc)

        oc2d= sc2dcentered(length(oc.Y), oc.Y, oc.Z, rad_2d, loc_2d , nbin , r2d, ρ2d, err2d )

        println(oc2d.densbin)
        # println(oc2d.densbinerr)
        if verbose println("## Density profiles computed...") end

        ## prior estimation
        ##
        prior, priorDisp= priorGuess(r2d, ρ2d)

        if verbose
            println("## Prior Guess:")
            println(prior)
            println("## prior estimated...")
        end

        C_mean= prior.C ; C_disp= priorDisp.C
        s_mean= prior.s ; s_disp= priorDisp.s
        m_mean= prior.m ; m_disp= priorDisp.m

        pinit= mcmcCauchy(C_mean, C_disp, s_mean, s_disp, m_mean, m_disp, nburnin, niter)

        ## if true it takes the mean as first step, otherwise random??

        if verbose println("## Running mcmc 2D fitting") end

        for i in 1:ntest
            mci ,  θbest , ηbest = main_mcmc(oc2d, pinit, likelihood2dbin, false, verbose)

            lags= collect(1:niter)
            acfC= StatsBase.autocor(mci.C, lags)[end]
            acfs= StatsBase.autocor(mci.s, lags)[end]
            acfm= StatsBase.autocor(mci.m, lags)[end]

            if verbose
                println("## test: $i")
                println(θbest)
                println("## Likelihood: $ηbest")
                println("## ACF(C,s,m): $acfC , $acfs , $acfm")
            end

            ## AFC < 0.1 to be considered as not correlated...
            if abs(acfC) < 0.1 && abs(acfs) < 0.1 && abs(acfm) < 0.1 && θbest != 0
              push!(solC,abs(θbest.C))
              push!(sols,abs(θbest.s))
              push!(solm,abs(θbest.m))

                if ηbest < ηBest
                    fitfound= true
                    ηBest= ηbest
                    θBest= θbest
                end
            end
        end

        if !fitfound
          θBest= modelCauchy(1,1,1)
        end
        if length(solC) > 1
          θerr= modelCauchy(std(solC),std(sols),std(solm))
        else
          θerr= modelCauchy(1,1,1)
        end

        if θBest.C<0 || θBest.s<0 || θBest.m<0
            θBest.C= abs(θBest.C)
            θBest.s= abs(θBest.s)
            θBest.m= abs(θBest.m)
        end
        println("## OC fitting done.")
        return(θBest, θerr, fitfound)
    end
end

## Estimate  prior parameters
function priorGuess(rad, dens)
    m= 3
    mdisp= 2

    C= maximum(dens)
    Cdisp= C/2

    s= StatsBase.mean(rad, StatsBase.weights(convert(Array{Float64,1}, dens)))
    sdisp= s/2

    priorMean= modelCauchy(C, s, m)
    priorDisp= modelCauchy(Cdisp, sdisp, mdisp)

    return(priorMean, priorDisp)
end

####################
## radial density functions
####################

function volumeSphere(r)
    return(4 / 3 * π * r^3)
end

function areaDisk(r)
    return(π * r^2)
end

function density2D(x , y , nbin=10)
    center= [mean(x) ; mean(y)]

    nxy= length(x)
    A= Array{Float64}(undef,2,nxy)
    for i in 1:nxy
        A[1,i]= x[i]
        A[2,i]= y[i]
    end

    d= Distances.Euclidean()
    r= Distances.colwise(d, A, center)

    ρ= [] ; radius= [] ; err= []
    rmax= maximum(r)
    dr= rmax/nbin
    r0=0

    for i in 1:nbin
        r1= i*dr
        indr= (r .> r0) .& (r .<= r1)
        nstar= count(indr)
        dens= nstar/(areaDisk(r1)-areaDisk(r0))
        if nstar > 0
            errDensity= dens / sqrt(nstar)
        else
            errDensity= 1e9
        end

        push!(ρ, dens)
        push!(radius, (r0+r1)/2)
        push!(err, errDensity)
        r0= r1
    end

    return(radius, ρ , err)
end

function density3D(x , y , z, nbin=10)
    center= [mean(x) ; mean(y) ; mean(z)]
    println(center)

    nxy= length(x)
    A= Array{Float64}(undef,3,nxy)
    for i in 1:nxy
        A[1,i]= x[i]
        A[2,i]= y[i]
        A[3,i]= z[i]
    end

    d= Distances.Euclidean()
    r= Distances.colwise(d, A, center)

    ρ= [] ; radius= [] ; err=[]
    rmax= maximum(r)
    dr= rmax/nbin
    r0=0
    for i in 1:nbin
        r1= i*dr
        indr= (r .> r0) .& (r .<= r1)
        nstar= count(indr)
        dens= nstar/(volumeSphere(r1)-volumeSphere(r0))
        errDensity= dens / sqrt(nstar)
        push!(ρ, dens)
        push!(radius, (r1+r0)/2)
        push!(err, errDensity)
        r0= r1
    end

    return(radius, ρ , err)
end

## Estimation of the local density around each point
## nei: number of neighbor to estimate the density
##

function locdens2d(x , y , nei=10)
    nxy= length(x)
    A= Array{Float64}(undef,2,nxy)
    center= [mean(x) ; mean(y)]

    for i in 1:nxy
        A[:,i]= [x[i] y[i]]
    end

    d= Euclidean()
    p= pairwise(d, A, A, dims=2)
    r= Distances.colwise(d, A, center)
    locdensrad= []

    for i in 1:nxy
        rsort= sort(p[i,:])
        rmax= rsort[nei+1]
        locdens= (nei+1)/areaDisk(rmax)
        push!(locdensrad, locdens)
    end

    return(r,locdensrad)
end

function locdens3d(x , y , z, nei=10)
    nxyz= length(x)
    A= Array{Float64}(undef,3,nxyz)
    center= [mean(x) ; mean(y) ; mean(z)]

    for i in 1:nxyz
        A[:,i]= [x[i] y[i] z[i]]
    end

    d= Euclidean()
    p= pairwise(d, A, A, dims=2)
    r= Distances.colwise(d, A, center)
    locdensrad= []

    for i in 1:nxyz
        rsort= sort(p[i,:])
        rmax= rsort[nei+1]
        locdens= (nei+1)/volumeSphere(rmax)
        push!(locdensrad, locdens)
    end

    return(r,locdensrad)
end

## Main MCMC to estimate parameters
## firstvalue: true(set to mean), false(random from prior)
function main_mcmc(oc::sc2dcentered, mcmc::mcmcCauchy, likelihood::Function, firstvalue=false , verbose=true)
    let
        Random.seed!()
        mci = mcCauchy(zeros(Float64,0),zeros(Float64,0),zeros(Float64,0))
        mi, probi = theta(mcmc, oc, likelihood, firstvalue)
        if verbose println("## Init done...") end

        niter = mcmc.niter
        nburn = mcmc.nburnin
        nchain = 1
        loopAgain = true
        burndone = false
        ηbest= 0 ; θbest= 0

        while loopAgain
            micurrent, probcurrent = thetaiter(mi , mcmc, oc, likelihood)

        ### Metropolis-Hasting
            α = probcurrent / probi
            if α > rand()
                mi = micurrent
                probi = probcurrent
                nchain += 1
                if (nchain%50000 == 0 && verbose ) println("### chain:",nchain) end
                if nchain > nburn && !burndone  nchain = 0 ; burndone = true end    ## burnin done
                if nchain > niter loopAgain = false end
                if burndone
                    push!(mci.C, mi.C)
                    push!(mci.s, mi.s)
                    push!(mci.m, mi.m)
                end
            else
                nchain += 1
                if (nchain%50000 == 0 && verbose) println("### chain:",nchain) end
                if nchain > nburn && !burndone   nchain = 0 ; burndone = true end    ## burnin done
                if nchain > niter loopAgain = false end
                if burndone
                    push!(mci.C, mi.C)
                    push!(mci.s, mi.s)
                    push!(mci.m, mi.m)
                end
            end

            η= likelihood(mi, oc)
            if η > ηbest
                ηbest= η
                θbest= mi
            end
        end

        if verbose println("## MCMC 2D fitting done.") end

        return(mci, θbest, ηbest)
    end
end

## density and likelihood functions
function fdens1(r, θ::modelCauchy)
    r2= r*r
    s2= θ.s*θ.s
    val = θ.C / (1+r2/s2)^θ.m
    return(val)
end

function fdens1log(r, θ::modelCauchy)
    r2= r*r
    s2= exp(θ.s)*exp(θ.s)
    val = θ.C - θ.m * log(1+r2/s2)
    return(val)
end

## using local density for all points
function likelihood2d(θ::modelCauchy, oc::sc2dcentered)
    χ2= 0
    for i in 1:oc.nxy
        desti= fdens1(oc.radius[i], θ)
        χ2+= ((desti-oc.dens[i])/  oc.densbinerr[i])^2/oc.nxy     ## σ=1
    end
    p = min(1, exp(-χ2))
    #p= min(1, 1/χ2)
    return(p)
end

### binnned 2D density
function likelihood2dbin(θ::modelCauchy, oc::sc2dcentered)
    χ2= 0
    for i in 1:oc.nbin
        desti= fdens1(oc.radbin[i], θ)
        χ2+= ((desti-oc.densbin[i]) / oc.densbinerr[i])^2
    end
    p = min(1, exp(-χ2))
    return(p)
end

### binnned 2D/3D density log
### LOG
function likelihood2dbinlog(θ::modelCauchy, oc::sc2dcentered)
    χ2= 0
    for i in 1:oc.nbin
        desti= fdens1log(oc.radbin[i], θ)
        χ2+= ((desti-oc.densbin[i]) / 1)^2
        #println(i)
        #println(θ)
        #println(oc.densbin[i])
        #println(desti)
    end
    p = min(1, exp(-χ2))

    return(p)
end

## Model iterations...

function theta(mcmc::mcmcCauchy, oc::sc2dcentered, likelihood::Function, firstvalue::Bool)
    pC  = Normal(mcmc.Cmean, mcmc.Cdisp)
    pm  = Normal(mcmc.mmean, mcmc.mdisp)
    ps  = Normal(mcmc.smean, mcmc.sdisp)

    if firstvalue
        C= mcmc.Cmean
        s= mcmc.smean
        m= mcmc.mmean
    else
        C = rand(pC)
        s = rand(ps)
        m = rand(pm)
    end

    pdfC = pdf(pC,C)
    pdfm = pdf(pm, m)
    pdfs = pdf(ps, s)

    θ= modelCauchy(C, s, m)

    pdfC = pdf(pC , C)
    pdfm = pdf(pm, m)
    pdfs = pdf(ps, s)
    ptotal = pdfC*pdfm*pdfs*likelihood(θ, oc)

    return(θ, ptotal)
end

## next iteration
function thetaiter(θi::modelCauchy, mcmc::mcmcCauchy, oc::sc2dcentered, likelihood::Function)
    pC  = Normal(mcmc.Cmean, mcmc.Cdisp)
    pm  = Normal(mcmc.mmean, mcmc.mdisp)
    ps  = Normal(mcmc.smean, mcmc.sdisp)

    C_rw= Normal(θi.C, 0.5)
    s_rw= Normal(θi.s  , 0.5)
    m_rw= Normal(θi.m , 0.5)

    new_C = 0.
    new_m = 0
    new_s = 0

    new_C  = rand(C_rw)
    new_m =  rand(m_rw)
    new_s =  rand(s_rw)

    θ= modelCauchy(new_C, new_s, new_m)

    pdfC = pdf(pC , new_C)
    pdfm = pdf(pm, new_m)
    pdfs = pdf(ps, new_s)
    ptotal = pdfC*pdfm*pdfs*likelihood(θ, oc)

    return(θ, ptotal)
end

function model_rad(radius, θ, densModel::Function)
    ρfit= []
    for i in 1:length(radius)
        ρ= densModel(radius[i], θ)
        push!(ρfit, ρ)
    end

    return(ρfit)
end

ext: jl

max_stars_repo_name: bosscha/gaia-shock

max_issues_repo_name: bosscha/gaia-shock

max_forks_repo_name: bosscha/gaia-shock


text:
## Payoffs
abstract type AbstractPayoff{numtype <: Number} end

struct ControlVariates{abstractPayoff <: AbstractPayoff{<:Number}} <: AbstractMonteCarloMethod
    variate::AbstractPayoff
    conf_variate::MonteCarloConfiguration{<:Integer, <:Integer, <:AbstractMonteCarloMethod, <:BaseMode}
    function ControlVariates(variate::T1, conf_variate::MonteCarloConfiguration) where {T1 <: AbstractPayoff{<:Number}}
        return new{T1}(variate, conf_variate)
    end
end

abstract type SingleNamePayoff{numtype <: Number} <: AbstractPayoff{numtype} end
abstract type EuropeanPayoff{numtype <: Number} <: SingleNamePayoff{numtype} end
abstract type AmericanPayoff{numtype <: Number} <: SingleNamePayoff{numtype} end
abstract type BermudanPayoff{numtype <: Number} <: SingleNamePayoff{numtype} end
abstract type PathDependentPayoff{numtype <: Number} <: SingleNamePayoff{numtype} end
abstract type BarrierPayoff{numtype <: Number} <: PathDependentPayoff{numtype} end
abstract type AsianPayoff{numtype <: Number} <: PathDependentPayoff{numtype} end

#No multiple inheritance, sigh
abstract type BasketPayoff{numtype <: Number} <: AbstractPayoff{numtype} end
abstract type EuropeanBasketPayoff{numtype <: Number} <: BasketPayoff{numtype} end
abstract type AmericanBasketPayoff{numtype <: Number} <: BasketPayoff{numtype} end
abstract type BarrierBasketPayoff{numtype <: Number} <: BasketPayoff{numtype} end
abstract type AsianBasketPayoff{numtype <: Number} <: BasketPayoff{numtype} end

function maturity(x::abstractPayoff) where {abstractPayoff <: AbstractPayoff{<:Number}}
    return x.T
end

"""
General Interface for Payoff computation from MonteCarlo paths

		Payoff=payoff(S,optionData,rfCurve,T1=optionData.T)

Where:\n
		S           = Paths of the Underlying.
		optionData  = Datas of the Option.
		rfCurve  = Datas of the Risk Free Curve.
		T1 [Optional, default to optionData.T]=Final Time of Spot Simulation (default equals Time to Maturity of the option)

		Payoff      = payoff of the Option.

"""
function payoff(S::AbstractMatrix, optionData::AbstractPayoff, rfCurve::AbstractZeroRateCurve, T1::Number = optionData.T)
    error("Function used just for documentation")
end

function payoff(S::Array{Array{num, 2}, 1}, optionData::SingleNamePayoff, rfCurve::abstractZeroRateCurve, T1::num2 = optionData.T) where {abstractZeroRateCurve <: AbstractZeroRateCurve, num <: Number, num2 <: Number}
    #switch to only in julia 1.4
    @assert length(S) == 1
    return payoff(S[1], optionData, rfCurve, T1)
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    abstract type ENSProposal

*Experimental feature, not part of stable public API.*

Abstract type for the algorithms to propose new live points used used by EllipsoidalNestedSampling.
"""
abstract type ENSProposal end

"""
    struct BAT.ENSUniformly <: ENSProposal

*Experimental feature, not part of stable public API.*

Each point in the bounding volume has an uniform chance to be proposed as a new live point.
"""
struct ENSUniformly <: ENSProposal end

"""
    struct BAT.ENSAutoProposal <: ENSProposal

*Experimental feature, not part of stable public API.*

Choose the proposal depending from the number of dimensions:
     ndims < 10: Proposals.Uniform,
     10 ≤ ndims ≤ 20: Proposals.RWalk,
     ndims > 20: Proposals.Slice
"""
struct ENSAutoProposal <: ENSProposal end

"""
    struct BAT.ENSRandomWalk <: ENSProposal

*Experimental feature, not part of stable public API.*

New live point is proposed by using a random walk away from an existing live point.

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:
$(TYPEDFIELDS)
"""
@with_kw struct ENSRandomWalk <: ENSProposal

    "Acceptance ratio for the random walk."
    ratio::Float64 = 0.5

    "Minimum number of random walk steps."
    walks::Int64 = 25

    "Scale of the proposal distribution."
    scale::Float64 = 1.0 # >= 0
end

"""
    struct BAT.ENSSlice <: ENSProposal

*Experimental feature, not part of stable public API.*

New live point is proposed by a serie of random slices from an existing live-point.

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)
"""
@with_kw struct ENSSlice <: ENSProposal

    "Minimum number of slices"
    slices::Int64 = 5

    "Scale of the proposal distribution."
    scale::Float64 = 1.0 # >= 0
end

function ENSprop(prop::ENSUniformly)
    return Proposals.Uniform()
end

function ENSprop(prop::ENSAutoProposal)
    return :auto     # :auto declaration: ndims < 10: Proposals.Uniform, 10 ≤ ndims ≤ 20: Proposals.RWalk, ndims > 20: Proposals.Slice
end

function ENSprop(prop::ENSRandomWalk)
    return Proposals.RWalk(;ratio=prop.ratio, walks=prop.walks, scale=prop.scale)
end

function ENSprop(prop::ENSSlice)
    return Proposals.Slice(;slices=prop.slices, scale=prop.scale)
end

function ENSprop(prop::ENSProposal) # if nothing is choosen
    return :auto
end


text:
function J2DifferentialEffect(arg0::Orbit, arg1::Orbit, arg2::jboolean, arg3::UnnormalizedSphericalHarmonicsProvider)
    return J2DifferentialEffect((Orbit, Orbit, jboolean, UnnormalizedSphericalHarmonicsProvider), arg0, arg1, arg2, arg3)
end

function J2DifferentialEffect(arg0::Orbit, arg1::Orbit, arg2::jboolean, arg3::jdouble, arg4::jdouble, arg5::jdouble)
    return J2DifferentialEffect((Orbit, Orbit, jboolean, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5)
end

function J2DifferentialEffect(arg0::SpacecraftState, arg1::AdapterPropagator_DifferentialEffect, arg2::jboolean, arg3::UnnormalizedSphericalHarmonicsProvider)
    return J2DifferentialEffect((SpacecraftState, AdapterPropagator_DifferentialEffect, jboolean, UnnormalizedSphericalHarmonicsProvider), arg0, arg1, arg2, arg3)
end

function J2DifferentialEffect(arg0::SpacecraftState, arg1::AdapterPropagator_DifferentialEffect, arg2::jboolean, arg3::jdouble, arg4::jdouble, arg5::jdouble)
    return J2DifferentialEffect((SpacecraftState, AdapterPropagator_DifferentialEffect, jboolean, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5)
end

function apply(obj::J2DifferentialEffect, arg0::Orbit)
    return jcall(obj, "apply", Orbit, (Orbit,), arg0)
end

function apply(obj::J2DifferentialEffect, arg0::SpacecraftState)
    return jcall(obj, "apply", SpacecraftState, (SpacecraftState,), arg0)
end


text:
import LinAlgExtensions: tolerance, tolreldefault, tolabsdefault
import Base: getproperty

srand(1)
let X = [spzeros(5,1) sprandn(5, 6, 0.5)], m = size(X, 1), n = size(X, 2)

stype(X) = X isa AbstractSparseMatrix ? "sparse" : "dense"

@testset "$(stype(A))" for A in (X, Matrix(X))
    @test normrows(A) ≈ sqrt.(diag(A*A'))
    @test normcols(A) ≈ sqrt.(diag(A'A))
end

ncx = maximum(normcols(X))

@test tolreldefault(1.0+0im) ≈ eps() * 20
@test tolabsdefault(X) ≈ ncx * eps() * n * 20

@test tolerance(X, true, 0.0, 0.0) == tolreldefault(X) * ncx
@test tolerance(X, true, 1e-5, 0.0)  == 1e-5 * ncx
@test tolerance(X, true, 0.0, 1e-5)  == 1e-5
@test tolerance(X, true, 1e-5, 1e-4)  == 1e-4
@test tolerance(X, true, 1e-5, 1e-20)  == 1e-5 * ncx
@test tolerance(X, false, 0.0, 0.0) == 0
@test tolerance(X, false, 1e-5, 0.0)  == 1e-5 * ncx
@test tolerance(X, false, 0.0, 1e-5)  == 1e-5
@test tolerance(X, false, 1e-5, 1e-4)  == 1e-4
@test tolerance(X, false, 1e-5, 1e-20)  == 1e-5 * ncx

@testset "$(stype(A))" for A in (X, Matrix(X))
    @test adjustsize(A, 5, 7) === A
    @test adjustsize(A, 7, 7) == [A; zeros(2, 7)]
    @test adjustsize(A, 5, 9) == [A zeros(5, 2)]
    @test adjustsize(A, 7, 10) == [A zeros(5, 3); zeros(2, 10)]
    @test adjustsize(A, 4, 10) == [A[1:4,:] zeros(4, 3)]
    @test adjustsize(A, 8, 2) == [A[:,1:2]; zeros(3, 2)]
    @test adjustsize(A, 4, 5) == A[1:4,1:5]
end

@testset "$(stype(A)) pivot=$piv" for A in (X, Matrix(X)), piv in (true, false)
    qrf = qrfactors(A, pivot=piv)
    @test size(qrf) == size(A)
    @test size(qrf.Q) == (size(A, 1), size(A, 1))
    @test size(qrf.R) == (min(size(A)...), size(A, 2))
    @test isperm(qrf.pcol)
    @test length(qrf.pcol) == n
    @test isperm(qrf.prow)
    @test length(qrf.prow) == m
end

@test_throws ArgumentError qrfactors(Float16.(X))
qrf1 = qrfactors(X')
qrf2 = qrfactors(copy(X'))
@testset "adjoint $d" for d in (:Q, :R, :pcol, :prow)
    @test getproperty(qrf1, d) == getproperty(qrf2, d)
end
end


text:
# addition
function (+)(f::SimpleRationalFunction, g::SimpleRationalFunction)
    top = f.num * g.den + f.den * g.num
    bot = f.den * g.den
    return SimpleRationalFunction(top, bot)
end

function (+)(f::SimpleRationalFunction, p::SimplePolynomial)
    return f + SimpleRationalFunction(p)
end
(+)(p::SimplePolynomial, f::SimpleRationalFunction) = f + p

(+)(f::SimpleRationalFunction, a::T) where {T<:CoefX} = f + SimpleRationalFunction(a)
(+)(a::T, f::SimpleRationalFunction) where {T<:CoefX} = f + a

# subtraction
(-)(f::SimpleRationalFunction) = SimpleRationalFunction(-f.num, f.den)

(-)(f::SimpleRationalFunction, g::SimpleRationalFunction) = f + (-g)

(-)(f::SimpleRationalFunction, p::SimplePolynomial) = f + (-p)
(-)(p::SimplePolynomial, f::SimpleRationalFunction) = p + (-f)

(-)(f::SimpleRationalFunction, a::T) where {T<:CoefX} = f + SimpleRationalFunction(-a)
(-)(a::T, f::SimpleRationalFunction) where {T<:CoefX} = SimpleRationalFunction(a) - f

# multiplication
(*)(f::SimpleRationalFunction, g::SimpleRationalFunction) =
    SimpleRationalFunction(f.num * g.num, f.den * g.den)

(*)(f::SimpleRationalFunction, p::SimplePolynomial) = f * SimpleRationalFunction(p)
(*)(p::SimplePolynomial, f::SimpleRationalFunction) = f * p

(*)(f::SimpleRationalFunction, a::T) where {T<:CoefX} = f * SimpleRationalFunction(a)
(*)(a::T, f::SimpleRationalFunction) where {T<:CoefX} = f * a

# division
inv(f::SimpleRationalFunction) = SimpleRationalFunction(f.den, f.num)

function inv(p::SimplePolynomial)
    T = eltype(p)
    SimpleRationalFunction(T(1), p)
end

SimpleThing = Union{SimplePolynomial,SimpleRationalFunction}

(/)(f::S, g::T) where {S<:SimpleThing,T<:SimpleThing} = f * inv(g)
(/)(f::SimpleRationalFunction, a::T) where {T<:CoefX} = f / SimplePolynomial(a)
(/)(a::T, f::S) where {T<:CoefX,S<:SimpleThing} = a * inv(f)

(//)(f::S, g::T) where {S<:SimpleThing,T<:SimpleThing} = f * inv(g)
(//)(f::SimpleRationalFunction, a::T) where {T<:CoefX} = f / SimplePolynomial(a)
(//)(a::T, f::S) where {T<:CoefX,S<:SimpleThing} = a * inv(f)

# exponentiation
function (^)(p::SimpleRationalFunction, k::S) where {S<:Integer}
    if k == 0
        return SimpleRationalFunction(1)
    end

    if k < 0
        return (inv(p))^(-k)
    end

    if k == 1
        return p
    end

    if k == 2
        return p * p
    end

    if k % 2 == 0  # even exponent
        j = k ÷ 2
        q = p^j
        return q * q
    end

    j = k ÷ 2
    q = p^j
    return q * q * p
end

# derivative
import SimplePolynomials: derivative
import Base: adjoint

function derivative(f::SimpleRationalFunction)
    p = numerator(f)
    q = denominator(f)
    return (p' * q - q' * p) / (q * q)
end
adjoint(f::SimpleRationalFunction) = derivative(f)


text:
function left_side(side)
    return HorizonSide(mod(Int(side) + 1 + 4, 4))
end
function right_side(side)
    return HorizonSide(mod(Int(side) - 1 + 4, 4))
end
function inverse(side)
    return HorizonSide(mod(Int(side) + 2, 4))
end
function mark_line(r, side, dif)
    res = 0
    while !isborder(r, side)
        res += dif
        move!(r, side)
        putmarker!(r)
    end
    return res
end

function get_to_start(r, side, x, y)
    for i in (1:x)
        move!(r, inverse(side))
    end
    y = abs(y)
    for i in (1:y)
        move!(r, inverse(right_side(side)))
    end
end

function complete_side(r, side)
    x = 0
    while(!isborder(r, side))
        x += 1
        move!(r, side)
    end
    putmarker!(r)
    y = 0
    y += mark_line(r, left_side(side), 1)
    y += mark_line(r, right_side(side), -1)
    get_to_start(r, side, x, y)
end

for i in (0:3)
    complete_side(r, HorizonSide(i))
end
show(r)


text:
"""
    LennardJones{T<:AbstractFloat} <: EmpiricalPotential{NamedTuple{(:ϵ, :σ)},NamedTuple{(:rcutoff,)}}

The Lennard-Jones potential is a simple two-body intermolecular potential with two typical parameters, primarily describing the interaction of neutral atoms. Formally, the interaction between two atoms at a distance, ``r``, is given by

```math
\\begin{equation}
V_{LJ}(r; \\epsilon, \\sigma, rcutoff, species) =
    \\begin{cases}
        0 & r > rcutoff \\\\
        4\\epsilon \\lbrack \\frac{\\sigma}{r})^{12} - (\\frac{\\sigma}{r})^6 \\rbrack & r < rcutoff.
    \\end{cases}
\\end{equation}
```

Users must supply two parameters, ``\\epsilon`` (units of energy), ``\\sigma`` (units of distance), and radial cutoff (units of distance).
"""
struct LennardJones{T<:AbstractFloat} <: EmpiricalPotential{NamedTuple{(:ϵ, :σ)},NamedTuple{(:rcutoff,)}}
    ϵ::T
    σ::T
    rcutoff::T
    species::Tuple
    function LennardJones(ϵ, σ, rcutoff, species)
        ϵ, σ, rcutoff = promote(austrip(ϵ), austrip(σ), austrip(rcutoff))
        new{typeof(rcutoff)}(ϵ, σ, rcutoff, Tuple(species))
    end
    function LennardJones{T}(; ϵ::T, σ::T, rcutoff::T, species::Tuple) where {T}
        new{T}(ϵ, σ, rcutoff, species)
    end
end

get_rcutoff(lj::LennardJones) = lj.rcutoff
get_species(lj::LennardJones) = lj.species

potential_energy(R::AbstractFloat, lj::LennardJones) = 4lj.ϵ * ((lj.σ / R)^12 - (lj.σ / R)^6)
force(R::AbstractFloat, r::SVector{3}, lj::LennardJones) = (24lj.ϵ * (2(lj.σ / R)^12 - (lj.σ / R)^6) / R^2)r


text:
module PlusPlus

export (++)

++(args::Union{AbstractString,AbstractChar}...) = *(args...)

++(args...) = vcat(args...)

end


text:
using Test
using Logging
logger = SimpleLogger(stderr, Logging.Info)
global_logger(logger)

import AoC2019
using OffsetArrays
input_data = AoC2019.read_numbers("input19.txt")

function check_beam(input_data)

  dims = (20, 20)
  beam_map = OffsetArray(zeros(Int8, dims...), 0:dims[1]-1, 0:dims[2]-1)
  println("Let's start")

  min_x = OffsetVector(zeros(Int32, 1000), 0:999)
  max_x = OffsetVector(zeros(Int32, 1000), 0:999)
  max_x[0] = 1
  @test AoC2019.run_program(input_data, input=[0,0]).output[1] == 1

  for x in 0:dims[1]-1
    for y in 0:dims[2]-1
      beam_map[y, x] = AoC2019.run_program(input_data, input=[x,y]).output[1]
    end
  end

  y = 5
  min_x[4] = findfirst(x -> x == 1, beam_map[y, :])
  max_x[4] = findlast(x -> x == 1, beam_map[y, :]) + 1
  while y < 1000
    println(stderr, y)

    current_left_border = min_x[y-1]
    current_right_border = max_x[y-1]
    println(stderr, "$current_left_border, $current_right_border")
    while AoC2019.run_program(input_data, input=[current_left_border,y]).output[1] == 0
      current_left_border += 1
    end
    min_x[y] = current_left_border
    while AoC2019.run_program(input_data, input=[current_right_border,y]).output[1] == 1
      println(stderr, "right border: $current_right_border")
      current_right_border += 1
    end
    max_x[y] = current_right_border

    #beam_map[y, current_left_border:current_right_border] .= 1
    #plt = heatmap(beam_map, yflip=true)
    #gui(plt)
    if y > 100
      println(stderr, "left border: $current_left_border")
      println(stderr, "Old y: $(y-99)")
      if max_x[y-99] - min_x[y] >= 100
        println("min_x: $(min_x[y]), max_x: $(max_x[y-99] - 1), min_y: $(y-99), max_y: $y")
        @assert AoC2019.run_program(input_data, input=[min_x[y],y-99]).output[1] == 1
        @assert AoC2019.run_program(input_data, input=[max_x[y-99]-1,y]).output[1] == 1
        @assert AoC2019.run_program(input_data, input=[max_x[y-99]-1,y-99]).output[1] == 1
        @assert AoC2019.run_program(input_data, input=[min_x[y],y]).output[1] == 1
        return (min_x[y] * 10000 + (y-99))
        break
      end
    end
    y += 1
  end
  return y
end
solution = check_beam(input_data)

max_stars_repo_name: buotex/aoc

max_issues_repo_name: buotex/aoc

max_forks_repo_name: buotex/aoc


text:
module BenchFib

import Schedulers
using Schedulers.Internal: WorkStealingScheduler, WorkerStates
using BenchmarkTools

module SchedulersFib
using Schedulers
function fib(n::Int)
    n <= 1 && return n
    t = Schedulers.spawn() do
        fib(n - 2)
    end
    return fib(n - 1) + fetch(t)::Int
end
end # module SchedulersFib

module BaseFib
function fib(n::Int)
    n <= 1 && return n
    t = Threads.@spawn fib(n - 2)
    return fib(n - 1) + fetch(t)::Int
end
end # module BaseFib

function unsafe_gc!(
    scheduler::WorkStealingScheduler = Schedulers.Internal.DEFAULT_WORKSTEALING_SCHEDULER_HANDLE[],
)
    nspins = 0
    for worker in scheduler.workers
        while (@atomic :monotonic worker.state) != WorkerStates.WAITING
            yield()
            nspins += 1
            if nspins > 1_000_000
                return
            end
        end
    end
    for deque in scheduler.queues
        data = deque.buffer.data
        n = length(data)
        empty!(data)
        resize!(data, n)
    end
end

function setup(; ns = [20])
    suite = BenchmarkGroup()
    for n in ns
        sn = suite["n=$n"] = BenchmarkGroup()
        sn["scheduler=:workstealing"] = @benchmarkable(
            Schedulers.global_workstealing() do
                SchedulersFib.fib($n)
            end,
            # Calling GC manually to avoid hitting copy-stack
            setup = begin
                unsafe_gc!()
                GC.gc(false)
            end
        )
        sn["scheduler=:prioritized"] = @benchmarkable(
            Schedulers.global_prioritized() do
                SchedulersFib.fib($n)
            end,
            # Calling GC manually to avoid hitting copy-stack
            setup = GC.gc(false)
        )
        sn["scheduler=:base"] = @benchmarkable(
            BaseFib.fib($n),
            # To be equivalent to what other benchmarks do:
            setup = GC.gc(false),
        )
    end
    return suite
end

function clear() end

end  # module


text:
using PlanerTransformer
using Plots

# define pcb
trace_edge_gap = 0.16e-3
trace_trace_gap = 0.13e-3
outer = copper_weight_to_meters(1.0)
inner = copper_weight_to_meters(1.0)
number_of_layers = 4
dielectric = (1.6e-3-2*outer-2*inner)/(number_of_layers-1)

stackup = Stackup([copper,fr408,copper,fr408,copper,fr408,copper],
                  [outer,dielectric,inner,dielectric,inner,dielectric,outer])
my_pcb = PCB_Specification(trace_edge_gap,
                        trace_trace_gap,
                        stackup)
turns_per_layer = 1:4
my_windings = [windings(my_pcb, core_geometry_dict["e32_plt"],
                        x,x,prisec"P-S-S-P",true,true) for x in turns_per_layer]
material = ferrite_dict["3f4"]
my_transformers = Transformer.(material,my_windings)

frequency = linspace(fmin(material),fmax(material),100)
v_in = 20.0
i_out = -0.5
my_converters = PushPull.(v_in, i_out, frequency)

my_analysis = [transformer_power_analysis.(t,c) for t in my_transformers, c in my_converters]
e = efficiency.(my_analysis)

p = plot()
xlabel!(p,"frequency")
ylabel!(p,"efficency")
for n in turns_per_layer
    plot!(p,frequency,e[n,:],label="$n turns per layer")
end
display(p)


text:
###
### Visualize the prior densities
###

cn = ConstrainedNormal(0.0,2.0,0.0,1e-5)

# Sample from ConstrainedNormal
vals = [inverse_parameter_transform(cn, x) for x in rand(convert_prior(cn), 10000000)]

f = Figure()
axtop = Axis(f[1, 1])
density!(axtop, vals)
# xlims!(0,2e-5)
save("visualize_ConstrainedNormal.png", f)
display(f)

density!(axright, ensemble[:, 2], direction = :y)

vlines!(axmain, [convective_κz], color=:red)
vlines!(axtop, [convective_κz], color=:red)
hlines!(axmain, [background_κz], color=:red)
hlines!(axright, [background_κz], color=:red)
colsize!(f.layout, 1, Fixed(300))
colsize!(f.layout, 2, Fixed(200))
rowsize!(f.layout, 1, Fixed(200))
rowsize!(f.layout, 2, Fixed(300))
leg = Legend(f[1, 2], scatters, ["Initial ensemble", "Iteration 1", "Iteration 2", "Iteration 10"], position = :lb)
hidedecorations!(axtop, grid = false)
hidedecorations!(axright, grid = false)
save("distributions_makie.png", f)


text:
Base.@kwdef mutable struct Benchmark
    b::Vector{Float64} # benchmark weights
    c::Vector{Float64} # Climate costs per $ mv
    σ::Vector{Float64} # Volatility
    ρ::Matrix # Correlation matrix
    Σ::Union{Matrix,Nothing} = nothing # Covariance Matrix
end

ρ = [
    1 0.8 0.7 0.6 0.7 0.5 0.7 0.6
    0.8 1 0.75 0.65 0.5 0.6 0.5 0.65
    0.7 0.75 1 0.8 0.7 0.7 0.7 0.7
    0.6 0.65 0.8 1 0.85 0.8 0.75 0.75
    0.7 0.5 0.7 0.85 1 0.6 0.8 0.65
    0.5 0.6 0.7 0.8 0.6 1 0.5 0.7
    0.7 0.5 0.7 0.75 0.8 0.5 1 0.8
    0.6 0.65 0.7 0.75 0.65 0.7 0.8 1
]

our_benchmark = Benchmark(b = [0.23, 0.19, 0.17, 0.13, 0.09, 0.08, 0.06, 0.05],
                        c= [0.06, 0.8, 0.39, 0.01, 0.05, 0.13, 0.5, 0.45],
                        σ = [0.22, 0.20, 0.25, 0.18, 0.35, 0.23, 0.13, 0.29],
                        ρ = ρ)

function get_benchmark_climate_costs(bench::Benchmark, mv::Int)::Float64
    return - round(bench.b' * bench.c * mv; digits= 2)
end

get_benchmark_climate_costs(our_benchmark, 1000)

Base.@kwdef mutable struct TiltedPortfolio
    x::Vector{Float64} # Optimal weights
end

other_portfolio = TiltedPortfolio(x = [0.30, 0.10, 0.08, 0.02, 0.48, 0.01, 0.005, 0.0055])

function get_climate_benefits(portfolio::TiltedPortfolio, bench::Benchmark, mv::Int)::Float64
    return -1 * round((portfolio.x - bench.b)' * bench.c * mv;digits = 5)
end

get_climate_benefits(other_portfolio, our_benchmark,1000)

check = other_portfolio.x' * our_benchmark.c * 1000

using LinearAlgebra
function get_cov_from_corr!(bench::Benchmark)::Benchmark
    D = diagm(bench.σ) # diagonal matrix of volatiltiy
    bench.Σ = D * bench.ρ * D # multiplying columns and rows of correlation matrix by volatiltiy
    return bench
end

get_cov_from_corr!(our_benchmark)

function tracking_error_volatility(portfolio::TiltedPortfolio, bench::Benchmark)::Float64
    return sqrt((portfolio.x - bench.b)' * bench.Σ * (portfolio.x - bench.b))
end

# get_excess_score(portfolio, our_benchmark)

# tracking_error_volatility(portfolio, our_benchmark)

using JuMP, COSMO # for the optimization

# tilt the benchmark according to gamma
function portfolio_tilting(bench::Benchmark, γ::Float64)::TiltedPortfolio
    n = length(bench.c) # number of assets
    model = JuMP.Model(COSMO.Optimizer)
    # the optimal weights we want to find
    @variable(model, x[1:n])
    @objective(model, Min, 1/2 * x' * bench.Σ * x - x' * (γ * -1 * bench.c + bench.Σ * bench.b))
    @constraint(model, zeros(n) .<= x .<= ones(n))
    @constraint(model, ones(n)' * x == 1)
    JuMP.optimize!(model)
    portfolio = TiltedPortfolio(x = JuMP.value.(x))
    return portfolio
end

# simulate for gammas between 0 and 10
using Plots
gammas = [i for i in 0.0:0.001:0.012]
climate_benefits = zeros(length(gammas))
excess_te = zeros(length(gammas))

for i in 1:length(gammas)
    new_portfolio = portfolio_tilting(our_benchmark, gammas[i])
    climate_benefits[i] = get_climate_benefits(new_portfolio, our_benchmark, 1000) # suppose 1000 billion of MV for the benchmark
    excess_te[i] = tracking_error_volatility(new_portfolio, our_benchmark) * 100
end

plot(excess_te, climate_benefits, xlabel = "Tracking error volatility (in %)",
                                ylabel = "Climate benefits for society (in billion USD)",
                                label = "")

max_stars_repo_name: tlorans/SustainableFinance

max_issues_repo_name: tlorans/SustainableFinance

max_forks_repo_name: tlorans/SustainableFinance


text:
using PyPlot

mutable struct point
    x::Float64
    y::Float64
end

function inCone(pOrg, pDeb, pFin, pCur)
    # pOrg : point origine du cone (la ou il est pointe)
    # pDeb : point depart du cone (point du rayon [pOrg,pDeb])
    # pFin : point final du cone (point du rayon [pOrg,pFin])
    # pCur : point courant a tester
    # retourne VRAI si pCur est dans le cone pDeb-pFin-pOrg, FAUX sinon

    cp_pDeb_pFin = (pDeb.x - pOrg.x) * (pFin.y - pOrg.y) - (pDeb.y - pOrg.y) * (pFin.x - pOrg.x)
    cp_pDeb_pCur = (pDeb.x - pOrg.x) * (pCur.y - pOrg.y) - (pDeb.y - pOrg.y) * (pCur.x - pOrg.x)
    cp_pFin_pCur = (pFin.x - pOrg.x) * (pCur.y - pOrg.y) - (pFin.y - pOrg.y) * (pCur.x - pOrg.x)

    if (cp_pDeb_pFin > 0)
        if ((cp_pDeb_pCur >= 0) && (cp_pFin_pCur <= 0))
            return true
        else
            return false
        end
    else
        if (!((cp_pDeb_pCur < 0) && (cp_pFin_pCur > 0)))
            return true
        else
            return false
        end
    end
end

function inCone1VersZ(pOrg, pDeb, pFin, pCur)
    return inCone(pOrg, pDeb, pFin, pCur)
end

function inCone2Vers0(pOrg, pDeb, pFin, pCur)
    return !inCone(pOrg, pDeb, pFin, pCur)
end

xlim(0,10)
ylim(0,10)
B=point(2.0,5.0); scatter(2.0,5.0,marker=".",color="black")
O=point(2.5,2.5); scatter(2.5,2.5,marker="x",color="red")
A=point(5.0,2.0); scatter(5.0,2.0,marker=".",color="black")
#N=point(5.0,5.0); scatter(5.0,5.0,marker="+")

x=round(rand()*8+2;digits=2); y=round(rand()*8+2;digits=2)
N=point(x,y); scatter(x,y,marker="x",color="blue")

x=round(rand()*10;digits=2); y=round(rand()*10;digits=2)
M=point(x,y); scatter(x,y,marker="o",color="magenta")

lx1=[O.x,B.x]; ly1=[O.y,B.y]; plot(lx1,ly1,color="red")
lx1=[O.x,0]; ly1=[O.y,(B.y-O.y)/(B.x-O.x)*(0-O.x)+O.y]; plot(lx1,ly1,color="red",linestyle="dotted")
lx2=[O.x,A.x]; ly2=[O.y,A.y]; plot(lx2,ly2,color="red")
lx2=[O.x,10]; ly2=[O.y,(A.y-O.y)/(A.x-O.x)*(10-O.x)+O.y]; plot(lx2,ly2,color="red",linestyle="dotted")

lx3=[N.x,B.x]; ly3=[N.y,B.y]; plot(lx3,ly3,color="blue",linestyle="dashed")
lx4=[N.x,A.x]; ly4=[N.y,A.y]; plot(lx4,ly4,color="blue",linestyle="dashed")

@show inCone1VersZ(O, A, B, M)
@show inCone2Vers0(N, A, B, M)
nothing

max_stars_repo_name: ducanh-le/GravityMachine

max_issues_repo_name: ducanh-le/GravityMachine

max_forks_repo_name: ducanh-le/GravityMachine


text:
Level1 = true
Level2 = true
Level3 = true
BenchNode = Dict{String,Vector{Float64}}()

function set_benchmark_levels(level1, level2, level3)
    global Level1 = level1
    global Level2 = level2
    global Level3 = level3
end

#=
macro level_bench(args...)
    return _level_bench(args[1], esc(args[2]))
end

function _level_bench(arg1, arg2)
    if eval(arg1) == true
        t0 = time_ns()
        val = $(arg2)
        println(val)
        t1 = time_ns()
        name = string(arg2.args[1])
        if haskey(BenchNode, name) != true
            BenchNode[name] = [(t1-t0)/1e9]
        else
            push!(BenchNode[name], (t1-t0)/1e9)
        end
        println("elapsed time in ", name, " ", (t1-t0)/1e9, " seconds")
        return val
    else
        return arg2
    end
end
=#

macro level_bench(arg1, arg2)
    name = string(arg2)
    quote
        while false; end
        if $(esc(arg1)) == true
            local t0 = time_ns()
            val = $(esc(arg2))
            local t1 = time_ns()
            name = $(esc(name))
            if haskey(BenchNode, name) != true
                BenchNode[name] = [(t1-t0)/1e9]
            else
                push!(BenchNode[name], (t1-t0)/1e9)
            end
            #println("elapsed time in ", name, " ", (t1-t0)/1e9, " seconds")
            val
        else
            $(esc(arg2))
        end
    end
end


text:
# This file was generated, do not modify it. # hide
tree_model = DecisionTreeClassifier()
tree = machine(tree_model, X, y)

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
# cd("c:/git/Game2048/benchmarks/")

include("move_combine_then_compact.jl")
include("move_ifs.jl")
include("move-loop.jl")
include("move-loop2.jl")
include("move_lookup.jl")

# generate all possible combinations
const VALS = Int8.(1:10)

function gen_all_comb()
    all_comb = Iterators.product([VALS for i in 1:4]...) |> collect ;
    all_comb = [[v...] for v in all_comb];
    all_comb
end

function test(all_comb)
    for v_orig in all_comb
        v_ifs = move_ifs!(copy(v_orig))
        v_combine = move_combine_compact!(copy(v_orig))
        v_loop = move_loop!(copy(v_orig))
        v_loop2 = move_loop2!(copy(v_orig))
        v_lookup = move_lookup!(copy(v_orig))

        if v_combine != v_ifs
            println("wrong ifs")
            println(v_orig)
            println(v_combine)
            println(v_ifs)
        end

        if v_combine != v_loop
            println("wrong loop")
            println(v_orig)
            println(v_combine)
            println(v_loop)
        end

        # if v_combine != v_lookup
        #     println("wrong lookup")
        #     println(v_orig)
        #     println(v_combine)
        #     println(v_lookup)
        # end

        # if v_combine != v_loop2
        #     println("wrong loop2")
        #     println(v_orig)
        #     println(v_combine)
        #     println(v_loop2)
        # end
    end
end

@time test(gen_all_comb())

using BenchmarkTools
# CONCLUSION: they all roughly perform the same so prefer the simpler implementation

@benchmark move_combine_compact!.(all_comb)  setup=(all_comb = gen_all_comb())
# BenchmarkTools.Trial:
#   memory estimate:  78.25 KiB
#   allocs estimate:  3
#   --------------
#   minimum time:     52.100 μs (0.00% GC)
#   median time:      66.101 μs (0.00% GC)
#   mean time:        74.256 μs (3.74% GC)
#   maximum time:     4.943 ms (97.62% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1
@benchmark move_ifs!.(all_comb)  setup=(all_comb = gen_all_comb())
# BenchmarkTools.Trial:
#   memory estimate:  78.25 KiB
#   allocs estimate:  3
#   --------------
#   minimum time:     50.900 μs (0.00% GC)
#   median time:      65.799 μs (0.00% GC)
#   mean time:        74.563 μs (3.42% GC)
#   maximum time:     2.497 ms (93.35% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1
@benchmark move_loop!.(all_comb)  setup=(all_comb = gen_all_comb())
# BenchmarkTools.Trial:
#   memory estimate:  78.25 KiB
#   allocs estimate:  3
#   --------------
#   minimum time:     50.900 μs (0.00% GC)
#   median time:      65.799 μs (0.00% GC)
#   mean time:        74.563 μs (3.42% GC)
#   maximum time:     2.497 ms (93.35% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1


text:
module Cliquing

using Memento
using LinearAlgebra

export AbstractClique, Clique
export GreedyClique, greedycliquing

const LOGGER = getlogger(@__MODULE__)

__init__() = Memento.register(LOGGER)

include("clique_type.jl")
include("greedycliquing.jl")

end


text:
function readAnnotate!(annotationFile::String, snpID::Vector{String},
                       nSNP::Int, grpInfo::Vector{Int64},
                       offsetSize::Vector{Int64}, geneName::Vector{String})

  annotatedata = readdlm(annotationFile, ',');
  annotatedata = convert(Array{String, 2}, annotatedata);
  annotatesnp = annotatedata[:, 2];
  nAnno = size(annotatedata, 1);
  nGrp = 0;
  annoIdx = 1;
  pregene = "";
  curgene = "";

  for i = 1 : nSNP

    while annoIdx <= nAnno
      if annotatesnp[annoIdx] == snpID[i]
        curgene = annotatedata[annoIdx, 1];
        annoIdx += 1;
        break;
      else
        annoIdx += 1;
      end
    end

    if curgene == pregene
      grpInfo[nGrp] += 1;
    else
      nGrp += 1;
      grpInfo[nGrp] += 1;
      geneName[nGrp] = curgene;
      if nGrp > 1
        offsetSize[nGrp] = offsetSize[nGrp - 1] + grpInfo[nGrp - 1];
      end
      pregene = curgene;
    end

  end

end


text:
if !@isdefined(Femshop)
    include("../Femshop.jl");
    using .Femshop
end
init_femshop("FVheat2d");

useLog("FVheat2dlog")

# Configuration setup
domain(2)
solverType(FV)
timeStepper(RK4)

# Mesh
n = 20 # number of elements
mesh(QUADMESH, elsperdim=n)

# Variables and BCs
u = variable("u", SCALAR, CELL)
boundary(u, 1, FLUX, "0")

# Time interval and initial condition
T = 0.1;
timeInterval(T)
initial(u, "(sin(pi*x)*sin(2*pi*y))^4")

# The flux and source terms of the conservation equation
# F and S in the following equation:
# Dt(int(u dx)/A) = int(S dx) - int(F.n ds)
coefficient("D", 0.1) # Diffusion rate
fluxAndSource(u, "-D*dot(grad(u),normal())", "0")

#exportCode("fvheat2dcode") # uncomment to export generated code to a file
#importCode("fvheat2dcode") # uncomment to import code from a file

solve(u)

@finalize()

##### Uncomment below to plot
# xy = Femshop.fv_info.cellCenters

# using Plots
# pyplot();
# display(plot(xy[1,:], xy[2,:], u.values[:], st=:surface))


text:
#= test/linop.jl
=#

@testset "linop" begin
  n = 10

  @time for (t, herm) in ((LinearOperator, false), (HermitianLinearOperator, true))
    for T in (Float32, Float64, ComplexF32, ComplexF64)
      println("  $t/$T")

      A = rand(T, n, n)
      if herm
        A += A'
      end

      L = LinearOperator(A)
      @test isa(L, t)

      F = Matrix(L)

      @test A ≈ F
      @test F' ≈ Matrix(L')
      @test A ≈ L*Matrix(1.0I,n,n)
      @test A ≈ Matrix(1.0I,n,n)*L

      x = rand(T, n)
      @test A  *x ≈ L  *x
      @test A' *x ≈ L' *x
      @test transpose(A)*x ≈ transpose(L)*x

      B = rand(T, n, n)
      @test A  *B   ≈ L  *B
      @test A  *B'  ≈ L  *B'
      @test A  *transpose(B) ≈ L  *transpose(B)
      @test A' *B   ≈ L' *B
      @test A' *B'  ≈ L' *B'
      @test transpose(A)*B   ≈ transpose(L)*B
      @test transpose(A)*transpose(B) ≈ transpose(L)*transpose(B)
      @test B  *A   ≈ B  *L
      @test B  *A'  ≈ B  *L'
      @test B  *transpose(A) ≈ B  *transpose(L)
      @test B' *A   ≈ B' *L
      @test B' *A'  ≈ B' *L'
      @test transpose(B)*A   ≈ transpose(B)*L
      @test transpose(B)*transpose(A) ≈ transpose(B)*transpose(L)

      c = rand()
      @test (c*A) *x ≈ (c*L) *x
      @test (c*A)'*x ≈ (c*L)'*x
      @test (c\A) *x ≈ (c\L) *x
      @test (c\A)'*x ≈ (c\L)'*x
      @test (A/c) *x ≈ (L/c) *x
      @test (A/c)'*x ≈ (L/c)'*x

      M = c*L
      M = L + M
      @test isa(M, t)
      @test (1 + c)*(A*x) ≈ M*x

      M = c*L
      M = L - M
      @test isa(M, t)
      @test (1 - c)*(A*x) ≈ M*x

      M = L^2
      @test A*(A*x) ≈ M*x

      M = convert(LinearOperator{T}, L)
      @test A*x == M*x
      M = convert(LinearOperator, L)
      @test A*x == M*x
    end
  end
end


text:
function get_dimension(obj::Space)
    return jcall(obj, "getDimension", jint, ())
end

function get_sub_space(obj::Space)
    return jcall(obj, "getSubSpace", Space, ())
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module IOTests

using Test
using Dates

@testset "string/show representation of Date" begin
    @test string(Dates.Date(1, 1, 1)) == "0001-01-01" # January 1st, 1 AD/CE
    @test sprint(show, Dates.Date(1, 1, 1)) == "Dates.Date(\"0001-01-01\")"
    @test string(Dates.Date(0, 12, 31)) == "0000-12-31" # December 31, 1 BC/BCE
    @test Dates.Date(1, 1, 1) - Dates.Date(0, 12, 31) == Dates.Day(1)
    @test Dates.Date(Dates.UTD(-306)) == Dates.Date(0, 2, 29)
    @test string(Dates.Date(0, 1, 1)) == "0000-01-01" # January 1st, 1 BC/BCE
    @test string(Dates.Date(-1, 1, 1)) == "-0001-01-01" # January 1st, 2 BC/BCE
    @test string(Dates.Date(-1000000, 1, 1)) == "-1000000-01-01"
    @test string(Dates.Date(1000000, 1, 1)) == "1000000-01-01"
    @test string(Dates.DateTime(2000, 1, 1, 0, 0, 0, 1)) == "2000-01-01T00:00:00.001"
    @test sprint(show, Dates.DateTime(2000, 1, 1, 0, 0, 0, 1)) == "Dates.DateTime(\"2000-01-01T00:00:00.001\")"
    @test string(Dates.DateTime(2000, 1, 1, 0, 0, 0, 2)) == "2000-01-01T00:00:00.002"
    @test string(Dates.DateTime(2000, 1, 1, 0, 0, 0, 500)) == "2000-01-01T00:00:00.5"
    @test string(Dates.DateTime(2000, 1, 1, 0, 0, 0, 998)) == "2000-01-01T00:00:00.998"
    @test string(Dates.DateTime(2000, 1, 1, 0, 0, 0, 999)) == "2000-01-01T00:00:00.999"
end

@testset "string/show representation of Time" begin
    tests = [
        Dates.Time(0) => ("00:00:00", "Dates.Time(0)"),
        Dates.Time(0, 1) => ("00:01:00", "Dates.Time(0, 1)"),
        Dates.Time(0, 1, 2) => ("00:01:02", "Dates.Time(0, 1, 2)"),
        Dates.Time(0, 1, 2, 3) => ("00:01:02.003", "Dates.Time(0, 1, 2, 3)"),
        Dates.Time(0, 1, 2, 3, 4) => ("00:01:02.003004", "Dates.Time(0, 1, 2, 3, 4)"),
        Dates.Time(0, 1, 2, 3, 4, 5) => ("00:01:02.003004005", "Dates.Time(0, 1, 2, 3, 4, 5)"),
        Dates.Time(0, 0, 0, 0, 0, 1) => ("00:00:00.000000001", "Dates.Time(0, 0, 0, 0, 0, 1)"),
        Dates.Time(0, 0, 0, 1) => ("00:00:00.001", "Dates.Time(0, 0, 0, 1)"),
    ]

    for (t, (printed, shown)) in tests
        @test sprint(print, t) == printed
        @test string(t) == printed
        @test sprint(show, t) == shown
        @test repr(t) == shown
    end
end

@testset "DateTime parsing" begin
    # Useful reference for different locales: http://library.princeton.edu/departments/tsd/katmandu/reference/months.html

    # Allow parsing of strings which are not representable as a TimeType
    str = "02/15/1996 25:00"
    df = Dates.DateFormat("mm/dd/yyyy HH:MM")
    parsed = Any[
        Dates.Month(2), Dates.Day(15), Dates.Year(1996), Dates.Hour(25), Dates.Minute(0)
    ]
    @test Dates.parse_components(str, df) == parsed
    @test_throws ArgumentError Dates.parse(DateTime, str, df)
end

@testset "DateFormat printing" begin
    @test sprint(show, DateFormat("yyyzzxmmdd\\MHH:MM:SS\\P")) == "dateformat\"yyyzzxmmdd\\MHH:MM:SSP\""
    @test sprint(show, DateFormat("yyy").tokens[1]) == "DatePart(yyy)"
    @test sprint(show, DateFormat("mmzzdd").tokens[2]) == "Delim(zz)"
    @test sprint(show, DateFormat("ddxmm").tokens[2]) == "Delim(x)"
    @test sprint(show, DateFormat("xxmmxx").tokens[2]) == "DatePart(mm)"
end
@testset "Common Parsing Patterns" begin
    #'1996-January-15'
    dt = Dates.DateTime(1996, 1, 15)
    f = "yy-mm-dd"
    a = "96-01-15"
    @test DateTime(a, f) + Dates.Year(1900) == dt
    @test Dates.format(dt, f) == a
    a1 = "96-1-15"
    @test Dates.DateTime(a1, f) + Dates.Year(1900) == dt
    @test Dates.format(dt, "yy-m-dd") == a1
    a2 = "96-1-1"
    @test Dates.DateTime(a2, f) + Dates.Year(1900) + Dates.Day(14) == dt
    @test Dates.format(dt - Dates.Day(14), "yy-m-d") == a2
    a3 = "1996-1-15"
    @test Dates.DateTime(a3, f) == dt
    @test Dates.format(dt, "yyyy-m-d") == a3
    a4 = "1996-Jan-15"
    @test_throws ArgumentError Dates.DateTime(a4, f) # Trying to use month name, but specified only "mm"

    f = "yy/uuu/dd"
    b = "96/Feb/15"
    @test Dates.DateTime(b, f) + Dates.Year(1900) == dt + Dates.Month(1)
    @test Dates.format(dt + Dates.Month(1), f) == b
    b1 = "1996/Feb/15"
    @test Dates.DateTime(b1, f) == dt + Dates.Month(1)
    @test Dates.format(dt + Dates.Month(1), "yyyy/uuu/dd") == b1
    b2 = "96/Feb/1"
    @test Dates.DateTime(b2, f) + Dates.Year(1900) + Dates.Day(14) == dt + Dates.Month(1)
    @test Dates.format(dt + Dates.Month(1) - Dates.Day(14), "yy/uuu/d") == b2
    # Here we've specified a text month name, but given a number
    b3 = "96/2/15"
    @test_throws ArgumentError Dates.DateTime(b3, f)
    try
        Dates.parse(DateTime, "2012/2/20T9:9:31.25i90", dateformat"yyyy/mm/ddTHH:MM:SS.s")
        @test false
    catch err
        @test isa(err, ArgumentError)
        @test err.msg == "Found extra characters at the end of date time string"
    end
    try
        Dates.parse(DateTime, "2012/2/20T9:9:3i90", dateformat"yyyy/mm/ddTHH:MM:SS.s")
        @test false
    catch err
        @test isa(err, ArgumentError)
        @test err.msg == "Unable to parse date time. Expected directive Delim(.) at char 16"
    end

    f = "yy:dd:mm"
    c = "96:15:01"
    @test Dates.DateTime(c, f) + Dates.Year(1900) == dt
    @test Dates.format(dt, f) == c
    c1 = "1996:15:01"
    @test Dates.DateTime(c1, f) == dt
    @test Dates.format(dt, "yyyy:dd:mm") == c1
    c2 = "96:15:1"
    @test Dates.DateTime(c2, f) + Dates.Year(1900) == dt
    @test Dates.format(dt, "yy:dd:m") == c2
    c3 = "96:1:01"
    @test Dates.DateTime(c3, f) + Dates.Year(1900) + Dates.Day(14) == dt
    @test Dates.format(dt - Dates.Day(14), "yy:m:dd") == c3
    c4 = "1996:15:01 # random comment"
    @test_throws ArgumentError Dates.DateTime(c4, f)

    f = "yyyy, uuu, dd"
    d = "1996, Jan, 15"
    @test Dates.DateTime(d, f) == dt
    @test Dates.format(dt, f) == d
    d1 = "96, Jan, 15"
    @test Dates.DateTime(d1, f) + Dates.Year(1900) == dt
    @test Dates.format(dt, "yy, uuu, dd") == d1
    d2 = "1996, Jan, 1"
    @test Dates.DateTime(d2, f) + Dates.Day(14) == dt
    @test Dates.format(dt - Dates.Day(14), "yyyy, uuu, d") == d2
    d3 = "1996, 2, 15"
    @test_throws ArgumentError Dates.DateTime(d3, f)

    f = "yyyy.U.dd"
    e = "1996.January.15"
    @test Dates.DateTime(e, f) == dt
    @test Dates.format(dt, f) == e
    e1 = "96.January.15"
    @test Dates.DateTime(e1, f) + Dates.Year(1900) == dt
    @test Dates.format(dt, "yy.U.dd") == e1

    fo = "yyyy m dd"
    f = "1996 1 15"
    @test Dates.DateTime(f, fo) == dt
    @test Dates.format(dt, fo) == f
    f1 = "1996 01 15"
    @test Dates.DateTime(f1, fo) == dt
    @test Dates.format(dt, "yyyy mm dd") == f1
    f2 = "1996 1 1"
    @test Dates.DateTime(f2, fo) + Dates.Day(14) == dt
    @test Dates.format(dt - Dates.Day(14), "yyyy m d") == f2

    j = "1996-01-15"
    f = "yyyy-mm-dd zzz"
    @test Dates.DateTime(j, f) == dt
    @test Dates.format(dt, f) == j * " zzz"
    k = "1996-01-15 10:00:00"
    f = "yyyy-mm-dd HH:MM:SS zzz"
    @test Dates.DateTime(k, f) == dt + Dates.Hour(10)
    @test Dates.format(dt + Dates.Hour(10), f) == k * " zzz"
    l = "1996-01-15 10:10:10.25"
    f = "yyyy-mm-dd HH:MM:SS.ss zzz"
    @test Dates.DateTime(l, f) == dt + Dates.Hour(10) + Dates.Minute(10) + Dates.Second(10) + Dates.Millisecond(250)
    @test Dates.format(dt + Dates.Hour(10) + Dates.Minute(10) + Dates.Second(10) + Dates.Millisecond(250), f) == l * " zzz"

    r = "1/15/1996" # Excel
    f = "m/dd/yyyy"
    @test Dates.DateTime(r, f) == dt
    @test Dates.format(dt, f) == r
    s = "19960115"
    f = "yyyymmdd"
    @test Dates.DateTime(s, f) == dt
    @test Dates.format(dt, f) == s
    v = "1996-01-15 10:00:00"
    f = "yyyy-mm-dd HH:MM:SS"
    @test Dates.DateTime(v, f) == dt + Dates.Hour(10)
    @test Dates.format(dt + Dates.Hour(10), f) == v
    w = "1996-01-15T10:00:00"
    f = "yyyy-mm-ddTHH:MM:SS zzz"
    @test Dates.DateTime(w, f) == dt + Dates.Hour(10)
    @test Dates.format(dt + Dates.Hour(10), f) == w * " zzz"

    f = "yyyy/m"
    y = "1996/1"
    @test Dates.DateTime(y, f) == dt - Dates.Day(14)
    @test Dates.format(dt, f) == y
    y1 = "1996/1/15"
    @test_throws ArgumentError Dates.DateTime(y1, f)
    y2 = "96/1"
    @test Dates.DateTime(y2, f) + Dates.Year(1900) == dt - Dates.Day(14)
    @test Dates.format(dt, "yy/m") == y2

    f = "yyyy"
    z = "1996"
    @test Dates.DateTime(z, f) == dt - Dates.Day(14)
    @test Dates.format(dt, f) == z
    z1 = "1996-3"
    @test_throws ArgumentError Dates.DateTime(z1, f)
    z2 = "1996-3-1"
    @test_throws ArgumentError Dates.DateTime(z2, f)

    aa = "1/5/1996"
    f = "m/d/yyyy"
    @test Dates.DateTime(aa, f) == dt - Dates.Day(10)
    @test Dates.format(dt - Dates.Day(10), f) == aa
    bb = "5/1/1996"
    f = "d/m/yyyy"
    @test Dates.DateTime(bb, f) == dt - Dates.Day(10)
    @test Dates.format(dt - Dates.Day(10), f) == bb
    cc = "01151996"
    f = "mmddyyyy"
    @test Dates.DateTime(cc, f) == dt
    @test Dates.format(dt, f) == cc
    dd = "15011996"
    f = "ddmmyyyy"
    @test Dates.DateTime(dd, f) == dt
    @test Dates.format(dt, f) == dd
    ee = "01199615"
    f = "mmyyyydd"
    @test Dates.DateTime(ee, f) == dt
    @test Dates.format(dt, f) == ee
    ff = "1996-15-Jan"
    f = "yyyy-dd-uuu"
    @test Dates.DateTime(ff, f) == dt
    @test Dates.format(dt, f) == ff
    gg = "Jan-1996-15"
    f = "uuu-yyyy-dd"
    @test Dates.DateTime(gg, f) == dt
    @test Dates.format(dt, f) == gg
    hh = "1996#1#15"
    f = "yyyy#m#d"
    @test Dates.DateTime(hh, f) == dt
    @test Dates.format(dt, f) == hh

    f = "ymd"
    @test Dates.Date("111", f) == Dates.Date(1)
    @test Dates.Date("1", f) == Dates.Date(1)

    @test Dates.DateTime("20140529 120000", "yyyymmdd HHMMSS") == Dates.DateTime(2014, 5, 29, 12)

    @test Dates.Date(string(Dates.Date(dt))) == Dates.Date(dt)
    @test Dates.DateTime(string(dt)) == dt
end
@testset "prefix." begin
    s = "/1996/1/15"
    f = "/yyyy/m/d"
    dt = Dates.DateTime(1996, 1, 15)
    @test Dates.DateTime(s, f) == dt
    @test Dates.format(dt, f) == s
    @test_throws ArgumentError Dates.DateTime("1996/1/15", f)
end
@testset "French and Chinese" begin
    # from Jiahao
    @test Dates.Date("2009年12月01日", "yyyy年mm月dd日") == Dates.Date(2009, 12, 1)
    @test Dates.format(Dates.Date(2009, 12, 1), "yyyy年mm月dd日") == "2009年12月01日"
    @test Dates.Date("2009-12-01", "yyyy-mm-dd") == Dates.Date(2009, 12, 1)

    # French: from Milan
    f = "dd/mm/yyyy"
    f2 = "dd/mm/yy"
    @test Dates.Date("28/05/2014", f) == Dates.Date(2014, 5, 28)
    @test Dates.Date("28/05/14", f2) + Dates.Year(2000) == Dates.Date(2014, 5, 28)
end

@testset "Customizing locale" begin
    Dates.LOCALES["french"] = Dates.DateLocale(
        ["janvier", "février", "mars", "avril", "mai", "juin",
         "juillet", "août", "septembre", "octobre", "novembre", "décembre"],
        ["janv", "févr", "mars", "avril", "mai", "juin",
         "juil", "août", "sept", "oct", "nov", "déc"],
        ["lundi", "mardi", "mercredi", "jeudi", "vendredi", "samedi", "dimanche"],
        [""],
    )

    f = "dd uuuuu yyyy"
    @test Dates.Date("28 mai 2014", f; locale="french") == Dates.Date(2014, 5, 28)
    @test Dates.format(Dates.Date(2014, 5, 28), f; locale="french") == "28 mai 2014"
    @test Dates.Date("28 févr 2014", f; locale="french") == Dates.Date(2014, 2, 28)
    @test Dates.format(Dates.Date(2014, 2, 28), f; locale="french") == "28 févr 2014"
    @test Dates.Date("28 août 2014", f; locale="french") == Dates.Date(2014, 8, 28)
    @test Dates.format(Dates.Date(2014, 8, 28), f; locale="french") == "28 août 2014"
    @test Dates.Date("28 avril 2014", f; locale="french") == Dates.Date(2014, 4, 28)
    @test Dates.format(Dates.Date(2014, 4, 28), f; locale="french") == "28 avril 2014"

    f = "dd u yyyy"
    @test Dates.Date("28 avril 2014", f; locale="french") == Dates.Date(2014, 4, 28)
    f = "dduuuuyyyy"
    # parses 3 and 4 character month names
    @test Dates.Date("28mai2014", f; locale="french") == Dates.Date(2014, 5, 28)
    @test Dates.Date("28août2014", f; locale="french") == Dates.Date(2014, 8, 28)
    # doesn't parse month name greater than 4 chars
    @test_throws ArgumentError Dates.Date("28avril2014", f; locale="french")
end
@testset "year digits parsing" begin
    # From Tony Fong
    f = "dduuuyy"
    @test Dates.Date("01Dec09", f) + Dates.Year(2000) == Dates.Date(2009, 12, 1)
    @test Dates.format(Dates.Date(2009, 12, 1), f) == "01Dec09"
    f = "dduuuyyyy"
    @test Dates.Date("01Dec2009", f) == Dates.Date(2009, 12, 1)
    @test Dates.format(Dates.Date(2009, 12, 1), f) == "01Dec2009"
    f = "duy"
    globex = ["f", "g", "h", "j", "k", "m", "n", "q", "u", "v", "x", "z"]
    locale = Dates.DateLocale(globex, map(uppercase, globex), globex[1:7], globex[1:7])
    @test Dates.Date("1F4", f; locale=locale) + Dates.Year(2010) == Dates.Date(2014, 1, 1)
    @test Dates.format(Dates.Date(2014, 1, 1), f; locale=locale) == "1F4"

    # From Matt Bauman
    f = "yyyy-mm-ddTHH:MM:SS"
    @test Dates.DateTime("2014-05-28T16:46:04", f) == Dates.DateTime(2014, 5, 28, 16, 46, 04)
end

@testset "Error handling" begin
    # Specified mm/dd, but date string has day/mm
    @test_throws ArgumentError Dates.DateTime("18/05/2009", "mm/dd/yyyy")
    @test_throws ArgumentError Dates.DateTime("18/05/2009 16", "mm/dd/yyyy hh")
    # Used "mm" for months AND minutes
    @test_throws ArgumentError Dates.DateTime("18/05/2009 16:12", "mm/dd/yyyy hh:mm")
    # Date string has different delimiters than format string
    @test_throws ArgumentError Dates.DateTime("18:05:2009", "mm/dd/yyyy")

    f = "y m d"
    @test Dates.Date("1 1 1", f) == Dates.Date(1)
    @test Dates.Date("10000000000 1 1", f) == Dates.Date(10000000000)
    @test_throws ArgumentError Dates.Date("1 13 1", f)
    @test_throws ArgumentError Dates.Date("1 1 32", f)
    @test_throws ArgumentError Dates.Date(" 1 1 32", f)
    @test_throws ArgumentError Dates.Date("# 1 1 32", f)
    @test Dates.Date("1", f) == Dates.Date(1)
    @test Dates.Date("1 ", f) == Dates.Date(1)
    @test Dates.Date("1 2", f) == Dates.Date(1, 2)
    # can't find space delimiter (finds '/'), so fails
    @test_throws ArgumentError Dates.Date("2000/1", f)
end

@testset "formerly vectorized Date/DateTime/format methods" begin
    dr = ["2000-01-01", "2000-01-02", "2000-01-03", "2000-01-04", "2000-01-05",
          "2000-01-06", "2000-01-07", "2000-01-08", "2000-01-09", "2000-01-10"]
    dr2 = [Dates.Date(2000) : Dates.Day(1) : Dates.Date(2000, 1, 10);]
    @test Dates.Date.(dr) == dr2
    @test Dates.Date.(dr, dateformat"yyyy-mm-dd") == dr2
    @test Dates.DateTime.(dr) == Dates.DateTime.(dr2)
    @test Dates.DateTime.(dr, dateformat"yyyy-mm-dd") == Dates.DateTime.(dr2)

    @test Dates.format.(dr2, "yyyy-mm-dd") == dr

    @test typeof(Dates.Date.(dr)) == Array{Date, 1}
end
@testset "Issue 13" begin
    t = Dates.DateTime(1, 1, 1, 14, 51, 0, 118)
    @test Dates.DateTime("[14:51:00.118]", "[HH:MM:SS.sss]") == t
    @test Dates.DateTime("14:51:00.118", "HH:MM:SS.sss") == t
    @test Dates.DateTime("[14:51:00.118?", "[HH:MM:SS.sss?") == t
    @test Dates.DateTime("?14:51:00.118?", "?HH:MM:SS.sss?") == t
    @test Dates.DateTime("x14:51:00.118", "xHH:MM:SS.sss") == t
    @test Dates.DateTime("14:51:00.118]", "HH:MM:SS.sss]") == t
end
@testset "RFC1123Format" begin
    dt = Dates.DateTime(2014, 8, 23, 17, 22, 15)
    @test Dates.format(dt, Dates.RFC1123Format) == "Sat, 23 Aug 2014 17:22:15"
    @test Dates.DateTime(Dates.format(dt, Dates.RFC1123Format), Dates.RFC1123Format) == dt
    @test Dates.format(dt, "yyyy-mm-ddTHH:MM:SS E") == "2014-08-23T17:22:15 Saturday"
    @test Dates.format(dt, "yyyy-mm-ddTHH:MM:SS e") == "2014-08-23T17:22:15 Sat"
    @test Dates.format(dt, "yyyy-mm-dd E") == "2014-08-23 Saturday"
    @test Dates.format(dt, "yyyy-mm-dd e") == "2014-08-23 Sat"
    @test Dates.format(dt, "yyyy-e-mm-dd") == "2014-Sat-08-23"

    @test Dates.format(Dates.DateTime(2014, 1, 2, 0, 0, 0, 999), Dates.RFC1123Format) == "Thu, 02 Jan 2014 00:00:00"
    @test Dates.format(Dates.DateTime(2014, 2, 18, 0, 0, 0, 9), Dates.RFC1123Format) == "Tue, 18 Feb 2014 00:00:00"
    @test Dates.format(Dates.DateTime(2014, 3, 8, 0, 0, 0, 9), Dates.RFC1123Format) == "Sat, 08 Mar 2014 00:00:00"
    @test Dates.format(Dates.DateTime(2014, 4, 28, 0, 0, 0, 9), Dates.RFC1123Format) == "Mon, 28 Apr 2014 00:00:00"
    @test Dates.format(Dates.DateTime(2014, 5, 10, 0, 0, 0, 9), Dates.RFC1123Format) == "Sat, 10 May 2014 00:00:00"
    @test Dates.format(Dates.DateTime(2014, 6, 4, 0, 0, 0, 9), Dates.RFC1123Format) == "Wed, 04 Jun 2014 00:00:00"
    @test Dates.format(Dates.DateTime(2014, 7, 13, 0, 0, 0, 9), Dates.RFC1123Format) == "Sun, 13 Jul 2014 00:00:00"
    @test Dates.format(Dates.DateTime(2014, 8, 17, 0, 0, 0, 9), Dates.RFC1123Format) == "Sun, 17 Aug 2014 00:00:00"
    @test Dates.format(Dates.DateTime(2014, 9, 20, 0, 0, 0, 9), Dates.RFC1123Format) == "Sat, 20 Sep 2014 00:00:00"
    @test Dates.format(Dates.DateTime(2014, 10, 31, 0, 0, 0, 9), Dates.RFC1123Format) == "Fri, 31 Oct 2014 00:00:00"
    @test Dates.format(Dates.DateTime(2014, 11, 2, 0, 0, 0, 9), Dates.RFC1123Format) == "Sun, 02 Nov 2014 00:00:00"
    @test Dates.format(Dates.DateTime(2014, 12, 5, 0, 0, 0, 9), Dates.RFC1123Format) == "Fri, 05 Dec 2014 00:00:00"

    dt = Dates.DateTime(2016, 11, 12, 7, 45, 36)
    @test parse(Dates.DateTime, "Sat, 12 Nov 2016 07:45:36", Dates.RFC1123Format) == dt
    @test parse(Dates.DateTime, "Mon, 12 Nov 2016 07:45:36", Dates.RFC1123Format) == dt  # Wrong day of week
    @test_throws ArgumentError parse(Date, "Foo, 12 Nov 2016 07:45:36", Dates.RFC1123Format)
end
@testset "Issue 15195" begin
    f = "YY"
    @test Dates.format(Dates.Date(1999), f) == "1999"
    @test Dates.format(Dates.Date(9), f) == "09"
    @test Dates.format(typemax(Dates.Date), f) == "252522163911149"
end

@testset "TimeZones.jl Issue #19" begin
    Zulu = String

    function Dates.tryparsenext(d::Dates.DatePart{'Z'}, str, i, len)
        return Dates.tryparsenext_word(str, i, len, Dates.min_width(d), Dates.max_width(d))
    end

    str = "2015-07-24T05:38:19.591Z"
    dt = Dates.DateTime(2015, 7, 24, 5, 38, 19, 591)
    parsed = Any[
        Dates.Year(2015), Dates.Month(7), Dates.Day(24),
        Dates.Hour(5), Dates.Minute(38), Dates.Second(19), Dates.Millisecond(591)
    ]

    format = "yyyy-mm-ddTHH:MM:SS.sssZ"
    escaped_format = "yyyy-mm-dd\\THH:MM:SS.sss\\Z"

    # Typically 'Z' isn't treated as a specifier so it doesn't have to be escaped
    @test Dates.parse_components(str, Dates.DateFormat(format)) == parsed
    @test Dates.parse_components(str, Dates.DateFormat(escaped_format)) == parsed

    try
        # Make 'Z' into a specifier
        Dates.CONVERSION_SPECIFIERS['Z'] = Zulu
        Dates.CONVERSION_DEFAULTS[Zulu] = ""

        @test Dates.parse_components(str, Dates.DateFormat(format)) == [parsed; Zulu("Z")]
        @test Dates.parse_components(str, Dates.DateFormat(escaped_format)) == parsed
    finally
        delete!(Dates.CONVERSION_SPECIFIERS, 'Z')
        delete!(Dates.CONVERSION_DEFAULTS, Zulu)
    end

    # Ensure that the default behaviour has been restored
    @test Dates.parse_components(str, Dates.DateFormat(format)) == parsed
    @test Dates.parse_components(str, Dates.DateFormat(escaped_format)) == parsed
end

@testset "Issue 10817" begin
    @test Dates.Date("Apr 01 2014", "uuu dd yyyy") == Dates.Date(2014, 4, 1)
    @test_throws ArgumentError Dates.Date("Apr 01 xx 2014", "uuu dd zz yyyy")
    @test_throws ArgumentError Dates.Date("Apr 01 xx 2014", "uuu dd    yyyy")
end
@testset "Issue 21001" begin
    for (ms, str) in zip([0, 1, 20, 300, 450, 678], ["0", "001", "02", "3", "45", "678"])
        local dt = DateTime(2000, 1, 1, 0, 0, 0, ms)
        @test Dates.format(dt, "s") == str
        @test Dates.format(dt, "ss") == rpad(str, 2, '0')
        @test Dates.format(dt, "sss") == rpad(str, 3, '0')
        @test Dates.format(dt, "ssss") == rpad(str, 4, '0')
    end
end
# Issue #21504
@test tryparse(Dates.Date, "0-1000") === nothing

@testset "parse milliseconds, Issue #22100" begin
    @test Dates.DateTime("2017-Mar-17 00:00:00.0000", "y-u-d H:M:S.s") == Dates.DateTime(2017, 3, 17)
    @test Dates.parse_components(".1", Dates.DateFormat(".s")) == [Dates.Millisecond(100)]
    @test Dates.parse_components(".12", Dates.DateFormat(".s")) == [Dates.Millisecond(120)]
    @test Dates.parse_components(".123", Dates.DateFormat(".s")) == [Dates.Millisecond(123)]
    @test Dates.parse_components(".1230", Dates.DateFormat(".s")) == [Dates.Millisecond(123)]
    @test_throws InexactError Dates.parse_components(".1234", Dates.DateFormat(".s"))

    # Ensure that no overflow occurs when using Int32 literals: Int32(10)^10
    @test Dates.parse_components("." * rpad(999, 10, '0'), Dates.DateFormat(".s")) == [Dates.Millisecond(999)]
end

@testset "Time Parsing" begin
    let t
        time_tuple(t::Dates.Time) = (
            Dates.hour(t), Dates.minute(t), Dates.second(t),
            Dates.millisecond(t), Dates.microsecond(t), Dates.nanosecond(t)
        )

        ## default ISOTimeFormat
        t = Dates.Time("01")
        @test time_tuple(t) == (1, 0, 0, 0, 0, 0)
        t = Dates.Time("01:23")
        @test time_tuple(t) == (1, 23, 0, 0, 0, 0)
        t = Dates.Time("01:23:45")
        @test time_tuple(t) == (1, 23, 45, 0, 0, 0)
        t = Dates.Time("01:23:45.678")
        @test time_tuple(t) == (1, 23, 45, 678, 0, 0)

        ## string format
        t = Dates.Time("23:56:12.1", "HH:MM:SS.s")
        @test time_tuple(t) == (23, 56, 12, 100, 0, 0)

        ## precomputed DateFormat
        t = Dates.Time("04:09:45.012", DateFormat("HH:MM:SS.s"))
        @test time_tuple(t) == (4, 9, 45, 12, 0, 0)
        t = Dates.Time("21 07", DateFormat("HH MM"))
        @test time_tuple(t) == (21, 7, 0, 0, 0, 0)
        t = Dates.Time("4.02", DateFormat("H.MM"))
        @test time_tuple(t) == (4, 2, 0, 0, 0, 0)
        t = Dates.Time("1725", DateFormat("HHMM"))
        @test time_tuple(t) == (17, 25, 0, 0, 0, 0)

        ## exceptions
        @test_throws ArgumentError Dates.Time("24:00")  # invalid hours
        @test_throws ArgumentError Dates.Time("00:60")  # invalid minutes
        @test_throws ArgumentError Dates.Time("00:00:60")  # invalid seconds
        @test_throws ArgumentError Dates.Time("20:03:20", DateFormat("HH:MM"))  # too much precision
        @test_throws ArgumentError Dates.Time("10:33:51", DateFormat("YYYY-MM-DD HH:MM:SS"))  # Time can't hold year/month/day
    end
end

@testset "midnight" begin
    # issue #28203: 24:00 is a valid ISO 8601 time
    @test DateTime("2018-01-01 24:00","yyyy-mm-dd HH:MM") == DateTime("2018-01-02T00:00:00") ==
          DateTime(2018, 1, 1, 24) == DateTime(2018, 1, 2)
    @test_throws ArgumentError DateTime("2018-01-01 24:01","yyyy-mm-dd HH:MM")
    @test_throws ArgumentError DateTime(2018, 1, 1, 24, 0, 1)
    @test_throws ArgumentError DateTime(2018, 1, 1, 24, 0, 0, 1)
end

@testset "AM/PM" begin
    # get the current locale
    LC_TIME = 2
    time_locale = ccall(:setlocale, Cstring, (Cint, Cstring), LC_TIME, C_NULL)
    try
        # set the locale
        ccall(:setlocale, Cstring, (Cint, Cstring), LC_TIME, "C")

        for (t12,t24) in (("12:00am","00:00"), ("12:07am","00:07"), ("01:24AM","01:24"),
                        ("12:00pm","12:00"), ("12:15pm","12:15"), ("11:59PM","23:59"))
            d = DateTime("2018-01-01T$t24:00")
            t = Time("$t24:00")
            for HH in ("HH","II")
                @test DateTime("2018-01-01 $t12","yyyy-mm-dd $HH:MMp") == d
                @test Time("$t12","$HH:MMp") == t
            end
            tmstruct = Libc.strptime("%I:%M%p", t12)
            @test Time(tmstruct) == t
            @test uppercase(t12) == Dates.format(t, "II:MMp") ==
                                    Dates.format(d, "II:MMp") ==
                Libc.strftime("%I:%M%p", tmstruct)
        end
        for bad in ("00:24am", "00:24pm", "13:24pm", "2pm", "12:24p.m.", "12:24 pm", "12:24pµ")
            @eval @test_throws ArgumentError Time($bad, "II:MMp")
        end
        # if am/pm is missing, defaults to 24-hour clock
        @eval Time("13:24", "II:MMp") == Time("13:24", "HH:MM")
    finally
        # recover the locale
        ccall(:setlocale, Cstring, (Cint, Cstring), LC_TIME, time_locale)
    end
end

end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
"""
    ϵ-Greedy Implementation
"""

type epsGreedy <: BanditAlgorithmBase
    noOfArms::Int64
    noOfSteps::Int64
    lastPlayedArm::Int64

    ϵ::Float64
    cummReward::Vector{Float64}
    count::Vector{Int64}
    avgValue::Vector{Float64}

    function epsGreedy( noOfArms::Int64, ϵ::Float64 )
        new( noOfArms,
             0,
             0,
             ϵ,
             zeros(Float64,noOfArms),
             zeros(Int64,noOfArms),
             zeros(Float64,noOfArms)
        )
    end
end

function get_arm_index( agent::epsGreedy )
    # if any(agent.count.==0)
    #     agent.lastPlayedArm = rand( find(agent.count.==0) )
    # else
        if rand() > agent.ϵ
            agent.lastPlayedArm = findmax(agent.avgValue)[2]
        else
            agent.lastPlayedArm = rand(1:agent.noOfArms)
        end
    # end
    return agent.lastPlayedArm
end

function update_reward!( agent::epsGreedy, r::Real )
    agent.cummReward[agent.lastPlayedArm] += r
    agent.count[agent.lastPlayedArm] += 1
    agent.noOfSteps += 1

    agent.avgValue[agent.lastPlayedArm] = agent.cummReward[agent.lastPlayedArm] ./
                                            agent.count[agent.lastPlayedArm]

    nothing
end

function reset!( agent::epsGreedy )
    agent.noOfSteps     = 0
    agent.lastPlayedArm = 0

    agent.cummReward    = zeros( Float64, agent.noOfArms )
    agent.count         = zeros( Int64, agent.noOfArms )
    agent.avgValue      = zeros( Float64, agent.noOfArms )

    nothing
end

function info_str( agent::epsGreedy, latex::Bool )
    if latex
        return @sprintf( "\$\\epsilon\$-Greedy (\$\\epsilon = %4.3f\$)", agent.ϵ )
    else
        return @sprintf( "ϵ-Greedy (ϵ = %4.3f)", agent.ϵ )
    end
end

"""
    ϵ_n Greedy Implementation
    Based on Auer, P., Bianchi, N. C., & Fischer, P. (2002). Finite time analysis of the multiarmed bandit problem. Machine Learning, 47, 235–256.
"""

type epsNGreedy <: BanditAlgorithmBase
    noOfArms::Int64
    noOfSteps::Int64
    lastPlayedArm::Int64

    param_c::Float64
    param_d::Float64
    ϵ::Float64

    cummReward::Vector{Float64}
    count::Vector{Int64}
    avgValue::Vector{Float64}

    function epsNGreedy( noOfArms::Int64, param_c::Real, param_d::Float64 )
        new( noOfArms,
             0,
             0,
             param_c,
             param_d,
             1,
             zeros(Float64,noOfArms),
             zeros(Int64,noOfArms),
             zeros(Float64,noOfArms)
        )
    end

    function epsNGreedy( noOfArms::Int64 )
        new( noOfArms,
             0,
             0,
             1/noOfArms,
             1,
             1,
             zeros(Float64,noOfArms),
             zeros(Int64,noOfArms),
             zeros(Float64,noOfArms)
        )
    end
end

function get_arm_index( agent::epsNGreedy )
    # if any(agent.count.==0)
    #     agent.lastPlayedArm = rand( find(agent.count.==0) )
    # else
        if rand() > agent.ϵ
            agent.lastPlayedArm = findmax(agent.avgValue)[2]
        else
            agent.lastPlayedArm = rand(1:agent.noOfArms)
        end
    # end
    return agent.lastPlayedArm
end

function update_reward!( agent::epsNGreedy, r::Real )
    # Book keeping
    agent.cummReward[agent.lastPlayedArm] += r
    agent.count[agent.lastPlayedArm] += 1
    agent.noOfSteps += 1

    # Update the observed reward to the corresponding arm
    agent.avgValue[agent.lastPlayedArm] = agent.cummReward[agent.lastPlayedArm] ./
                                            agent.count[agent.lastPlayedArm]

    # Also change the exploration rate
    agent.ϵ     = min( 1, (agent.param_c*agent.noOfArms)/(agent.param_d*agent.noOfSteps) )

    nothing
end

function reset!( agent::epsNGreedy )
    agent.noOfSteps     = 0
    agent.lastPlayedArm = 0
    agent.ϵ             = 1

    agent.cummReward    = zeros( Float64, agent.noOfArms )
    agent.count         = zeros( Int64, agent.noOfArms )
    agent.avgValue      = zeros( Float64, agent.noOfArms )

    nothing
end

function info_str( agent::epsNGreedy, latex::Bool )
    if latex
        return @sprintf( "\$\\epsilon_n\$-Greedy (\$c = %4.3f, d = %4.3f\$)", agent.param_c, agent.param_d )
    else
        return @sprintf( "ϵₙ-Greedy (c = %4.3f, d = %4.3f)", agent.param_c, agent.param_d )
    end
end


text:
# Autogenerated wrapper script for SLICOT_jll for x86_64-apple-darwin-libgfortran3
export libslicot

using libblastrampoline_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("SLICOT")
JLLWrappers.@declare_library_product(libslicot, "@rpath/libslicot.dylib")
function __init__()
    JLLWrappers.@generate_init_header(libblastrampoline_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libslicot,
        "lib/libslicot.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function get_methodinstance(@nospecialize(sig);
                            world=Base.get_world_counter(),
                            interp=Core.Compiler.NativeInterpreter(world))
    ms = Base._methods_by_ftype(sig, 1, Base.get_world_counter())
    @assert length(ms) == 1
    m = ms[1]
    display(m)
    mi = ccall(:jl_specializations_get_linfo,
               Ref{Core.MethodInstance}, (Any, Any, Any),
               m[3], m[1], m[2])
    return mi
end

function code_ircode_by_signature(@nospecialize(sig);
                                  world=Base.get_world_counter(),
                                  interp=Core.Compiler.NativeInterpreter(world))
    return [code_ircode(ccall(:jl_specializations_get_linfo,
                              Ref{Core.MethodInstance},
                              (Any, Any, Any),
                              data[3], data[1], data[2]);
                        world=world, interp=interp)
            for data in Base._methods_by_ftype(sig, -1, world)]
end

function code_ircode(@nospecialize(f), @nospecialize(types=Tuple);
                     world=Base.get_world_counter(),
                     interp=Core.Compiler.NativeInterpreter(world))
    return [code_ircode(mi; world=world, interp=interp)
            for mi in Base.method_instances(f, types, world)]
end

function code_ircode(mi::Core.Compiler.MethodInstance;
                     world=Base.get_world_counter(),
                     interp=Core.Compiler.NativeInterpreter(world))
    ccall(:jl_typeinf_begin, Cvoid, ())
    result = Core.Compiler.InferenceResult(mi)
    frame = Core.Compiler.InferenceState(result, false, interp)
    frame === nothing && return nothing
    if Core.Compiler.typeinf(interp, frame)
        opt_params = Core.Compiler.OptimizationParams(interp)
        opt = Core.Compiler.OptimizationState(frame, opt_params, interp)
        ir = Core.Compiler.run_passes(opt.src, opt.nargs - 1, opt)
        opt.src.inferred = true
    end
    ccall(:jl_typeinf_end, Cvoid, ())
    frame.inferred || return nothing
    resize!(ir.argtypes, opt.nargs)
    return ir => Core.Compiler.widenconst(result.result)
end

emit_optimized(f, tt...) =
    emit(typeof(f), tt,
         emit_fptr=false,
         dump_options=[Brutus.DumpCanonicalized])

emit_lowered(f, tt...) =
    emit(typeof(f), tt,
         emit_fptr=false, # TODO: change to true when ready
         dump_options=[Brutus.DumpLoweredToStd,
                       Brutus.DumpLoweredToLLVM])

emit_translated(f, tt...) =
    emit(typeof(f), tt,
         emit_fptr=false,
         dump_options=[Brutus.DumpTranslated])

emit_llvm(f, tt...) =
    emit(typeof(f), tt,
         emit_fptr=false, # TODO: change to true when ready
         dump_options=[Brutus.DumpLoweredToLLVM,
                       Brutus.DumpTranslateToLLVM])

function lit(emit_sym)
    @assert length(ARGS) == 1
    temp = Module(:Main)
    expr = quote
            include(fname) = Core.include(@__MODULE__, fname)
            import Brutus
            emit(args...;kwargs...) = Brutus.$(emit_sym)(args...;kwargs...)
            include($(ARGS[1]))
    end
    Base.eval(temp, Expr(:toplevel, expr))
end

max_stars_repo_name: stephenneuendorffer/brutus

max_issues_repo_name: stephenneuendorffer/brutus

max_forks_repo_name: stephenneuendorffer/brutus


text:
struct IsingOptions
    n_iters::Int
    n_save_iters::Int
    n_save::Int
end

"""
    `IsingOptions(;n_iters = 10^4, n_save_iters=1)`

Set options for samplers.

### Fields

* `n_iters`       - Set the number of iterations of the sampler
* `n_save_iters`  - Set the frequency at which iterations are saved.  If
                  n_save_iters=1, every iteration is saved.  If n_save_iters=n_iters,
                  only the final iteration is saved.
"""
function IsingOptions(;n_iters = 10^4, n_save_iters=1)
        return IsingOptions(n_iters, n_save_iters, n_iters÷n_save_iters)
end

"""
`RandomLattice` - Generate a random square lattice

### Fields
* `N` - Lattice size, N×N
"""
function RandomLattice(N)
    return Int8.(2*rand(Bool, N, N).-1)
end

"""
`RandomMLattice` - Generate a random square lattice with magnetization M.  This
assumes that M is compatible with the allowable magenetizations of the lattice

### Fields
* `N` - Lattice size, N×N
* `M` - Target magnetization, M
"""
function RandomMLattice(N,M)

    if M ≥ 0
        σ = Int8(1);
    else
        σ = Int8(-1);
    end
    x = σ  * ones(Int8, N, N);

    # compute number of sites to flip
    n_flip = (N^2 -abs(M)) ÷ 2;
    # sample the sites
    flip_sites = sample(collect(Iterators.product(1:N,1:N))[:],n_flip,replace=false)
    # flip
    for site in flip_sites
        x[site[1],site[2]] *=Int8(-1);
    end

    return x

end

"""
`GetNeighbors!` - Generate a list of neighbors

### Fields
* `𝒩` - Neighbor set to be populated.  Must be a mutable type
* `i` - i-index
* `j` - j-index
* `N` - Lattice size
"""
function GetNeighbors!(𝒩, i, j, N)
    𝒩[1][1] = mod1(i+1,N);
    𝒩[1][2] = j;

    𝒩[2][1] = mod1(i-1,N);
    𝒩[2][2] = j;

    𝒩[3][1] = i;
    𝒩[3][2] = mod1(j+1,N);

    𝒩[4][1] = i;
    𝒩[4][2] = mod1(j-1,N);

    𝒩
end


text:
export build_groups, build_groupgm, parse_config, build_network

using Clustering
using GZip

"Build GroupGM groups from a correlation matrix and a set of ids."
function build_groups(C, header)

    # force C to be perfectly symmetric (rounding errors make hclust unhappy)
    symC = deepcopy(C)
    for i in 1:size(C)[1], j in 1:size(C)[2]
        symC[j,i] = symC[i,j]
    end

    h = hclust(1 .- abs(symC), :complete)

    # build a list of all internal nodes in the hclust tree
    internalNodes = ASCIIString[]
    internalScores = Float64[]
    for i in 1:length(header)-1
        v1 = h.merge[i,1] < 0 ? header[-h.merge[i,1]] : internalNodes[h.merge[i,1]]
        v2 = h.merge[i,2] < 0 ? header[-h.merge[i,2]] : internalNodes[h.merge[i,2]]
        push!(internalNodes, "$v1 $v2") # note this keeps sub group names next to each other in larger group names
    end

    # build the results
    groups = Tuple{Float64, ASCIIString}[]
    for (i,v) in enumerate(header)
        push!(groups, (0.0, v))
    end
    for (i,v) in enumerate(internalNodes)
        push!(groups, (h.height[i], v))
    end

    groups
end

"Build a GroupGM score matrix form an inverse correlation matrix and the given groups."
function build_groupgm(I, header, groups; combiner=sum)

    # create a dense array and a map to its indexes
    G = zeros(length(groups), length(groups))
    indexMap = Dict{ASCIIString,Int64}()
    for i in 1:length(header)
        indexMap[header[i]] = i
    end

    # fill in our group matrix
    indexGroups = collect(map(g->collect(map(x->indexMap[x], split(g[2]))), groups))
    for i in 1:length(indexGroups)
        for j in i:length(indexGroups)
            G[j,i] = G[i,j] = combiner(I[indexGroups[j],indexGroups[i]])
        end
    end

    G,[g[2] for g in groups]
end

# each line is tab separated with the follow format (tab separated)
# BED_FILE_NAME SHORT_TITLE LONG_TITLE CELL_TYPE LAB EXPERIMENT_ID ANTIBODY_ID TREATMENTS ORGANISM LIFE_STAGE LINK
function parse_config(stream, fileRoot)
    metadata = Dict()
    for line in eachline(stream)
        parts = split(strip(line), '\t')
        l = length(parts)

        obj = Dict()
        if ismatch(r".*\.bam$", parts[1])
            obj["fileType"] = "bam"
        elseif ismatch(r".*\.bed$", parts[1])
            obj["fileType"] = "bed"
        elseif ismatch(r".*\.bed\.gz$", parts[1])
            obj["fileType"] = "bed.gz"
        elseif ismatch(r".*\.narrowPeak$", parts[1])
            obj["fileType"] = "narrowPeak"
        elseif ismatch(r".*\.narrowPeak\.gz$", parts[1])
            obj["fileType"] = "narrowPeak.gz"
        else
            obj["fileType"] = "unknown"
        end
        obj["dataFile"] = joinpath(fileRoot, parts[1])
        obj["name"] = l >= 2 ? parts[2] : basename(obj["dataFile"])
        obj["description"] = l >= 3 ? parts[3] : obj["name"]
        obj["cellType"] = l >= 4 ? parts[4] : "Unlabeled custom upload"
        obj["lab"] = l >= 5 ? parts[5] : "Unlabeled custom upload"
        obj["id"] = l >= 6 ? parts[6] : obj["name"]*"_"*randstring(15)
        obj["antibody"] = l >= 7 ? parts[7] : "Unlabeled custom upload"
        obj["treatments"] = l >= 8 ? parts[8] : "None"
        obj["organism"] = l >= 9 ? parts[9] : "Homo sapiens"
        obj["lifeStage"] = l >= 10 ? parts[10] : "None"
        obj["link"] = l >= 11 ? parts[11] : ""
        metadata[obj["id"]] = obj
    end
    metadata
end

"Build a network from a GroupGM matrix, groups, and metadata."
function build_network(G, groups, metadata; threshold=0.2, groupScoreThreshold=0.7, quiet=false)

    function find_parent(groups, id, ind)
        for i in ind:length(groups)
            if searchindex(groups[i][2], id) != 0
                return i
            end
        end
        0
    end

    # build the node data
    nodes = Any[]
    for (i,id) in enumerate(map(x->x[2], groups))
        if search(id, ' ') == 0
            md = metadata[id]
            push!(nodes, Dict{Any,Any}(
                "name" => get(md, "name", id),
                "group" => 0,
                "groupScore" => 1.0,
                "parent" => find_parent(groups, id, i+1)-1,
                "data" => Dict{Any,Any}(
                    "id" => id,
                    "description" => get(md, "description", get(md, "name", id)),
                    "cellType" => get(md, "cellType", "Unknown"),
                    "lab" => get(md, "lab", "Unknown"),
                    "organism" => get(md, "organism", "Unknown"),
                    "lifeStage" => get(md, "lifeStage", "Unknown"),
                    "treatments" => get(md, "treatments", "None"),
                    "antibody" => get(md, "antibody", "Unknown")
                )
            ))
        else
            push!(nodes, Dict{Any,Any}(
                "name" => "$i",
                "group" => 1,
                "groupScore" => 1-groups[i][1],
                "parent" => find_parent(groups, id, i+1)-1,
                "data" => Dict{Any,Any}(
                    "id" => "$i"
                )
            ))
        end
    end

    allNodes = Array{ASCIIString,1}[]
    for i in 1:length(groups)
        push!(allNodes, ASCIIString[])
    end
    for (i,(v,id)) in enumerate(groups)

        # if we are a leaf node then add our id to all ancestors (including ourselves)
        if !(' ' in id)
            ind = i
            while ind != 0
                push!(allNodes[ind], id)
                ind = nodes[ind]["parent"]+1
            end
        end
    end

    # build the edge array
    links = Any[]
    for i in 1:length(nodes)
        idi = nodes[i]["data"]["id"]
        for j in i+1:length(nodes)
            idj = nodes[j]["data"]["id"]

            # make sure we pass the threshold and are not between nested groups
            if (abs(G[i,j]) > threshold && searchindex(groups[i][2], groups[j][2]) == 0
                && searchindex(groups[j][2], groups[i][2]) == 0
                && groups[i][1] < 1-groupScoreThreshold && groups[j][1] < 1-groupScoreThreshold)

                d = Dict{Any,Any}(
                    "source" => i-1,
                    "target" => j-1,
                    "coeff" => -G[i,j],
                    "labels" => get(metadata, "$idi,$idj", ASCIIString[])
                )

                # find all the labels that are represented by the edges between these groups
                labels = ASCIIString[]
                for id1 in allNodes[i]
                    for id2 in allNodes[j]
                        if haskey(metadata, "$id1,$id2")
                            push!(labels, metadata["$id1,$id2"]...)
                        end
                    end
                end

                if length(labels) > 0
                    d["labels"] = unique(labels)
                end

                push!(links, d)
            end
        end
    end

    if !quiet
        println(STDERR, "Created network with $(length(links)) links and $(length(nodes)) nodes.")
    end
    Dict("nodes" => nodes, "links" => links)
end


text:
module Dates

import ..Clarus

function fixingdates(;params...)
  return Clarus.api_request("Dates","FixingDates";params...)
end

function fxforwarddate(;params...)
  return Clarus.api_request("Dates","FXForwardDate";params...)
end

function fxoptiondate(;params...)
  return Clarus.api_request("Dates","FXOptionDate";params...)
end

function irdspotdate(;params...)
  return Clarus.api_request("Dates","IRDSpotDate";params...)
end

function maintenanceperiods(;params...)
  return Clarus.api_request("Dates","MaintenancePeriods";params...)
end

function schedule(;params...)
  return Clarus.api_request("Dates","Schedule";params...)
end

end #MODULE-END


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Prjtrellis_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Prjtrellis")
JLLWrappers.@generate_main_file("Prjtrellis", UUID("6a9aa77d-afea-5b9d-b858-ecc31eb61f73"))
end  # module Prjtrellis_jll


text:
const mathml_ns = "http://www.w3.org/1998/Math/MathML"

"""
    mathml_to_nums()

given a filename or an `EzXML.Document` or `EzXML.Node`,
finds all of the <ci>s and defines them as Symbolics.jl Nums
returns a Vector{Num}.
Note, the root namespace needs to be MathML
"""
function mathml_to_nums end

function mathml_to_nums(fn::AbstractString)
    doc = readxml(fn)
    mathml_to_nums(doc)
end

function mathml_to_nums(xml::EzXML.Document)
    doc_root = EzXML.root(xml)
    mathml_to_nums(doc_root)
end

"""
    extract_mathml()

given a filename, `EzXML.Document`, or `EzXML.Node`
returns all of the MathML nodes.
"""
function extract_mathml end

function extract_mathml(fn::AbstractString)
    extract_mathml(readxml(fn))
end

function extract_mathml(doc::EzXML.Document)
    extract_mathml(EzXML.root(doc))
end

function extract_mathml(node::EzXML.Node)
    disambiguate_equality!(node)
    findall("//x:math", node, ["x" => mathml_ns])
end

"""
    @disambiguate_equality!

utility function to replace <eq> inside piecewise subtrees to
disambiguate from the assignement <eq>
"""
function disambiguate_equality!(node)
    nodes = findall("//x:piecewise//x:eq", node, ["x" => mathml_ns])
    for n in nodes
        setnodename!(n, "equal")
    end
end

"""
    @xml_str(s)

utility macro for parsing xml strings into node
"""
macro xml_str(s)
    parsexml(s).root
end

"""
    @MathML_str(s)

utility macro for parsing xml strings into symbolics
"""
macro MathML_str(s)
    MathML.parse_str(s)
end

# bindings for viewing call tree
AbstractTrees.children(n::EzXML.Node) = elements(n)
AbstractTrees.printnode(io::IO, node::EzXML.Node) = print(io, getproperty(node, :name))
AbstractTrees.nodetype(::EzXML.Node) = EzXML.Node

function custom_root(x)
    length(x) == 1 ? sqrt(x...) : Base.:^(x[2], x[1])
end

"ensure theres only one independent variable, returns false if more than one iv"
function check_ivs(node)
    x = findall("//x:bvar", node, ["x" => MathML.mathml_ns])
    all(y -> y.content == x[1].content, x)
end

# conditional and rounding hacks
H(x) = IfElse.ifelse(x >= 0, one(x), zero(x))
const ϵ = eps(Float64)
frac(x) = 0.5 - atan(cot(π * x)) / π
heaviside_or(x) = length(x) == 1 ? x[1] : x[1] + heaviside_or(x[2:end]) - x[1] * heaviside_or(x[2:end])