text:
function gettargetDF(VAF; fmin = 0.05, fmax = 0.75)
  targetdataCDF = getCDF(VAF, 0.001, fmin = fmin, fmax = fmax)
  return targetdataCDF, VAF
end

function getmodel(abcres, model)
    indeces = map(x -> x.model, abcres.particles) .== model
    abcresnew = deepcopy(abcres)
    abcresnew.particles = abcresnew.particles[indeces]
    abcresnew.parameters = abcresnew.parameters[model]
    abcresnew.weights = abcresnew.weights[model]
    return abcresnew
end

function getCDF(VAF::Array, step_size::Float64; fmin = 0.05, fmax = 0.7)
  #fast way to calculate CDF
  out = cumsum(fit(Histogram, VAF, fmax:-step_size:fmin,closed=:left).weights[1:end - 1])
  out = out .- out[1]
  return out
end

function tumourABCneutral(parameters, constants, targetdata)
    #function to simulate neutral tumour growth, outputs the euclidean distance between targetdata and simulation and
    #and histogram values which can be used to plot the output

    cst = constants
    simdata = simulate(nclones = 0,
                ploidy = cst[1],
                read_depth = cst[2],
                μ = parameters[1],
                clonalmutations = round(parameters[2]),
                d = cst[3],
                b = cst[4],
                ρ = cst[5],
                Nmax = 500,
                timefunction = cst[7],
                detectionlimit = cst[8],
                s = Float64[],
                tevent = Float64[],
                cellularity = parameters[3])

    simdatacdf = getCDF(simdata.sampleddata.VAF, 0.001,
                    fmin = cst[9],
                    fmax = cst[10])

    return euclidean(simdatacdf, targetdata), [simdata.sampleddata.DF], true
end

function tumourABCselection(parameters, constants, targetdata)

    # we only consider the tumour to have a clone if it is observeable in the data, above frequency 5%.

    cst = constants
    simdata = simulate(nclones = 1,
                ploidy = cst[1],
                read_depth = cst[2],
                μ = parameters[1],
                clonalmutations = round(parameters[2]),
                d = cst[3],
                b = cst[4],
                ρ = cst[5],
                Nmax = cst[6],
                s = [parameters[3]],
                tevent = [parameters[4]],
                timefunction = cst[7],
                detectionlimit = cst[8],
                cellularity = parameters[5])

    simdatacdf = getCDF(simdata.sampleddata.VAF, 0.001,
                    fmin = cst[9],
                    fmax = cst[10])

    if length(simdata.output.subclonemutations) == 0
        out = [simdata.sampleddata.DF, 0.0, 0.0, 0.0]
    else
        out = [simdata.sampleddata.DF, simdata.output.subclonemutations[1], simdata.output.Ndivisions[1], simdata.output.clonefreq[1]]
    end

    #return true if clone is between 0.05 and 0.95
    c = false
    if ((sum(simdata.output.clonefreq.<0.95).==1) & (sum(simdata.output.clonefreq.>0.05).==1))[1] == true
        c = true
    end

    euclidean(simdatacdf, targetdata), out, c
end

