text:
using QuBase
using QuDynamics
using BenchmarkLite

type SampleProc{Alg} <: Proc end

# Jaynes-Cummings type which takes the configuration parameters.
# These parameters are used for benchmarking, creating various
# JaynesCummings types with various parameters and using
# the solvers to get the timing results.
type JaynesCummings
    N::Int
    wc::Float64
    wa::Float64
    g::Float64
    kappa::Float64
    gamma::Float64
    use_rwa::Bool
    ntraj
    tlist
end

# Name of the solver, used in the creating the table of timing analysis
Base.string{Alg}(::SampleProc{Alg}) = lowercase(string(Alg))

# Constructing a valid configuration, this is set to `true` as we have
# set the type of the required fields in the construct `JaynesCummings` itself.
Base.isvalid{Alg}(p::SampleProc{Alg}, cfg::JaynesCummings) = true

# This returns the `hamitlonian`, `collapse operators`, `state drawn from QuMCWFEnsemble`,
# `array initialized for storing the evolved states` taking in the `cfg`
# configuration of the system.
function operator_set(p::SampleProc, cfg::JaynesCummings)
    idc = QuArray(eye(cfg.N))
    ida = QuArray(eye(2))
    a  = tensor(lowerop(cfg.N), ida)
    sm = tensor(idc, lowerop(2))
    if cfg.use_rwa
        H = cfg.wc * a' * a + cfg.wa * sm' * sm + cfg.g * (a' * sm + a * sm')
    else
        H = cfg.wc * a' * a + cfg.wa * sm' * sm + cfg.g * (a' + a) * (sm + sm')
    end
    c_op_list = Array(QuBase.AbstractQuMatrix, 0)
    n_th_a = 0.0
    rate_1 = cfg.kappa * (1 + n_th_a)
    if rate_1 > 0.0
        push!(c_op_list, full(sqrt(rate_1) * a))
    end
    rate = cfg.kappa * n_th_a
    if rate > 0.0
        push!(c_op_list, full(sqrt(rate) * a'))
    end
    rate = cfg.gamma
    if rate > 0.0
        push!(c_op_list, full(sqrt(rate) * sm))
    end
    psi = complex(tensor(statevec(1, FiniteBasis(cfg.N)), statevec(2, FiniteBasis(2))))
    rho = psi*psi'
    qumcwfen = QuMCWFEnsemble(complex(psi), cfg.ntraj)
    rhos = Array(typeof(rho), length(cfg.tlist)-1)
    for i=1:length(cfg.tlist)-1
        rhos[i] = complex(zeros(rho))
    end
    return full(H), c_op_list, qumcwfen, rhos
end

# Start states are pre-allocated, that is they are constructed before the evaluation
# here we return the `hamiltonian`, `collapse operator list`, `initial density matrix`,
# `observable whose expectation is to be calculated` using the above function
# i.e., operator_set. The first three remain fixed, the fourth can be varied as required.
Base.start(p::SampleProc, cfg::JaynesCummings) = (operator_set(p, cfg))

Base.length(p::SampleProc, cfg::JaynesCummings) = cfg.N

# This function gets the timing done for calculation of trace, using the pre-allocated
# constructions from the start method i.e., the hamiltonain, collapse operator list,
# initial density matrix and observable whose expectation is to be calculated.
function Base.run{Alg}(p::SampleProc{Alg}, cfg::JaynesCummings, s::(QuBase.QuArray, Array, QuDynamics.QuMCWFEnsemble, Array))
    for psi0 in s[3]
        i = 1
        for (t,psi) in QuPropagator(s[1], s[2], psi0, cfg.tlist, Alg())
            s[4][i] = s[4][i] + (psi*psi')/length(s[3])/norm(psi)^2
            i = i + 1
        end
    end
end

Base.done(p::SampleProc, cfg, s) = nothing

# Creating the array of methods on which benchmarks have to be performed
procs = Proc[ SampleProc{QuMCWF}()]

# Setting up various configurations.
cfgs = JaynesCummings[JaynesCummings(2, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, true, 200, 0.:0.1:2*pi),
                      JaynesCummings(3, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, true, 300, 0.:0.1:2*pi),
                      JaynesCummings(6, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, true, 400, 0.:0.1:2*pi),
                      JaynesCummings(10, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, true, 500, 0.:0.1:2*pi),
                      JaynesCummings(2, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, false, 200, 0.:0.1:2*pi),
                      JaynesCummings(3, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, false, 300, 0.:0.1:2*pi),
                      JaynesCummings(6, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, false, 400, 0.:0.1:2*pi),
                      JaynesCummings(10, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, false, 500, 0.:0.1:2*pi),
                      JaynesCummings(2, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, true, 200, 0.0:0.20134228187919462:30.0),
                      JaynesCummings(3, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, true, 300, 0.0:0.20134228187919462:30.0),
                      JaynesCummings(6, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, true, 400, 0.0:0.20134228187919462:30.0),
                      JaynesCummings(10, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, true, 500, 0.0:0.20134228187919462:30.0),
                      JaynesCummings(2, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, false, 200, 0.0:0.20134228187919462:30.0),
                      JaynesCummings(3, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, false, 300, 0.0:0.20134228187919462:30.0),
                      JaynesCummings(6, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, false, 400, 0.0:0.20134228187919462:30.0),
                      JaynesCummings(10, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, false, 500, 0.0:0.20134228187919462:30.0)]

# Generating the timing analysis.
rtable = run(procs, cfgs)

# Table display from the above method.
show(rtable; unit=:usec)


text:
__precompile__()
module MovingWeightedLeastSquares

const Point = Vector{T} where {T <: Real}

using DynamicPolynomials
using MultivariatePolynomials
using NearestNeighbors
using DataStructures

include("polynomial-generator.jl")
export WlsObject, wls, calcWlsCoefficients, plotWls
include("wls-object.jl")
include("wls-ops.jl")
export CellLinkedList, cllAdd!, cll_add!, cllRemove!, cll_remove!, cllModify!, cll_modify!, cllInrange, cll_inrange
include("cell-linked-list.jl")
export MwlsObject, mwlsNaive, mwls_naive, mwlsKd, mwls_kd, mwlsCll, mwls_cll
include("mwls-object.jl")
export calcMwlsCoefficients, mwls_coefficients,
       approximate, calcDiffMwlsPolys,
       mwls_diff_polys, mwlsDiff, mwls_diff
include("mwls-ops.jl")

end

ext: jl


text:
function RK=MethodRKIMEX(name)
switch name
  case 'ARK2'
    RK.nStage=3;
    RK.aE=zeros(RK.nStage,RK.nStage);
    RK.aI=zeros(RK.nStage,RK.nStage);
    RK.cE=zeros(RK.nStage,1);
    RK.cI=zeros(RK.nStage,1);
    RK.bE=zeros(1,RK.nStage);
    RK.bI=zeros(1,RK.nStage);
    RK.gamma=1-1/sqrt(2);
    alpha=(3+2*sqrt(2))/6;
    delta=1/(2*sqrt(2));
    gamma=1-1/sqrt(2);

    RK.gamma=gamma;
    RK.aE(2,1)=2-sqrt(2);
    RK.aE(3,1)=1-alpha;
    RK.aE(3,2)=alpha;
    RK.cE(2)=2-sqrt(2);
    RK.cE(3)=1;
    RK.bE(1)=delta;
    RK.bE(2)=delta;
    RK.bE(3)=gamma;

    RK.aI(2,1)=gamma;
    RK.aI(2,2)=gamma;
    RK.aI(3,1)=delta;
    RK.aI(3,2)=delta;
    RK.aI(3,3)=gamma;
    RK.cI=RK.cE;
    RK.bI=RK.bE;

    case 'HSDIRK2(2,2,2)'
    RK.nStage=2;
    RK.aE=zeros(RK.nStage,RK.nStage);
    RK.aI=zeros(RK.nStage,RK.nStage);
    RK.cE=zeros(RK.nStage,1);
    RK.cI=zeros(RK.nStage,1);
    RK.bE=zeros(1,RK.nStage);
    RK.bI=zeros(1,RK.nStage);
    RK.gamma=0.5;

    RK.aE(2,1)=1;
    RK.cE(2)=1;
    RK.bE(1)=.5;
    RK.bE(2)=.5;

    RK.aI(1,1)=.5;
    RK.aI(2,2)=.5;
    RK.cI(1)=.5;
    RK.cI(2)=.5;
    RK.bI=RK.bE;
  case default
    a='Keine Methode'
end
end


text:
##
# sample test program to produce monthly aggregate count of smileys
# used in tweets between year 2006 to 2009, based on data from infochimps.

using DataFrames
using Gaston
using HDFS
using HDFS.MapReduce
using Blocks

import HDFS.MapReduce.find_rec

function find_rec(rdr, iter_status, t::Type{DataFrame}, rec_sep::Char='\n', col_sep::Char=',')
    (iter_status != nothing) && return (nothing, true, iter_status)
    ios = HDFS.MapReduce.get_stream(rdr)
    sz = nb_available(ios)
    # frequent peek done by DataFrame reader causes libhdfs to misbehave
    buff = Array(Uint8, sz)
    read(ios, buff)
    pb = PipeBuffer(buff)
    #println("bytes available: $sz")
    #df = readtable(pb, header=false, separator=col_sep, nrows=sz, buffersize=sz)
    df = as_dataframe(pb, header=false, separator=col_sep)
    return (df, false, position(ios))

    #(dlmarr, iseof, iter_status) = find_rec(rdr, iter_status, Matrix, rec_sep, col_sep)
    #df = (nothing == dlmarr) ? nothing : DataFrame(dlmarr)
    #println("hdfs_find_rec_table got $dlmarr")
    #(df, iseof, iter_status)
end

##
# find smiley records from HDFS CSV file
find_smiley_df(r, next_rec_pos) = find_rec(r, next_rec_pos, DataFrame, '\n', '\t')

##
# reduce smiley counts or array of counts
reduce_smiley_df(reduced, results...) = +((nothing == reduced) ? 0 : reduced, filter(x->(x!=nothing),results)...)

##
# for finding total counts across all years
function map_total(rec)
    #println("map got $rec")
    (nothing == rec) && return []
    [sum(rec[3])]
end

function collect_total(results, rec)
    #println("collect got $rec")
    (length(rec) == 0) && return results
    (results == nothing) ? rec : (results+rec)
end

function do_dataframe_test(furl::String, use_hdfs=true)
    println("starting dmapreduce job...")
    finp = use_hdfs ? MRHdfsFileInput([furl], find_smiley_df) : MRFsFileInput([furl], find_smiley_df, 1024*16)
    j = dmapreduce(finp, map_total, collect_total, reduce_smiley_df)
    println("waiting for jobs to finish...")
    loopstatus = true
    while(loopstatus)
        sleep(2)
        jstatus,jstatusinfo = status(j,true)
        ((jstatus == "error") || (jstatus == "complete")) && (loopstatus = false)
        (jstatus == "running") && println("$(j): $(jstatusinfo)% complete...")
    end

    wait(j)
    println("time taken (total time, wait time, run time): $(times(j))")
    println("")
    jstatus,total_count = results(j)
    println("completed with status: $jstatus")
    if(jstatus == "complete")
        println("results:")
        println("total number of smileys: $total_count")
    else
        println(status(j, true))
    end
    unload(j) # ensure file handles are closed cleanly
end


text:
###Cargar datos
using TopicModelsVB
import Distributions.sample
#Arreglar direccion
corp = readcorp(:mac)

corp.docs = vcat([sample(filter(doc -> round(doc.stamp / 100) == y, corp.docs), 400, replace=false) for y in 1984:2005]...)

fixcorp!(corp, abr=100, len=10) # Remove words which appear < 100 times and documents of length < 10.

basemodel = LDA(corp, 9)
train!(basemodel, iter=150, chkelbo=151)

# training...

model = DTM(corp, 9, 200, basemodel)
train!(model, iter=10) # This will likely take about an hour on a personal computer.
                       # Convergence for all other models is worst-case quadratic,
                       # while DTM convergence is linear or at best super-linear.
# training...


text:
using .AlgebraOfGraphics

include("AoG/plot_chains.jl")


text:
"""
    multisurf(data::Array{<:Real,2}, target::Array{<:Integer,1}, m::Signed=-1,
              dist_func::Any=(e1, e2) -> sum(abs.(e1 .- e2), dims=2);
              f_type::String="continuous")::Array{Float64,1}

Compute feature weights using MultiSURF algorithm. The f_type argument specifies whether the features are continuous or discrete
and can either have the value of "continuous" or "discrete".

---
# Reference:
- Ryan Urbanowicz, Randal Olson, Peter Schmitt, Melissa Meeker, and
Jason Moore. Benchmarking Relief-based feature selection methods for
bioinformatics data mining. Journal of Biomedical Informatics, 85, 2018.
"""
function multisurf(data::Array{<:Real,2}, target::Array{<:Integer,1}, m::Signed=-1,
                   dist_func::Any=(e1, e2) -> sum(abs.(e1 .- e2), dims=2);
                   f_type::String="continuous")::Array{Float64,1}

    # Initialize feature weights vector.
    weights = zeros(Float64, 1, size(data, 2))

    # Compute vectors of maximum and minimum feature values.
    max_f_vals = maximum(data, dims=1)
    min_f_vals = minimum(data, dims=1)

    # Sample m samples without replacement. If m has signal value -1, use all samples.
    sample_idxs = StatsBase.sample(1:size(data, 1), if (m==-1) size(data,1) else m end, replace=false)
    if (m == -1) m = size(data, 1) end # If m has signal value -1, set m to total number of examples.

    # Compute pairwise distances between samples (vector form).
    dists = Array{Float64}(undef, Int64((size(data, 1)^2 - size(data, 1))/2 + 1))
    dists[1] = 0  # Set first value of distances vector to 0 - accessed when i == j in square form indices.

    # Construct pairwise distances vector using vectorized distance function.
    top_ptr = 2
    @inbounds for idx = 1:size(data,1)-1
        upper_lim = top_ptr + size(data, 1) - idx - 1
        dists[top_ptr:upper_lim] = dist_func(data[idx:idx, :], data[idx+1:end, :])
        top_ptr = upper_lim + 1
    end

    # Go over sampled samples.
    @inbounds for idx = sample_idxs

        # Row and column indices for querying pairwise distance vector.
        row_idxs = repeat([idx - 1], size(data, 1))
        col_idxs = collect(0:size(data, 1)-1)

        # Get indices in distance vector (from square form indices).
        neigh_idx = square_to_vec(row_idxs, col_idxs, size(data, 1)) .+ 2

        # Query distances to neighbors to get masks for both zones.
        dists_neighbors = dists[neigh_idx[neigh_idx .!= 1]]
        mu = Statistics.mean(dists_neighbors)
        sig = Statistics.std(dists_neighbors)
        thresh_near = mu - sig/2.0
        neigh_mask_near = dists_neighbors .< thresh_near
        insert!(neigh_mask_near, idx, 0)

        # Get class values of miss neighbors.
        miss_classes_near = target[neigh_mask_near .& (target .!= target[idx])]

        # Get masks for considered regions.
        hit_neigh_mask_near = neigh_mask_near .& (target .== target[idx])
        miss_neigh_mask_near = neigh_mask_near .& (target .!= target[idx])

        # Ignore samples that have no hits or misses in radius.
        if sum(hit_neigh_mask_near) == 0 || sum(miss_neigh_mask_near) == 0
            continue
        end

        # Compute weights for near misses and compute weighting vector.
        weights_mult = Array{Float64}(undef, length(miss_classes_near))  # Allocate weights multiplier vector.
        cm = countmap(miss_classes_near)                                 # Count unique values.
        u = collect(keys(cm))
        c = collect(values(cm))
        neighbor_weights = c ./ length(miss_classes_near)  # Compute misses' weights.

        @inbounds for (w, val) = zip(neighbor_weights, u)  # Build multiplier vector.
            find_res = findall(miss_classes_near .== val)
            weights_mult[find_res] .= w
        end

        ### Weights Update ###

        if f_type == "continuous"
            # If features continuous.

            # Penalty term
            penalty = sum(abs.(data[idx:idx, :] .- data[hit_neigh_mask_near, :])./((max_f_vals .- min_f_vals) .+ eps(Float64)), dims=1)

            # Reward term
            reward = sum(weights_mult .* (abs.(data[idx:idx, :] .- data[miss_neigh_mask_near, :])./((max_f_vals .- min_f_vals .+ eps(Float64)))), dims=1)

            # Weights update
            weights = weights .- penalty./(size(data, 1)*sum(hit_neigh_mask_near) + eps(Float64)) .+
                reward./(size(data, 1)*sum(miss_neigh_mask_near) + eps(Float64))

        elseif f_type == "discrete"
            # If features discrete.

            # Penalty term
            penalty = sum(Int64.(data[idx:idx, :] .!= data[hit_neigh_mask_near, :]), dims=1)

            # Reward term
            reward = sum(weights_mult .* (Int64.(data[idx:idx, :] .!= data[miss_neigh_mask_near, :])), dims=1)

            # Weights update
            weights = weights .- penalty./(size(data, 1)*sum(hit_neigh_mask_near) + eps(Float64)) .+
                reward./(size(data, 1)*sum(miss_neigh_mask_near) + eps(Float64))

        else
            throw(DomainError(f_type, "f_type can only be equal to \"continuous\" or \"discrete\"."))
        end

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

    end

    # Return computed feature weights.
    return vec(weights)

end


text:
# Autogenerated wrapper script for SoapyLMS7_jll for aarch64-linux-gnu-cxx11
export LimeUtil, libLMS7Support, libLimeSuite

using soapysdr_jll
using libusb_jll
JLLWrappers.@generate_wrapper_header("SoapyLMS7")
JLLWrappers.@declare_executable_product(LimeUtil)
JLLWrappers.@declare_library_product(libLMS7Support, "libLMS7Support.so")
JLLWrappers.@declare_library_product(libLimeSuite, "libLimeSuite.so.20.10-1")
function __init__()
    JLLWrappers.@generate_init_header(soapysdr_jll, libusb_jll)
    JLLWrappers.@init_executable_product(
        LimeUtil,
        "bin/LimeUtil",
    )

    JLLWrappers.@init_library_product(
        libLMS7Support,
        "lib/SoapySDR/modules0.8/libLMS7Support.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libLimeSuite,
        "lib/libLimeSuite.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# box_test.jl
# Box-Pierce and Ljung-Box tests for autocorrelation
#
# Copyright (C) 2017   Benjamin Born
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

export BoxPierceTest, LjungBoxTest

# Box-Pierce test

struct BoxPierceTest <: HypothesisTest
    y::Vector{Float64}  # time series vector
    n::Int              # number of observations
    lag::Int            # number of lags in test statistic
    dof::Int            # degrees-of-freedom correction (if y is residual)
    Q::Float64          # test statistic
end

"""
    BoxPierceTest(y, lag, dof=0)

Compute the Box-Pierce `Q` statistic to test the null hypothesis of
independence in a time series `y`.

`lag` specifies the number of lags used in the construction of `Q`. When
testing the residuals of an estimated model, `dof` has to be set to the number
of estimated parameters. E.g., when testing the residuals of an ARIMA(p,0,q)
model, set `dof=p+q`.

# External links

  * [Box-Pierce test on Wikipedia](https://en.wikipedia.org/wiki/Ljung–Box_test#Box-Pierce_test)
"""
function BoxPierceTest(y::AbstractVector{T}, lag::Int, dof::Int=0) where T<:Real
    if dof>=lag
        throw(ArgumentError("Number of lags has to be larger than degrees of" *
        " freedom correction"))
    end
    n = size(y,1)
    Q = n * first(sum(k -> autocor(y,k:k).^2, 1:lag))
    BoxPierceTest(y,n,lag,dof,Q)
end

testname(::BoxPierceTest) = "Box-Pierce autocorrelation test"
population_param_of_interest(x::BoxPierceTest) = ("autocorrelations up to lag k",
    "all zero", NaN)
default_tail(test::BoxPierceTest) = :right

function show_params(io::IO, x::BoxPierceTest, ident)
    println(io, ident, "number of observations:         ", x.n)
    println(io, ident, "number of lags:                 ", x.lag)
    println(io, ident, "degrees of freedom correction:  ", x.dof)
    println(io, ident, "Q statistic:                    ", x.Q)
end

pvalue(x::BoxPierceTest) = pvalue(Chisq(x.lag-x.dof), x.Q; tail=:right)

#Ljung-Box test

struct LjungBoxTest <: HypothesisTest
    y::Vector{Float64}  # time series vector
    n::Int              # number of observations
    lag::Int            # number of lags in test statistic
    dof::Int            # degrees-of-freedom correction (if y is residual)
    Q::Float64          # test statistic
end

"""
    LjungBoxTest(y, lag, dof=0)

Compute the Ljung-Box `Q` statistic to test the null hypothesis of
independence in a time series `y`.

`lag` specifies the number of lags used in the construction of `Q`. When
testing the residuals of an estimated model, `dof` has to be set to the number
of estimated parameters. E.g., when testing the residuals of an ARIMA(p,0,q)
model, set `dof=p+q`.

# External links

  * [Ljung-Box test on Wikipedia](https://en.wikipedia.org/wiki/Ljung–Box_test)
"""
function LjungBoxTest(y::AbstractVector{T}, lag::Int, dof::Int=0) where T<:Real
    if dof>=lag
        throw(ArgumentError("Number of lags has to be larger than degrees of" *
        " freedom correction"))
    end
    n = size(y,1)
    Q = n*(n+2)* first(sum(k -> autocor(y,k:k).^2/(n-k), 1:lag))
    LjungBoxTest(y,n,lag,dof,Q)
end

testname(::LjungBoxTest) = "Ljung-Box autocorrelation test"
population_param_of_interest(x::LjungBoxTest) = ("autocorrelations up to lag k",
    "all zero", NaN)
default_tail(test::LjungBoxTest) = :right

function show_params(io::IO, x::LjungBoxTest, ident)
    println(io, ident, "number of observations:         ", x.n)
    println(io, ident, "number of lags:                 ", x.lag)
    println(io, ident, "degrees of freedom correction:  ", x.dof)
    println(io, ident, "Q statistic:                    ", x.Q)
end

pvalue(x::LjungBoxTest) = pvalue(Chisq(x.lag-x.dof), x.Q; tail=:right)


text:
module TestPotentials1D
# a module of simple potentials for testing optimization and sampling methods

"""
SymmetricDoubleWell - The classic symmetric doublewell potential.
"""
function SymmetricDoubleWell(x)
    return (x^2 -1)^2;
end

"""
ASymmetricDoubleWell - An asymmetric double well potential requires |δ|< 1 for
multiple minima.
"""
function ASymmetricDoubleWell(x;δ=0.5)
    return x^4 - 1.5 * x^2 - δ * x;
end

"""
FatSkinnyDoubleWell10 - The fat-skinny double well of degree 10. This introduces
entropic effects in dimension 1.
"""
function FatSkinnyDoubleWell10(x)
    return ((8-5*x)^8 * (2+5*x)^2)/(2^26);
end

end # end module

max_stars_repo_name: gideonsimpson/TestLandscapes

max_issues_repo_name: gideonsimpson/TestLandscapes

max_forks_repo_name: gideonsimpson/TestLandscapes


text:
@testset "training.jl inside the solver (with backtracking)" begin

    println("Training inside the solver -------- ")

    generator = SeaPearl.HomogenousGraphColoringGenerator(10, 0.1)

    numInFeatures = 3
    # Model definition
    approximator_GNN = SeaPearl.GraphConv(64 => 64, Flux.leakyrelu)
    target_approximator_GNN = SeaPearl.GraphConv(64 => 64, Flux.leakyrelu)
    gnnlayers = 10

    approximator_model = SeaPearl.CPNN(
        graphChain = Flux.Chain(
            SeaPearl.GraphConv(numInFeatures => 64, Flux.leakyrelu),
            [approximator_GNN for i = 1:gnnlayers]...
        ),
        nodeChain = Flux.Chain(
            Flux.Dense(64, 32, Flux.leakyrelu),
            Flux.Dense(32, 32, Flux.leakyrelu),
            Flux.Dense(32, 16, Flux.leakyrelu),
        ),
        outputChain = Flux.Dense(16, generator.nb_nodes),
    )
    target_approximator_model = SeaPearl.CPNN(
        graphChain = Flux.Chain(
            SeaPearl.GraphConv(numInFeatures => 64, Flux.leakyrelu),
            [target_approximator_GNN for i = 1:gnnlayers]...
        ),
        nodeChain = Flux.Chain(
            Flux.Dense(64, 32, Flux.leakyrelu),
            Flux.Dense(32, 32, Flux.leakyrelu),
            Flux.Dense(32, 16, Flux.leakyrelu),
        ),
        outputChain = Flux.Dense(16, generator.nb_nodes),
    )

    # Agent definition
    agent = RL.Agent(
        policy = RL.QBasedPolicy(
            learner = RL.DQNLearner(
                approximator = RL.NeuralNetworkApproximator(
                    model = approximator_model,
                    optimizer = ADAM(0.0005f0)
                ),
                target_approximator = RL.NeuralNetworkApproximator(
                    model = target_approximator_model,
                    optimizer = ADAM(0.0005f0)
                ),
                loss_func = Flux.Losses.huber_loss,
                stack_size = nothing,
                γ = 0.99f0,
                batch_size = 32,
                update_horizon = 4,
                min_replay_history = 32,
                update_freq = 4,
                target_update_freq = 20,
            ),
            explorer = RL.EpsilonGreedyExplorer(
                ϵ_stable = 0.001,
                kind = :exp,
                ϵ_init = 1.0,
                warmup_steps = 0,
                decay_steps = 5000,
                step = 1,
                is_break_tie = false,
                #is_training = true,
                rng = MersenneTwister(33)
            )
        ),
        trajectory = RL.CircularArraySARTTrajectory(
            capacity = 800,
            state = SeaPearl.DefaultTrajectoryState[] => (),
        )
    )

    learnedHeuristic = SeaPearl.LearnedHeuristic(agent)

    initial_params = deepcopy(params(learnedHeuristic.agent.policy.learner.approximator.model))

    metrics,evalmetrics = SeaPearl.train!(
        valueSelectionArray=learnedHeuristic,
        generator=generator,
        nbEpisodes= 10, #3,
        strategy=SeaPearl.DFSearch(),
        variableHeuristic=SeaPearl.MinDomainVariableSelection{true}(),
        out_solver = false,
        verbose = false,
        evaluator=SeaPearl.SameInstancesEvaluator([learnedHeuristic],generator)  #need to be a vector of Heuristic
    )

    final_params = params(learnedHeuristic.agent.policy.learner.approximator.model)

    @test final_params != initial_params

end

@testset "training.jl outside the solver (without backtracking)" begin

    println("Let's try without backtracking   :   ------------")

    generator = SeaPearl.HomogenousGraphColoringGenerator(10, 0.1)

    numInFeatures = 3

    # Model definition
    approximator_GNN = SeaPearl.GraphConv(64 => 64, Flux.leakyrelu)
    target_approximator_GNN = SeaPearl.GraphConv(64 => 64, Flux.leakyrelu)
    gnnlayers = 10

    approximator_model = SeaPearl.CPNN(
        graphChain = Flux.Chain(
            SeaPearl.GraphConv(numInFeatures => 64, Flux.leakyrelu),
            [approximator_GNN for i = 1:gnnlayers]...
        ),
        nodeChain = Flux.Chain(
            Flux.Dense(64, 32, Flux.leakyrelu),
            Flux.Dense(32, 32, Flux.leakyrelu),
            Flux.Dense(32, 16, Flux.leakyrelu),
        ),
        outputChain = Flux.Dense(16, generator.nb_nodes),
    )
    target_approximator_model = SeaPearl.CPNN(
        graphChain = Flux.Chain(
            SeaPearl.GraphConv(numInFeatures => 64, Flux.leakyrelu),
            [target_approximator_GNN for i = 1:gnnlayers]...
        ),
        nodeChain = Flux.Chain(
            Flux.Dense(64, 32, Flux.leakyrelu),
            Flux.Dense(32, 32, Flux.leakyrelu),
            Flux.Dense(32, 16, Flux.leakyrelu),
        ),
        outputChain = Flux.Dense(16, generator.nb_nodes),
    )

    # Agent definition
    agent = RL.Agent(
        policy = RL.QBasedPolicy(
            learner = RL.DQNLearner(
                approximator = RL.NeuralNetworkApproximator(
                    model = approximator_model,
                    optimizer = ADAM(0.0005f0)
                ),
                target_approximator = RL.NeuralNetworkApproximator(
                    model = target_approximator_model,
                    optimizer = ADAM(0.0005f0)
                ),
                loss_func = Flux.Losses.huber_loss,
                stack_size = nothing,
                γ = 0.99f0,
                batch_size = 32,
                update_horizon = 4,
                min_replay_history = 32,
                update_freq = 4,
                target_update_freq = 20,
            ),
            explorer = RL.EpsilonGreedyExplorer(
                ϵ_stable = 0.001,
                kind = :exp,
                ϵ_init = 1.0,
                warmup_steps = 0,
                decay_steps = 5000,
                step = 1,
                is_break_tie = false,
                #is_training = true,
                rng = MersenneTwister(33)
            )
        ),
        trajectory = RL.CircularArraySLARTTrajectory(
            capacity = 800,
            state = SeaPearl.DefaultTrajectoryState[] => (),
            legal_actions_mask = Vector{Bool} => (generator.nb_nodes, ),

        )
    )

    learnedHeuristic = SeaPearl.LearnedHeuristic(agent)
    initial_params = deepcopy(params(learnedHeuristic.agent.policy.learner.approximator.model))

    metrics,evalmetrics = SeaPearl.train!(
        valueSelectionArray=learnedHeuristic,
        generator=generator,
        nbEpisodes=10,
        strategy=SeaPearl.DFSearch(),
        variableHeuristic=SeaPearl.MinDomainVariableSelection{true}(),
        out_solver = true,
        verbose = false,
        evaluator=nothing
    )

    final_params = params(learnedHeuristic.agent.policy.learner.approximator.model)

    @test final_params != initial_params

end


text:
using VLKeggSDK
using DataFrames

# get codes -
ca = (get_kegg_organism_codes() |> check) |> DataFrame


text:
#
# Functions that return the atoms of the second set that are closer to
# each molecule of the first set (only one list is returned).
#

function update_list_cross!(
    i,
    j,
    d2,
    mol_index_i::F,
    list::List,
) where {F<:Function}
    d = sqrt(d2)
    imol = mol_index_i(i)
    if d < list[imol].d
        list[imol] = MinimumDistance(true, i, j, d)
    end
    return list
end

"""

```
function minimum_distances!(
    mol_index_i::F,
    x_list::AbstractVector{<:MinimumDistance},
    y::AbstractVector,
    box, cl;
    parallel = true
) where {F<:Function}
```

Compute the list of minimum distances given the precomputed cell lists, auxiliary vectors, and
the functions that associate each atom to the corresponding molecule. Should be preferred for multiple calls of this function,
with the outer update of the cell lists.

## Example

```
julia> x = [ rand(SVector{3,Float64}) for _ in 1:12 ];

julia> y = [ rand(SVector{2,Float64}) for _ in 1:800 ];

julia> using CellListMap

julia> box = Box([1,1],0.2);

julia> cl = CellList(x,y,box);

julia> x_list = init_list(x, i -> mol_index(i,4)) # 4 atoms per molecule
3-element Vector{MinimumDistance{Float64}}:
 MinimumDistance{Float64}(false, 0, 0, Inf)
 MinimumDistance{Float64}(false, 0, 0, Inf)
 MinimumDistance{Float64}(false, 0, 0, Inf)

julia> x_list = init_list(y, i -> mol_index(i,10)); # 10 atoms per molecule

julia> x_list_threaded = [ copy(y_list) for _ in 1:nbatches(cl) ]

julia> minimum_distances!(
           i -> mol_index(i,4),
           x_list,
           y,
           box,
           cl;
           list_threaded = y_list_threaded
       );

julia> x_list
3-element Vector{MinimumDistance{Float64}}:
 MinimumDistance{Float64}(true, 1, 516, 0.007047201657257146)
 MinimumDistance{Float64}(true, 7, 415, 0.007128399163692169)
 MinimumDistance{Float64}(true, 9, 511, 0.0033074993141050577)

```

"""
function minimum_distances!(
    mol_index_i::F,
    x_list::List,
    box::Box,
    cl::CellListMap.CellListPair;
    parallel = true,
    list_threaded = nothing,
) where {F<:Function}
    reset!(x_list)
    reset!(list_threaded)
    map_pairwise!(
        (x, y, i, j, d2, x_list) -> update_list_cross!(i, j, d2, mol_index_i, x_list),
        x_list,
        box,
        cl;
        parallel = parallel,
        reduce = reduce_list!,
        output_threaded = list_threaded
    )
    return x_list
end

"""

```
function minimum_distances(
    x, y,
    n_atoms_per_molecule_x,
    box;
    parallel = true
)
```

Compute the list of minimum distances given the vector of atomic coordinates `x` and `y`,
the number of atoms of the molecules, the box (of type `CellListMap.Box`). This is the simplest
interface, but it allocates new arrays. For iterative computations use `minimum_distances!`
instead.

## Example

```
julia> x = [ rand(SVector{3,Float64}) for _ in 1:100 ];

julia> y = [ rand(SVector{3,Float64}) for _ in 1:90 ];

julia> box =  Box([1,1,1],0.2);

julia> x_list, y_list = minimum_distances(x, y, 5, 3, box);

julia> x_list
20-element Vector{MinimumDistance{Float64}}:
 MinimumDistance{Float64}(true, 4, 39, 0.11311002375318155)
 MinimumDistance{Float64}(true, 9, 42, 0.0752872244199588)
 MinimumDistance{Float64}(true, 12, 8, 0.06411501358230641)
 ⋮
 MinimumDistance{Float64}(true, 92, 88, 0.0468461712318602)
 MinimumDistance{Float64}(true, 96, 50, 0.09478144847155845)

```

"""
function minimum_distances(
    x::AbstractVector{<:AbstractVector{T}},
    y::AbstractVector{<:AbstractVector{T}},
    n_atoms_per_molecule_x::Int,
    box::Box{UnitCellType,N,T};
    parallel = true,
) where {UnitCellType,N,T}
    cl = CellList(x, y, box, parallel = parallel)
    x_list = init_list(x, i -> mol_index(i, n_atoms_per_molecule_x))
    minimum_distances!(
        i -> mol_index(i, n_atoms_per_molecule_x),
        x_list,
        box,
        cl;
        parallel = parallel,
    )
    return x_list
end

"""

```
function minimum_distances(
    x::AbstractVector{<:AbstractVector{T}},
    y::AbstractVector{<:AbstractVector{T}},
    n_atoms_per_molecule_x::Int,
    cutoff::T;
    parallel = true,
) where {T}
```

Input only the cutoff, assuming that no periodic boundary conditions are used.

"""
function minimum_distances(
    x::AbstractVector{<:AbstractVector{T}},
    y::AbstractVector{<:AbstractVector{T}},
    n_atoms_per_molecule_x::Int,
    cutoff::T;
    parallel = true,
) where {T}
    box = Box(limits(x,y),cutoff)
    cl = CellList(x, y, box, parallel = parallel)
    x_list = init_list(x, i -> mol_index(i, n_atoms_per_molecule_x))
    minimum_distances!(
        i -> mol_index(i, n_atoms_per_molecule_x),
        x_list,
        box,
        cl;
        parallel = parallel,
    )
    return x_list
end


text:
module Gallery

using CellwiseValues
using Base: @propagate_inbounds, @pure

export CellValueFromArray
import Base: iterate
import Base: length
import Base: eltype
import Base: size
import Base: getindex
import Base: IndexStyle

struct CellValueFromArray{T,N,V<:AbstractArray{T,N}} <: IndexCellValue{T,N}
  v::V
end

@propagate_inbounds function getindex(
  self::CellValueFromArray{T,N},I::Vararg{Int,N}) where {T,N}
  @inbounds self.v[I...]
end

size(self::CellValueFromArray) = size(self.v)

@pure IndexStyle(::Type{CellValueFromArray{T,N,V}}) where {T,N,V} = IndexStyle(V)

end # module Gallery


text:
function scale!(
  df::DataFrame,
  vars::Vector{Symbol},
  ext="_s")

  for var in vars
    mean_var = mean(df[!, var])
    std_var = std(df[!, var])

    df[!, Symbol("$(String(var))$ext")] =
      (df[:, var] .- mean_var)/std_var
  end
  df
end

function scale!(
  df::DataFrame,
  var::Symbol,
  ext="_s")

  mean_var = mean(df[!, var])
  std_var = std(df[!, var])

  df[!, Symbol("$(String(var))$ext")] =
    (df[:, var] .- mean_var)/std_var

end


text:
"""
    upload(bucket::Bucket, filename::AbstractString, source::AbstractMongoStream;
        options::Union{Nothing, BSON}=nothing,
        file_id=BSONObjectId())

Uploads data from `source` to a GridFS file `filename`.

`bson_file_id` is a BSON document with a `_id` field.
If `_id` field is not present, a `BSONObjectId` will be generated.
"""
function upload(bucket::Bucket, filename::AbstractString, source::AbstractMongoStream;
        file_id=BSONObjectId(),
        chunk_size::Union{Nothing, Integer}=nothing,
        metadata::Union{Nothing, BSON}=nothing)

    # Julia does not support C unions.
    # libmongoc represents bson_value_t as a C union.
    # To make things work, here we construct a BSON document
    # and set the id. Later, we retrieve a bson_value_t
    # using get_as_bson_value, which uses bson_iter_value.
    # new_bson_with_id returns a tuple with the BSON doc
    # along with the id to ensure that the lifetime of the BSON
    # doc exceeds the lifetime of the bson_value_t
    # that we'll pass to mongoc_gridfs_bucket_upload_from_stream_with_id.

    function new_bson_with_id(id)
        doc = BSON("_id" => id)
        return doc, get_as_bson_value(doc, "_id")
    end
    doc, file_id_as_bson_value = new_bson_with_id(file_id)

    function new_bson_upload_opts(chunk_size::Union{Nothing, Integer}, metadata::Union{Nothing, BSON}) :: Union{Nothing, BSON}
        if chunk_size == nothing && metadata == nothing
            return nothing
        else
            result = BSON()

            if chunk_size != nothing
                result["chunkSizeBytes"] = Int32(chunk_size)
            end

            if metadata != nothing
                result["metadata"] = metadata
            end

            return result
        end
    end
    options = new_bson_upload_opts(chunk_size, metadata)
    options_handle = options == nothing ? C_NULL : options.handle

    err_ref = Ref{BSONError}()

    ok = mongoc_gridfs_bucket_upload_from_stream_with_id(
            bucket.handle,
            file_id_as_bson_value.handle,
            filename,
            source.handle,
            options_handle,
            err_ref
        )

    if !ok
        throw(err_ref[])
    end

    nothing
end

"""
    upload(bucket::Bucket, remote_filename::AbstractString, local_source_filepath::AbstractString;
        options::Union{Nothing, BSON}=nothing,
        file_id=BSONObjectId())

High-level interface to upload a local file to a GridFS bucket.
"""
function upload(bucket::Bucket, remote_filename::AbstractString, local_source_filepath::AbstractString;
        file_id=BSONObjectId(),
        chunk_size::Union{Nothing, Integer}=nothing,
        metadata::Union{Nothing, BSON}=nothing)

    @assert isfile(local_source_filepath) "Couldn't find $local_source_filepath."
    source_stream = MongoStreamFile(local_source_filepath)

    try
        upload(bucket, remote_filename, source_stream,
            file_id=file_id, chunk_size=chunk_size, metadata=metadata)
    finally
        close(source_stream)
    end
end

"""
    download(bucket::Bucket, file_id, target::AbstractMongoStream)

Download a GridFS file identified by `file_id` to `target` stream.

`file_id` may be either a `BSONValue` or a `BSON` document with an `_id` field.
"""
function download(bucket::Bucket, file_id::BSONValue, target::AbstractMongoStream)
    err_ref = Ref{BSONError}()
    ok = mongoc_gridfs_bucket_download_to_stream(
            bucket.handle, file_id.handle, target.handle, err_ref
        )

    if !ok
        throw(err_ref[])
    end
    nothing
end

function download(bucket::Bucket, file_info::BSON, target::AbstractMongoStream)
    @assert haskey(file_info, "_id")
    download(bucket, get_as_bson_value(file_info, "_id"), target)
end

function find_file_info(bucket::Bucket, filename::AbstractString) :: BSON
    bson_filter = BSON()
    bson_filter["filename"] = String(filename)

    local found::Bool=false
    local result::BSON

    for doc in find(bucket, bson_filter)
        @assert !found "More than 1 file was found with name $filename."
        found = true

        result = doc
    end
    @assert found "Coudln't find file with name $filename."

    return result
end

function download(bucket::Bucket, remote_file_id_or_info::Union{BSONValue, BSON}, local_filepath::AbstractString;
        flags::Integer=(Base.Filesystem.JL_O_CREAT | Base.Filesystem.JL_O_RDWR), mode::Integer=0o600)

    download_file_stream = Mongoc.MongoStreamFile(local_filepath, flags=(Base.Filesystem.JL_O_CREAT | Base.Filesystem.JL_O_RDWR), mode=0o600)

    try
        Mongoc.download(bucket, remote_file_id_or_info, download_file_stream)
    finally
        Mongoc.close(download_file_stream)
    end
end

"""
    download(bucket::Bucket, remote_file, local_filepath::AbstractString;
        flags::Integer=(Base.Filesystem.JL_O_CREAT | Base.Filesystem.JL_O_RDWR), mode::Integer=0o600

Downloads a GridFS file `remote_file` to `local_filepath`.

`remote_file` can be either a unique filename, or a `file_id::BSONValue`, or `file_info::BSON`.
"""
function download(bucket::Bucket, remote_filename::AbstractString, local_filepath::AbstractString;
        flags::Integer=(Base.Filesystem.JL_O_CREAT | Base.Filesystem.JL_O_RDWR), mode::Integer=0o600)

    file_info = find_file_info(bucket, remote_filename)
    download(bucket, file_info, local_filepath, flags=flags, mode=mode)
end

"""
    find(bucket::Bucket, bson_filter::BSON=BSON();
        options::BSON=BSON()) :: Cursor

Looks for files in GridFS bucket.
"""
function find(bucket::Bucket, bson_filter::BSON=BSON();
        options::BSON=BSON()) :: Cursor

    cursor_handle = mongoc_gridfs_bucket_find(bucket.handle, bson_filter.handle, options.handle)

    if cursor_handle == C_NULL
        error("Couldn't execute query.")
    end
    return Cursor(bucket, cursor_handle)
end

"""
    delete(bucket::Bucket, file_id)
    delete(bucket::Bucket, file_metadata::BSON)

Deletes a file from a GridFS Bucket.
"""
function delete(bucket::Bucket, file_id::BSONValue)
    err_ref = Ref{BSONError}()
    ok = mongoc_gridfs_bucket_delete_by_id(bucket.handle, file_id.handle, err_ref)
    if !ok
        throw(err_ref[])
    end
    nothing
end

function delete(bucket::Bucket, file_metadata::BSON)
    oid_as_bson_value = get_as_bson_value(file_metadata, "_id")
    delete(bucket, oid_as_bson_value)
end

function delete(bucket::Bucket, file_id)
    file_metadata = BSON("_id" => file_id)
    delete(bucket, file_metadata)
end

function Base.empty!(bucket::Bucket)
    for doc in find(bucket)
        delete(bucket, doc)
    end
end

function Base.isempty(bucket::Bucket)
    for _ in find(bucket)
        return false
    end
    return true
end

function open_download_stream(bucket::Bucket, file_id::BSONValue;
        timeout_msec::Integer=DEFAULT_TIMEOUT_MSEC,
        chunk_size::Integer=DEFAULT_CHUNK_SIZE) :: MongoIOStream

    err_ref = Ref{BSONError}()
    stream_handle = mongoc_gridfs_bucket_open_download_stream(bucket.handle, file_id.handle, err_ref)
    if stream_handle == C_NULL
        throw(err_ref[])
    end
    return MongoIOStream(bucket, stream_handle; timeout_msec=timeout_msec, chunk_size=chunk_size)
end

function open_download_stream(bucket::Bucket, filename::AbstractString) :: MongoIOStream
    local result::MongoIOStream
    file_info = find_file_info(bucket, filename)
    return open_download_stream(bucket, file_info)
end

function open_download_stream(bucket::Bucket, file_info::BSON;
        timeout_msec::Integer=DEFAULT_TIMEOUT_MSEC,
        chunk_size::Integer=DEFAULT_CHUNK_SIZE) :: MongoIOStream

    @assert haskey(file_info, "_id")
    return open_download_stream(bucket, get_as_bson_value(file_info, "_id"),
                timeout_msec=timeout_msec,
                chunk_size=chunk_size)
end

"""
    open_download_stream(f, bucket, filename)
    open_download_stream(bucket, filename) :: MongoIOStream

Opens a stream for reading a remote file identified by `filename`.

# Example

Open a stream, do some work, then close it.

```julia
io = Mongoc.open_download_stream(bucket, file)
try
    tmp_str = read(io, String)
finally
    close(io)
end
```

Use *do-syntax* to ensure that the `io` stream will
be closed in case something goes wrong.

```julia
Mongoc.open_download_stream(bucket, remote_filename) do io
    tmp_str = read(io, String)
end
```
"""
function open_download_stream(f::Function, bucket::Bucket, file::Union{BSONValue, AbstractString})
    stream = open_download_stream(bucket, file)
    try
        f(stream)
    finally
        close(stream)
    end
end

function open_upload_stream(bucket::Bucket, file_id::BSONValue, filename::AbstractString;
        options::Union{Nothing, BSON}=nothing,
        timeout_msec::Integer=DEFAULT_TIMEOUT_MSEC,
        chunk_size::Integer=DEFAULT_CHUNK_SIZE) :: MongoIOStream

    err_ref = Ref{BSONError}()
    options_handle = options == nothing ? C_NULL : options.handle

    stream_handle = mongoc_gridfs_bucket_open_upload_stream_with_id(
            bucket.handle, file_id.handle, filename, options_handle,
            err_ref)

    if stream_handle == C_NULL
        throw(err_ref[])
    end

    return MongoIOStream(
            bucket,
            stream_handle,
            timeout_msec=timeout_msec,
            chunk_size=chunk_size
        )
end

function open_upload_stream(bucket::Bucket, file_info::BSON;
        options::Union{Nothing, BSON}=nothing,
        timeout_msec::Integer=DEFAULT_TIMEOUT_MSEC,
        chunk_size::Integer=DEFAULT_CHUNK_SIZE) :: MongoIOStream

    @assert haskey(file_info, "filename") "`file_info` must have a `filename` field."

    # generates a new _id if not present in `file_info`
    if !haskey(file_info, "_id")
        file_info["_id"] = BSONObjectId()
    end

    return open_upload_stream(
            bucket, get_as_bson_value(file_info, "_id"), file_info["filename"],
            options=options, timeout_msec=timeout_msec, chunk_size=chunk_size
        )
end

"""
    open_upload_stream(bucket, file_id, filename; [options], [timeout_msec], [chunk_size]) :: MongoIOStream
    open_upload_stream(bucket, filename; [options], [timeout_msec], [chunk_size]) :: MongoIOStream
    open_upload_stream(bucket, file_info; [options], [timeout_msec], [chunk_size]) :: MongoIOStream

Opens a stream to upload a file to a GridFS Bucket.

`file_info` is a BSON document with the following fields:

* `_id` as an optional identifier.

* `filename` as the name of the file in the remote bucket.

If `_id` is not provided, a `BSONObjectId` will be generated.

# Example

```julia
data = rand(UInt8, 3_000_000)
remote_filename = "uploaded.data"
io = Mongoc.open_upload_stream(bucket, remote_filename)
write(io, data)
close(io)
```
"""
function open_upload_stream(bucket::Bucket, file_id, filename::AbstractString;
        options::Union{Nothing, BSON}=nothing,
        timeout_msec::Integer=DEFAULT_TIMEOUT_MSEC,
        chunk_size::Integer=DEFAULT_CHUNK_SIZE) :: MongoIOStream

    file_info = BSON("_id" => file_id, "filename" => filename)
    return open_upload_stream(bucket, file_info,
            options=options, timeout_msec=timeout_msec, chunk_size=chunk_size)
end

function open_upload_stream(bucket::Bucket, filename::AbstractString;
        options::Union{Nothing, BSON}=nothing,
        timeout_msec::Integer=DEFAULT_TIMEOUT_MSEC,
        chunk_size::Integer=DEFAULT_CHUNK_SIZE) :: MongoIOStream

    file_info = BSON("filename" => filename)
    return open_upload_stream(bucket, file_info,
            options=options, timeout_msec=timeout_msec, chunk_size=chunk_size)
end


text:
function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    solution_info = model.inner.solution_info
    if solution_info === nothing
        return MOI.OPTIMIZE_NOT_CALLED
    end
    solver_status = solution_info.solver_status
    model_status = solution_info.model_status

    if solver_status == NORMAL_COMPLETION
        if model_status == OPTIMAL
            return MOI.OPTIMAL
        elseif model_status == INFEASIBLE
            return MOI.INFEASIBLE
        elseif model_status == UNBOUNDED
            return MOI.DUAL_INFEASIBLE
        elseif model_status == INTERMEDIATE_FEASIBLE
            return LOCALLY_SOLVED
        elseif model_status == UNKNOWN
            return MOI.OTHER_ERROR
        end
    elseif solver_status == INSUFFICIENT_MEMORY_FOR_NODES
        return MOI.MEMORY_LIMIT
    elseif solver_status == ITERATION_LIMIT
        return MOI.ITERATION_LIMIT
    elseif solver_status == TIME_LIMIT
        return MOI.TIME_LIMIT
    elseif solver_status == NUMERICAL_SENSITIVITY
        return MOI.NUMERICAL_ERROR
    elseif solver_status == INSUFFICIENT_MEMORY_FOR_SETUP
        return MOI.MEMORY_LIMIT
    elseif solver_status == RESERVED
        return MOI.OTHER_ERROR
    elseif solver_status == TERMINATED_BY_BARON
        return MOI.OTHER_ERROR
    elseif solver_status == SYNTAX_ERROR
        return MOI.INVALID_MODEL
    elseif solver_status == LICENSING_ERROR
        return MOI.OTHER_ERROR
    elseif solver_status == USER_HEURISTIC_TERMINATION
        return MOI.OTHER_LIMIT
    end

    error("Unrecognized Baron status: $solver_status, $model_status")
end

function MOI.get(model::Optimizer, ::MOI.ResultCount)
    return (model.inner.solution_info.feasible_point === nothing) ? 0 : 1
end

function MOI.get(model::Optimizer, ::MOI.PrimalStatus)
    solution_info = model.inner.solution_info
    if solution_info === nothing || solution_info.feasible_point === nothing
        return MOI.NO_SOLUTION
    else
        return solution_info.model_status == UNBOUNDED ? MOI.INFEASIBILITY_CERTIFICATE : MOI.FEASIBLE_POINT
    end
end

MOI.get(model::Optimizer, ::MOI.ObjectiveValue) = model.inner.solution_info.objective_value

function MOI.get(model::Optimizer, ::MOI.VariablePrimal, vi::VI)
    solution_info = model.inner.solution_info
    if solution_info === nothing || solution_info.feasible_point === nothing
        error("VariablePrimal not available.")
    end
    check_variable_indices(model, vi)
    return solution_info.feasible_point[vi.value]
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveBound)
    solution_info = model.inner.solution_info
    return solution_info.dual_bound
end

function MOI.get(model::Optimizer, ::MOI.SolveTime)
    solution_info = model.inner.solution_info
    return solution_info.wall_time
end

# TODO: desirable?
function MOI.get(model::MOIU.CachingOptimizer{BARON.Optimizer}, attr::MOI.ConstraintPrimal, ci::MOI.ConstraintIndex)
    return MOIU.get_fallback(model, attr, ci)
end

# TODO: MOI getter for solvetime


text:
using Weave

weave("01_WhyJulia.jmd", doctype=:"github", out_path=:pwd)
weave("02_Functions.jmd", doctype=:"github", out_path=:pwd)
weave("03_Graphs.jmd", doctype=:"github", out_path=:pwd)
weave("04_Transformations.jmd", doctype=:"github", out_path=:pwd)
weave("05_Inverse.jmd", doctype=:"github", out_path=:pwd)
weave("06_Polynomials.jmd", doctype=:"github", out_path=:pwd)
weave("07_PolynomialRoots.jmd", doctype=:"github", out_path=:pwd)
weave("08_Rational-Functions.jmd", doctype=:"github", out_path=:pwd)
weave("09_Exp-Log.jmd", doctype=:"github", out_path=:pwd)
weave("10_Trigonometric-Functions.jmd", doctype=:"github", out_path=:pwd)
weave("11_Limits.jmd", doctype=:"github", out_path=:pwd)
weave("12_Limit-Extensions.jmd", doctype=:"github", out_path=:pwd)
weave("13_Continuity.jmd", doctype=:"github", out_path=:pwd)
weave("13_Continuity.jmd", doctype=:"github", out_path=:pwd)

max_stars_repo_path: Deep Into to Julia/build.jl

max_stars_repo_name: bmoretz/Statistical-Computing

max_issues_repo_path: Deep Into to Julia/build.jl

max_issues_repo_name: bmoretz/Statistical-Computing

max_forks_repo_path: Deep Into to Julia/build.jl

max_forks_repo_name: bmoretz/Statistical-Computing


text:
using Revise

using Connectomes
using Test
using SparseArrays
using DataFrames
using LinearAlgebra
using FileIO

@testset "Connectomes.jl" begin
    connectome_path = "/"*relpath((@__FILE__)*"/../..","/") * "/assets/connectomes/hcp-scale1-standard-master.graphml"

    connectome = Connectome(connectome_path)

    @test typeof(connectome) == Connectome
    @test typeof(connectome.A) == SparseMatrixCSC{Float64, Int64}

    @test typeof(connectome.parc) == DataFrame
    @test typeof(connectome.parc.x) == typeof(connectome.parc.y) == typeof(connectome.parc.z)

    @test length(connectome.parc.ID) == 83

end


text:
import Base: iterate, indexed_iterate

### ITERATION: We define iterate in terms of getindex which does the necessary tracking.

# Iterate is used in `for x in a` loops
iterate(a::Value,i=nothing)=throw(MethodError(iterate,(a,i)))
iterate(a::Value{T},i=1) where {T<:AbstractArray} = i > length(a) ? nothing : (a[i],i+1)
iterate(a::Value{T},i=1) where {T<:Tuple}         = i > length(a) ? nothing : (a[i],i+1)
iterate(a::Value{T})   where {T<:AbstractDict} = (v=iterate(value(a));   v==nothing ? v : (((k,_),j)=v;(k=>a[k],j)))
iterate(a::Value{T},i) where {T<:AbstractDict} = (v=iterate(value(a),i); v==nothing ? v : (((k,_),j)=v;(k=>a[k],j)))
iterate(a::Value{T})   where {T<:Number} = (a,nothing)
iterate(a::Value{T},i) where {T<:Number} = nothing

# indexed_iterate for `(x,y)=a` multiple assignments.
indexed_iterate(a::Value{T},i::Int,state=1) where {T<:Tuple} = (a[i],i+1)
indexed_iterate(a::Value{T},i::Int,state=1) where {T<:Array} = (a[i],i+1)


text:
using LaserLab
using Test
using Unitful
import Unitful:
    nm, μm, mm, cm, m, km,
    mg, g, kg,
    ps, ns, μs, ms, s, minute, hr, d, yr, Hz, kHz, MHz, GHz,
    eV,
    μJ, mJ, J,
	μW, mW, W

function hnu_eV(lambda_nm)
    return 1240.0 * (1.0/lambda_nm)
end

dspot(λ::Real, NA::Real) = 1.83*λ/(2*NA)

include("dffunctions_test.jl")
include("setup_test.jl")
include("glaser_test.jl")

max_stars_repo_name: proscrite/LaserLab

max_issues_repo_name: proscrite/LaserLab

max_forks_repo_name: proscrite/LaserLab


text:
using MatpowerCases, Test

println("- extract_bus check")
for c in casenames()
    mpc = loadcase(c; describe=false)
    @assert isa(extract_bus(mpc), MatpowerCases.Bus)
    @assert isa(extract_bus(c), MatpowerCases.Bus)
end


text:
using Base.Test
using Clustering

srand(34568)

m = 3
n = 1000
k = 5

x = rand(m,n)

fuzziness = 2.0
srand(34568)
r = fuzzy_cmeans(x, k, fuzziness)
@test isa(r, FuzzyCMeansResult{Float64})
@test size(r.centers) == (m,k)
@test size(r.weights) == (n,k)
@test all([s ≈ 1.0 for s in sum(r.weights,2)])
@test all(0 .<= r.weights .<= 1)

fuzziness = 3.0
srand(34568)
r = fuzzy_cmeans(x, k, fuzziness)
@test isa(r, FuzzyCMeansResult{Float64})
@test size(r.centers) == (m,k)
@test size(r.weights) == (n,k)
@test all([s ≈ 1.0 for s in sum(r.weights,2)])
@test all(0 .<= r.weights .<= 1)


text:
push!(LOAD_PATH,"../src/")

using Documenter, EarthEngine

pages = [
    "Home" => "index.md",
    "Usage" => "usage.md",
    "Examples" => [
        "Landsat 8 Cloud Masking" => "examples/landsat8_cloudmasking.md",
        "Sentinel 1 Composite" => "examples/sentinel1_composite.md",
        "Trendy Lights" => "examples/trendy_lights.md",
        "Classification" => "examples/classification.md",
    ],
    "API" => "api.md",
]

makedocs(;
    modules = [EarthEngine],
    authors = "Kel Markert",
    repo = "https://github.com/KMarkert/EarthEngine.jl/blob/{commit}{path}#L{line}",
    sitename = "EarthEngine.jl",
    # format = Documenter.HTML(;
    #     prettyurls = get(ENV, "CI", "false") == "true",
    #     canonical = "https://deltares.github.io/Wflow.jl",
    #     assets = String[],
    # ),
    pages = pages,
)

deploydocs(;
    deps = Deps.pip("earthengine-api"),
    repo = "github.com/KMarkert/EarthEngine.jl.git",
    devbranch = "main"
)


text:
using LazyArrays: ApplyArray

"""
    parent_coefficients(R::RecombinedBSplineBasis, coefs::AbstractVector)

Returns the coefficients associated to the parent B-spline basis, from the
coefficients `coefs` in the recombined basis.

Note that this function doesn't allocate, since it returns a lazy concatenation
(via LazyArrays.jl) of two StaticArrays and a view of the `coefs` vector.
"""
function parent_coefficients(R::RecombinedBSplineBasis, coefs::AbstractVector)
    length(coefs) == length(R) ||
        throw(DimensionMismatch("wrong number of coefficients"))

    M = recombination_matrix(R)

    # Upper-left and lower-right submatrices
    A = M.ul
    C = M.lr

    N = length(coefs)
    Na = size(A, 2)
    Nc = size(C, 2)
    H = N - Nc

    Base.require_one_based_indexing(coefs)

    xA = @inbounds view(coefs, 1:Na)
    xC = @inbounds view(coefs, (H + 1):N)

    yA = A * xA
    yB = @inbounds view(coefs, (Na + 1):H)
    yC = C * xC

    ApplyArray(vcat, yA, yB, yC)
end

# Returns spline written in the parent basis (usually a regular BSplineBasis).
# Those splines can be evaluated, differentiated, etc... using known algorithms.
function Splines.parent_spline(R::RecombinedBSplineBasis, S::Spline)
    @assert basis(S) === R
    B = parent(R)
    coefs = parent_coefficients(R, coefficients(S))
    Spline(B, coefs)
end

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
"""
    CategoricalTree{T}
    where T <: Real

A Datastructure designed for efficiently sampling and mutating a categorical distribution with many categories.

Given an array of values, we store these, and their sums in a tree where the leaf nodes are the values, and each parent node contains the sum of up to two children, all the way up to the root node which contains the sum of all values.

This allows for ``O(1)`` read, ``O(\\log n)`` write, ``O(1)`` sum, and ``O(\\log n)`` sampling of a random index (when the values are describe a categorical distribution).
"""
struct CategoricalTree{T <: Real}
    ns::Vector{Int} # number of nodes at each layer of the tree (`ns[end]` should be 1, i.e. the root of the tree
    as::Vector{Vector{T}} # store each layer into a separate array. Each parent node contains the sum of its children
end

function CategoricalTree(a::Vector{T}) where T<:Real
    n = length(a)
    ns::Vector{Int64} = [n]
    while n > 1
        if isodd(n)
            n = Int((n+1)/2)
            push!(ns,n)
        else
            n = Int(n/2)
            push!(ns,n)
        end
    end
    as = [zeros(T,k) for k in ns]
    as[1] = a
    for i in 1:(length(ns)-1)
        for j in 2:2:ns[i]
            as[i+1][Int(j/2)] = as[i][j-1] + as[i][j]
        end
        if isodd(ns[i])
            as[i+1][ns[i+1]] = as[i][ns[i]]
        end
    end
    return CategoricalTree{T}(ns,as)
end

Base.iterate(ct::CategoricalTree) = iterate(ct.as[1])
Base.iterate(ct::CategoricalTree, state) = iterate(ct.as[1], state)

Base.isempty(ct::CategoricalTree) = isempty(ct.as[1])
Base.length(ct::CategoricalTree) = length(ct.as[1])

Base.firstindex(::CategoricalTree) = 1
Base.lastindex(ct::CategoricalTree) = length(ct)

Base.in(ct::CategoricalTree) = in(ct.as[1])

Base.getindex(ct::CategoricalTree, i) = ct.as[1][i]

function Base.setindex!(ct::CategoricalTree{T}, x::T, i) where {T}
    ct.as[1][i] = x
    j = i
    for k in 2:length(ct.ns)
        j = isodd(j) ? Int((j+1)/2) : Int(j/2)
        if isodd(ct.ns[k-1]) && j == ct.ns[k]
                ct.as[k][j] = ct.as[k-1][2*j-1]
        else
            ct.as[k][j] = ct.as[k-1][2*j-1] + ct.as[k-1][2*j]
        end
    end
end

# since we already computed the sum, we might as well make use of it
Base.sum(ct::CategoricalTree) = ct.as[end][1]

"""
    rand_categorical(ct::CategoricalTree{<:Real}, a0)

Pick one index according to a categorical distribution described by `ct`.

The leaf values `a[k]` are assumed to all be non-negative. The index `k` is picked with probability `a[k]/sum(ct)`.
"""
function rand_categorical(ct::CategoricalTree, a0 = sum(ct))
    if length(ct.ns) == 1
        return 1
    end
    r = rand()*sum(ct) # ignore a0
    i = length(ct.ns) - 1
    j = 1
    while i > 1
        if r <= ct.as[i][j]
            j = 2*j - 1
        else
            r -= ct.as[i][j]
            j = 2*j + 1
        end
        i -= 1
    end
    if r <= ct.as[i][j]
        return j
    else
        return j+1
    end
end


text:
package java.util;

public class ArrayList implements Collection {
    private Object[] objects;
    private int size;

    public ArrayList() {
	objects = new Object[1];
	size = 0;
    }

    public int indexOf(Object o) {
	for (int i = 0; i < objects.length; i++) {
	    if (objects[i].equals(o)) {
		return i;
	    }
	}
	return -1;
    }

    public int size() {
	return size;
    }

    public Object get(int i) {
	return objects[i];
    }

    public Object remove(int i) {
	Object o = objects[i];
	objects[i] = null;
	size -= 1;
	return o;
    }

    public boolean remove(Object o) {
	if (o.equals(objects[0])) {
	    remove(0);
	    return true;
	} else {
	    return false;
	}
    }

    public boolean add(Object o) {
	boolean present = o.equals(objects[0]);
	objects[0] = o;
	size += 1;
	return present;
    }

    public boolean isEmpty() {
	return size > 0;
    }

    public boolean contains(Object o) {
	return o.equals(objects[0]);
    }

    public boolean containsAll(Collection c) {
	return contains(c);
    }

    public boolean addAll(Collection c) {
	return add(c);
    }

    public boolean retainAll(Collection c) {
	return remove(c);
    }

    public boolean removeAll(Collection c) {
	return remove(c);
    }

    public void clear() {
	removeAll(this);
    }

    private class ALIterator implements Iterator {
	private ArrayList al;

	public ALIterator(ArrayList al) {
	    this.al = al;
	}

	public boolean hasNext() {
	    return al.isEmpty();
	}

	public Object next() {
	    return al.get(0);
	}

	public void remove() {
	    al.remove(0);
	}
    }

    public Object[] toArray(Object[] array) {
	for (int i = 0; i < size; i++) {
	    array[i] = objects[i];
	}
	return array;
    }

    public Object[] toArray() {
	return objects;
    }

    public Iterator iterator() {
	return new ALIterator(this);
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
using SeisRequests, Test, Dates

@testset "Common" begin
    @testset "POST" begin
        let r1 = FDSNStation()
    end
end


text:
########################################################################
# Gradient oracles
########################################################################
using Printf

# run SVRG to obtain an exact/ineact dual gradient
function svrg(X::SparseMatrixCSC{Float64, Int64}, Y::Vector{Int64}, w::Matrix{Float64}, y::Matrix{Float64})
    1
end

# run Newton's method to obtain an exact/ineact dual gradient
function newton!(X::SparseMatrixCSC{Float64, Int64}, Xt::SparseMatrixCSC{Float64, Int64}, Y::Vector{Int64}, λ::Float64, W::Matrix{Float64}, y::Matrix{Float64}, η::Float64; T::Int64 = 50, tol::Float64=1e-6)
    for t = 1:T
        g = getGradient(X, Xt, Y, W, λ) - y
        gnorm = norm(g)
        # @printf("   gnorm: %4.4e\n", gnorm)
        if gnorm < tol
            break
        end
        t1 = time()
        D = ComputeNewtonDirection( X, Xt, W, λ, g)
        t2 = time()
        # @printf("Time spent on computing Newton direction: %4.4f\n", t2 - t1)
        η = lineSearch2(X, Y, y, D, W, g, λ)
        # η = 1.0
        W .-= η*D
        if t == T
            @warn("Reached maximum iteration in Newton's method.")
        end
    end
end

# run SGD to obtain an exact/ineact dual gradient
function sgd!(X::SparseMatrixCSC{Float64, Int64}, Xt::SparseMatrixCSC{Float64, Int64}, Y::Vector{Int64}, λ::Float64, W::Matrix{Float64}, y::Matrix{Float64}, η::Float64; T::Int64 = 20, is_lazy::Bool = true)
    n, d = size(X)
    _, k = size(W)
    if is_lazy
        hitTime = zeros(Int64, d, k)
    end
    perm = collect(1:n)
    for t = 1:T
        shuffle!(perm)
        if is_lazy
            fill!(hitTime, 0)
            timeStep = 1
        end
        for i in perm
            if is_lazy
                g = getStochasticGrad(Xt, Y, W, i)
                I, J, V = findnz(g)
                for j = 1:length(I)
                    idx1 = I[j]
                    idx2 = J[j]
                    # Lazy update
                    W[idx1, idx2] *= (1 - η*λ)^(timeStep-hitTime[idx1, idx2])
                    W[idx1, idx2] += (timeStep-hitTime[idx1, idx2])*η*y[idx1, idx2]
                    W[idx1, idx2] -= η*V[j]
                    # Update hitTime to the current time
                    hitTime[idx1, idx2] = timeStep
                end
                timeStep += 1
            else
                g = getStochasticGrad(Xt, Y, W, i) + λ*W - y
                W .-= η*g
            end
        end
        # Lazy update for staled coordinates
        if is_lazy
            timeStep -= 1
            for i = 1:d
                for j = 1:k
                    if hitTime[i, j] < timeStep
                        W[i, j] *= (1 - η*λ)^(timeStep-hitTime[i, j])
                        W[i, j] += (timeStep-hitTime[i, j])*η*y[i, j]
                    end
                end
            end
        end
    end
    return nothing
end

# run adam to obtain an exact/ineact dual gradient
function adam!(Xt::SparseMatrixCSC{Float64, Int64}, Y::Flux.OneHotArray, W::Flux.Chain, y::Zygote.Params, λ::Float64; num_epoches::Int64=5)
    # data
    data = Flux.Data.DataLoader((Xt, Y), batchsize=128, shuffle=true)
    # loss
    sqnorm(w) = sum(abs2, w)
    loss(x, l) = Flux.crossentropy(W(x), l) + (λ/2) * sum(sqnorm, params(W)) - dot(params(W), y)
    # optimizer
    # opt = ADAMW(0.001, (0.89, 0.995), 1e-8)
    opt = Descent(1e-1)
    # train
    for t = 1:num_epoches
        Flux.train!(loss, params(W), data, opt)
        # @printf "epoch: %d, obj: %.2f, acc: %.2f\n" t obj(Xt, Y, W, λ) accuracy(Xt', Y, W)
    end
end

# Run SGD for each individual problem with neural nets as their models.
function sgd!(Xt::SparseMatrixCSC{Float64, Int64}, Y::Flux.OneHotArray, W::Flux.Chain, λ::Float64, μ::Float64, lr::Float64, num_epoches::Int64)
    # data
    data = Flux.Data.DataLoader((Xt, Y), batchsize=128, shuffle=true)
    # Store Wold
    Wold = deepcopy( params(W) )
    # loss
    sqnorm(w) = sum(abs2, w)
    loss(x, l) = Flux.crossentropy(W(x), l) + ((λ+μ)/2) * sum(sqnorm, params(W)) - μ*dot(params(W), Wold)
    # optimizer
    opt = Descent(lr)
    # train
    for t = 1:num_epoches
        # @printf "   epoch: %d, obj: %4.4e, acc: %4.4e\n" t obj(Xt, Y, W, λ) accuracy(Xt', Y, W)
        Flux.train!(loss, params(W), data, opt)
    end
end


text:
# Autogenerated wrapper script for S2Geometry_jll for powerpc64le-linux-gnu-cxx11
export libs2

using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("S2Geometry")
JLLWrappers.@declare_library_product(libs2, "libs2.so.0")
function __init__()
    JLLWrappers.@generate_init_header(OpenSSL_jll)
    JLLWrappers.@init_library_product(
        libs2,
        "lib/libs2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using PtFEM, Test

data = Dict(
  # Plane(ndim, nst, nxe, nye, nip, direction, finite_element(nod, nodof), axisymmetric)
  :struc_el => Plane(2, 3, 2, 2, 1, :x, Triangle(3, 2), false),
  :properties => [1.0e6 0.3;],
  :x_coords => [0.0,  0.5,  1.0],
  :y_coords => [0.0,  -0.5,  -1.0],
  :support => [
    (1, [0 1]),
    (4, [0 1]),
    (7, [0 0]),
    (8, [1 0]),
    (9, [1 0])
    ],
  :loaded_nodes => [
    (1, [0.0 -0.25]),
    (2, [0.0 -0.50]),
    (3, [0.0 -0.25])
    ]
)

@time m, dis_df, fm_df = p51(data)

arr = Vector{Float64}(m.loads[0:12])
@test  arr ≈ [0.0, -9.100000000000005e-7, 1.950000000000001e-7,
-9.100000000000006e-7, 3.900000000000002e-7, -9.1e-7, -4.5500000000000025e-7,
1.950000000000002e-7, -4.550000000000004e-7, 3.900000000000004e-7,
-4.5499999999999993e-7, 1.9500000000000038e-7, 3.9000000000000045e-7]


text:
using DataFrames
using Printf
using PrettyTables

struct TrainingBatch
    game_type::Type{<:AbstractEnv}
    action_space_size::Integer
    samples::Vector{TrainingSample}
end

samples(x::TrainingBatch) = x.samples
# @todo change to get_tensor
tensor(x::TrainingBatch) = TrainingBatchTensor(x.game_type, x.action_space_size, x.samples)

function Base.show(io::IO, m::MIME"text/plain", x::TrainingBatch)
    action_hl = Highlighter((data, i, j)->(j==3) && data[i,j] == nothing , crayon"#666666");
    pretty_table(
        io,
        get_df(x);
        show_row_number=true,
        nosubheader=true,
        title="Training Batch",
        highlighters=(action_hl),
    )
end

Base.getindex(x::TrainingBatch, idx::Integer) = x.samples[idx]
Base.length(x::TrainingBatch) = length(x.samples)
Base.lastindex(x::TrainingBatch) = length(x)
Base.iterate(x::TrainingBatch, state=1) = state > length(x) ? nothing : (x[state], state+1)

function Base.rand(x::TrainingBatch)
    if isempty(x.samples)
        return nothing
    end
    rand(x.samples)
end

function get_state(x::TrainingBatch, row::Integer)
    state_type = x.samples[1].image_state |> typeof
    df = x |> get_df
    state_key = df."State Key"[row]
    state_type(state_key)
end

function get_df(x::TrainingBatch)
    sampleidx_col = "Sample #"
    df = reduce(
        append!,
        map(
            (x)->begin
                idx, sample = x
                df = getdf(sample)
                df[!, sampleidx_col] = [idx for _ in 1:size(df, 1)]
                df
            end,
            enumerate(x.samples)
        )
    )
    ncols = size(df,2)
    df[!, [ncols, vec(1:(ncols-1))...]]
end


text:
# #-------------------------------------------------------------------------------------------------------
# #-------------------------------------------------------------------------------------------------------
# # This file contains functions used to calculate the log-posterior.
# #-------------------------------------------------------------------------------------------------------
# #-------------------------------------------------------------------------------------------------------

#######################################################################################################################
# CALCULATE TOTAL (LOG) PRIOR PROBABILITY.
#######################################################################################################################
# Description: This creates a function that will calculates the total (log) prior probability of the uncertain model,
#              initial condition, and statistical process parameters specific to the CH₄ cycle model selected.
#
# Function Arguments:
#
#       climate_model = A symbol identifying the specific version of SNEASY+CH4 (options are :sneasy_fair, :sneasy_fund,
#                       :sneasy_hector, and :sneasy_magicc).
#----------------------------------------------------------------------------------------------------------------------

function construct_log_prior(climate_model::Symbol)

    # Declare all prior distributions for all uncertain parameters found in the four versions of SNEASY+CH4.

    # -----------------------------------------
    # Statistical Process Priors.
    # -----------------------------------------
    prior_σ_temperature      = Uniform(0, 0.2)
    prior_σ_ocean_heat       = Uniform(0, 4)
    prior_σ²_white_noise_CO₂ = Uniform(0, 200)
    prior_σ²_white_noise_CH₄ = Uniform(0, 800)

    prior_ρ_temperature      = Uniform(-0.99, 0.99)
    prior_ρ_ocean_heat       = Uniform(-0.99, 0.99)
    prior_α₀_CO₂             = Uniform(0.01, 11.5)
    prior_α₀_CH₄             = Uniform(0.01, 11.5)

    # -----------------------------------------
    # Initial Condition Priors.
    # -----------------------------------------
    prior_temperature_0     = Normal()
    prior_ocean_heat_0      = Uniform(-100, 0)
    prior_CO₂_0             = Uniform(275, 281)
    prior_CH₄_0             = Uniform(691, 753)
    prior_N₂O_0             = Uniform(264, 282)

    # -----------------------------------------
    # Climate & Radiative Forcing Priors.
    # -----------------------------------------
    prior_ECS               = Truncated(Cauchy(3.0,2.0), 0.0, 10.0)
    prior_heat_diffusivity  = LogNormal(1.1, 0.3)
    prior_rf_scale_aerosol  = TriangularDist(0., 3., 1.)
    prior_rf_scale_CH₄      = Uniform(0.72, 1.28)
    prior_F2x_CO₂           = Uniform(2.968, 4.452)

    # -----------------------------------------
    # Carbon Cycle Priors.
    # -----------------------------------------
    prior_Q10               = Uniform(1.0, 5)
    prior_CO₂_fertilization = Uniform(0., 1)
    prior_CO₂_diffusivity   = Uniform(0., 200)

    # -----------------------------------------
    # Methane Cycle Priors.
    # -----------------------------------------
    prior_τ_troposphere     = Uniform(6.7, 13.5)
    prior_CH₄_natural       = Uniform(150, 484)
    prior_τ_soil            = Uniform(96, 144)
    prior_τ_stratosphere    = Uniform(100, 200)

    #----------------------------------------------------------------------------------------------------------------
    # Create a function unique to each CH₄ cycle model to calculate prior distribution for uncertain CH₄ parameters.
    #----------------------------------------------------------------------------------------------------------------

    # CH₄ prior function (requires the particular ordering of uncertain parameters used below).
    ch4_model_priors =

        # Create function based on uncertain SNEASY-FAIR CH₄ cycle parameters.
        if climate_model == :sneasy_fair
            function(p)
                τ_troposphere = p[22]
                CH₄_natural   = p[23]
                return logpdf(prior_τ_troposphere, τ_troposphere) + logpdf(prior_CH₄_natural, CH₄_natural)
            end

        # Create function based on uncertain SNEASY-FUND CH₄ cycle parameters.
        elseif climate_model == :sneasy_fund
            function(p)
                τ_troposphere = p[22]
                return logpdf(prior_τ_troposphere, τ_troposphere)
            end

        # Create function based on uncertain SNEASY-Hector or SNEASY-MAGICC CH₄ cycle parameters.
        elseif climate_model == :sneasy_hector || climate_model == :sneasy_magicc
            function(p)
                τ_troposphere  = p[22]
                CH₄_natural    = p[23]
                τ_soil         = p[24]
                τ_stratosphere = p[25]
                return logpdf(prior_τ_troposphere, τ_troposphere) + logpdf(prior_CH₄_natural, CH₄_natural) + logpdf(prior_τ_soil, τ_soil) + logpdf(prior_τ_stratosphere, τ_stratosphere)
            end

        # Throw an error for incorrect model selection.
        else
            error("Incorrect climate model selected. Options include :sneasy_fair, :sneasy_fund, :sneasy_hector, or :sneasy_magicc.")
        end

    #------------------------------------------------------------------------------------
    # Create function that returns the log-prior of all uncertain model parameters.
    #------------------------------------------------------------------------------------

    function total_log_prior(p::Array{Float64,1})

        # Assign parameter values (common to all four models) names for convenience.
        σ_temperature      = p[1]
        σ_ocean_heat       = p[2]
        σ²_white_noise_CO₂ = p[3]
        σ²_white_noise_CH₄ = p[4]
        ρ_temperature      = p[5]
        ρ_ocean_heat       = p[6]
        α₀_CO₂             = p[7]
        α₀_CH₄             = p[8]
        temperature_0      = p[9]
        ocean_heat_0       = p[10]
        CO₂_0              = p[11]
        CH₄_0              = p[12]
        N₂O_0              = p[13]
        ECS                = p[14]
        heat_diffusivity   = p[15]
        rf_scale_aerosol   = p[16]
        rf_scale_CH₄       = p[17]
        F2x_CO₂            = p[18]
        Q10                = p[19]
        CO₂_fertilization  = p[20]
        CO₂_diffusivity    = p[21]

        # Return total log-prior of all uncertain parameters (using "ch4_model_priors" function to account for different versions of SNEASY+CH4).
        log_prior = logpdf(prior_σ_temperature, σ_temperature) + logpdf(prior_σ_ocean_heat, σ_ocean_heat) + logpdf(prior_σ²_white_noise_CO₂, σ²_white_noise_CO₂) + logpdf(prior_σ²_white_noise_CH₄, σ²_white_noise_CH₄) +
                    logpdf(prior_ρ_temperature, ρ_temperature) + logpdf(prior_ρ_ocean_heat, ρ_ocean_heat) + logpdf(prior_α₀_CO₂, α₀_CO₂) + logpdf(prior_α₀_CH₄, α₀_CH₄) +
                    logpdf(prior_temperature_0, temperature_0) + logpdf(prior_ocean_heat_0, ocean_heat_0) + logpdf(prior_CO₂_0, CO₂_0) + logpdf(prior_CH₄_0, CH₄_0) + logpdf(prior_N₂O_0, N₂O_0) +
                    logpdf(prior_ECS, ECS) + logpdf(prior_heat_diffusivity, heat_diffusivity) + logpdf(prior_rf_scale_aerosol, rf_scale_aerosol) + logpdf(prior_rf_scale_CH₄, rf_scale_CH₄) + logpdf(prior_F2x_CO₂, F2x_CO₂) +
                    logpdf(prior_Q10, Q10) + logpdf(prior_CO₂_fertilization, CO₂_fertilization) + logpdf(prior_CO₂_diffusivity, CO₂_diffusivity) +
                    ch4_model_priors(p)

        return log_prior
    end

    # Return total log-prior function for specified version of SNEASY+CH4.
    return total_log_prior
end

#######################################################################################################################
# CALCULATE LOG POSTERIOR.
#######################################################################################################################
# Description: This creates a function that calculates the log-posterior probability of the uncertain model, initial
#              condition, and statistical process parameters specific to the CH₄ cycle model selected.
#
# Function Arguments:
#
#       f_run_model   = A function that runs the specific version of SNEASY+CH4 and returns the output being calibrated to observations.
#       climate_model = A symbol identifying the specific version of SNEASY+CH4 (options are :sneasy_fair, :sneasy_fund,
#                       :sneasy_hector, and :sneasy_magicc).
#       end_year      = The final year to run the model calibration (defaults to 2017).
#----------------------------------------------------------------------------------------------------------------------

function construct_log_posterior(f_run_model, climate_model::Symbol; end_year::Int=2017)

    # Initialize model to start in 1765.
    start_year = 1765

    # Get log-prior function based on version of SNEASY+CH4 being used.
    total_log_prior = construct_log_prior(climate_model)

    # Create a vector of calibration years and calculate total number of years to run model.
    calibration_years = collect(start_year:end_year)
    n = length(calibration_years)

    # Load calibration data/observations.
    calibration_data = load_calibration_data(end_year)

    # Calculate indices for each year that has an observation in calibration data sets.
    indices_temperature_data   = findall(x-> !ismissing(x), calibration_data.hadcrut_temperature_obs)
    indices_oceanheat_data     = findall(x-> !ismissing(x), calibration_data.ocean_heat_obs)
    indices_oceanco2_flux_data = findall(x-> !ismissing(x), calibration_data.oceanco2_flux_obs)
    indices_maunaloa_co2_data  = findall(x-> !ismissing(x), calibration_data.maunaloa_co2_obs)
    indices_lawdome_co2_data   = findall(x-> !ismissing(x), calibration_data.lawdome_co2_obs)
    indices_noaa_ch4_data      = findall(x-> !ismissing(x), calibration_data.noaa_ch4_obs)
    indices_lawdome_ch4_data   = findall(x-> !ismissing(x), calibration_data.lawdome_ch4_obs)

    # Combine CO₂ and CH₄ indices from Law Dome and more recent observations.
    indices_co2_data = sort(vcat(indices_lawdome_co2_data, indices_maunaloa_co2_data))
    indices_ch4_data = sort(vcat(indices_noaa_ch4_data, indices_lawdome_ch4_data))

    # Combine CO₂ and CH₄ measurement errors from Law Dome and more recent observations (just for convenience).
    calibration_data.co2_combined_sigma = calibration_data.lawdome_co2_sigma
    calibration_data.co2_combined_sigma[indices_maunaloa_co2_data] = calibration_data.maunaloa_co2_sigma[indices_maunaloa_co2_data]

    calibration_data.ch4_combined_sigma = calibration_data.lawdome_ch4_sigma
    calibration_data.ch4_combined_sigma[indices_noaa_ch4_data] = calibration_data.noaa_ch4_sigma[indices_noaa_ch4_data]

    # Calculate number of ice core observations for CO₂ and CH₄ (used for indexing).
    n_lawdome_co2 = length(indices_lawdome_co2_data)
    n_lawdome_ch4 = length(indices_lawdome_ch4_data)

    # Allocate arrays to store data-model residuals.
    temperature_residual = zeros(length(indices_temperature_data))
    ocean_heat_residual  = zeros(length(indices_oceanheat_data))
    co2_residual         = zeros(length(indices_co2_data))
    ch4_residual         = zeros(length(indices_ch4_data))

    # Allocate array to store likelihoods for individual ocean CO₂ flux data points that need to be summed up to get a total likelihood (assuming iid error structure).
    oceanco2_flux_single_llik   = zeros(length(indices_oceanco2_flux_data))

    # Allocate vectors to store model output being calibrated to the observations.
    modeled_CO₂           = zeros(n)
    modeled_CH₄           = zeros(n)
    modeled_oceanCO₂_flux = zeros(n)
    modeled_temperature   = zeros(n)
    modeled_ocean_heat    = zeros(n)

    #---------------------------------------------------------------------------------------------------------------------------------------
    # Create a function to calculate the log-likelihood for the observations, assuming residual independence across calibration data sets.
    #---------------------------------------------------------------------------------------------------------------------------------------

    function total_log_likelihood(p::Array{Float64,1})

        # Assign names to uncertain statistical process parameters used in log-likelihood calculations.
        σ_temperature      = p[1]
        σ_ocean_heat       = p[2]
        σ²_white_noise_CO₂ = p[3]
        σ²_white_noise_CH₄ = p[4]
        ρ_temperature      = p[5]
        ρ_ocean_heat       = p[6]
        α₀_CO₂             = p[7]
        α₀_CH₄             = p[8]

        # Run an instance of SNEASY+CH4 with sampled parameter set and return model output being compared to observations.
        f_run_model(p, modeled_CO₂, modeled_CH₄, modeled_oceanCO₂_flux, modeled_temperature, modeled_ocean_heat)

        #---------------------------------------------------------------------------
        # Global Surface Temperature (normalized to 1861-1880 mean) Log-Likelihood.
        #---------------------------------------------------------------------------

        llik_temperature = 0.0

        # Calculate temperature residuals.
        for (i, index)=enumerate(indices_temperature_data)
            temperature_residual[i] = calibration_data[index, :hadcrut_temperature_obs] - modeled_temperature[index]
        end

        # Calculate temperature log-likelihood.
        llik_temperature = hetero_logl_ar1(temperature_residual, σ_temperature, ρ_temperature, calibration_data[indices_temperature_data,:hadcrut_temperature_sigma])

        #-----------------------------------------------------------------------
        # Ocean Heat Content Log-Likelihood
        #-----------------------------------------------------------------------

        llik_ocean_heat = 0.0

        # Calculate ocean heat residuals.
        for (i, index)=enumerate(indices_oceanheat_data)
            ocean_heat_residual[i] = calibration_data[index, :ocean_heat_obs] - modeled_ocean_heat[index]
        end

        # Calculate ocean heat log-likelihood.
        llik_ocean_heat = hetero_logl_ar1(ocean_heat_residual, σ_ocean_heat, ρ_ocean_heat, calibration_data[indices_oceanheat_data, :ocean_heat_sigma])

        #-----------------------------------------------------------------------
        # Atmospheric CO₂ Concentration Log-Likelihood
        #-----------------------------------------------------------------------

        llik_co2 = 0.

        # Calculate CO₂ concentration (Law Dome) residuals (assuming 8 year model mean centered on year of ice core observation).
        for (i, index)=enumerate(indices_lawdome_co2_data)
            co2_residual[i] = calibration_data[index, :lawdome_co2_obs] - mean(modeled_CO₂[index .+ (-4:3)])
        end

        # Calculate CO₂ concentration (Mauna Loa) residuals.
        for (i, index)=enumerate(indices_maunaloa_co2_data)
            co2_residual[i+n_lawdome_co2] = calibration_data[index, :maunaloa_co2_obs] - modeled_CO₂[index]
        end

        # Calculate atmospheric CO₂ concentration log-likelihood.
        llik_co2 = hetero_logl_car1(co2_residual, indices_co2_data, σ²_white_noise_CO₂, α₀_CO₂, calibration_data[indices_co2_data, :co2_combined_sigma])

        #-----------------------------------------------------------------------
        # Atmospheric CH₄ Concentration Log-Likelihood
        #-----------------------------------------------------------------------

        llik_ch4 = 0.

        # Calculate CH₄ concentration (Law Dome) residuals for individual data points (assuming 8 year model mean centered on year of ice core observation).
        for (i,index) = enumerate(indices_lawdome_ch4_data)
            ch4_residual[i] = calibration_data[index, :lawdome_ch4_obs] - mean(modeled_CH₄[index .+ (-4:3)])
        end

        # Calculate CH₄ concentration (NOAA) residuals.
        for (i, index)=enumerate(indices_noaa_ch4_data)
            ch4_residual[i+n_lawdome_ch4] = calibration_data[index, :noaa_ch4_obs] - modeled_CH₄[index]
        end

        # Calculate atmospheric CH₄ concentration log-likelihood.
        llik_ch4 = hetero_logl_car1(ch4_residual, indices_ch4_data, σ²_white_noise_CH₄, α₀_CH₄, calibration_data[indices_ch4_data, :ch4_combined_sigma])

        #-----------------------------------------------------------------------
        # Ocean Carbon Flux Log-Likelihood
        #-----------------------------------------------------------------------

        llik_oceanco2_flux = 0.

        # Calculate ocean CO₂ flux log-likelihood for individual data points.
        for (i,index) = enumerate(indices_oceanco2_flux_data)
            oceanco2_flux_single_llik[i] = logpdf(Normal(modeled_oceanCO₂_flux[index], calibration_data[index, :oceanco2_flux_sigma]), calibration_data[index, :oceanco2_flux_obs])
        end

        # Calculate ocean CO₂ flux total log-likelihood as sum of individual data point likelihoods.
        llik_oceanco2_flux = sum(oceanco2_flux_single_llik)

        #-----------------------------------------------------------------------
        # Total Log-Likelihood
        #-----------------------------------------------------------------------

        # Calculate the total log-likelihood (assuming residual independence across data sets).
        llik = llik_temperature + llik_ocean_heat + llik_co2 + llik_ch4 + llik_oceanco2_flux

        return llik
    end

    #---------------------------------------------------------------------------------------------------------------
    # Create a function to calculate the log-posterior of uncertain parameters 'p' (posterior ∝ likelihood * prior)
    #---------------------------------------------------------------------------------------------------------------

    function log_posterior(p)

        # Calculate log-prior
        log_prior = total_log_prior(p)

        # In case a parameter sample leads to non-physical model outcomes, return -Inf rather than erroring out.
        try
            log_post = isfinite(log_prior) ? total_log_likelihood(p) + log_prior : -Inf
        catch
            log_post = - Inf
        end
    end

    # Return log posterior function given user specifications.
    return log_posterior
end


text:
struct pyconvert_rule_numpysimplevalue{R,S} <: Function end

function (::pyconvert_rule_numpysimplevalue{R,SAFE})(::Type{T}, x::Py) where {R,SAFE,T}
    ans = GC.@preserve x C.PySimpleObject_GetValue(R, getptr(x))
    if SAFE
        pyconvert_return(convert(T, ans))
    else
        pyconvert_tryconvert(T, ans)
    end
end

const NUMPY_SIMPLE_TYPES = [
    ("int8", Int8),
    ("int16", Int16),
    ("int32", Int32),
    ("int64", Int64),
    ("uint8", UInt8),
    ("uint16", UInt16),
    ("uint32", UInt32),
    ("uint64", UInt64),
    ("float16", Float16),
    ("float32", Float32),
    ("float64", Float64),
    ("complex32", ComplexF16),
    ("complex64", ComplexF32),
    ("complex128", ComplexF64),
]

function init_numpy()
    for (t,T) in NUMPY_SIMPLE_TYPES
        isint = occursin("int", t)
        isuint = occursin("uint", t)
        isfloat = occursin("float", t)
        iscomplex = occursin("complex", t)
        isreal = isint || isfloat
        isnumber = isreal || iscomplex

        name = "numpy/$t"
        rule = pyconvert_rule_numpysimplevalue{T, false}()
        saferule = pyconvert_rule_numpysimplevalue{T, true}()

        pyconvert_add_rule(name, T, saferule, 100)
        isuint && pyconvert_add_rule(name, UInt, sizeof(T) ≤ sizeof(UInt) ? saferule : rule)
        isuint && pyconvert_add_rule(name, Int, sizeof(T) < sizeof(Int) ? saferule : rule)
        isint && !isuint && pyconvert_add_rule(name, Int, sizeof(T) ≤ sizeof(Int) ? saferule : rule)
        isint && pyconvert_add_rule(name, Integer, rule)
        isfloat && pyconvert_add_rule(name, Float64, saferule)
        isreal && pyconvert_add_rule(name, Real, rule)
        iscomplex && pyconvert_add_rule(name, ComplexF64, saferule)
        iscomplex && pyconvert_add_rule(name, Complex, rule)
        isnumber && pyconvert_add_rule(name, Number, rule)
    end
end


text:
"""
object which holds a transformation matrix and a model reference
"""
type Transform
  model::Mat4x4f
  mvp::Mat4x4f
end

"""
array which holds raw data of mesh
"""
type MeshArray
  elements::UInt32
  count::UInt32
  bufferID::GLuint
  bufferType::GLenum
  loaded::Bool

  data::AbstractArray

  MeshArray() = new(1,0,0,GL_ARRAY_BUFFER,false,[])
end

"""
data which holds various arrays of mesh
"""
type MeshData
  vao::GLuint
  arrays::Dict{Symbol,MeshArray}
  draw::Union{Void,MeshArray}

  MeshData() = new(0,Dict(),nothing)
end

gltypes=Dict(Float32=>GL_FLOAT,Float64=>GL_DOUBLE,UInt32=>GL_UNSIGNED_INT,Int32=>GL_INT)

"""
set attributes for shader
"""
function setAttributes(this::MeshArray, program, attrb)
  if this.count == 0 return end

  const STRIDE = GLsizei(length(attrb)<=1?0:reduce(+, (x->sizeof(x[2])*x[3]).(attrb)))
  OFFSET =  C_NULL

  glBindBuffer(this.bufferType, this.bufferID)

  for (name,typ,elems,inst) in attrb
    atr = glGetAttribLocation(program, name)
    glCheckError("glGetAttribLocation")
    if atr > -1
      glEnableVertexAttribArray(atr)
      glCheckError("glEnableVertexAttribArray")
      glVertexAttribPointer(atr, elems, gltypes[typ], GL_FALSE, STRIDE, OFFSET)
      glCheckError("glVertexAttribPointer")
      if inst>0 glVertexAttribDivisor(atr, inst); glCheckError("glVertexAttribDivisor") end
      OFFSET +=  sizeof(typ)*elems
    else warn("Could not load Attribute \"$name\"")
    end
  end

  glBindBuffer(this.bufferType, 0)
end

"""
creates gpu buffer
"""
function createBuffers(this::MeshData)
  if this.vao == 0 this.vao = glGenVertexArray() end
  for (s,a) in this.arrays if a.bufferID == 0 && a.count > 0 a.bufferID = glGenBuffer() end end
end

"""
sets attributes for shader
"""
function setAttributes(this::MeshData, program)
  glBindVertexArray(this.vao)
  glCheckError("glBindVertexArray")
  if haskey(this.arrays,:vertices) setAttributes(this.arrays[:vertices], program, [("iVertex",Float32,3,0)]) end
  if haskey(this.arrays,:instances) setAttributes(this.arrays[:instances], program, [("iInstancePos",Float32,3,1),("iInstanceFlags",Float32,2,1)]) end
  glBindVertexArray(0)
end

"""
sets array which will be drawn by opengl
"""
function setDrawArray(this::MeshData, key::Symbol)
  this.draw = this.arrays[key]
end

"""
sets data
"""
function setData(this::MeshArray, data, elems=0)
  if elems != 0 this.elements = elems end
  this.data = data
  this.count = length(this.data)/this.elements
end

"""
links data
"""
function linkData(this::MeshData, args...)
  d=Dict(args)

  this.arrays = Dict()
  this.draw = nothing

  for (s,x) in d
    l=isa(x,Tuple)?length(x):0

    data = l>0?x[1]:x
    dtyp = eltype(data)
    elems = l>1?x[2]:1
    btyp = l>2?x[3]:GL_ARRAY_BUFFER
    draw = l>3?x[4]:false

    if haskey(this.arrays,s)
      a = this.arrays[s]
      a.loaded = false
    else
      this.arrays[s] = a = MeshArray()
      a.bufferType = btyp
    end

    setData(a, data, elems)

    if draw this.draw = a
    elseif this.draw == nothing this.draw = a
    end
  end

  createBuffers(this)
  upload(this)
end

"""
uploads data
"""
function upload(this::MeshArray)
  if this.loaded || this.count == 0 return end
  const data = this.data
  glBindBuffer(this.bufferType, this.bufferID)
  glCheckError("glBindBuffer")
  glBufferData(this.bufferType, sizeof(data), data, GL_STATIC_DRAW)
  glCheckError("glBufferData")
  this.loaded=true
end

"""
uploads data
"""
function upload(this::MeshData)
  glBindVertexArray(this.vao)
  glCheckError("glBindVertexArray")
  for (s,a) in this.arrays upload(a) end
  glBindVertexArray(0)
end

"""
uploads data
"""
function upload(this::MeshData, key::Symbol, data::AbstractArray)
   if !haskey(this.arrays,key) warn("Could not find key $key"); return end
  a=this.arrays[key]
  setData(a, data)
  a.loaded = false
  glBindVertexArray(this.vao)
  glCheckError("glBindVertexArray")
  upload(a)
  glBindVertexArray(0)
end

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
function ContinuousSpace(
    extent,
    spacing;
    kwargs...,
)
    @warn "Giving `spacing` as a positional argument to `ContinuousSpace` is "*
    "deprecated, provide it as a keyword instead."
    return ContinuousSpace(extend; spacing, kwargs...)
end

function __init__()
    # Plot recipes
    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
        include("visualization/plot-recipes.jl")
    end
    # Workaround for Documenter.jl, so we don't need to include
    # heavy dependencies to build documentation
    @require Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" begin
        include("visualization/plot-recipes.jl")
    end
end

# 4.0 Depreciations
@deprecate space_neighbors nearby_ids
@deprecate node_neighbors nearby_positions
@deprecate get_node_contents ids_in_position
@deprecate get_node_agents agents_in_position
@deprecate pick_empty random_empty
@deprecate find_empty_nodes empty_positions
@deprecate has_empty_nodes has_empty_positions
@deprecate nodes positions

# 4.2 Deprecations
@deprecate fastest Schedulers.fastest
@deprecate by_id Schedulers.by_id
@deprecate as_added Schedulers.by_id
@deprecate random_activation Schedulers.randomly
@deprecate partial_activation Schedulers.partially
@deprecate property_activation Schedulers.by_property
@deprecate by_type Schedulers.by_type

@deprecate osm_random_road_position OSM.random_road_position
@deprecate osm_plan_route OSM.plan_route
@deprecate osm_map_coordinates OSM.map_coordinates
@deprecate osm_road_length OSM.road_length
@deprecate osm_random_route! OSM.random_route!
@deprecate osm_latlon OSM.latlon
@deprecate osm_intersection OSM.intersection
@deprecate osm_road OSM.road
@deprecate move_agent!(a, model::ABM{<:OpenStreetMapSpace}, distance) move_along_route!(a, model, distance)


text:
using Test

using JuLIP, NBodyIPs
using NBodyIPs.Polys: NBPoly
using NBodyIPs: BondLengthDesc
const Env = NBodyIPs.EnvIPs

println("-------------------")
println("Tests for EnvIPs")
println("-------------------")
println("Setting up the test systems ...")
r0 = rnn(:Cu)
rcut3 = 2.1 * r0
at = rattle!(bulk(:Cu, cubic=true) * 2, 0.02)
D3 = BondLengthDesc(ExpTransform(3, r0),
                    CosCut(0.66*rcut3, rcut3) )

random_3body(ntup=1) = NBPoly( [tuple( [rand(1:4, 3); 0]... ) for n = 1:ntup],
                                 (0.1.+rand(ntup))/factorial(3), D3 )

Vn = ("exp(- 3 * ((r/$r0)-1))", 1.8*r0)

V3env_0 = Env.EnvIP(0, random_3body(), Vn...)
V3env_1 = Env.EnvIP(1, random_3body(), Vn...)
V3env_2 = Env.EnvIP(2, random_3body(), Vn...)

println("Finite-difference test on configurations")
println("------------------------------------------------")
set_constraint!(at, VariableCell(at))
for Venv in [V3env_0, V3env_1, V3env_2]
   (@test JuLIP.Testing.fdtest(Venv, at)) |> println
end

println("Testing combined assembly of basis")
B = Env.envpolys(3, D3, 5, Vn, 2)

# nb: the , false is to set typewarn=false, since eltype(B) is not a leaftype
print("energy: ")
E1 = energy(B, at, false)
E2 = [ energy(b, at) for b in B ]
println(@test E1 ≈ E2)

print("forces: ")
F1 = forces(B, at, false)
F2 = [ forces(b, at) for b in B ]
println(@test F1 ≈ F2)

print("virial: ")
V1 = virial(B, at, false)
V2 = [ virial(b, at) for b in B ]
println(@test V1 ≈ V2)

println("Check saving and loading of an EnvIP")
IP = NBodyIP(B, rand(length(B)))
fname = tempname() * ".json"
save_ip(fname, IP)
IP1, _ = load_ip(fname)
println(@test IP1 == IP)
run(`rm $fname`)

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs