text:
#MC Hammer Charts for simulation results

using Gadfly

function density_chrt(Data, x_label)
      if x_label ==""
            x_label="x"
      end
chart = plot(Data, x=Data, Guide.xlabel(x_label), Guide.ylabel("Frequency"), Geom.density)
describe(Data)
println("")
print("Mean: ", mean(ProjectNPV),"\n")
print("Std.Dev: ", std(ProjectNPV),"\n")
print("Prob. of Neg.: ", GetCertainty(ProjectNPV,0,0),"\n")
println("")
print("p10, p50, p90 : ", quantile(collect(Float64, ProjectNPV),[0.1,0.5,0.9]),"\n")
return chart
end

function histogram_chrt(Data, x_label)
      if x_label ==""
            x_label="x"
      end
chart = plot(Data, x=Data, Guide.xlabel(x_label), Guide.ylabel("Frequency"), Geom.histogram)
describe(Data)
println("")
print("Mean: ", mean(ProjectNPV),"\n")
print("Std.Dev: ", std(ProjectNPV),"\n")
print("Prob. of Neg.: ", GetCertainty(ProjectNPV,0,0),"\n")
println("")
print("p10, p50, p90 : ", quantile(collect(Float64, ProjectNPV),[0.1,0.5,0.9]),"\n")
return chart
end

#SENSITIVITY CHART (ArrayName = Sensitivity_Tbl; TargetCol=1; RankOrder=1)

function sensitivity_chrt(ArrayName, TargetCol, RankOrder)
cor_mat = []
M_Size = size(ArrayName,2)
if RankOrder == 1
      for i=TargetCol
      cor_vector = []
            for i2=1:M_Size
             cor_i = cor(tiedrank(ArrayName[i]),tiedrank(ArrayName[i2]));
             push!(cor_vector,cor_i)
             #print(cov_i)
            end
       push!(cor_mat,cor_vector)
      end
else
      for i=1:TargetCol
      cor_vector = []
            for i2=1:M_Size
             cor_i = cor(ArrayName[i], ArrayName[i2]);
             push!(cor_vector,cor_i)
             #print(cov_i)
      end
      push!(cor_mat,cor_vector)
      end
end

correl_vals = hcat(cor_mat...)

#This removes all entries of 1 (self)
#correl_vals = correl_vals[(correl_vals[:].!=1)]

#setup input names, next version will get variable names

ystr = names(ArrayName)
color_code = correl_vals .<0

impact=[]
for i in 1:size(correl_vals,1)
      if color_code[i] == true
            value="Negative"
      else
            value="Positive"
      end
      push!(impact, value)
end

#Contribution to Variance
cont_var = correl_vals .^2
cont_var = cont_var ./ (sum(cont_var)-1)

graph_tbl = DataFrame(hcat(ystr,correl_vals, abs.(correl_vals), cont_var, color_code, impact))
names!(graph_tbl, [:name, :correlation, :abs_cor, :cont_var, :color_code, :impact])
graph_tbl = graph_tbl[graph_tbl[:correlation].!=1,:]
#graph_tbl = vcat(graph_tbl, cont_var)

#This removes all entries of 1 (self) in DataTable
graph_tbl = graph_tbl[graph_tbl[:correlation].!=1,:]

#Sort
graph_tbl = sort(graph_tbl,3,rev=false)

println(graph_tbl)

return plot(graph_tbl, y=:name, x=:correlation, Guide.Title("Variables with Biggest Impact"), Guide.xlabel("Correlation Score"), Guide.ylabel("Input"), Geom.bar(orientation=:horizontal), Theme(bar_spacing=10pt), color=:impact, Scale.color_discrete_manual(colorant"red", colorant"deep sky blue"))

end

