text:
using OPFInstances, BenchmarkTools, DataStructures

println("Loading module MathProgComplex")
tic();
using MathProgComplex
toc();

function main(instance = case30pwl)
    instancepath = getinstancepath("Matpower", "QCQP", instance)
    # instancepath = getinstancepath("Matpower", "QCQP", "case89pegase")
    # instancepath = getinstancepath("Matpower", "QCQP", "case300")
    println("\n\nWorking on $(splitdir(instancepath))")

    setlevel!(getlogger(MathProgComplex), "info")
    @show MathProgComplex.DictType

    pb_c, pt = import_from_dat(instancepath)
    problem = pb_cplx2real(pb_c)

    relax_ctx = set_relaxation(problem; hierarchykind=:Real,
                                        # symmetries=[PhaseInvariance],
                                        d = 1,
                                        params = Dict(:opt_outlev=>0,
                                                      :opt_outmode=>0,
                                                      :opt_outcsv=>0))

    max_cliques = get_maxcliques(relax_ctx, problem)
    relax_ctx.relaxparams[:opt_nb_cliques] = length(max_cliques)

    momentmat_param, localizingmat_param = build_sparsity(relax_ctx, problem, max_cliques)

    println("--- time build_sparsity")
    @btime build_sparsity($relax_ctx, $problem, $max_cliques)

    mmtrel_pb = build_momentrelaxation(relax_ctx, problem, momentmat_param, localizingmat_param, max_cliques)

    println("--- time build_momentrelaxation")
    @btime build_momentrelaxation($relax_ctx, $problem, $momentmat_param, $localizingmat_param, $max_cliques)

    sdpinstance = build_SOSrelaxation(relax_ctx, mmtrel_pb)

    println("--- time build_SOSrelaxation")
    @btime build_SOSrelaxation($relax_ctx, $mmtrel_pb)

    path = joinpath(pwd(), "Mosek_runs", "worksdp")
    mkpath(path)
    export_SDP(sdpinstance, path)
    sdp_instance = read_SDPPrimal(path)

    println("--- time read_SDPPrimal")
    @btime read_SDPPrimal($path)

    # sdp = SDP_Problem()

    # set_constraints!(sdp, sdp_instance)
    # set_vartypes!(sdp, sdp_instance)
    # set_blocks!(sdp, sdp_instance)
    # set_linvars!(sdp, sdp_instance)

    # set_matrices!(sdp, sdp_instance)
    # set_linear!(sdp, sdp_instance)
    # set_const!(sdp, sdp_instance)

    # # println(sdp)

    # primal = SortedDict{Tuple{String,String,String}, Float64}()
    # dual = SortedDict{Tuple{String, String, String}, Float64}()

    # primobj, dualobj = solve_mosek(sdp::SDP_Problem, primal, dual, sol_info=relax_ctx.relaxparams)

    # # final_output(relax_ctx)

    return
end

function main_work()
    instancepath = getinstancepath("Matpower", "QCQP", "case300")
    println("\nWorking on $(splitdir(instancepath))")

    setlevel!(getlogger(MathProgComplex), "debug")
    @show MathProgComplex.DictType

    pb_c, pt = import_from_dat(instancepath)

    pb_cplx2real(pb_c)
    setlevel!(getlogger(MathProgComplex), "info")

    return
end

# main_work()
main("WB2")
main("case30pwl")
main("case89pegase")
main("case300")


text:
module MultivariateStats_

export RidgeRegressor, PCA, KernelPCA, ICA, LDA, BayesianLDA

import MLJModelInterface
import MLJModelInterface: @mlj_model, metadata_pkg, metadata_model,
    Table, Continuous, Finite

using MLJBase # need the UnivariateFinite *type* not just the method

const MMI = MLJModelInterface

import StatsBase: proportions, CovarianceEstimator
using Distances
using LinearAlgebra
using Tables

import MultivariateStats

const MS = MultivariateStats

struct LinearFitresult{F} <: MMI.MLJType
    coefficients::Vector{F}
    intercept::F
end

const RIDGE_DESCR = "Ridge regressor with regularization parameter lambda. Learns a linear regression with a penalty on the l2 norm of the coefficients."
const PCA_DESCR = "Principal component analysis. Learns a linear transformation to project the data  on a lower dimensional space while preserving most of the initial variance."
const KPCA_DESCR = "Kernel principal component analysis."
const ICA_DESCR = "Independent component analysis."
const LDA_DESCR = "Multiclass linear discriminant analysis. The algorithm learns a projection matrix `W` that projects the feature matrix `Xtrain` onto a lower dimensional space of dimension `out_dim` such that the between-class variance in the transformed space is maximized relative to the within-class variance."
const BayesianLDA_DESCR = "Bayesian Multiclass linear discriminant analysis. The algorithm learns a projection matrix `W` that projects the feature matrix `Xtrain` onto a lower dimensional space of dimension `out_dim` such that the between-class variance in the transformed space is maximized relative to the within-class variance and classifies using Bayes rule."
const BayesianSubspaceLDA_DESCR = "Bayesian Multiclass linear discriminant analysis. Suitable for high dimensional data. The algorithm learns a projection matrix `W` (without computing covariance matrices Sw ,Sb) in the subspace spanned by the within class scatter matrix, that projects the feature matrix `Xtrain` onto a lower dimensional space of dimension `out_dim` such that the between-class variance in the transformed space is maximized relative to the within-class variance and classifies using Bayes rule."
const SubspaceLDA_DESCR = "Multiclass linear discriminant analysis. Suitable for high dimensional data. The algorithm learns a projection matrix `W` (without computing covariance matrices Sw ,Sb) in the subspace spanned by the within class scatter matrix, that projects the feature matrix `Xtrain` onto a lower dimensional space of dimension `out_dim` such that the between-class variance in the transformed space is maximized relative to the within-class variance."

####
#### RIDGE
####

"""
    RidgeRegressor(; lambda=1.0)

$RIDGE_DESCR

## Parameters

* `lambda=1.0`: non-negative parameter for the regularization strength.
"""
@mlj_model mutable struct RidgeRegressor <: MMI.Deterministic
    lambda::Real = 1.0::(_ ≥ 0)
end

function MMI.fit(model::RidgeRegressor, verbosity::Int, X, y)
    Xmatrix   = MMI.matrix(X)
    features  = Tables.schema(X).names
    θ         = MS.ridge(Xmatrix, y, model.lambda)
    coefs     = θ[1:end-1]
    intercept = θ[end]

    fitresult = LinearFitresult(coefs, intercept)
    report    = NamedTuple()
    cache     = nothing

    return fitresult, cache, report
end

MMI.fitted_params(::RidgeRegressor, fr) =
    (coefficients=fr.coefficients, intercept=fr.intercept)

function MMI.predict(::RidgeRegressor, fr, Xnew)
    Xmatrix = MMI.matrix(Xnew)
    return Xmatrix * fr.coefficients .+ fr.intercept
end

####
#### PCA
####

const PCAFitResultType = MS.PCA

"""
    PCA(; maxoutdim=nothing, method=:auto, pratio=0.99, mean=nothing)

$PCA_DESCR

## Parameters

* `maxoutdim=nothing`: maximum number of output dimensions, unconstrained if
                       nothing.
* `method=:auto`:      method to use to solve the problem, one of `:auto`,
                       `:cov` or `:svd`
* `pratio=0.99`:       ratio of variance preserved
* `mean=nothing`:      if set to nothing centering will be computed and
                       applied, if set to `0` no centering (assumed pre-
                       centered), if a vector is passed, the centering is done
                       with that vector.
"""
@mlj_model mutable struct PCA <: MMI.Unsupervised
    maxoutdim::Union{Nothing,Int} = nothing::(_ === nothing || _ ≥ 1)
    method::Symbol  = :auto::(_ in (:auto, :cov, :svd))
    pratio::Float64 = 0.99::(0.0 < _ ≤ 1.0)
    mean::Union{Nothing, Real, Vector{Float64}} = nothing::(_ === nothing || (_ isa Real && iszero(_)) || true)
end

function MMI.fit(model::PCA, verbosity::Int, X)
    Xarray = MMI.matrix(X)
    mindim = minimum(size(Xarray))

    maxoutdim = model.maxoutdim === nothing ? mindim : model.maxoutdim

    # NOTE: copy/transpose
    fitresult = MS.fit(MS.PCA, permutedims(Xarray);
                       method=model.method,
                       pratio=model.pratio,
                       maxoutdim=maxoutdim,
                       mean=model.mean)

    cache = nothing
    report = (indim=MS.indim(fitresult),
              outdim=MS.outdim(fitresult),
              mean=MS.mean(fitresult),
              principalvars=MS.principalvars(fitresult),
              tprincipalvar=MS.tprincipalvar(fitresult),
              tresidualvar=MS.tresidualvar(fitresult),
              tvar=MS.tvar(fitresult))

    return fitresult, cache, report
end

MMI.fitted_params(::PCA, fr) = (projection=fr,)

function MMI.transform(::PCA, fr::PCAFitResultType, X)
    # X is n x d, need to transpose and copy twice...
    Xarray = MMI.matrix(X)
    Xnew   = permutedims(MS.transform(fr, permutedims(Xarray)))
    return MMI.table(Xnew, prototype=X)
end

####
#### KernelPCA
####

const KernelPCAFitResultType = MS.KernelPCA

default_kernel = (x, y) -> x'y

"""
    KernelPCA(; maxoutdim=nothing, kernel=(x,y)->x'y, solver=:auto, pratio=0.99, mean=nothing)

$KPCA_DESCR

## Parameters

* `maxoutdim=nothing`: maximum number of output dimensions, unconstrained if
                       nothing.
* `kernel=nothing`:    kernel function of 2 vector arguments x and y, returns a
                       scalar value, (x,y)->x'y if nothing
* `solver=:auto`:      solver to use for the eigenvalues, one of `:eig`
                       (default), `:eigs`
* `inverse=false`:     perform calculation for inverse transform
* `beta=1.0`:          strength of the ridge regression that learns the inverse
                       transform when inverse is true
* `tol=0.0`:           Convergence tolerance for eigs solver
* `maxiter=300`:       maximum number of iterations for eigs solver
"""
@mlj_model mutable struct KernelPCA <: MMI.Unsupervised
    maxoutdim::Union{Nothing,Int} = nothing::(_ === nothing || _ ≥ 1)
    kernel::Function = default_kernel
    solver::Symbol   = :eig::(_ in (:eig, :eigs))
    inverse::Bool    = false
    beta::Real       = 1.0::(_ ≥ 0.0)
    tol::Real        = 1e-6::(_ ≥ 0.0)
    maxiter::Int     = 300::(_ ≥ 1)
end

function MMI.fit(model::KernelPCA, verbosity::Int, X)
    Xarray = MMI.matrix(X)
    mindim = minimum(size(Xarray))
    # default max out dim if not given
    maxoutdim = model.maxoutdim === nothing ? mindim : model.maxoutdim

    fitresult = MS.fit(MS.KernelPCA, permutedims(Xarray);
                       kernel=model.kernel,
                       maxoutdim=maxoutdim,
                       solver=model.solver,
                       inverse=model.inverse,
                       β=model.beta,
                       tol=model.tol,
                       maxiter=model.maxiter)

    cache  = nothing
    report = (indim=MS.indim(fitresult),
              outdim=MS.outdim(fitresult),
              projection=MS.projection(fitresult),
              principalvars=MS.principalvars(fitresult))

    return fitresult, cache, report
end

MMI.fitted_params(::KernelPCA, fr) = (projection=fr,)

function MMI.transform(::KernelPCA, fr::KernelPCAFitResultType, X)
    # X is n x d, need to transpose and copy twice...
    Xarray = MMI.matrix(X)
    Xnew   = permutedims(MS.transform(fr, permutedims(Xarray)))
    return MMI.table(Xnew, prototype=X)
end

####
#### ICA
####

const ICAFitResultType = MS.ICA

"""
    ICA(; k=0, alg=:fastica, fun=:tanh, do_whiten=true, maxiter=100,
      tol=1e-6, mean=nothing, winit=zeros(0,0))

$ICA_DESCR

## Parameters

* `k=0`:              number of independent components to recover, set
                      automatically if `0`
* `alg=:fastica`:     algorithm to use (only `:fastica` is supported at the
                      moment)
* `fun=:tanh`:        approximate neg-entropy functor, via the function
                      `MultivariateStats.icagfun`, one of `:tanh` and `:gaus`
* `do_whiten=true`:   whether to perform pre-whitening
* `maxiter=100`:      maximum number of iterations
* `tol=1e-6`:         convergence tolerance for change in matrix W
* `mean=nothing`:     mean to use, if nothing centering is computed and
                      applied, if zero, no centering, a vector of means can be
                      passed
* `winit=zeros(0,0)`: initial guess for matrix `W` either an empty matrix
                      (random init), a matrix of size `k × k` (if `do_whiten`
                      is true), a matrix of size `m × k` otherwise
"""
@mlj_model mutable struct ICA <: MMI.Unsupervised
    k::Int          = 0::(_ ≥ 0)
    alg::Symbol     = :fastica::(_ in (:fastica,))
    fun::Symbol     = :tanh::(_ in (:tanh, :gaus))
    do_whiten::Bool = true
    maxiter::Int    = 100
    tol::Real       = 1e-6::(_ ≥ 0.0)
    winit::Union{Nothing,Matrix{<:Real}} = nothing
    mean::Union{Nothing, Real, Vector{Float64}} = nothing::(_ === nothing || (_ isa Real && iszero(_)) || true)
end

function MMI.fit(model::ICA, verbosity::Int, X)
    Xarray = MMI.matrix(X)
    n, p   = size(Xarray)

    m = min(n, p)
    k = ifelse(model.k ≤ m, model.k, m)

    fitresult = MS.fit(MS.ICA, permutedims(Xarray), k;
                       alg=model.alg,
                       fun=MS.icagfun(model.fun),
                       do_whiten=model.do_whiten,
                       maxiter=model.maxiter,
                       tol=model.tol,
                       mean=model.mean,
                       winit=ifelse(model.winit === nothing, zeros(0,0), model.winit))

    cache = nothing
    report = (indim=MS.indim(fitresult),
              outdim=MS.outdim(fitresult),
              mean=MS.mean(fitresult))

    return fitresult, cache, report
end

MMI.fitted_params(::ICA, fr) = (projection=fr,)

function MMI.transform(::ICA, fr::ICAFitResultType, X)
    # X is n x d, need to transpose and copy twice...
    Xarray = MMI.matrix(X)
    Xnew   = permutedims(MS.transform(fr, permutedims(Xarray)))
    return MMI.table(Xnew, prototype=X)
end

####
#### MulticlassLDA
####

"""
    LDA(; kwargs...)

$LDA_DESCR

## Parameters

* `method=:gev`:      choice of solver, one of `:gevd` or `:whiten` methods
* `cov_w=SimpleCovariance()`: an estimator for the within-class covariance, by
                              default set to the standard
                              `MultivariateStats.CovarianceEstimator` but could
                              be set to any robust estimator from `CovarianceEstimation.jl`.
* `cov_b=SimpleCovariance()`: same as `cov_w` but for the between-class
                              covariance.
* `out_dim`:          the output dimension, i.e dimension of the transformed
                      space, automatically set if 0 is given (default).
* `regcoef`:          regularization coefficient (default value 0.0). A
                      positive value `regcoef * eigmax(Sw)` where `Sw` is the
                      within-class covariance estimator, is added to the
                      diagonal of Sw to improve numerical stability. This can
                      be useful if using the standard covariance estimator.
* `dist=SqEuclidean`: the distance metric to use when performing classification
                      (to compare the distance between a new point and
                      centroids in the transformed space), an alternative
                      choice can be the `CosineDist`.

See also the [package documentation](https://multivariatestatsjl.readthedocs.io/en/latest/lda.html).
For more information about the algorithm, see the paper by Li, Zhu and Ogihara, [Using Discriminant Analysis for Multi-class Classification: An Experimental Investigation](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.89.7068&rep=rep1&type=pdf).
"""
@mlj_model mutable struct LDA <: MMI.Probabilistic
    method::Symbol   = :gevd::(_ in (:gevd, :whiten))
    cov_w::CovarianceEstimator = MS.SimpleCovariance()
    cov_b::CovarianceEstimator = MS.SimpleCovariance()
    out_dim::Int     = 0::(_ ≥ 0)
    regcoef::Float64    = 1e-6::(_ ≥ 0)
    dist::SemiMetric = SqEuclidean()
end

function MMI.fit(model::LDA, ::Int, X, y)
    class_list = MMI.classes(y[1]) #class list containing unique entries in y
    nclasses   = length(class_list)

    # NOTE: copy/transpose
    Xm_t   = MMI.matrix(X, transpose=true) # now p x n matrix
    yplain = MMI.int(y) # vector of n ints in {1,..., nclasses}
    p      = size(Xm_t, 1)

     #check to make sure we have more than one class
    nclasses >= 2 ||  throw(ArgumentError("The number of classes has to be greater than one"))

    # check output dimension default is min(p, nc-1)
    def_outdim = min(p, nclasses - 1)
    # if unset (0) use the default; otherwise try to use the provided one
    out_dim = ifelse(model.out_dim == 0, def_outdim, model.out_dim)
    # check if the given one is sensible
    out_dim ≤ def_outdim || throw(ArgumentError("`out_dim` must not be larger than `min(p, nc-1)` where `p` is the dimension of `X` and `nc` is the number of classes."))

    core_res = MS.fit(MS.MulticlassLDA, nclasses, Xm_t, Int.(yplain);
                      method=model.method,
                      outdim=out_dim,
                      regcoef=model.regcoef,
                      covestimator_within=model.cov_w,
                      covestimator_between=model.cov_b)

    cache     = nothing
    report    = NamedTuple{}()
    fitresult = (core_res, class_list)

    return fitresult, cache, report
end

function MMI.fitted_params(::LDA, (core_res, class_list))
    return (class_means       = MS.classmeans(core_res),
            projection_matrix = MS.projection(core_res))
end

function MMI.predict(m::LDA, (core_res, class_list), Xnew)
    # projection of Xnew XWt is n x o  where o = number of out dims
    XWt = MMI.matrix(Xnew) * core_res.proj
    # centroids in the transformed space, nc x o
    centroids = permutedims(core_res.pmeans)

    # compute the distances in the transformed space between pairs of rows
    # the probability matrix is `n x nc` and normalised accross rows
    P = pairwise(m.dist, XWt, centroids, dims=1)
    # apply a softmax transformation
    P .-= maximum(P, dims=2)
    P  .= exp.(-P)
    P ./= sum(P, dims=2)

    return MMI.UnivariateFinite(class_list, P)
end

####
#### BayesianLDA
####

"""
    BayesianLDA(; kwargs...)

$BayesianLDA_DESCR

## Parameters

* `method=:gev`:    choice of solver, one of `:gevd` or `:whiten`
                    methods
* `cov_w=SimpleCovariance()`: an estimator for the within-class covariance, by
                              default set to the standard
                              `MultivariateStats.CovarianceEstimator` but could
                              be set to any robust estimator from
                              `CovarianceEstimation.jl`.
* `cov_b=SimpleCovariance()`: same as `cov_w` but for the between-class
                              covariance.
* `out_dim`:                  the output dimension, i.e dimension of the
                    transformed space, automatically set if 0 is
                    given (default).
* `regcoef`:        regularization coefficient (default value 1e-6).
                    A positive value `regcoef * eigmax(Sw)` where `Sw` is the
                    within-class covariance estimator, is added to the diagonal
                    of Sw to improve numerical stability. This can be useful if
                    using the standard covariance estimator.
* `priors=nothing`: if `priors = nothing` estimates the prior probabilities
                    from the data else it uses the user specified
                    `UnivariateFinite` prior probabilities.

See also the [package documentation](https://multivariatestatsjl.readthedocs.io/en/latest/lda.html).
For more information about the algorithm, see the paper by Li, Zhu and Ogihara, [Using Discriminant Analysis for Multi-class Classification: An Experimental Investigation](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.89.7068&rep=rep1&type=pdf).
"""
@mlj_model mutable struct BayesianLDA <: MMI.Probabilistic
    method::Symbol   = :gevd::(_ in (:gevd, :whiten))
    cov_w::CovarianceEstimator = MS.SimpleCovariance()
    cov_b::CovarianceEstimator = MS.SimpleCovariance()
    out_dim::Int     = 0::(_ ≥ 0)
    regcoef::Float64    = 1e-6::(_ ≥ 0)
    priors::Union{Nothing,MLJBase.UnivariateFinite} = nothing
end

function MMI.fit(model::BayesianLDA, ::Int, X, y)
    class_list = MMI.classes(y[1]) #class list containing unique entries in y
    nclasses   = length(class_list)

    # NOTE: copy/transpose
    Xm_t   = MMI.matrix(X, transpose=true) # now p x n matrix
    yplain = MMI.int(y) # vector of n ints in {1,..., nclasses}
    p, n   = size(Xm_t)

     #check to make sure we have more than one class
    nclasses >= 2 ||  throw(ArgumentError("The number of classes has to be greater than one"))

    # check output dimension default is min(p, nc-1)
    def_outdim = min(p, nclasses - 1)
    # if unset (0) use the default; otherwise try to use the provided one
    out_dim = ifelse(model.out_dim == 0, def_outdim, model.out_dim)
    # check if the given one is sensible
    out_dim ≤ def_outdim || throw(ArgumentError("`out_dim` must not be larger than `min(p, nc-1)` where `p` is the dimension of `X` and `nc` is the number of classes."))

    ## Estimates prior probabilities if unspecified by user.
    if model.priors == nothing
        priors = proportions(yplain)
    else
        #check if the length of priors is same as nclasses
        length(MMI.levels(model.priors)) == nclasses || throw(ArgumentError("Invalid size of `priors`"))
        priors = MMI.pdf.(model.priors, class_list)
    end

    core_res = MS.fit(MS.MulticlassLDA, nclasses, Xm_t, Int.(yplain);
                      method=model.method,
                      outdim=out_dim,
                      regcoef= model.regcoef,
                      covestimator_within=model.cov_w,
                      covestimator_between=model.cov_b)

    ## The original projection matrix satisfies Pᵀ*Sw*P=I
    ## scaled projection_matrix and core_res.proj by multiplying by sqrt(n - nclasses) this ensures Pᵀ*Σ*P=I
    ## where covariance estimate Σ = Sw / (n - nclasses)
    core_res.proj   .*= sqrt(n - nclasses)
    core_res.pmeans .*= sqrt(n - nclasses)

    cache     = nothing
    report    = NamedTuple{}()
    fitresult = (core_res, class_list, priors)

    return fitresult, cache, report
end

function MMI.fitted_params(::BayesianLDA, (core_res, class_list, priors))
    return (class_means       = MS.classmeans(core_res),
            projection_matrix = MS.projection(core_res),
            priors            = priors)
end

function MMI.predict(m::BayesianLDA, (core_res, class_list, priors), Xnew)
    # projection of Xnew XWt is n x o  where o = number of out dims
    XWt = MMI.matrix(Xnew) * core_res.proj
    # centroids in the transformed space, nc x o
    centroids = permutedims(core_res.pmeans)

    n  = core_res.stats.tweight # n is the Number of training examples
    nclasses = size(class_list)

    # compute the distances in the transformed space between pairs of rows
    # The discriminant matrix `P` is of dimension `n x nc`
    #  P[i,k] = -0.5*(xᵢ −  µₖ)ᵀΣ⁻¹(xᵢ −  µₖ) + log(priorsₖ) and Σ⁻¹ = I due to the nature of the projection_matrix
    P = pairwise(SqEuclidean(), XWt, centroids, dims=1)
    P .*= -0.5
    P .+= log.(priors)'

    # apply a softmax transformation to convert P to a probability matrix
    P .-= maximum(P, dims=2)
    P  .= exp.(P)
    P ./= sum(P, dims=2)

    return [MMI.UnivariateFinite(class_list, P[j, :]) for j in 1:size(P, 1)]
end

####
#### BayesianSubspaceLDA
####

"""
    BayesianSubspaceLDA(; kwargs...)

$BayesianSubspaceLDA_DESCR

## Parameters

* `normalize=true`: Option to normalize the between class variance by number of
                    observations in each class, one of `true` or `false`.
* `out_dim`:        the output dimension, i.e dimension of the transformed
                    space, automatically set if 0 is given (default).
* `priors=nothing`: if `priors = nothing` estimates the prior probabilities
                    from the data, else it uses the user specified
                    `UnivariateFinite` prior probabilities.

See also the [package documentation](https://multivariatestatsjl.readthedocs.io/en/latest/lda.html).
For more information about the algorithm, see the paper by Howland & Park (2006), "Generalizing discriminant analysis using the generalized singular value decomposition", IEEE Trans. Patt. Anal. & Mach. Int., 26: 995-1006.
"""
@mlj_model mutable struct BayesianSubspaceLDA <: MMI.Probabilistic
    normalize::Bool=true
    out_dim::Int   = 0::(_ ≥ 0)
    priors::Union{Nothing,MLJBase.UnivariateFinite} = nothing
end

