text:
using Distributed
@everywhere include("./elastic_gd.jl")
@everywhere begin
	using LinearAlgebra
	using Random, Distributions
  using .ElasticGD
end

using DelimitedFiles, Formatting

@everywhere function make_data_oa(seed, rho1, rho2, p, sigma, n, beta_var)
	Random.seed!(seed)

	Ss11=(1-rho1)*I+rho1*ones(p,p)
	Ss12=rho2*ones(p,p)
	Ss=vcat(hcat(Ss11,Ss12),hcat(Ss12',Ss11))

  beta_star=vcat(rand(Normal(2,beta_var),p),zeros(p))

	X=rand(MvNormal(Ss),n)'
	y=X*beta_star+rand(Normal(0,sigma),n)
  X_train, y_train, X_val, y_val, X_test, y_test = train_val_test(X,y,splits=[0.6, 0.2, 0.2])

	return X_train, y_train, X_val, y_val, X_test, y_test, beta_star
end

function print_log(LOG_NAME, data_mat)
  n_obs=size(data_mat)[1]

	sens_mse_eg = data_mat[:,1]
  sens_mse_en = data_mat[:,2]
  spec_mse_eg = data_mat[:,3]
  spec_mse_en = data_mat[:,4]
  sens_cv_eg =  data_mat[:,5]
  sens_cv_en =  data_mat[:,6]
  spec_cv_eg =  data_mat[:,7]
  spec_cv_en =  data_mat[:,8]
  path_frac_eg = data_mat[:,9]
  path_frac_en = data_mat[:,10]
  mse_mse_eg =  data_mat[:,11]
  mse_mse_en =  data_mat[:,12]
  mse_cv_eg =  data_mat[:,13]
  mse_cv_en =  data_mat[:,14]
  alpha_eg = data_mat[:,15]
  alpha_en = data_mat[:,16]

  sens_mse_diff=sens_mse_eg.-sens_mse_en
  spec_mse_diff=spec_mse_eg.-spec_mse_en
  mse_mse_diff= mse_mse_eg.- mse_mse_en
  sens_cv_diff= sens_cv_eg.- sens_cv_en
  spec_cv_diff= spec_cv_eg.- spec_cv_en
  mse_cv_diff=  mse_cv_eg.-  mse_cv_en
  path_frac_diff=path_frac_eg.-path_frac_en

  sens_mc_eg =  data_mat[:,17]
  sens_mc_en =  data_mat[:,18]
  spec_mc_eg =  data_mat[:,19]
  spec_mc_en =  data_mat[:,20]
  sens_mc1_eg = data_mat[:,21]
  sens_mc1_en = data_mat[:,22]
  spec_mc1_eg = data_mat[:,23]
  spec_mc1_en = data_mat[:,24]

  sens_mc_diff= sens_mc_eg.-sens_mc_en
  spec_mc_diff= spec_mc_eg.-spec_mc_en
  sens_mc1_diff=sens_mc_eg.-sens_mc_en
  spec_mc1_diff=spec_mc_eg.-spec_mc_en

	data_mat_save = hcat(mse_mse_eg, mse_mse_en, mse_cv_eg, mse_cv_en, path_frac_eg, path_frac_en, alpha_eg, alpha_en)
	writedlm(string("data/",LOG_NAME,".txt"), data_mat_save, ' ')
	println1(LOG_NAME, format("Observations: {1:d}", n_obs))
	println1(LOG_NAME, "MSE")
	println1(LOG_NAME, format("eg. sens: {1:.3f} ({2:.3f}). spec: {3:.3f} ({4:.3f}). mse: {5:.3f} ({6:.3f}).", mean(sens_mse_eg), 3/sqrt(n_obs)*std(sens_mse_eg), mean(spec_mse_eg), 3/sqrt(n_obs)*std(spec_mse_eg), mean(mse_mse_eg), 3/sqrt(n_obs)*std(mse_mse_eg)))
	println1(LOG_NAME, format("en. sens: {1:.3f} ({2:.3f}). spec: {3:.3f} ({4:.3f}). mse: {5:.3f} ({6:.3f}).", mean(sens_mse_en), 3/sqrt(n_obs)*std(sens_mse_en), mean(spec_mse_en), 3/sqrt(n_obs)*std(spec_mse_en), mean(mse_mse_en), 3/sqrt(n_obs)*std(mse_mse_en)))
	println1(LOG_NAME, format("diff. sens: {1:.3f} ({2:.3f}). spec: {3:.3f} ({4:.3f}). mse: {5:.3f} ({6:.3f}).", mean(sens_mse_diff), 3/sqrt(n_obs)*std(sens_mse_diff), mean(spec_mse_diff), 3/sqrt(n_obs)*std(spec_mse_diff), mean(mse_mse_diff), 3/sqrt(n_obs)*std(mse_mse_diff)))
	println1(LOG_NAME, "CV")
	println1(LOG_NAME, format("eg. sens: {1:.3f} ({2:.3f}). spec: {3:.3f} ({4:.3f}). mse: {5:.3f} ({6:.3f}).", mean(sens_cv_eg), 3/sqrt(n_obs)*std(sens_cv_eg), mean(spec_cv_eg), 3/sqrt(n_obs)*std(spec_cv_eg), mean(mse_cv_eg), 3/sqrt(n_obs)*std(mse_cv_eg)))
	println1(LOG_NAME, format("en. sens: {1:.3f} ({2:.3f}). spec: {3:.3f} ({4:.3f}). mse: {5:.3f} ({6:.3f}).", mean(sens_cv_en), 3/sqrt(n_obs)*std(sens_cv_en), mean(spec_cv_en), 3/sqrt(n_obs)*std(spec_cv_en), mean(mse_cv_en), 3/sqrt(n_obs)*std(mse_cv_en)))
	println1(LOG_NAME, format("diff. sens: {1:.3f} ({2:.3f}). spec: {3:.3f} ({4:.3f}). mse: {5:.3f} ({6:.3f}).", mean(sens_cv_diff), 3/sqrt(n_obs)*std(sens_cv_diff), mean(spec_cv_diff), 3/sqrt(n_obs)*std(spec_cv_diff), mean(mse_cv_diff), 3/sqrt(n_obs)*std(mse_cv_diff)))

	println1(LOG_NAME, "PATH FRAC")
	println1(LOG_NAME, format("eg: {1:.3f} ({2:.3f}).", mean(path_frac_eg), 3/sqrt(n_obs)*std(path_frac_eg)))
	println1(LOG_NAME, format("en: {1:.3f} ({2:.3f}).", mean(path_frac_en), 3/sqrt(n_obs)*std(path_frac_en)))
	println1(LOG_NAME, format("diff: {1:.3f} ({2:.3f}).", mean(path_frac_diff), 3/sqrt(n_obs)*std(path_frac_diff)))

	println1(LOG_NAME, "MC")
	println1(LOG_NAME, format("eg. sens: {1:.3f} ({2:.3f}). spec: {3:.3f} ({4:.3f}).", mean(sens_mc_eg), 3/sqrt(n_obs)*std(sens_mc_eg), mean(spec_mc_eg), 3/sqrt(n_obs)*std(spec_mc_eg)))
	println1(LOG_NAME, format("en. sens: {1:.3f} ({2:.3f}). spec: {3:.3f} ({4:.3f}).", mean(sens_mc_en), 3/sqrt(n_obs)*std(sens_mc_en), mean(spec_mc_en), 3/sqrt(n_obs)*std(spec_mc_en)))
	println1(LOG_NAME, format("diff. sens: {1:.3f} ({2:.3f}). spec: {3:.3f} ({4:.3f}).", mean(sens_mc_diff), 3/sqrt(n_obs)*std(sens_mc_diff), mean(spec_mc_diff), 3/sqrt(n_obs)*std(spec_mc_diff)))

	println1(LOG_NAME, "MC1")
	println1(LOG_NAME, format("eg. sens: {1:.3f} ({2:.3f}). spec: {3:.3f} ({4:.3f}).", mean(sens_mc1_eg), 3/sqrt(n_obs)*std(sens_mc1_eg), mean(spec_mc1_eg), 3/sqrt(n_obs)*std(spec_mc1_eg)))
	println1(LOG_NAME, format("en. sens: {1:.3f} ({2:.3f}). spec: {3:.3f} ({4:.3f}).", mean(sens_mc1_en), 3/sqrt(n_obs)*std(sens_mc1_en), mean(spec_mc1_en), 3/sqrt(n_obs)*std(spec_mc1_en)))
	println1(LOG_NAME, format("diff. sens: {1:.3f} ({2:.3f}). spec: {3:.3f} ({4:.3f}).", mean(sens_mc1_diff), 3/sqrt(n_obs)*std(sens_mc1_diff), mean(spec_mc1_diff), 3/sqrt(n_obs)*std(spec_mc1_diff)))

	println1(LOG_NAME, "")
end

@everywhere function eg_en_iter(seed, rho1, rho2, p, sigma, n, beta_var, alg, n_folds)
	X_train, y_train, X_val, y_val, X_test, y_test, beta_star = make_data_oa(seed, rho1, rho2, p, sigma, n, beta_var)
  beta_mse_eg, beta_mse_en, beta_cv_eg, beta_cv_en, beta_path_eg, beta_path_en, alpha_eg, alpha_en, beta_mc_eg, beta_mc_en, beta_mc1_eg, beta_mc1_en = sweep_alpha(X_train, y_train, X_val, y_val, alg; n_folds=n_folds)
	sens_mse_eg, sens_mse_en, spec_mse_eg, spec_mse_en = get_sens_spec(beta_mse_eg, beta_mse_en, beta_star)
	sens_cv_eg,  sens_cv_en,  spec_cv_eg,  spec_cv_en  = get_sens_spec(beta_cv_eg,  beta_cv_en,  beta_star)
  path_frac_eg = get_path_frac(beta_path_eg, beta_star)
  path_frac_en = get_path_frac(beta_path_en, beta_star)
  mse_mse_eg=get_mse(y_test, X_test*beta_mse_eg)
  mse_mse_en=get_mse(y_test, X_test*beta_mse_en)
  mse_cv_eg=get_mse(y_test, X_test*beta_cv_eg)
  mse_cv_en=get_mse(y_test, X_test*beta_cv_en)

	sens_mc_eg,  sens_mc_en,  spec_mc_eg,  spec_mc_en  = get_sens_spec(beta_mc_eg,  beta_mc_en,  beta_star)
	sens_mc1_eg, sens_mc1_en, spec_mc1_eg, spec_mc1_en = get_sens_spec(beta_mc1_eg, beta_mc1_en, beta_star)

	return sens_mse_eg, sens_mse_en, spec_mse_eg, spec_mse_en, sens_cv_eg, sens_cv_en, spec_cv_eg, spec_cv_en, path_frac_eg, path_frac_en, mse_mse_eg, mse_mse_en, mse_cv_eg, mse_cv_en, alpha_eg, alpha_en, sens_mc_eg, sens_mc_en, spec_mc_eg, spec_mc_en, sens_mc1_eg, sens_mc1_en, spec_mc1_eg, spec_mc1_en
end

N_OBS=1000
RHO1=0.8
RHO2=0.5
p=10
sigma=10
n=100
BETA_VAR=1
N_FOLDS=10
GPU=""
ALG="norm_ns" #unnorm, norm_ns, norm_s, gs_ns, gs_s
for arg in ARGS
  eval(Meta.parse(arg))
end

LOG_NAME=string("run_synth_",RHO1,"_",RHO2,"_",ALG)
LOG_NAME=string(LOG_NAME,GPU)

nw=nworkers()
n_loops=N_OBS ÷ nw
seed_ranges=Any[]
for i_loop in 1:n_loops
	push!(seed_ranges,(i_loop-1)*nw+1:i_loop*nw)
end
n_last_loop=N_OBS- n_loops*nw
if n_last_loop>0
  push!(seed_ranges,n_loops*nw+1:n_loops*nw+n_last_loop)
end

eg_en_iter(1, RHO1, RHO2, p, sigma, n, BETA_VAR, ALG, N_FOLDS)
data_mat = zeros(0,24)

for seeds in seed_ranges
	data_mat_loc = pmap(seed -> eg_en_iter(seed, RHO1, RHO2, p, sigma, n, BETA_VAR, ALG, N_FOLDS), seeds)
  global data_mat=vcat(data_mat,hcat(map(collect, data_mat_loc)...)')

  print_log(LOG_NAME, data_mat)
end


text:
using StructJuMP

m = StructuredModel()

n = 10
@variable(m, x[1:n] >= 0)
@variable(m, y[1:n] >= 0)
for i in 1:n
    @constraint(m, x[i] - y[i] <= i)
end
@objective(m, :Min, x[1] - y[3] + 2x[6] + x[1]*x[2] - 2y[4]*x[3])

numScen = 8

for i in 1:numScen
    bl = StructuredModel(parent=m,id=i)
    @variable(bl, w[1:n] >= 0)
    @constraint(bl, sum(i*w[i] for i=2:2:n) == 1)
    @constraint(bl, sum(x[i] - 2w[n-i+1] for i=1:n) <= 0)
    @objective(bl, :Min, w[5] - w[1]*w[1] + 2w[2]*w[3])
end

ext: jl


text:
module Metrics

using ..Trajectories
using ..Trajectories:genrate_v!
using ..Force

export WorkHeat, get_workheat

mutable struct WorkHeat
    work::Real
    heat::Real
end

function get_workheat(tj::Trajectory, workrate::Function, potential::Function, step::Real)::WorkHeat
    if isnothing(tj.v)
        genrate_v!(tj)
    end
    work = sum(workrate.(tj.x[2:end-1], tj.v, tj.t[2:end-1])) * step
    dp = potential(tj.x[end-1], tj.t[end-1]) - potential(tj.x[2], tj.t[2])
    heat = work - dp
    WorkHeat(work, heat)
end

function get_workheat(tj::Trajectory, f::AbstractForce, step::Real)::WorkHeat
    get_workheat(tj,
                 (x, v, t) -> virtual_workrate(f, x, v, t),
                 (x, t) -> potential(f, x, t),
                 step
                )
end

end


text:
include("src/integration_helpers.jl")

using Test

# @testset ".jl" begin
    x0 = 2 * rand(10) .- 1
    x = 2 * rand(10) .- 1
    ϵ = rand(10)
    for _x0 in x0
        for _x in x
            for _ϵ in ϵ
                result = _x0 - _ϵ < _x < _x0 + _ϵ
                 @show is_in_Ball(_x, _x0, _ϵ) == result
            end
        end
    end
# end


text:
function perform_linprog(A,b,c,maxit=-1,tol=1e-10)
    """
        LINPROG  This code uses the revised simplex method to solve the linear
               programming problem: Minimize the cost c'x subject to
               equations Ax=b and nonnegativity x >=0:

                        f =  Min { c'x; Ax = b, x >= 0 },

               You must define an m by n matrix A , a column vector b with
               m components, and a column vector c with n components.

               The output vector x gives the minimum cost, which is the output f.

                        [x,f,itn] = linprog(A,b,c,k,maxit,tol)

               At most "maxit" iterations (default 10*length(b)) are applied
               and the actual number of iterations is returned in "itn".

               If the optimal solution is unbounded or the constraints are
               inconsistent then a diagnostic is displayed.
               Bland's rule is used to resolve degeneracies.  In exact
               arithmetic cycling is not possible.  But in real life!!!
               If x has more than 20 components it is returned in the
               sparse format. Note that if A has many zeros it is worth
               passing it to linprog in sparse format.

               Although written for teaching purposes this routine has
               successfully solved some problems with size(A) = [50,100000]!

               Please report any difficulties to: idc@math.canterbury.ac.nz

               New version                  (c) I.D.Coope, 1988, 1993
    """

    c = vec(c)
    b = vec(b)
    m,n = size(A)
    if maxit == -1
        maxit = 10*m
    end

    it=0

    D= sign(sign(b)+.5)
    D = diagm(D)                       # initial (inverse) basis matrix
    A = [A D]                      # incorporate slack/artificial variables
    B = collect(n+1:n+m)                     # initial basis
    N = collect(1:n)                    # non-basis

    phase = 1
    xb = abs(b)
    s =[zeros(n); ones(m)]    # supercost

    while phase < 3
        df = -1
        t = Inf
        yb= D'*s[B]        # multipliers for Ax=b
        while (it < maxit)
            if isempty(N)
                break
            end                              # no freedom for minimization
            r = s[N] - A[:,N]'*yb          # reduced costs
            rmin,q = findmin(r)                                    # determine new basic variable
            if rmin >= -tol*(norm(s[N],Inf)+1)
                break # optimal!
            end
            it = it+1
            if df >= 0                         # apply Bland's rule to avoid cycling
                J = find(r.<0)
                Nq = minimum(N[J])
                q = find(N.==Nq)[1]
            end
            d = vec(D*A[:,N[q]])
            I = find(d .> tol)
            if isempty(I)
                print("Solution is unbounded")
                it = -it
                break
            end
            xbd=xb[I]./d[I]
            r,p = findmin(xbd)
            p = I[p]
            if df >= 0                         # apply Bland's rule to avoid cycling
                J = find(xbd .== r)
                Bp = minimum(B[I[J]])
                p = find(B .== Bp)[1]
            end
            xb= xb - r*d

            xb[p] = r                            # update x
            df=r*rmin;
            # change in f
            v = vec(D[p,:])./d[p] # row vector
            yb= yb + v.*(s[N[q]] - dot(d',s[B]))
            d[p] = d[p] - 1
            D = D - d*v'                 # update inverse basis matrix
            t = B[p]
            B[p] = N[q]

            if t > n
                deleteat!(N,q)
            else
                N[q] = t
            end

        end
        xb=xb+D*(b-A[:,B]*xb);  # iterative refinement
        I = find(xb .< 0)                  # must be due to rounding error
        if length(I) > 0
            xb[I]=xb[I]-xb[I]                      # so correct
        end

        if phase == 2 || it < 0
            break                                # B, xb,n,m,res=A(:,B)*xb-b
        end
        if sum(xb.*s[B]) > tol
            it=-it
            print("No feasible solution")
            break
        end
        phase=phase+1                            # re-initialise for Phase 2
        s=1e6*norm(c, Inf)*s
        s[1:n]=c
    end
    x = zeros(2*n)
    x[B]=xb
    x=x[1:n]
    f=c'*x
    if it >= maxit
        print("Too many iterations")
        it=-it
    end
    return x
end


text:
# https://arxiv.org/abs/1309.3458

typealias Verts Array{Float64, 2}
typealias Cells SparseMatrixCSC{Int8, Int}
typealias Cell SparseVector{Int8, Int}
using IntervalTrees

function normalizer(V::Verts)
    d = size(V, 2)
    upper = mapslices(x->max(x...), V, 1)
    lower = mapslices(x->min(x...), V, 1)
    diff = upper-lower

    T = eye(d+1)
    T[d+1, 1:d] = -lower

    S = eye(d+1)
    S[1:d, 1:d] = diagm(vec(map(inv, diff)))

    mat = T*S
    mat, inv(mat)
end
function submanifold_mapping(V, face)
    p1, p2, p3 = map(i->V[face.nzind[i], :], 1:3)
    u1 = p2-p1
    u2 = p3-p1
    u3 = cross(u1, u2)
    T = eye(4)
    T[4, 1:3] = -p1
    M = eye(4)
    M[1:3, 1:3] = [u1 u2 u3]
    return T*M
end
function spatial_index(V::Verts, CV::Cells)
    d = size(V,2)
    cell_num = size(CV, 1)
    IntervalsType = IntervalValue{Float64, Int64}
    boxes1D = Array{IntervalsType, 2}(0, d)
    for ci in 1:cell_num
        intervals = map((l,u)->IntervalsType(l,u,ci), bbox(V,CV[ci, :])...)
        boxes1D = vcat(boxes1D, intervals)
    end
    trees = mapslices(IntervalTree{Float64, IntervalsType}, sort(boxes1D, 1), 1)

    function intersect_intervals(intervals)
        cells = Array{Int64,1}[]
        for axis in 1:d
            vs = map(i->i.value, intersect(trees[axis], intervals[axis]))
            push!(cells, vs)
        end
        mapreduce(x->x, intersect, cells)
    end

    mapping = Array{Int64,1}[]
    for ci in 1:cell_num
        cell_indexes = setdiff(intersect_intervals(boxes1D[ci, :]), [ci])
        push!(mapping, cell_indexes)
    end

    mapping
end

max_stars_repo_name: furio/lar-parallel-tests

max_issues_repo_name: furio/lar-parallel-tests

max_forks_repo_name: furio/lar-parallel-tests


text:
function Netflix()::Persa.CFDataset
    filename = "$(defdir)/netflix/netflix.csv"
    file = readtable(filename, separator = ',', header = false)

    df = DataFrame()

    df[:user] = labelencode(labelmap(file[:,1]), file[:,1])
    df[:item] = file[:,2]
    df[:rating] = file[:,3]

    return Persa.Dataset(df)
end

function MovieTweeting()::Persa.TimeCFDataset
    filename = "$(defdir)/Movie-Tweeting-200k/ratings.dat"
    file = readtable(filename, separator = ':', header = false)

    df = DataFrame()

    df[:user] = file[:,1]
    df[:item] = labelencode(labelmap(file[:,2]), file[:,2])
    df[:rating] = file[:,3]
    df[:timestamp] = file[:,4]

    return Persa.Dataset(df)
end

function MovieTweeting10k()::Persa.TimeCFDataset
    filename = "$(defdir)/mt-snapshot-10k/ratings.dat"
    file = readtable(filename, separator = ':', header = false)

    df = DataFrame()

    df[:user] = file[:,1]
    df[:item] = labelencode(labelmap(file[:,2]), file[:,2])
    df[:rating] = file[:,3]
    df[:timestamp] = file[:,4]

    return Persa.Dataset(df)
end

function CiaoDVD()::Persa.TimeCFDataset
    filename = "$(defdir)/CiaoDVD/movie-ratings.txt"
    file = readtable(filename, separator = ',', header = false)

    df = DataFrame()

    df[:user] = file[:,1]
    df[:item] = file[:,2]
    df[:rating] = file[:,5]
    df[:timestamp] = convert(Array{Int}, Dates.datetime2unix.(Dates.DateTime(file[:,6])))

    return Persa.Dataset(df)
end

function FilmTrust()::Persa.CFDataset
    filename = "$(defdir)/FilmTrust/ratings.txt"
    file = readtable(filename, separator = ' ', header = false)

    df = DataFrame()

    df[:user] = file[:,1]
    df[:item] = file[:,2]
    df[:rating] = file[:,3]

    return Persa.Dataset(df)
end

function YahooMusic()::Persa.CFDataset
    filename = "$(defdir)/yahoo-music-r3/ymusic-r3-dummy-time.dat"
    file = readtable(filename, separator = ' ', header = false)

    df = DataFrame()

    df[:user] = file[:,1]
    df[:item] = file[:,2]
    df[:rating] = file[:,3]

    return Persa.Dataset(df)
end

function LastFM()::Persa.CFDataset
    filename = "$(defdir)/lastfm/last_fm.dat"
    file = readtable(filename, separator = ',', header = false)

    df = DataFrame()

    df[:user] = labelencode(labelmap(file[:,1]), file[:,1])
    df[:item] = labelencode(labelmap(file[:,2]), file[:,2])
    df[:rating] = file[:,3]

    return Persa.Dataset(df)
end


text:
export stage, cached_stage, compute, debug_compute, cached, free!
using Compat

compute(x) = compute(Context(), x)
compute(ctx, c::Chunk) = c

collect(ctx::Context, c) = collect(ctx, compute(ctx, c))
collect(d::Union{Chunk,Thunk}) = collect(Context(), d)

abstract type Computation end

compute(ctx, c::Computation) = compute(ctx, stage(ctx, c))
collect(c::Computation) = collect(Context(), c)

function finish_task!(state, node, node_order; free=true)
    if istask(node) && node.cache
        node.cache_ref = Nullable{Any}(state.cache[node])
    end
    immediate_next = false
    for dep in sort!(collect(state.dependents[node]), by=node_order)
        set = state.waiting[dep]
        pop!(set, node)
        if isempty(set)
            pop!(state.waiting, dep)
            push!(state.ready, dep)
            immediate_next = true
        end
        # todo: free data
    end
    for inp in inputs(node)
        if inp in keys(state.waiting_data)
            s = state.waiting_data[inp]
            if node in s
                pop!(s, node)
            end
            if free && isempty(s)
                if haskey(state.cache, inp)
                    _node = state.cache[inp]
                    free!(_node, force=false, cache=(istask(inp) && inp.cache))
                    pop!(state.cache, inp)
                end
            end
        end
    end
    push!(state.finished, node)
    pop!(state.running, node)
    immediate_next
end

###### Scheduler #######
"""
Compute a Thunk - creates the DAG, assigns ranks to
nodes for tie breaking and runs the scheduler.
"""
function compute(ctx, d::Thunk)
    master = OSProc(myid())
    @dbg timespan_start(ctx, :scheduler_init, 0, master)

    ps = procs(ctx)
    chan = Channel{Any}(32)
    deps = dependents(d)
    ord = order(d, noffspring(deps))

    node_order = x -> -get(ord, x, 0)
    state = start_state(deps, node_order)
    # start off some tasks
    for p in ps
        isempty(state.ready) && break
        task = pop_with_affinity!(ctx, state.ready, p, false)
        if task !== nothing
            fire_task!(ctx, task, p, state, chan, node_order)
        end
    end
    @dbg timespan_end(ctx, :scheduler_init, 0, master)

    while !isempty(state.ready) || !isempty(state.running)

        if isempty(state.running) && !isempty(state.ready)
            for p in ps
                isempty(state.ready) && break
                task = pop_with_affinity!(ctx, state.ready, p, false)
                if task !== nothing
                    fire_task!(ctx, task, p, state, chan, node_order)
                end
            end
        end

        if isempty(state.running)
            # the block above fired only meta tasks
            continue
        end

        proc, thunk_id, res = take!(chan)
        if isa(res, CapturedException) || isa(res, RemoteException)
            rethrow(res)
        end
        node = _thunk_dict[thunk_id]
        @logmsg("WORKER $(proc.pid) - $node ($(node.f)) input:$(node.inputs)")
        state.cache[node] = res

        @dbg timespan_start(ctx, :scheduler, thunk_id, master)
        immediate_next = finish_task!(state, node, node_order)
        if !isempty(state.ready)
            thunk = pop_with_affinity!(Context(ps), state.ready, proc, immediate_next)
            if thunk !== nothing
                fire_task!(ctx, thunk, proc, state, chan, node_order)
            end
        end
        @dbg timespan_end(ctx, :scheduler, thunk_id, master)
    end
    state.cache[d]
end

function pop_with_affinity!(ctx, tasks, proc, immediate_next)
    if immediate_next
        # fast path
        if proc in first.(affinity(tasks[end]))
            return pop!(tasks)
        end
    end
    parent_affinities = affinity.(tasks)
    for i=length(tasks):-1:1
        # TODO: use the size
        if proc in first.(parent_affinities[i])
            t = tasks[i]
            deleteat!(tasks, i)
            return t
        end
    end
    for i=length(tasks):-1:1
        # use up tasks without affinitites
        # let the procs with the respective affinities pick up
        # other tasks
        if isempty(parent_affinities[i])
            t = tasks[i]
            deleteat!(tasks, i)
            return t
        end
        aff = first.(parent_affinities[i])
        if all(!(p in aff) for p in procs(ctx)) # no proc is ever going to ask for it
            t = tasks[i]
            deleteat!(tasks, i)
            return t
        end
    end
    return nothing
end

function fire_task!(ctx, thunk, proc, state, chan, node_order)
    @logmsg("W$(proc.pid) + $thunk ($(showloc(thunk.f, length(thunk.inputs)))) input:$(thunk.inputs) cache:$(thunk.cache) $(thunk.cache_ref)")
    push!(state.running, thunk)
    if thunk.cache && !isnull(thunk.cache_ref)
        # the result might be already cached
        data = unrelease(get(thunk.cache_ref)) # ask worker to keep the data around
                                          # till this compute cycle frees it
        if !isnull(data)
            @logmsg("cache hit: $(get(thunk.cache_ref))")
            state.cache[thunk] = get(data)
            immediate_next = finish_task!(state, thunk, node_order; free=false)
            if !isempty(state.ready)
                thunk = pop_with_affinity!(ctx, state.ready, proc, immediate_next)
                if thunk !== nothing
                    fire_task!(ctx, thunk, proc, state, chan, node_order)
                end
            end
            return
        else
            thunk.cache_ref = Nullable{Any}()
            @logmsg("cache miss: $(thunk.cache_ref) recomputing $(thunk)")
        end
    end

    if thunk.meta
        # Run it on the parent node
        # do not _move data.
        p = OSProc(myid())
        @dbg timespan_start(ctx, :comm, thunk.id, p)
        fetched = map(thunk.inputs) do x
            istask(x) ? state.cache[x] : x
        end
        @dbg timespan_end(ctx, :comm, thunk.id, p)

        @dbg timespan_start(ctx, :compute, thunk.id, p)
        res = thunk.f(fetched...)
        @dbg timespan_end(ctx, :compute, thunk.id, p)

        #push!(state.running, thunk)
        state.cache[thunk] = res
        immediate_next = finish_task!(state, thunk, node_order; free=false)
        if !isempty(state.ready)
            if immediate_next
                thunk = pop!(state.ready)
            else
                thunk = pop_with_affinity!(ctx, state.ready, proc, immediate_next)
            end
            if thunk !== nothing
                fire_task!(ctx, thunk, proc, state, chan, node_order)
            end
        end
        return
    end

    data = map(thunk.inputs) do x
        istask(x) ? state.cache[x] : x
    end
    async_apply(ctx, proc, thunk.id, thunk.f, data, chan, thunk.get_result, thunk.persist, thunk.cache)
end

##### Scheduling logic #####

"find the set of direct dependents for each task"
function dependents(node::Thunk, deps=Dict{Thunk, Set{Thunk}}())
    if !haskey(deps, node)
        deps[node] = Set{Thunk}()
    end
    for inp = inputs(node)
        if isa(inp, Thunk)
            s::Set{Thunk} = Base.@get!(deps, inp, Set{Thunk}())
            push!(s, node)
            dependents(inp, deps)
        end
    end
    deps
end

"""
recursively find the number of taks dependent on each task in the DAG.
Input: dependents dict
"""
function noffspring(dpents::Dict)
    Dict(node => noffspring(node, dpents) for node in keys(dpents))
end

function noffspring(n, dpents)
    if haskey(dpents, n)
        ds = dpents[n]
        reduce(+, length(ds), noffspring(d, dpents) for d in ds)
    else
        0
    end
end

"""
Given a root node of the DAG, calculates a total order for tie-braking

  * Root node gets score 1,
  * rest of the nodes are explored in DFS fashion but chunks
    of each node are explored in order of `noffspring`,
    i.e. total number of tasks depending on the result of the said node.

Args:
    - node: root node
    - ndeps: result of `noffspring`
"""
function order(node::Thunk, ndeps)
    function recur(nodes, s)
        for n in nodes
            output[n] = s += 1
            parents = collect(Iterators.filter(istask, inputs(n)))
            s = recur(sort!(parents, by=k->get(ndeps,k,0)), s)
        end
        return s
    end
    output = Dict{Thunk,Int}()
    recur([node], 0)
    return output
end

const OneToMany = Dict{Thunk, Set{Thunk}}
struct ComputeState
    dependents::OneToMany
    finished::Set{Thunk}
    waiting::OneToMany
    waiting_data::OneToMany
    ready::Vector{Thunk}
    cache::Dict{Thunk, Any}
    running::Set{Thunk}
end

function start_state(deps::Dict, node_order)
    state = ComputeState(
                  deps,
                  Set{Thunk}(),
                  OneToMany(),
                  OneToMany(),
                  Vector{Thunk}(0),
                  Dict{Thunk, Any}(),
                  Set{Thunk}()
                 )

    nodes = sort(collect(keys(deps)), by=node_order)
    merge!(state.waiting_data, deps)
    for k in nodes
        if istask(k)
            waiting = Set{Thunk}(Iterators.filter(istask,
                                                  inputs(k)))
            if isempty(waiting)
                push!(state.ready, k)
            else
                state.waiting[k] = waiting
            end
        end
    end
    state
end

_move(ctx, to_proc, x) = x
_move(ctx, to_proc::OSProc, x::Union{Chunk, Thunk}) = collect(ctx, x)

function do_task(ctx, proc, thunk_id, f, data, send_result, persist, cache)
    @dbg timespan_start(ctx, :comm, thunk_id, proc)
    time_cost = @elapsed fetched = map(x->_move(ctx, proc, x), data)
    @dbg timespan_end(ctx, :comm, thunk_id, proc)

    @dbg timespan_start(ctx, :compute, thunk_id, proc)
    result_meta = try
        res = f(fetched...)
        (proc, thunk_id, send_result ? res : tochunk(res, persist=persist, cache=persist ? true : cache)) #todo: add more metadata
    catch ex
        bt = catch_backtrace()
        (proc, thunk_id, RemoteException(myid(), CapturedException(ex, bt)))
    end
    @dbg timespan_end(ctx, :compute, thunk_id, proc)
    result_meta
end

function async_apply(ctx, p::OSProc, thunk_id, f, data, chan, send_res, persist, cache)
    @schedule begin
        try
            put!(chan, Base.remotecall_fetch(do_task, p.pid, ctx, p, thunk_id, f, data, send_res, persist, cache))
        catch ex
            bt = catch_backtrace()
            put!(chan, (p, thunk_id, CapturedException(ex, bt)))
        end
        nothing
    end
end

function debug_compute(ctx::Context, args...; profile=false)
    @time res = compute(ctx, args...)
    get_logs!(ctx.log_sink), res
end

function debug_compute(arg; profile=false)
    ctx = Context()
    dbgctx = Context(procs(ctx), LocalEventLog(), profile)
    debug_compute(dbgctx, arg)
end

Base.@deprecate gather(ctx, x) collect(ctx, x)
Base.@deprecate gather(x) collect(x)


text:
function solve()
    s, w = [parse(Int, x) for x in split(readline())]
    if w  ≥ s
        "unsafe"
    else
        "safe"
    end
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
"""
    riemann_integration(f::Function, a::Real, b::Real, n::Int, approx::Symbol = :midpoint)

a Riemann sum is a certain kind of approximation of an integral by a finite sum.
The sum is calculated by partitioning the region into shapes (rectangles, trapezoids, parabolas, or cubics) that together form a region that is similar to the region being measured, then calculating the area for each of these shapes, and finally adding all of these small areas together.

Because the region filled by the small shapes is usually not exactly the same shape as the region being measured, the Riemann sum will differ from the area being measured.
This error can be reduced by dividing up the region more finely, using smaller and smaller shapes.
As the shapes get smaller and smaller, the sum approaches the Riemann integral.

# Arguments
- `f`: the function to integrate. (at the momment only single variable is suported)
- `a`: Start of the integration limit.
- `b`: End of the integration limit.
- `n`: Number of points to sample. (as n increase, error decrease)
- `approx`: Indicate the method of approximation (midpoint, left or right)

# Examples
```julia
julia> riemann_integration(x -> x, 1, 3, 1_000, :midpoint)  # 4.0
julia> riemann_integration(x -> x, 1, 3, 1_000, :left)      # 3.997997997997998
julia> riemann_integration(x -> x, 1, 3, 1_000, :right)     # 4.002002002002002
julia> riemann_integration(x -> 3*x^2, 0, 1, 100000)        # integrate a polynomial
0.9999999999750021
julia> riemann_integration(x -> sin(x), 0, pi, 1000)          # integrate the sin function
2.0000008241146774
```

# Refereces
- https://www.khanacademy.org/math/ap-calculus-ab/ab-integration-new/ab-6-2/a/riemann-sums-review
- https://math.libretexts.org/Courses/Mount_Royal_University/MATH_2200%3A_Calculus_for_Scientists_II/2%3A_Techniques_of_Integration/2.5%3A_Numerical_Integration_-_Midpoint%2C_Trapezoid%2C_Simpson's_rule
- https://abel.math.harvard.edu/~knill/teaching/math1a_2011/handouts/40-numerical.pdf
- https://en.wikipedia.org/wiki/Riemann_integral

Contributed By:- [AugustoCL](https://github.com/AugustoCL)
"""
function riemann_integration(f::Function, a::Real, b::Real, n::Int, approx::Symbol = :midpoint)
    # width of the rectangles
    Δₓ = (b - a) / n

    # methods of approximation (:midpoint, :left, :right)
    if approx == :midpoint
        sum_range = (a + Δₓ/2):Δₓ:(b - Δₓ/2)
    elseif approx == :left
        sum_range = a:Δₓ:(b - Δₓ)
    elseif approx == :right
        sum_range = (a + Δₓ):Δₓ:b
    else
        throw(ArgumentError("The symbol :$approx is not a valid argument. Insert :midpoint, :left or :right"))
    end

    # sum of the height of the rectangles
    Σ = 0.0
    for i in 2:length(sum_range)
        a  = sum_range[i-1]
        b  = sum_range[i]
        xᵢ = a + (b-a)*rand()   # draw a uniform(a,b) for each subinterval [a,b]
        Σ += f(xᵢ)
    end

    # approximate integral of f
    return Δₓ * Σ
end

max_stars_repo_name: frankschmitt/Julia

max_issues_repo_name: frankschmitt/Julia

max_forks_repo_name: frankschmitt/Julia


text:
using Requires, FillArrays

instead(::Operations, args, ::AType{Type{Fill}}, x, sz) =
  ([fill, args[2], args[3]], (Const(fill), x, sz))

nop() = return

xlogy(x, y) = x*log(y)

function requires()
  @require Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" begin
    instead(::Operations, args, ::AType{typeof(Zygote.accum_global)}, _...) =
      [nop], [Const(nop)]
    instead(::Operations, args, ::AType{typeof(Zygote.accum_param)}, _...) =
      [nop], [Const(nop)]
    instead(::Operations, args, ::AType{typeof(Zygote._push!)}, xs, x) =
      [push!, args[2], args[3]], [Const(push!), xs, x]
  end
  @require Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" begin
    if isdefined(Flux, :xlogy)
      @eval instead(::Operations, args, ::AType{typeof(Flux.xlogy)}, x, y) =
              [xlogy, args[2], args[3]], [Const(xlogy), x, y]
    end
    constructor(::Type{<:Flux.Chain}) = xs -> Flux.Chain(xs...)
  end
end


text:
#Table 6.2
#Run uility functions at bottom
#Change the setting of Lambda if necessary
using NMF, LinearAlgebra, MultiscaleGraphSignalTransforms, ScikitLearn
@sk_import linear_model:Lasso

########################### utility functions ###########################
function create_scaling_vectors(GProws::GraphPart)
    GP = GProws
    ### creating scaling vectors
    temp = GP.tag.== 0;
    scaling_index = findall(temp);# index of the scaling vectors .i.e. vectors with tag == 0

    num_scaling = length(scaling_index);
    scaling_vectors = fill(0., (size(GP.tag,1), num_scaling));

    for i = 1:num_scaling
        ind = scaling_index[i]
        dvec = fill(0., (size(GP.tag,1),1))
        dvec[ind[1],1] = 1.;
        BS = bs_level(GP, ind[2] - 1)
        scaling_vectors[:,i] = ghwt_synthesis(dvec, GP, BS)
    end
    return scaling_vectors
end

function create_PHI(scaling_vectors_rows, scaling_vectors_cols, m, n)
    PHI_right = repeat(scaling_vectors_cols[:,1:n], outer = (1, m))
    PHI_left = repeat(scaling_vectors_rows[:,1:m], inner = (1, n))
    return PHI_left, PHI_right
end

function ghwtinit(matrix, k)
    GProws, GPcols = partition_tree_matrixDhillon(matrix)
    dmatrix = ghwt_analysis_2d(matrix, GProws, GPcols)
    scaling_vectors_rows = create_scaling_vectors(GProws)
    scaling_vectors_cols = create_scaling_vectors(GPcols)

    PHI_left, PHI_right = create_PHI(scaling_vectors_rows, scaling_vectors_cols, 10, 10)
    scaling_vectors = fill(0., (length(matrix[:]), size(PHI_left,2)))
    for i = 1:size(PHI_left,2)
        scaling_vectors[:,i] = (PHI_left[:,i]*PHI_right[:,i]')[:]
    end
    ##########################
    ####################################### perform lasso

    f = matrix[:]
    #Lambda = range(0.00001,stop = 0.0001, length = 10)
    Lambda = range(0.001,stop = 0.01, length = 10)# m=125,n=25,k=5,sigma = 0.5,
    #Lambda = range(0.001, stop = 0.01, length = 10)
    B = fill(0., (size(scaling_vectors,2), length(Lambda))) ### initialize coefficients
    for i = 1:length(Lambda)
        lasso_out = Lasso(alpha = Lambda[i], positive = true, fit_intercept = false)
        ScikitLearn.fit!(lasso_out, scaling_vectors, f)
        B[:,i] = lasso_out.coef_
    end

    j = 1
    if sum((B[:,j].!=0)) < k
        return false, false
    else
        while j <10 && sum(B[:,j+1].!=0) >= k
            j += 1
        end
    end
    index = findall(B[:,j].!=0)[1:k]
    rescaleW = repeat(sqrt.(B[index,j]'),size(PHI_left,1),1)
    rescaleH = repeat(sqrt.(B[index,j]'),size(PHI_right,1),1)
    Wtilda = PHI_left[:,index].*rescaleW
    Htilda = Array{Float64,2}((PHI_right[:,index].*rescaleH)')
    return Wtilda, Htilda
end

function improved_NNDSVD(matrix, k)
    m,n = size(matrix)
    U,S,V = svd(matrix)
    W = zeros(m,k)
    H = zeros(k,n)
    Y = U*sqrt.(Diagonal(S))
    Z = sqrt.(Diagonal(S))*transpose(V)

    W[:,1] = abs.(Y[:,1])
    H[1,:] = abs.(Z[1,:])
    j = 2
    for i = 2:k
        if mod(i,2) == 0
            W[:,i] = Y[:,j]
            H[i,:] = Z[i,:]
        else
            W[:,i] = .-Y[:,j]
            H[i,:] = .-Z[i,:]
            j += 1
        end
        W[W[:,i] .< 0,i] .= 0
        H[i,H[i,:] .< 0] .= 0
    end
    return W, H
end
########################### end of utility functions ###########################

### Now the start of the actual NMF experiments.

k = 5
m = 125
n = 25
sigma = 0.5 #noise parameter
#methodid = 1 corresponds to ALSPGRAD.
#methodid = 2 corresponds to HALS.
methodid = 2

num_pos_eghwt_svd = 0
num_neg_eghwt_svd = 0
num_pos_eghwt_rand = 0
num_neg_eghwt_rand = 0
num_pos_svd_rand = 0
num_neg_svd_rand = 0
num_pos_svd_newsvd = 0
num_neg_svd_newsvd = 0
num_pos_eghwt_newsvd = 0
num_neg_eghwt_newsvd = 0
for i = 1:50
    Wtrue = rand(m,k)
    Htrue = rand(k,n)
    X = Wtrue*Htrue + sigma*rand(m,n)

    #X = rand(m,n)

    ################
    ######################################

    Wtilda, Htilda = ghwtinit(X, k)
    if Wtilda == false
        continue
    end
    Wsvd, Hsvd = NMF.nndsvd(X, k)
    Wrand, Hrand = NMF.randinit(X, k)
    Wnewsvd, Hnewsvd = improved_NNDSVD(X, k)

    iters = 1000000
    method = [NMF.ALSPGrad{Float64}(maxiter=iters, tolg=1.0e-6), NMF.CoordinateDescent{Float64}(maxiter=iters, α=0.5, l₁ratio=0.5)]

    a = NMF.solve!(method[methodid], X, copy(Wtilda), copy(Htilda))
    b = NMF.solve!(method[methodid], X, copy(Wsvd), copy(Hsvd))
    c = NMF.solve!(method[methodid], X, copy(Wrand), copy(Hrand))
    d = NMF.solve!(method[methodid], X, copy(Wnewsvd), copy(Hnewsvd))

    #println("initial norm ",norm(X - Wtilda*Htilda)," ",norm(X-Wsvd*Hsvd)," ",norm(X - Wrand*Hrand)," ", norm(X - Wnewsvd*Hnewsvd))
    #println(" ",a.niters," ",b.niters," ",c.niters," ",d.niters)
    #println("final norm",norm(X - a.W*a.H)," ",norm(X-b.W*b.H)," ",norm(X - c.W*c.H)," ", norm(X - d.W*d.H))
    if a.niters < b.niters
        global num_pos_eghwt_svd += 1
    else
        global num_neg_eghwt_svd += 1
    end
    if a.niters < c.niters
        global num_pos_eghwt_rand += 1
    else
        global num_neg_eghwt_rand += 1
    end
    if b.niters < c.niters
        global num_pos_svd_rand += 1
    else
        global num_neg_svd_rand += 1
    end
    if a.niters < d.niters
        global num_pos_eghwt_newsvd += 1
    else
        global num_neg_eghwt_newsvd += 1
    end
    if b.niters < d.niters
        global num_pos_svd_newsvd += 1
    else
        global num_neg_svd_newsvd += 1
    end
end

println("compare eGHWT with NNDSVD ", num_pos_eghwt_svd/(num_pos_eghwt_svd + num_neg_eghwt_svd))
println("compare eGHWT with Random ", num_pos_eghwt_rand/(num_pos_eghwt_rand + num_neg_eghwt_rand))
println("compare eGHWT with NNSVD-LRC ", num_pos_eghwt_newsvd/(num_pos_eghwt_newsvd + num_neg_eghwt_newsvd))


text:
"""
    A = JopPermutedims(sp, perm)

where `sp::JetSpace` is the domain of the operator, and `perm` is an array
specifying the permutation of array dimensions.  The range of the operator
is inferred from `sp` and `perm`.

# Example:
```julia
A = JopPermutedims(JetSpace(Float64,3,2),[2;1])
m = [1 2 3 ; 4 5 6]
d = A*m # d = [1 4 ; 2 5 ; 3 6]
```
"""
function JopPermutedims(dom, perm)
    shpdom = size(dom)
    shprng = ntuple(i->shpdom[perm[i]], length(shpdom))
    rng = JetSpace(eltype(dom), shprng)

    P = zeros(Int,length(perm),length(perm))
    for i = 1:length(perm)
        P[i,perm[i]] = 1
    end
    iperm = P'*[1:length(shpdom);]

    JopLn(dom = dom, rng = rng, df! =  JopPermutedims_df!, df′! = JopPermutedims_df′!, s = (perm=perm, iperm=iperm))
end

export JopPermutedims

JopPermutedims_df!(d, m; perm, kwargs...) = permutedims!(d, m, perm)
JopPermutedims_df′!(m, d; iperm, kwargs...) = permutedims!(m, d, iperm)


text:
using Burgers
using Test

@testset "Burgers" begin
    include("data.jl")
    include("deeponet.jl")
end


text:
## EXPOSE OPTIMISERS TO MLJ (for eg, tuning)

# Here we make the optimiser structs "transparent" so that their
# field values are exposed by calls to MLJ.params

for opt in (:Descent,
            :Momentum,
            :Nesterov,
            :RMSProp,
            :ADAM,
            :RADAM,
            :AdaMax,
            :OADAM,
            :ADAGrad,
            :ADADelta,
            :AMSGrad,
            :NADAM,
            :AdaBelief,
            :Optimiser,
            :InvDecay, :ExpDecay, :WeightDecay,
            :ClipValue,
            :ClipNorm) # last updated: Flux.jl 0.12.3

    @eval begin
        MLJModelInterface.istransparent(m::Flux.$opt) = true
    end
end

## GENERAL METHOD TO OPTIMIZE A CHAIN

struct Mover{R<:AbstractResource}
    resource::R
end

(::Mover{<:CPU1})(data) = Flux.cpu(data)
(::Mover{<:CUDALibs})(data) = Flux.gpu(data)

"""
Custom training loop. Here, `loss_func` is the objective
function to optimise, `parameters` are the model parameters,
`optimiser` is the optimizer to be used, `X` (input features)is a
vector of arrays where the last dimension is the batch size. `y`
is the target observation vector.
"""
function train!(loss_func, parameters, optimiser, X, y)
    for i=1:length(X)
        gs = Flux.gradient(parameters) do
            training_loss = loss_func(X[i], y[i])
            return training_loss
        end
        Flux.update!(optimiser, parameters, gs)
    end
end

"""
    fit!(chain,
         optimiser,
         loss,
         epochs,
         lambda,
         alpha,
         verbosity,
         acceleration,
         X,
         y)

Optimize a Flux model `chain` using the regularization parameters
`lambda` (strength) and `alpha` (l2/l1 mix), where `loss(yhat, y) ` is
the supervised loss for instances (or vectors of instances) of the
target predictions `yhat` and target observations `y`.

Here `chain` is a `Flux.Chain` object, or other "Flux model" such that
`Flux.params(chain)` returns the parameters to be optimised.

The `X` argument is the training features and `y` argument is the
target:

- `X` and `y` have type `Vector{<:Array{<:AbstractFloat}}`

- the shape of each elment of `X` is `(n1, n2, ..., nk, batch_size)`
  where `(n1, n2, ..., nk)` is the shape of the inputs of `chain`

- the shape of each element of `y` is `(m1, m2, ..., mk, batch_size)`
  where `(m1, m2, ..., mk)` is the shape of the `chain` outputs (even
  if `batch_size == 1`).

- the vectors `X` and `y` have the same length, coinciding with the
  total number of training batches.

The contribution to the objective function of a single training batch
`(X[i], y[i])` is

    loss(chain(X[i]), y[i]) + lambda*(model.alpha*l1) + (1 - model.alpha)*l2

where `l1 = sum(norm, params(chain)` and `l2 = sum(norm, params(chain))`.

One must have `acceleration isa CPU1` or `acceleration isa CUDALibs`
(for running on a GPU) where `CPU1` and `CUDALibs` are types defined
in `ComputationalResources.jl`.

### Return value

`(chain_trained, history)`, where `chain_trained` is a trained version
of `chain` (possibly moved to a gpu) and `history` is a vector of
losses - one intial loss, and one loss per epoch. The method may
mutate the argument `chain`, depending on cpu <-> gpu movements.

"""
function  fit!(chain, optimiser, loss, epochs,
               lambda, alpha, verbosity, acceleration, X, y)

    # intitialize and start progress meter:
    meter = Progress(epochs+1, dt=0, desc="Optimising neural net:",
                     barglyphs=BarGlyphs("[=> ]"), barlen=25, color=:yellow)
    verbosity != 1 || next!(meter)

    move = Mover(acceleration)
    X = move(X)
    y = move(y)
    chain = move(chain)

    loss_func(x, y) = loss(chain(x), y)

    # initiate history:
    prev_loss = mean(loss_func(X[i], y[i]) for i=1:length(X))
    history = [prev_loss,]

    for i in 1:epochs
        # We're taking data in a Flux-fashion.
#        @show i rand()
        train!(loss_func, Flux.params(chain), optimiser, X, y)
        current_loss =
            mean(loss_func(X[i], y[i]) for i=1:length(X))
        verbosity < 2 ||
            @info "Loss is $(round(current_loss; sigdigits=4))"
        push!(history, current_loss)

        # Early stopping is to be externally controlled.
        # So @ablaom has commented next 5 lines :
        # if current_loss == prev_loss
        #     @info "Model has reached maximum possible accuracy."*
        #     "More training won't increase accuracy"
        #     break
        # end

        prev_loss = current_loss
        verbosity != 1 || next!(meter)

    end

    return Flux.cpu(chain), history

end

# TODO: add callback functionality to above.

## BUILDING CHAINS A FROM HYPERPARAMETERS + INPUT/OUTPUT SHAPE

# We introduce chain builders as a way of exposing neural network
# hyperparameters (describing, architecture, dropout rates, etc) to
# the MLJ user. These parameters generally exclude the input/output
# dimensions/shapes, as the MLJ fit methods will determine these from
# the training data. A `Builder` object stores the parameters and an
# associated `fit` method builds a corresponding chain given the
# input/output dimensions/shape.

# Below n or (n1, n2) etc refers to network inputs, while m or (m1,
# m2) etc refers to outputs.

abstract type Builder <: MLJModelInterface.MLJType end

"""
    Linear(; σ=Flux.relu)

MLJFlux builder that constructs a fully connected two layer network
with activation function `σ`. The number of input and output nodes is
determined from the data.

"""
mutable struct Linear <: Builder
    σ
end
Linear(; σ=Flux.relu) = Linear(σ)
build(builder::Linear, n::Integer, m::Integer) =
    Flux.Chain(Flux.Dense(n, m, builder.σ))

"""
    Short(; n_hidden=0, dropout=0.5, σ=Flux.sigmoid)

MLJFlux builder that constructs a full-connected three-layer network
using `n_hidden` nodes in the hidden layer and the specified `dropout`
(defaulting to 0.5). An activation function `σ` is applied between the
hidden and final layers. If `n_hidden=0` (the default) then `n_hidden`
is the geometric mean of the number of input and output nodes.  The
number of input and output nodes is determined from the data.

"""
mutable struct Short <: Builder
    n_hidden::Int     # if zero use geometric mean of input/output
    dropout::Float64
    σ
end
Short(; n_hidden=0, dropout=0.5, σ=Flux.sigmoid) = Short(n_hidden, dropout, σ)
function build(builder::Short, n, m)
    n_hidden =
        builder.n_hidden == 0 ? round(Int, sqrt(n*m)) : builder.n_hidden
    return Flux.Chain(Flux.Dense(n, n_hidden, builder.σ),
                      Flux.Dropout(builder.dropout),
                       Flux.Dense(n_hidden, m))
end

## HELPERS

"""
    gpu_isdead()

Returns `true` if `acceleration=CUDALibs()` option is unavailable, and
false otherwise.

"""
gpu_isdead() = Flux.gpu([1.0, ]) isa Array

"""
    nrows(X)

Find the number of rows of `X`, where `X` is an `AbstractVector or
Tables.jl table.
"""
function nrows(X)
    Tables.istable(X) || throw(ArgumentError)
    Tables.columnaccess(X) || return length(collect(X))
    # if has columnaccess
    cols = Tables.columntable(X)
    !isempty(cols) || return 0
    return length(cols[1])
end
nrows(y::AbstractVector) = length(y)

reformat(X) = reformat(X, scitype(X))

# ---------------------------------
# Reformatting tables

reformat(X, ::Type{<:Table}) = MLJModelInterface.matrix(X)'

# ---------------------------------
# Reformatting images

reformat(X, ::Type{<:GrayImage}) =
    reshape(Float32.(X), size(X)..., 1)

function reformat(X, ::Type{<:AbstractVector{<:GrayImage}})
    ret = zeros(Float32, size(first(X))..., 1, length(X))
    for idx=1:size(ret, 4)
        ret[:, :, :, idx] .= reformat(X[idx])
    end
    return ret
end

function reformat(X, ::Type{<:ColorImage})
    ret = zeros(Float32, size(X)... , 3)
    for w = 1:size(X)[1]
        for h = 1:size(X)[2]
            ret[w, h, :] .= Float32.([X[w, h].r, X[w, h].g, X[w, h].b])
        end
    end
    return ret
end

function reformat(X, ::Type{<:AbstractVector{<:ColorImage}})
    ret = zeros(Float32, size(first(X))..., 3, length(X))
    for idx=1:size(ret, 4)
        ret[:, :, :, idx] .= reformat(X[idx])
    end
    return ret
end

# ------------------------------------------------------------
# Reformatting vectors of length n into matrices of dimension n * 1
# This enables compatibility with Flux's BatchNorm. This is currently
# used only in `predict`. In the future, when MLJ's "data front end"
# is implemented, `tomat` and the `reformat` of continuous vectors
# that follows will be collapsed and there will be some
# simplification.

function tomat end
tomat(x::Matrix) = x
tomat(x::Vector) = reshape(x, size(x, 1), 1)

# ------------------------------------------------------------
# Reformatting vectors of "scalar" types

reformat(y, ::Type{<:AbstractVector{<:Continuous}}) = reshape(y, 1, length(y))
function reformat(y, ::Type{<:AbstractVector{<:Finite}})
    levels = y |> first |> MLJModelInterface.classes
    return hcat([Flux.onehot(ele, levels) for ele in y]...,)
end

function reformat(y, ::Type{<:AbstractVector{<:Count}})
    levels = y |> first |> MLJModelInterface.classes
    return hcat([Flux.onehot(ele, levels) for ele in y]...,)
end

_get(Xmatrix::AbstractMatrix, b) = Xmatrix[:, b]
_get(y::AbstractVector, b) = y[b]

# each element in X is a single image of size (w, h, c)
_get(X::AbstractArray{<:Any, 4}, b) = X[:, :, :, b]

"""
    collate(model, X, y)

Return the Flux-friendly data object required by `MLJFlux.fit!`, given
input `X` and target `y` in the form required by
`MLJModelInterface.input_scitype(X)` and
`MLJModelInterface.target_scitype(y)`. (The batch size used is given
by `model.batch_size`.)

"""
function collate(model, X, y)
    row_batches = Base.Iterators.partition(1:nrows(y), model.batch_size)
    Xmatrix = reformat(X)
    ymatrix = reformat(y)
    return [_get(Xmatrix, b) for b in row_batches], [_get(ymatrix, b) for b in row_batches]
end


text:
"""
    get_mechanism(model; kwargs...)

    constructs mechanism

    model: name of mechanism
    kwargs: mechanism specific parameters
"""
function get_mechanism(model; kwargs...)
    # convert potential strings to symbols
    mech = eval(Symbol(:get, :_, string_to_symbol(model)))(; string_to_symbol(kwargs)...)
    return mech
end

"""
    initialize!(mechanism, model; kwargs)

    state initialization for mechanism

    mechanism: Mechanism
    model: name of mechanism
    kwargs: mechanism specific parameters
"""
function initialize!(mechanism::Mechanism, model; kwargs...)
    eval(Symbol(:initialize, :_, string_to_symbol(model), :!))(mechanism; kwargs...)
end


text:
using ExtremeFloats
using Test

@testset "conversion to Float64" begin
    numbers = rand(100)
    for num in numbers
        @test (num |> ExtremeFloat |> Float64) ≈ num
    end
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export split

"""
    split(;kwargs...)
    split(children::Any;kwargs...)
    split(children_maker::Function;kwargs...)

A Split component.
ExampleComponent is an example component.
It takes a property, `label`, and
displays it.
It renders an input with the property `value`
which is editable by the user.
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): The children to render in the split
- `id` (String; optional): The ID used to identify this component in Dash callbacks.
- `className` (String; optional): The class of the container (div)
- `cursor` (String; optional): Cursor to display while dragging.
- `direction` (a value equal to: "horizontal", "vertical"; optional): Direction to split: horizontal or vertical.
- `dragInterval` (Real; optional): Number of pixels to drag.
- `expandToMin` (Bool; optional): Grow initial sizes to minSize
- `gutterAlign` (String; optional): Gutter alignment between elements.
- `gutterSize` (Real; optional): Gutter size in pixels.
- `maxSize` (Real | Array of Reals; optional): Maximum size of each element.
- `minSize` (Real | Array of Reals; optional): Minimum size of each element.
- `sizes` (Array of Reals; optional): Initial sizes of each element in percents or CSS values.
- `snapOffset` (Real | Array of Reals; optional): Snap to minimum size offset in pixels.
- `style` (Dict; optional): The style of the container (div)
"""
function split(; kwargs...)
        available_props = Symbol[:children, :id, :className, :cursor, :direction, :dragInterval, :expandToMin, :gutterAlign, :gutterSize, :maxSize, :minSize, :sizes, :snapOffset, :style]
        wild_props = Symbol[]
        return Component("split", "Split", "dash_split", available_props, wild_props; kwargs...)
end

split(children::Any; kwargs...) = split(;kwargs..., children = children)
split(children_maker::Function; kwargs...) = split(children_maker(); kwargs...)

max_stars_repo_name: TCatshoek/dash-split

max_issues_repo_name: TCatshoek/dash-split

max_forks_repo_name: TCatshoek/dash-split


text:
"""
    alpha(p::Int)

Returns coefficients for the spline kernel of order `p`.
"""
function alpha(p::Int)
    if p < 1
        throw(DomainError("Spline order has to be strictly positive"))
    end
    a = zeros(p)
    for i = 0:p-1
        for j = 0:i
            for k = i:p-1
                a[i+1] = a[i+1] + (-1.0)^(k-j)/((k+j+1)*
                        factorial(j)*
                        factorial(p-1)*
                        factorial(p-1-k)*
                        factorial(i-j)*
                        factorial(k-i))
            end
        end
    end
    return a
end

"""
    spline_kernel(t, p)

Returns generators `U` and `V` for the spline kernel of order `p` given a list of knots `t`.
"""
function spline_kernel(t::AbstractArray, p::Int)

    fp = factorial.(p-1:-1:0)
    a = alpha(p).*fp
    Ut = (repeat(t,p,1).^Vector(p-1:-1:0))./fp
    Vt = (repeat(t,p,1).^Vector(p:2*p-1)).*a

    return Ut,Vt
end

"""
    spline_kernel_matrix(U, V)

Returns the dense spline kernel matrix with generators `U` and `V`.
"""
function spline_kernel_matrix(U, V)
    return tril(U*V') + triu(V*U',1)
end

"""
    spline_kernel_matrix(U, V, d)

Returns the dense spline kernel matrix with generators `U` and `V` plus `Diagonal(d)`.
"""
function spline_kernel_matrix(U, V, d )
    return spline_kernel_matrix(U, V) + Diagonal(d)
end

"""
    logp(Σ, T, y, σf, σn)

Generalized Log-likelihood evalation for the spline kernel.
"""
function logp(Σ, T, y, σf, σn)
    m = rank(T)
    n = length(y)
    Ki = cholesky(σf^2*Σ  + σn^2*I)
    if LinearAlgebra.issuccess(Ki) == false
        throw(DomainError("Cholesky Factorization failed"))
    end
    A = T'*(Ki\T)
    C = (Ki\T)*(A\(T'*inv(Ki)))
    lp = -0.5*y'*(Ki\y)  +
          0.5*y'*(C*y)   +
         -0.5*logdet(Ki) +
         -0.5*logdet(A)  +
         -(n-m)*0.5*log(2*π)
    return lp
end

"""
    dlogp(Σ, T, y, σf, σn)

Derivatives of the generalized Log-likelihood for the spline kernel.
"""
function dlogp(Σ, T, y, σf, σn)
    Ki = cholesky(σf^2*Σ  + σn^2*I)
    if LinearAlgebra.issuccess(Ki) == false
        throw(DomainError("Cholesky Factorization failed"))
    end
    A = T'*(Ki\T)
    # TO_DO: Solving for the identity?
    dKf = 2.0*σf*Σ
    dAf = -T'*(Ki\dKf)*(Ki\T)
    dKn = 2.0*σn*I
    dAn = -T'*(inv(Ki)*dKn)*(Ki\T)
    Pf = Ki\dKf*inv(Ki)
    Pn = Ki\(dKn*inv(Ki))
    G = T*(A\(T'*(Ki\y)))
    dlpf =  0.5*y'*(Pf*y)                +
            0.5*(-2*(y'*Pf)*G + G'*Pf*G) +
           -0.5*tr(Ki\dKf)               +
           -0.5*tr(A\dAf)
   dlpn =   0.5*y'*(Pn*y)                   +
            0.5*(-2*(y'*Pn)*G + G'*Pn*G)    +
           -0.5*tr(inv(Ki)*dKn)             +
           -0.5*tr(A\dAn)
    return [dlpf; dlpn]
end

"""
    create_T(t, m)

Creates `m`-order Taylor-polynomial basis from the knots in `t`.
"""
function create_T(t, m)
    T = ones(length(t), m)
    for ν = 2:m
        T[:, ν] = t.^(ν-1)/factorial(ν - 1)
    end
    return T
end


text:
using Jl2Py
using Test

@testset "Jl2Py.jl" begin
    # Literal constants
    @test jl2py("2") == "2"
    @test jl2py("2.0") == "2.0"
    @test jl2py("\"hello\"") == "'hello'"
    @test jl2py("false") == "False"
    @test jl2py("true") == "True"
    @test jl2py("foo") == "foo"

    # UAdd & USub
    @test jl2py("+3") == "3"
    @test jl2py("-3") == "-3"
    @test jl2py("+a") == "+a"
    @test jl2py("-a") == "-a"

    # Add
    @test jl2py("1 + 1") == "1 + 1"
    @test jl2py("1 + 1 + 1") == "1 + 1 + 1"
    @test jl2py("1 + 1 + 1 + 1") == "1 + 1 + 1 + 1"
    @test jl2py("a + 2") == "a + 2"
    @test jl2py("a + b") == "a + b"

    # Sub
    @test jl2py("1 - 1") == "1 - 1"
    @test jl2py("a - 2") == "a - 2"
    @test jl2py("a - b") == "a - b"

    # Mult
    @test jl2py("1 * 2") == "1 * 2"
    @test jl2py("1 * 2 * 3") == "1 * 2 * 3"

    # Div & FloorDiv
    @test jl2py("1 / 2") == "1 / 2"
    @test jl2py("1 ÷ 3") == "1 // 3"
    @test jl2py("div(3, 1)") == "3 // 1"

    # Mod
    @test jl2py("1 % 2") == "1 % 2"
    @test jl2py("mod(3, 1)") == "3 % 1"

    # Pow
    @test jl2py("10 ^ 5") == "10 ** 5"

    # Complex arithmetic
    @test jl2py("(1 + 5) * (2 - 5) / (3 * 6)") == "(1 + 5) * (2 - 5) / (3 * 6)"

    # Binary comparisons
    @test jl2py("1 == 2") == "1 == 2"
    @test jl2py("1 === 2") == "1 == 2"
    @test jl2py("1 != 2") == "1 != 2"
    @test jl2py("1 < 2") == "1 < 2"
    @test jl2py("1 <= 2") == "1 <= 2"
    @test jl2py("1 > 2") == "1 > 2"
    @test jl2py("1 >= 2") == "1 >= 2"

    # Multiple comparisons
    @test jl2py("1 < 2 < 3") == "1 < 2 < 3"
    @test jl2py("1 < 2 < 3 < 4") == "1 < 2 < 3 < 4"
    @test jl2py("3 > 2 != 3 == 3") == "3 > 2 != 3 == 3"

    # List
    @test jl2py("[1, 2, 3]") == "[1, 2, 3]"
    @test jl2py("[1, 2 + 3, 3]") == "[1, 2 + 3, 3]"
    @test jl2py("[1, 2, [3, 4], []]") == "[1, 2, [3, 4], []]"

    # Assign
    @test jl2py("a = 2") == "a = 2"
    @test jl2py("a = 2 + 3") == "a = 2 + 3"
    @test jl2py("a = b = 2") == "a = b = 2"
    @test jl2py("a = b = c = 23 + 3") == "a = b = c = 23 + 3"

    # Function call
    @test jl2py("print(2)") == "print(2)"
end


text:
module ParamLevelSet
using LinearAlgebra
using SparseArrays
using jInv.Mesh
using jInv.Utils
using NearestNeighbors

export centerHeaviside,deltaHeaviside,getDefaultHeaviside
centerHeavySide = 0.3;
deltaHeavySide = 0.05;

function getDefaultHeaviside()
	return (u,dsu)->heaviside!(u,dsu,ParamLevelSet.centerHeavySide,ParamLevelSet.deltaHeavySide);
end

## Other things to include: bf type in psi.
## Thresholding parameters for the Jacobians.

function getIdentityMat(n::Int64)
return Matrix(1.0I,n,n);

end
include("SpMatBuilder.jl");
include("RBFs.jl");
include("MeshFreeRBFs.jl");
include("RBFs_extended.jl");
include("heaviside.jl");
include("RBFRotation3D.jl");
include("SPD_regularization.jl");
include("rotation3D.jl")

end


text:
@doc doc"Hyperbolic tangent transfer layer" ->
type TanH
    output::DenseRealTensor
    grad_input::DenseRealTensor
end

function forward!(m::TanH, input::DenseMatrix)
    self.output = tanh(input)
end

function compute_inputgradient!(m::TanH, grad_output::DenseMatrix)
    @simd for i = eachindex(grad_output.values)
        self.grad_input[i] = (one(F) - m.output[i] * m.output[i]) * grad_output[i]
    end
    self.grad_input
end


text:
function custatevecApplyGeneralizedPermutationMatrix(handle, sv, svDataType, nIndexBits, permutation, diagonals, diagonalsDataType, adjoint, basisBits, nBasisBits, maskBitString, maskOrdering, maskLen, extraWorkspace, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecApplyGeneralizedPermutationMatrix, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{custatevecIndex_t}, Ptr{Cvoid}, cudaDataType_t, Int32, Ptr{Int32}, UInt32, Ptr{Int32}, Ptr{Int32}, UInt32, Ptr{Cvoid}, Csize_t), handle, sv, svDataType, nIndexBits, permutation, diagonals, diagonalsDataType, adjoint, basisBits, nBasisBits, maskBitString, maskOrdering, maskLen, extraWorkspace, extraWorkspaceSizeInBytes)
end

function custatevecGetErrorName(status)
    initialize_context()
    ccall((:custatevecGetErrorName, libcustatevec()), Cstring, (custatevecStatus_t,), status)
end

function custatevecAccessor_create(handle, sv, svDataType, nIndexBits, accessor, bitOrdering, bitOrderingLen, maskBitString, maskOrdering, maskLen, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecAccessor_create, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{custatevecAccessorDescriptor_t}, Ptr{Int32}, UInt32, Ptr{Int32}, Ptr{Int32}, UInt32, Ptr{Csize_t}), handle, sv, svDataType, nIndexBits, accessor, bitOrdering, bitOrderingLen, maskBitString, maskOrdering, maskLen, extraWorkspaceSizeInBytes)
end

function custatevecCreate(handle)
    initialize_context()
    ccall((:custatevecCreate, libcustatevec()), custatevecStatus_t, (Ptr{custatevecHandle_t},), handle)
end

function custatevecExpectationsOnPauliBasis(handle, sv, svDataType, nIndexBits, expectationValues, pauliOperatorsArray, basisBitsArray, nBasisBitsArray, nPauliOperatorArrays)
    initialize_context()
    ccall((:custatevecExpectationsOnPauliBasis, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{Cdouble}, Ptr{Ptr{custatevecPauli_t}}, Ptr{Ptr{Int32}}, Ptr{UInt32}, UInt32), handle, sv, svDataType, nIndexBits, expectationValues, pauliOperatorsArray, basisBitsArray, nBasisBitsArray, nPauliOperatorArrays)
end

function custatevecApplyMatrix(handle, sv, svDataType, nIndexBits, matrix, matrixDataType, layout, adjoint, targets, nTargets, controls, nControls, controlBitValues, computeType, extraWorkspace, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecApplyMatrix, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, CuPtr{Cvoid}, cudaDataType_t, UInt32, PtrOrCuPtr{Cvoid}, cudaDataType_t, custatevecMatrixLayout_t, Int32, Ptr{Int32}, UInt32, Ptr{Int32}, UInt32, Ptr{Int32}, custatevecComputeType_t, CuPtr{Cvoid}, Csize_t), handle, sv, svDataType, nIndexBits, matrix, matrixDataType, layout, adjoint, targets, nTargets, controls, nControls, controlBitValues, computeType, extraWorkspace, extraWorkspaceSizeInBytes)
end

function custatevecAccessor_setExtraWorkspace(handle, accessor, extraWorkspace, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecAccessor_setExtraWorkspace, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, Ptr{custatevecAccessorDescriptor_t}, Ptr{Cvoid}, Csize_t), handle, accessor, extraWorkspace, extraWorkspaceSizeInBytes)
end

function custatevecSampler_create(handle, sv, svDataType, nIndexBits, sampler, nMaxShots, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecSampler_create, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, CuPtr{Cvoid}, cudaDataType_t, UInt32, Ptr{custatevecSamplerDescriptor_t}, UInt32, Ptr{Csize_t}), handle, sv, svDataType, nIndexBits, sampler, nMaxShots, extraWorkspaceSizeInBytes)
end

function custatevecCollapseOnZBasis(handle, sv, svDataType, nIndexBits, parity, basisBits, nBasisBits, norm)
    initialize_context()
    ccall((:custatevecCollapseOnZBasis, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Int32, Ptr{Int32}, UInt32, Cdouble), handle, sv, svDataType, nIndexBits, parity, basisBits, nBasisBits, norm)
end

function custatevecExpectation(handle, sv, svDataType, nIndexBits, expectationValue, expectationDataType, residualNorm, matrix, matrixDataType, layout, basisBits, nBasisBits, computeType, extraWorkspace, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecExpectation, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, PtrOrCuPtr{Cvoid}, cudaDataType_t, UInt32, Ptr{Cvoid}, cudaDataType_t, Ptr{Cdouble}, Ptr{Cvoid}, cudaDataType_t, custatevecMatrixLayout_t, Ptr{Int32}, UInt32, custatevecComputeType_t, CuPtr{Cvoid}, Csize_t), handle, sv, svDataType, nIndexBits, expectationValue, expectationDataType, residualNorm, matrix, matrixDataType, layout, basisBits, nBasisBits, computeType, extraWorkspace, extraWorkspaceSizeInBytes)
end

function custatevecBatchMeasure(handle, sv, svDataType, nIndexBits, bitString, bitOrdering, bitStringLen, randnum, collapse)
    initialize_context()
    ccall((:custatevecBatchMeasure, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{Int32}, Ptr{Int32}, UInt32, Cdouble, custatevecCollapseOp_t), handle, sv, svDataType, nIndexBits, bitString, bitOrdering, bitStringLen, randnum, collapse)
end

function custatevecExpectation_bufferSize(handle, svDataType, nIndexBits, matrix, matrixDataType, layout, nBasisBits, computeType, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecExpectation_bufferSize, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, cudaDataType_t, UInt32, Ptr{Cvoid}, cudaDataType_t, custatevecMatrixLayout_t, UInt32, custatevecComputeType_t, Ptr{Csize_t}), handle, svDataType, nIndexBits, matrix, matrixDataType, layout, nBasisBits, computeType, extraWorkspaceSizeInBytes)
end

function custatevecApplyMatrix_bufferSize(handle, svDataType, nIndexBits, matrix, matrixDataType, layout, adjoint, nTargets, nControls, computeType, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecApplyMatrix_bufferSize, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, cudaDataType_t, UInt32, Ptr{Cvoid}, cudaDataType_t, custatevecMatrixLayout_t, Int32, UInt32, UInt32, custatevecComputeType_t, Ptr{Csize_t}), handle, svDataType, nIndexBits, matrix, matrixDataType, layout, adjoint, nTargets, nControls, computeType, extraWorkspaceSizeInBytes)
end

function custatevecLoggerSetMask(mask)
    ccall((:custatevecLoggerSetMask, libcustatevec()), custatevecStatus_t, (Int32,), mask)
end

function custatevecAccessor_createReadOnly(handle, sv, svDataType, nIndexBits, accessor, bitOrdering, bitOrderingLen, maskBitString, maskOrdering, maskLen, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecAccessor_createReadOnly, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{custatevecAccessorDescriptor_t}, Ptr{Int32}, UInt32, Ptr{Int32}, Ptr{Int32}, UInt32, Ptr{Csize_t}), handle, sv, svDataType, nIndexBits, accessor, bitOrdering, bitOrderingLen, maskBitString, maskOrdering, maskLen, extraWorkspaceSizeInBytes)
end

function custatevecGetStream(handle, streamId)
    initialize_context()
    ccall((:custatevecGetStream, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, Ptr{CUstream}), handle, streamId)
end

function custatevecAccessor_get(handle, accessor, externalBuffer, _begin, _end)
    initialize_context()
    ccall((:custatevecAccessor_get, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, Ptr{custatevecAccessorDescriptor_t}, Ptr{Cvoid}, custatevecIndex_t, custatevecIndex_t), handle, accessor, externalBuffer, _begin, _end)
end

function custatevecGetVersion()
    ccall((:custatevecGetVersion, libcustatevec()), Csize_t, ())
end

function custatevecLoggerForceDisable()
    ccall((:custatevecLoggerForceDisable, libcustatevec()), custatevecStatus_t, ())
end

function custatevecGetDefaultWorkspaceSize(handle, workspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecGetDefaultWorkspaceSize, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, Ptr{Csize_t}), handle, workspaceSizeInBytes)
end

function custatevecSetStream(handle, streamId)
    initialize_context()
    ccall((:custatevecSetStream, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, CUstream), handle, streamId)
end

function custatevecSetWorkspace(handle, workspace, workspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecSetWorkspace, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, Csize_t), handle, workspace, workspaceSizeInBytes)
end

function custatevecApplyExp(handle, sv, svDataType, nIndexBits, theta, paulis, targets, nTargets, controls, controlBitValues, nControls)
    initialize_context()
    ccall((:custatevecApplyExp, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Cdouble, Ptr{custatevecPauli_t}, Ptr{Int32}, UInt32, Ptr{Int32}, Ptr{Int32}, UInt32), handle, sv, svDataType, nIndexBits, theta, paulis, targets, nTargets, controls, controlBitValues, nControls)
end

function custatevecAbs2SumArray(handle, sv, svDataType, nIndexBits, abs2sum, bitOrdering, bitOrderingLen, maskBitString, maskOrdering, maskLen)
    initialize_context()
    ccall((:custatevecAbs2SumArray, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{Cdouble}, Ptr{Int32}, UInt32, Ptr{Int32}, Ptr{Int32}, UInt32), handle, sv, svDataType, nIndexBits, abs2sum, bitOrdering, bitOrderingLen, maskBitString, maskOrdering, maskLen)
end

function custatevecGetProperty(type, value)
    initialize_context()
    ccall((:custatevecGetProperty, libcustatevec()), custatevecStatus_t, (libraryPropertyType, Ptr{Int32}), type, value)
end

function custatevecAbs2SumOnZBasis(handle, sv, svDataType, nIndexBits, abs2sum0, abs2sum1, basisBits, nBasisBits)
    initialize_context()
    ccall((:custatevecAbs2SumOnZBasis, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Int32}, UInt32), handle, sv, svDataType, nIndexBits, abs2sum0, abs2sum1, basisBits, nBasisBits)
end

function custatevecLoggerOpenFile(logFile)
    ccall((:custatevecLoggerOpenFile, libcustatevec()), custatevecStatus_t, (Cstring,), logFile)
end

function custatevecDestroy(handle)
    initialize_context()
    ccall((:custatevecDestroy, libcustatevec()), custatevecStatus_t, (custatevecHandle_t,), handle)
end

function custatevecLoggerSetLevel(level)
    initialize_context()
    ccall((:custatevecLoggerSetLevel, libcustatevec()), custatevecStatus_t, (Int32,), level)
end

function custatevecApplyGeneralizedPermutationMatrix_bufferSize(handle, svDataType, nIndexBits, permutation, diagonals, diagonalsDataType, basisBits, nBasisBits, maskLen, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecApplyGeneralizedPermutationMatrix_bufferSize, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, cudaDataType_t, UInt32, Ptr{custatevecIndex_t}, Ptr{Cvoid}, cudaDataType_t, Ptr{Int32}, UInt32, UInt32, Ptr{Csize_t}), handle, svDataType, nIndexBits, permutation, diagonals, diagonalsDataType, basisBits, nBasisBits, maskLen, extraWorkspaceSizeInBytes)
end

function custatevecLoggerSetFile(file)
    ccall((:custatevecLoggerSetFile, libcustatevec()), custatevecStatus_t, (Ptr{Cvoid},), file)
end

function custatevecGetErrorString(status)
    ccall((:custatevecGetErrorString, libcustatevec()), Cstring, (custatevecStatus_t,), status)
end

function custatevecLoggerSetCallback(callback)
    ccall((:custatevecLoggerSetCallback, libcustatevec()), custatevecStatus_t, (custatevecLoggerCallback_t,), callback)
end

function custatevecCollapseByBitString(handle, sv, svDataType, nIndexBits, bitString, bitOrdering, bitStringLen, norm)
    initialize_context()
    ccall((:custatevecCollapseByBitString, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{Int32}, Ptr{Int32}, UInt32, Cdouble), handle, sv, svDataType, nIndexBits, bitString, bitOrdering, bitStringLen, norm)
end

function custatevecSampler_preprocess(handle, sampler, extraWorkspace, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecSampler_preprocess, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, Ptr{custatevecSamplerDescriptor_t}, CuPtr{Cvoid}, Csize_t), handle, sampler, extraWorkspace, extraWorkspaceSizeInBytes)
end

function custatevecMeasureOnZBasis(handle, sv, svDataType, nIndexBits, parity, basisBits, nBasisBits, randnum, collapse)
    initialize_context()
    ccall((:custatevecMeasureOnZBasis, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{Int32}, Ptr{Int32}, UInt32, Cdouble, custatevecCollapseOp_t), handle, sv, svDataType, nIndexBits, parity, basisBits, nBasisBits, randnum, collapse)
end

function custatevecSampler_sample(handle, sampler, bitStrings, bitOrdering, bitStringLen, randnums, nShots, output)
    initialize_context()
    ccall((:custatevecSampler_sample, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, Ptr{custatevecSamplerDescriptor_t}, Ptr{custatevecIndex_t}, Ptr{Int32}, UInt32, Ptr{Cdouble}, UInt32, custatevecSamplerOutput_t), handle, sampler, bitStrings, bitOrdering, bitStringLen, randnums, nShots, output)
end

function custatevecAccessor_set(handle, accessor, externalBuffer, _begin, _end)
    initialize_context()
    ccall((:custatevecAccessor_set, libcustatevec()), custatevecStatus_t, (custatevecHandle_t, Ptr{custatevecAccessorDescriptor_t}, Ptr{Cvoid}, custatevecIndex_t, custatevecIndex_t), handle, accessor, externalBuffer, _begin, _end)
end


text:
# Useful utilities from Base
using Base: sign_mask, exponent_mask, significand_mask, exponent_one, exponent_bias, significand_bits

using Base.Math: @horner

# Similar to @horner, but converts coefficients to same type as x
macro horner_oftype(x, p...)
    ex = :(oftype($x,$(esc(p[end]))))
    for i = length(p)-1:-1:1
        ex = :(muladd(t, $ex, oftype($x,$(esc(p[i])))))
    end
    Expr(:block, :(t = $(esc(x))), ex)
end

"""
    highword(x::Union{Float32,Float64}) -> UInt32

Get the most significant 32 bits as a `UInt32` from `x`.
Corresponds to `GET_HIGH_WORD` in musl
"""
@inline highword(x::Float64) = UInt32(reinterpret(UInt64, x) >> 32)
@inline highword(x::Float32) =  reinterpret(UInt32, x)

"""
    trunclo(x::Union{Float32,Float64})

Truncates the lower order bits of `x` so that the result of the multiplication `trunclo(x)
* trunclo(y)` is exact, assuming no underflow or overflow occurs.

This relies on the following property: if `a` has `n` significant bits, and `b` has `m`
significant bits, then exact product `a*b` has either `n+m-1` or `n+m` significant bits:

* `Float64` has 53 significant bits (including implicit leading bit), so we need to
  truncate the last 27 significant bits (leaving 26 bits).

* `Float32` has 24 significant bits (including implicit leading bit), so we need to
  truncate the last 12 significant bits (leaving 12 bits).

This is typically faster than other methods of truncating lower order bits (such as
Veltkamp splitting, or converting `Float64`s to `Float32`s and back again). For LLVM 3.8
or greater, this should give optimal `ANDPD`/`ANDSD` instructions on supported x86
architectures, which doesn't require moving registers (Julia issue #9868).

NOTE: For odd significand sizes, such as `Float64`, when used as a replacement Veltkamp
splitting for computing extended precision multiplications with a Dekker-style `mul12`
algorithm, it can lose the last bit of precision. For example the function:

    function incorrect_mul2(x,y)
        hx = trunclo(x); tx = x-hx
        hy = trunclo(y); ty = y-hy
        p = hx*hy
        q = hx*ty + tx*hy
        z = p+q
        zz = p-z+q+tx*ty
        z, zz
    end

will give an incorrect result for the case `x = y = 0x1.800000e000001p+0`.
"""
@inline trunclo(x::Float64) =
    reinterpret(Float64, reinterpret(UInt64,x) & 0xffff_ffff_f800_0000)
@inline trunclo(x::Float32) =
    reinterpret(Float32, reinterpret(UInt32,x) & 0xffff_f000)

# determine if hardware FMA is available
# should probably check with LLVM, see #9855.
"""
    is_fma_fast(T)

Checks if the `fma` function is fast for the floating point type `T`: typically is it a
native instruction (`true`) or does it fall back on a software implementation (`false`).
"""
function is_fma_fast end
for T in (Float32, Float64)
    @eval is_fma_fast(::Type{$T}) = $(muladd(nextfloat(one(T)),nextfloat(one(T)),-nextfloat(one(T),2)) != zero(T))
end


text:
# ---
# title: 129. Sum Root to Leaf Numbers
# id: problem129
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/sum-root-to-leaf-numbers/description/>
# hidden: true
# ---
#
# Given a binary tree containing digits from `0-9` only, each root-to-leaf path
# could represent a number.
#
# An example is the root-to-leaf path `1->2->3` which represents the number
# `123`.
#
# Find the total sum of all root-to-leaf numbers.
#
# **Note:**  A leaf is a node with no children.
#
# **Example:**
#
#
#
#     Input: [1,2,3]
#         1
#        / \
#       2   3
#     Output: 25
#     Explanation:
#     The root-to-leaf path 1->2 represents the number 12.
#     The root-to-leaf path 1->3 represents the number 13.
#     Therefore, sum = 12 + 13 = 25.
#
# **Example 2:**
#
#
#
#     Input: [4,9,0,5,1]
#         4
#        / \
#       9   0
#      / \
#     5   1
#     Output: 1026
#     Explanation:
#     The root-to-leaf path 4->9->5 represents the number 495.
#     The root-to-leaf path 4->9->1 represents the number 491.
#     The root-to-leaf path 4->0 represents the number 40.
#     Therefore, sum = 495 + 491 + 40 = 1026.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
import Base:length

length(heap::Heap) = length(heap.cells) - heap.total_used

isless(hc1::HeapCell, hc2::HeapCell)  =
    isless(hc2.score, hc1.score) # because we want maximum and not minimum

function push!(heap::Heap, action_ix::Int32, parent_id::Int32, score::Float32)
    heap_cell = HeapCell(
        false, # not used
        score,
        action_ix,
        parent_id,
    )
    heappush!(heap.cells, heap_cell)
end

function pop_max!(heap::Heap)
    length(heap) == 0 && error("empty heap in pop_max!")
    if heap.total_used / length(heap) > HEAP_GC_FAC
        garbage_collect!(heap)
    end
    while true
        hc = heappop!(heap.cells)
        if hc.is_used
            heap.total_used -= 1
            continue
        end
        return hc.action_ix, hc.parent_id
    end
end

function pop_rand!(heap::Heap)
    length(heap) == 0 && error("empty heap in pop_rand!")
    if heap.total_used / length(heap) > HEAP_GC_FAC
        garbage_collect!(heap)
    end
    while true
        ix = rand(1:length(heap.cells))
        hc = heap.cells[ix]
        hc.is_used && continue
        # mark as used:
        heap.cells[ix] = HeapCell(
            true,
            hc.score,
            hc.action_ix,
            hc.parent_id,
        )
        heap.total_used += 1
        return hc.action_ix, hc.parent_id
    end
end

function garbage_collect!(heap::Heap)
    tmp = HeapCell[]
    for hc in heap.cells
        hc.is_used && continue
        push!(tmp, hc)
    end
    empty!(heap.cells)
    for hc in tmp
        heappush!(heap.cells, hc)
    end
    heap.total_used = 0
end


text:
# try to emulate the output of $ ./lasertag -r 7
using LaserTag
using POMDPs
using Test

let
    model = LaserTagPOMDP{DESPOTEmu, DMeas}(
        obstacles = Set{Coord}(
            Coord(xy...) for xy in [
                (3,1),
                (6,4),
                (10,4),
                (10,5),
                (2,6),
                (6,6),
                (1,7),
                (5,7)
            ]
        )
    )

    is = LTState(Coord(1,1), Coord(10,7), false)

    show(stdout, MIME("text/plain"), LaserTagVis(model, s=is))

    a = Dict{Symbol, Int}(:North=>1, :South=>2, :East=>3, :West=>4, :Tag=>5)

    aspor = [
        (a[:South], LTState([1,2], [11,7], false), DMeas(4, 9, 0, 0, 6, 2, 0, 0), -1.0),
        (a[:East], LTState([2,2], [11,7], false), DMeas(1, 8, 0, 0, 2, 0, 1, 1), -1.0),
        (a[:East], LTState([3,2], [11,7], false), DMeas(5, 8, 0, 0, 7, 1, 1, 0), -1.0),
        (a[:East], LTState([4,2], [11,7], false), DMeas(4, 7, 0, 0, 0, 1, 0, 3), -1.0),
        (a[:East], LTState([5,2], [11,7], false), DMeas(2, 3, 0, 2, 4, 0, 0, 5), -1.0),
        (a[:East], LTState([6,2], [11,7], false), DMeas(0, 3, 1, 3, 2, 0, 2, 5), -1.0),
        (a[:East], LTState([7,2], [11,7], false), DMeas(4, 4, 0, 6, 4, 2, 0, 6), -1.0),
        (a[:South], LTState([7,3], [11,7], false), DMeas(4, 2, 0, 6, 4, 4, 2, 0), -1.0),
        (a[:South], LTState([7,4], [11,7], false), DMeas(2, 2, 3, 0, 3, 4, 3, 4), -1.0),
        (a[:South], LTState([7,5], [11,7], false), DMeas(1, 2, 4, 4, 2, 6, 0, 0), -1.0),
        (a[:East], LTState([8,5], [11,7], false), DMeas(1, 1, 2, 5, 2, 4, 6, 1), -1.0),
        (a[:South], LTState([8,6], [11,7], false), DMeas(1, 2, 4, 0, 0, 0, 1, 2), -1.0),
        (a[:East], LTState([9,6], [11,7], false), DMeas(0, 2, 4, 0, 1, 0, 6, 2), -1.0),
        (a[:East], LTState([10,6], [11,7], false), DMeas(0, 0, 0, 1, 0, 1, 4, 0), -1.0),
        (a[:South], LTState([10,7], [11,7], false), DMeas(0, 0, 1, 2, 1, 0, 9, 0), -1.0),
        (a[:East], LTState([11,7], [11,7], false), LaserTag.D_SAME_LOC, -1.0),
        (a[:Tag], LTState([11,7], [11,7], true), LaserTag.D_SAME_LOC, 10.0),
    ]

    oprobs = [
        0.00258604,
        0.000232729,
        0.000177984,
        0.000319759,
        1.58596e-05,
        1.2223e-05,
        3.82869e-05,
        3.78645e-05,
        0.000100641,
        0.000288588,
        1.19013e-05,
        0.00023327,
        0.000348013,
        0.000925637,
        0.000897613,
        1,
        1
    ]

    discounted_rsum = 0.0
    disc = 1.0
    @show discount(model)
    for i in 1:length(aspor)
        a, sp, o, r = aspor[i]
        if i == 1
            s = is
        else
            s = aspor[i-1][2]
        end
        od = observation(model, s, a, sp)
        @test isapprox(pdf(od, o), oprobs[i], rtol=0.001)
        @test isapprox(reward(model, s, a, sp), r, atol=0.0001)
        @test !isterminal(model, s)
        discounted_rsum = discounted_rsum + disc*r
        @show discounted_rsum
        disc*=discount(model)
    end
    @show discounted_rsum

    @test isterminal(model, aspor[end][2])
    @test isapprox(discounted_rsum, -6.7962, atol=0.0001)

    println("Matches CPP!")
end