function tumourABCselection2(parameters, constants, targetdata)

    cst = constants
    simdata = simulate(nclones = 2,
                ploidy = cst[1],
                read_depth = cst[2],
                μ = parameters[1],
                clonalmutations = round(parameters[2]),
                d = cst[3],
                b = cst[4],
                ρ = cst[5],
                Nmax = cst[6],
                s = [parameters[3], parameters[5]],
                tevent = [parameters[4], parameters[6]],
                timefunction = cst[7],
                detectionlimit = cst[8],
                cellularity = parameters[7])

    simdatacdf = getCDF(simdata.sampleddata.VAF, 0.001,
                    fmin = cst[9],
                    fmax = cst[10])

    if length(simdata.output.subclonemutations) < 2
        out = [simdata.sampleddata.DF, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    else
        out = [simdata.sampleddata.DF, simdata.output.subclonemutations[1],
        simdata.output.subclonemutations[2],
        simdata.output.Ndivisions[1],
        simdata.output.Ndivisions[2],
        simdata.output.clonefreq[1],
        simdata.output.clonefreq[2]]
    end

    #return true if clones are between 0.05 and 0.95
    c = false
    if ((sum(simdata.output.clonefreq.<0.95).==2) & (sum(simdata.output.clonefreq.>0.05).==2))[1] == true
        c = true
    end

    euclidean(simdatacdf, targetdata), out, c
end

timefunc() = 1
timefuncrand() = -log(rand())

function getsetup(maxclones; nparticles = 100, maxiterations = 10^4, convergence = 0.1, ϵT = 1.0, read_depth = 100.0, Nmax = 10^3, detectionlimit = 0.05, modelkern = 0.4, ϵ1 = 10^6, mincellularity = 0.1, maxcellularity = 1.1, ρ = 0.0, maxclonalmutations = 10000.0, maxmu = 620.0, timefunction = timefunc, ploidy = 2, d = 0.0, b = log(2), fmin = 0.05, fmax = 0.05, Nmaxinf = 10^10)

  #function to define priors, constants and create ABCSMC model type

  read_depth = read_depth
  Nmax = Nmax
  ρ = ρ

  cst = [ploidy, read_depth, d, b, ρ, Nmax, timefunction, detectionlimit, fmin, fmax];

  priormu = [0.01, maxmu]
  priorcm = [0.0, Float64(maxclonalmutations)]
  priorcellularity = [mincellularity, maxcellularity]

  #need to create Prior type which has a distribution type array with a corresponding distribution specific parameter array
  priorneutral = Prior([Uniform(priormu...),
     Uniform(priorcm...),
     Uniform(priorcellularity...)])

  priort = [2.0, log(Nmax)/log(2) + (MathConstants.eulergamma / log(2))]
  maxsel = selection(log(2), 0.99, priort[2], priort[2] - 1)
  priorsel = [0.0, maxsel]

priorselection = Prior([Uniform(priormu...),
   Uniform(priorcm...),
   Uniform(priorsel...),
   Uniform(priort...),
   Uniform(priorcellularity...)])

priorselection2 = Prior([Uniform(priormu...),
   Uniform(priorcm...),
   Uniform(priorsel...),
   Uniform(priort...),
   Uniform(priorsel...),
   Uniform(priort...),
   Uniform(priorcellularity...)])

  nparamatersneutral = 3
  nparametersselection = 5
  nparametersselection2 = 7

  #for the ABCSMC Model selection, each input is an array of model specific parameters,

  if maxclones == 1
    setup = ABCSMCModel(
    [tumourABCneutral, tumourABCselection], #simulation functions
    [nparamatersneutral, nparametersselection], # number of parameters for each model
    ϵT, #desired final tolerence
    [priorneutral, priorselection], #priors
    constants = [cst, cst], #constants
    nparticles = nparticles,
    maxiterations = maxiterations,
    convergence = convergence,
    modelkern = modelkern,
    ϵ1 = ϵ1,
    other = Nmaxinf);
  elseif maxclones == 2
    setup = ABCSMCModel(
    [tumourABCneutral, tumourABCselection, tumourABCselection2], #simulation functions
    [nparamatersneutral, nparametersselection, nparametersselection2], # number of parameters for each model
    ϵT, #desired final tolerence
    [priorneutral, priorselection, priorselection2], #priors
    constants = [cst, cst, cst], #constants
    nparticles = nparticles,
    maxiterations = maxiterations,
    convergence = convergence,
    modelkern = modelkern,
    ϵ1 = ϵ1,
    other = Nmaxinf);
  end

  return setup
end

"""
    fitABCmodels(data::Array{Float64, 1}, sname::String; <keyword arguments>)

Fit a stochastic model of cancer evolution to cancer sequencing data using Approximate Bayesian computation and infer the number of subclones (up to 2) and their relative fitness and time they emerge.
...
## Arguments
- `read_depth = 200.0`: Mean read depth of the target data set
- `minreads = 5`: Minimum number of reads needed to identify a mutation, this should be the same as was set in the sequencing data analysis
- `minvaf = 0.0`: Minimum VAF to identify a mutation, will override minreads if > 0.0
- `fmin = 0.01`: Minimum range of VAF to perform inference, if minvaf > 0.0 then fmin = minvaf
- `fmax = 0.75`: Maximum range of VAF to perform inference
- `maxiterations = 10^6`: Maximum number of iterations before ABC terminates
- `maxclones = 2`: Maximum number of clones, can be 1 or 2
- `nparticles = 500`: Number of particles (ie samples) in the ABC output
- `Nmax = 10^4`: Maximum population size used to fit data, increase if suspect that there is a late arising subclone that dominates the distribution
- `resultsdirectory = "output" `: Directory where posterior distributions will be saved to, this directory will be created if it does not already exist.
- `progress = true`: Show progress of ABC sampler with `ProgressMeter` package
- `verbose = true`: Print out summary at each ABC population step.
- `save = false `: Save output or not
- `ϵ1 = 10^6 `: Target ϵ for first ABC step
- `firstpass = false`: If set to true will run a limited first pass of the algorithm to determine a good starting ϵ1 if this unknown, otherwise ϵ1 is effectively infinity.
- `Nmaxinf = 10^10`: Scales selection coefficient value assuming the tumour size is Nmaxinf. Once value >10^9 has limited effect.
- `ρ = 0.0`: Overdispersion parameter for beta-binomial model of sequencing data. ρ = 0.0 means model is binomial sampling
- `adaptpriors = true`: If true priors on μ and clonalmutations are adapted based on the number of mutations in the data set which provide an upper an lower limit, this is an experimental feature that needs further validation, although initial tests suggest it performs well and does not skew inferences
- `timefunction = timefunc`: Function for KMC algorithm timestep. timefunc returns 1 meaning the timestep is the average of stochastic process. Alternatively timefuncrand can be specified which uses `-log(rand())` to increase the time step, so it is exponentially distributed rather than the mean of the exponential distribution. We use population doublings as our unit of time so this does not change the method and is an algorithmic choice.
- `ploidy = 2`: ploidy of the genome
- `d = 0.0`: Death rate of the the host population in the tumour, it is the compound parameter μ/β (where β = (b-d)/b) that can be inferred from the VAF distribution hence setting d = 0.0 means there is one fewer parameter in the ABC and hence increases efficiency.
- `b = log(2)`: Birth rate of the population. Default is set to `log(2)` so that tumour doubles with each unit increase in t in the absence of cell death.
- `mincellularity = 0.1`: If some prior knowledge on cellularity is known this can be modifed
- `maxcellularity = 1.1`: If some prior knowledge on cellularity is known this can be modifed. This is set to > 1.0 so that there is some flexibility in the inference.
- `convergence = 0.07`: Convergence for ABC. If new population ϵ is within convergence % of previous population then ABC stops.
- `savepopulations = false`: Save results from intermediate populations.
...
"""
function fitABCmodels(data::Array{Float64, 1}, sname::String;
  read_depth = 200.0, minreads = 5, minvaf = 0.0, fmin = 0.05,
  fmax = 0.75, maxiterations = 10^4, maxclones = 2,
  nparticles = 500, Nmax = 10^4, resultsdirectory::String = "output",
  progress = true, verbose = false, save = false,
  ϵ1 = 10^6, mincellularity = 0.1, maxcellularity = 1.1, firstpass = false,
  Nmaxinf = 10^10, ρ = 0.0,
  adaptpriors = true, timefunction = timefunc, ploidy = 2, d = 0.0, b = log(2),
  maxmu = 500, maxclonalmutations = 5000, convergence = 0.07,
  savepopulations = false)

  println("It is recommended that you use mobster for this type of analysis moving forward. mobster is an R package that provides similar functionality with orders of magnitude increases in speed and has many other features. SubClonalSelection.jl will remain here but I am unlikely to actively develop the package or be able to provide much ongoing support.")

  #make output directories
  if save
    makedirectories(joinpath(resultsdirectory, sname, "finalpopulation"))
  end
  #sequencing error is 1%, otherwise detection limit is controlled by minimum number of reads to call a variant
  detectionlimit = maximum([0.01, minreads/read_depth])

  #if a hard cutoff for minvaf is specified detectionlimit is changed also changes fmin in this case as its assumed all mutations < minvaf have been removed or shouldn't be used
  if minvaf > 0.0
    detectionlimit = minvaf
  end

  targetdata, VAF = gettargetDF(data, fmin = fmin, fmax = fmax)
  targetdataCDF = targetdata
  if save != false
    writedlm(joinpath(joinpath(resultsdirectory, sname, "finalpopulation"),
    "data", "$(sname).txt"), VAF)
  end

  dl = detectionlimit
  eps1 = ϵ1

  #adaptively set priors on μ and clonalmutations based on number of mutations
  # two extreme cases are mutations are all subclonal which defines a maximum mu,
  #or mutations are all clonal hich defines a maximum clonalmutational load
  if adaptpriors == true
    maxmu = 1.5 * (length(data) / ((1/detectionlimit) - 1))
    maxclonalmutations = length(data)
  end

  nparts = nparticles
  if (firstpass == true) && (length(VAF) < 3000) #I found this not to be useful when number of mutations is large
    println("################################################")
    println("Running first pass to get starting point for ABC")
    println("################################################")
    println("")
    ABCsetup = getsetup(1, detectionlimit = detectionlimit,
    read_depth = read_depth,
    nparticles = 100,
    maxiterations = 3 * nparts,
    modelkern = 0.5,
    Nmax = Nmax,
    mincellularity = mincellularity,
    maxcellularity = maxcellularity,
    maxclonalmutations = maxclonalmutations,
    maxmu = maxmu,
    timefunction = timefunction,
    ploidy = ploidy,
    d = d,
    b = b,
    fmin = fmin,
    fmax = fmax,
    Nmaxinf = Nmaxinf
    )
    abcres = runabcCancer(ABCsetup, targetdataCDF, verbose = verbose, progress = progress);
    eps1 = abcres.ϵ[maximum([1, length(abcres.ϵ) - 1])]
    println("################################################")
    println("Now running inference with ϵ1 = $(eps1)")
    println("################################################")
    println("")
  end

  ABCsetup = getsetup(maxclones, detectionlimit = dl,
  read_depth = read_depth,
  maxiterations = maxiterations,
  nparticles = nparticles,
  modelkern = 0.5,
  convergence = convergence,
  ϵ1 = eps1,
  Nmax = Nmax,
  mincellularity = mincellularity,
  ρ = ρ,
  maxclonalmutations = maxclonalmutations,
  maxmu = maxmu,
  timefunction = timefunction,
  ploidy = ploidy,
  d = d,
  b = b,
  fmin = fmin,
  fmax = fmax,
  Nmaxinf = Nmaxinf
  )

  abcres = runabcCancer(ABCsetup, targetdataCDF, verbose = verbose,
  progress = progress, Nmaxinf = Nmaxinf,
  savepopulations = savepopulations, sname = sname, resultsdirectory = resultsdirectory, VAF = VAF);

  posteriors, DFmp = getresults(abcres, joinpath(resultsdirectory), sname, VAF, save = save, Nmaxinf = Nmaxinf, savepopulations = false);

  finalresults = Results(ABCsetup, abcres, VAF, posteriors, DFmp, sname);

  return finalresults
end

"""
    fitABCmodels(data::String, sname::String; <keyword arguments>)

If data is a string will read in file. File should be a 1 column text file with VAF values in the rows.
"""
function fitABCmodels(data::String, sname::String;
  read_depth = 200.0, minreads = 5, fmin = 0.05, minvaf = 0.0,
  fmax = 0.75, maxiterations = 10^4, maxclones = 2,
  nparticles = 500, Nmax = 10^4, resultsdirectory::String = "output",
  progress = true, verbose = false, save = false,
  ϵ1 = 10^6, mincellularity = 0.1, maxcellularity = 1.1,
  firstpass = false, Nmaxinf = 10^10, ρ = 0.0,
  adaptpriors = true, timefunction = timefunc, ploidy = 2, d = 0.0, b = log(2),
  convergence = 0.07,
  savepopulations = false)

  VAF = readdlm(data)[:, 1]

  return fitABCmodels(VAF, sname;
  fmin = fmin, fmax = fmax, minreads = minreads, minvaf = minvaf,
  read_depth = read_depth, maxiterations = maxiterations,
  maxclones = maxclones, nparticles = nparticles,
  Nmax = Nmax, resultsdirectory = resultsdirectory,
  progress = progress, verbose = verbose, save = save,
  ϵ1 = ϵ1, mincellularity = mincellularity, maxcellularity = maxcellularity,
  firstpass = firstpass, Nmaxinf = Nmaxinf,
  ρ = ρ, adaptpriors = adaptpriors,
  timefunction = timefunction, ploidy = ploidy, d = d, b = b,
  convergence = convergence,
  savepopulations = savepopulations)
end


text:
time_modeling(model::Model, srcGeometry, srcData, recGeometry, recData, perturbation, srcnum::UnitRange{Int64}, op::Char, mode::Int64) =
    time_modeling(model, srcGeometry, srcData, recGeometry, recData, perturbation, srcnum, op, mode, Options())

function time_modeling(model::Model, srcGeometry, srcData, recGeometry, recData, perturbation, srcnum::UnitRange{Int64}, op::Char, mode::Int64, options)
# time_modeling function for multiple sources. Depending on the operator and mode, this function distributes the sources
# and if applicable the input data amongst the available workers.

    p = default_worker_pool()
    time_modeling_par = remote(TimeModeling.time_modeling)
    time_modeling = retry(time_modeling_par)

    numSources = length(srcnum)
    results = Array{Any}(undef, numSources)

    # Process shots from source channel asynchronously
    @sync begin
        for j=1:numSources

            # local geometry for current position
            if srcGeometry == nothing
                srcGeometryLocal = nothing
            else
                srcGeometryLocal = subsample(srcGeometry,j)
            end
            if recGeometry == nothing
                recGeometryLocal = nothing
            else
                recGeometryLocal = subsample(recGeometry,j)
            end
            opt_local = subsample(options,j)
            numSources > 1 && (opt_local.save_wavefield_to_disk=true)    # don't collect wavefields on master

            # Parallelization
            if op=='F' && mode==1
                srcDataLocal = Array{Any}(undef, 1)
                srcDataLocal[1] = srcData[j]
                @async results[j] = time_modeling(model, srcGeometryLocal, srcDataLocal, recGeometryLocal, nothing, nothing, j, op, mode, opt_local)
            elseif op=='F' && mode==-1
                recDataLocal = Array{Any}(undef, 1)
                recDataLocal[1] = recData[j]
                @async results[j] = time_modeling(model, srcGeometryLocal, nothing, recGeometryLocal, recDataLocal, nothing, j, op, mode, opt_local)
            elseif op=='J' && mode==1
                srcDataLocal = Array{Any}(undef, 1)
                srcDataLocal[1] = srcData[j]
                @async results[j] = time_modeling(model, srcGeometryLocal, srcDataLocal, recGeometryLocal, nothing, perturbation, j, op, mode, opt_local)
            elseif op=='J' && mode==-1
                srcDataLocal = Array{Any}(undef, 1)
                srcDataLocal[1] = srcData[j]
                recDataLocal = Array{Any}(undef, 1)
                recDataLocal[1] = recData[j]
                @async results[j] = time_modeling(model, srcGeometryLocal, srcDataLocal, recGeometryLocal, recDataLocal, nothing, j, op, mode, opt_local)
            end
        end
    end

    if op=='F' || (op=='J' && mode==1)
        argout1 = results[1]
        for j=2:numSources
            argout1 = [argout1; results[j]]
        end
    elseif op=='J' && mode==-1
        argout1 = results[1]
        for j=2:numSources
            argout1 += results[j]
        end
    else
        error("operation no defined")
    end
    return argout1
end

ext: jl


text:
module TodoMVC

using Redux

include("Todo.jl")
using .Todo

include("TodoFilter.jl")
using .TodoFilter

# state
struct State <: AbstractImmutableState
    todos::Vector{Todo.State}
    visibility::TodoFilter.State
end

const INITIAL_STATE = State(Todo.INITIAL_STATE, TodoFilter.State(TodoFilter.SHOW_ALL))

# reducers
todo_mvc(state::AbstractState, action::AbstractAction) = state
todo_mvc(state::Vector{<:AbstractState}, action::AbstractAction) = state
function todo_mvc(state::State, action::AbstractSyncAction)
    next_todos = Todo.todo(state.todos, action)
    next_visibility = TodoFilter.todo_filter(state.visibility, action)
    return State(next_todos, next_visibility)
end

end # module


text:
@inline function first_cache_buffer(::Val{T}) where {T}
  first_cache_buffer(Val{T}(), first_cache_size(Val(T)))
end
@static if Sys.WORD_SIZE == 32
  @inline function first_cache_buffer(::Val{T}, N) where {T}
    reinterpret(Ptr{T}, ACACHEPTR[] + ((Threads.threadid()-1) * N) * static_sizeof(T))
  end
else
  @inline function first_cache_buffer(::Val{T}, ::StaticInt{N}) where {T,N}
    MemoryBuffer{N,T}(undef)
  end
end

BCache(i::Integer) = BCache(BCACHEPTR[]+cld_fast(second_cache_size()*i, Threads.nthreads()), i % UInt)
BCache(::Nothing) = BCache(BCACHEPTR[], nothing)

@inline Base.pointer(b::BCache) = b.p
@inline Base.unsafe_convert(::Type{Ptr{T}}, b::BCache) where {T} = Base.unsafe_convert(Ptr{T}, b.p)

function _use_bcache()
  while Threads.atomic_cas!(BCACHE_LOCK, zero(UInt), typemax(UInt)) != zero(UInt)
    pause()
  end
  return BCache(nothing)
end
@inline _free_bcache!(::BCache{Nothing}) = reset_bcache_lock!()

_use_bcache(::Nothing) = _use_bcache()
function _use_bcache(i)
  f = one(UInt) << i
  while (Threads.atomic_or!(BCACHE_LOCK, f) & f) != zero(UInt)
    pause()
  end
  BCache(i)
end
_free_bcache!(b::BCache{UInt}) = (Threads.atomic_xor!(BCACHE_LOCK, one(UInt) << b.i); nothing)

"""
  reset_bcache_lock!()

Currently not using try/finally in matmul routine, despite locking.
So if it errors for some reason, you may need to manually call `reset_bcache_lock!()`.
"""
@inline reset_bcache_lock!() = (BCACHE_LOCK[] = zero(UInt); nothing)


text:
### This is a validation test
### the objective is to show that KPLS can easily fit any curve using the same matrix
### If not, the regressor could be wrong

### Example extracted from a python KPLS implementation: https://github.com/jhumphry/regressions/blob/master/examples/kpls_example.py

using PLSRegressor
using Gadfly
import Random

Random.seed!(1)

z(x)     = 4.26 * (exp.(-x) - 4 * exp.(-2.0*x) + 3 * exp.(-3.0*x))
x_values = range(0.0,step=3.5,length=100)
z_pure   = z(x_values)
noise    = randn(100)
z_noisy  = z_pure + noise
X        = collect(x_values)[:,:]
Y        = z_noisy[:,:] #z_pure

global min_mae = 10
global best_pred
global best_w = 10
global best_g = 10

for g in [1,2], w in range(0.01,step=3,length=10)
    print(".")
    model      = PLSRegressor.fit(X,Y,centralize=true,nfactors=g,kernel="rbf",width=w)
    Y_pred     = PLSRegressor.predict(model,X)
    mae = mean(abs.(Y .- Y_pred))
    if mae < min_mae
       min_mae = mae
       best_pred = Y_pred[:]
       best_g    = g
       best_w    = w
   end
end

print("[KPLS] min mae error : $(min_mae)")
print("[KPLS] best factor : $(best_g)")
print("[KPLS] best width : $(best_w)")

plot([Y,best_pred ], x=Row.index, y=Col.value, color=Col.index, Geom.line)


text:
"""
## bmat_nonaxi!

This subroutine forms the strain-displacement matrix for
axisymmetric solids subjected to non-axisymmetric loading.

### Method
```julia
bmat_nonaxi!(bee, radius, coord, deriv, fun, iflag, lth)
```

### Arguments
```julia
* bee::Matrix{Float64}         : Bee matrix (Updated)
* radius::Float64              : r coordinate of the Gauss point
* coord::Matrix{Float64}       : Nodal coordinate matrix
* deriv::Matrix{Float64}       : Derivative
* fun::Vector{Float64}         : Shape function
* iflag::Int                 : 1 = symmetric, -1 = anti-symmetric
* lth::Int                   : Harmonic
```
"""
function bmat_nonaxi!(bee::Matrix{Float64}, radius::Float64,
  coord::LinearAlgebra.Adjoint{Float64,Array{Float64,2}}, deriv::Matrix{Float64},
  fun::Vector{Float64}, iflag::Int, lth::Int)
  #
  # This subroutine forms the strain-displacement matrix for
  # axisymmetric solids subjected to non-axisymmetric loading.
  #
  # Updates bee, radius
  #
 nod = size(deriv, 2)
 radius = sum(fun.*coord[:,1])
 for m in 1:nod
   n = 3 * m
   k = n-1
   l = k-1
   bee[1,l] = deriv[1, m]
   bee[2,k] = deriv[2, m]
   bee[3,l] = fun[m]/radius
   bee[3,n] = iflag*lth*bee[3, l]
   bee[4,l] = deriv[2, m]
   bee[4,k] = deriv[1, m]
   bee[5,k] = -iflag*lth*fun[m]/radius
   bee[5,n] = deriv[2, m]
   bee[6,l] = bee[5,k]
   bee[6,n] = deriv[1, m]-fun[m]/radius
 end
 (radius, bee)
end


text:
### A Pluto.jl notebook ###
# v0.14.1

using Markdown
using InteractiveUtils

# ╔═╡ b37139ed-ca54-4a70-a4aa-190fc4ff223a
using PlutoUI

# ╔═╡ 48988162-68d1-11eb-22dc-2bbc0e41043b
using Test

# ╔═╡ cb0ae3b2-42b8-485b-9c04-c1190e5eed94
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ b6e6f16a-68ce-11eb-320d-6b292123f85a
md"""
### Functions in julia
"""

# ╔═╡ d8456146-68ce-11eb-2872-3d66933b3385
begin
	DT = Float64

	function δ(a::DT, b::DT, c::DT)::DT
    	b * b - 4 * a * c
	end
end

# ╔═╡ fbec986e-68ce-11eb-1598-619db8bca4cc
function root(a::DT, b::DT, c::DT)::NamedTuple{(:x1, :x2),Tuple{DT, DT}}
    Δ = δ(a, b, c)
    Δ < 0 && throw(ArgumentError("No solution in ℜ"))

	x₁ = (-b + √ Δ) / 2a
   	x₂ = (-b - √ Δ) / 2a
    (x1 = x₁, x2 = x₂) # named tuple
end

# ╔═╡ 3750513e-68cf-11eb-0e76-7169551337af
with_terminal() do
	println(typeof(δ(1., 1., 1.)))
end

# ╔═╡ 50a85e6a-68cf-11eb-0e5e-2f2cbbfcce73
with_terminal() do
	r = root(2., -3., 1.)
	println("r: $(r) / typeof(r): $(typeof(r))")
end

# ╔═╡ 80512f66-68cf-11eb-3f58-498f113425d4
with_terminal() do
	r = root(1., -2., 1.)
	println("r: $(r) / typeof(r): $(typeof(r))")
end

# ╔═╡ 94a3ce9c-68cf-11eb-1a26-e57d02e78f19
with_terminal() do
	try
    	println(root(1., -1., 1.))
	catch ex
    	println("Intercepted: ", string(ex))
	end
end

# ╔═╡ a5f7141a-68cf-11eb-3dd1-9b04e8b62317
md"""
#### Function composition
"""

# ╔═╡ b2f4f754-68cf-11eb-3d27-dd9df6feb65f
function f(x::Real)::Real  ## Real is an abstract type
    2. * x
end

# ╔═╡ c06c4d8a-68cf-11eb-133e-ed1e6ae4c889
function g(x::Real)::Real
    x - 5.
end

# ╔═╡ cadfceca-68cf-11eb-1b1b-7fcc66d4c632
(g ∘ f)(4.)  # type: \circ<tab>: g(f(4.)) = 2 * 4. - 5, = 3

# ╔═╡ d97425bc-68cf-11eb-20c0-2dcaf9ede5f4
(f ∘ g)(7.)

# ╔═╡ e0553894-68cf-11eb-0e1f-5be65d986ed8
ary = [x for x in 0:5:100]

# ╔═╡ e4d32200-68cf-11eb-2ff7-0b227d2de510
nary = map(x -> 3x + 2, ary) |>
    x -> filter(isodd, x)

# ╔═╡ f250750e-68cf-11eb-3e4e-b5fcbcdf2320
map(x -> 3x + 2, ary) |>
    x -> filter(isodd, x) |>
    x -> reduce(+, x, init=0)

# ╔═╡ fcf27c16-68cf-11eb-3ae6-9b5a3d98c6ab
reduce(+, [x for x in 1:10], init=0)

# ╔═╡ 07a84580-68d0-11eb-0dcc-2dfc25b93d8b
reduce(+, 1:10, init=11)

# ╔═╡ 0b7b3b18-68d0-11eb-0e84-57cbc7e344d1
md"""
#### Function and Closure
"""

# ╔═╡ 1c25a7be-68d0-11eb-29f3-4dd237a2ef2e
adder = function(α::Number=2)
    x -> x + α
end

# ╔═╡ 44e600ba-68d0-11eb-122f-e9973916ffc8
add2 = adder()   # ≡ ad2(x) = x + 2

# ╔═╡ 53b8de8a-68d0-11eb-34e0-c1bbb3824ede
add2(2)

# ╔═╡ 7174059e-68d0-11eb-01fe-e712691890a4
begin
	add3 = adder(3)
	add3(10)
end

# ╔═╡ 840923b0-68d0-11eb-039f-05cd7abb2c04
map(add2, 1:10)

# ╔═╡ 8c579696-68d0-11eb-1321-ad2e8290b3bd
counter = function(; init=0)
    cnt = init

    inc() = cnt += 1
	dec() = cnt -= 1
	reset() = cnt = 0

	## Declare public:
    ()->(inc, dec, reset)
end

# ╔═╡ 2542f116-68d1-11eb-0eb2-0fa4398666ee
begin
	xc = counter()

	@test xc.inc() == 1
	@test xc.inc() == 2
	@test xc.inc() == 3

	@test xc.dec() == 2

	@test xc.reset() == 0
end

# ╔═╡ 6d948d44-68d1-11eb-283f-95c7209184f3
begin
	yc = counter(init=100)

	@test yc.inc() == 101
	@test yc.inc() == 102
	@test yc.inc() == 103

	@test yc.dec() == 102

	@test yc.reset() == 0
end

# ╔═╡ 8d553214-68d1-11eb-29d7-d791c9dff9f4
counter_ng = function(;init=0)
    cnt = init

	# closure
    λ = function(action=:inc)
        if action == :inc
            cnt += 1
        elseif action == :dec
            cnt -= 1
        elseif action == :reset
            cnt = 0
        else
            throw(ArgumentError("action $(action) not yet supported..."))
        end
    end

    λ
end

# ╔═╡ b2a41b98-68d1-11eb-21b3-9b48c657591d
begin
	ix = counter_ng()

	@test ix(:inc) == 1
	@test ix(:inc) == 2
	@test ix(:inc) == 3

	@test ix(:dec) == 2
	@test ix(:dec) == 1

	@test ix(:reset) == 0
end

# ╔═╡ eadfe19a-68d1-11eb-05c0-43915ae2d1b1
begin
	jy = counter_ng(init=100)

	@test jy(:inc) == 101
	@test jy(:inc) == 102
	@test jy(:inc) == 103

	@test jy(:dec) == 102
	@test_throws ArgumentError jy(:rec)

	@test jy(:reset) == 0
end

# ╔═╡ 3dc862b0-a58a-4731-a3c8-4eb44ded1443
html"""
<style>
  main {
        max-width: calc(800px + 25px + 6px);
  }
  .plutoui-toc.aside {
    background: linen;
  }
  h3, h4 {
        background: wheat;
  }
</style>
"""

# ╔═╡ Cell order:
# ╠═b37139ed-ca54-4a70-a4aa-190fc4ff223a
# ╠═cb0ae3b2-42b8-485b-9c04-c1190e5eed94
# ╟─b6e6f16a-68ce-11eb-320d-6b292123f85a
# ╠═d8456146-68ce-11eb-2872-3d66933b3385
# ╠═fbec986e-68ce-11eb-1598-619db8bca4cc
# ╠═48988162-68d1-11eb-22dc-2bbc0e41043b
# ╠═3750513e-68cf-11eb-0e76-7169551337af
# ╠═50a85e6a-68cf-11eb-0e5e-2f2cbbfcce73
# ╠═80512f66-68cf-11eb-3f58-498f113425d4
# ╠═94a3ce9c-68cf-11eb-1a26-e57d02e78f19
# ╟─a5f7141a-68cf-11eb-3dd1-9b04e8b62317
# ╠═b2f4f754-68cf-11eb-3d27-dd9df6feb65f
# ╠═c06c4d8a-68cf-11eb-133e-ed1e6ae4c889
# ╠═cadfceca-68cf-11eb-1b1b-7fcc66d4c632
# ╠═d97425bc-68cf-11eb-20c0-2dcaf9ede5f4
# ╠═e0553894-68cf-11eb-0e1f-5be65d986ed8
# ╠═e4d32200-68cf-11eb-2ff7-0b227d2de510
# ╠═f250750e-68cf-11eb-3e4e-b5fcbcdf2320
# ╠═fcf27c16-68cf-11eb-3ae6-9b5a3d98c6ab
# ╠═07a84580-68d0-11eb-0dcc-2dfc25b93d8b
# ╟─0b7b3b18-68d0-11eb-0e84-57cbc7e344d1
# ╠═1c25a7be-68d0-11eb-29f3-4dd237a2ef2e
# ╠═44e600ba-68d0-11eb-122f-e9973916ffc8
# ╠═53b8de8a-68d0-11eb-34e0-c1bbb3824ede
# ╠═7174059e-68d0-11eb-01fe-e712691890a4
# ╠═840923b0-68d0-11eb-039f-05cd7abb2c04
# ╠═8c579696-68d0-11eb-1321-ad2e8290b3bd
# ╠═2542f116-68d1-11eb-0eb2-0fa4398666ee
# ╠═6d948d44-68d1-11eb-283f-95c7209184f3
# ╠═8d553214-68d1-11eb-29d7-d791c9dff9f4
# ╠═b2a41b98-68d1-11eb-21b3-9b48c657591d
# ╠═eadfe19a-68d1-11eb-05c0-43915ae2d1b1
# ╟─3dc862b0-a58a-4731-a3c8-4eb44ded1443

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
# include packages and other codes -
include("Include.jl")

# make an executable script -
function main()

    # ask the user for thier net id -
    print("What's your netid? ")

    # Calling rdeadline() function
    netid = readline()
    uuid_value = uuid4()

    # create a message, then display it -
    message = "Hello World from $(netid)! Your secret (random) 5440/7770 code is: $(uuid_value)"
    println(message)
end

# execute -
main()


text:
using Test

using PSFs
using NDTools

pp = PSFParams()
sz = (128,128,128)
sampling=(0.100,0.100,0.150)

@testset "Compare scalar aPSFs" begin
    a_prop = apsf(PSFs.MethodPropagate, sz,pp, sampling=sampling);
    a_sincR = apsf(PSFs.MethodSincR, sz,pp, sampling=sampling);
    a_shell = apsf(PSFs.MethodShell, sz,pp, sampling=sampling);
    a_iter = apsf(PSFs.MethodPropagateIterative, sz,pp, sampling=sampling);
    sz_big = (512,512,128)
    a_prop2 = NDTools.select_region(apsf(PSFs.MethodPropagate, sz_big,pp, sampling=sampling), new_size=sz);
    # @vt a_prop2 a_iter a_sincR  a_prop a_shell
    # @vtp ft2d(a_prop2) ft2d(a_iter) ft2d(a_sincR) ft2d(a_prop) ft2d(a_shell)
end

pp = PSFParams(pol=pol_x)
sz = (128,128,128)
sampling=(0.100,0.100,0.150)

@testset "Compare scalar aPSFs" begin
    a_prop = apsf(PSFs.MethodPropagate, sz,pp, sampling=sampling);
    a_sincR = apsf(PSFs.MethodSincR, sz,pp, sampling=sampling);
    a_shell = apsf(PSFs.MethodShell, sz,pp, sampling=sampling);
    a_iter = apsf(PSFs.MethodPropagateIterative, sz,pp, sampling=sampling);
    # @vt a_prop a_sincR a_iter
end

@testset "vectorial pupil" begin
    ppv = PSFParams(pol=pol_x)
    pupil = pupil_xyz(sz, ppv, sampling)
end

return


text:
# functions that return arrays
get_current_x(relation) = value.(collect(keys(relation)))
get_current_F(relation) = value.(collect(values(relation)))
get_variables(relation) = collect(keys(relation))

#=
function createSolutionDict(m::Model)
    relation = getVIPData(m).relation

    sol_dict = Dict()
    for variable in keys(relation)
        xstar = value.(variable)
        Fstar = value.(relation[variable])
        sol_dict[variable] = value.(variable)
    end

    return sol_dict
end
=#

function _fixed_point(m, step_size, tolerance, max_iter)
    relation = getVIPData(m).relation

    # initialization
    var = get_variables(relation)
    xk = get_current_x(relation)
    fk = get_current_F(relation)
    tau = step_size

    for i=1:max_iter
        @objective(m, Min,
            sum( ( var[j] - (xk[j]-tau*fk[j]) )^2 for j in 1:length(var) )
        )

        optimize!(m)

        new_xk = get_current_x(relation)
        new_fk = get_current_F(relation)
        rel_err = norm(xk-new_xk)  #/ norm(xk)

        @printf("iteration %4d: rel_err=%15.12f\n", i, rel_err)

        if  rel_err < tolerance
            break
        end

        xk = new_xk
        fk = new_fk
    end

    return xk, fk
end

# Extragradient Section 12.1.2 of Facchinei and Pang (2003)
function _extra_gradient(m, step_size, tolerance, max_iter)
    relation = getVIPData(m).relation

    # initialization
    var = get_variables(relation)
    xk = get_current_x(relation)
    fk = get_current_F(relation)
    tau = step_size

    for i=1:max_iter
        @objective(m, Min,
            sum( ( var[j] - (xk[j]-tau*fk[j]) )^2 for j in 1:length(var) )
        )
        optimize!(m)
        status = termination_status(m)

        mid_fk = get_current_F(relation)

        @objective(m, Min,
            sum( ( var[j] - (xk[j]-tau*mid_fk[j]) )^2 for j in 1:length(var) )
        )
        optimize!(m)
        status = termination_status(m)

        new_xk = get_current_x(relation)
        new_fk = get_current_F(relation)
        rel_err = norm(xk-new_xk) / norm(xk)

        @printf("iteration %4d: rel_err=%15.12f\n", i, rel_err)

        if  rel_err < tolerance
            break
        end

        xk = new_xk
        fk = new_fk
    end

    return xk, fk
end

#=
# Hyperplane Projection Algorithm Section 12.1.3 of Facchinei and Pang (2003)
function _hyperplane(m, x0, step_size, tolerance, max_iter)

    F = getVIPData(m).F
    var = getVIPData(m).var
    @assert length(F) == length(var)
    dim = length(F)

    # initialization
    xk = x0
    setvalue(var, x0)
    fk = value.(F)
    tau = step_size
    sigma = 0.5

    for k=1:max_iter

        # yk = P_X ( xk - tau F(xk) )
        @NLobjective(m, Min, sum{ ( var[i] - (xk[i]-tau*fk[i]) )^2, i=1:dim} )
        optimize!(m)
        status = termination_status(m)
        yk = value.(var)

        # find the smallest nonnegative integer i such that...
        # F(2^(-i)yk + (1-2^(-i))xk) ^T (xk - yk) >= sigma / tau || xk - yk ||^2
        ik = Inf
        for i=0:10000
            # println("------------------------")
            # println("i=", i)
            # println("var=", 2.0^(-i)*yk + (1-2.0^(-i))*xk)

            setvalue(var, 2.0^(-i)*yk + (1-2.0^(-i))*xk)
            FF = value.(F)

            # println("FF=", FF)
            # println("dot(FF, xk-yk)=", dot(FF, xk-yk))
            # println("sigma / tau * dot(xk-yk, xk-yk)=", sigma / tau * dot(xk-yk, xk-yk))

            if dot(FF, xk-yk) >= sigma / tau * dot(xk-yk, xk-yk)
                ik = i
                break
            end
        end

        zk = 2.0^(-ik) * yk + (1-2.0^(-ik)) * xk
        setvalue(var, zk)
        Fzk = value.(F)

        wk = xk - dot(Fzk, xk-zk) / dot(Fzk, Fzk) * Fzk

        @NLobjective(m, Min, sum{ ( var[i] - wk[i] )^2, i=1:dim} )
        optimize!(m)
        status = termination_status(m)
        new_xk = value.(var)
        new_fk = value.(F)

        rel_err = norm(xk-new_xk) / norm(xk)

        @printf("iteration %4d: rel_err=%15.12f\n", k, rel_err)

        if  rel_err < tolerance
            break
        end

        xk = new_xk
        fk = new_fk
    end

    return xk, fk
end
=#


text:
module MultivariateMixtures

using LinearAlgebra
using Statistics: mean, var, cov, covm
using StatsBase: fit, ConvergenceException
using MultivariateStats: PPCA, FactorAnalysis, outdim
using Distributions
using Distributions: log2π
using Clustering: kmeans, assignments

import Base: show, length, convert
import Distributions: probs, ncomponents, components, component, logpdf, logpdf!
import Statistics: mean

export fit_mm, FinateMixtureModel

include("types.jl")
include("utils.jl")
include("gmm.jl")
include("mfa.jl")

# move to MultivariateStats
length(fa::FactorAnalysis) = outdim(fa)

# auxiliary functions
convert(::Type{MvNormal}, fa::FactorAnalysis) = MvNormal(mean(fa), cov(fa))
logpdf(fa::FactorAnalysis{T}, x::AbstractVector{T}) where {T<:Real} =
    logpdf(convert(MvNormal, fa), x)
logpdf(fa::FactorAnalysis{T}, X::AbstractMatrix{T}) where {T<:Real} =
    logpdf!(Vector{T}(undef, size(X,2)), fa, X)
logpdf!(r::AbstractArray{T}, fa::FactorAnalysis{T}, X::AbstractMatrix{T}) where {T<:Real} =
    logpdf!(r, convert(MvNormal, fa), X)
predict(classifier::Union{MixtureModel,FinateMixtureModel}, data::AbstractVector) =
    findmax(componentwise_logpdf(classifier, data))[end]
predict(classifier::Union{MixtureModel,FinateMixtureModel}, data::AbstractMatrix; dims=2) =
    map(ci->ci[2], findmax(componentwise_logpdf(classifier, (dims == 1 ? data' : data)), dims=2)[end])[:]

end


text:
@testset "Loo" begin
    result = @suppress Escape.run(
        Escape.HLAModel{4}(), Escape.HLADataset("Test").data[1],
        iter = 10, warmup = 10, chains = 4, seed = 123
    )

    sf = Escape.stanfit(result)
    p = extract(sf)

    @test Escape.theta_i(sf, p, 1) isa Vector{Float64}
    @test Escape.pointwise_loglikelihoods(sf, p, 1) isa Vector{Vector{Float64}}
    @test @suppress Escape.loo(result) isa Loo.LooResult
end


text:
abstract PCModification

abstract Rotation <: PCModification
abstract Reflection <: PCModification

abstract CWRotation <: Rotation
abstract counterCWRotation <: Rotation

abstract NoModification <: PCModification

#######################################################
## resolving pdf and cdf arguments for modifications ##
#######################################################

function resolve_args(mod::PCModification,
                      u1::FloatVec, u2::FloatVec)
    error("pdf modification still needs to be implemented for this
    type of copula modification")
end

function resolve_args(mod::PCModification, u1::Float64, u2::Float64)
    return resolve_args(mod, [u1], [u2])
end

function resolve_pval(mod::PCModification,
                      u1::FloatVec, u2::FloatVec)
    error("cdf modification still needs to be implemented for this
    type of copula modification")
end

function resolve_pval(mod::PCModification, u1::Float64, u2::Float64)
    return resolve_pval(mod, [u1], [u2])
end


text:
## Auxiliary functions
export acc

@inline acc(yt,yp) = count(x->x==true, yt .== yp)/length(yt)

# used in linear and voted perceptron
@inline sinal(x) = ( x>=0 ? 1.0 : 0.0 )
# used in kernel perceptron
@inline sign(val) = ( val >=0 ? 1.0 : -1.0 )

## checks PLS input data and params
function check_data(X::Matrix{T},Y::Union{Vector{T},Matrix{T}}) where T<:AbstractFloat
    !isempty(X) ||
        throw(DimensionMismatch("Empty input data (X)."))
    !isempty(Y) ||
        throw(DimensionMismatch("Empty target data (Y)."))
    size(X, 1) == size(Y, 1) ||
        throw(DimensionMismatch("Incompatible number of rows of input data (X) and target data (Y)."))
end

function check_data(X::Matrix{T},nfeatures::Int) where T<:AbstractFloat
    !isempty(X) ||
        throw(DimensionMismatch("Empty input data (X)."))
    size(X, 2) == nfeatures ||
        throw(DimensionMismatch("Incompatible number of columns of input data (X) and original training X columns."))
end

function check_params(kernel::AbstractString,mode::AbstractString)
    kernel in ["rbf","linear"] || error("kernel must be 'linear' or 'rbf'")
    mode   in ["kernel","linear","voted","averaged"] || error("mode must be 'linear' or 'kernel' or 'voted' or 'averaged'")
end

## checks constant columns
check_constant_cols(X::Matrix{T}) where {T<:AbstractFloat} = size(X,1)>1 && !any(all(X .== X[1,:]',dims=1)) || error("You must remove constant columns of input data (X) before train")
check_constant_cols(Y::Vector{T}) where {T<:AbstractFloat} = length(Y)>1 && length(unique(Y)) > 1 || error("Your target values are constant. All values are equal to $(Y[1])")

## Preprocessing data using z-score statistics. this is due to the fact that if X and Y are z-scored, than X'Y returns for W vector a pearson correlation for each element! :)
centralize_data(D::AbstractArray{T}, m::AbstractArray{T}, s::AbstractArray{T}) where {T<:AbstractFloat}  = (D .-m)./s
centralize_data(D::Vector{T}, m::T, s::T)  where {T<:AbstractFloat} = (D -m)/s

decentralize_data(D::AbstractArray{T}, m::AbstractArray{T}, s::AbstractArray{T}) where {T<:AbstractFloat} = D .*s .+m
decentralize_data(D::Vector{T}, m::T, s::T)  where {T<:AbstractFloat}               = D *s +m

check_linear_binary_labels(Y::Vector{T}) where {T<:Number} = length(setdiff(Y,[0.,1.]))==0 || error("Your target values must be 1 and 0 only. You have $(unique(Y)) distinct labels")


text:
using GLFW, ModernGL

# whenever the window size changed (by OS or user resize) this callback function executes
function framebuffer_size_callback(window::GLFW.Window, width::Cint, height::Cint)
    # make sure the viewport matches the new window dimensions; note that width and
    # height will be significantly larger than specified on retina displays.
	glViewport(0, 0, width, height)
end

# process all input: query GLFW whether relevant keys are pressed/released this frame and react accordingly
function processInput(window::GLFW.Window)
    GLFW.GetKey(window, GLFW.KEY_ESCAPE) == GLFW.PRESS && GLFW.SetWindowShouldClose(window, true)
end

# misc. config
const SCR_WIDTH = 800
const SCR_HEIGHT = 600

const vertexShaderSource = """
	#version 330 core
	layout (location = 0) in vec3 aPos;
	void main(void)
	{
	    gl_Position = vec4(aPos.x, aPos.y, aPos.z, 1.0);
	}"""

const fragmentShader1Source = """
	#version 330 core
	out vec4 FragColor;
	void main(void)
	{
	    FragColor = vec4(1.0f, 0.5f, 0.2f, 1.0f);
	}"""

const fragmentShader2Source = """
	#version 330 core
	out vec4 FragColor;
	void main(void)
	{
	    FragColor = vec4(1.0f, 1.0f, 0.0f, 1.0f);
	}"""

GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 3)
GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE)

# create window
window = GLFW.CreateWindow(SCR_WIDTH, SCR_HEIGHT, "LearnOpenGL")
window == C_NULL && error("Failed to create GLFW window.")
GLFW.MakeContextCurrent(window)
GLFW.SetFramebufferSizeCallback(window, framebuffer_size_callback)

# build and compile shader program
vertexShader = glCreateShader(GL_VERTEX_SHADER)
fragmentShaderOrange = glCreateShader(GL_FRAGMENT_SHADER)
fragmentShaderYellow = glCreateShader(GL_FRAGMENT_SHADER)
shaderProgramOrange = glCreateProgram()
shaderProgramYellow = glCreateProgram()

glShaderSource(vertexShader, 1, Ptr{GLchar}[pointer(vertexShaderSource)], C_NULL)
glCompileShader(vertexShader)
glShaderSource(fragmentShaderOrange, 1, Ptr{GLchar}[pointer(fragmentShader1Source)], C_NULL)
glCompileShader(fragmentShaderOrange)
glShaderSource(fragmentShaderYellow, 1, Ptr{GLchar}[pointer(fragmentShader2Source)], C_NULL)
glCompileShader(fragmentShaderYellow)

# link shaders
glAttachShader(shaderProgramOrange, vertexShader)
glAttachShader(shaderProgramOrange, fragmentShaderOrange)
glLinkProgram(shaderProgramOrange)

glAttachShader(shaderProgramYellow, vertexShader)
glAttachShader(shaderProgramYellow, fragmentShaderYellow)
glLinkProgram(shaderProgramYellow)

# set up vertex data (and buffer(s)) and configure vertex attributes
firstTriangle = GLfloat[-0.9,  0.5, 0.0,    # left
                         0.0, -0.5, 0.0,    # right
                         0.45, 0.5, 0.0]    # top
secondTriangle = GLfloat[0.0, -0.5, 0.0,    # left
						 0.9, -0.5, 0.0,    # right
					     0.45, 0.5, 0.0]    # top
vbos = zeros(GLuint, 2)
vaos = zeros(GLuint, 2)
glGenVertexArrays(2, vaos)
glGenBuffers(2, vbos)
# first triangle setup
glBindVertexArray(vaos[1])
glBindBuffer(GL_ARRAY_BUFFER, vbos[1])
glBufferData(GL_ARRAY_BUFFER, sizeof(firstTriangle), firstTriangle, GL_STATIC_DRAW)
glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 3 * sizeof(GLfloat), Ptr{Cvoid}(0))
glEnableVertexAttribArray(0)
# second triangle setup
glBindVertexArray(vaos[2])
glBindBuffer(GL_ARRAY_BUFFER, vbos[2])
glBufferData(GL_ARRAY_BUFFER, sizeof(secondTriangle), secondTriangle, GL_STATIC_DRAW)
glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 3 * sizeof(GLfloat), Ptr{Cvoid}(0))
glEnableVertexAttribArray(0)

# uncomment this call to draw in wireframe polygons.
# glPolygonMode(GL_FRONT_AND_BACK, GL_LINE)

# render loop
while !GLFW.WindowShouldClose(window)
    processInput(window)

    # render
    glClearColor(0.2, 0.3, 0.3, 1.0)
    glClear(GL_COLOR_BUFFER_BIT)

	# draw first triangle
	glUseProgram(shaderProgramOrange)
    glBindVertexArray(vaos[1])
	glDrawArrays(GL_TRIANGLES, 0, 3)

	# draw  second triangle
	glUseProgram(shaderProgramYellow)
    glBindVertexArray(vaos[2])
	glDrawArrays(GL_TRIANGLES, 0, 3)

    # swap buffers and poll IO events
    GLFW.SwapBuffers(window)
    GLFW.PollEvents()
end

# optional: de-allocate all resources once they've outlived their purpose
glDeleteVertexArrays(2, vaos)
glDeleteBuffers(2, vbos)

GLFW.DestroyWindow(window)

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_name: Gnimuc/Videre


text:
using Documenter, AdventOfCode

makedocs(sitename="AdventOfCode.jl")
deploydocs(
    repo = "github.com/SebRollen/AdventOfCode.jl.git",
)


text:
"""
    SparseTabularMDP

An MDP object where states and actions are integers and the transition is represented by a list of sparse matrices.
This data structure can be useful to exploit in vectorized algorithm (e.g. see SparseValueIterationSolver).
The recommended way to access the transition and reward matrices is through the provided accessor functions: `transition_matrix` and `reward_vector`.

# Fields
- `T::Vector{SparseMatrixCSC{Float64, Int64}}` The transition model is represented as a vector of sparse matrices (one for each action). `T[a][s, sp]` the probability of transition from `s` to `sp` taking action `a`.
- `R::Array{Float64, 2}` The reward is represented as a matrix where the rows are states and the columns actions: `R[s, a]` is the reward of taking action `a` in sate `s`.
- `terminal_states::Set{Int64}` Stores the terminal states
- `discount::Float64` The discount factor

# Constructors

- `SparseTabularMDP(mdp::MDP)` : One can provide the matrices to the default constructor or one can construct a `SparseTabularMDP` from any discrete state MDP defined using the explicit interface.
Note that constructing the transition and reward matrices requires to iterate over all the states and can take a while.
To learn more information about how to define an MDP with the explicit interface please visit https://juliapomdp.github.io/POMDPs.jl/latest/explicit/ .
- `SparseTabularMDP(smdp::SparseTabularMDP; transition, reward, discount)` : This constructor returns a new sparse MDP that is a copy of the original smdp except for the field specified by the keyword arguments.

"""
struct SparseTabularMDP <: MDP{Int64, Int64}
    T::Vector{SparseMatrixCSC{Float64, Int64}} # T[a][s, sp]
    R::Array{Float64, 2} # R[s, a]
    terminal_states::Set{Int64}
    discount::Float64
end

function SparseTabularMDP(mdp::MDP)
    T = transition_matrix_a_s_sp(mdp)
    R = reward_s_a(mdp)
    ts = terminal_states_set(mdp)
    return SparseTabularMDP(T, R, ts, discount(mdp))
end

@POMDP_require SparseTabularMDP(mdp::MDP) begin
    P = typeof(mdp)
    S = statetype(P)
    A = actiontype(P)
    @req discount(::P)
    @subreq ordered_states(mdp)
    @subreq ordered_actions(mdp)
    @req transition(::P,::S,::A)
    @req reward(::P,::S,::A,::S)
    @req stateindex(::P,::S)
    @req actionindex(::P, ::A)
    @req actions(::P, ::S)
    as = actions(mdp)
    ss = states(mdp)
    @req length(::typeof(as))
    @req length(::typeof(ss))
    a = first(as)
    s = first(ss)
    dist = transition(mdp, s, a)
    D = typeof(dist)
    @req support(::D)
    @req pdf(::D,::S)
end

function SparseTabularMDP(mdp::SparseTabularMDP;
                          transition::Union{Nothing, Vector{SparseMatrixCSC{Float64, Int64}}} = nothing,
                          reward::Union{Nothing, Array{Float64, 2}} = nothing,
                          discount::Union{Nothing, Float64} = nothing,
                          terminal_states::Union{Nothing, Set{Int64}} = nothing)
    T = transition != nothing ? transition : mdp.T
    R = reward != nothing ? reward : mdp.R
    d = discount != nothing ? discount : mdp.discount
    ts = terminal_states != nothing ? terminal_states : mdp.terminal_states
    return SparseTabularMDP(T, R, ts, d)
end

"""
    SparseTabularPOMDP

A POMDP object where states and actions are integers and the transition and observation distributions are represented by lists of sparse matrices.
This data structure can be useful to exploit in vectorized algorithms to gain performance (e.g. see SparseValueIterationSolver).
The recommended way to access the transition, reward, and observation matrices is through the provided accessor functions: `transition_matrix`, `reward_vector`, `observation_matrix`.

# Fields
- `T::Vector{SparseMatrixCSC{Float64, Int64}}` The transition model is represented as a vector of sparse matrices (one for each action). `T[a][s, sp]` the probability of transition from `s` to `sp` taking action `a`.
- `R::Array{Float64, 2}` The reward is represented as a matrix where the rows are states and the columns actions: `R[s, a]` is the reward of taking action `a` in sate `s`.
- `O::Vector{SparseMatrixCSC{Float64, Int64}}` The observation model is represented as a vector of sparse matrices (one for each action). `O[a][sp, o]` is the probability of observing `o` from state `sp` after having taken action `a`.
- `terminal_states::Set{Int64}` Stores the terminal states
- `discount::Float64` The discount factor

# Constructors

- `SparseTabularPOMDP(pomdp::POMDP)` : One can provide the matrices to the default constructor or one can construct a `SparseTabularPOMDP` from any discrete state MDP defined using the explicit interface.
Note that constructing the transition and reward matrices requires to iterate over all the states and can take a while.
To learn more information about how to define an MDP with the explicit interface please visit https://juliapomdp.github.io/POMDPs.jl/latest/explicit/ .
- `SparseTabularPOMDP(spomdp::SparseTabularMDP; transition, reward, observation, discount)` : This constructor returns a new sparse POMDP that is a copy of the original smdp except for the field specified by the keyword arguments.

"""
struct SparseTabularPOMDP <: POMDP{Int64, Int64, Int64}
    T::Vector{SparseMatrixCSC{Float64, Int64}} # T[a][s, sp]
    R::Array{Float64, 2} # R[s,sp]
    O::Vector{SparseMatrixCSC{Float64, Int64}} # O[a][sp, o]
    terminal_states::Set{Int64}
    discount::Float64
end

function SparseTabularPOMDP(pomdp::POMDP)
    T = transition_matrix_a_s_sp(pomdp)
    R = reward_s_a(pomdp)
    O = observation_matrix_a_sp_o(pomdp)
    ts = terminal_states_set(pomdp)
    return SparseTabularPOMDP(T, R, O, ts, discount(pomdp))
end

@POMDP_require SparseTabularPOMDP(pomdp::POMDP) begin
    P = typeof(pomdp)
    S = statetype(P)
    A = actiontype(P)
    O = obstype(P)
    @req discount(::P)
    @subreq ordered_states(pomdp)
    @subreq ordered_actions(pomdp)
    @subreq ordered_observations(pomdp)
    @req transition(::P,::S,::A)
    @req reward(::P,::S,::A,::S)
    @req observation(::P, ::A, ::S)
    @req stateindex(::P,::S)
    @req actionindex(::P, ::A)
    @req actions(::P, ::S)
    @req observations(::P)
    @req obsindex(::P, ::O)
    as = actions(pomdp)
    ss = states(pomdp)
    @req length(::typeof(as))
    @req length(::typeof(ss))
    a = first(as)
    s = first(ss)
    dist = transition(pomdp, s, a)
    D = typeof(dist)
    @req support(::D)
    @req pdf(::D,::S)
    odist = observation(pomdp, a, s)
    OD = typeof(odist)
    @req support(::OD)
    @req pdf(::OD, ::O)
end

function SparseTabularPOMDP(pomdp::SparseTabularPOMDP;
                          transition::Union{Nothing, Vector{SparseMatrixCSC{Float64, Int64}}} = nothing,
                          reward::Union{Nothing, Array{Float64, 2}} = nothing,
                          observation::Union{Nothing, Vector{SparseMatrixCSC{Float64, Int64}}} = nothing,
                          discount::Union{Nothing, Float64} = nothing,
                          terminal_states::Union{Nothing, Set{Int64}} = nothing)
    T = transition != nothing ? transition : pomdp.T
    R = reward != nothing ? reward : pomdp.R
    d = discount != nothing ? discount : pomdp.discount
    O = observation != nothing ? transition : pomdp.O
    ts = terminal_states != nothing ? terminal_states : pomdp.terminal_states
    return SparseTabularPOMDP(T, R, O, ts, d)
end

const SparseTabularProblem = Union{SparseTabularMDP, SparseTabularPOMDP}

function transition_matrix_a_s_sp(mdp::Union{MDP, POMDP})
    # Thanks to zach
    na = length(actions(mdp))
    state_space = states(mdp)
    ns = length(state_space)
    transmat_row_A = [Int[] for _ in 1:n_actions(mdp)]
    transmat_col_A = [Int[] for _ in 1:n_actions(mdp)]
    transmat_data_A = [Float64[] for _ in 1:n_actions(mdp)]

    for s in states(mdp)
        si = stateindex(mdp, s)
        for a in actions(mdp, s)
            ai = actionindex(mdp, a)
            if isterminal(mdp, s) # if terminal, there is a probability of 1 of staying in that state
                push!(transmat_row_A[ai], si)
                push!(transmat_col_A[ai], si)
                push!(transmat_data_A[ai], 1.0)
            else
                td = transition(mdp, s, a)
                for (sp, p) in weighted_iterator(td)
                    if p > 0.0
                        spi = stateindex(mdp, sp)
                        push!(transmat_row_A[ai], si)
                        push!(transmat_col_A[ai], spi)
                        push!(transmat_data_A[ai], p)
                    end
                end
            end
        end
    end
    transmats_A_S_S2 = [sparse(transmat_row_A[a], transmat_col_A[a], transmat_data_A[a], ns, ns) for a in 1:na]
    # if an action is not valid from a state, the transition is 0.0 everywhere
    # @assert all(all(sum(transmats_A_S_S2[a], dims=2) .≈ ones(ns)) for a in 1:na) "Transition probabilities must sum to 1"
    return transmats_A_S_S2
end

function reward_s_a(mdp::Union{MDP, POMDP})
    state_space = states(mdp)
    action_space = actions(mdp)
    reward_S_A = fill(-Inf, (length(state_space), length(action_space))) # set reward for all actions to -Inf unless they are in actions(mdp, s)
    for s in state_space
        if isterminal(mdp, s)
            reward_S_A[stateindex(mdp, s), :] .= 0.0
        else
            for a in actions(mdp, s)
                td = transition(mdp, s, a)
                r = 0.0
                for (sp, p) in weighted_iterator(td)
                    if p > 0.0
                        r += p*reward(mdp, s, a, sp)
                    end
                end
                reward_S_A[stateindex(mdp, s), actionindex(mdp, a)] = r
            end
        end
    end
    return reward_S_A
end

function terminal_states_set(mdp::Union{MDP, POMDP})
    ts = Set{Int64}()
    for s in states(mdp)
        if isterminal(mdp, s)
            si = stateindex(mdp, s)
            push!(ts, si)
        end
    end
    return ts
end

function observation_matrix_a_sp_o(pomdp::POMDP)
    state_space = states(pomdp)
    action_space = actions(pomdp)
    obs_space = observations(pomdp)
    na = length(action_space)
    ns = length(state_space)
    no = length(obs_space)
    obsmat_row_A = [Int[] for _ in 1:na]
    obsmat_col_A = [Int[] for _ in 1:na]
    obsmat_data_A = [Float64[] for _ in 1:na]

    for sp in states(pomdp)
        spi = stateindex(pomdp, sp)
        for a in actions(pomdp)
            ai = actionindex(pomdp, a)
            od = observation(pomdp, a, sp)
            for (o, p) in weighted_iterator(od)
                if p > 0.0
                    oi = obsindex(pomdp, o)
                    push!(obsmat_row_A[ai], spi)
                    push!(obsmat_col_A[ai], oi)
                    push!(obsmat_data_A[ai], p)
                end
            end
        end
    end
    obsmats_A_SP_O = [sparse(obsmat_row_A[a], obsmat_col_A[a], obsmat_data_A[a], ns, ns) for a in 1:na]
    @assert all(all(sum(obsmats_A_SP_O[a], dims=2) .≈ ones(no)) for a in 1:na) "Observation probabilities must sum to 1"
    return obsmats_A_SP_O
end

# MDP and POMDP common methods

POMDPs.states(p::SparseTabularProblem) = 1:size(p.T[1], 1)
POMDPs.actions(p::SparseTabularProblem) = 1:size(p.T, 1)
POMDPs.actions(p::SparseTabularProblem, s::Int64) = [a for a in actions(p) if sum(transition_matrix(p, a)) ≈ size(p.T[1], 1)]

POMDPs.stateindex(::SparseTabularProblem, s::Int64) = s
POMDPs.actionindex(::SparseTabularProblem, a::Int64) = a

POMDPs.discount(p::SparseTabularProblem) = p.discount

POMDPs.transition(p::SparseTabularProblem, s::Int64, a::Int64) = SparseCat(findnz(p.T[a][s, :])...) # XXX not memory efficient

POMDPs.reward(p::SparseTabularProblem, s::Int64, a::Int64) = p.R[s, a]

POMDPs.isterminal(p::SparseTabularProblem, s::Int64) = s ∈ p.terminal_states

"""
    transition_matrix(p::SparseTabularProblem, a)
Accessor function for the transition model of a sparse tabular problem.
It returns a sparse matrix containing the transition probabilities when taking action a: T[s, sp] = Pr(sp | s, a).
"""
transition_matrix(p::SparseTabularProblem, a) = p.T[a]

"""
    transition_matrices(p::SparseTabularProblem)
Accessor function for the transition model of a sparse tabular problem.
It returns a list of sparse matrices for each action of the problem.
"""
transition_matrices(p::SparseTabularProblem) = p.T

"""
    reward_vector(p::SparseTabularProblem, a)
Accessor function for the reward function of a sparse tabular problem.
It returns a vector containing the reward for all the states when taking action a: R(s, a).
The length of the return vector is equal to the number of states.
"""
reward_vector(p::SparseTabularProblem, a) = view(p.R, :, a)

"""
    reward_matrix(p::SparseTabularProblem)
Accessor function for the reward matrix R[s, a] of a sparse tabular problem.
"""
reward_matrix(p::SparseTabularProblem) = p.R

# POMDP only methods

POMDPs.observations(p::SparseTabularPOMDP) = 1:size(p.O[1], 2)

POMDPs.observation(p::SparseTabularPOMDP, a::Int64, sp::Int64) = SparseCat(findnz(p.O[a][sp, :])...)

POMDPs.obsindex(p::SparseTabularPOMDP, o::Int64) = o

"""
    observation_matrix(p::SparseTabularPOMDP, a::Int64)
Accessor function for the observation model of a sparse tabular POMDP.
It returns a sparse matrix containing the observation probabilities when having taken action a: O[sp, o] = Pr(o | sp, a).
"""
observation_matrix(p::SparseTabularPOMDP, a::Int64) = p.O[a]

"""
    observation_matrices(p::SparseTabularPOMDP)
Accessor function for the observation model of a sparse tabular POMDP.
It returns a list of sparse matrices for each action of the problem.
"""
observation_matrices(p::SparseTabularPOMDP) = p.O


text:
###################
### OnlineStats ###
###################
"""
$(TYPEDEF)

# Usage

    Skip(b::Int, stat::OnlineStat)

Skips the first `b` observations before passing them on to `stat`.
"""
mutable struct Skip{T, O<:OnlineStat{T}} <: OnlineStat{T}
    b::Int
    current_index::Int
    stat::O
end

Skip(b::Int, stat) = Skip(b, 0, stat)

OnlineStats.nobs(o::Skip) = OnlineStats.nobs(o.stat)
OnlineStats.value(o::Skip) = OnlineStats.value(o.stat)
function OnlineStats._fit!(o::Skip, x::Real)
    if o.current_index > o.b
        OnlineStats._fit!(o.stat, x)
    end
    o.current_index += length(x)

    return o
end

Base.show(io::IO, o::Skip) = print(
    io,
    "Skip ($(o.b)): current_index=$(o.current_index) | stat=$(o.stat)`"
)

"""
$(TYPEDEF)

# Usage

    Thin(b::Int, stat::OnlineStat)

Thins `stat` with an interval `b`, i.e. only passes every b-th observation to `stat`.
"""
mutable struct Thin{T, O<:OnlineStat{T}} <: OnlineStat{T}
    b::Int
    current_index::Int
    stat::O
end

Thin(b::Int, stat) = Thin(b, 0, stat)

OnlineStats.nobs(o::Thin) = OnlineStats.nobs(o.stat)
OnlineStats.value(o::Thin) = OnlineStats.value(o.stat)
function OnlineStats._fit!(o::Thin, x::Real)
    if (o.current_index % o.b) == 0
        OnlineStats._fit!(o.stat, x)
    end
    o.current_index += length(x)

    return o
end

Base.show(io::IO, o::Thin) = print(
    io,
    "Thin ($(o.b)): current_index=$(o.current_index) | stat=$(o.stat)`"
)

"""
$(TYPEDEF)

# Usage

    WindowStat(b::Int, stat::O) where {O <: OnlineStat}

"Wraps" `stat` in a `MovingWindow` of length `b`.

`value(o::WindowStat)` will then return an `OnlineStat` of the same type as
`stat`, which is *only* fitted on the batched data contained in the `MovingWindow`.

"""
struct WindowStat{T, O} <: OnlineStat{T}
    window::MovingWindow{T}
    stat::O
end

WindowStat(b::Int, T::Type, o) = WindowStat{T, typeof(o)}(MovingWindow(b, T), o)
WindowStat(b::Int, o::OnlineStat{T}) where {T} = WindowStat{T, typeof(o)}(
    MovingWindow(b, T), o
)

# Proxy methods to the window
OnlineStats.nobs(o::WindowStat) = OnlineStats.nobs(o.window)
OnlineStats._fit!(o::WindowStat, x) = OnlineStats._fit!(o.window, x)

function OnlineStats.value(o::WindowStat{<:Any, <:OnlineStat})
    stat_new = deepcopy(o.stat)
    fit!(stat_new, OnlineStats.value(o.window))
    return stat_new
end

function OnlineStats.value(o::WindowStat{<:Any, <:Function})
    stat_new = o.stat()
    fit!(stat_new, OnlineStats.value(o.window))
    return stat_new
end


text:
module DCCA

include("calculs.jl")

"""
    rhoDCCA(x,y; box_start = 3, box_stop = div(length(x),10), nb_pts = 30)

Performs the DCCA analysis of `x` and `y`. The default analysis starts with a window size of 3 up to one tenth of the total length of `x` for statistical reasons.

returns the different window sizes used for the analysis, and the associated dcca coefficients.
"""
function rhoDCCA(data1::Array{Float64,1},data2::Array{Float64,1}; box_start = 3, box_stop = div(length(data1),10), nb_pts = 30)
    x = log_space(box_start,box_stop,nb_pts)
    y = dcca(data1,data2; box_start = box_start, box_stop = box_stop, nb_pts = nb_pts)
    return x, y
end

function rhoDCCA_CI(x,y)
    lower, upper = zeros(500), zeros(500)
    for i in 1:500
        tmp = dcca(shuffle(x), shuffle(y))
        lower[i], upper[i] = minimum(tmp), maximum(tmp)
    end
    sort!(lower)
    sort!(upper)
    return lower[25], upper[475]
end

export rhoDCCA, rhoDCCA_CI

end

max_stars_repo_name: Datseris/VelocityMTD

max_issues_repo_name: Datseris/VelocityMTD

max_forks_repo_name: Datseris/VelocityMTD


text:
function to_roman(number)

end


text:
# This Map has non-trivial domain, thus we need the define testargs
"""
    Reindex(values) -> Map
"""
struct Reindex{A} <: Map
  values::A
end

function testargs(k::Reindex,i)
  @check length(k.values) !=0 "This map has empty domain"
  (one(i),)
end
function testargs(k::Reindex,i::Integer...)
  @check length(k.values) !=0 "This map has empty domain"
  map(one,i)
end
function return_value(k::Reindex,i...)
  length(k.values)!=0 ? evaluate(k,testargs(k,i...)...) : testitem(k.values)
end
return_cache(k::Reindex,i...) = array_cache(k.values)
evaluate!(cache,k::Reindex,i...) = getindex!(cache,k.values,i...)
evaluate(k::Reindex,i...) = k.values[i...]

#"""
#    reindex(i_to_v::AbstractArray, j_to_i::AbstractArray)
#"""
#function reindex(i_to_v::AbstractArray, j_to_i::AbstractArray)
#  lazy_map(Reindex(i_to_v),j_to_i)
#end

function lazy_map(k::Reindex{<:Fill},::Type{T}, j_to_i::AbstractArray) where T
  v = k.values.value
  Fill(v,size(j_to_i)...)
end

function lazy_map(k::Reindex{<:CompressedArray},::Type{T}, j_to_i::AbstractArray) where T
  i_to_v = k.values
  values = i_to_v.values
  ptrs = lazy_map(Reindex(i_to_v.ptrs),j_to_i)
  CompressedArray(values,ptrs)
end

function lazy_map(k::Reindex{<:LazyArray},::Type{T},j_to_i::AbstractArray) where T
  i_to_maps = k.values.maps
  i_to_args = k.values.args
  j_to_maps = lazy_map(Reindex(i_to_maps),eltype(i_to_maps),j_to_i)
  j_to_args = map(i_to_fk->lazy_map(Reindex(i_to_fk),eltype(i_to_fk),j_to_i), i_to_args)
  LazyArray(T,j_to_maps,j_to_args...)
end

# This optimization is important for surface-coupled problems
function lazy_map(k::Reindex{<:LazyArray{<:Fill{<:PosNegReindex}}},::Type{T},j_to_i::AbstractArray) where T
  i_to_iposneg = k.values.args[1]
  ipos_to_value = k.values.maps.value.values_pos
  ineg_to_value = k.values.maps.value.values_neg
  if aligned_with_pos(i_to_iposneg,j_to_i,length(ipos_to_value))
    ipos_to_value
  elseif aligned_with_neg(i_to_iposneg,j_to_i,length(ineg_to_value))
    ineg_to_value
  elseif all_in_pos(i_to_iposneg,j_to_i)
    j_to_ipos = lazy_map(Reindex(get_array(i_to_iposneg)),j_to_i)
    j_to_value = lazy_map(Reindex(ipos_to_value),j_to_ipos)
  elseif all_in_neg(i_to_iposneg,j_to_i)
    j_to_ineg = lazy_map(Reindex(get_array(i_to_iposneg)),j_to_i)
    j_to_value = lazy_map(Reindex(ineg_to_value),lazy_map(ineg->-ineg,j_to_ineg))
  else
    j_to_iposneg = lazy_map(Reindex(get_array(i_to_iposneg)),j_to_i)
    j_to_value = lazy_map(PosNegReindex(ipos_to_value,ineg_to_value),j_to_iposneg)
  end
end

function _lazy_map_evaluate_posneg(
  i_to_jposneg::LazyArray{<:Fill{<:Reindex}},::Type{T},i_a,i_x::Fill) where T
  i_y = CompressedArray([i_x.value],Fill(Int32(1),length(i_x)))
  _lazy_map_evaluate_posneg(i_to_jposneg,T,i_a,i_y)
end

function _lazy_map_evaluate_posneg(i_to_jposneg::LazyArray{<:Fill{<:Reindex}},::Type{T},i_a,i_x::CompressedArray) where T
  jpos_a = i_a.maps.value.values_pos
  jneg_a = i_a.maps.value.values_neg
  j_jposneg = i_to_jposneg.maps.value.values
  jpos_j, jneg_j = pos_and_neg_indices(j_jposneg)
  k_x = deepcopy(i_x.values)
  i_k = i_x.ptrs
  push!(k_x,testitem(i_x))
  i_j = i_to_jposneg.args[1]
  j_k = fill(Int32(length(k_x)),length(j_jposneg))
  j_k[i_j] .= i_k
  jpos_k = j_k[jpos_j]
  jneg_k = j_k[jneg_j]
  jpos_x = CompressedArray(k_x,jpos_k)
  jneg_x = CompressedArray(k_x,jneg_k)
  jpos_b = lazy_map(evaluate,jpos_a,jpos_x)
  jneg_b = lazy_map(evaluate,jneg_a,jneg_x)
  j_b = lazy_map(PosNegReindex(jpos_b,jneg_b),j_jposneg)
  i_b = lazy_map(Reindex(j_b),i_j)
  i_b
end

function lazy_map(k::Reindex{<:LazyArray{<:Fill{<:PosNegReindex}}},::Type{T},indices::IdentityVector) where T
  @check length(k.values) == length(indices)
  k.values
end

function lazy_map(k::Reindex{<:LazyArray},::Type{T},indices::IdentityVector) where T
  @check length(k.values) == length(indices)
  k.values
end

function lazy_map(k::Reindex{<:AbstractArray},::Type{T},indices::IdentityVector) where T
  @check length(k.values) == length(indices)
  k.values
end

function lazy_map(k::Reindex{<:Fill},::Type{T},indices::IdentityVector) where T
  @check length(k.values) == length(indices)
  k.values
end

function lazy_map(k::Reindex{<:CompressedArray},::Type{T},indices::IdentityVector) where T
  @check length(k.values) == length(indices)
  k.values
end

@propagate_inbounds function Base.setindex!(a::LazyArray{<:Fill{<:Reindex}},v,j::Integer)
  k = a.maps.value
  i_to_v = k.values
  j_to_i, = a.args
  i = j_to_i[j]
  i_to_v[i]=v
end


text:
# You reach another volcanically active part of the cave. It would be nice if
# you could do some kind of thermal imaging so you could tell ahead of time
# which caves are too hot to safely enter.
#
# Fortunately, the submarine seems to be equipped with a thermal camera! When
# you activate it, you are greeted with:
#
# Congratulations on your purchase! To activate this infrared thermal imaging
# camera system, please enter the code found on page 1 of the manual.
#
# Apparently, the Elves have never used this feature. To your surprise, you
# manage to find the manual; as you go to open it, page 1 falls out. It's a
# large sheet of transparent paper! The transparent paper is marked with random
# dots and includes instructions on how to fold it up (your puzzle input). For
# example:
#
# 6,10
# 0,14
# 9,10
# 0,3
# 10,4
# 4,11
# 6,0
# 6,12
# 4,1
# 0,13
# 10,12
# 3,4
# 3,0
# 8,4
# 1,10
# 2,14
# 8,10
# 9,0
#
# fold along y=7
# fold along x=5
#
# The first section is a list of dots on the transparent paper. 0,0 represents
# the top-left coordinate. The first value, x, increases to the right. The
# second value, y, increases downward. So, the coordinate 3,0 is to the right
# of 0,0, and the coordinate 0,7 is below 0,0. The coordinates in this example
# form the following pattern, where # is a dot on the paper and . is an empty,
# unmarked position:
#
# ...#..#..#.
# ....#......
# ...........
# #..........
# ...#....#.#
# ...........
# ...........
# ...........
# ...........
# ...........
# .#....#.##.
# ....#......
# ......#...#
# #..........
# #.#........
#
# Then, there is a list of fold instructions. Each instruction indicates a line
# on the transparent paper and wants you to fold the paper up (for horizontal
# y=... lines) or left (for vertical x=... lines). In this example, the first
# fold instruction is fold along y=7, which designates the line formed by all
# of the positions where y is 7 (marked here with -):
#
# ...#..#..#.
# ....#......
# ...........
# #..........
# ...#....#.#
# ...........
# ...........
# -----------
# ...........
# ...........
# .#....#.##.
# ....#......
# ......#...#
# #..........
# #.#........
#
# Because this is a horizontal line, fold the bottom half up. Some of the dots
# might end up overlapping after the fold is complete, but dots will never
# appear exactly on a fold line. The result of doing this fold looks like this:
#
# #.##..#..#.
# #...#......
# ......#...#
# #...#......
# .#.#..#.###
# ...........
# ...........
#
# Now, only 17 dots are visible.
#
# Notice, for example, the two dots in the bottom left corner before the
# transparent paper is folded; after the fold is complete, those dots appear in
# the top left corner (at 0,0 and 0,1). Because the paper is transparent, the
# dot just below them in the result (at 0,3) remains visible, as it can be seen
# through the transparent paper.
#
# Also notice that some dots can end up overlapping; in this case, the dots
# merge together and become a single dot.
#
# The second fold instruction is fold along x=5, which indicates this line:
#
# #.##.|#..#.
# #...#|.....
# .....|#...#
# #...#|.....
# .#.#.|#.###
# .....|.....
# .....|.....
#
# Because this is a vertical line, fold left:
#
# #####
# #...#
# #...#
# #...#
# #####
# .....
# .....
#
# The instructions made a square!
#
# The transparent paper is pretty big, so for now, focus on just completing the
# first fold. After the first fold in the example above, 17 dots are visible -
# dots that end up overlapping after the fold is completed count as a single
# dot.
#
# How many dots are visible after completing just the first fold instruction on
# your transparent paper?

# read data
iter = eachline("input.txt")
values = map(l -> parse.(Int, split(l, ',')), Iterators.takewhile(l -> l != "", iter))
folds = map(iter) do l
  (f, i) = split(l, '=')
  (f[end], parse(Int, i))
end

# build sheet - I really hate 1-based indexing
width = maximum(first, values) + 1
height = maximum(last, values) + 1
sheet = falses(width, height)
foreach(((x,y),) -> sheet[x+1,y+1] = true, values)

# fold
fold_along = folds[1][2]
if folds[1][1] == 'x'
  sheet = sheet[1:fold_along,:] .| sheet[end:-1:(fold_along+2),:]
else
  sheet = sheet[:,1:fold_along] .| sheet[:,end:-1:(fold_along+2)]
end

println("Count: ", count(sheet))


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module GBDTs

export GBDTParams, GBDTResult, induce_tree, apply_expr, classify

using ExprSearch, DecisionTrees
using RLESUtils, DataFrameSets, GitUtils, CPUTimeUtils, LogSystems, Observers,
    TreeIterators
using JLD

import RLESTypes.SymbolTable
import DecisionTrees: get_truth, get_splitter, get_split_labels, classify

#user should implemented these for their problem
apply_expr(problem::ExprProblem, ids::Vector{Int64}, expr) =
    error("user should override apply_expr()")

include("logdefs.jl")

type GBDTParams
    problem::ExprProblem
    num_data::Int64
    split_params::SearchParams
    maxdepth::Int64
    predict_type::Type #T1
    label_type::Type #T2
    logsys::LogSystem
end
GBDTParams(problem::ExprProblem, num_data::Int64, split_params::SearchParams,
    maxdepth::Int64, predict_type::Type, label_type::Type) =
    GBDTParams(problem, split_params,
    maxdepth, predict_type, label_type, logsystem())

type GBDTParamsSerial
    num_data::Int64
    split_params::SearchParams
    maxdepth::Int64
    predict_type::Type #T1
    label_type::Type #T2
end

function JLD.writeas(p::GBDTParams)
    GBDTParamsSerial(p.num_data, p.split_params, p.maxdepth, p.predict_type, p.label_type)
end

type GBDTResult
    params::GBDTParams
    dt_params::DTParams #need this for classify to work
    decision_tree::DecisionTree
end

function induce_tree(p::GBDTParams)
    @notify_observer(p.logsys.observer, "verbose1", ["Starting decision tree induction..."])
    @notify_observer(p.logsys.observer, "computeinfo", ["starttime", string(now())])

    tstart = CPUtime_start()
    dt_userargs = SymbolTable(:params=>p, :decision_id=>0)
    dt_params = DTParams(p.num_data, p.maxdepth, p.predict_type, p.label_type, dt_userargs)
    dtree = build_tree(dt_params)
    result = GBDTResult(p, dt_params, dtree)

    @notify_observer(p.logsys.observer, "result", [result])

    for node in tree_iter(dtree.root)
        @notify_observer(p.logsys.observer, "members", [node.id, string(node.members)])
    end

    #meta info
    @notify_observer(p.logsys.observer, "computeinfo", ["endtime",  string(now())])
    @notify_observer(p.logsys.observer, "computeinfo", ["hostname", gethostname()])
    @notify_observer(p.logsys.observer, "computeinfo", ["gitSHA",  get_SHA(dirname(@__FILE__))])
    @notify_observer(p.logsys.observer, "computeinfo", ["cpu_time",
        string(CPUtime_elapsed_s(tstart))])
    @notify_observer(p.logsys.observer, "parameters", ["num_data", p.num_data])
    @notify_observer(p.logsys.observer, "parameters", ["maxdepth", p.maxdepth])
    @notify_observer(p.logsys.observer, "parameters", ["predict_type", string(p.predict_type)])
    @notify_observer(p.logsys.observer, "parameters", ["label_type", string(p.label_type)])

    result
end

function get_truth(members::Vector{Int64}, dt_userargs::SymbolTable)
    p = dt_userargs[:params]
    labels(p.problem.data, members)
end

function get_splitter(members::Vector{Int64}, dt_userargs::SymbolTable)
    tstart = CPUtime_start()

    p = dt_userargs[:params]
    dt_userargs[:decision_id] += 1
    decision_id = dt_userargs[:decision_id]

    p.split_params.userargs[:ids] = members
    split_result = exprsearch(p.split_params, p.problem)

    @notify_observer(p.logsys.observer, "elapsed_cpu_s", [decision_id, CPUtime_elapsed_s(tstart)])
    @notify_observer(p.logsys.observer, "split_result_raw", [decision_id, split_result])
    @notify_observer(p.logsys.observer, "split_result", [decision_id,
        string(get_expr(split_result)), get_fitness(split_result)])

    split_result
end

function get_split_labels(result::SearchResult, members::Vector{Int64},
    dt_userargs::SymbolTable)

    p = dt_userargs[:params]
    expr = get_expr(result)
    apply_expr(p.problem, members, expr)
end
function get_split_labels(result::SearchResult, records::AbstractVector{DataFrame},
    dt_userargs::SymbolTable)

    p = dt_userargs[:params]
    expr = get_expr(result)
    apply_expr(p.problem, records, expr)
end

function classify{T}(result::GBDTResult, v::Vector{T})
   [classify(result, v[i]) for i=1:length(v)]
end
function classify(result::GBDTResult, x)
    classify(result.dt_params, result.decision_tree, x)
end

end #module

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
using ImGuiOpenGLBackend
using Test

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


text:
# SingleDecision #
# ========================== #
mutable struct SingleDecisionConstraintBridge{T, S <: MOI.AbstractScalarSet} <: MOIB.Constraint.AbstractBridge
    constraint::CI{MOI.SingleVariable, S}
    decision::SingleDecision
end

function MOIB.Constraint.bridge_constraint(::Type{SingleDecisionConstraintBridge{T,S}},
                                           model,
                                           f::SingleDecision,
                                           set::S) where {T, S <: MOI.AbstractScalarSet}
    # Perform the bridge mapping manually
    g = MOIB.bridged_variable_function(model, f.decision)
    mapped_variable = MOI.SingleVariable(only(g.terms).variable_index)
    # Add the bridged constraint
    constraint = MOI.add_constraint(model,
                                    mapped_variable,
                                    set)
    # Save the constraint index and the decision to allow modifications
    return SingleDecisionConstraintBridge{T,S}(constraint, f)
end

function MOI.supports_constraint(::Type{<:SingleDecisionConstraintBridge{T}},
                                 ::Type{SingleDecision},
                                 ::Type{<:MOI.AbstractScalarSet}) where T
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:SingleDecisionConstraintBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(::Type{<:SingleDecisionConstraintBridge{T,S}}) where {T,S}
    return [(MOI.SingleVariable, S), (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T})]
end
function MOIB.Constraint.concrete_bridge_type(::Type{<:SingleDecisionConstraintBridge{T}},
                              ::Type{SingleDecision},
                              S::Type{<:MOI.AbstractScalarSet}) where T
    return SingleDecisionConstraintBridge{T,S}
end

MOI.get(b::SingleDecisionConstraintBridge{T,S}, ::MOI.NumberOfConstraints{MOI.SingleVariable, S}) where {T,S} = 1
MOI.get(b::SingleDecisionConstraintBridge{T}, ::MOI.ListOfConstraintIndices{MOI.SingleVariable, S}) where {T,S} = [b.constraint]

function MOI.get(model::MOI.ModelLike, attr::MOI.AbstractConstraintAttribute,
                 bridge::SingleDecisionConstraintBridge)
    return MOI.get(model, attr, bridge.constraint)
end

function MOI.get(model::MOI.ModelLike, ::DecisionIndex,
                 bridge::SingleDecisionConstraintBridge{T}) where T
    return bridge.constraint
end

function MOI.delete(model::MOI.ModelLike, bridge::SingleDecisionConstraintBridge)
    if bridge.constraint.value != 0
        MOI.delete(model, bridge.constraint)
    end
    return nothing
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintSet,
                 bridge::SingleDecisionConstraintBridge{T,S}, change::S) where {T, S <: MOI.AbstractScalarSet}
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint, change)
    return nothing
end


text:
using FactCheck
using LLVM

facts("Context") do
    context("simple constructor test") do
        ctx = LLVM.create_ctx()
        @fact isa(ctx, LLVM.Context) => true
        LLVM.dispose!(ctx)
        @fact LLVM.isnull(ctx.handle) => true
    end

    context("simple global context constructor test") do
        ctx1 = LLVM.global_ctx()
        ctx2 = LLVM.global_ctx()
        @fact ctx1 => ctx2
    end

    context("simple with context constructor test") do
        LLVM.with_ctx() do ctx
            @fact isa(ctx, LLVM.Context) => true
        end
    end

    context("simple with global context constructor test") do
        LLVM.with_global_ctx() do ctx
            @fact isa(ctx, LLVM.Context) => true
            @fact LLVM.isglobalctx(ctx) => true
        end
    end
end


text:
@testset "symbolic" begin

# @testset "Dolang.eq_expr" begin
#     ex1 = :(z[0] = x + y(1))
#     ex2 = :(z[0] == x + y(1))

#     for ex in (ex1, ex2)
#         @test Dolang.stringify(ex) == :(_x_ + _y__1_ - _z__0_)
#         @test Dolang.stringify(ex, targets=[:_z__0_]) == :(_z__0_ = _x_ + _y__1_)
#         @test Dolang.stringify(ex, targets=[(:z, 0)]) == :(_z__0_ = _x_ + _y__1_)
#     end

# end

@testset "Dolang.stringify" begin
    @testset "Dolang.stringify(::Union{Symbol,String}, Integer)" begin
        @test Dolang.stringify(:x, 0) == :x__0_
        @test Dolang.stringify(:x, 1) == :x__1_
        @test Dolang.stringify(:x, -1) == :x_m1_
        @test Dolang.stringify(:x, -100) == :x_m100_

        # @test Dolang.stringify("x", 0) == :x__0_
        # @test Dolang.stringify("x", 1) == :x__1_
        # @test Dolang.stringify("x", -1) == :x_m1_
        # @test Dolang.stringify("x", -100) == :x_m100_

        @test Dolang.stringify((:x, 0)) == :x__0_
        @test Dolang.stringify((:x, 1)) == :x__1_
        @test Dolang.stringify((:x, -1)) == :x_m1_
        @test Dolang.stringify((:x, -100)) == :x_m100_
    end

    @testset "numbers" begin
        for T in (Float16, Float32, Float64, Int8, Int16, Int32, Int64)
            x = rand(T)
            @test Dolang.stringify(x) == x
        end
    end

    @testset "symbols" begin
        for i=1:10
            s = gensym()
            want = Symbol(s, "_")
            @test Dolang.stringify(s) == want
        end
    end

    @testset "x_(shift_Integer)" begin
        # for i=1:10, T in (Int8, Int16, Int32, Int64)
        for i=1:10, T in (Int64,)
            @test Dolang.stringify(string("x[t+", T(i), "]")) == Symbol("x__$(i)_")
            @test Dolang.stringify(string("x[t", T(-i), "]")) == Symbol("x_m$(i)_")
        end
    end

    @testset "other function calls" begin
        @testset "one argument" begin
            @test Dolang.stringify("sin(x)") == :(sin(x_))
            @test Dolang.stringify("sin(x[t-1])") == :(sin(x_m1_))
            @test Dolang.stringify("foobar(x[t+2])") == :(foobar(x__2_))
        end

        @testset "two arguments" begin
            @test Dolang.stringify("dot(x, y[t+1])") == :(dot(x_, y__1_))
            @test Dolang.stringify("plot(x[t-1], y)") == :(plot(x_m1_, y_))
            @test Dolang.stringify("bingbong(x[t+2], y)") == :(bingbong(x__2_, y_))
        end

        @testset "more args" begin
            for i=3:10
                ex = Expr(:call, :my_func, [:(x[t+$j]) for j in 1:i]...)
                want = Expr(:call, :my_func, [Symbol("x__", j, "_") for j in 1:i]...)
                @test Dolang.stringify(ex) == want
            end
        end

        @testset "arithmetic" begin
            @test Dolang.stringify(:(a[t+1] + b + c[t+2] + d[t-1])) == :(a__1_ + b_ + c__2_ + d_m1_)
            @test Dolang.stringify(:(a[t+1] * b * c[t+2] * d[t-1])) == :(a__1_ * b_ * c__2_ * d_m1_)
            @test Dolang.stringify(:(a[t+1] - b - c[t+2] - d[t-1])) == :(((a__1_ - b_) - c__2_) - d_m1_)
            @test Dolang.stringify(:(a[t+ 1]^ b)) == :(a__1_ ^ b_)
        end

    end

    # @testset "Expr(:(=), ...)" begin
    #     @testset "without targets" begin
    #         @test Dolang.stringify(:(x = y)) == :(_y_ - _x_)
    #     end

    #     @testset "with targets" begin
    #         @test Dolang.stringify(:(x = log(y(-1))); targets=[:x]) == :(_x_ = log(_y_m1_))
    #         @test Dolang.stringify(:(x == log(y(-1))); targets=[:x]) == :(_x_ = log(_y_m1_))
    #         @test_throws Dolang.stringifyError Dolang.stringify(:(x = y); targets=[:y])
    #     end
    # end

    @testset "stringify(::Tuple{Symbol,Int})" begin
        @test Dolang.stringify((:x, 0)) == :x__0_
        @test Dolang.stringify((:x, 1)) == :x__1_
        @test Dolang.stringify((:x, -1)) == :x_m1_
        @test Dolang.stringify((:x, -100)) == :x_m100_
    end

end

@testset "Dolang.time_shift" begin
    defs = Dict(:a=>:(b[t-1]/c))
    defs_sanitized = Dict(:a=>:(b[t-1]/c[t]))
    funcs = [:foobar]
    for shift in [-1, 0, 1]
        have = Dolang.time_shift(:(a+b[t+1] + c), shift)
        @test have == :(a + $(Dolang.create_variable(:b,shift+1)) + c)

        have = Dolang.time_shift(:(a+b[t+1] + c[t]), shift)
        @test have == :(a + $(Dolang.create_variable(:b,shift+1)) + $(Dolang.create_variable(:c,shift)) )
    end
end

@testset "Dolang.steady_state" begin
    @test Dolang.steady_state(:(a+b[t+1] + c)) == :(a+b[t]+c)
end

@testset "Dolang.list_symbols" begin
    ex = :(a+b[t+1]+c)
    out = Dolang.list_symbols(ex)
    @test out.variables == [(:b, 1)]
    @test out.parameters == [:a, :c]
    @test out.variables == @inferred Dolang.list_variables(ex)
    @test out.parameters == @inferred Dolang.list_parameters(ex)

    ex = :(a+b[t+1]+c[t+0])
    out = Dolang.list_symbols(ex)
    @test out.variables == [(:b, 1), (:c, 0)]
    @test out.parameters == [:a]
    @test out.variables == @inferred Dolang.list_variables(ex)
    @test out.parameters == @inferred Dolang.list_parameters(ex)

end

@testset "arg_name, arg_time, arg_name_time, arg_names" begin
    for i in 1:5
        s = gensym()
        s_vec = [gensym() for xxxxxx in 1:3]
        s_dict = OrderedDict(:m => s_vec, :M => [s])
        @test Dolang.arg_name(s) == s
        @test Dolang.arg_time(s) == 0
        @test Dolang.arg_name_time(s) == (s, 0)
        @test Dolang.arg_names(s_vec) == s_vec
        @test Dolang.arg_names(s_dict) == vcat(s_vec, s)
        for t in 1:10
            @test Dolang.arg_name((s, t)) == s
            @test Dolang.arg_time((s, t)) == t
            @test Dolang.arg_name_time((s, t)) == (s, t)
            @test Dolang.arg_name(Expr(:call, s, t)) == s
            @test Dolang.arg_time(Expr(:call, s, t)) == t
            @test Dolang.arg_name_time(Expr(:call, s, t)) == (s, t)
        end
    end

    for f in (Dolang.arg_name, Dolang.arg_time, Dolang.arg_name_time, Dolang.arg_names)
        @test_throws MethodError f(1)        # Number
        @test_throws MethodError f([1])      # Array of number
    end

    @test Dolang.arg_name_time(:_x__1_) == (:x, 1)
    @test Dolang.arg_name_time(:_x_m100_) == (:x, -100)
    @test Dolang.arg_name_time(:_this_is_x_m100_) == (:this_is_x, -100)
end

end  # @testset "symbolic"

max_issues_repo_name: EconForge/Dolang


text:
push!(LOAD_PATH,SRC_PATH)

using InputParser
using BuildBasis
using BuildOverlap
using NuclearRepulsion
using KineticEnergy
using NuclearAttraction
using ElecElecRepulsion
using SCF


text:
function nuts(θ0, δ, L, ∇L, L∇L=θ->(L(θ),∇L(θ)), M=100, Madapt=M, verbose=true, Δ_max=1000)
  doc"""
    - θ0      : initial model parameter
    - δ       : desirable average accept rate
    - L       : likelihood function
    - ∇L      : gradient of L
    - L∇L     : function to compute both L and ∇L, ideally in one pass
    - M       : sample number
    - Madapt  : number of samples for step size adaptation
    - verbose : whether to show log
  """

  function leapfrog(θ, r, ϵ)
    doc"""
      - θ : model parameter
      - r : momentum variable
      - ϵ : leapfrog step size
    """
    r̃ = r + 0.5 * ϵ * ∇L(θ)
    θ̃ = θ + ϵ * r̃
    r̃ = r̃ + 0.5 * ϵ * ∇L(θ̃)
    return θ̃, r̃
  end

  function find_reasonable_ϵ(θ)
    ϵ, r = 1, randn(length(θ))
    θ′, r′ = leapfrog(θ, r, ϵ)

    # This trick prevents the log-joint or its graident from being infinte
    # Ref: code start from Line 111 in https://github.com/mfouesneau/NUTS/blob/master/nuts.py
    # QUES: will this lead to some bias of the sampler?
    (l,g) = L∇L(θ′)
    while isinf(l) || any(isinf.(g))
      ϵ = ϵ * 0.5
      θ′, r′ = leapfrog(θ, r, ϵ)
      (l,g) = L∇L(θ′)
    end

    a = 2 * (exp(L(θ′) - 0.5 * dot(r′, r′)) / exp(L(θ) - 0.5 * dot(r, r)) > 0.5) - 1
    while (exp(L(θ′) - 0.5 * dot(r′, r′)) / exp(L(θ) - 0.5 * dot(r, r)))^float(a) > 2^float(-a)
      ϵ = 2^float(a) * ϵ
      θ′, r′ = leapfrog(θ, r, ϵ)
    end
    return ϵ
  end

  function build_tree(θ, r, u, v, j, ϵ, θ0, r0)
    doc"""
      - θ   : model parameter
      - r   : momentum variable
      - u   : slice variable
      - v   : direction ∈ {-1, 1}
      - j   : depth of tree
      - ϵ   : leapfrog step size
      - θ0  : initial model parameter
      - r0  : initial mometum variable
    """
    if j == 0
      # Base case - take one leapfrog step in the direction v.
      θ′, r′ = leapfrog(θ, r, v * ϵ)
      # NOTE: this trick prevents the log-joint or its graident from being infinte
      (l,g) = L∇L(θ′)
      while isinf(l) || any(isinf.(g))
        ϵ = ϵ * 0.5
        θ′, r′ = leapfrog(θ, r, ϵ)
        (l,g) = L∇L(θ′)
      end

      n′ = u <= exp(L(θ′) - 0.5 * dot(r′, r′))
      s′ = u < exp(Δ_max + L(θ′) - 0.5 * dot(r′, r′))
      return θ′, r′, θ′, r′, θ′, n′, s′, min(1, exp(L(θ′) - 0.5 * dot(r′, r′) - L(θ0) + 0.5 * dot(r0, r0))), 1
    else
      # Recursion - build the left and right subtrees.
      θm, rm, θp, rp, θ′, n′, s′, α′, n′_α = build_tree(θ, r, u, v, j - 1, ϵ, θ0, r0)
      if s′ == 1
        if v == -1
          θm, rm, _, _, θ′′, n′′, s′′, α′′, n′′_α = build_tree(θm, rm, u, v, j - 1, ϵ, θ0, r0)
        else
          _, _, θp, rp, θ′′, n′′, s′′, α′′, n′′_α = build_tree(θp, rp, u, v, j - 1, ϵ, θ0, r0)
        end
        if rand() < n′′ / (n′ + n′′)
          θ′ = θ′′
        end
        α′ = α′ + α′′
        n′_α = n′_α + n′′_α
        s′ = s′′ & (dot(θp - θm, rm) >= 0) & (dot(θp - θm, rp) >= 0)
        n′ = n′ + n′′
      end
      return θm, rm, θp, rp, θ′, n′, s′, α′, n′_α
    end
  end

  θs = Array{Array}(M + 1)

  θs[1], ϵ = θ0, find_reasonable_ϵ(θ0)
  μ, γ, t_0, κ = log(10 * ϵ), 0.05, 10, 0.75
  ϵ̄, H̄ = 1, 0

  if verbose println("[NUTS] start sampling for $M samples with inital ϵ=$ϵ") end

  for m = 1:M
    if verbose print('.') end
    r0 = randn(length(θ0))
    u = rand() * exp(L(θs[m]) - 0.5 * dot(r0, r0)) # Note: θ^{m-1} in the paper corresponds to
                                                   #       `θs[m]` in the code
    θm, θp, rm, rp, j, θs[m + 1], n, s = θs[m], θs[m], r0, r0, 0, θs[m], 1, 1
    α, n_α = NaN, NaN
    while s == 1
      v = rand([-1, 1])
      if v == -1
        θm, rm, _, _, θ′, n′, s′, α, n_α = build_tree(θm, rm, u, v, j, ϵ, θs[m], r0)
      else
        _, _, θp, rp, θ′, n′, s′, α, n_α = build_tree(θp, rp, u, v, j, ϵ, θs[m], r0)
      end
      if s′ == 1
        if rand() < min(1, n′ / n)
          θs[m + 1] = θ′
        end
      end
      n = n + n′
      s = s′ & (dot(θp - θm, rm) >= 0) & (dot(θp - θm, rp) >= 0)
      j = j + 1
    end
    if m + 1 <= Madapt + 1
      # NOTE: H̄ goes to negative when δ - α / n_α < 0
      H̄ = (1 - 1 / (m + t_0)) * H̄ + 1 / (m + t_0) * (δ - α / n_α)
      ϵ = exp(μ - sqrt(m) / γ * H̄)
      ϵ̄ = exp(m^float(-κ) * log(ϵ) + (1 - m^float(-κ)) * log(ϵ̄))
    else
      ϵ = ϵ̄
    end
  end

  if verbose println() end
  if verbose println("[NUTS] sampling complete with final apated ϵ = $ϵ") end

  return θs
end


text:
module PolynomialBasesTests

using Test
using Gridap

p1 = Point(1,2)
p2 = Point(0,3)
p = [p1,p2]

T = VectorValue{3,Float64}
orders = (1,1)
basis = MonomialBasis(T,orders)

v = evaluate(basis,p)
g = evaluate(∇(basis),p)

test_basis(basis,p,v,g)

T = Float64
orders = (1,1)
basis = MonomialBasis(T,orders)

v = evaluate(basis,p)
g = evaluate(∇(basis),p)

test_basis(basis,p,v,g)

T = VectorValue{2,Float64}
filter(e,order) = sum(e) <= order
order = 1
basis = MonomialBasis(2,T,filter,order)

@test length(basis) == 6

T = VectorValue{2,Float64}
basis = GradMonomialBasis(T,2)
v = evaluate(basis,p)
g = evaluate(∇(basis),p)

@test length(basis) == 12

p = Point{3,Int}[(1,2,3),(3,2,4)]

T = VectorValue{3,Float64}
order = 1
basis = CurlGradMonomialBasis(T,order)
v = evaluate(basis,p)
g = evaluate(∇(basis),p)

@test length(basis) == 6

end # module


text:
using QuPhys
using Test

using LinearAlgebra
using SparseArrays

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

@testset "Eigenvalues" begin
    A = sprand(ComplexF64, 15, 15, 0.1)
    A *= A'
    eigs1 = eigen(collect(A)).values
    eigs2 = eigensystem(A, k = 8, sigma = -5, krylovdim = 10)[1]
    eigs3 = eigensystem(A, k = 8, krylovdim = 10)[1]
    @test sum(abs.(eigs1[1:8] .- eigs2)) < 0.1
    @test sum(reverse(sort(abs.(eigs1))[end-7:end]) .- abs.(eigs3)) < 0.1
end

@testset "Time Evolution" begin
    N = 10

    a = kron(destroy(N), eye(2))
    a_d = a'
    sp = kron(eye(N), destroy(2))
    sm = sp'
    sx = sm + sp
    sy = 1im * (sm - sp)
    sz = sp * sm - sm * sp
    η = 0.01
    H = a_d * a + 0.5 * sz - 1im * η * (a - a_d) * sx
    psi0 = kron(fock(N, 0), fock(2, 0))
    t_l = LinRange(0, 1000, 1000)
    e_ops = [a_d * a]
    sol, expect_se = sesolve(H, psi0, t_l, e_ops = e_ops)
    @test sum(abs.(expect_se[1:end, 1] .- sin.(η * t_l).^2)) / length(t_l) < 0.1

    a = destroy(N)
    a_d = a'
    H = a_d * a
    c_ops = [sqrt(0.1) * a]
    e_ops = [a_d * a]
    psi0 = fock(N, 3)
    t_l = LinRange(0, 100, 1000)
    sol_me, expect_me = mesolve(H, psi0, t_l, c_ops, e_ops = e_ops);
    sol_mc, expect_mc = mcsolve(H, psi0, t_l, c_ops, n_traj = 500, e_ops = e_ops, ensemble_method = EnsembleSerial());
    @test sum(abs.(expect_mc .- expect_me)) / length(t_l) < 0.1

    sp1 = kron(destroy(2), eye(2))
    sm1 = sp1'
    sx1 = sm1 + sp1
    sy1 = 1im * (sm1 - sp1)
    sz1 = sp1 * sm1 - sm1 * sp1
    sp2 = kron(eye(2), destroy(2))
    sm2 = sp2'
    sx2 = sm2 + sp2
    sy2 = 1im * (sm2 - sp2)
    sz2 = sp2 * sm2 - sm2 * sp2
    ωq1, ωq2 = 1, 1
    γ1, γ2 = 0.05, 0.1
    H = 0.5 * ωq1 * sz1 + 0.5 * ωq2 * sz2
    c_ops = [sqrt(γ1) * sm1, sqrt(γ2) * sm2]
    psi0_1 = normalize(fock(2, 0) + fock(2, 1))
    psi0_2 = normalize(fock(2, 0) + fock(2, 1))
    psi0 = kron(psi0_1, psi0_2)
    t_l = LinRange(0, 10 / (γ1 + γ2), 1000)
    sol_me, expect_me = mesolve(H, psi0, t_l, c_ops, e_ops = [sp1 * sm1, sp2 * sm2]);
    sol_mc, expect_mc = mcsolve(H, psi0, t_l, c_ops, n_traj = 500, e_ops = [sp1 * sm1, sp2 * sm2], ensemble_method = EnsembleSerial());
    @test sum(abs.(expect_mc[1:end, 1:2] .- expect_me[1:end, 1:2])) / length(t_l) < 0.1
end

@testset "Entanglement" begin
    g = fock(2, 1)
    e = fock(2, 0)
    state = normalize(kron(g, e) + kron(e, g))
    rho = state * state'
    @test entanglement(state, [1], (2, 2)) / log(2) ≈ 1
end