function MMI.fit(model::BayesianSubspaceLDA, ::Int, X, y)
    class_list = MMI.classes(y[1]) #class list containing unique entries in y
    nclasses   = length(class_list)

    # NOTE: copy/transpose
    Xm_t   = MMI.matrix(X, transpose=true) # now p x n matrix
    yplain = MMI.int(y) # vector of n ints in {1,..., nclasses}
    p, n   = size(Xm_t)

     #check to make sure we have more than one class
    nclasses >= 2 ||  throw(ArgumentError("The number of classes has to be greater than one"))

    # check output dimension default is min(p, nc-1)
    def_outdim = min(p, nclasses - 1)
    # if unset (0) use the default; otherwise try to use the provided one
    out_dim = ifelse(model.out_dim == 0, def_outdim, model.out_dim)
    # check if the given one is sensible
    out_dim ≤ def_outdim || throw(ArgumentError("`out_dim` must not be larger than `min(p, nc-1)` where `p` is the dimension of `X` and `nc` is the number of classes."))

    # check to make sure sample size is greater or equals nclasses
    n ≥ nclasses || throw(ArgumentError("The number of samples is less than the number of classes"))

    ## Estimates prior probabilities if unspecified by user.
    if model.priors == nothing
        priors = proportions(yplain)
    else
        #check if the length of priors is same as nclasses
        length(MMI.levels(model.priors)) == nclasses || throw(ArgumentError("Invalid size of `priors`"))
        priors = MMI.pdf.(model.priors, class_list)
    end

    # # Code gotten from JuliaStats/MultivariateStats.jl
    # # Compute centroids, class weights, and deviation from centroids
    # # Note Sb = Hb*Hb', Sw = Hw*Hw'
    cmeans, cweights, Hw = MS.center(Xm_t, Int.(yplain), nclasses) # Hw is p x n matrix , cmeans is p x nclasses matrix, cweights is nclasses x 1 vector
    dmeans = cmeans .- (model.normalize ? MS.mean(cmeans, dims=2) : cmeans * (cweights / (n))) #dmeans is p x nclasses
    Hb = model.normalize ? dmeans .* sqrt(n)  : dmeans .* sqrt.(cweights) #Hb is p x nclasses

    # Project to the subspace spanned by the within-class scatter
    # (essentially, PCA before LDA)
    Uw, Σw, _ = svd(Hw, full = false)
    keep = Σw .> sqrt(eps()) * maximum(Σw)
    projw = Uw[:, keep]
    pHb = projw' * Hb
    pHw = projw' * Hw

    # compute LDA projection G in the subspace spanned by the within-class variance.
    λ, G = MS.lda_gsvd(pHb, pHw, cweights) # λ is a (nclasses -1) x 1 vector containing the eigen values sorted in descending order.

    G = G[:, 1:out_dim]

    core_res = MS.SubspaceLDA(projw, G, λ, cmeans, cweights)

    ## scaled the overall projection_matrix (P = projw * G) by multiplying by sqrt(n - nclasses) this ensures Pᵀ*Σ*P=I
    ## where covariance estimate Σ = Sw / (n - nclasses)
    core_res.projLDA   .*= sqrt(n - nclasses)

    vproportions = λ ./ sum(λ) #proportions of variance
    cache     = nothing
    report    = (vproportions = vproportions,)
    fitresult = (core_res, class_list, priors)

    return fitresult, cache, report
end

function MMI.fitted_params(::BayesianSubspaceLDA, (core_res, class_list, priors))
    return (class_means       = MS.classmeans(core_res),
            projection_matrix = MS.projection(core_res),
            priors            = priors)
end