#TREND CHART
function trend_chrt(SimTimeArray, PeriodRange, quantiles)
#In order to stack vector entries, use the mapping function ****
AC_DF = vcat(map(x->x',SimTimeArray)...)
AC_DF = DataFrame(AC_DF)

trend_chart = []

for i = 1:size(AC_DF,2)
    push!(trend_chart, quantile!(AC_DF[i], quantiles))
end

#this allows us to join arrays generated by the quantile function
trend_chart = vcat(map(x->x',trend_chart)...)

#convert to time array to join dates and reconvert to DataFrame
trend_chart = TimeArray(PeriodRange, trend_chart)
trend_chart = DataFrame(trend_chart)

#Build Charts
names!(trend_chart, [:timestamp, :LowerBound, :p50, :UpperBound])
return plot(stack(trend_chart), y=:value, x=:timestamp, color=:variable, Geom.line)

end

#Error Bar Example

# x = sort(rand(10)); xmin = 0.9x; xmax = 1.1x;
# ystr = map(i -> string(i) * " as a string", 1:10)
# plot(y = ystr, x = ex, xmin = xmin, xmax = xmax, Geom.errorbar)


text:
"""
constants = get_constants(freq).\n
freq is the frequency of interest of your simulation.\n
contants include:
c  : lightspeed\n
μ₀ : permeability of vacuum\n
ϵ₀ : permettivity of vacuum\n
f₀ : frequency in Hz\n
ω  : frequency in rad/s\n
kb : wavenumber of vacuum\n
λ  : wavelength in vacuum\n
"""
function get_constants(freq)
    c  = 299792458 #apparently this is what sim4life uses
    μ₀ = 4f-7*π
    ϵ₀ = 1/(c^2*μ₀)
    f₀ = freq #round is for the fact that i set the 3T and 7T freq to 128 and 298 MHz in sim4life
    ω  = 2*(π*f₀) |> Float32
    kb = ω/c
    λ  = c/f₀
    return (;c, μ₀, ϵ₀, f₀, ω, kb, λ)
end

ext: jl


text:
using Celeste.SDSSIO
import Celeste.SDSSIO: SkyIntensity

@testset "sky interpolations" begin
    small_sky = [1.  2.  3.  4.;
                 5.  6.  7.  8.;
                 9. 10. 11. 12]
    sky_x = [0.1, 2.5]
    sky_y = [0.5, 2.5, 4.]
    sky = SkyIntensity(small_sky, sky_x, sky_y, ones(2))

    @test sky[1, 1] ≈ 1.0
    @test sky[2, 1] ≈ 7.0
    @test sky[1, 2] ≈ 2.5
    @test sky[2, 2] ≈ 8.5
    @test sky[1, 3] ≈ 4.0
    @test sky[2, 3] ≈ 10.0
end

@testset "test coordinates out of bounds" begin
    small_sky = [1.  2.  3.  4.;
            5.  6.  7.  8.;
            9. 10. 11. 12]
    sky_x = [-5.0, 4.0]
    sky_y = [-4.0, 5.0]
    sky = SkyIntensity(small_sky, sky_x, sky_y, ones(2))

    @test sky[1,1] == 1.0
    @test sky[1,2] == 4.0
    @test sky[2,1] == 9.0
    @test sky[2,2] == 12.0
end

@testset "read_photoobj handles missing table extensions" begin
    # get an example of a photoobj file with a missing table
    strategy = PlainFITSStrategy(joinpath(Pkg.dir("Celeste"), "test", "data"))
    rcf = RunCamcolField(6597, 4, 25)
    fname = SDSSIO.compute_fname(strategy, SDSSIO.PhotoObj(rcf))[1]
    if !isfile(fname)
        run(`curl --create-dirs -o $fname https://data.sdss.org/sas/dr12/boss/photoObj/301/6597/4/photoObj-006597-4-0025.fits`)
    end

    catalog = SDSSIO.read_photoobj(strategy, rcf)
end


text:
using SymPy
if VERSION >= v"0.5.0-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

@testset "Matrix" begin
    ## matrices
    (x,) = @syms x
    A = [x 1; 1 x]
    b = [x, 2]

    ## These fail for older installations of SymPy
    @test simplify(det(A)) == x^2 - 1

    ## we use inverse for A[:inv]()
    inv(A) # aliased to use inverse
    @test simplify((A[:inv]() - inv(A))[1,1]) == 0

    adjoint(A)
    dual(A)
    cholesky(A)
    ## other functions, could wrap
    b = subs(A, x, 2)
    QRdecomposition(b)

    # is_lower, is_square, is_symmetric much slower than julia only counterparts. May deprecate, but for now they are here
    @test is_lower(A) == istril(A)
    @test is_square(A) == true
    @compat @test is_symmetric(A) == issymmetric(A)
    @test Set(eigvals(A)) == Set([x-1, x+1])

    A = [1 0 0; 0 1 0; 0 0 x]
    evecs = eigvecs(A)
    @test evecs[1] == [1, 0, 0]

    ##
    M = Sym[1 2 3; 3 6 2; 2 0 1]
    q,r = QRdecomposition(M)
    @test (q * r - M)[1,1] == 0
    if VERSION >= v"0.5"
        L = [Sym[2,3,5], Sym[3,6,2], Sym[8,3,6]]
        out = GramSchmidt(L)
    end
    A = Sym[4 3; 6 3]
    L, U, _ = LUdecomposition(A)
    @test L == Sym[1 0; 3//2 1]

end


text:
using Antennas
using Documenter

DocMeta.setdocmeta!(Antennas, :DocTestSetup, :(using Antennas); recursive=true)

makedocs(;
    modules=[Antennas],
    authors="Kiran Shila <me@kiranshila.com> and contributors",
    repo="https://github.com/kiranshila/Antennas.jl/blob/{commit}{path}#{line}",
    sitename="Antennas.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://kiranshila.github.io/Antennas.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/kiranshila/Antennas.jl",
)


text:
using Documenter, ImageFeatures

makedocs(format   = :html,
         sitename = "ImageFeatures",
         pages    = ["Home" => "index.md",
                     "Tutorials" => [
                         "BRIEF" => "tutorials/brief.md",
                         "ORB" => "tutorials/orb.md",
                         "BRISK" => "tutorials/brisk.md",
                         "FREAK" => "tutorials/freak.md",
                         "Gray level co-occurence matrix" => "tutorials/glcm.md",
                         "Local binary patterns" => "tutorials/lbp.md",
                         "Object Detection using HOG" => "tutorials/object_detection.md"
                     ],
                     "Function reference" => "function_reference.md",
                     ],
         )

deploydocs(repo   = "github.com/JuliaImages/ImageFeatures.jl.git",
           julia  = "0.6",
           target = "build",
           deps   = nothing,
           make   = nothing,
           )


text:
using NeuralPDE
using Quadrature, Cubature, Cuba
using Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux
using Plots
using DelimitedFiles
using QuasiMonteCarlo
import ModelingToolkit: Interval, infimum, supremum

function nernst_planck(strategy, minimizer, maxIters)

    ##  DECLARATIONS
    @parameters t x y z
    @variables c(..)

    Dt = Differential(t)
    Dx = Differential(x)
    Dy = Differential(y)
    Dz = Differential(z)
    Dxx = Differential(x)^2
    Dyy = Differential(y)^2
    Dzz = Differential(z)^2

    ## DOMAINS AND OPERATORS

    # Discretization
    xwidth      = 1.0
    ywidth      = 1.0
    zwidth      = 1.0
    tmax        = 1.0
    xMeshNum    = 10
    yMeshNum    = 10
    zMeshNum    = 10
    tMeshNum    = 10

    dx  = xwidth/xMeshNum
    dy  = ywidth/yMeshNum
    dz  = zwidth/zMeshNum
    dt  = tmax/tMeshNum

    domains = [t ∈ Interval(0.0,tmax),
               x ∈ Interval(0.0,xwidth),
               y ∈ Interval(0.0,ywidth),
               z ∈ Interval(0.0,zwidth)]

    xs = 0.0 : dx : xwidth
    ys = 0.0 : dy : ywidth
    zs = 0.0 : dz : zwidth
    ts = 0.0 : dt : tmax

    # Constants
    D = 1  #dummy
    ux = 10 #dummy
    uy = 10 #dummy
    uz = 10 #dummy

    # Operators
    div = - D*(Dxx(c(t,x,y,z)) + Dyy(c(t,x,y,z)) + Dzz(c(t,x,y,z)))
          + (ux*Dx(c(t,x,y,z)) + uy*Dy(c(t,x,y,z)) + uz*Dz(c(t,x,y,z)))

    # Equation
    eq = Dt(c(t,x,y,z)) + div ~ 0      #NERNST-PLANCK EQUATION

    # Boundary conditions
    bcs = [c(0,x,y,z) ~ 0]

    ## NEURAL NETWORK
    n = 16   #neuron number

    chain = FastChain(FastDense(4,n,Flux.σ),FastDense(n,n,Flux.σ),FastDense(n,1))   #Neural network from Flux library

    indvars = [t,x,y,z]   #independent variables
    depvars = [c(t,x,y,z)]       #dependent (target) variable

    dim = length(domains)

    losses = []
    error = []
    times = []

    dx_err = 0.2

    error_strategy = GridTraining(dx_err)

    initθ = Float64.(DiffEqFlux.initial_params(chain))
    eltypeθ = eltype(initθ)
    parameterless_type_θ = DiffEqBase.parameterless_type(initθ)

    phi = NeuralPDE.get_phi(chain,parameterless_type_θ)
    derivative = NeuralPDE.get_numeric_derivative()

    _pde_loss_function = NeuralPDE.build_loss_function(eq,indvars,depvars,
                                             phi,derivative,nothing,chain,initθ,error_strategy)

    bc_indvars = NeuralPDE.get_variables(bcs,indvars,depvars)
    _bc_loss_functions = [NeuralPDE.build_loss_function(bc,indvars,depvars,
                                              phi,derivative,nothing,chain,initθ,error_strategy,
                                              bc_indvars = bc_indvar) for (bc,bc_indvar) in zip(bcs,bc_indvars)]

    train_sets = NeuralPDE.generate_training_sets(domains,dx_err,[eq],bcs,eltypeθ,indvars,depvars)
    train_domain_set, train_bound_set = train_sets

    pde_loss_functions = [NeuralPDE.get_loss_function(_pde_loss_function,
                                          train_domain_set[1],eltypeθ,
                                          parameterless_type_θ,error_strategy)]

    bc_loss_functions = [NeuralPDE.get_loss_function(bc,
                                          train_set,eltypeθ,
                                          parameterless_type_θ,error_strategy) for (train_set,bc) in zip(train_bound_set, _bc_loss_functions)]
    loss_functions =  [pde_loss_functions;bc_loss_functions]
    function loss_function_(θ,p)
        return sum(map(l->l(θ) ,loss_functions))
    end

    cb_ = function (p,l)
        deltaT_s = time_ns() #Start a clock when the callback begins, this will evaluate questo misurerà anche il calcolo degli uniform error

        ctime = time_ns() - startTime - timeCounter #This variable is the time to use for the time benchmark plot
        append!(times, ctime/10^9) #Conversion nanosec to seconds
        append!(losses, l)
        loss_ = loss_function_(p,nothing)
        append!(error, loss_)

        timeCounter = timeCounter + time_ns() - deltaT_s #timeCounter sums all delays due to the callback functions of the previous iterations

        return false
    end

    @named pde_system = PDESystem(eq, bcs, domains, indvars, depvars)

    discretization = NeuralPDE.PhysicsInformedNN(chain,strategy; init_params =initθ)
    prob = NeuralPDE.discretize(pde_system,discretization)

    timeCounter = 0.0
    startTime = time_ns() #Fix initial time (t=0) before starting the training
    res = GalacticOptim.solve(prob, minimizer, cb = cb_, maxiters=maxIters)

    phi = discretization.phi

    params = res.minimizer

    # Model prediction
    domain = [ts, xs, ys, zs]

    u_predict  = [reshape([phi([t,x,y,z],res.minimizer) for x in xs for y in ys for z in zs],
                 (length(xs),length(ys),length(zs))) for t in ts]

    return [error, params, domain, times]
end

maxIters = [(1,1,1,1000,1000,1000,1000),(1,1,1,300,300,300,300)] #iters for ADAM/LBFGS
# maxIters = [(1,1,1,10,10,10,10),(1,1,1,3,3,3,3)] #iters for ADAM/LBFGS

strategies = [NeuralPDE.QuadratureTraining(quadrature_alg = CubaCuhre(), reltol = 1e-4, abstol = 1e-4, maxiters = 50),
              NeuralPDE.QuadratureTraining(quadrature_alg = HCubatureJL(), reltol = 1e-4, abstol = 1e-4, maxiters = 50, batch = 0),
              NeuralPDE.QuadratureTraining(quadrature_alg = CubatureJLh(), reltol = 1e-4, abstol = 1e-4, maxiters = 50),
              NeuralPDE.QuadratureTraining(quadrature_alg = CubatureJLp(), reltol = 1e-4, abstol = 1e-4, maxiters = 50),
              NeuralPDE.GridTraining(0.2),
              NeuralPDE.StochasticTraining(400 ; bcs_points= 50),
              NeuralPDE.QuasiRandomTraining(400 ; bcs_points= 50)]

strategies_short_name = ["CubaCuhre",
                        "HCubatureJL",
                        "CubatureJLh",
                        "CubatureJLp",
                        "GridTraining",
                        "StochasticTraining",
                        "QuasiRandomTraining"]

minimizers = [ADAM(0.005),
              #BFGS()]
              LBFGS()]

minimizers_short_name = ["ADAM",
                         "LBFGS"]
                        # "BFGS"]

# Run models
error_res =  Dict()
domains = Dict()
params_res = Dict()  #to use same params for the next run
times = Dict()

## Convergence

for strat=1:length(strategies) # strategy
      for min =1:length(minimizers) # minimizer
            # println(string(strategies_short_name[strat], "  ", minimizers_short_name[min]))
            res = nernst_planck(strategies[strat], minimizers[min], maxIters[min][strat])
            push!(error_res, string(strat,min)     => res[1])
            push!(params_res, string(strat,min) => res[2])
            push!(domains, string(strat,min)        => res[3])
            push!(times, string(strat,min)        => res[4])
      end
end

#Plotting the first strategy with the first minimizer out from the loop to initialize the canvas
current_label = string(strategies_short_name[1], " + " , minimizers_short_name[1])
error = Plots.plot(times["11"], error_res["11"], yaxis=:log10, label = current_label)#, xlims = (0,10))#legend = true)#, size=(1200,700))
plot!(error, times["21"], error_res["21"], yaxis=:log10, label = string(strategies_short_name[2], " + " , minimizers_short_name[1]))
plot!(error, times["31"], error_res["31"], yaxis=:log10, label = string(strategies_short_name[3], " + " , minimizers_short_name[1]))
plot!(error, times["41"], error_res["41"], yaxis=:log10, label = string(strategies_short_name[4], " + " , minimizers_short_name[1]))
plot!(error, times["51"], error_res["51"], yaxis=:log10, label = string(strategies_short_name[5], " + " , minimizers_short_name[1]))
plot!(error, times["61"], error_res["61"], yaxis=:log10, label = string(strategies_short_name[6], " + " , minimizers_short_name[1]))
plot!(error, times["71"], error_res["71"], yaxis=:log10, label = string(strategies_short_name[7], " + " , minimizers_short_name[1]))

plot!(error, times["12"], error_res["12"], yaxis=:log10, label = string(strategies_short_name[1], " + " , minimizers_short_name[2]))
plot!(error, times["22"], error_res["22"], yaxis=:log10, label = string(strategies_short_name[2], " + " , minimizers_short_name[2]))
plot!(error, times["32"], error_res["32"], yaxis=:log10, label = string(strategies_short_name[3], " + " , minimizers_short_name[2]))
plot!(error, times["42"], error_res["42"], yaxis=:log10, label = string(strategies_short_name[4], " + " , minimizers_short_name[2]))
plot!(error, times["52"], error_res["52"], yaxis=:log10, label = string(strategies_short_name[5], " + " , minimizers_short_name[2]))
plot!(error, times["62"], error_res["62"], yaxis=:log10, label = string(strategies_short_name[6], " + " , minimizers_short_name[2]))
plot!(error, times["72"], error_res["72"], yaxis=:log10, title = string("Nernst Planck convergence ADAM/LBFGS"), ylabel = "log(error)", xlabel = "t",label = string(strategies_short_name[7], " + " , minimizers_short_name[2]))

using SciMLBenchmarks
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

max_stars_repo_name: SciML/SciMLBenchmarksOutput

max_issues_repo_name: SciML/SciMLBenchmarksOutput

max_forks_repo_name: SciML/SciMLBenchmarksOutput


text:
"""
Subroutine dspcg

This subroutine generates a sequence of approximate minimizers
for the subproblem

    min { q(x) : xl <= x <= xu }.

The quadratic is defined by

    q(x[0]+s) = 0.5*s'*A*s + g'*s,

where x[0] is a base point provided by the user, A is a symmetric
matrix in compressed column storage, and g is a vector.

At each stage we have an approximate minimizer x[k], and generate
a direction p[k] by using a preconditioned conjugate gradient
method on the subproblem

    min { q(x[k]+p) : || L'*p || <= delta, s(fixed) = 0 },

where fixed is the set of variables fixed at x[k], delta is the
trust region bound, and L is an incomplete Cholesky factorization
of the submatrix

    B = A(free:free),

where free is the set of free variables at x[k]. Given p[k],
the next minimizer x[k+1] is generated by a projected search.

The starting point for this subroutine is x[1] = x[0] + s, where
x[0] is a base point and s is the Cauchy step.

The subroutine converges when the step s satisfies

    || (g + A*s)[free] || <= rtol*|| g[free] ||

In this case the final x is an approximate minimizer in the face
defined by the free variables.

The subroutine terminates when the trust region bound does
not allow further progress, that is, || L'*p[k] || = delta.
In this case the final x satisfies q(x) < q(x[k]).

MINPACK-2 Project. March 1999.
Argonne National Laboratory.
Chih-Jen Lin and Jorge J. More'.

March 2000

Clarified documentation of nv variable.
Eliminated the nnz = max(nnz,1) statement.
"""
function dspcg(n,x,xl,xu,A,g,delta,
               rtol,s,nv,itermax,
               B, L,
               indfree,gfree,w,wa,iwa)
    zero = 0.0
    one = 1.0

    # Compute A*(x[1] - x[0]) and store in w.

    dssyax(n, A, s, w)

    # Compute the Cauchy point.

    daxpy(n,one,s,1,x,1)
    dmid(n,x,xl,xu)

    # Start the main iteration loop.
    # There are at most n iterations because at each iteration
    # at least one variable becomes active.

    info = 3
    iters = 0
    for nfaces=1:n

        # Determine the free variables at the current minimizer.
        # The indices of the free variables are stored in the first
        # n free positions of the array indfree.
        # The array iwa is used to detect free variables by setting
        # iwa[i] = nfree if the ith variable is free, otherwise iwa[i] = 0.

        nfree = 0
        # TODO
        @inbounds for j=1:n
            if xl[j] < x[j] && x[j] < xu[j]
                nfree = nfree + 1
                indfree[nfree] = j
                iwa[j] = nfree
            else
                iwa[j] = 0
            end
        end

        # Exit if there are no free constraints.

        if nfree == 0
            info = 1
            return info, iters
        end

        # Obtain the submatrix of A for the free variables.
        # Recall that iwa allows the detection of free variables.
        nnz = reorder!(B, A, indfree, nfree, iwa)

        # Compute the incomplete Cholesky factorization.
        alpha = zero
        # TODO
        dicfs(nfree, nnz, B, L,
              nv, alpha,
              iwa, view(wa,1:n), view(wa,n+1:2*n))

        # Compute the gradient grad q(x[k]) = g + A*(x[k] - x[0]),
        # of q at x[k] for the free variables.
        # Recall that w contains A*(x[k] - x[0]).
        # Compute the norm of the reduced gradient Z'*g.

        # TODO
        @inbounds for j=1:nfree
            gfree[j] = w[indfree[j]] + g[indfree[j]]
            wa[j] = g[indfree[j]]
        end
        gfnorm = dnrm2(nfree,wa,1)

        # Save the trust region subproblem in the free variables
        # to generate a direction p[k]. Store p[k] in the array w.

        tol = rtol*gfnorm
        stol = zero

        infotr,itertr = dtrpcg(nfree,B,gfree,delta, L,
                               tol,stol,itermax,w,
                               view(wa,1:n),view(wa,n+1:2*n),view(wa,2*n+1:3*n),
                               view(wa,3*n+1:4*n),view(wa,4*n+1:5*n))

        iters = iters + itertr
        dtsol(nfree, L, w)

        # Use a projected search to obtain the next iterate.
        # The projected search algorithm stores s[k] in w.

        # TODO
        @inbounds for j=1:nfree
            wa[j] = x[indfree[j]]
            wa[n+j] = xl[indfree[j]]
            wa[2*n+j] = xu[indfree[j]]
        end

        dprsrch(nfree,view(wa,1:n),view(wa,n+1:2*n),view(wa,2*n+1:3*n),
                B,gfree,w,
                view(wa,3*n+1:4*n), view(wa,4*n+1:5*n))

        # Update the minimizer and the step.
        # Note that s now contains x[k+1] - x[0].

        # TODO
        @inbounds for j=1:nfree
            x[indfree[j]] = wa[j]
            s[indfree[j]] = s[indfree[j]] + w[j]
        end

        # Compute A*(x[k+1] - x[0]) and store in w.

        dssyax(n, A, s, w)

        # Compute the gradient grad q(x[k+1]) = g + A*(x[k+1] - x[0])
        # of q at x[k+1] for the free variables.

        # TODO
        @inbounds for j=1:nfree
            gfree[j] = w[indfree[j]] + g[indfree[j]]
        end
        gfnormf = dnrm2(nfree, gfree, 1)

        # Convergence and termination test.
        # We terminate if the preconditioned conjugate gradient
        # method encounters a direction of negative curvature, or
        # if the step is at the trust region bound.

        if gfnormf <= rtol*gfnorm
            info = 1
            return info, iters
        elseif infotr == 3 || infotr == 4
            info = 2
            return info, iters
        elseif iters > itermax
            info = 3
            return info, iters
        end
    end

    return info, iters
end

@inline function dspcg(n::Int, delta::Float64, rtol::Float64, itermax::Int,
               x::CuDeviceArray{Float64,1}, xl::CuDeviceArray{Float64,1},
               xu::CuDeviceArray{Float64,1}, A::CuDeviceArray{Float64,2},
               g::CuDeviceArray{Float64,1}, s::CuDeviceArray{Float64,1},
               B::CuDeviceArray{Float64,2}, L::CuDeviceArray{Float64,2},
               indfree::CuDeviceArray{Int,1}, gfree::CuDeviceArray{Float64,1},
               w::CuDeviceArray{Float64,1}, iwa::CuDeviceArray{Int,1},
               wa1::CuDeviceArray{Float64,1}, wa2::CuDeviceArray{Float64,1},
               wa3::CuDeviceArray{Float64,1}, wa4::CuDeviceArray{Float64,1},
               wa5::CuDeviceArray{Float64,1})

    tx = threadIdx().x
    ty = threadIdx().y

    zero = 0.0
    one = 1.0

    # Compute A*(x[1] - x[0]) and store in w.

    dssyax(n, A, s, w)

    # Compute the Cauchy point.

    daxpy(n,one,s,1,x,1)
    dmid(n,x,xl,xu)

    # Start the main iteration loop.
    # There are at most n iterations because at each iteration
    # at least one variable becomes active.

    info = 3
    iters = 0
    for nfaces=1:n

        # Determine the free variables at the current minimizer.
        # The indices of the free variables are stored in the first
        # n free positions of the array indfree.
        # The array iwa is used to detect free variables by setting
        # iwa[i] = nfree if the ith variable is free, otherwise iwa[i] = 0.

        # Use a single thread to avoid multiple branch divergences.
        # XXX: Would there be any gain in employing multiple threads?
        nfree = 0
        if tx == 1 && ty == 1
            @inbounds for j=1:n
                if xl[j] < x[j] && x[j] < xu[j]
                    nfree = nfree + 1
                    indfree[nfree] = j
                    iwa[j] = nfree
                else
                    iwa[j] = 0
                end
            end
        end
        nfree = CUDA.shfl_sync(0xffffffff, nfree, 1)

        # Exit if there are no free constraints.

        if nfree == 0
            info = 1
            return info, iters
        end

        # Obtain the submatrix of A for the free variables.
        # Recall that iwa allows the detection of free variables.
        reorder!(n, nfree, B, A, indfree, iwa)

        # Compute the incomplete Cholesky factorization.
        alpha = zero
        dicfs(nfree, alpha, B, L, wa1, wa2)

        # Compute the gradient grad q(x[k]) = g + A*(x[k] - x[0]),
        # of q at x[k] for the free variables.
        # Recall that w contains A*(x[k] - x[0]).
        # Compute the norm of the reduced gradient Z'*g.

        if tx <= nfree && ty == 1
            @inbounds begin
                gfree[tx] = w[indfree[tx]] + g[indfree[tx]]
                wa1[tx] = g[indfree[tx]]
            end
        end
        CUDA.sync_threads()
        gfnorm = dnrm2(nfree,wa1,1)

        # Save the trust region subproblem in the free variables
        # to generate a direction p[k]. Store p[k] in the array w.

        tol = rtol*gfnorm
        stol = zero

        infotr,itertr = dtrpcg(nfree,B,gfree,delta,L,
                               tol,stol,itermax,w,
                               wa1,wa2,wa3,wa4,wa5)

        iters += itertr
        dtsol(nfree, L, w)

        # Use a projected search to obtain the next iterate.
        # The projected search algorithm stores s[k] in w.

        if tx <= nfree && ty == 1
            @inbounds begin
                wa1[tx] = x[indfree[tx]]
                wa2[tx] = xl[indfree[tx]]
                wa3[tx] = xu[indfree[tx]]
            end
        end
        CUDA.sync_threads()

        dprsrch(nfree,wa1,wa2,wa3,B,gfree,w,wa4,wa5)

        # Update the minimizer and the step.
        # Note that s now contains x[k+1] - x[0].

        if tx <= nfree && ty == 1
            @inbounds begin
                x[indfree[tx]] = wa1[tx]
                s[indfree[tx]] += w[tx]
            end
        end
        CUDA.sync_threads()

        # Compute A*(x[k+1] - x[0]) and store in w.

        dssyax(n, A, s, w)

        # Compute the gradient grad q(x[k+1]) = g + A*(x[k+1] - x[0])
        # of q at x[k+1] for the free variables.

        if tx == 1 && ty == 1
            @inbounds for j=1:nfree
                gfree[j] = w[indfree[j]] + g[indfree[j]]
            end
        end
        CUDA.sync_threads()

        gfnormf = dnrm2(nfree, gfree, 1)

        # Convergence and termination test.
        # We terminate if the preconditioned conjugate gradient
        # method encounters a direction of negative curvature, or
        # if the step is at the trust region bound.

        if gfnormf <= rtol*gfnorm
            info = 1
            return info, iters
        elseif infotr == 3 || infotr == 4
            info = 2
            return info, iters
        elseif iters > itermax
            info = 3
            return info, iters
        end
    end

    return info, iters
end


text:
module VDPTag2Experiments

import ARDESPOT

using POMDPs
using MCTS
using ParticleFilters
using StaticArrays
using VDPTag2

using POMDPModelTools: GenerativeBeliefMDP
using Random: Random, MersenneTwister, AbstractRNG, randperm

export
    RootToNextMLFirst,
    ObsAdaptiveParticleFilter,
    VDPUpper

struct RootToNextMLFirst
    rng::MersenneTwister
end

function MCTS.next_action(gen::RootToNextMLFirst, p::VDPTagPOMDP, b, node)
    if isroot(node) && n_children(node) < 1
        target_sum=MVector(0.0, 0.0)
        agent_sum=MVector(0.0, 0.0)
        for s in particles(b::ParticleCollection)
            target_sum += s.target
            agent_sum += s.agent
        end
        next = VDPTag2.next_ml_target(mdp(p), target_sum/n_particles(b))
        diff = next-agent_sum/n_particles(b)
        return TagAction(false, atan(diff[2], diff[1]))
    else
        return rand(gen.rng, actions(p))
    end
end

MCTS.next_action(gen::RootToNextMLFirst, p::GenerativeBeliefMDP, b, node) = next_action(gen, p.pomdp, b, node)

struct ObsAdaptiveParticleFilter{P<:POMDP,S,R,RNG<:AbstractRNG} <: Updater
    pomdp::P
    resample::R
    max_frac_replaced::Float64
    n_init::Int
    rng::RNG
    _pm::Vector{S}
    _wm::Vector{Float64}
end

function ObsAdaptiveParticleFilter(p::POMDP, resample, max_frac_replaced, n_init, rng::AbstractRNG)
    S = statetype(p)
    return ObsAdaptiveParticleFilter(p, resample, max_frac_replaced, n_init, rng, S[], Float64[])
end

POMDPs.initialize_belief(up::ObsAdaptiveParticleFilter, d::Any) = ParticleCollection([rand(up.rng, d) for i in 1:up.n_init])
POMDPs.update(up::ObsAdaptiveParticleFilter, b, a, o) = update(up, resample(up.resample, b, up.rng), a, o)

function POMDPs.update(up::ObsAdaptiveParticleFilter, b::ParticleCollection, a, o)
    if n_particles(b) > 2*up.resample.n
        b = resample(up.resample, b, up.rng)
    end

    ps = particles(b)
    pm = up._pm
    wm = up._wm
    resize!(pm, 0)
    resize!(wm, 0)

    all_terminal = true
    for i in 1:n_particles(b)
        s = ps[i]
        if !isterminal(up.pomdp, s)
            all_terminal = false
            sp = @gen(:sp)(up.pomdp, s, a, up.rng)
            push!(pm, sp)
            od = observation(up.pomdp, s, a, sp)
            push!(wm, pdf(od, o))
        end
    end
    ws = sum(wm)
    if all_terminal || ws < eps(1.0/length(wm))
        # warn("All states in particle collection were terminal.")
        return initialize_belief(up, reset_distribution(up.pomdp, b, a, o))
    end

    pc = resample(up.resample, WeightedParticleBelief{statetype(up.pomdp)}(pm, wm, ws, nothing), up.rng)
    ps = particles(pc)

    mpw = max_possible_weight(up.pomdp, a, o)
    frac_replaced = up.max_frac_replaced*max(0.0, 1.0 - maximum(wm)/mpw)
    n_replaced = floor(Int, frac_replaced*length(ps))
    is = randperm(up.rng, length(ps))[1:n_replaced]
    for i in is
        ps[i] = new_particle(up.pomdp, b, a, o, up.rng)
    end
    return pc
end

max_possible_weight(pomdp::VDPTagPOMDP, a::TagAction, o) = 0.0
new_particle(pomdp::VDPTagPOMDP, a::TagAction, o) = error("shouldn't get here")
reset_distribution(p::POMDP, b, a, o) = initialstate(p)

struct VDPUpper end

function ARDESPOT.ubound(ub::VDPUpper, pomdp::POMDP, b::ARDESPOT.ScenarioBelief)
    if all(isterminal(pomdp, s) for s in particles(b))
        return 0.0
    else
        return mdp(cproblem(pomdp)).tag_reward
    end
end

Random.seed!(p::Policy, ::Any) = nothing

end


text:
a = rand(Float32, 10, 10)
ad = AFArray(a)

# Issue https://github.com/gaika/ArrayFire.jl/issues/26
@test im * ad == complex(zeros(ad), ad)

# Indexing - issue #96
let
    b = Float32[1.,2.,3.]
    bd = AFArray(b)
    ind = AFArray([false, true, true])
    @test Array(bd[ind]) == Float32[2., 3.]
    ind[3] = false
    @test Array(bd[ind]) == Float32[2.]
end

# Broadcast
if VERSION >= v"0.5.0"
    @test sin(ad) == sin.(ad)
end

# Sign - issue #109
if VERSION >= v"0.5.0"
    let
        a = randn(Float32, 10)
        ad = AFArray(a)
        @test Array(signbit.(ad)) == signbit.(a)
        @test Array(sign.(ad)) == sign.(a)
    end
end

# Indexing - issue #115
let
    x = rand(Float32, 3,3)
    xd = AFArray(x)
    y = x[:, [1,3]]
    yd = xd[:, [1,3]]
    @test x ≈ Array(xd)
end

# Return types of sum reduction
let
    for T in (Float32, Complex{Float32}, Int32,
              Int64, UInt32, UInt8, UInt64, Bool)
        s = AFArray(rand(T, 10))
        if T == UInt8
            @test typeof(sum(s)) == UInt32
        elseif T == Bool
            @test typeof(sum(s)) == Int64
        else
            @test typeof(sum(s)) == T
        end
    end
end

# Issue #125
let
    a = rand(Float32, 5, 5)
    ad = AFArray(a)
    if VERSION >= v"0.6.0"
        for f in (>, >=, <, <=, ==)
            for val in (0.5f0, 1)
                b = broadcast(f, val, a)
                bd = broadcast(f, val, ad)
                @test b ≈ Array(bd)
                b = broadcast(f, a, val)
                bd = broadcast(f, ad, val)
                @test b ≈ Array(bd)
            end
        end
    else
        for f in (.>, .>=, .<, .<=, .==)
            for val in (0.5f0, 1)
                b = f(val, a)
                bd = f(val, ad)
                @test b ≈ Array(bd)
                b = f(a, val)
                bd = f(ad, val)
                @test b ≈ Array(bd)
            end
        end
    end
end

# Issue #131
let
    set_device(0)
    for i = 1:10
        a = rand(AFArray{Float32}, 10)
        a + 1
    end
    sync(0)
end

# Issue #133
let
    for sz in ((10,), (10, 10), (10, 10, 10))
        a = rand(Complex{Float32}, sz...)
        ad = AFArray(a)
        fft!(a)
        fft!(ad)
        @test a ≈ Array(ad)
    end
end

# Issue #143
let
    AD = AFArray(ones(Float32,3,3)+im*ones(Float32,3,3))
    AH = Array(AD)
    @test Array(abs(AD)) ≈ abs.(AH)
end

# Issue #103
let
    a = rand(AFArray{ComplexF64}, 2, 2)
    a[1,1]  = Float32(1)+Float32(2)im
    @test a[1,1] == 1.0f0 + 2.0f0im
end


text:
## version range

struct VersionRange
    lower::VersionNumber
    upper::VersionNumber
end

Base.in(v::VersionNumber, r::VersionRange) = (v >= r.lower && v <= r.upper)

import Base.(:)
(:)(a::VersionNumber, b::VersionNumber) = VersionRange(a, b)

Base.intersect(v::VersionNumber, r::VersionRange) =
    v < r.lower ? (r.lower:v) :
    v > r.upper ? (v:r.upper) : (v:v)


text:
using LinearAlgebra, OrdinaryDiffEq, Test
f = (du,u,p,t) -> du .= u ./ t
jac = (J,u,p,t) -> (J[1,1] = 1/t; J[2,2] = 1/t; J)

jp_diag = Diagonal(zeros(2))
fun = ODEFunction(f; jac=jac, jac_prototype=jp_diag)
prob = ODEProblem(fun,ones(2),(1.0,10.0))
sol = solve(prob,Rosenbrock23(autodiff=false))
@test sol[end] ≈ [10.0,10.0]
@test length(sol) < 60

jp = Tridiagonal(jp_diag)
fun = ODEFunction(f; jac=jac, jac_prototype=jp)
prob = ODEProblem(fun,ones(2),(1.0,10.0))
sol = solve(prob,Rosenbrock23(autodiff=false))
@test sol[end] ≈ [10.0,10.0]
@test length(sol) < 60

jp = SymTridiagonal(jp_diag)
fun = ODEFunction(f; jac=jac, jac_prototype=jp)
prob = ODEProblem(fun,ones(2),(1.0,10.0))
sol = solve(prob,Rosenbrock23(autodiff=false))
@test sol[end] ≈ [10.0,10.0]
@test length(sol) < 60

@test_broken begin
    jp = Hermitian(jp_diag)
    fun = ODEFunction(f; jac=jac, jac_prototype=jp)
    prob = ODEProblem(fun,ones(2),(1.0,10.0))
    sol = solve(prob,Rosenbrock23(autodiff=false))
    @test sol[end] ≈ [10.0,10.0]
    @test length(sol) < 60
end

@test_broken begin
    jp = Symmetric(jp_diag)
    fun = ODEFunction(f; jac=jac, jac_prototype=jp)
    prob = ODEProblem(fun,ones(2),(1.0,10.0))
    sol = solve(prob,Rosenbrock23(autodiff=false))
    @test sol[end] ≈ [10.0,10.0]
    @test length(sol) < 60
end


text:
function OC_doDeepOcnCorrectionOfT!(
    ocn :: Ocean,
    i   :: Integer,
    j   :: Integer;
    Δt  :: Float64,
    τ   :: Float64 = ocn.Ts_clim_relax_time,
    start_layer :: Integer,
)

    return doDeepOcnCorrection!(
        h_ML    = ocn.h_ML[i, j],
        qs      = ocn.cols.Ts[i, j],
        qs_clim = ocn.cols.Ts_clim[i, j],
        FLDO    = ocn.FLDO[i, j],
        Nz      = ocn.Nz[i, j],
        τ       = τ,
        Δt      = Δt,
        hs      = ocn.cols.hs[i, j],
        zs      = ocn.cols.zs[i, j],
        start_layer = start_layer,
    )

end

function OC_doDeepOcnCorrectionOfS!(
    ocn :: Ocean,
    i   :: Integer,
    j   :: Integer;
    Δt  :: Float64,
    start_layer :: Integer,
)

    return doDeepOcnCorrection!(
        h_ML    = ocn.h_ML[i, j],
        qs      = ocn.cols.Ss[i, j],
        qs_clim = ocn.cols.Ss_clim[i, j],
        FLDO    = ocn.FLDO[i, j],
        Nz      = ocn.Nz[i, j],
        τ       = ocn.Ss_clim_relax_time,
        Δt      = Δt,
        hs      = ocn.cols.hs[i, j],
        zs      = ocn.cols.zs[i, j],
        start_layer = start_layer,
    )

end

"""

    This function newtonian-relaxes `qs` to `qs_clim` with e-folding time `τ`.
    It has to be noted that mixed-layer is not relaxed.

    Also, Euler backward integration scheme is used.

"""
function doDeepOcnCorrection!(;
    h_ML       :: Float64,
    qs         :: AbstractArray{Float64, 1},
    qs_clim    :: AbstractArray{Float64, 1},
    FLDO       :: Integer,
    Nz         :: Integer,
    τ          :: Float64,
    Δt         :: Float64,
    hs         :: AbstractArray{Float64, 1},
    zs         :: AbstractArray{Float64, 1},
    start_layer:: Integer,
)

    src_and_sink = 0.0

    if τ > 0.0

        r = Δt / τ
        if FLDO != -1 && start_layer <= Nz

            for i = max(FLDO, start_layer):Nz
                dq = r * (qs_clim[i] - qs[i]) / (1.0 + r)
                src_and_sink += ((i == FLDO) ? - zs[FLDO+1] - h_ML : hs[i]) * dq
                #qs[i] = (qs[i] + r * qs_clim[i]) / (1+r)
                qs[i] += dq
            end

        end

    elseif τ == 0.0

        if FLDO != -1 && start_layer <= Nz
            for i = max(FLDO, start_layer):Nz
                src_and_sink += ((i == FLDO) ? - zs[FLDO+1] - h_ML : hs[i]) * (qs_clim[i] - qs[i])
                qs[i] = qs_clim[i]
            end
        end

    end

    return src_and_sink / Δt
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using Sobol, Random, Distributions

export ParticleSampler

"""
    ParticleSampler{D,V}( sampling_type, symmetric, dims, n_particles)

Particle initializer class with various functions to initialize a particle.

- `sampling_type` : `:random` or `:sobol`
- `symmetric` : `true` or `false`
- `n_particles` : number of particles
"""
struct ParticleSampler{D,V}
    sampling_type::Symbol
    dims::Tuple{Int,Int}
    n_particles::Int
    symmetric::Bool
    seed::Int

    function ParticleSampler{D,V}(
        sampling_type::Symbol, symmetric::Bool, n_particles::Int, seed::Int=1234
    ) where {D,V}
        if !(sampling_type in [:random, :sobol])
            throw(ArgumentError("Sampling type $sampling_type
                      not implemented"))
        end

        dims = (D, V)

        # Make sure that the particle number is conforming
        # with symmetric sampling (if necessary)
        if symmetric
            ncopies = 2^(sum(dims))
            np = mod(n_particles, ncopies)
            if np != 0
                n_particles += np
            end
        else
            ncopies = 1
        end

        return new(sampling_type, dims, n_particles, symmetric, seed)
    end

    function ParticleSampler{D,V}(
        sampling_type::Symbol, n_particles::Int, seed::Int=1234
    ) where {D,V}
        if !(sampling_type in [:random, :sobol])
            throw(ArgumentError("Sampling type $sampling_type
                      not implemented"))
        end

        dims = (D, V)
        symmetric = false

        return new(sampling_type, dims, n_particles, symmetric, seed)
    end
end

export sample!

"""
    sample!( pg, ps, df, mesh)

Sample from a Particle sampler

- `pg`   : Particle group
- `ps`   : Particle sampler
- `df`   : Distribution function
- `xmin` : lower bound of the domain
- `dimx` : length of the domain.
"""
function sample!(
    pg::ParticleGroup{1,2}, ps::ParticleSampler, df::AbstractCosGaussian, mesh::AbstractGrid
)
    if ps.symmetric
        sample_sym(ps, pg, df, mesh)
    else
        sample_all(ps, pg, df, mesh)
    end
end

"""
    sample_all( ps, pg, df, mesh )

Helper function for pure sampling
"""
function sample_all(ps, pg::ParticleGroup{1,2}, df::AbstractCosGaussian, mesh)
    ndx, ndv = df.dims

    x = zeros(ndx)
    v = zeros(ndv)

    n_rnds = 0
    if df.params.n_gaussians > 1
        n_rnds = 1
    end

    δ = zeros(df.params.n_gaussians)
    for i_v in 1:(df.params.n_gaussians)
        δ[i_v] = sum(df.params.δ[1:i_v])
    end

    n_rnds += ndx + ndv
    rdn = zeros(ndx + ndv + 1)

    if ps.sampling_type == :sobol
        rng_sobol = SobolSeq(ndx)
    else
        rng_random = MersenneTwister(ps.seed)
    end

    d = Normal()

    for i_part in 1:(pg.n_particles)
        if ps.sampling_type == :sobol
            x .= mesh.xmin .+ Sobol.next!(rng_sobol) * mesh.dimx
        else
            x .= mesh.xmin .+ rand(rng_random, ndx) * mesh.dimx
        end

        # Set weight according to value of perturbation
        w = eval_x_density(df, x) .* mesh.dimx

        v .= rand!(d, v)

        # For multiple Gaussian, draw which one to take
        rnd_no = rdn[ndx + ndv + 1]
        i_gauss = 1
        while (rnd_no > δ[i_gauss])
            i_gauss += 1
        end
        v .= v .* df.params.σ[i_gauss] .+ df.params.μ[i_gauss]

        # Copy the generated numbers to the particle
        set_x!(pg, i_part, x)
        set_v!(pg, i_part, v)

        # Set weights.
        set_weights!(pg, i_part, w)
    end
end

"""
    sample_sym( ps, pg, df, mesh )

Helper function for antithetic sampling in 1d2v
"""
function sample_sym(ps, pg, df, mesh)
    ndx, ndv = df.dims

    x = zeros(Float64, ndx)
    v = zeros(Float64, ndv)

    n_rnds = 0
    (df.params.n_gaussians > 1) && (n_rnds = 1)

    δ = zeros(Float64, df.params.n_gaussians)
    for i_v in 1:(df.params.n_gaussians)
        δ[i_v] = sum(df.params.δ[1:i_v])
    end

    n_rnds = n_rnds + ndx + ndv
    rdn = zeros(ndx + ndv + 1)

    if ps.sampling_type == :sobol
        rng_sobol = SobolSeq(ndx + ndv + 1)
    else
        rng_random = MersenneTwister(ps.seed)
    end

    dnormal = Normal()

    i_gauss = 1::Int
    wi = 0.0::Float64

    for i_part in 1:(pg.n_particles)
        ip = i_part % 8

        if ip == 1

            # Generate Random or Sobol numbers on [0,1]
            if ps.sampling_type == :sobol
                rdn .= Sobol.next!(rng_sobol)
            else
                rdn .= rand!(rng_random, rdn)
            end

            # Transform rdn to the interval
            x[1:ndx] .= mesh.xmin .+ mesh.dimx .* rdn[1:ndx]

            # Set weight according to value of perturbation
            wi = eval_x_density(df, x[1:ndx]) * mesh.dimx

            # Maxwellian distribution of the temperature

            v .= rand!(dnormal, v)

            # For multiple Gaussian, draw which one to take
            rnd_no = rdn[ndx + ndv + 1]
            i_gauss = 1
            while i_gauss < df.params.n_gaussians && rnd_no > δ[i_gauss]
                i_gauss += 1
            end

            v .= v .* df.params.σ[i_gauss] .+ df.params.μ[i_gauss]

        elseif ip == 5
            x[1] = mesh.dimx - x[1] + 2.0 * mesh.xmin

        elseif ip % 2 == 0
            v[1] = -v[1] + 2.0 * df.params.μ[i_gauss][1]

        else
            v[2] = -v[2] + 2.0 * df.params.μ[i_gauss][2]
        end

        # Copy the generated numbers to the particle

        set_x!(pg, i_part, x)
        set_v!(pg, i_part, v)
        set_weights!(pg, i_part, wi)
    end
end


text:
"""
    vertices(g)

Returns an iterator to the vertices of a graph (i.e. 1:nv(g))
"""
vertices(g::AGraphOrDiGraph) = one(vertextype(g)):nv(g)

"""
    adjacency_list(g::AGraph)
    adjacency_list(g::ADiGraph, dir=:out)

Returns the adjacency list `a` of a graph (a vector of vector of ints). The `i`-th
element of the adjacency list is a vector containing the neighbors of `i` in `g`.

For directed graphs a second optional argument can be specified (`:out` or `:in`).
The neighbors in the returned adjacency list are considered accordingly as those
related through outgoing or incoming edges.

The elements of  `a[i]` have the same order as in the iterator `(out_/in_)neighbors(g,i)`.

*Attention*: For some graph types it returns a reference, not a copy,
therefore the returned object should not be modified.
"""
adjacency_list(g::AGraph) = out_adjlist(g)
adjacency_list(g::ADiGraph, dir=:out) = dir == :out ? out_adjlist(g) :
                                        dir == :in ? in_adjlist(g) :
                                        error("Second argument has to be :out or :in.")

"""
    add_vertices!(g, n)

Add `n` new vertices to the graph `g`. Returns the final number
of vertices.
"""
function add_vertices!(g::AGraphOrDiGraph, n)
    added = true
    for i = 1:n
        add_vertex!(g)
    end
    return nv(g)
end

"""
    has_vertex(g, v)

Return true if `v` is a vertex of `g`.
"""
has_vertex(g::AGraphOrDiGraph, v) = v in vertices(g)

show(io::IO, g::AGraphOrDiGraph) = shortshow(io, g)
show(io::IO, g::ANetOrDiNet) = longshow(io, g)

function shortshow(io::IO, g::G) where G<:AGraphOrDiGraph
    print(io, split("$G",'.')[end],
        "($(nv(g)), $(ne(g)))")
        # is_directed(g) ? " undirected graph" : " directed graph")
end

function longshow(io::IO, g::G) where G<:ANetOrDiNet
    print(io, split("$G",'.')[end], "($(nv(g)), $(ne(g)))")
    print(io, " with ")
    _printstrvec(io, gprop_names(g))
    print(io, " graph, ")
    _printstrvec(io, vprop_names(g))
    print(io," vertex, ")
    _printstrvec(io, eprop_names(g))
    print(io, " edge properties")
end

function _printstrvec(io::IO, vs::Vector{String})
    print(io,"[")
    if length(vs) > 0
        for s in vs[1:end-1]
            print(io, "\"" * s * "\", ")
        end
        print(io, "\"" * last(vs) * "\"")
    end
    print(io,"]")
end

"""
    is_directed(g)

Check if `g` a graph with directed edges.
"""
is_directed(g::AGraph) = false
is_directed(g::ADiGraph) = true
is_directed(::Type{G}) where {G<:AGraph} = false
is_directed(::Type{G}) where {G<:ADiGraph} = true

"""
    in_degree(g, v)

Returns the number of edges which start at vertex `v`.
"""
in_degree(g::AGraphOrDiGraph, v::Integer) = length(in_neighbors(g,v))

"""
    out_degree(g, v)

Returns the number of edges which end at vertex `v`.
"""
out_degree(g::AGraphOrDiGraph, v::Integer) = length(out_neighbors(g,v))

"""
    degree(g, v)

Return the number of edges  from the vertex `v`.
"""
degree(g::AGraphOrDiGraph, v::Integer) = out_degree(g, v)

in_degree(g::AGraphOrDiGraph, v::AbstractVector{<:Integer} = vertices(g)) = [in_degree(g,x) for x in v]
out_degree(g::AGraphOrDiGraph, v::AbstractVector{<:Integer} = vertices(g)) = [out_degree(g,x) for x in v]
degree(g::AGraphOrDiGraph, v::AbstractVector{<:Integer} = vertices(g)) = [degree(g,x) for x in v]

"""
    neighbors(g, v)

Returns a list of all neighbors from vertex `v` in `g`.

For directed graph, this is equivalent to [`out_neighbors`](@ref)(g, v).

NOTE: it may return a reference, not a copy. Do not modify result.
"""
neighbors(g::AGraphOrDiGraph, v::Integer) = out_neighbors(g, v)
in_neighbors(g::AGraph, v::Integer) = out_neighbors(g, v)

"""
    all_neighbors(g, v)

Iterates over all distinct in/out neighbors of vertex `v` in `g`.
"""
all_neighbors(g::AGraph, v::Integer) = out_neighbors(g, v)

all_neighbors(g::ADiGraph, v::Integer) =
    distinct(flatten((out_neighbors(g, v), in_neighbors(g, v))))

"""
    density(g)

Density is defined as the ratio of the number of actual edges to the
number of possible edges. This is ``|v| |v-1|`` for directed graphs and
``(|v| |v-1|) / 2`` for undirected graphs.
"""
density(g::AGraph) = (2*ne(g)) / (nv(g) * (nv(g)-1))
density(g::ADiGraph) = ne(g) / (nv(g) * (nv(g)-1))

"""
    clean_vertex!(g, v)

Remove all incident edges on vertex `v` in `g`.
"""
function clean_vertex!(g::AGraphOrDiGraph, v::Integer)
    edgs = collect(all_edges(g, v))
    for e in edgs
        rem_edge!(g, e)
    end
end

copy(g::AGraphOrDiGraph) = deepcopy(g)

graphtype(::Type{G}) where {G<:AGraph} = G
digraphtype(::Type{G}) where {G<:ADiGraph} = G

edgetype(g::G) where {G<:AGraphOrDiGraph} = edgetype(G)
graphtype(g::G) where {G<:AGraphOrDiGraph} = graphtype(G)
digraphtype(g::G) where {G<:AGraphOrDiGraph} = digraphtype(G)
vertextype(g::G) where {G<:AGraphOrDiGraph} = vertextype(G)

graph(g::AGraph) = g
digraph(g::ADiGraph) = g

#### FALLBACKS #################

function (::Type{G})(g::AGraph) where G <: AGraph
    gnew = G(nv(g))
    for e in edges(g)
        add_edge!(gnew, src(e), dst(e))
    end
    return gnew
end

function (::Type{G})(g::ADiGraph) where G <: ADiGraph
    gnew = G(nv(g))
    for e in edges(g)
        add_edge!(gnew, src(e), dst(e))
    end
    return gnew
end

function (::Type{G})(adjmx::AbstractMatrix{T}
        ; selfedges=true) where {G<:ADiGraph, T<:Number}
    op =  selfedges ? ((x,y) -> true)  : (!=)
    return _graph_from_matr!(G, op, adjmx)
end

function (::Type{G})(adjmx::AbstractMatrix{T}
        ; upper=false, selfedges=true) where {G<:AGraph, T<:Number}

  op =  upper && !selfedges  ? (<) :
        upper && selfedges   ? (<=)  :
        !upper && !selfedges ? (!=) : ((x,y) -> true)

  return _graph_from_matr!(G, op, adjmx)
end

function _graph_from_matr!(G, op, adjmx::AbstractMatrix)
    n, m = size(adjmx)
    n == m || error("Adjacency / distance matrices must be square")
    g = G(n)
    for i in eachindex(adjmx)
        adjmx[i] == 0 && continue
        u, v = Tuple(CartesianIndices(adjmx)[i])
        op(u,v) || continue
        add_edge!(g, u, v)
    end
    return g
end

function _graph_from_matr!(G, op, adjmx::SparseMatrixCSC)
    n, m = size(adjmx)
    n == m || error("Adjacency / distance matrices must be square")
    g = G(n)
    rows = rowvals(adjmx)
    for v=1:n
        for j in nzrange(adjmx, v)
           u = rows[j]
           op(u,v) || continue
           add_edge!(g, u, v)
        end
    end
    return g
end

(::Type{G})(n::Integer, m::Integer; seed = -1) where {G<:AGraphOrDiGraph} = erdos_renyi(n, m, G; seed=seed)

function digraph(g::AGraph)
    G = digraphtype(g)
    h = G(nv(g))
    for e in edges(g)
        add_edge!(h, src(e), dst(e))
        add_edge!(h, dst(e), src(e))
    end
    return h
end

function graph(g::ADiGraph)
    G = graphtype(g)
    h = G(nv(g))
    for e in edges(g)
        add_edge!(h, src(e), dst(e))
    end
    return h
end

function ==(g::G, h::G) where G<:AGraphOrDiGraph
    nv(g) != nv(h) && return false
    ne(g) != ne(h) && return false
    for i=1:nv(g)
        if sort(collect(out_neighbors(g, i))) != sort(collect(out_neighbors(h, i)))
            return false
        end
    end
    return true
end

#not exported
"""
    in_adjlist(g)

Returns the backward adjacency list of a graph.
For each vertex the vector of neighbors though incoming edges.

    in_adjlist(g) == [collect(in_neighbors(i)) for i=1:nv(g)]

It is the same as [`adjlist`](@ref) and [`out_adjlist`](@ref) for
undirected graphs.

NOTE: returns a reference, not a copy. Do not modify result.
"""
in_adjlist(g::ADiGraph) = Vector{Int}[collect(in_neighbors(g, i)) for i=1:nv(g)]

#not exported
"""
    out_adjlist(g)

Returns the forward adjacency list of a graph, i.e. a vector of vectors
containing for each vertex the neighbors trhough outgoing edges.

    out_adjlist(g) == [collect(out_neighbors(i)) for i=1:nv(g)]

The adjacency list is be pre-calculated for most graph types.
It is the same as [`adjlist`](@ref) and [`in_adjlist`](@ref) for
undirected graphs and the same as [`adjlist`](@ref) for directed ones.

NOTE: It may return a reference, not a copy. Do not modify result.

"""
out_adjlist(g::AGraphOrDiGraph) = Vector{Int}[collect(out_neighbors(g, i)) for i=1:nv(g)]

"""
    has_edge(g, e)
    has_edge(g, u, v)

Returns true if the graph `g` has an edge `e` (from `u` to `v`).
"""
function has_edge(g::AGraph, u, v)
    (u > nv(g) || v > nv(g)) && return false
    if degree(g, u) > degree(g, v)
        u, v = v, u
    end
    return v ∈ neighbors(g, u)
end

function has_edge(g::ADiGraph, u, v)
    (u > nv(g) || v > nv(g)) && return false
    if out_degree(g, u) < in_degree(g, v)
        return v ∈ out_neighbors(g, u)
    else
        return u ∈ in_neighbors(g, v)
    end
end

"""
    in_edges(g, v)

Returns an iterator to the edges in `g` going to vertex `v`.
`v == dst(e)` for each returned edge `e`.
"""
in_edges(g::AGraphOrDiGraph, v) = (edge(g, x, v) for x in in_neighbors(g, v))

"""
    out_edges(g, v)

Returns an iterator to the edges in `g` coming from vertex `v`.
`v == src(e)` for each returned edge `e`.
"""
out_edges(g::AGraphOrDiGraph, v) = (edge(g, v, x) for x in out_neighbors(g, v))

"""
    edges(g, v)

Returns an iterator to the edges in `g` coming from vertex `v`.
`v == src(e)` for each returned edge `e`.

It is equivalent to [`out_edges`](@ref).

For digraphs, use [`all_edges`](@ref) to iterate over
both in and out edges.
"""
edges(g::AGraphOrDiGraph, v) = out_edges(g, v)

"""
    all_edges(g, v)

Iterates over all in and out edges of vertex `v` in `g`.
"""
all_edges(g::AGraph, v) = out_edges(g, v)
all_edges(g::ADiGraph, v) = flatten((out_edges(g, v), in_edges(g, v)))
#TODO fix chain eltype, since collect gives Any[...]

"""
    reverse(g::ADiGraph)

Produces a graph where all edges are reversed from the
original.
"""
function reverse(g::G) where G<:ADiGraph
    h = G(nv(g))
    for e in edges(g)
        add_edge!(h, reverse(e))
    end
    return h
end

"""
    reverse!(g::DiGraph)

In-place reverse (modifies the original graph).
"""
reverse!(g::ADiGraph) = nothing

add_edge!(g::AGraphOrDiGraph, e::AEdge) = add_edge!(g, src(e), dst(e))
rem_edge!(g::AGraphOrDiGraph, e::AEdge) = rem_edge!(g, src(e), dst(e))
has_edge(g::AGraphOrDiGraph, e::AEdge) = has_edge(g, src(e), dst(e))

"""
    unsafe_add_edge!(g, u, v)

Possibly faster and unsafer version of [`add_edge!`](@ref), which doesn't guarantee
some graph invariant properties.

For example, some graph types (e.g. `Graph`) assume sorted adjacency lists as members.
In this case order is not preserved while inserting new edges, resulting in a
faster construction of the graph. As a consequence though, some functions such
`has_edge(g, u, v)` could give incorrect results.

To restore the correct behaviour, call [`rebuild!`](@ref)(g) after the last
call to `unsafe_add_edge!`.
"""
unsafe_add_edge!(g::AGraphOrDiGraph, u, v) = add_edge!(g, u, v)
unsafe_add_edge!(g::AGraphOrDiGraph, e::AEdge) = unsafe_add_edge!(g, src(e), dst(e))

"""
    rebuild!(g)

Check and restore the structure of `g`, which could be corrupted by the
use of unsafe functions (e. g. [`unsafe_add_edge!`](@ref))
"""
rebuild!(g::AGraphOrDiGraph) = nothing

"""
    rem_vertex!(g, v)

Remove the vertex `v` from graph `g`.
It will change the label of the last vertex of the old graph to `v`.

See also [`rem_vertices!`](@ref).
"""
function rem_vertex!(g::AGraphOrDiGraph, v)
    v in vertices(g) || return false
    clean_vertex!(g, v)
    swap_vertices!(g, v, nv(g))
    pop_vertex!(g)
    return true
end

"""
    rem_vertices!(g, vs)
    rem_vertices!(g, v1, v2, ....)

Remove the vertices in `vs` from graph `g`.

Some vertices of `g` may be reindexed during the removal. To keep track of the reindexing,
a vertex map is returned, associating vertices with changed indexes to their old indexes.
"""
rem_vertices!(g, vs::Set) = _rem_vertices!(g, sort!(collect(vs)))
rem_vertices!(g, vs) = _rem_vertices!(g, sort!(union(vs)))
rem_vertices!(g, vs::Integer...) = _rem_vertices!(g, sort!(union(vs)))

# assume `vlist` sorted and not containing duplicates
function _rem_vertices!(g::AGraphOrDiGraph, vlist::Vector)
    n = nv(g)
    nrem = length(vlist)
    # TODO use default vertex map when available
    vmap = VertexMap(g, vertextype(g))
    vswap = n-nrem+1
    for v in vlist
        if v <= n-nrem
            while vswap in vlist
                vswap += 1
            end
            clean_vertex!(g, v)
            vmap[v] = vswap
            vswap += 1
        end
    end

    for v in vlist
        if v <= n-nrem
            swap_vertices!(g, v, vmap[v])
        end
    end
    for i=1:nrem
        pop_vertex!(g)
    end
    return vmap
end


text:
# surrogate_ground_truth_radon-mala-huberloss
#
# This script runs a long chain based on the
# Metropolis-adjusted Langevin algorithm (MALA). This is used as a gold
# standard as this procedure yields an asymptotically unbiased chain
# and thus is used to compare the efficiency of the biased chains produced
# in the experiment bias_variance-arwmh-huberloss.

using DataFrames
using Distributions: Normal

using SteinDistributions: SteinHuberLossRegressionGaussianPrior
using SteinSamplers: runmala

include("experiment_utils.jl")

# Set random number generator seed
@parseintcli seed "s" "seed" 7

# Set up the data
data_set = "radon"
rdat = readtable("src/experiments/radon/srrs2.dat")
cdat = readtable("src/experiments/radon/cty.dat")
rename!(cdat, :ctfips, :cntyfips)
# join the two datasets
dat = join(rdat, cdat, on = [:stfips, :cntyfips])
# subset only to Minnesota! (Gelman does this)
dat = dat[dat[:st] .== "MN",:]
# log radon
dat[:logradon] = log(dat[:activity])
dat[dat[:activity] .== 0, :logradon] = log(0.1)
# log uranium
dat[:loguranium] = log(dat[:Uppm])
dat[dat[:Uppm] .== 0, :loguranium] = log(0.1)
# make up X and y
X = Array(Float64, nrow(dat), 3)
X[:,1] = 1.0
X[:,2] = dat[:floor]
X[:,3] = dat[:loguranium]
y = convert(Array, dat[:logradon])

# N is number of samples, p is dimension of parameter
N, d = size(X)
# the max number of samples to run in the solver
n = 100000000
# name of prior + posterior
distname = "huberloss-gaussianprior"
# setup distribution
target = SteinHuberLossRegressionGaussianPrior(X, y)
# intiial guess for beta
beta0 = zeros(d)
# the epsilons to run the experiemnt
epsilons = [1e-4]

for epsilon = epsilons
    # Draw MALA samples
    (betas_mala, numgrad_mala) =
        @setseed runmala(target, beta0;
                         epsilonfunc=t -> epsilon,
                         numiter=n)
    # estimate the final beta moments
    betahat = vec(mean(betas_mala, 1))
    betavcov = cov(betas_mala)
    betasd = sqrt(diag(betavcov))
    # compute the running estimated yhat E{Y | X = x}
    # when using the first k values from the chain
    ks = round(Int, linspace(0, n, 1001))[2:end]
    meanyhats = Array(Float64, N, length(ks))
    for (i,k) in enumerate(ks)
        meanbeta = mean(betas_mala[1:k,:], 1)
        yhats = X * meanbeta'
        meanyhats[:,i] = vec(yhats)
    end
    # scale the mu's by the l_infty norm of each X_i
    scaledyhats = meanyhats ./ vec(maximum(X, 2))
    # fetch constants for computing confidence intervals of beta
    phi = Normal()
    zscore = quantile(phi, 1 - 0.05/2)
    # Package and save results
    summary_data = Dict{Any, Any}(
        "n" => n,
        "seed" => seed,
        "distname" => distname,
        "epsilon" => epsilon,
        "beta.hat" => betahat,
        "beta.vcov" => betavcov,
        "beta.sd" => betasd,
        "beta.hat.low" => betahat - zscore .* betasd,
        "beta.hat.high" => betahat + zscore .* betasd,
        "scaled.yhat" => scaledyhats,
        "scaled.yhat.ii" => ks,
        "numgrad" => numgrad_mala,
        "sampler" => "radon-mala"
    )
    save_json(
        summary_data;
        dir="bias_variance-arwmh-huberloss-MALA",
        distname=distname,
        dataset=data_set,
        n=n,
        epsilon=epsilon,
        d=d,
        seed=seed
    )
end


text:
macro horner(x, p...)
ex = esc(p[end])
for i = length(p)-1:-1:1
ex = :(muladd(t, $ex, $(esc(p[i]))))
end
Expr(:block, :(t = $(esc(x))), ex)
end

max_stars_repo_path: Module 2/Chapter04/macro_horner.jl

max_issues_repo_path: Module 2/Chapter04/macro_horner.jl

max_forks_repo_path: Module 2/Chapter04/macro_horner.jl


text:
platforms_and_sources = Dict(
    Platform("x86_64", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.2.4/11.4_20210831/cudnn-11.4-linux-x64-v8.2.4.15.tgz",
                      "0e5d2df890b9967efa6619da421310d97323565a79f05a1a8cb9b7165baad0d7")],
    Platform("powerpc64le", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.2.4/11.4_20210831/cudnn-11.4-linux-ppc64le-v8.2.4.15.tgz",
                      "af8749ca83fd6bba117c8bee31b787b7f204946e864294030ee0091eb7d3577e")],
    Platform("aarch64", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.2.4/11.4_20210831/cudnn-11.4-linux-aarch64sbsa-v8.2.4.15.tgz",
                      "48b11f19e9cd3414ec3c6c357ad228aebbd43282aae372d42cab2af67c32a08b")],
    Platform("x86_64", "windows") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.2.4/11.4_20210831/cudnn-11.4-windows-x64-v8.2.4.15.zip",
                      "f01594639de35c380b4e360673ccaf04cdb238578e4b284935ee3d5a45f51f3c")],
)

max_stars_repo_name: cmhamel/Yggdrasil

max_issues_repo_name: cmhamel/Yggdrasil

max_forks_repo_name: cmhamel/Yggdrasil


text:
@coeff function a(reactor::Reactor)
  epsilon * R_0
end


text:
__precompile__(false)

module Libfreenect2

export Freenect2,enumerateDevices, getDeviceSerialNumber,
    getDefaultDeviceSerialNumber,openDevice, openDefaultDevice,
    CpuPacketPipeline, OpenGLPacketPipeline, OpenCLPacketPipeline,
    OpenCLKdePacketPipeline, getSerialNumber, getFirmwareVersion,
    getColorCameraParams, getIrCameraParams, setColorFrameListener,
    setIrAndDepthFrameListener, start, startStreams, stop, close,
    SyncMultiFrameListenerPtr, hasNewFrame, waitForNewFrame, release,
    FrameMap, FrameType, FramePtr, Registration, apply, undistortDepth

# Load dependency
deps = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(deps)
    include(deps)
else
    error("Libfreenect2 not properly installed. Please run Pkg.build(\"Libfreenect2\")")
end

import Base: convert, getindex, start, close, gamma

using Cxx

Libdl.dlopen_e(libfreenect2, Libdl.RTLD_GLOBAL)

function include_headers(top)
    addHeaderDir(top, kind=C_System)
    addHeaderDir(joinpath(top, "libfreenect2"), kind=C_System)
    for name in [
        "frame_listener.hpp",
        "frame_listener_impl.h",
        "libfreenect2.hpp",
        "packet_pipeline.h",
        "registration.h",
        ]
        cxxinclude(joinpath(top, "libfreenect2/$name"))
    end
end

const system_include_top = "/usr/local/include"
const local_include_top = joinpath(Pkg.dir("Libfreenect2", "deps", "usr", "include"))

if isdir(local_include_top)
    # info("Including headers from local path: $local_include_top")
    include_headers(local_include_top)
elseif isdir(joinpath(system_include_top, "libfreenect2"))
    # info("Including headers from system path: $system_include_top")
    include_headers(system_include_top)
else
    error("Cannot find libfreenect2 headers")
end

CpuPacketPipeline() = icxx"new libfreenect2::CpuPacketPipeline();"
OpenGLPacketPipeline() = icxx"new libfreenect2::OpenGLPacketPipeline();"
OpenCLPacketPipeline() = icxx"new libfreenect2::OpenCLPacketPipeline();"
OpenCLKdePacketPipeline() = icxx"new libfreenect2::OpenCLKdePacketPipeline();"

import Cxx: CppEnum
const Libfreenect2FrameType = CppEnum{Symbol("libfreenect2::Frame::Type"),Int32}

# This is used in the SyncMultiFrameListener constructor, since it requires frame types
# to be specified by unsigned int
const FRAME_COLOR = Cuint(1)
const FRAME_IR    = Cuint(2)
const FRAME_DEPTH = Cuint(4)

module FrameType

import ..Libfreenect2: Libfreenect2FrameType

const COLOR = Libfreenect2FrameType(1)
const IR = Libfreenect2FrameType(2)
const DEPTH = Libfreenect2FrameType(4)

end # module FrameType
end

### FrameMap ###

type FrameMap
    handle::cxxt"libfreenect2::FrameMap"
end

function release(jlframes::FrameMap)
    frames = jlframes.handle
    icxx"""
    for (libfreenect2::FrameMap::iterator it = $frames.begin();
         it != $frames.end(); ++it) {
      if (it->second != nullptr) {
        delete it->second;
        it->second = nullptr;
      }
    }
    $frames.clear();
    """
end

### Frame ###

type FramePtr
    handle::pcpp"libfreenect2::Frame"
    frame_type::Cuint

    # Allocated Frame should be set from
    # SyncMultiFrameListener::waitForNewFrame
    FramePtr(frame::pcpp"libfreenect2::Frame", key=0) = new(frame, key)

    function FramePtr(width, height, bytes_per_pixel; key=0)
        frame = icxx"new libfreenect2::Frame($width, $height, $bytes_per_pixel);"
        p = new(frame, key)
        finalizer(p, destroy)
        p
    end
end

function destroy(frame::FramePtr)
    icxx"""
    if ($(frame.handle) != nullptr) {
        delete $(frame.handle);
        $(frame.handle) = nullptr;
    }
    """
end

using CxxStd
function getindex(map::CxxStd.StdMap, key::Libfreenect2FrameType)
    icxx"""
        libfreenect2::Frame* f = $map[$key];
        return f;
    """
end

function getindex(frames::FrameMap, key::Libfreenect2FrameType)
    @assert isa(key, Libfreenect2FrameType)
    FramePtr(getindex(frames.handle, key), key.val)
end

for member in [
    :timestamp,
    :sequence,
    :data,
    :exposure,
    :gain,
    :gamma
    ]
    body = "\$(frame.handle)->$(member);"
    ex = Expr(:macrocall, Symbol("@icxx_str"), body)
    @eval $member(frame::FramePtr) = $ex
end

for member in [
    :width,
    :height,
    :bytes_per_pixel,
    ]
    body = "\$(frame.handle)->$(member);"
    ex = Expr(:macrocall, Symbol("@icxx_str"), body)
    @eval $member(frame::FramePtr) = convert(Int, $ex)
end

### Frame to Array conversion ###

"""Convenient function to convert Frame* to Array
"""
@inline function _asarray(frame::FramePtr; do_reshape::Bool=true)
    data_ptr = data(frame)::Ptr{UInt8}

    if frame.frame_type == FRAME_COLOR
        total_length = width(frame) * height(frame) * 4
        array = unsafe_wrap(Array, data_ptr, total_length)
        if do_reshape
            array = reshape(array, 4, width(frame), height(frame))
        end
        return array
    elseif frame.frame_type == FRAME_IR || frame.frame_type == FRAME_DEPTH
        total_length = width(frame) * height(frame) * 1
        array = unsafe_wrap(Array, convert(Ptr{Float32}, data_ptr), total_length)
        if do_reshape
            array = reshape(array, width(frame), height(frame))
        end
        return array
    else
        error("Cannnot determine type of raw data")
    end
end

function convert(::Type{Array}, frame::FramePtr)
    _asarray(frame, do_reshape=true)
end
function convert(::Type{Vector}, frame::FramePtr)
    _asarray(frame, do_reshape=false)
end

# TODO: should use multiple dispatch instead?
function convert{T,N}(::Type{Array{T,N}}, frame::FramePtr)
    typ = frame.frame_type
    if typ == FRAME_COLOR && T == UInt8 && N == 3
        return _asarray(frame, do_reshape=true)
    elseif typ == FRAME_COLOR && T == UInt8 && N == 1
        return _asarray(frame, do_reshape=false)
    elseif (typ == FRAME_IR || typ == FRAME_DEPTH) && T == Float32 && N == 2
        return _asarray(frame, do_reshape=true)
    elseif (typ == FRAME_IR || typ == FRAME_DEPTH) && T == Float32 && N == 1
        return _asarray(frame, do_reshape=false)
    else
        error("unsupported conversion")
    end
end

### SyncMultiFrameListener ###

type SyncMultiFrameListenerPtr
    handle::cxxt"std::shared_ptr<libfreenect2::SyncMultiFrameListener>"

    function SyncMultiFrameListenerPtr(
            frame_types=FRAME_COLOR | FRAME_IR | FRAME_DEPTH)
        new(icxx"""
            std::shared_ptr<libfreenect2::SyncMultiFrameListener>(
                new libfreenect2::SyncMultiFrameListener($frame_types));
        """)
    end
end

# TODO: remove this inefficient glue code
cxx"""
namespace legacy {
libfreenect2::FrameMap getFrameMap(
    libfreenect2::SyncMultiFrameListener* listener) {
    libfreenect2::FrameMap frames;
    listener->waitForNewFrame(frames);
    return frames;
}
}
"""
function waitForNewFrame(listener::SyncMultiFrameListenerPtr)
    frames = icxx"""return legacy::getFrameMap($(listener.handle).get());"""
    return FrameMap(frames)
end

function release(listener::SyncMultiFrameListenerPtr, jlframes::FrameMap)
    frames = jlframes.handle
    # really ugly
    icxx"""
    libfreenect2::FrameMap cxxframes;
    for (libfreenect2::FrameMap::iterator it = $frames.begin();
         it != $frames.end(); ++it) {
        cxxframes[it->first] = it->second;
    }
    $(listener.handle)->release(cxxframes);
    """
end

hasNewFrame(listener::SyncMultiFrameListenerPtr) =
    icxx"$(listener.handle)->hasNewFrame()"

"""libfreenect2::Freenect2"""
typealias Freenect2 cxxt"libfreenect2::Freenect2"
(::Type{Freenect2})() = icxx"libfreenect2::Freenect2();"

for name in [
    :enumerateDevices,
    ]
    ex = Expr(:macrocall, Symbol("@icxx_str"), "\$f.$(name)();")
    @eval $name(f::Freenect2) = $ex
end

function openDefaultDevice(f::Freenect2, pipeline=Union{})
    if is(pipeline, Union{})
        icxx"$f.openDefaultDevice();"
    else
        icxx"$f.openDefaultDevice($pipeline);"
    end
end

function getDeviceSerialNumber(f::Freenect2, idx)
    unsafe_string(icxx"$f.getDeviceSerialNumber($idx);")
end

function getDefaultDeviceSerialNumber(f::Freenect2)
    unsafe_string(icxx"$f.getDefaultDeviceSerialNumber();")
end

function openDevice(f::Freenect2, name, pipeline=Union{})
    arg1 = isa(name, AbstractString) ? pointer(name) : name

    if is(pipeline, Union{})
        icxx"$f.openDevice($arg1);"
    else
        icxx"$f.openDevice($arg1, $pipeline);"
    end
end

"""libfreenect2::Freenect2Device*"""
const Freenect2DevicePtr = pcpp"libfreenect2::Freenect2Device"

function getSerialNumber(device::Freenect2DevicePtr)
    unsafe_string(icxx"$device->getSerialNumber();")
end

function getFirmwareVersion(device::Freenect2DevicePtr)
    unsafe_string(icxx"$device->getFirmwareVersion();")
end

import Base: start, close

for f in [
    :getColorCameraParams,
    :getIrCameraParams,
    ]
    ex = Expr(:macrocall, Symbol("@icxx_str"), "\$device->$f();")
    @eval $f(device::Freenect2DevicePtr) = $ex
end

for f in [
    :start,
    :stop,
    :close,
    ]
    cxxcall = Expr(:macrocall, Symbol("@icxx_str"), "\$device->$f();")
    @eval begin
        function $f(device::Freenect2DevicePtr)
            r = $cxxcall
            r || error("problem happens in device operation")
        end
    end
end

function startStreams(device::Freenect2DevicePtr, rgb::Bool, depth::Bool)
    icxx"$(device)->startStreams($rgb, $depth);"
end

function setColorFrameListener(device::Freenect2DevicePtr,
        listener::SyncMultiFrameListenerPtr)
    icxx"$device->setColorFrameListener($(listener.handle).get());"
end

function setIrAndDepthFrameListener(device::Freenect2DevicePtr,
        listener::SyncMultiFrameListenerPtr)
    icxx"$device->setIrAndDepthFrameListener($(listener.handle).get());"
end

const ColorCameraParams = cxxt"libfreenect2::Freenect2Device::ColorCameraParams"
const IrCameraParams = cxxt"libfreenect2::Freenect2Device::IrCameraParams"

type Registration
    handle::cxxt"std::shared_ptr<libfreenect2::Registration>"

    function Registration(irparams::IrCameraParams, cparams::ColorCameraParams)
        p = icxx"std::shared_ptr<libfreenect2::Registration>(
            new libfreenect2::Registration($irparams, $cparams));"
    new(p)
    end
end

function apply{T<:FramePtr}(registration::Registration, color::T, depth::T,
    undistorted::T, registered::T; enable_filter::Bool=true)
    icxx"$(registration.handle)->apply($(color.handle), $(depth.handle),
        $(undistorted.handle), $(registered.handle), $enable_filter);"
end

function undistortDepth{T<:FramePtr}(registration::Registration, depth::T,
    undistorted::T)
    icxx"""
        $(registration.handle)->undistortDepth(
            $(depth.handle), $(undistorted.handle));
    """
end

# deprecates
@deprecate SyncMultiFrameListener SyncMultiFrameListenerPtr
@deprecate FrameMapContainer FrameMap
@deprecate FrameContainer FramePtr

end # module Libfreenect2


text:
const PROJ_IDER = "iJR904"
const BIOMASS_IDER = "BiomassEcoli"
const ATPM_IDER = "ATPM"
const COST_IDER = "tot_cost"
const EX_GLC_IDER = "EX_glc_LPAREN_e_RPAREN_"
const MAX_ABS_BOUND = 100


text:
#!/usr/bin/env julia

include("../src/read.jl")

vcf  = ARGS[1] # Variant Calling file e.g: chr3.vcf.gz

"""
Computing consensus fo variants with at least 3 reads
"""

function consensus(vcf_line)
	cols = split(vcf_line)
	alts = []
	for sample in cols[5:end]
		altRead = parse(Int64, split(sample, ':')[1])
		if altRead >= 5
			push!(alts, 1)
		else
			push!(alts, 0)
		end
	end
	consensus = (sum(alts)/length(alts))*100
end

"""
Filtering each position from input vcf with the consensus
value of 10% and returning filtered variants in .vcf format
"""
function filterVCF(vcf_path)
	lines = read(vcf_path)
	println(chomp(lines[1]))
	for line = lines[2:end]
		cons = consensus(line)
		# minimum 3 reads having mutation in at least 10% of thae samples
		if cons >= 10
			print(line)
		end
	end
end

filterVCF(vcf)


text:
using ModelingToolkit, OrdinaryDiffEq
using Plots

@parameters t
D = Differential(t)

@connector function Pin(;name)
    sts = @variables v(t)=1.0 i(t)=1.0
    ODESystem(Equation[], t, sts, []; name=name)
end

# tell MTK how to connect equations
function ModelingToolkit.connect(::Type{Pin}, ps...)
    # junction laws
    eqs = [0 ~ sum(p->p.i, ps)]

    for i ∈1:length(ps)-1
        push!(eqs, ps[i].v ~ ps[i+1].v)
    end
    eqs
end

# define a ground component
function Ground(;name)
    @named g = Pin()
    eqs = [0 ~ g.v]
    compose(ODESystem(eqs, t, [], []; name=name), g)
end

# abstract component with b.c.
function OnePort(;name)
    @named p = Pin()
    @named n = Pin()
    sts = @variables v(t)=1.0 i(t)=1.0
    eqs = [v ~ p.v - n.v,
           0 ~ p.i + n.i,
           i ~ p.i]
    compose(ODESystem(eqs, t, sts, []; name=name), p, n)
end

# extend OnePort for Capacitor, Resistor, Voltage Source components
function Resistor(;name, R=1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters R=R
    eqs = [v ~ i * R] # Ohm's law
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Capacitor(;name, C=1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps =@parameters C=C
    eqs = [D(v) ~  i / C]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function ConstantVoltage(;name, V=1.0)  # i.e. a battery
    @named oneport = OnePort()
    @unpack v = oneport
    ps = @parameters V=V
    eqs = [v ~ V]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

R = 1.0
C = 1.0
V = 1.0

@named resistor = Resistor(R=R)
@named capacitor = Capacitor(C=C)
@named source = ConstantVoltage(V=V)
@named ground =Ground()

# define the system by connecting the components!
rc_eqs = [
    connect(source.p, resistor.p)
    connect(resistor.n, capacitor.p)
    connect(capacitor.n, ground.g, source.n)
]
@named rc_model = compose(ODESystem(rc_eqs, t), [resistor, capacitor, source, ground])

# We have 20 equations but we should only have 1!!!
sys = structural_simplify(rc_model)

u0 = [capacitor.v => 0.0]
prob = ODAEProblem(sys, u0, (0.0, 10.0))

sol = solve(prob, Tsit5())

plot(sol, vars=[capacitor.v, resistor.v])


text:
using KUparser, KUnet, HDF5, JLD
@date @load "zn11oparse.jld"
@show KUnet.atype(Array)
@show niter = 0
@show net = Layer[Perceptron(size(ydev,1))]
@date train(net, xdev, ydev; iters=niter)
@show net = Layer[Perceptron(size(ydev,1))]
@date train(net, xdev, ydev; iters=niter)
@show net = Layer[Perceptron(size(ydev,1))]
@date @profile train(net, xdev, ydev; iters=niter)
Profile.print()

# for i=1:100
#     @date train(net, xtst, ytst)
#     @date println((i, accuracy(ytst, predict(net, xtst)),
#                    accuracy(ydev, predict(net, xdev))))
# end
# :ok


text:
using Statistics
using Plots
using Flux
using OceanParameterizations
using Oceananigans.Grids
using BSON
using OrdinaryDiffEq, DiffEqSensitivity
using WindMixing
using JLD2
using FileIO

PATH = joinpath(pwd(), "extracted_training_output")
# PATH = "D:\\University Matters\\Massachusetts Institute of Technology\\CLiMA Project\\OceanParameterizations.jl\\training_output"

DATA_NAME = "NDE_training_mpp_8sim_wind_mixing_cooling_diffusivity_1e-1_Ri_1e-1_weights_divide1f5_gradient_smallNN_scale_5e-3_rate_1e-4"
DATA_PATH = joinpath(PATH, "$(DATA_NAME)_extracted.jld2")
ispath(DATA_PATH)
                    # FILE_PATH = "D:\\University Matters\\Massachusetts Institute of Technology\\CLiMA Project\\OceanParameterizations.jl\\training_output"
FILE_PATH = joinpath(pwd(), "Output")
# VIDEO_NAME = "u_v_T_mpp_2sim_wind_mixing_-1e-3_cooling_5e-8_diffusivity_1e-1_Ri_1e-1_weights_divide1f5_gradient_smallNN_scale_1e-2_rate_1e-4_test_5e-8"
VIDEO_NAME = "test"
# SIMULATION_NAME = "NN Smoothing Wind-Mixing, Testing Data"
SIMULATION_NAME = "4 Simulation Training (Wind Mixing, Free Convection)"

# file = jldopen(DATA_PATH, "r")
file = jldopen(DATA_PATH, "r")

losses = file["losses"]

minimum(losses)

train_files = file["training_info/train_files"]
train_parameters = file["training_info/parameters"]

if haskey(file["training_info"], "loss_scalings")
    loss_scalings = file["training_info/loss_scalings"]
elseif haskey(train_parameters, "gradient_scaling")
    gradient_scaling = train_parameters["gradient_scaling"]
    loss_scalings = (u=1f0, v=1f0, T=1f0, ∂u∂z=gradient_scaling, ∂v∂z=gradient_scaling, ∂T∂z=gradient_scaling)
end

# Plots.plot(1:1:length(losses), losses, yscale=:log10)
# Plots.xlabel!("Iteration")
# Plots.ylabel!("Loss mse")
# savefig(joinpath(PATH, "Output", "NDE_training_modified_pacanowski_philander_1sim_-1e-3_smaller_learning_rate_loss.pdf"))
# train_files = ["-1e-3"]
𝒟train = WindMixing.data(train_files, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=false)

test_files = ["wind_-1e-3_heating_-4e-8"]
𝒟test = WindMixing.data(test_files, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=false)
uw_NN = file["neural_network/uw"]
vw_NN = file["neural_network/vw"]
wT_NN = file["neural_network/wT"]

close(file)
[uw_NN(rand(96)) uw_NN(rand(96))]

# N_inputs = 96
# hidden_units = 400
# N_outputs = 31

# weights, re = Flux.destructure(Chain(Dense(N_inputs, hidden_units, relu), Dense(hidden_units, N_outputs)))

# uw_NN = re(zeros(Float32, length(weights)))
# vw_NN = re(zeros(Float32, length(weights)))
# wT_NN = re(zeros(Float32, length(weights)))

𝒟test = WindMixing.data(test_files, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=false)
trange = 1:1:1153
plot_data = NDE_profile(uw_NN, vw_NN, wT_NN, test_files[1], 𝒟test, 𝒟train, trange,
                        modified_pacanowski_philander=train_parameters["modified_pacanowski_philander"],
                        # ν₀=1f-4, ν₋=0.1f0, ΔRi=1f-1,
                        ν₀=train_parameters["ν₀"], ν₋=train_parameters["ν₋"], ΔRi=train_parameters["ΔRi"],
                        Riᶜ=train_parameters["Riᶜ"], convective_adjustment=train_parameters["convective_adjustment"],
                        smooth_NN=train_parameters["smooth_NN"], smooth_Ri=train_parameters["smooth_Ri"],
                        zero_weights=train_parameters["zero_weights"],
                        loss_scalings=loss_scalings)

# WindMixing.animate_profiles_fluxes(plot_data, joinpath(FILE_PATH, VIDEO_NAME), dimensionless=false, SIMULATION_NAME=SIMULATION_NAME)

animation_type = "Training"
n_trainings = length(train_files)
training_types = "Wind Mixing, Free Convection"
VIDEO_NAME = "test"
animate_profiles_fluxes_comparison(plot_data, plot_data, plot_data, joinpath(FILE_PATH, VIDEO_NAME), fps=30,
                                                animation_type=animation_type, n_trainings=n_trainings, training_types=training_types)

# VIDEO_NAME = "u_v_T_modified_pacanowski_philander_1sim_-1e-3_test2"

# keys(plot_data)

# plot_data["truth_T"][:,1]

# uvT_truth = [plot_data["truth_u"]; plot_data["truth_v"]; plot_data["truth_T"]]
# Ris = local_richardson(uvT_truth, 𝒟test, unscale=true)

# animate_local_richardson_profile(uvT_truth, 𝒟test, joinpath(FILE_PATH, "Ris_convective_adjustment_1sim_-1e-3_2_test"), unscale=true)

# plot(Ris[:,3], plot_data["depth_flux"])
# xlabel!("Ri")
# ylabel!("z")

# animate_profile_flux(plot_data, "u", "uw", joinpath(FILE_PATH, "u_uw_modified_pacanowski_philander_1sim_-1e-3_test"), gif=true, dimensionless=false)
# animate_profile_flux(plot_data, "v", "vw", joinpath(FILE_PATH, "v_vw_modified_pacanowski_philander_1sim_-1e-3_test"), gif=true, dimensionless=false)
# animate_profile_flux(plot_data, "T", "wT", joinpath(FILE_PATH, "w_wT_modified_pacanowski_philander_1sim_-1e-3_test"), gif=true, dimensionless=false)

# animate_profiles(plot_data, joinpath(FILE_PATH, VIDEO_NAME), dimensionless=false)

# animate_profile(plot_data, "u", "uw", joinpath(FILE_PATH, "u_uw_convective_adjustment_viscosity_empty"), gif=true)
# animate_profile(plot_data, "v", "vw", joinpath(FILE_PATH, "v_vw_convective_adjustment_viscosity_empty"), gif=true)
# animate_profile(plot_data, "T", "wT", joinpath(FILE_PATH, "w_wT_convective_adjustment_viscosity_empty"), gif=true)

# animate_flux(plot_data, "uw", joinpath(FILE_PATH, "uw_test"))
# animate_flux(plot_data, "vw", joinpath(FILE_PATH, "vw_test"))
# animate_flux(plot_data, "wT", joinpath(FILE_PATH, "wT_test"))

# train_files = ["-1e-3", "-9e-4", "-8e-4", "-7e-4", "-6e-4", "-5e-4", "-4e-4", "-3e-4", "-2e-4"]
train_files = ["wind_-1e-3_heating_-4e-8"]
# train_files = ["cooling_1e-8"]

𝒟train = WindMixing.data(train_files, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=false)
𝒟train.T.coarse
𝒟train.t
VIDEO_NAME = "u_v_T_3sim_-1e-3_-8e-4_-5e-4_diffusivity_1e-1_Ri_1e-1_weights_divide1f5_gradient_smallNN_scale_1e-2_rate_2e-4_test_-1e-3"
VIDEO_NAME = "test_video"

animate_training_data_profiles_fluxes(train_files, joinpath(FILE_PATH, VIDEO_NAME))

test_files = ["-1e-3", "cooling_5e-8", "-8e-4", "cooling_4e-8"]

animate_training_results(test_files, DATA_NAME, trange=1:1:1153)


text:
struct RibbonGeometry <: Gadfly.GeometryElement
    default_statistic::Gadfly.StatisticElement
    fill::Bool
    tag::Symbol
end
RibbonGeometry(default_statistic=Gadfly.Stat.identity(); fill=true, tag=empty_tag) =
        RibbonGeometry(default_statistic, fill, tag)

"""
    Geom.ribbon[(; fill=true)]

Draw a ribbon at the positions in `x` bounded above and below by `ymax` and
`ymin`, respectively.  Optionally draw multiple ribbons by grouping with `color` and `alpha` (for `fill=true`),
 or `color` and `linestyle` (for `fill=false`).
"""
const ribbon = RibbonGeometry

default_statistic(geom::RibbonGeometry) = geom.default_statistic

element_aesthetics(::RibbonGeometry) = [:x, :ymin, :ymax, :y, :xmin, :xmax, :color, :linestyle, :alpha]

function render(geom::RibbonGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    orientation = if isempty(Gadfly.undefined_aesthetics(aes, :x, :ymin, :ymax))
        :vertical
    elseif isempty(Gadfly.undefined_aesthetics(aes, :y, :xmin, :xmax))
        :horizontal
    else
        error("For `Geom.ribbon`, define either `x`, `ymin`, `ymax` or `y`, `xmin`, `xmax`")
    end

    var, minvar, maxvar = if orientation==:vertical
        :x, :ymin, :ymax
    else
        :y, :xmin, :xmax
    end

    Gadfly.assert_aesthetics_defined("Geom.ribbon", aes, var, minvar, maxvar)
    Gadfly.assert_aesthetics_equal_length("Geom.ribbon", aes, var, minvar, maxvar)

    vals = getfield(aes, var)
    maxvals = getfield(aes, maxvar)
    minvals = getfield(aes, minvar)

    default_aes = Gadfly.Aesthetics()
    default_aes.linestyle = [1]
    default_aes.color = Colorant[theme.default_color]
    default_aes.alpha = [1]
    aes = inherit(aes, default_aes)

    aes_vals, aes_minvals, aes_maxvals, aes_color, aes_linestyle, aes_alpha =
         concretize(vals, minvals, maxvals, aes.color, aes.linestyle, aes.alpha)
    XT, CT, LST, AT = eltype(aes_vals), eltype(aes_color), eltype(aes_linestyle), eltype(aes_alpha)
    YT = eltype(aes_minvals)
    groups = collect((Tuple{CT, LST, AT}), Compose.cyclezip(aes_color, aes_linestyle, aes_alpha))
    ugroups = unique(groups)
    nugroups = length(ugroups)

    V = Vector{Tuple{XT, YT}}
    K = Tuple{CT, LST, AT}

    max_points = Dict{K, V}(g=>V[] for g in ugroups)
    min_points = Dict{K, V}(g=>V[] for g in ugroups)
    orderf = first
    if orientation==:vertical
        for (x, ymin, ymax, c, ls, a) in Compose.cyclezip(aes_vals, aes_minvals, aes_maxvals, aes_color, aes_linestyle, aes_alpha)
            push!(max_points[(c,ls,a)], (x, ymax))
            push!(min_points[(c,ls,a)], (x, ymin))
        end
    else
        for (x, ymin, ymax, c, ls, a) in Compose.cyclezip(aes_vals, aes_minvals, aes_maxvals, aes_color, aes_linestyle, aes_alpha)
            push!(max_points[(c,ls,a)], (ymax, x))
            push!(min_points[(c,ls,a)], (ymin, x))
        end
        orderf = last
    end
    kys = keys(max_points)
    for k in kys
        sort!(max_points[k], by=orderf)
        sort!(min_points[k], by=orderf, rev=true)
    end

    polys = [collect(Tuple{XT, YT}, Iterators.flatten((min_points[k], max_points[k]))) for k in kys]
    lines = [collect(Tuple{XT, Union{YT, Missing}}, Iterators.flatten((min_points[k], [(last(min_points[k])[1], missing)], max_points[k]))) for k in kys]

    classes = Vector{String}(undef, nugroups)
    colors = Vector{Colorant}(undef, nugroups)
    linestyles = Vector{Vector{Measure}}(undef, nugroups)
    linestyle_palette_length = length(theme.line_style)
    alphas = Vector{Float64}(undef, nugroups)
    alpha_discrete  = AT <: Int
    linestyle_discrete = LST <: Int

    for (i, (c,ls,a)) in enumerate(kys)
        classes[i] = svg_color_class_from_label(aes.color_label([c])[1])
        colors[i] = parse_colorant(theme.lowlight_color(c))
        linestyles[i] = linestyle_discrete ? get_stroke_vector(theme.line_style[mod1(ls, linestyle_palette_length)]) : get_stroke_vector(ls)
        alphas[i] = alpha_discrete ? theme.alphas[a] : a
    end

    ctx = context()
    geom.fill ?
        compose!(ctx, (context(), Compose.polygon(polys, geom.tag),
                       fill(colors), fillopacity(alphas), svgclass(classes))) :
        compose!(ctx, (context(), Compose.line(lines, geom.tag),
                       stroke(colors), strokedash(linestyles), svgclass(classes)))

    return compose!(ctx, svgclass("geometry"), linewidth(theme.line_width))
end


text:
using FastGaussQuadrature
using Test, LinearAlgebra, Random, SpecialFunctions

@testset "FastGaussQuadrature.jl" begin
    include("test_gausschebyshev.jl")
    include("test_gausslegendre.jl")
    include("test_gaussjacobi.jl")
    include("test_gaussradau.jl")
    include("test_gausslobatto.jl")
    include("test_besselroots.jl")
    include("test_gausslaguerre.jl")
    include("test_gausshermite.jl")
end


text:
#AM to MPCC
export AMtoMPCC

using AmplNLReader

function AMtoMPCC(AM::AmplModel)

 nc=find(x->x == 0,AM.cvar)
 ncc=find(x->x > 0,AM.cvar)
 n_cc=length(ncc)
 ncon=length(nc)
 nvar=AM.meta.nvar

 #pas optimal comme technique :
 Temp=sparse(zeros(n_cc,nvar))
 for i=1:n_cc
  Temp[i,AM.cvar[ncc[i]]]=1.0
 end
 maj(t)=Temp'*t
 Temp1=sparse(zeros(n_cc,ncon+n_cc))
 for i=1:n_cc
  Temp1[i,ncc[i]]=1.0
 end
 maj1(t)=Temp1'*t
 Temp2=sparse(zeros(ncon,ncon+n_cc))
 for i=1:ncon
  Temp2[i,nc[i]]=1.0
 end
 maj2(t)=Temp2'*t

#doublon de la contrainte de positivité...
nlp=SimpleNLPModel(x->obj(AM,x),AMnlp.meta.x0,
                   lvar=AM.meta.lvar,uvar=AM.meta.uvar,
                   c=x->cons(AM,x)[nc],lcon=AM.meta.lcon[nc],
                   ucon=AM.meta.ucon[nc],g=x->grad(AM,x),
                   H=(x;y=zeros(ncon),obj_weight=1.0)->hess(AM,x;y=maj2(y),obj_weight=obj_weight),
                   J=x->jac(AM,x)[nc,1:nvar],
                   Jtp=(x,v)->jac(AM,x)[nc,1:nvar]'*v)

G=SimpleNLPModel(()->(),AM.meta.x0,c=x->cons(AM,x)[ncc],
                   lcon=AM.meta.lcon[ncc],
                   ucon=AM.meta.ucon[ncc],J=x->jac(AM,x)[ncc,1:nvar],
                   Jtp=(x,v)->(jac(AM,x)[ncc,1:nvar])'*v,
                   H=(x;y=zeros(n_cc),obj_weight=1.0)->hess(AM,x;y=maj1(y),obj_weight=0.0))

H=SimpleNLPModel(()->(),AM.meta.x0,c=x->x[AM.cvar[ncc]],
                   lcon=AM.meta.lvar[AM.cvar[ncc]],
                   ucon=AM.meta.uvar[AM.cvar[ncc]],
                   J=x->Temp,Jtp=(x,v)->maj(v),
                   H=(x;y=zeros(n_cc),obj_weight=1.0)->zeros(nvar,nvar))

 return nlp,G,H
end


text:
# ---
# title: 871. Minimum Number of Refueling Stops
# id: problem871
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming, Heap
# link: <https://leetcode.com/problems/minimum-number-of-refueling-stops/description/>
# hidden: true
# ---
#
# A car travels from a starting position to a destination which is `target`
# miles east of the starting position.
#
# Along the way, there are gas stations.  Each `station[i]` represents a gas
# station that is `station[i][0]` miles east of the starting position, and has
# `station[i][1]` liters of gas.
#
# The car starts with an infinite tank of gas, which initially has `startFuel`
# liters of fuel in it.  It uses 1 liter of gas per 1 mile that it drives.
#
# When the car reaches a gas station, it may stop and refuel, transferring all
# the gas from the station into the car.
#
# What is the least number of refueling stops the car must make in order to
# reach its destination?  If it cannot reach the destination, return `-1`.
#
# Note that if the car reaches a gas station with 0 fuel left, the car can still
# refuel there.  If the car reaches the destination with 0 fuel left, it is
# still considered to have arrived.
#
#
#
# **Example 1:**
#
#
#
#     Input: target = 1, startFuel = 1, stations = []
#     Output: 0
#     Explanation: We can reach the target without refueling.
#
#
# **Example 2:**
#
#
#
#     Input: target = 100, startFuel = 1, stations = [[10,100]]
#     Output: -1
#     Explanation: We can't reach the target (or even the first gas station).
#
#
# **Example 3:**
#
#
#
#     Input: target = 100, startFuel = 10, stations = [[10,60],[20,30],[30,30],[60,40]]
#     Output: 2
#     Explanation:
#     We start with 10 liters of fuel.
#     We drive to position 10, expending 10 liters of fuel.  We refuel from 0 liters to 60 liters of gas.
#     Then, we drive from position 10 to position 60 (expending 50 liters of fuel),
#     and refuel from 10 liters to 50 liters of gas.  We then drive to and reach the target.
#     We made 2 refueling stops along the way, so we return 2.
#
#
#
#
# **Note:**
#
#   1. `1 <= target, startFuel, stations[i][1] <= 10^9`
#   2. `0 <= stations.length <= 500`
#   3. `0 < stations[0][0] < stations[1][0] < ... < stations[stations.length-1][0] < target`
#
#
## @lc code=start
using LeetCode

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


text:
function bandwidth(A::SparseMatrixCSC{Tv, Ti}) where {Tv, Ti}
    m, n = size(A)
    bot = zeros(Ti, n + 1)
    bot[n + 1] = m + 1
    for j = n:-1:1
        if A.colptr[j] < A.colptr[j + 1]
            bot[j] = min(bot[j + 1], A.rowval[A.colptr[j]])
        else
            bot[j] = bot[j + 1]
        end
    end
    top = zeros(Ti, n + 1)
    top[1] = 0
    for j = 1:n
        if A.colptr[j] < A.colptr[j + 1]
            top[j + 1] = max(top[j], A.rowval[A.colptr[j + 1] - 1])
        else
            top[j + 1] = top[j]
        end
    end

    return maximum(top[2:n + 1] .- bot[1:n])
end

@testset "CuthillMcKee" begin
    A = matrices["HB/can_292"]
    (σ, τ) = permute_plaid(A, CuthillMcKeePermuter())
    @test σ.prm === τ.prm
    (σ′, τ′) = permute_plaid(A, CuthillMcKeePermuter(reverse=true))
    @test σ.prm == reverse(σ′.prm)
    @test τ.prm == reverse(τ′.prm)
    A = SparseMatrixCSC(A)
    (σ, τ) = permute_plaid(A, CuthillMcKeePermuter())
    @test σ.prm === τ.prm
    (σ′, τ′) = permute_plaid(A, CuthillMcKeePermuter(reverse=true))
    @test σ.prm == reverse(σ′.prm)
    @test τ.prm == reverse(τ′.prm)
    @test bandwidth(A[perm(σ), perm(τ)]) < bandwidth(A)
    @test bandwidth(permutedims(A[perm(σ), perm(τ)])) < bandwidth(permutedims(A))
    A = matrices["HB/west0132"]
    (σ, τ) = permute_plaid(A, CuthillMcKeePermuter())
    (σ′, τ′) = permute_plaid(A, CuthillMcKeePermuter(reverse=true))
    @test σ.prm == reverse(σ′.prm)
    @test τ.prm == reverse(τ′.prm)
    @test bandwidth(A[perm(σ), perm(τ)]) < bandwidth(A)
    @test bandwidth(permutedims(A[perm(σ), perm(τ)])) < bandwidth(permutedims(A))
    A = matrices["LPnetlib/lp_etamacro"]
    (σ, τ) = permute_plaid(A, CuthillMcKeePermuter())
    (σ′, τ′) = permute_plaid(A, CuthillMcKeePermuter(reverse=true))
    @test σ.prm == reverse(σ′.prm)
    @test τ.prm == reverse(τ′.prm)
    @test bandwidth(A[perm(σ), perm(τ)]) < bandwidth(A)
    @test bandwidth(permutedims(A[perm(σ), perm(τ)])) < bandwidth(permutedims(A))
end