function MMI.predict(m::BayesianSubspaceLDA, (core_res, class_list, priors), Xnew)
    # projection of Xnew XWt is n x o  where o = number of out dims
    proj = core_res.projw * core_res.projLDA #proj is the projection_matrix
    XWt = MMI.matrix(Xnew) * proj
    # centroids in the transformed space, nc x o
    centroids = permutedims(proj' * core_res.cmeans)

    # compute the distances in the transformed space between pairs of rows
    # The discriminant matrix `P` is of dimension `n x nc`
    #  P[i,k] = -0.5*(xᵢ −  µₖ)ᵀΣ⁻¹(xᵢ −  µₖ) + log(priorsₖ) and Σ⁻¹ = I due to the nature of the projection_matrix (P = projw * G)
    P = pairwise(SqEuclidean(), XWt, centroids, dims=1)
    P .*= -0.5
    P .+= log.(priors)'

    # apply a softmax transformation to convert P to a probability matrix
    P .-= maximum(P, dims=2)
    P  .= exp.(P)
    P ./= sum(P, dims=2)

    return MMI.UnivariateFinite(class_list, P)
end

####
#### SubspaceLDA
####

"""
    SubspaceLDA(; kwargs...)

$SubspaceLDA_DESCR

## Parameters

* `normalize=true`:   Option to normalize the between class variance by number
                      of observations in each class, one of `true` or `false`.
* `out_dim`:          the output dimension, i.e dimension of the transformed
                      space, automatically set if 0 is given (default).
* `dist=SqEuclidean`: the distance metric to use when performing classification
                      (to compare the distance between a new point and
                      centroids in the transformed space), an alternative
                      choice can be the `CosineDist`.

See also the [package documentation](https://multivariatestatsjl.readthedocs.io/en/latest/lda.html).
For more information about the algorithm, see the paper by Howland & Park (2006), "Generalizing discriminant analysis using the generalized singular value decomposition", IEEE Trans. Patt. Anal. & Mach. Int., 26: 995-1006.
"""
@mlj_model mutable struct SubspaceLDA <: MMI.Probabilistic
    normalize::Bool=true
    out_dim::Int   = 0::(_ ≥ 0)
    dist::SemiMetric = SqEuclidean()
end

function MMI.fit(model::SubspaceLDA, ::Int, X, y)
    features   = Tables.schema(X).names
    class_list = MMI.classes(y[1]) #class list containing unique entries in y
    nclasses   = length(class_list)

    # NOTE: copy/transpose
    Xm_t   = MMI.matrix(X, transpose=true) # now p x n matrix
    yplain = MMI.int(y) # vector of n ints in {1,..., nclasses}
    p, n   = size(Xm_t)

     #check to make sure we have more than one class
    nclasses >= 2 ||  throw(ArgumentError("The number of classes has to be greater than one"))

    # check output dimension default is min(p, nc-1)
    def_outdim = min(p, nclasses - 1)
    # if unset (0) use the default; otherwise try to use the provided one
    out_dim = ifelse(model.out_dim == 0, def_outdim, model.out_dim)
    # check if the given one is sensible
    out_dim ≤ def_outdim || throw(ArgumentError("`out_dim` must not be larger than `min(p, nc-1)` where `p` is the dimension of `X` and `nc` is the number of classes."))

    # check to make sure sample size is greater or equals nclasses
    n ≥ nclasses || throw(ArgumentError("The number of samples is less than the number of classes"))

    # # Code gotten from JuliaStats/MultivariateStats.jl
    # # Compute centroids, class weights, and deviation from centroids
    # # Note Sb = Hb*Hb', Sw = Hw*Hw'
    cmeans, cweights, Hw = MS.center(Xm_t, Int.(yplain), nclasses) # Hw is p x n matrix , cmeans is p x nclasses matrix, cweights is nclasses x 1 vector
    dmeans = cmeans .- (model.normalize ? MS.mean(cmeans, dims=2) : cmeans * (cweights / (n))) #dmeans is p x nclasses
    Hb = model.normalize ? dmeans .* sqrt(n)  : dmeans .* sqrt.(cweights)

    # Project to the subspace spanned by the within-class scatter
    # (essentially, PCA before LDA)
    Uw, Σw, _ = svd(Hw, full = false)
    keep = Σw .> sqrt(eps()) * maximum(Σw)
    projw = Uw[:, keep]
    pHb = projw' * Hb
    pHw = projw' * Hw

    # compute LDA projection G in the subspace spanned by the within-class variance.
    λ, G = MS.lda_gsvd(pHb, pHw, cweights) # λ is a (nclasses -1) x 1 vector containing the eigen values sorted in descending order.

    G = G[:, 1:out_dim]

    core_res = MS.SubspaceLDA(projw, G, λ, cmeans, cweights)
    vproportions = λ ./ sum(λ) #proportions of variance

    cache     = nothing
    report    = (vproportions = vproportions,)
    fitresult = (core_res, class_list, features)

    return fitresult, cache, report
end

function MMI.fitted_params(::SubspaceLDA, (core_res, class_list))
    return (class_means       = MS.classmeans(core_res),
            projection_matrix = MS.projection(core_res))
end

function MMI.predict(m::SubspaceLDA, (core_res, class_list), Xnew)
    # projection of Xnew, XWt is n x o  where o = number of out dims
    proj = core_res.projw * core_res.projLDA #proj is the projection_matrix
    XWt = MMI.matrix(Xnew) * proj
    # centroids in the transformed space, nc x o
    centroids = permutedims(proj' * core_res.cmeans)

    # compute the distances in the transformed space between pairs of rows
    # the probability matrix is `n x nc` and normalised accross rows
    P = pairwise(m.dist, XWt, centroids, dims=1)
    # apply a softmax transformation
    P .-= maximum(P, dims=2)
    P  .= exp.(-P)
    P ./= sum(P, dims=2)

    return MMI.UnivariateFinite(class_list, P)
end

function MMI.transform(m::T, (core_res,), X) where T<:Union{LDA, SubspaceLDA, BayesianLDA, BayesianSubspaceLDA}
    # projection of X, XWt is n x o  where o = number of out dims
    proj = core_res.projw * core_res.projLDA #proj is the projection_matrix
    XWt = MMI.matrix(X) * proj
    return MMI.table(XWt, prototype = X)
end

####
#### METADATA
####

metadata_pkg.(
    (RidgeRegressor, PCA, KernelPCA, ICA, LDA,
     BayesianLDA, SubspaceLDA,  BayesianSubspaceLDA),
    name       = "MultivariateStats",
    uuid       = "6f286f6a-111f-5878-ab1e-185364afe411",
    url        = "https://github.com/JuliaStats/MultivariateStats.jl",
    license    = "MIT",
    julia      = true,
    is_wrapper = false)

metadata_model(RidgeRegressor,
    input   = Table(Continuous),
    target  = AbstractVector{Continuous},
    weights = false,
    descr   = RIDGE_DESCR)

metadata_model(PCA,
    input   = Table(Continuous),
    target  = Table(Continuous),
    weights = false,
    descr   = PCA_DESCR)

metadata_model(KernelPCA,
    input   = Table(Continuous),
    target  = Table(Continuous),
    weights = false,
    descr   = KPCA_DESCR)

metadata_model(ICA,
    input   = Table(Continuous),
    target  = Table(Continuous),
    weights = false,
    descr   = ICA_DESCR)

metadata_model(LDA,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    output  = Table(Continuous),
    descr   = LDA_DESCR)

metadata_model(BayesianLDA,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    output  = Table(Continuous),
    descr   = BayesianLDA_DESCR)

metadata_model(SubspaceLDA,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    output  = Table(Continuous),
    descr   = SubspaceLDA_DESCR)

metadata_model(BayesianSubspaceLDA,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    output  = Table(Continuous),
    descr   = BayesianSubspaceLDA_DESCR)

end

ext: jl


text:
using PyPlot, HDF5, MRIReco, LinearAlgebra

filename = "rawdata_heart_radial_55proj_34ch.h5"
data = permutedims(h5read(filename, "rawdata"),[3,2,1,4])
traj = permutedims(h5read(filename, "trajectory"),[3,2,1])
N = 240
Nc = 34

#############################################
# load data and form Acquisition data object
##############################################
tr = Trajectory(reshape(traj[1:2,:,:],2,:) ./ N, 55, 320, circular=false)
dat = Array{Array{Complex{Float64},2},3}(undef,1,1,1)
dat[1,1,1] = 1.e8.*reshape(data,:,Nc)
acqData = AcquisitionData(tr, dat, encodingSize=[N,N,1])

################################
# generate coil sensitivity maps
################################
@info "Espirit"
acqDataCart = regrid2d(acqData, (N,N); cgnr_iter=3)
sensitivity = espirit(acqDataCart,(6,6), 30, eigThresh_1=0.02, eigThresh_2=0.96)

##############################
# undersampled reconstructions
##############################
@info "undersampled reco"
rf = [1, 5/3, 5/2, 5]
img_cg = Vector{Array{ComplexF64,5}}(undef,4)

params = Dict{Symbol, Any}()
params[:reconSize] = (N,N)
params[:regularization] = "L2"
params[:solver] = "cgnr"
params[:senseMaps] = reshape(sensitivity, N, N, 1, Nc)
params[:iterations] = 20
params[:relTol] = 1.e-3
params[:λ] = 1.e-2

for (i,r) in enumerate(rf)
  @info "r=$r"
  # undersample profiles
  acqDataSub = convertUndersampledData(sample_kspace(acqData, Float64.(r), "regular",step=1))
  # SENSE reconstruction while monitoring error
  params[:reco] = "multiCoil"
  img_cg[i] = reconstruction(acqDataSub, params).data
end

##############
# plot figures
##############

figure("cardiac reconstructions",figsize=(6,1.9))
clf()
subplot(1,4,1); title("11 shots")
subplot(1,4,2); title("22 shots")
subplot(1,4,3); title("33 shots")
subplot(1,4,4); title("55 shots")
for (i,r) in enumerate(rf)

  subplot(1,4,5-i)
  imshow(abs.(img_cg[i][end:-1:1,end:-1:1,1,1,1]),cmap="gray")
  xticks([]);yticks([])
end
subplots_adjust(wspace=0.04,hspace=0.04,left=0.01,bottom=0.02,right=0.99,top=0.85)
savefig("Fig6.png",dpi=300)


text:
using TravelingSalesmanExact, GLPK, Test
using TravelingSalesmanExact: format_time
set_default_optimizer!(GLPK.Optimizer)

function test_valid_tour(tour, L)
    @test length(tour) == L
    @test isperm(tour)
    @test tour isa Vector{Int}
end

function tour_cost(tour, cost::AbstractMatrix)
    s = zero(eltype(cost))
    for (index, city) in enumerate(tour)
        nextcity = tour[mod1(index+1, end)]
        s += cost[city, nextcity]
    end
    return s
end

function tour_cost(
    tour,
    cities::AbstractVector;
    distance = TravelingSalesmanExact.euclidean_distance,
)
    n = length(cities)
    cost = [distance(cities[i], cities[j]) for i = 1:n, j = 1:n]
    return tour_cost(tour, cost)
end

function test_tour(input, opt = TravelingSalesmanExact.get_default_optimizer(); kwargs...)
    t, c = get_optimal_tour(input, opt; kwargs...)
    test_valid_tour(t, size(input, 1))
    if :distance in keys(kwargs)
        @test c ≈ tour_cost(t, input; distance = kwargs[:distance])
    else
        @test c ≈ tour_cost(t, input)
    end
    return t, c
end

@testset "`format_time`" begin
    @test format_time(0.052342) == "0.0523 seconds"
    @test format_time(0.452342) == "0.452 seconds"
    @test format_time(1000.452342) == "1000 seconds"
    @test format_time(1.452342) == "1.45 seconds"
end

@testset "Small random asymmetric" begin
    # chosen via rand(5,5)
    cost = 10 *
           [
        0.474886 0.350983 0.262651 0.138455 0.904042;
        0.683586 0.922968 0.278874 0.408406 0.0224372;
        0.513651 0.778167 0.140392 0.981211 0.891122;
        0.20529 0.976361 0.784706 0.98504 0.385203;
        0.489131 0.783738 0.538762 0.998821 0.0324331
    ]
    t1, c1 = test_tour(cost; verbose = true)
    t2, c2 = test_tour(cost; verbose = false)
    t3, c3 = test_tour(cost; verbose = false, lazy_constraints = true)
    @test c1 ≈ c2
    @test c2 ≈ c3

    # incorrect `symmetric` should give the wrong answers
    t4, c4 = get_optimal_tour(cost; verbose = false, symmetric = true)
    test_valid_tour(t4, 5)
    @test !(c2 ≈ c4)
    @test !(c4 ≈ tour_cost(t4, cost))
end

@testset "Medium random asymmetric" begin
    cost = 10 * rand(15, 15)
    t1, c1 = test_tour(cost; verbose = false)
    t2, c2 = test_tour(cost; verbose = true, lazy_constraints = true)
    @test c1 ≈ c2
end

@testset "Small random symmetric" begin
    cost = 10 *
           [
        0.474886 0.350983 0.262651 0.138455 0.904042;
        0.683586 0.922968 0.278874 0.408406 0.0224372;
        0.513651 0.778167 0.140392 0.981211 0.891122;
        0.20529 0.976361 0.784706 0.98504 0.385203;
        0.489131 0.783738 0.538762 0.998821 0.0324331
    ]
    cost_sym = cost + transpose(cost)
    t4, c4 = test_tour(cost_sym; verbose = false)
    t5, c5 = test_tour(cost_sym; symmetric = true, verbose = false)
    t6, c6 = test_tour(cost_sym; symmetric = false, verbose = false)
    t7, c7 = test_tour(cost_sym; symmetric = false, verbose = true, lazy_constraints = true)
    @test c4 ≈ c5
    @test c5 ≈ c6
    @test c6 ≈ c7
end

@testset "Small random cities" begin
    # cities = [ 100*rand(2) for _ in 1:5]
    cities = Array{Float64,1}[
        [48.8885, 41.0517],
        [35.6635, 12.1844],
        [95.6122, 15.9847],
        [67.5772, 9.54407],
        [16.6325, 51.9001],
    ]
    t7, c7 = test_tour(cities; verbose = true)
    t8, c8 = @inferred test_tour(cities; verbose = false, symmetric = false)
    cost = [TravelingSalesmanExact.euclidean_distance(c1, c2) for c1 in cities, c2 in cities]
    t9, c9 = test_tour(cost; verbose = false, symmetric = false)
    t10, c10 = test_tour(cost; verbose = false, symmetric = true)
    @test c7 ≈ c8
    @test c8 ≈ c9
    @test c9 ≈ c10
    @test c8 isa Float64
end

@testset "Medium random cities" begin
    # cities = [ 100*rand(2) for _ in 1:15]
    cities = Array{Float64,1}[
        [53.80996646037832, 19.19022341362766],
        [66.31428678667224, 61.80858818555463],
        [36.353148448445616, 41.255152847278474],
        [90.27061139767463, 6.978822544781482],
        [53.764677734646504, 29.81193167526588],
        [99.44803233011115, 66.55866205554817],
        [38.928236331570986, 18.22266760966855],
        [21.194033394576728, 56.14389699364914],
        [6.430489908372161, 92.08608470530129],
        [67.01950679275805, 53.82196541323463],
        [57.04078593130859, 88.74875672836768],
        [40.91801899608445, 23.47673907020007],
        [8.889246634653713, 44.183494522157886],
        [82.70468776751652, 74.73476674725042],
        [8.637366704185245, 87.23819026270408]
    ]
    t7, c7 = test_tour(cities; verbose = true, lazy_constraints=true)
    t8, c8 = @inferred test_tour(cities; verbose = false, symmetric = false)
    cost = [TravelingSalesmanExact.euclidean_distance(c1, c2) for c1 in cities, c2 in cities]
    t9, c9 = test_tour(cost; verbose = false, symmetric = false)
    t10, c10 = test_tour(cost; verbose = false, symmetric = true)
    @test c7 ≈ c8
    @test c8 ≈ c9
    @test c9 ≈ c10
    @test c8 isa Float64
end

@testset "Exceptions" begin
    cost = rand(5, 4)
    @test_throws ArgumentError test_tour(cost)

    cost = rand(5, 5)
    TravelingSalesmanExact.reset_default_optimizer!()
    @test_throws ArgumentError test_tour(cost)
    set_default_optimizer!(GLPK.Optimizer)
end

@testset "att48.tsp" begin
    cities = simple_parse_tsp(joinpath(@__DIR__, "..", "data", "att48.tsp"))

    @test cities == TravelingSalesmanExact.get_ATT48_cities()

    sym_tour, sym_cost = test_tour(
        cities;
        distance = TravelingSalesmanExact.ATT,
        verbose = true,
        slow = true,
    )
    @test sym_cost ≈ 10628

    asym_tour, asym_cost = test_tour(
        cities;
        distance = TravelingSalesmanExact.ATT,
        symmetric = false,
    )
    @test asym_cost ≈ 10628
end


text:
function DateBasedManeuverTriggers(arg0::AbsoluteDate, arg1::jdouble)
    return DateBasedManeuverTriggers((AbsoluteDate, jdouble), arg0, arg1)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_duration(obj::DateBasedManeuverTriggers)
    return jcall(obj, "getDuration", jdouble, ())
end

function get_end_date(obj::DateBasedManeuverTriggers)
    return jcall(obj, "getEndDate", AbsoluteDate, ())
end

function get_events_detectors(obj::DateBasedManeuverTriggers)
    return jcall(obj, "getEventsDetectors", Stream, ())
end

function get_field_events_detectors(obj::DateBasedManeuverTriggers, arg0::Field)
    return jcall(obj, "getFieldEventsDetectors", Stream, (Field,), arg0)
end

function get_name(obj::ManeuverTriggers)
    return jcall(obj, "getName", JString, ())
end

function get_parameters_drivers(obj::ManeuverTriggers)
    return jcall(obj, "getParametersDrivers", Vector{ParameterDriver}, ())
end

function get_start_date(obj::DateBasedManeuverTriggers)
    return jcall(obj, "getStartDate", AbsoluteDate, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function init(obj::DateBasedManeuverTriggers, arg0::SpacecraftState, arg1::AbsoluteDate)
    return jcall(obj, "init", void, (SpacecraftState, AbsoluteDate), arg0, arg1)
end

function is_firing(obj::DateBasedManeuverTriggers, arg0::AbsoluteDate)
    return jcall(obj, "isFiring", jboolean, (AbsoluteDate,), arg0)
end

function is_firing(obj::DateBasedManeuverTriggers, arg0::AbsoluteDate, arg1::Vector{jdouble})
    return jcall(obj, "isFiring", jboolean, (AbsoluteDate, Vector{jdouble}), arg0, arg1)
end

function is_firing(obj::DateBasedManeuverTriggers, arg0::FieldAbsoluteDate, arg1::Vector{RealFieldElement})
    return jcall(obj, "isFiring", jboolean, (FieldAbsoluteDate, Vector{RealFieldElement}), arg0, arg1)
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
# here we assume, that groups is a list of cnt of each value
# of the target variable
function information_gain(group, den)
    res = 0.0
    for x in group
        p = x/den
        res -= p * log2(p)
    end
    return res
end

function information_impurity(information_before, left, right, ll, lr, lt)
    information_before -
    ll/lt * information_gain(left, ll) -
    lr/lt * information_gain(right, lr)
end

function gini_index(group, den)
    den == 0 && return 0.0
    res = 0.0
    for x in group
        res += (x/den)^2
    end
    return res
end

function gini_impurity(gini_before, left, right, ll, lr, lt)
    ll/lt * gini_index(left, ll) +
    lr/lt * gini_index(right, lr) -
    gini_before
end


text:
module cfgrib

const cfgrib_jl_version = "0.0.0"

include("constants.jl")
include("cfmessage.jl")
include("indexing.jl")
include("dataset.jl")

end


text:
using LinearAlgebra
using CVChannel
using Test
using DelimitedFiles

"""
A qubit to qubit channel can be parameterized (up to a unitary)
by 6 parameters (t1,t2,t3,λ1,λ2,λ3). It is only non-unital if
it satisfies certain properties (Generalized Fujiwara-Algoet conditions).
These parameters may be shown all have absolute value bounded by 1.

In this script, we scan over the non-unital qubit to qubit channels
using the generalized FA conditions and check for multiplicativity for
such channels. You can hardcode the range you go over. You choose the
stepsize during the run. All data from the run is saved to a csv file
that you choose. The code tells you if non-multiplicativity was detected
regardless, in case you don't care about the specific data.
"""

#This function checks if the generalized Fujiwara-Algoet conditions hold.
#Assumes the input is of ordering (t1,t2,t3,λ1,λ2,λ3) and that it is not the
#case that t1=t2=t3=0 (i.e. it needs to be non-unital)
#See Theorem 4 of https://arxiv.org/abs/1306.0495
function genFAConditons(t_vec :: Vector, λ_vec :: Vector) :: Bool
    t = norm(t_vec); u = t_vec / t;

    q0 = (1 + λ_vec[1] + λ_vec[2] + λ_vec[3])/4;
    q1 = (1 + λ_vec[1] - λ_vec[2] - λ_vec[3])/4;
    q2 = (1 - λ_vec[1] + λ_vec[2] - λ_vec[3])/4;
    q3 = (1 - λ_vec[1] - λ_vec[2] - λ_vec[3])/4;
    q_vec = [q0,q1,q2,q3];
    q = 256 * prod(q_vec);                          #Eqn. (32)
    r = 1 - sum(λ_vec.^2) + 2*sum((λ_vec .* u).^2)  #Eqn. (31)

    if findfirst(x -> x < 0, q_vec) !== nothing
        return false
    elseif r^2 - q  < 0
        return false
    elseif t^2 - r + sqrt(r^2 - q) > 0              #Eqn. (30)
        return false
    else
        return true
    end
end

print("\nFirst we verify our function for generating the Choi matrix from the")
print(" 6 parameters does what it ought.")
#This function converts the parameters into the Choi matrix    #Eqn. (27)
#Assumes the input is of ordering (t1,t2,t3,λ1,λ2,λ3)
function tMatrixChoi(t :: Vector, λ :: Vector) :: Matrix{<:Number}
    choi_matrix = [
        0.5*(1+λ[3]+t[3])    0                   0.5*(t[1]+1im*t[2]) 0.5*(λ[1]+λ[2])
        0                    0.5*(1-λ[3]+t[3])   0.5*(λ[1]-λ[2])     0.5*(t[1]+1im*t[2])
        1/2*(t[1]-1im*t[2])  0.5*(λ[1]-λ[2])     0.5*(1-λ[3]-t[3])   0
        1/2*(λ[1]+λ[2])      0.5*(t[1]-1im*t[2]) 0                   1/2*(1+λ[3]-t[3])
    ]
    return choi_matrix
end

@testset "Verifying tMatrixChoi" begin
    t = [0,0,0]; λ=[0,0,0];
    target = [1/2 0 0 0 ; 0 1/2 0 0  ; 0 0 1/2 0 ; 0 0 0 1/2]
    @test tMatrixChoi(t,λ) == target

    t = [0,0,0]; λ=[1/3,1/4,1/2];
    target = [3/4 0 0 7/24 ; 0 1/4 1/24 0 ; 0 1/24 1/4 0 ; 7/24 0 0 3/4]
    @test isapprox(tMatrixChoi(t,λ),target, atol= 1e-9)

    t = [1/4,1/3,1/2]; λ=[0,0,0];
    target = [3/4 0 0.5*(1/4+1im*1/3) 0 ;
            0 3/4 0 0.5*(1/4+1im*1/3) ;
            0.5*(1/4-1im*1/3) 0 1/4 0 ;
            0 0.5*(1/4-1im*1/3) 0 1/4 ]
    @test isapprox(tMatrixChoi(t,λ),target, atol=1e-9)

    t = [1/4,1/3,1/2]; λ=[1/3,1/4,1/2];
    target = [ 1 0 1/8+1im/6 7/24;
               0 1/2 1/24 1/8+1im/6;
               1/8-1im/6 1/24 0 0;
               7/24 1/8-1im/6 0 1/2]
    @test isapprox(tMatrixChoi(t,λ),target,atol=1e-9)
end

#Here you set the range you want to run over
function stepSize() :: Float64
    print("\nNow one decides how to scan over the space.")
    print("\nIf you want to alter the step-size from non-uniform, you must hard-code it in the script. \n")
    print("\nPlease enter your choice of stepsize in (10^(-5),1): \n")
    acceptable_input = false
    while !acceptable_input
        step_size = parse(Float64,readline())
        if step_size ≥ 1 || step_size ≤ 1e-5
            print("\nPlease enter a number in (10^(-5),1): ")
        else
            return step_size
        end
    end
end

step_size = stepSize()

step_sizes = step_size*[1,1,1,1,1,1]
#|λ_i | ≤ 1 for all i ; See between (36) and (37)
#From this, it's easy to at least numerically see that |t_i| ≤ 1 ∀i
min_vals = [-1,-1,-1,-1,-1,-1]
max_vals = [1,1,1,1,1,1]

function verifyScanningRange(step_sizes :: Vector, min_vals :: Vector, max_vals :: Vector) :: Bool
    print("\n\n To be safe, we check how many points we'll need to evaluate.")
    total_num_of_scan_points = 1;
    for i in 1:6
        total_num_of_scan_points = total_num_of_scan_points * length(min_vals[i]:step_sizes[i]:max_vals[i])
    end

    print("\n Warning: if you check this range for evaluable_points, you will check ", total_num_of_scan_points," points. Do you wish to continue? (y/n) \n")
    proper_response = false
    while !proper_response
        user_response = readline()
        if user_response=="y"
            return true
        elseif user_response == "n"
            print("\n See you later!")
            return false
        else
            print("Please respond with `y` or `n`. \n ")
        end
    end
end

want_to_scan_range = verifyScanningRange(step_sizes, min_vals, max_vals)

if want_to_scan_range
    function numOfEvaluations(step_sizes :: Vector, min_vals :: Vector ,max_vals ::Vector) :: Matrix{Float16}
        #Here we check how many points you'll actually be running the SDP for
        #if you check this range. We also return what those parameters are
        nonUnitalParams = Float16[0 0 0 0 0 0]
        for it1 in min_vals[1]:step_sizes[1]:max_vals[1]
            for it2 in min_vals[2]:step_sizes[2]:max_vals[2]
                for it3 in min_vals[3]:step_sizes[3]:max_vals[3]
                    for it4 in min_vals[4]:step_sizes[4]:max_vals[4]
                        for it5 in min_vals[5]:step_sizes[5]:max_vals[5]
                            for it6 in min_vals[6]:step_sizes[6]:max_vals[6]
                                t = [it1,it2,it3];
                                λ = [it4,it5,it6]
                                if genFAConditons(t,λ) && !(it1==it2==it3==0)
                                    x = hcat(t,λ)
                                    newRow = copy(transpose(convert(Matrix{Float16},reshape(x,6,1))))
                                    nonUnitalParams = vcat(nonUnitalParams,newRow)
                                end
                            end
                        end
                    end
                end
            end
        end
        nonUnitalParams = nonUnitalParams[setdiff(1:end, 1), :]
        return nonUnitalParams
    end

    nonUnitalParams = numOfEvaluations(step_sizes, min_vals,max_vals)
else
    "\nSee you later!"
    exit()
end

function verifyProceeding(evaluable_points :: Int64) :: Bool
    print("Warning: if you run over this range, you will check ", evaluable_points," points. Do you wish to continue? (y/n) \n")
    proper_response = false
    while !proper_response
        user_response = readline()
        if user_response=="y"
            return true
        elseif user_response == "n"
            return false
        else
            print("Please respond with `y` or `n`. \n ")
        end
    end
end

evaluable_points = size(nonUnitalParams)[1]
want_to_proceed = verifyProceeding(evaluable_points)

#We write all the data to a csv, but we have the computer keep track
#of if non-multiplcativity has been detected
function scanRange(nonUnitalParams) :: Tuple{String,Bool}
    evaluable_points = size(nonUnitalParams)[1]
    diffMatrix = zeros(Float16,evaluable_points,1)
    nonMultDetected = false
    print("Please name the file you'd like to write to: \n")
    file_name = readline()
    file_to_open = string(file_name,".csv")
    evaluable_points = size(nonUnitalParams)[1]
    println("\nRunning the SDP solver for the first time. It will run faster after the first point.")
    for ctr = 1 : evaluable_points
        if (ctr % 5 == 0)
            println("\nNow evaluating evaluable point ", ctr," of ", evaluable_points,".")
        end

        t = nonUnitalParams[ctr,1:3]; λ = nonUnitalParams[ctr,4:6]
        orig_choi = tMatrixChoi(t,λ)
        test1 = pptCVDual(orig_choi,2,2)
        kron_par_choi = kron(orig_choi,orig_choi)
        par_choi = permuteSubsystems(kron_par_choi,[1,3,2,4],[2,2,2,2])
        test2 = pptCVPrimal(par_choi,4,4)
        diffMatrix[ctr] = test2[1] - test1[1]^2
        if abs(diffMatrix[ctr]) > 1e-5 && diffMatrix[ctr] > 0
            if eigmin(orig_choi) >= 0 #This should always be true, but I don't trust numerical error, so it's a final sanity check
                nonMultDetected = true
            end
        end
    end

    data_to_save = hcat(nonUnitalParams,diffMatrix)
    writedlm(file_to_open, data_to_save, ',')
    return file_to_open, nonMultDetected
end

if want_to_proceed
    @testset "Multiplicative Over Range" begin
    (file_to_open,non_mult_detected)=scanRange(nonUnitalParams)
    print("\nIf it was multiplicative over the range, the test passes.")
    @test !non_mult_detected
    end
    print("\nAll of the data from the run was saved to your file, if you'd like to look at it.")
else
    print("\nSee you later!")
end

max_forks_repo_name: ChitambarLab/cv-channel


text:
#                                                             convert A.U. -> Å
ase_cell(lattice) = pyimport("ase").cell.Cell(Array(lattice)' / austrip(1u"Å"))
ase_cell(model::Model) = ase_cell(model.lattice)

function ase_atoms(lattice_or_model, atoms, positions, magnetic_moments=[])
    if isempty(magnetic_moments)  # Collect Z magnetic moments
        magmoms = nothing
    else
        magmoms = [normalize_magnetic_moment(mom)[3] for mom in magnetic_moments]
    end
    cell    = ase_cell(lattice_or_model)
    symbols = string.(atomic_symbol.(atoms))
    scaled_positions = reduce(hcat, positions)'
    pyimport("ase").Atoms(;symbols, cell, pbc=true, scaled_positions, magmoms)
end
function ase_atoms(model::Model, magnetic_moments=[])
    ase_atoms(model.lattice, model.atoms, model.positions, magnetic_moments)
end

function load_lattice_ase(pyobj::PyObject)
    ase = pyimport("ase")
    if pyisinstance(pyobj, ase.Atoms)
        if !all(pyobj.pbc)
            error("DFTK only supports calculations with periodic boundary conditions.")
        end
        load_lattice_ase(pyobj.cell)
    elseif pyisinstance(pyobj, ase.cell.Cell)
        lattice = zeros(3, 3)
        cell_julia = convert(Array, pyobj)  # Array of arrays
        for i = 1:3, j = 1:3
            lattice[i, j] = austrip(cell_julia[j][i] * u"Å")
        end
        Mat3(lattice)
    else
        error("load_lattice_ase not implemented for python type $pyobj")
    end
end

function load_atoms_ase(pyobj::PyObject)
    @assert pyisinstance(pyobj, pyimport("ase").Atoms)
    # TODO Be smarter and look at the calculator to determine the psps
    [ElementCoulomb(number) for number in pyobj.get_atomic_numbers()]
end

function load_positions_ase(pyobj::PyObject)
    @assert pyisinstance(pyobj, pyimport("ase").Atoms)
    [Vec3(pos) for pos in eachrow(pyobj.get_scaled_positions())]
end

function load_magnetic_moments_ase(pyobj::PyObject)
    @assert pyisinstance(pyobj, pyimport("ase").Atoms)
    [normalize_magnetic_moment(magmom)
     for magmom in pyobj.get_initial_magnetic_moments()]
end


text:
include("../tools/play.jl")

big_fish = BigFish()

width = RES_W * ENLARGEMENT
height = RES_H * ENLARGEMENT
number_iterations = 5000

# Set the random seeds
rand_seed = 0
game_level_seed_gen = RandGen()
seed(game_level_seed_gen, rand_seed)
seed(big_fish.basic_game.game.level_seed_rand_gen, randint(game_level_seed_gen))
big_fish.basic_game.game.level_seed_low = 0
big_fish.basic_game.game.level_seed_high = typemax(UInt32)

# Options

# Game will only use assets from a single theme
#big_fish.basic_game.game.options.restrict_themes = true

# Game will use monochromatic rectangles instead of human designed assets. Best used with option restrict_themes set to true.
#big_fish.basic_game.game.options.use_monochrome_assets = true

# Game will use pure black background
#big_fish.basic_game.game.options.use_backgrounds = false

# Reaching the end of a level does not end the episode, and the seed for the new level is derived from the current level seed
#big_fish.basic_game.game.options.use_sequential_levels = true

# Game will be easier
#big_fish.basic_game.game.options.distribution_mode = EasyMode

# Only affects some procgen games, doesn't affect bigfish
# Observations will be centered on the agent and not display the full level
#big_fish.basic_game.game.options.center_agent = true

game_init(big_fish)
reset(big_fish)

first_iteration = true

@play width height number_iterations begin

    if first_iteration
        observe(big_fish)
        global first_iteration = false
    else
        # The Framerate is 15 Hz
        sleep(1/15)

        set_action(big_fish, current_key)
        step(big_fish)
    end
end

println("Finished")


text:
# example:
#   df=wdi("NY.GNP.PCAP.CD", ["US","BR"], 1980, 2012, true)
function wdi(indicators::Union(ASCIIString,Array{ASCIIString,1}),
             countries::Union(ASCIIString,Array{ASCIIString,1}),
             startyear::Integer=1800,endyear::Integer=3000;
             extra::Bool = false,
             format::Type = AbstractTimedata)
    ## download multiple indicators / countries

    if countries == "all"
        countries = all_countries
    end

    if typeof(countries) == ASCIIString
        countries = [countries]
    end

    for c in countries
        if ! (c in all_countries)
            error("country ",c," not found")
        end
    end

    if ! (startyear < endyear)
        error("startyear has to be < endyear. startyear=",startyear,". endyear=",endyear)
    end

    df = DataFrame()

    if typeof(indicators) == ASCIIString
        indicators=[indicators]
    end

    firstInd = true
    for ind in indicators
        dfn = wdi_download(ind, countries, startyear, endyear)
        if firstInd
            df = vcat(df, dfn)
            firstInd = false
        else
            nams = names(df)
            nNams = length(nams)
            ## if nNams > 4
                ## allNamesExceptIndicator = nams[[[1:3], [5:end]]]
            ## else
            allNamesExceptIndicator = nams[1:3]
            ## end
            df = join(df, dfn, on = allNamesExceptIndicator,
                      kind = :outer)
        end
    end

    if extra
        cntdat = getWBMeta("countries")
        df = join(df, cntdat, on = :iso2c)
    end

    ## sort data for plotting
    sort!(df, cols = [:iso2c, :time])

    if format == AbstractTimedata

        dats = df[:time]
        delete!(df, :time)
        return Timedata(df, dats)
    elseif format == DataFrame
        return df
    end
end

function wdi_download(indicator::ASCIIString,
                      country::Union(ASCIIString,Array{ASCIIString,1}),
                      startyear::Integer,
                      endyear::Integer)
    # download single indicator for possibly multiple countries and
    # years
    if typeof(country) == ASCIIString # single country
        countryString = country
    else
        countryString = join(country, ";")
    end

    url = string("http://api.worldbank.org/countries/", countryString, "/indicators/", indicator,
                 "?date=", startyear,":", endyear, "&per_page=25000", "&format=json")
    json = [download_parse_json(url)[2]]

    parse_wdi(indicator,json, startyear, endyear)
end

function download_parse_json(url::ASCIIString)
    ## try to download data and process JSON to Dict
    println("download: ",url)
    request = HTTPC.get(url)
    if request.http_code != 200
        error("download failed")
    end
    JSON.parse(bytestring(request.body))
end

function parse_wdi(indicator::ASCIIString, json, startyear::Integer,
                   endyear::Integer)
    country_id = ASCIIString[] # no NA
    country_name = UTF8String[] # no NA
    value = DataArray(Float64, 0)
    ## value = ASCIIString[]
    date = ASCIIString[] # no NA

    for d in json
        append!(country_id, [d["country"]["id"]]) # iso-2
        append!(country_name, [d["country"]["value"]]) # name
        push!(value, tofloat(d["value"]))
        append!(date, [d["date"]])
    end

    date = Date[formatDate(dat) for dat in date]

    df = DataFrame(time = date, iso2c = country_id,
                   country = country_name)
    df[symbol(indicator)] = value

    # filter missing/wrong data
    complete_cases!(df)

    startY = formatDate(string(startyear))
    endY = formatDate(string(endyear))
    checkyear(x) = (x >= startY) & (x <= endY)
    yind = map(checkyear,df[:time])
    yind = convert(DataArray{Bool, 1}, yind)
    df[yind, :]
end

all_countries = ["AW","AF","A9","AO","AL","AD","L5","1A","AE","AR","AM","AS","AG","AU","AT","AZ","BI","BE","BJ","BF","BD","BG","BH","BS","BA","BY","BZ","BM","BO","BR","BB","BN","BT","BW","C9","CF","CA","C4","B8","C5","CH","JG","CL","CN","CI","C6","C7","CM","CD","CG","CO","KM","CV","CR","C8","S3","CU","CW","KY","CY","CZ","DE","DJ","DM","DK","DO","DZ","4E","Z4","7E","Z7","EC","EG","XC","ER","ES","EE","ET","EU","F1","FI","FJ","FR","FO","FM","GA","GB","GE","GH","GN","GM","GW","GQ","GR","GD","GL","GT","GU","GY","XD","HK","HN","XE","HR","HT","HU","ID","IM","IN","XY","IE","IR","IQ","IS","IL","IT","JM","JO","JP","KZ","KE","KG","KH","KI","KN","KR","KV","KW","XJ","LA","LB","LR","LY","LC","ZJ","L4","XL","XM","LI","LK","XN","XO","LS","LT","LU","LV","MO","MF","MA","L6","MC","MD","MG","MV","ZQ","MX","MH","XP","MK","ML","MT","MM","XQ","ME","MN","MP","MZ","MR","MU","MW","MY","XU","M2","NA","NC","NE","NG","NI","NL","XR","NO","NP","NZ","XS","OE","OM","S4","PK","PA","PE","PH","PW","PG","PL","PR","KP","PT","PY","PS","S2","PF","QA","RO","RU","RW","8S","SA","L7","SD","SN","SG","SB","SL","SV","SM","SO","RS","ZF","SS","ZG","S1","ST","SR","SK","SI","SE","SZ","SX","A4","SC","SY","TC","TD","TG","TH","TJ","TM","TL","TO","TT","TN","TR","TV","TZ","UG","UA","XT","UY","US","UZ","VC","VE","VI","VN","VU","1W","WS","A5","YE","ZA","ZM","ZW"]

## all_countries = ["AW", "AF", "A9", "AO", "AL", "AD", "1A", "AE", "AR", "AM", "AS", "AG", "AU", "AT", "AZ", "BI", "BE", "BJ", "BF", "BD", "BG", "BH", "BS", "BA", "BY", "BZ", "BM", "BO", "BR", "BB", "BN", "BT", "BW", "C9", "CF", "CA", "C4", "C5", "CH", "JG", "CL", "CN", "CI", "C6", "C7", "CM", "CD", "CG", "CO", "KM", "CV", "CR", "C8", "S3", "CU", "CW", "KY", "CY", "CZ", "DE", "DJ", "DM", "DK", "DO", "DZ", "4E", "Z4", "7E", "Z7", "EC", "EG", "XC", "ER", "ES", "EE", "ET", "EU", "FI", "FJ", "FR", "FO", "FM", "GA", "GB", "GE", "GH", "GN", "GM", "GW", "GQ", "GR", "GD", "GL", "GT", "GU", "GY", "XD", "HK", "HN", "XE", "HR", "HT", "HU", "ID", "IM", "IN", "XY", "IE", "IR", "IQ", "IS", "IL", "IT", "JM", "JO", "JP", "KZ", "KE", "KG", "KH", "KI", "KN", "KR", "KV", "KW", "XJ", "LA", "LB", "LR", "LY", "LC", "ZJ", "XL", "XM", "LI", "LK", "XN", "XO", "LS", "LT", "LU", "LV", "MO", "MF", "MA", "MC", "MD", "MG", "MV", "ZQ", "MX", "MH", "XP", "MK", "ML", "MT", "MM", "XQ", "ME", "MN", "MP", "MZ", "MR", "MU", "MW", "MY", "XU", "M2", "NA", "NC", "NE", "NG", "NI", "NL", "XR", "NO", "NP", "NZ", "XS", "OE", "OM", "S4", "PK", "PA", "PE", "PH", "PW", "PG", "PL", "PR", "KP", "PT", "PY", "PS", "S2", "PF", "QA", "RO", "RU", "RW", "8S", "SA", "SD", "SN", "SG", "SB", "SL", "SV", "SM", "SO", "RS", "ZF", "SS", "ZG", "S1", "ST", "SR", "SK", "SI", "SE", "SZ", "SX", "A4", "SC", "SY", "TC", "TD", "TG", "TH", "TJ", "TM", "TL", "TO", "TT", "TN", "TR", "TV", "TZ", "UG", "UA", "XT", "UY", "US", "UZ", "VC", "VE", "VI", "VN", "VU", "1W", "WS", "A5", "YE", "ZA", "ZM", "ZW" ]


text:
mutable struct BSONIter
    _wrap_::Vector{UInt8}
    done::Bool

    function BSONIter(bsonObject::BSONObject)
        bsonIter = new(alloc_buf(128), false)
        ccall((:bson_iter_init, libbson), Bool, (BytePtr, BytePtr),
              bsonIter._wrap_, bsonObject._wrap_) || error("BSONIter(): failure")
        bsonIter.done = !ccall((:bson_iter_next, libbson), Bool, (BytePtr, ), bsonIter._wrap_)
        bsonIter
    end

    function BSONIter(bsonObject::BSONObject, key::AbstractString)
        bsonIter = new(alloc_buf(128), false)
        ccall((:bson_iter_init, libbson), Bool, (BytePtr, BytePtr),
              bsonIter._wrap_, bsonObject._wrap_) || error("BSONIter(): failure")
        keyCStr = string(key)
        bsonIter.done = !ccall((:bson_iter_find, libbson), Bool, (BytePtr, BytePtr),
                               bsonIter._wrap_, keyCStr)
        bsonIter
    end

    function BSONIter(bsonArray::BSONArray)
        bsonIter = new(alloc_buf(128), false)
        ccall((:bson_iter_init, libbson), Bool, (BytePtr, BytePtr),
              bsonIter._wrap_, bsonArray._wrap_) || error("BSONIter(): failure")
        bsonIter.done = !ccall((:bson_iter_next, libbson), Bool, (BytePtr, ), bsonIter._wrap_)
        bsonIter
    end

    function BSONIter(bsonArray::BSONArray, key::Integer)
        bsonIter = new(alloc_buf(128), false)
        ccall((:bson_iter_init, libbson), Bool, (BytePtr, BytePtr),
              bsonIter._wrap_, bsonArray._wrap_) || error("BSONIter(): failure")
        keyCStr = string(key)
        bsonIter.done = !ccall((:bson_iter_find, libbson), Bool, (BytePtr, BytePtr),
                               bsonIter._wrap_, keyCStr)
        bsonIter
    end
end
export BSONIter

# Index

function Base.getindex(bsonObject::BSONObject, key::AbstractString)
    bsonIter = BSONIter(bsonObject, key)
    bsonIter.done && error("key not found: $(repr(key))")
    value(bsonIter)
end

function Base.getindex(bsonArray::BSONArray, key::Integer)
    bsonIter = BSONIter(bsonArray, key - 1)
    bsonIter.done && error("key not found: $(repr(key))")
    value(bsonIter)
end

# Associative collection

function Base.get!(bsonObject::BSONObject, key::AbstractString, default)
    iter = BSONIter(bsonObject, key)
    iter.done || return value(iter)
    append(bsonObject, key, default)
    default
end

# Iterator

@static if VERSION < v"0.7-"
Base.start(bsonObject::BSONObject) = BSONIter(bsonObject)

Base.next(bsonObject::BSONObject, bsonIter::BSONIter) =
    ((key(bsonIter), value(bsonIter)), next!(bsonIter))

Base.done(bsonObject::BSONObject, bsonIter::BSONIter) = bsonIter.done

Base.start(bsonArray::BSONArray) = BSONIter(bsonArray)

Base.next(bsonArray::BSONArray, bsonIter::BSONIter) = (value(bsonIter), next!(bsonIter))

Base.done(bsonArray::BSONArray, bsonIter::BSONIter) = bsonIter.done
else
Base.iterate(bsonObject::BSONObject) = BSONIter(bsonObject)
Base.iterate(bsonArray::BSONArray)   = BSONIter(bsonArray)

Base.iterate(bsonObject::BSONObject, iter::BSONIter) = ((key(iter), value(iter)), next!(iter))
Base.iterate(bsonArray::BSONArray,   iter::BSONIter) = (value(iter), next!(iter))
end

# Private

function key(bsonIter::BSONIter)
    bsonIter.done && error("already done iteration")
    cstr = ccall((:bson_iter_key, libbson), BytePtr, (BytePtr, ), bsonIter._wrap_)
    cstr == C_NULL && error("bson_iter_key: failure")
    String(unsafe_string(cstr))
end

function value_type(bsonIter::BSONIter)
    bsonIter.done && error("already done iteration")
    ccall((:bson_iter_type, libbson), BSONType, (BytePtr, ), bsonIter._wrap_)
end

function value(bsonIter::BSONIter)
    ty = value_type(bsonIter)
    if ty == BSON_TYPE_DOUBLE
        ccall((:bson_iter_double, libbson), Cdouble, (BytePtr, ), bsonIter._wrap_)
    elseif ty == BSON_TYPE_INT32
        ccall((:bson_iter_int32, libbson), Int32, (BytePtr, ), bsonIter._wrap_)
    elseif ty == BSON_TYPE_INT64
        ccall((:bson_iter_int64, libbson), Int64, (BytePtr, ), bsonIter._wrap_)
    elseif ty == BSON_TYPE_BOOL
        ccall((:bson_iter_bool, libbson), Bool, (BytePtr, ), bsonIter._wrap_)
    elseif ty == BSON_TYPE_DATE_TIME
        ts = ccall((:bson_iter_date_time, libbson), Int64, (BytePtr, ), bsonIter._wrap_)
        Dates.unix2datetime(ts / 1000)
    elseif ty == BSON_TYPE_NULL
        nothing
    elseif ty == BSON_TYPE_MINKEY
        :minkey
    elseif ty == BSON_TYPE_MAXKEY
        :maxkey
    elseif ty == BSON_TYPE_UTF8
        String(unsafe_string(ccall((:bson_iter_utf8, libbson), BytePtr,
                                   (BytePtr, BytePtr), bsonIter._wrap_, C_NULL)))
    elseif ty == BSON_TYPE_OID
        data = alloc_buf(12)
        ptr = ccall((:bson_iter_oid, libbson), BytePtr, (BytePtr,), bsonIter._wrap_)
        unsafe_copyto!(pointer(data), ptr, 12)
        BSONOID(data)
    elseif ty == BSON_TYPE_DOCUMENT
        rlen = Ref{UInt32}(0)
        data = Ref{BytePtr}(C_NULL)
        ccall((:bson_iter_document, libbson), VoidPtr,
              (BytePtr, Ref{UInt32}, Ref{BytePtr}),
              bsonIter._wrap_, rlen, data)
        BSONObject(data[], rlen[], bsonIter)
    elseif ty == BSON_TYPE_ARRAY
        rlen  = Ref{UInt32}(0)
        data = Ref{BytePtr}(C_NULL)
        ccall((:bson_iter_array, libbson), VoidPtr,
              (BytePtr, Ref{UInt32}, Ref{BytePtr}),
              bsonIter._wrap_, rlen, data)
        BSONArray(data[], rlen[], bsonIter)
    elseif ty == BSON_TYPE_BINARY
        rlen = Ref{UInt32}(0)
        data = Ref{BytePtr}(C_NULL)
        ccall((:bson_iter_binary, libbson), VoidPtr,
              (BytePtr, VoidPtr, Ref{UInt32}, Ref{BytePtr}),
              bsonIter._wrap_, C_NULL, rlen, data)
        len = Int(rlen[])
        dataArray = alloc_buf(len)
        unsafe_copyto!(pointer(dataArray), data[], len)
        dataArray
    else
        error("unhandled BSONType $ty")
    end
end
export value

function next!(bsonIter::BSONIter)
    bsonIter.done = !ccall((:bson_iter_next, libbson), Bool, (BytePtr, ), bsonIter._wrap_)
    bsonIter
end


text:
const __bodyfunction__ = Dict{Method,Any}()

# Find keyword "body functions" (the function that contains the body
# as written by the developer, called after all missing keyword-arguments
# have been assigned values), in a manner that doesn't depend on
# gensymmed names.
# `mnokw` is the method that gets called when you invoke it without
# supplying any keywords.
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{MinusEq{typeof(^)},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}},GVar{Float64, Float64}})   # time: 0.088883124
    Base.precompile(Tuple{PlusEq{typeof(^)},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}},GVar{Float64, Float64}})   # time: 0.08728426
    Base.precompile(Tuple{MinusEq{typeof(/)},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}}})   # time: 0.07497767
    Base.precompile(Tuple{PlusEq{typeof(exp)},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}}})   # time: 0.06419293
    Base.precompile(Tuple{PlusEq{typeof(/)},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}}})   # time: 0.041938376
    Base.precompile(Tuple{MinusEq{typeof(exp)},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}}})   # time: 0.04189861
    Base.precompile(Tuple{typeof(ROT),Float64,Float64,Irrational{:π}})   # time: 0.03481159
    Base.precompile(Tuple{MinusEq{typeof(exp)},ComplexF64,ComplexF64})   # time: 0.034116935
    Base.precompile(Tuple{PlusEq{typeof(exp)},ComplexF64,ComplexF64})   # time: 0.034112282
    Base.precompile(Tuple{typeof(PUSH!),Tuple{Float64, Float64, Vector{Float64}}})   # time: 0.031773232
    Base.precompile(Tuple{typeof(auto_expand),Expr})   # time: 0.031103289
    Base.precompile(Tuple{MinusEq{typeof(log)},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}}})   # time: 0.02943123
    Base.precompile(Tuple{PlusEq{typeof(log)},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}}})   # time: 0.028707368
    Base.precompile(Tuple{PlusEq{typeof(log)},ComplexF64,ComplexF64})   # time: 0.02753613
    Base.precompile(Tuple{MinusEq{typeof(inv)},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}}})   # time: 0.025810735
    Base.precompile(Tuple{PlusEq{typeof(inv)},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}}})   # time: 0.0257753
    Base.precompile(Tuple{PlusEq{typeof(^)},ComplexF64,ComplexF64,Float64})   # time: 0.019021787
    Base.precompile(Tuple{typeof(POP!),Vector{Any}})   # time: 0.01667686
    Base.precompile(Tuple{PlusEq{typeof(/)},ComplexF64,ComplexF64,ComplexF64})   # time: 0.016225988
    Base.precompile(Tuple{typeof(COPYPOP!),Vector{Any}})   # time: 0.014254036
    Base.precompile(Tuple{MinusEq{typeof(/)},ComplexF64,ComplexF64,ComplexF64})   # time: 0.014187833
    Base.precompile(Tuple{typeof(IROT),Float64,Float64,Irrational{:π}})   # time: 0.013347673
    Base.precompile(Tuple{typeof(auto_alloc),Expr})   # time: 0.011528118
    Base.precompile(Tuple{MinusEq{typeof(log)},ComplexF64,ComplexF64})   # time: 0.008786706
    Base.precompile(Tuple{typeof(COPYPUSH!),Tuple{Float64, Float64, Vector{Float64}}})   # time: 0.007922663
    Base.precompile(Tuple{PlusEq{typeof(inv)},ComplexF64,ComplexF64})   # time: 0.00758488
    Base.precompile(Tuple{typeof(NEG),Complex{GVar{Float64, Float64}}})   # time: 0.00560381
    Base.precompile(Tuple{typeof(POP!),Vector{Any},Vector{Float64}})   # time: 0.005197541
    Base.precompile(Tuple{PlusEq{typeof(*)},ComplexF64,ComplexF64,ComplexF64})   # time: 0.004626822
    Base.precompile(Tuple{MinusEq{typeof(^)},ComplexF64,ComplexF64,Float64})   # time: 0.004009681
    Base.precompile(Tuple{typeof(NEG),ComplexF64})   # time: 0.002459327
    Base.precompile(Tuple{MinusEq{typeof(*)},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}}})   # time: 0.002440985
    Base.precompile(Tuple{typeof(PUSH!),Vector{Any}})   # time: 0.002278265
    Base.precompile(Tuple{typeof(PUSH!),Vector{Any},Vector{Int64}})   # time: 0.002102671
    Base.precompile(Tuple{PlusEq{typeof(+)},ComplexF64,ComplexF64,ComplexF64})   # time: 0.001917077
    Base.precompile(Tuple{PlusEq{typeof(-)},ComplexF64,ComplexF64,ComplexF64})   # time: 0.001910042
    Base.precompile(Tuple{MinusEq{typeof(*)},ComplexF64,ComplexF64,ComplexF64})   # time: 0.001726165
    Base.precompile(Tuple{Type{BennettLog}})   # time: 0.001334291
    Base.precompile(Tuple{PlusEq{typeof(identity)},ComplexF64,ComplexF64})   # time: 0.001304064
    Base.precompile(Tuple{PlusEq{typeof(-)},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}}})   # time: 0.001289984
    Base.precompile(Tuple{MinusEq{typeof(+)},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}}})   # time: 0.001287298
    Base.precompile(Tuple{PlusEq{typeof(+)},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}}})   # time: 0.0012845
    Base.precompile(Tuple{MinusEq{typeof(-)},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}},Complex{GVar{Float64, Float64}}})   # time: 0.001269163
    Base.precompile(Tuple{typeof(-),NoGrad{Float64}})   # time: 0.001207526
    Base.precompile(Tuple{MinusEq{typeof(+)},ComplexF64,ComplexF64,ComplexF64})   # time: 0.001151932
    Base.precompile(Tuple{MinusEq{typeof(inv)},ComplexF64,ComplexF64})   # time: 0.001151206
    Base.precompile(Tuple{MinusEq{typeof(-)},ComplexF64,ComplexF64,ComplexF64})   # time: 0.001148769
    Base.precompile(Tuple{typeof(chfield),ComplexF64,typeof(real),Float64})   # time: 0.001119287
    Base.precompile(Tuple{typeof(POP!),Tuple{Float64, GVar{Float64, Float64}, Vector{GVar{Float64, Float64}}}})   # time: 0.001055845
    Base.precompile(Tuple{PlusEq{typeof(complex)},Complex{GVar{Float64, Float64}},GVar{Float64, Float64},GVar{Float64, Float64}})   # time: 0.001005807
    Base.precompile(Tuple{PlusEq{typeof(complex)},ComplexF64,Float64,Float64})   # time: 0.001005075
    Base.precompile(Tuple{typeof(chfield),ComplexF64,typeof(imag),Float64})   # time: 0.001002244
end


text:
function inverted_vector(vector_in::Array{<:Real, 1})
    num_points = size(vector_in)[1]
    return reshape([vector_in[num_points + 1 - i]
                    for i in 1:num_points], (num_points, ))
end # inverted_vector


text:
@testset "example kernel run" begin
    #make system
    perbead_param_keys= (:k,:kθ)
    global_param_keys= (:d₀,)
    function chain_pe(p,b,n,perbead_params,global_params)
        rp= (p .- b)
        rn= (n .- b)
        dp= sqrt(sum(rp.^2))
        dn= sqrt(sum(rn.^2))
        cosθ= sum(rp .* rn)/dp/dn
        perbead_params.k*1//2*((dn-global_params.d₀)^2) + perbead_params.kθ*(cosθ+1)
    end
    beaddef1= BeadDefinition(chain_pe, perbead_param_keys, global_param_keys, Val(3))
    globalparams= (d₀=1.0, )
    system= ChainSystem(0.0,beaddef1,globalparams)
    #create some chain
    params= fill((k=3.5, kθ=4.4),10)
    params[10]= (k=0.0, kθ= 0.0)# turn off interation to split ends
    params[1]= (k=3.5, kθ= 0.0)
    pos= zeros(10,3)
    pos[:,3] .= range(0.0, step= 1.0, length=10)
    chain1= Chain(beaddef1,params, globalparams, pos)
    system= append(system,chain1)
    #create chain2 translated [1 0 0] from chain1
    chain2= Chain(beaddef1,params, globalparams, pos .+ [1 0 0])
    system= append(system,chain2)
    #specific forces
    params= (k=10.0,)
    pe(x,time,params)= 1//2 * params.k * sum((x[1]).^2)
    interactions = [((1,1),)
                    ((1,2),)]
    fdef= SpecificForce(pe, params, interactions, Val(3))
    system= append(system,fdef)

    params2= (k₁=10.0, k₂=2.0)
    pe2(x,time,params)= 1//2 * sum((params.k₁, params.k₂, params.k₂*time) .* ((x[1] .- x[2]).^2))
    interactions2 = [((1,1),(2,1))
                    ((1,2),(1,1))]
    fdef2= SpecificForce(pe2, params2, interactions2, Val(3))
    system= append(system,fdef2)

    simpos= copy(system.init_pos)
    simvel= copy(system.init_vel)
    invmass= 1.0
    γ= 0.0
    invβ= 1.0
    VNbeads= Val(size(simpos)[1])
    VNdims= Val(size(simpos)[2])
    step0= UInt64(0)
    time0= 0.0
    Vrngkey= Val(UInt64(123))
    Δt= 0.01
    VChainbounds= Val(system.chainbounds)
    VNsteps= Val(2)

    refcpukernel!(simpos,
                simvel,
                system.sp_forces,
                system.sp_interactions,
                system.sp_params,
                system.beaddef.chain_force,
                copy(system.perbead_params),
                system.global_params,
                invmass,
                γ,
                invβ,
                VNbeads,
                step0,
                time0,
                Vrngkey,
                Δt,
                VNdims,
                VChainbounds,
                VNsteps)
end


text:
julia> 0b1011010111 == 0o1327 == 0x2d7 == 727
true

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text: @show (asin ∘ sin)(0.5)


text:
module Dyno

    export import_modfile, import_data, solve, simulate, deterministic, Der, Model

    using Dolang
    using DataFrames
    using DataStructures
    using Distributions
    using Dolang
    import JSON
    import Dolang
    import Dolang: Der

    include("sanitize.jl")

    include("trisystem.jl")
    include("modfile_parsing.jl")
    include("model_import.jl")
    include("perturbation.jl")
    include("deterministic.jl")

end # module


text:
immutable RibbonGeometry <: Gadfly.GeometryElement
    default_statistic::Gadfly.StatisticElement

    function RibbonGeometry(default_statistic=Gadfly.Stat.identity())
        new(default_statistic)
    end
end

const ribbon = RibbonGeometry

function default_statistic(geom::RibbonGeometry)
    return geom.default_statistic
end

function element_aesthetics(::RibbonGeometry)
    return [:x, :ymin, :ymax, :color]
end

function render(geom::RibbonGeometry, theme::Gadfly.Theme,
                aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.ribbon", aes, :x, :ymin, :ymax)
    Gadfly.assert_aesthetics_equal_length("Geom.ribbon", aes,
                                          element_aesthetics(geom)...)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = PooledDataArray(ColorValue[theme.default_color])
    aes = inherit(aes, default_aes)

    if length(aes.color) == 1 &&
        !(isa(aes.color, PooledDataArray) && length(levels(aes.color)) > 1)
        max_points = collect((Any, Any), zip(aes.x, aes.ymax))
        sort!(max_points, by=first)

        min_points = collect((Any, Any), zip(aes.x, aes.ymin))
        sort!(min_points, by=first, rev=true)

        form = compose(polygon(chain(min_points, max_points)...),
                       fill(theme.lowlight_color(aes.color[1])),
                       svgclass("geometry"))
    else
        max_points = Dict{ColorValue, Array{(Any, Any), 1}}()
        for (x, y, c) in zip(aes.x, aes.ymax, aes.color)
            if !haskey(max_points, c)
                max_points[c] = Array((Any, Any),0)
            end
            push!(max_points[c], (x, y))
        end

        min_points = Dict{ColorValue, Array{(Any, Any), 1}}()
        for (x, y, c) in zip(aes.x, aes.ymin, aes.color)
            if !haskey(min_points, c)
                min_points[c] = Array((Any, Any),0)
            end
            push!(min_points[c], (x, y))
        end

        forms = Array(Any, length(max_points))
        for (i, c) in enumerate(keys(max_points))
            c_max_points = max_points[c]
            c_min_points = min_points[c]
            sort!(c_max_points, by=first)
            sort!(c_min_points, by=first, rev=true)

            forms[i] =
                compose(polygon(chain(c_min_points, c_max_points)...),
                        fill(theme.lowlight_color(c)),
                        opacity(theme.lowlight_opacity),
                        svgclass(@sprintf("geometry color_%s",
                                          escape_id(aes.color_label([c])[1]))))
        end
        form = combine(forms...)
    end

    compose(form, stroke(nothing), linewidth(theme.line_width))
end


text:
@testset for llog in (Libm.log_tang,)
    for T in (Float32, Float64)
        @test isnan(llog(T(NaN)))
        @test llog(T(Inf)) == Inf
        @test llog(T(0)) == -Inf
        s = logspace(T(-10),T(10),1000)
        @test Base.log.(s) ≈ llog.(s)
        s = linspace(T(0),realmin(T),1000)
        @test Base.log.(s) ≈ llog.(s)
        @test_throws DomainError llog(T(-1))

        @test Base.log(prevfloat(T(Inf))) == llog(prevfloat(T(Inf)))
        @test Base.log(nextfloat(T(0))) == llog(nextfloat(T(0)))
    end
end

@testset for llog1p in (Libm.log1p_tang,)
    for T in (Float32, Float64)
        @test isnan(llog1p(T(NaN)))
        @test llog1p(T(Inf)) == Inf
        @test llog1p(T(0)) == 0
        @test llog1p(T(-1)) == -Inf
        s = logspace(T(-10),T(10),1000)
        @test Base.log1p.(s) ≈ llog1p.(s)
        s = -logspace(T(-10),T(0),1000)
        @test Base.log1p.(s) ≈ llog1p.(s)
        @test_throws DomainError llog1p(T(-2))
        @test Base.log1p(prevfloat(T(Inf))) == llog1p(prevfloat(T(Inf)))
        @test Base.log1p(nextfloat(T(-1))) == llog1p(nextfloat(T(-1)))
    end
end


text:
SensitivityAlg(args...;kwargs...) = @error("The SensitivtyAlg choice mechanism was completely overhauled. Please consult the local sensitivity documentation for more information")

abstract type AbstractForwardSensitivityAlgorithm{CS,AD,FDT} <: DiffEqBase.AbstractSensitivityAlgorithm{CS,AD,FDT} end
abstract type AbstractAdjointSensitivityAlgorithm{CS,AD,FDT} <: DiffEqBase.AbstractSensitivityAlgorithm{CS,AD,FDT} end
abstract type AbstractSecondOrderSensitivityAlgorithm{CS,AD,FDT} <: DiffEqBase.AbstractSensitivityAlgorithm{CS,AD,FDT} end

struct ForwardSensitivity{CS,AD,FDT} <: AbstractForwardSensitivityAlgorithm{CS,AD,FDT}
  autojacvec::Bool
end
Base.@pure function ForwardSensitivity(;
                                       chunk_size=0,autodiff=true,
                                       diff_type=Val{:central},
                                       autojacvec=autodiff)
  ForwardSensitivity{chunk_size,autodiff,diff_type}(autojacvec)
end

struct ForwardDiffSensitivity{CS,CTS} <: AbstractForwardSensitivityAlgorithm{CS,Nothing,Nothing}
end
Base.@pure function ForwardDiffSensitivity(;chunk_size=0,convert_tspan=nothing)
  ForwardDiffSensitivity{chunk_size,convert_tspan}()
end

"""
ODE:
 Rackauckas, C. and Ma, Y. and Martensen, J. and Warner, C. and Zubov, K. and Supekar,
 R. and Skinner, D. and Ramadhana, A. and Edelman, A., Universal Differential Equations
 for Scientific Machine Learning,	arXiv:2001.04385

 Hindmarsh, A. C. and Brown, P. N. and Grant, K. E. and Lee, S. L. and Serban, R.
 and Shumaker, D. E. and Woodward, C. S., SUNDIALS: Suite of nonlinear and
 differential/algebraic equation solvers, ACM Transactions on Mathematical
 Software (TOMS), 31, pp:363–396 (2005)

 Chen, R.T.Q. and Rubanova, Y. and Bettencourt, J. and Duvenaud, D. K.,
 Neural ordinary differential equations. In Advances in neural information processing
 systems, pp. 6571–6583 (2018)

 Pontryagin, L. S. and Mishchenko, E.F. and Boltyanskii, V.G. and Gamkrelidze, R.V.
 The mathematical theory of optimal processes. Routledge, (1962)

 Rackauckas, C. and Ma, Y. and Dixit, V. and Guo, X. and Innes, M. and Revels, J.
 and Nyberg, J. and Ivaturi, V., A comparison of automatic differentiation and
 continuous sensitivity analysis for derivatives of differential equation solutions,
 arXiv:1812.01892

DAE:
 Cao, Y. and Li, S. and Petzold, L. and Serban, R., Adjoint sensitivity analysis
 for differential-algebraic equations: The adjoint DAE system and its numerical
 solution, SIAM journal on scientific computing 24 pp: 1076-1089 (2003)

SDE:
 Gobet, E. and Munos, R., Sensitivity Analysis Using Ito-Malliavin Calculus and
 Martingales, and Application to Stochastic Optimal Control,
 SIAM Journal on control and optimization, 43, pp. 1676-1713 (2005)

 Li, X. and Wong, T.-K. L.and Chen, R. T. Q. and Duvenaud, D.,
 Scalable Gradients for Stochastic Differential Equations,
 PMLR 108, pp. 3870-3882 (2020), http://proceedings.mlr.press/v108/li20i.html
"""
struct BacksolveAdjoint{CS,AD,FDT,VJP,NOISE} <: AbstractAdjointSensitivityAlgorithm{CS,AD,FDT}
  autojacvec::VJP
  checkpointing::Bool
  noise::NOISE
  noisemixing::Bool
end
Base.@pure function BacksolveAdjoint(;chunk_size=0,autodiff=true,
                                      diff_type=Val{:central},
                                      autojacvec=autodiff,
                                      checkpointing=true, noise=true,noisemixing=false)
  BacksolveAdjoint{chunk_size,autodiff,diff_type,typeof(autojacvec),typeof(noise)}(autojacvec,checkpointing,noise,noisemixing)
end

"""
 Rackauckas, C. and Ma, Y. and Martensen, J. and Warner, C. and Zubov, K. and Supekar,
 R. and Skinner, D. and Ramadhana, A. and Edelman, A., Universal Differential Equations
 for Scientific Machine Learning,	arXiv:2001.04385

 Hindmarsh, A. C. and Brown, P. N. and Grant, K. E. and Lee, S. L. and Serban, R.
 and Shumaker, D. E. and Woodward, C. S., SUNDIALS: Suite of nonlinear and
 differential/algebraic equation solvers, ACM Transactions on Mathematical
 Software (TOMS), 31, pp:363–396 (2005)

 Rackauckas, C. and Ma, Y. and Dixit, V. and Guo, X. and Innes, M. and Revels, J.
 and Nyberg, J. and Ivaturi, V., A comparison of automatic differentiation and
 continuous sensitivity analysis for derivatives of differential equation solutions,
 arXiv:1812.01892
"""
struct InterpolatingAdjoint{CS,AD,FDT,VJP,NOISE} <: AbstractAdjointSensitivityAlgorithm{CS,AD,FDT}
  autojacvec::VJP
  checkpointing::Bool
  noise::NOISE
  noisemixing::Bool
end
Base.@pure function InterpolatingAdjoint(;chunk_size=0,autodiff=true,
                                         diff_type=Val{:central},
                                         autojacvec=autodiff,
                                         checkpointing=false, noise=true,noisemixing=false)
  InterpolatingAdjoint{chunk_size,autodiff,diff_type,typeof(autojacvec),typeof(noise)}(autojacvec,checkpointing,noise,noisemixing)
end

struct QuadratureAdjoint{CS,AD,FDT,VJP} <: AbstractAdjointSensitivityAlgorithm{CS,AD,FDT}
  autojacvec::VJP
  abstol::Float64
  reltol::Float64
  compile::Bool
end
Base.@pure function QuadratureAdjoint(;chunk_size=0,autodiff=true,
                                         diff_type=Val{:central},
                                         autojacvec=autodiff,abstol=1e-6,
                                         reltol=1e-3,compile=false)
  QuadratureAdjoint{chunk_size,autodiff,diff_type,typeof(autojacvec)}(autojacvec,abstol,reltol,compile)
end

struct TrackerAdjoint <: AbstractAdjointSensitivityAlgorithm{nothing,true,nothing} end
struct ReverseDiffAdjoint <: AbstractAdjointSensitivityAlgorithm{nothing,true,nothing} end
struct ZygoteAdjoint <: AbstractAdjointSensitivityAlgorithm{nothing,true,nothing} end

"""
Wang, Q., Hu, R., and Blonigan, P. Least squares shadowing sensitivity analysis of
chaotic limit cycle oscillations. Journal of Computational Physics, 267, 210-224 (2014).
"""
struct ForwardLSS{CS,AD,FDT,aType} <: AbstractForwardSensitivityAlgorithm{CS,AD,FDT}
  alpha::aType # alpha: weight of the time dilation term in LSS.
end
Base.@pure function ForwardLSS(;
                                chunk_size=0,autodiff=true,
                                diff_type=Val{:central},
                                alpha=CosWindowing())
  ForwardLSS{chunk_size,autodiff,diff_type,typeof(alpha)}(alpha)
end

"""
Wang, Q., Hu, R., and Blonigan, P. Least squares shadowing sensitivity analysis of
chaotic limit cycle oscillations. Journal of Computational Physics, 267, 210-224 (2014).
"""
struct AdjointLSS{CS,AD,FDT,aType} <: AbstractAdjointSensitivityAlgorithm{CS,AD,FDT}
  alpha::aType # alpha: weight of the time dilation term in LSS.
end
Base.@pure function AdjointLSS(;
                                chunk_size=0,autodiff=true,
                                diff_type=Val{:central},
                                alpha=10.0)
  AdjointLSS{chunk_size,autodiff,diff_type,typeof(alpha)}(alpha)
end

abstract type WindowingChoice end
struct CosWindowing <: WindowingChoice end
struct Cos2Windowing <: WindowingChoice end

"""
Ni, A., and Wang, Q. Sensitivity analysis on chaotic dynamical systems by Non-Intrusive
Least Squares Shadowing (NILSS). Journal of Computational Physics 347, 56-77 (2017).
"""
struct NILSS{CS,AD,FDT,RNG} <: AbstractAdjointSensitivityAlgorithm{CS,AD,FDT}
  rng::RNG
  nseg::Int
  nstep::Int
  autojacvec::Bool
end
Base.@pure function NILSS(nseg, nstep; rng = Xorshifts.Xoroshiro128Plus(rand(UInt64)),
                                chunk_size=0,autodiff=true,
                                diff_type=Val{:central},
                                autojacvec = autodiff
                                )
  NILSS{chunk_size,autodiff,diff_type,typeof(rng)}(rng, nseg, nstep, autojacvec)
end

"""
 Johnson, S. G., Notes on Adjoint Methods for 18.336, Online at
 http://math.mit.edu/stevenj/18.336/adjoint.pdf (2007)
"""
struct SteadyStateAdjoint{CS,AD,FDT,VJP,LS} <: AbstractAdjointSensitivityAlgorithm{CS,AD,FDT}
  autojacvec::VJP
  linsolve::LS
end
Base.@pure function SteadyStateAdjoint(;chunk_size=0,autodiff=true,diff_type=Val{:central},autojacvec=autodiff,linsolve=DEFAULT_LINSOLVE)
  SteadyStateAdjoint{chunk_size,autodiff,diff_type,typeof(autojacvec),typeof(linsolve)}(autojacvec,linsolve)
end

abstract type VJPChoice end
struct ZygoteVJP <: VJPChoice end
struct EnzymeVJP <: VJPChoice end
struct TrackerVJP <: VJPChoice end
struct ReverseDiffVJP{compile} <: VJPChoice
  ReverseDiffVJP(compile=false) = new{compile}()
end

abstract type NoiseChoice end
struct ZygoteNoise <: NoiseChoice end
struct ReverseDiffNoise{compile} <: NoiseChoice
  ReverseDiffNoise(compile=false) = new{compile}()
end

@inline convert_tspan(::ForwardDiffSensitivity{CS,CTS}) where {CS,CTS} = CTS
@inline convert_tspan(::Any) = nothing
@inline alg_autodiff(alg::DiffEqBase.AbstractSensitivityAlgorithm{CS,AD,FDT}) where {CS,AD,FDT} = AD
@inline get_chunksize(alg::DiffEqBase.AbstractSensitivityAlgorithm{CS,AD,FDT}) where {CS,AD,FDT} = CS
@inline diff_type(alg::DiffEqBase.AbstractSensitivityAlgorithm{CS,AD,FDT}) where {CS,AD,FDT} = FDT
@inline function get_jacvec(alg::DiffEqBase.AbstractSensitivityAlgorithm)
  alg.autojacvec isa Bool ? alg.autojacvec : true
end
@inline ischeckpointing(alg::DiffEqBase.AbstractSensitivityAlgorithm, sol=nothing) = false
@inline ischeckpointing(alg::InterpolatingAdjoint) = alg.checkpointing
@inline ischeckpointing(alg::InterpolatingAdjoint, sol) = alg.checkpointing || !sol.dense
@inline ischeckpointing(alg::BacksolveAdjoint, sol=nothing) = alg.checkpointing

@inline isnoise(alg::DiffEqBase.AbstractSensitivityAlgorithm) = false
@inline isnoise(alg::InterpolatingAdjoint) = alg.noise
@inline isnoise(alg::BacksolveAdjoint) = alg.noise

@inline isnoisemixing(alg::DiffEqBase.AbstractSensitivityAlgorithm) = false
@inline isnoisemixing(alg::InterpolatingAdjoint) = alg.noisemixing
@inline isnoisemixing(alg::BacksolveAdjoint) = alg.noisemixing

@inline compile_tape(vjp::ReverseDiffVJP{compile}) where compile = compile
@inline compile_tape(noise::ReverseDiffNoise{compile}) where compile = compile
@inline compile_tape(autojacvec::Bool) = false
@inline compile_tape(sensealg::QuadratureAdjoint) = sensealg.compile

struct ForwardDiffOverAdjoint{A} <: AbstractSecondOrderSensitivityAlgorithm{nothing,true,nothing}
  adjalg::A
end


text:
module DiscordAPI

using Discord

include("./ZoomAPI.jl")

const PREFIX = "/zoomcord "

function create_cmd(client::Client, msg::Message)
    response = create_meeting("RCazfv6rSo6eBdoO4PXnWA")
    reply(client, msg, String(response))
end

function client(token::String, prefix::String)
    client = Client(token; prefix=prefix, presence=(game = (name = prefix, type = AT_LISTENING),),)
    add_command!(
        client, :create_meeting, create_cmd;
        pattern=r"^create$",
        help="Create a new meeting"
    )
    return client
end

function start()
    c = client(ENV["DISCORD_TOKEN"], PREFIX)
    while true
        open(c)
        wait(c)
        isopen(c) && close(c)
    end
end

end

max_stars_repo_name: MLH-Fellowship/ZoomCord

max_issues_repo_name: MLH-Fellowship/ZoomCord

max_forks_repo_name: MLH-Fellowship/ZoomCord