text:
"""
    eposvd(D; nlv)
Pre-processing X-data by external parameter orthogonalization (EPO) (Roger et al 2003).
* `D` : Data (m, p) containing detrimental information.
* `nlv` : Nb. of first loadings vectors of D considered for the orthogonalization.

The objective is to remove from a dataset X (n, p) some "detrimental"
information (e.g. humidity patterns in signals) defined by a dataset `D` (m, p).
The method orthogonalizes the observations (rows of X) to the
detrimental sub-space defined by the first `nlv` loadings vectors
computed from a (non-centered) PCA of D.

Function `eposvd` makes a SVD factorization of D and returns
two matrices:
* M (p, p) : The orthogonalization matrix that can be used to correct the X-data.
* P (p, `nlv`) : The matrix of the loading vectors of D.

The data corrected from the detrimental information D
can be computed by X_corrected = X * M.

# References

Roger, J.-M., Chauchard, F., Bellon-Maurel, V., 2003. EPO-PLS external parameter
orthogonalisation of PLS application to temperature-independent measurement
of sugar content of intact fruits.
Chemometrics and Intelligent Laboratory Systems 66, 191-204.
https://doi.org/10.1016/S0169-7439(03)00051-0

Roger, J.-M., Boulet, J.-C., 2018. A review of orthogonal projections for calibration.
Journal of Chemometrics 32, e3045. https://doi.org/10.1002/cem.3045

# Example

```julia
n = 4 ; p = 8 ;
X = rand(n, p)
m = 3 ;
D = rand(m, p)    # Detrimental information

nlv = 2
res = eposvd(D; nlv = nlv)
res.M      # orthogonalization matrix
res.P      # detrimental directions (columns of P = loadings of D)
```

The matrix corrected from D can be computed by:
```julia
X_corr = X * res.M
```
Rows of the corrected matrix X_corr
are orthogonal to the loadings vectors (columns of P):
```julia
X_corr * res.P
```
"""
function eposvd(D; nlv)
    D = ensure_mat(D)
    m, p = size(D)
    nlv = min(nlv, m, p)
    Id = Diagonal(I, p)
    if nlv == 0
        M = Id
        P = nothing
    else
        P = svd(D).V[:, 1:nlv]
        M = Id - P * P'
    end
    (M = M, P = P)
end

"""
    detrend(X)
Linear de-trend transformation of each row of X-data.
* `X` : X-data.

The function fits a univariate linear regression to each observation
and returns the residuals.
"""
function detrend(X)
    M = copy(X)
    detrend!(M)
    M
end

function detrend!(X)
    n, p = size(X)
    xmean = mean(1:p)
    x = collect(1:p)
    xc = x .- xmean
    xtx = dot(xc, xc)
    yc = similar(xc)
    @inbounds for i = 1:n
        vy = vrow(X, i)
        ymean = mean(vy)
        yc .= vy .- ymean
        b = dot(xc, yc) / xtx
        int = ymean - xmean * b
        X[i, :] .= vy .- (int .+ x * b)
    end
end

"""
    fdif(X; f = 2)
    fdif!(M, X; f = 2)
Finite differences of each row of a matrix X.
* `X` : X-data (n, p).
* `M` : Pre-allocated output matrix (n, p - f + 1).
* `f` : Size of the window (nb. points involved) for the finite differences.
    The range of the window (= nb. intervals of two successive colums) is f - 1.

The finite differences can be used for computing discrete derivates.
The method reduces the column-dimension: (n, p) --> (n, p - f + 1).

The in-place function stores the output in `M`.
"""
function fdif(X; f = 2)
    n, p = size(X)
    M = similar(X, n, p - f + 1)
    fdif!(M, X; f)
    M
end

function fdif!(M, X; f = 2)
    p = size(X, 2)
    zp = p - f + 1
    @inbounds for j = 1:zp
        M[:, j] .= vcol(X, j + f - 1) .- vcol(X, j)
    end
end

"""
    interpl(X, wl; wlfin,
        fun = CubicSpline)
Sampling signals by interpolation.
* `X` : Matrix (n, p) of signals (rows).
* `wl` : Values representing the column "names" of `X`.
    Must be a numeric vector of length p, or an AbstractRange.
* `wlfin` : Final values where to interpolate within the range of `wl`.
    Must be a numeric vector, or an AbstractRange.
* `fun` : Function defining the interpolation method.

The function uses package DataInterpolations.jl.

Possible values of `fun` (methods) are:
- `linear_int`: A linear interpolation (LinearInterpolation).
- `quadratic_int`: A quadratic interpolation (QuadraticInterpolation).
- `quadratic_spline`: A quadratic spline interpolation(QuadraticSpline).
- `cubic_spline`: A cubic spline interpolation (CubicSpline)

## References

Package Interpolations.jl
https://github.com/PumasAI/DataInterpolations.jl
https://htmlpreview.github.io/?https://github.com/PumasAI/DataInterpolations.jl/blob/v2.0.0/example/DataInterpolations.html
"""
function interpl(X, wl; wlfin,
        fun = cubic_spline)
    X = ensure_mat(X)
    n = size(X, 1)
    q = length(wlfin)
    zX = similar(X, n, q)
    @inbounds for i = 1:n
        itp = fun(vrow(X, i), wl)
        zX[i, :] .= itp.(wlfin)
    end
    zX
end
# Due to conflicts with Interpolations.jl
linear_int(y, x) = DataInterpolations.LinearInterpolation(y, x)
quadratic_int(y, x) = DataInterpolations.QuadraticInterpolation(y, x)
quadratic_spline(y, x) = DataInterpolations.QuadraticSpline(y, x)
cubic_spline(y, x) = DataInterpolations.CubicSpline(y, x)

"""
    interpl_mon(X, wl; wlfin,
        fun = FritschCarlsonMonotonicInterpolation)
Sampling signals by monotonic interpolation.
* `X` : Matrix (n, p) of signals (rows).
* `wl` : Values representing the column "names" of `X`.
    Must be a numeric vector of length p, or an AbstractRange.
* `wlfin` : Values where to interpolate within the range of `wl`.
    Must be a numeric vector, or an AbstractRange.
* `fun` : Function defining the interpolation method.

See e.g. https://en.wikipedia.org/wiki/Monotone_cubic_interpolation.

The function uses package Interpolations.jl.

Possible values of `fun` (methods) are:
- `LinearMonotonicInterpolation`
- `FiniteDifferenceMonotonicInterpolation` : Classic cubic
- `CardinalMonotonicInterpolation`
- `FritschCarlsonMonotonicInterpolation`
- `FritschButlandMonotonicInterpolation`
- `SteffenMonotonicInterpolation`
See https://github.com/JuliaMath/Interpolations.jl/pull/243/files#diff-92e3f2a374c9a54769084bad1bbfb4ff20ee50716accf008074cda7af1cd6149

## References

Package Interpolations.jl
https://github.com/JuliaMath/Interpolations.jl

Fritsch & Carlson (1980), "Monotone Piecewise Cubic Interpolation",
doi:10.1137/0717021.

Fritsch & Butland (1984), "A Method for Constructing Local Monotone Piecewise
Cubic Interpolants", doi:10.1137/0905021.

Steffen (1990), "A Simple Method for Monotonic Interpolation
in One Dimension", http://adsabs.harvard.edu/abs/1990A%26A...239..443S
"""
function interpl_mon(X, wl; wlfin,
        fun = FritschCarlsonMonotonicInterpolation)
    X = ensure_mat(X)
    n = size(X, 1)
    q = length(wlfin)
    zX = similar(X, n, q)
    @inbounds for i = 1:n
        itp = interpolate(wl, vrow(X, i), fun())
        zX[i, :] .= itp.(wlfin)
    end
    zX
end

"""
    mavg(X; f)
Moving averages smoothing of each row of X-data.
* `X` : X-data.
* `f` : Size (nb. points involved) of the filter.

The smoothing is computed by convolution (with padding), with function
imfilter of package ImageFiltering.jl. The centered kernel is ones(`f`) / `f`.
Each returned point is located on the center of the kernel.
"""
function mavg(X; f)
    M = copy(X)
    mavg!(M; f)
    M
end

function mavg!(X; f)
    n, p = size(X)
    kern = ImageFiltering.centered(ones(f) / f) ;
    out = similar(X, p)
    @inbounds for i = 1:n
        imfilter!(out, vrow(X, i), kern)
        X[i, :] .= out
    end
end

"""
    mavg_runmean(X, f)
    mavg_runmean!(M, X, f)
Moving average smoothing of each row of a matrix X.
* `X` : X-data (n, p).
* `M` : Pre-allocated output matrix (n, p - f + 1).
* `f` : Size (nb. points involved) of the filter.

The smoothing is computed by convolution, without padding (which reduces
the column dimension). The function is an adaptation/simplification of function
runmean (V. G. Gumennyy) of package Indicators.jl. See
https://github.com/dysonance/Indicators.jl/blob/a449c1d68487c3a8fea0008f7abb3e068552aa08/src/run.jl.
The kernel is ones(`f`) / `f`. Each returned point is located on the 1st unit of the kernel.
In general, this function can be faster than mavg, especialy for in-place versions.

The in-place function stores the output in `M`.
"""
function mavg_runmean(X; f)
    n, p = size(X)
    M = similar(X, n, p - f + 1)
    mavg_runmean!(M, X; f)
    M
end

function mavg_runmean!(M, X; f)
    n, zp = size(M)
    out = similar(M, zp)
    @inbounds for i = 1:n
        Jchemo.runmean!(out, vrow(X, i); f)
        M[i, :] .= out
    end
end

## adaptation from function runmean  (V. G. Gumennyy)
## of package Indicators.jl
function runmean!(out, x; f)
    ## x : (n,)
    ## out : (n - f + 1,)
    ## f : integer
    ## There is no padding
    ## Location on the 1st unit of the kernel
    n = length(x)
    zsum = 0.
    @inbounds for i = 1:f
        zsum += x[i]
    end
    out[1] = zsum / f
    @inbounds for i = (f + 1):n
        zsum += x[i] - x[i - f]
        out[i - f + 1] = zsum / f
    end
end

"""
    savgk(m, pol, d)
Compute the kernel of the Savitzky-Golay filter.
* `m` : Nb. points of the half window (m >= 1)
    --> the size of the kernel is odd (f = 2 * m + 1):
    x[-m], x[-m+1], ..., x[0], ...., x[m-1], x[m].
* `pol` : Polynom order (1 <= pol <= 2 * m).
    The case "pol = 0" (simple moving average) is not allowed by the funtion.
* `d` : Derivation order (0 <= d <= pol).
    If `d = 0`, there is no derivation (only polynomial smoothing).

## References

Luo, J., Ying, K., Bai, J., 2005. Savitzky–Golay smoothing and differentiation
filter for even number data. Signal Processing 85, 1429–1434.
https://doi.org/10.1016/j.sigpro.2005.02.002
"""
function savgk(m, pol, d)
    @assert m >= 1 "m must be >= 1"
    @assert pol >= 1 && pol <= 2 * m "pol must agree with: 1 <= pol <= 2 * m"
    @assert 0 <= d && d <= pol "d must agree with: 0 <= d <= pol"
    f = 2 * m + 1
    S = zeros(Int64(f), Int64(pol) + 1) ;
    u = collect(-m:m)
    @inbounds for j in 0:pol
        S[:, j + 1] .= u.^j
    end
    G = S * inv(S' * S)
    kern = factorial(d) * vcol(G, d + 1)
    (S = S, G = G, kern = kern)
end

"""
    savgol(X; f, pol, d)
Savitzky-Golay smoothing of each row of a matrix `X`.
* `X` : X-data (n, p).
* `f` : Size of the filter (nb. points involved in the kernel). Must be odd and >= 3.
    The half-window size is m = (f - 1) / 2.
* `pol` : Polynom order (1 <= pol <= f - 1).
* `d` : Derivation order (0 <= d <= pol).

The smoothing is computed by convolution (with padding), with function
imfilter of package ImageFiltering.jl. Each returned point is located on the center
of the kernel. The kernel is computed with function `savgk`.
"""
function savgol(X; f, pol, d)
    M = copy(X)
    savgol!(M; f, pol, d)
    M
end

function savgol!(X; f, pol, d)
    @assert isodd(f) && f >= 3 "f must be odd and >= 3"
    n, p = size(X)
    m = (f - 1) / 2
    kern = savgk(m, pol, d).kern
    zkern = ImageFiltering.centered(kern)
    out = similar(X, p)
    @inbounds for i = 1:n
        ## convolution with "replicate" padding
        imfilter!(out, vrow(X, i), reflect(zkern))
        X[i, :] .= out
    end
end

"""
    snv(X)
Standard-normal-variate (SNV) transformation of each row of X-data.
* `X` : X-data.

Each row of `X` is centered and scaled.
"""
function snv(X)
    M = copy(X)
    snv!(M)
    M
end

function snv!(X)
    p = size(X, 2)
    mu = vec(Statistics.mean(X; dims = 2))
    s = vec(Statistics.std(X; corrected = false, dims = 2))
    @inbounds for j = 1:p
        X[:, j] .= (vcol(X, j) .- mu) ./ s
    end
end


text:
using Distributions
using ProgressMeter
using Base.LinAlg.BLAS
using Compat
using Devectorize
using HDF5
using PyCall
@pyimport matplotlib.pyplot as plt
@pyimport numpy as np

import StatsBase.fit

function calculate_weight_gradient!(rbm::RBM, h_pos::Mat{Float64}, v_pos::Mat{Float64}, h_neg::Mat{Float64}, v_neg::Mat{Float64}, lr::Float64; approx="CD")
    ## Load step buffer with negative-phase
    gemm!('N', 'T', lr, h_neg, v_neg, 0.0, rbm.dW)          # dW <- LearRate*<h_neg,v_neg>
    ## Stubtract step buffer from positive-phase to get gradient
    gemm!('N', 'T', lr, h_pos, v_pos, -1.0, rbm.dW)         # dW <- LearnRate*<h_pos,v_pos> - dW
    ## Second-Order EMF Correction (for EMF-TAP2, EMF-TAP3)
    if contains(approx,"tap")
        buf2 = gemm('N', 'T', h_neg-abs2(h_neg), v_neg-abs2(v_neg)) .* rbm.W
        axpy!(-lr, buf2, rbm.dW)
    end
    ## Third-Order EMF Correction (for EMF-TAP3)
    if approx == "tap3"
        buf3 = gemm('N','T', (h_neg-abs2(h_neg)) .* (0.5-h_neg), (v_neg-abs2(v_neg)) .* (0.5-v_neg)) .* rbm.W2
        axpy!(-2.0*lr, buf3, rbm.dW)
    end
    ## Apply Momentum (adding last gradient to this one)
    axpy!(rbm.momentum, rbm.dW_prev, rbm.dW)    # rbm.dW += rbm.momentum * rbm.dW_prev
end

function update_weights!(rbm::RBM,approx::AbstractString)
    axpy!(1.0,rbm.dW,rbm.W)             # Take step: W = W + dW
    copy!(rbm.dW_prev, rbm.dW)          # Save the current step for future use
    # if contains(approx,"tap")
    rbm.W2 = rbm.W  .* rbm. W       # Update Square [for EMF-TAP2]
    # end
    if approx == "tap3"
        rbm.W3 = rbm.W2 .* rbm.W        # Update Cube   [for EMF-TAP3]
    end
 end

function regularize_weight_gradient!(rbm::RBM,LearnRate::Float64;L2Penalty::Float64=NaN,L1Penalty::Float64=NaN,DropOutRate::Float64=NaN)
    ## Quadratic penalty on weights (Energy shrinkage)
    if !isnan(L2Penalty)
        axpy!(-LearnRate*L2Penalty,rbm.W,rbm.dW)
    end
    ## Linear penalty on weights (Sparsifying)
    if !isnan(L1Penalty)
        axpy!(-LearnRate*L1Penalty,sign(rbm.W),rbm.dW)
    end
    ## Dropout Regularization (restricted set of updates)
    if !isnan(DropOutRate)
        # Not yet implemented, so we do nothing.
        # TODO: Implement Drop-out, here.
    end
end

function get_negative_samples(rbm::RBM,vis_init::Mat{Float64},hid_init::Mat{Float64},approx::AbstractString, iterations::Int)
    if approx=="naive" || contains(approx,"tap")
        v_neg, h_neg = equilibrate(rbm,vis_init,hid_init; iterations=iterations, approx=approx)
    end

    if approx=="CD"
        # In the case of Gibbs/MCMC sampling, we will take the binary visible samples as the negative
        # visible samples, and the expectation (means) for the negative hidden samples.
        v_neg, _, _, h_neg = MCMC(rbm, hid_init; iterations=iterations, StartMode="hidden")
    end

    return v_neg, h_neg
end

function generate(rbm::RBM,vis_init::Mat{Float64},approx::AbstractString,SamplingIterations::Int)
    Nsamples = size(vis_init,2)
    Nhid     = size(rbm.hbias,1)
    h_init  = zeros(Nsamples,Nhid)

    if approx=="naive" || contains(approx,"tap")
        _, hid_mag = equilibrate(rbm,vis_init,hid_init; iterations=SamplingIterations, approx=approx)
    end

    if approx=="CD"
        _, hid_mag, _, _ = MCMC(rbm, vis_init; iterations=SamplingIterations, StartMode="visible")
    end

    samples,_ = sample_visibles(rbm,hid_mag)

    return reshape(samples,rbm.VisShape)
end

function fit_batch!(rbm::RBM, vis::Mat{Float64};
                    persistent=true, lr=0.1, NormalizationApproxIter=1,
                    weight_decay="none",decay_magnitude=0.01, approx="CD")

    # Determine how to acquire the positive samples based upon the persistence mode.
    v_pos = vis
    h_samples, h_pos = sample_hiddens(rbm,v_pos)
    # Set starting points in teh case of persistence
    if persistent
        if approx=="naive" || contains(approx,"tap")
            v_init = copy(rbm.persistent_chain_vis)
            h_init = copy(rbm.persistent_chain_hid)
        end
        if approx=="CD"
            v_init = vis               # A dummy setting
            h_init,_ = sample_hiddens(rbm,rbm.persistent_chain_vis)
        end
    else
        if approx=="naive" || contains(approx,"tap")
            v_init = vis
            h_init = h_pos
        end
        if approx=="CD"
            v_init = vis               # A dummy setting
            h_init = h_samples
        end
    end

    # Calculate the negative samples according to the desired approximation mode
    v_neg, h_neg = get_negative_samples(rbm,v_init,h_init,approx,NormalizationApproxIter)

    # If we are in persistent mode, update the chain accordingly
    if persistent
        copy!(rbm.persistent_chain_vis,v_neg)
        copy!(rbm.persistent_chain_hid,h_neg)
    end

    # Update on weights
    calculate_weight_gradient!(rbm,h_pos,v_pos,h_neg,v_neg,lr,approx=approx)
    if weight_decay == "l2"
        regularize_weight_gradient!(rbm,lr;L2Penalty=decay_magnitude)
    end
    if weight_decay == "l1"
        regularize_weight_gradient!(rbm,lr;L1Penalty=decay_magnitude)
    end
    update_weights!(rbm,approx)

    # Gradient update on biases
    rbm.hbias += vec(lr * (sum(h_pos, 2) - sum(h_neg, 2)))
    rbm.vbias += vec(lr * (sum(v_pos, 2) - sum(v_neg, 2)))

    return rbm
end

"""
    # Boltzmann.fit (training.jl)
    ## Function Call
        `fit(rbm::RBM, X::Mat{Float64}[, persistent, lr, batch_size, NormalizationApproxIter, weight_decay,
                                         decay_magnitude, validation,monitor_ever, monitor_vis,
                                         approx, persistent_start, save_params])`
    ## Description
    The core RBM training function. Learns the weights and biasings using
    either standard Contrastive Divergence (CD) or Persistent CD, depending on
    the user options.

    - *rbm:* RBM object, initialized by `RBM()`/`GRBM()`
    - *X:* Set of training vectors

    ### Optional Inputs
     - *persistent:* Whether or not to use persistent-CD [default=true]
     - *persistent_start:* At which epoch to start using the persistent chains. Only
                           applicable for the case that `persistent=true`.
                           [default=1]
     - *lr:* Learning rate [default=0.1]
     - *n_iter:* Number of training epochs [default=10]
     - *batch_size:* Minibatch size [default=100]
     - *NormalizationApproxIter:* Number of Gibbs sampling steps on the Markov Chain [default=1]
     - *weight_decay:* A string value representing the regularization to add to apply to the
                       weight magnitude during training {"none","l1","l2"}. [default="none"]
     - *decay_magnitude:* Relative importance assigned to the weight regularization. Smaller
                          values represent less regularization. Should be in range (0,1).
                          [default=0.01]
     - *validation:* An array of validation samples, e.g. a held out set of training data.
                     If passed, `fit` will also track generalization progress during training.
                     [default=empty-set]
     - *score_every:* Controls at which epoch the progress of the fit is monitored. Useful to
                      speed up the fit procedure if detailed progress monitoring is not required.
                      [default=5]
     - *save_progress:* Controls the saving of RBM parameters throughout the course of the training.
                     Should be passed as a tuple in the following manner:
                        (::AbstractString,::Int)
                     where the first field is the filename for the HDF5 used to save results and
                     the second field specifies how often to write the parameters to disk. All
                     results will be stored in the specified HDF5 file under the root headings
                        `Epochxxxx___weight`
                        `Epochxxxx___vbias`
                        `Epochxxxx___bias`
                     where `xxxx` specifies the epoch number in the `%04d` format.
                     [default=nothing]

    ## Returns
     - *::RBM* -- A trained RBM model.
     - *::Monitor* -- A Monitor structure containing information on the training progress over
                      epochs.
"""
function fit(rbm::RBM, X::Mat{Float64};
             persistent=true, lr=0.1, n_iter=10, batch_size=100, NormalizationApproxIter=1,
             weight_decay="none",decay_magnitude=0.01,validation=[],
             monitor_every=5,monitor_vis=false, approx="CD",
             persistent_start=1, save_progress=nothing)

    # TODO: This line needs to be changed to accomodate real-valued units
    @assert minimum(X) >= 0 && maximum(X) <= 1

    n_valid=0
    n_features = size(X, 1)
    n_samples = size(X, 2)
    n_hidden = size(rbm.W,1)
    n_batches = @compat Int(ceil(n_samples / batch_size))
    N = n_hidden+n_features

    # Check for the existence of a validation set
    flag_use_validation=false
    if length(validation)!=0
        flag_use_validation=true
        n_valid=size(validation,2)
    end

    # Create the historical monitor
    ProgressMonitor = Monitor(n_iter,monitor_every;monitor_vis=monitor_vis,
                                                   validation=flag_use_validation)

    # Print info to user
    m_ = rbm.momentum
    info("=====================================")
    info("RBM Training")
    info("=====================================")
    info("  + Training Samples:     $n_samples")
    info("  + Features:             $n_features")
    info("  + Hidden Units:         $n_hidden")
    info("  + Epochs to run:        $n_iter")
    info("  + Persistent ?:         $persistent")
    info("  + Training approx:      $approx")
    info("  + Momentum:             $m_")
    info("  + Learning rate:        $lr")
    info("  + Norm. Approx. Iters:  $NormalizationApproxIter")
    info("  + Weight Decay?:        $weight_decay")
    info("  + Weight Decay Mag.:    $decay_magnitude")
    info("  + Validation Set?:      $flag_use_validation")
    info("  + Validation Samples:   $n_valid")
    info("=====================================")

    # Scale the learning rate by the batch size
    lr=lr/batch_size

    # Random initialization of the persistent chains
    rbm.persistent_chain_vis,_ = random_columns(X,batch_size)
    rbm.persistent_chain_hid = ProbHidCondOnVis(rbm, rbm.persistent_chain_vis)

    use_persistent = false
    for itr=1:n_iter
        # Check to see if we can use persistence at this epoch
        use_persistent = itr>=persistent_start ? persistent : false

        tic()

        # Mini-batch fitting loop.
        @showprogress 1 "Fitting Batches..." for i=1:n_batches
            batch = X[:, ((i-1)*batch_size + 1):min(i*batch_size, end)]
            batch = full(batch)

            fit_batch!(rbm, batch; persistent=use_persistent,
                                   NormalizationApproxIter=NormalizationApproxIter,
                                   weight_decay=weight_decay,
                                   decay_magnitude=decay_magnitude,
                                   lr=lr, approx=approx)

        end

        # Get the average wall-time in µs
        walltime_µs=(toq()/n_batches/N)*1e6

        UpdateMonitor!(rbm,ProgressMonitor,X,itr;bt=walltime_µs,validation=validation,lr=lr,mo=rbm.momentum)
        ShowMonitor(rbm,ProgressMonitor,X,itr)

        # Attempt to save parameters if need be
        if save_progress != nothing
            if itr%save_progress[2]==0
                rootName = @sprintf("Epoch%04d",itr)
                if isfile(save_progress[1])
                    info("Appending Params...")
                    append_params(save_progress[1],rbm,rootName)
                else
                    info("Creating file and saving params...")
                    save_params(save_progress[1],rbm,rootName)
                end
            end
        end
    end

    if monitor_vis
        plt.close()
    end

    return rbm, ProgressMonitor
end

ext: jl


text:
using ToyAD
using Base.Test

# ======================
# scalar
# ======================
x1 = ToyAD.AD(2.0)
y1 = ToyAD.AD(1)
r1 = (x1+y1)*(y1 + ToyAD.AD(1))
bp1 = ToyAD.backprop(r1, true)
@test r1.value == ones(1,1)*6
@test x1.grad == ones(1,1)*2
@test y1.grad == ones(1,1)*5

# ======================
# matrix
# ======================
m1 = ToyAD.AD(zeros(3,2))
m1.value[1,1] = 1.; m1.value[1,2] = 2.
m1.value[2,1] = 3.; m1.value[2,2] = 4.
m1.value[3,1] = 5.; m1.value[3,2] = 6.

m2 = ToyAD.AD(zeros(2,3))
m2.value[1,1] = 2.; m2.value[1,2] = 3.; m2.value[1,3] = 4.
m2.value[2,1] = 5.; m2.value[2,2] = 6.; m2.value[2,3] = 7.

# addition
m3 = m1 + m1
m3.grad[1,1] = .1; m3.grad[1,2] = .2
m3.grad[2,1] = .3; m3.grad[2,2] = .4
m3.grad[3,1] = .5; m3.grad[3,2] = .6

rmAdd = ToyAD.backprop(m3, false)
@test m1.grad[1,1] == 0.2
@test m1.grad[1,2] == 0.4
@test m1.grad[2,1] == 0.6
@test m1.grad[2,2] == 0.8
@test m1.grad[3,1] == 1.0
@test m1.grad[3,2] == 1.2

# multiplication
# mul test
m1.grad[:] = 0. # reset gradient matrices
m2.grad[:] = 0. # reset  gradient matrices
m3 = m1* m2
@test m3.value[1,1] == 12.
@test m3.value[1,2] == 15.
@test m3.value[1,3] == 18.
@test m3.value[2,1] == 26.
@test m3.value[2,2] == 33.
@test m3.value[2,3] == 40.
@test m3.value[3,1] == 40.
@test m3.value[3,2] == 51.
@test m3.value[3,3] == 62.

m3.grad[1,1] = .1; m3.grad[1,2] = .2; m3.grad[1,3] = .3
m3.grad[2,1] = .4; m3.grad[2,2] = .5; m3.grad[2,3] = .6
m3.grad[3,1] = .7; m3.grad[3,2] = .8; m3.grad[3,3] = .9

rmMul = ToyAD.backprop(m3, false)
# m1 gradient tests
@test m1.grad[1,1] == 2.
@test m1.grad[1,2] == 3.8000000000000003
@test m1.grad[2,1] == 4.699999999999999
@test m1.grad[2,2] == 9.2
@test m1.grad[3,1] == 7.4
@test m1.grad[3,2] == 14.600000000000001

# m2 gradient tests
@test m2.grad[1,1] == 4.800000000000001
@test m2.grad[1,2] == 5.7
@test m2.grad[1,3] == 6.6
@test m2.grad[2,1] == 5.999999999999999
@test m2.grad[2,2] == 7.200000000000001
@test m2.grad[2,3] == 8.4

# reul() tests
m4 = ToyAD.AD(zeros(3,2))
m4.value[1,1] = 1.; m4.value[1,2] =-2.
m4.value[2,1] =-3.; m4.value[2,2] = 4.
m4.value[3,1] = 5.; m4.value[3,2] =-6.
m5 = ToyAD.relu(m4)
@test m5.value[1,1] == 1.
@test m5.value[1,2] == 0.
@test m5.value[2,1] == 0.
@test m5.value[2,2] == 4.
@test m5.value[3,1] == 5.
@test m5.value[3,2] == 0.

m5.grad[1,1] =-.1; m5.grad[1,2] = .2
m5.grad[2,1] = .3; m5.grad[2,2] = .4
m5.grad[3,1] = .5; m5.grad[3,2] = .6

rmRelu = ToyAD.backprop(m5, false)

@test m4.grad[1,1] == -0.1
@test m4.grad[1,2] == 0.
@test m4.grad[2,1] == 0.
@test m4.grad[2,2] == 0.4
@test m4.grad[3,1] == 0.5
@test m4.grad[3,2] == 0.


text:
include("fprop.jl")
include("load_data.jl")
include("helper_functions.jl")
include("cost.jl")
using Devectorize

function train(epochs)
#=
Inputs:
  epochs: Number of epochs to run
Output:
  model: A structure containing the learned weights and biases and vocabulary
=#

# SET HYPERPARAMETERS
batchsize = 100
learning_rate = 0.1
momentum = 0.9
numhid1 = 50
numhid2 = 200
init_wt = 0.01  # This is the standard deviation of the normal distributios

# LOAD DATA
(train_input, train_target, valid_input, valid_target, test_input, test_target,
vocab) = load_data(batchsize)
numwords, batchsize, numbatches = size(train_input)
vocab_size = size(vocab, 2)

# Get initial weights and deltas
(word_embedding_weights,
 embed_to_hid_weights,
 hid_to_output_weights,
 hid_bias,
 output_bias,
 word_embedding_weights_delta,
 word_embedding_weights_gradient,
 embed_to_hid_weights_delta,
 hid_to_output_weights_delta,
 hid_bias_delta,
 output_bias_delta) = network_initializer(init_wt, numhid1, numhid2, numwords, vocab_size)

# Create expansion matrix for output layer error calculation
expansion_matrix = eye(vocab_size)

# TRAIN
for epoch = 1:epochs
  println("Epoch: ",epoch)
  for m = 1:numbatches
    input_batch = train_input[:, :, m]
    target_batch = train_target[:, :, m]

    # FORWARD PROPAGATE
    (embedding_layer_state, hidden_layer_state, output_layer_state) =
      fprop(input_batch,
            word_embedding_weights, embed_to_hid_weights,
            hid_to_output_weights, hid_bias, output_bias)
    # Compute the error derivative of the output layer
    expanded_target_batch = output_mapper(expansion_matrix, target_batch)
    @devec error_deriv = output_layer_state - expanded_target_batch

    # Measure the loss function
    CE = cost(expanded_target_batch, output_layer_state, batchsize)

    # BACK PROPAGATION

    # OUTPUT LAYER
    hid_to_output_weights_gradient = hidden_layer_state * error_deriv'
    output_bias_gradient = sum_rows(error_deriv)
    back_propogated_deriv_1 = hid_to_output_weights * error_deriv
    @devec back_propogated_deriv_1 = back_propogated_deriv_1 .* hidden_layer_state .* (1 - hidden_layer_state)

    # HIDDEN LAYER
    embed_to_hid_weights_gradient = embedding_layer_state * back_propogated_deriv_1'
    hid_bias_gradient = sum_rows(back_propogated_deriv_1)
    back_propogated_deriv_2 = embed_to_hid_weights * back_propogated_deriv_1

    word_embedding_weights_gradient[:, :] = 0
    # EMBEDDING LAYER
    for w = 1:numwords
      word_embedding_weights_gradient +=
        expansion_matrix[:, vec(input_batch[w, :])] * back_propogated_deriv_2[1 + (w-1) * numhid1: w * numhid1,:]'
    end

    # UPDATE NETWORK WEIGHTS AND BIASES
    @devec word_embedding_weights_delta = momentum .* word_embedding_weights_delta + word_embedding_weights_gradient ./ batchsize
    @devec word_embedding_weights = word_embedding_weights - learning_rate .* word_embedding_weights_delta

    @devec embed_to_hid_weights_delta = momentum .* embed_to_hid_weights_delta + embed_to_hid_weights_gradient ./ batchsize
    @devec embed_to_hid_weights = embed_to_hid_weights- learning_rate .* embed_to_hid_weights_delta

    @devec hid_to_output_weights_delta = momentum .* hid_to_output_weights_delta + hid_to_output_weights_gradient ./ batchsize
    @devec hid_to_output_weights = hid_to_output_weights - learning_rate .* hid_to_output_weights_delta;

    @devec hid_bias_delta = momentum .* hid_bias_delta + hid_bias_gradient ./ batchsize
    @devec hid_bias = hid_bias - learning_rate .* hid_bias_delta;

    @devec output_bias_delta = momentum .* output_bias_delta + output_bias_gradient ./ batchsize
    @devec output_bias = output_bias - learning_rate .* output_bias_delta
  end
    (embedding_layer_state, hidden_layer_state, output_layer_state) =
    fprop(valid_input,
          word_embedding_weights, embed_to_hid_weights,
          hid_to_output_weights, hid_bias, output_bias)
    datasetsize = size(valid_input, 2)
    expanded_valid_target = output_mapper(expansion_matrix, valid_target)
    CE = cost(expanded_valid_target, output_layer_state, datasetsize)
    println("Validation error after epoch: ", epoch, "is ", CE)
end

end
train(5)

max_stars_repo_name: kirnap/Machine-Learning-in-Julia

max_issues_repo_name: kirnap/Machine-Learning-in-Julia

max_forks_repo_name: kirnap/Machine-Learning-in-Julia


text:
###
# Test by simulating OU with known transition density
###
using Makie, ZigZagBoomerang, SparseArrays, LinearAlgebra
# using CairoMakie
include("faberschauder.jl")

# Drift
const β = 1.0
b(x) = -β*x
# First derivative
b′(x) = -β
# Second derivative
b″(x) = zero(x)

L = 7
n = (2 << L) + 1
T = 2.0 # length diffusion bridge
ξ0 = 0randn(n)
u, v = 2.5, -0.5  # initial and fianl point
ξ0[1] = u/sqrt(T)
ξ0[end] = v/sqrt(T)
c = 2 .^ (lvl.(1:n, L))
#c[end] = 1.0
c[1] = 0.0
θ0 = rand((-1.0, 1.0), n)
#θ0[end] =
θ0[1] =  0.0 # fix final point
T′ = 4000.0 # final clock of the pdmp

Γ = sparse(1.0I, n, n)
trace, (t, ξ, θ), (acc, num), c = @time spdmp(∇ϕmoving, 0.0, ξ0,
    θ0, T′, c, ZigZag(Γ, ξ0*0), SelfMoving(), L, T, adapt=false);
ts, ξs = splitpairs(discretize(trace, T′/n))
S = T*(0:n)/(n+1)

using Makie
p1 = lines(S, [dotψ(ξ, s, L, T) for s in S], linewidth=0.3)
for ξ in ξs[1:5:end]
    lines!(p1, S, [dotψ(ξ, s, L, T) for s in S], linewidth=0.3)
end
display(p1)

p3 = hbox([lines(ts, getindex.(ξs, i)) for i in [2,3,5,9,17,(n+1)÷2 + 1, n]]...)

vbox(p1, p3)

using Statistics
p4 = lines(S, [mean(dotψ(ξ, s, L, T) for ξ in ξs) for s in S])
p4 = lines!(p4, S, @.(u*exp(-β*S)))
p0 = vbox(hbox(p1, p4), p3)
save("figures/diffbridges2.png", p0)


text:
using FranklinUtils

@env function box(md; title="", color="purple")
    return html("""

$title

        """) * md * html("""

        """)
end

@lx function avatar(; link="", img="", name="", affil="")
    return html("""

            $name

              $affil

        """)
end

# NAV_YEAR tools

curyear() = parse(Int, match(r"^[\/|\\]?(20\d{2})[\/||\\]", locvar(:fd_rpath)::String).captures[1])

function hfun_curyear()
    return "$(curyear())"
end

function hfun_previous_editions()
    year = curyear()
    pre_years = year-1:-1:2014
    io = IOBuffer()
    for (i, y) in enumerate(pre_years)
        write(io, """$y""")
        i == length(pre_years) || write(io, "/")
    end
    prevs = String(take!(io))
    return """

          Previously: $prevs

        """
end

@env function centered(md; title="Gold")
    return html("""

$title

        """) * md * html("""

        """)
end

@lx function sponsor(; link="", img="", name="", level=1)
    return html("""

            """)
end


text:
using Pkg
Pkg.activate(".")

using Genie, Genie.Router, Genie.Renderer, Genie.Renderer.Html

form = """
<form action="/" method="POST" enctype="multipart/form-data">
  <input type="text" name="greeting" value="hello genie" />

  <input type="file" name="fileupload" />

  <input type="submit" value="Submit" />
</form>
"""

route("/") do
  html(form)
end

route("/", method = POST) do
  for (name,file) in @params(:FILES)
    write(file.name, IOBuffer(file.data))
  end

  write("-" * @params(:FILES)["fileupload"].name, IOBuffer(@params(:FILES)["fileupload"].data))

  @show @params(:greeting)

  @params(:greeting)
end

Genie.AppServer.startup(; open_browser = false, async = false)


text:
@testset "Edit distance" begin

    gtri = random_regular_graph(3,2)
    gquad = random_regular_graph(4,2)
    gpent = random_regular_graph(5,2)

    for triangle in testgraphs(gtri), quadrangle in testgraphs(gquad), pentagon in testgraphs(gpent)
        d, λ = @inferred(edit_distance(triangle, quadrangle, subst_cost=MinkowskiCost(1:3,1:4)))
        @test d == 1.0
        @test λ == Tuple[(1,1),(2,2),(3,3),(0,4)]

        d, λ = @inferred(edit_distance(quadrangle, triangle, subst_cost=MinkowskiCost(1:4,1:3)))
        @test d == 1.0
        @test λ == Tuple[(1,1),(2,2),(3,3),(4,0)]

        d, λ = @inferred(edit_distance(triangle, pentagon, subst_cost=MinkowskiCost(1:3,1:5)))
        @test d == 2.0
        @test λ == Tuple[(1,1),(2,2),(3,3),(0,4),(0,5)]

        d, λ = @inferred(edit_distance(pentagon, triangle, subst_cost=MinkowskiCost(1:5,1:3)))
        @test d == 2.0
        @test λ == Tuple[(1,1),(2,2),(3,3),(4,0),(5,0)]
    end
    cost = @inferred(MinkowskiCost(1:3,1:3))
    bcost = @inferred(BoundedMinkowskiCost(1:3,1:3))
    for i=1:3
      @test cost(i,i) == 0.
      @test bcost(i,i) == 2/3
    end

    g1c = CompleteGraph(4)
    g2c = CompleteGraph(4)
    rem_edge!(g2c, 1, 2)
    for g1 in testgraphs(g1c), g2 in testgraphs(g2c)
        d, λ = @inferred(edit_distance(g1, g2))
        @test d == 2.0
        @test λ == Tuple[(1,1),(2,2),(3,3),(4,4)]
    end
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
# Julia wrapper for header: /Users/daly/.julia/v0.6/Homebrew/deps/usr/Cellar/sdl2_mixer/2.0.2_3/include/SDL2/SDL_mixer.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

include("SDL_mixer_h.jl")

function Mix_Linked_Version()
    ccall((:Mix_Linked_Version, libsdl2_mixer), Ptr{SDL_version}, ())
end

function Mix_Init(flags::Cint)
    ccall((:Mix_Init, libsdl2_mixer), Cint, (Cint,), flags)
end

function Mix_Quit()
    ccall((:Mix_Quit, libsdl2_mixer), Cvoid, ())
end

function Mix_OpenAudio(frequency::Cint, format::Uint16, channels::Cint, chunksize::Cint)
    ccall((:Mix_OpenAudio, libsdl2_mixer), Cint, (Cint, Uint16, Cint, Cint), frequency, format, channels, chunksize)
end

function Mix_OpenAudioDevice(frequency::Cint, format::Uint16, channels::Cint, chunksize::Cint, device, allowed_changes::Cint)
    ccall((:Mix_OpenAudioDevice, libsdl2_mixer), Cint, (Cint, Uint16, Cint, Cint, Cstring, Cint), frequency, format, channels, chunksize, device, allowed_changes)
end

function Mix_AllocateChannels(numchans::Cint)
    ccall((:Mix_AllocateChannels, libsdl2_mixer), Cint, (Cint,), numchans)
end

function Mix_QuerySpec(frequency, format, channels)
    ccall((:Mix_QuerySpec, libsdl2_mixer), Cint, (Ptr{Cint}, Ptr{Uint16}, Ptr{Cint}), frequency, format, channels)
end

function Mix_LoadWAV_RW(src, freesrc::Cint)
    ccall((:Mix_LoadWAV_RW, libsdl2_mixer), Ptr{Mix_Chunk}, (Ptr{RWops}, Cint), src, freesrc)
end

function Mix_LoadMUS(file)
    ccall((:Mix_LoadMUS, libsdl2_mixer), Ptr{Mix_Music}, (Cstring,), file)
end

function Mix_LoadMUS_RW(src, freesrc::Cint)
    ccall((:Mix_LoadMUS_RW, libsdl2_mixer), Ptr{Mix_Music}, (Ptr{RWops}, Cint), src, freesrc)
end

function Mix_LoadMUSType_RW(src, _type::Mix_MusicType, freesrc::Cint)
    ccall((:Mix_LoadMUSType_RW, libsdl2_mixer), Ptr{Mix_Music}, (Ptr{RWops}, Mix_MusicType, Cint), src, _type, freesrc)
end

function Mix_QuickLoad_WAV(mem)
    ccall((:Mix_QuickLoad_WAV, libsdl2_mixer), Ptr{Mix_Chunk}, (Ptr{Uint8},), mem)
end

function Mix_QuickLoad_RAW(mem, len::Uint32)
    ccall((:Mix_QuickLoad_RAW, libsdl2_mixer), Ptr{Mix_Chunk}, (Ptr{Uint8}, Uint32), mem, len)
end

function Mix_FreeChunk(chunk)
    ccall((:Mix_FreeChunk, libsdl2_mixer), Cvoid, (Ptr{Mix_Chunk},), chunk)
end

function Mix_FreeMusic(music)
    ccall((:Mix_FreeMusic, libsdl2_mixer), Cvoid, (Ptr{Mix_Music},), music)
end

function Mix_GetNumChunkDecoders()
    ccall((:Mix_GetNumChunkDecoders, libsdl2_mixer), Cint, ())
end

function Mix_GetChunkDecoder(index::Cint)
    ccall((:Mix_GetChunkDecoder, libsdl2_mixer), Cstring, (Cint,), index)
end

function Mix_HasChunkDecoder(name)
    ccall((:Mix_HasChunkDecoder, libsdl2_mixer), bool, (Cstring,), name)
end

function Mix_GetNumMusicDecoders()
    ccall((:Mix_GetNumMusicDecoders, libsdl2_mixer), Cint, ())
end

function Mix_GetMusicDecoder(index::Cint)
    ccall((:Mix_GetMusicDecoder, libsdl2_mixer), Cstring, (Cint,), index)
end

function Mix_HasMusicDecoder(name)
    ccall((:Mix_HasMusicDecoder, libsdl2_mixer), bool, (Cstring,), name)
end

function Mix_GetMusicType(music)
    ccall((:Mix_GetMusicType, libsdl2_mixer), Mix_MusicType, (Ptr{Mix_Music},), music)
end

function Mix_SetPostMix(mix_func, arg)
    ccall((:Mix_SetPostMix, libsdl2_mixer), Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}), mix_func, arg)
end

function Mix_HookMusic(mix_func, arg)
    ccall((:Mix_HookMusic, libsdl2_mixer), Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}), mix_func, arg)
end

function Mix_HookMusicFinished(music_finished)
    ccall((:Mix_HookMusicFinished, libsdl2_mixer), Cvoid, (Ptr{Cvoid},), music_finished)
end

function Mix_GetMusicHookData()
    ccall((:Mix_GetMusicHookData, libsdl2_mixer), Ptr{Cvoid}, ())
end

function Mix_ChannelFinished(channel_finished)
    ccall((:Mix_ChannelFinished, libsdl2_mixer), Cvoid, (Ptr{Cvoid},), channel_finished)
end

function Mix_RegisterEffect(chan::Cint, f::Mix_EffectFunc_t, d::Mix_EffectDone_t, arg)
    ccall((:Mix_RegisterEffect, libsdl2_mixer), Cint, (Cint, Mix_EffectFunc_t, Mix_EffectDone_t, Ptr{Cvoid}), chan, f, d, arg)
end

function Mix_UnregisterEffect(channel::Cint, f::Mix_EffectFunc_t)
    ccall((:Mix_UnregisterEffect, libsdl2_mixer), Cint, (Cint, Mix_EffectFunc_t), channel, f)
end

function Mix_UnregisterAllEffects(channel::Cint)
    ccall((:Mix_UnregisterAllEffects, libsdl2_mixer), Cint, (Cint,), channel)
end

function Mix_SetPanning(channel::Cint, left::Uint8, right::Uint8)
    ccall((:Mix_SetPanning, libsdl2_mixer), Cint, (Cint, Uint8, Uint8), channel, left, right)
end

function Mix_SetPosition(channel::Cint, angle::Sint16, distance::Uint8)
    ccall((:Mix_SetPosition, libsdl2_mixer), Cint, (Cint, Sint16, Uint8), channel, angle, distance)
end

function Mix_SetDistance(channel::Cint, distance::Uint8)
    ccall((:Mix_SetDistance, libsdl2_mixer), Cint, (Cint, Uint8), channel, distance)
end

function Mix_SetReverseStereo(channel::Cint, flip::Cint)
    ccall((:Mix_SetReverseStereo, libsdl2_mixer), Cint, (Cint, Cint), channel, flip)
end

function Mix_ReserveChannels(num::Cint)
    ccall((:Mix_ReserveChannels, libsdl2_mixer), Cint, (Cint,), num)
end

function Mix_GroupChannel(which::Cint, tag::Cint)
    ccall((:Mix_GroupChannel, libsdl2_mixer), Cint, (Cint, Cint), which, tag)
end

function Mix_GroupChannels(from::Cint, to::Cint, tag::Cint)
    ccall((:Mix_GroupChannels, libsdl2_mixer), Cint, (Cint, Cint, Cint), from, to, tag)
end

function Mix_GroupAvailable(tag::Cint)
    ccall((:Mix_GroupAvailable, libsdl2_mixer), Cint, (Cint,), tag)
end

function Mix_GroupCount(tag::Cint)
    ccall((:Mix_GroupCount, libsdl2_mixer), Cint, (Cint,), tag)
end

function Mix_GroupOldest(tag::Cint)
    ccall((:Mix_GroupOldest, libsdl2_mixer), Cint, (Cint,), tag)
end

function Mix_GroupNewer(tag::Cint)
    ccall((:Mix_GroupNewer, libsdl2_mixer), Cint, (Cint,), tag)
end

function Mix_PlayChannelTimed(channel::Cint, chunk, loops::Cint, ticks::Cint)
    ccall((:Mix_PlayChannelTimed, libsdl2_mixer), Cint, (Cint, Ptr{Mix_Chunk}, Cint, Cint), channel, chunk, loops, ticks)
end

function Mix_PlayMusic(music, loops::Cint)
    ccall((:Mix_PlayMusic, libsdl2_mixer), Cint, (Ptr{Mix_Music}, Cint), music, loops)
end

function Mix_FadeInMusic(music, loops::Cint, ms::Cint)
    ccall((:Mix_FadeInMusic, libsdl2_mixer), Cint, (Ptr{Mix_Music}, Cint, Cint), music, loops, ms)
end

function Mix_FadeInMusicPos(music, loops::Cint, ms::Cint, position::Cdouble)
    ccall((:Mix_FadeInMusicPos, libsdl2_mixer), Cint, (Ptr{Mix_Music}, Cint, Cint, Cdouble), music, loops, ms, position)
end

function Mix_FadeInChannelTimed(channel::Cint, chunk, loops::Cint, ms::Cint, ticks::Cint)
    ccall((:Mix_FadeInChannelTimed, libsdl2_mixer), Cint, (Cint, Ptr{Mix_Chunk}, Cint, Cint, Cint), channel, chunk, loops, ms, ticks)
end

function Mix_Volume(channel::Cint, volume::Cint)
    ccall((:Mix_Volume, libsdl2_mixer), Cint, (Cint, Cint), channel, volume)
end

function Mix_VolumeChunk(chunk, volume::Cint)
    ccall((:Mix_VolumeChunk, libsdl2_mixer), Cint, (Ptr{Mix_Chunk}, Cint), chunk, volume)
end

function Mix_VolumeMusic(volume::Cint)
    ccall((:Mix_VolumeMusic, libsdl2_mixer), Cint, (Cint,), volume)
end

function Mix_HaltChannel(channel::Cint)
    ccall((:Mix_HaltChannel, libsdl2_mixer), Cint, (Cint,), channel)
end

function Mix_HaltGroup(tag::Cint)
    ccall((:Mix_HaltGroup, libsdl2_mixer), Cint, (Cint,), tag)
end

function Mix_HaltMusic()
    ccall((:Mix_HaltMusic, libsdl2_mixer), Cint, ())
end

function Mix_ExpireChannel(channel::Cint, ticks::Cint)
    ccall((:Mix_ExpireChannel, libsdl2_mixer), Cint, (Cint, Cint), channel, ticks)
end

function Mix_FadeOutChannel(which::Cint, ms::Cint)
    ccall((:Mix_FadeOutChannel, libsdl2_mixer), Cint, (Cint, Cint), which, ms)
end

function Mix_FadeOutGroup(tag::Cint, ms::Cint)
    ccall((:Mix_FadeOutGroup, libsdl2_mixer), Cint, (Cint, Cint), tag, ms)
end

function Mix_FadeOutMusic(ms::Cint)
    ccall((:Mix_FadeOutMusic, libsdl2_mixer), Cint, (Cint,), ms)
end

function Mix_FadingMusic()
    ccall((:Mix_FadingMusic, libsdl2_mixer), Mix_Fading, ())
end

function Mix_FadingChannel(which::Cint)
    ccall((:Mix_FadingChannel, libsdl2_mixer), Mix_Fading, (Cint,), which)
end

function Mix_Pause(channel::Cint)
    ccall((:Mix_Pause, libsdl2_mixer), Cvoid, (Cint,), channel)
end

function Mix_Resume(channel::Cint)
    ccall((:Mix_Resume, libsdl2_mixer), Cvoid, (Cint,), channel)
end

function Mix_Paused(channel::Cint)
    ccall((:Mix_Paused, libsdl2_mixer), Cint, (Cint,), channel)
end

function Mix_PauseMusic()
    ccall((:Mix_PauseMusic, libsdl2_mixer), Cvoid, ())
end

function Mix_ResumeMusic()
    ccall((:Mix_ResumeMusic, libsdl2_mixer), Cvoid, ())
end

function Mix_RewindMusic()
    ccall((:Mix_RewindMusic, libsdl2_mixer), Cvoid, ())
end

function Mix_PausedMusic()
    ccall((:Mix_PausedMusic, libsdl2_mixer), Cint, ())
end

function Mix_SetMusicPosition(position::Cdouble)
    ccall((:Mix_SetMusicPosition, libsdl2_mixer), Cint, (Cdouble,), position)
end

function Mix_Playing(channel::Cint)
    ccall((:Mix_Playing, libsdl2_mixer), Cint, (Cint,), channel)
end

function Mix_PlayingMusic()
    ccall((:Mix_PlayingMusic, libsdl2_mixer), Cint, ())
end

function Mix_SetMusicCMD(command)
    ccall((:Mix_SetMusicCMD, libsdl2_mixer), Cint, (Cstring,), command)
end

function Mix_SetSynchroValue(value::Cint)
    ccall((:Mix_SetSynchroValue, libsdl2_mixer), Cint, (Cint,), value)
end

function Mix_GetSynchroValue()
    ccall((:Mix_GetSynchroValue, libsdl2_mixer), Cint, ())
end

function Mix_SetSoundFonts(paths)
    ccall((:Mix_SetSoundFonts, libsdl2_mixer), Cint, (Cstring,), paths)
end

function Mix_GetSoundFonts()
    ccall((:Mix_GetSoundFonts, libsdl2_mixer), Cstring, ())
end

function Mix_EachSoundFont(_function, data)
    ccall((:Mix_EachSoundFont, libsdl2_mixer), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), _function, data)
end

function Mix_GetChunk(channel::Cint)
    ccall((:Mix_GetChunk, libsdl2_mixer), Ptr{Mix_Chunk}, (Cint,), channel)
end

function Mix_CloseAudio()
    ccall((:Mix_CloseAudio, libsdl2_mixer), Cvoid, ())
end

max_stars_repo_name: Kyjor/julia-game

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: Kyjor/julia-game

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: Kyjor/julia-game

max_forks_repo_licenses: MIT Unlicense


text:
# using Base.Test
# using MIToS.MSA

print("""

Tests for Annotations
=====================
""")

print("""

Annotations() for empty annotations
""")
let annot = Annotations()
  @test length(annot.file) == 0
  @test length(annot.sequences) == 0
  @test length(annot.columns) == 0
  @test length(annot.residues) == 0

  print("""
ncolums
""")
  @test ncolumns(annot) == -1

  print("""

Getters & Setters
-----------------
""")
  setannotresidue!(annot, "O31698/18-71", "SS", "CCCHHHHHHHHHHHHHHHEEEEEEEEEEEEEEEEHHH")

  @test_throws ErrorException setannotresidue!(annot, "O31698/18-71", randstring(51),    "CCCHHHHHHHHHHHHHHHEEEEEEEEEEEEEEEEHHH")
  @test_throws ErrorException setannotresidue!(annot, "O31698/18-71", "My Feature Name", "CCCHHHHHHHHHHHHHHHEEEEEEEEEEEEEEEEHHH")

  @test ncolumns(annot) == 37

  setannotfile!(annot, "AC", "PF00571")
  setannotcolumn!(annot, "SS_cons", "CCCCCHHHHHHHHHHHHHEEEEEEEEEEEEEEEEEEH")
  setannotsequence!(annot, "O31698/88-139", "OS", "Bacillus subtilis")

  @test getannotfile(annot, "AC") == "PF00571"
  @test getannotcolumn(annot, "SS_cons") == "CCCCCHHHHHHHHHHHHHEEEEEEEEEEEEEEEEEEH"
  @test getannotsequence(annot, "O31698/88-139", "OS") == "Bacillus subtilis"
  @test getannotresidue(annot, "O31698/18-71", "SS") == "CCCHHHHHHHHHHHHHHHEEEEEEEEEEEEEEEEHHH"

  @test getannotfile(annot, "An", "Default") == "Default"
  @test getannotcolumn(annot, "Other", "Default") == "Default"
  @test getannotsequence(annot, "O31698/1-88", "OS", "Default") == "Default"
  @test getannotresidue(annot, "O31698/1-88", "SS", "Default") == "Default"

  @test ncolumns(annot) == 37
  @test_throws DimensionMismatch setannotresidue!(annot, "O31698/18-71", "AS", "__*__")
  @test_throws DimensionMismatch setannotcolumn!(annot, "SS_cons", "---CCCCCHHHHHHHHHHHHHEEEEEEEEEEEEEEEEEEH---")

  print("""

copy, empty, isempty & empty!
-----------------------------
""")
  copy_annot = copy(annot)
  empty!(copy_annot)
  @test ncolumns(annot) == 37
  @test ncolumns(copy_annot) == -1

  @test !isempty(annot)
  @test isempty(copy_annot)

  no_annotations = empty(Annotations)
  @test ncolumns(no_annotations) == -1
  @test isempty(no_annotations)

  print("""

Filters
-------
""")
  filtersequences!(annot, IndexedArray(["O31698/88-139", "O31698/18-71"]), [false, true])
  @test length( getannotsequence(annot) ) == 0
  filtersequences!(annot, IndexedArray(["O31698/88-139", "O31698/18-71"]), [true, false])
  @test length( getannotresidue(annot) ) == 0

  mask = collect("CCCCCHHHHHHHHHHHHHEEEEEEEEEEEEEEEEEEH") .!= 'E'
  filtercolumns!(annot, mask)
  @test ncolumns(annot) == 19
  @test getannotcolumn(annot, "SS_cons") == "CCCCCHHHHHHHHHHHHHH"

  filtercolumns!(annot, [1,2,19])
  @test ncolumns(annot) == 3
  @test getannotcolumn(annot, "SS_cons") == "CCH"

end


text:
# Contains all the data relative to quadrature in the sphere
struct SphereQuadraturervec
    r1Vec::Array{Float64}
    r2Vec::Array{Float64}
    r3Vec::Array{Float64}
    wVec::Array{Float64}
end


text:
"""
    Image

An image, taken though a particular filter band.
"""
mutable struct Image{B<:AbstractMatrix, P<:AbstractPSFMap}
    H::Int  # image height in pixels
    W::Int  # image width in pixels

    # An HxW matrix of pixel intensities, in raw electron counts (nelec).
    pixels::Matrix{Float32}

    b::Int  # band id (values from 1 to 5).
    wcs::WCSTransform  # Transform between pixels and world coordinates
    psf::Vector{PsfComponent}  # The PSF at the center of the image

    # The background intensity in nanomaggies. (varies by position)
    sky::B

    # The expected number of photons contributed to this image
    # by a source 1 nanomaggie in brightness. (varies by row)
    nelec_per_nmgy::Array{Float32, 1}

    # The image PSF map: a callable f(x, y) -> Matrix that takes a pixel
    # coordinate and returns the rasterized PSF image at that coordinate,
    # with the PSF centered in the image.
    psfmap::P

    Image{B,P}(pixels::Matrix{Float32},
               b::Int,
               wcs::WCSTransform,
               psf::Vector{PsfComponent},
               sky::B,
               nelec_per_nmgy::Array{Float32, 1},
               psfmap::P) where {B<:AbstractMatrix, P<:AbstractPSFMap} =
        new(size(pixels, 1), size(pixels, 2), pixels, b, wcs, psf, sky,
            nelec_per_nmgy, psfmap)
end

Image(pixels::Matrix{Float32},
      b::Int,
      wcs::WCSTransform,
      psf::Vector{PsfComponent},
      sky::B,
      nelec_per_nmgy::Array{Float32, 1},
      psfmap::P) where {B<:AbstractMatrix, P<:AbstractPSFMap} =
    Image{B,P}(pixels, b, wcs, psf, sky, nelec_per_nmgy, psfmap)

# TODO: better name for this.
"""
    calibrated_pixels(im::Image)

Calibrated, sky-subtracted pixel values in nmgy.
"""
calibrated_pixels(im::Image) = im.pixels ./ im.nelec_per_nmgy .- im.sky

# The code below lets us JLD serialize images instances.
# Without this code we get an error for trying to serialize C pointers from WCS
# and some problems for StaticArrays too.
struct SerializedWCS
    header::String
end

JLD.writeas(wcs::WCSTransform) = SerializedWCS(WCS.to_header(wcs))
JLD.readas(s_wcs::SerializedWCS) = WCS.from_header(s_wcs.header)[1]

struct PsfComponentSerial
    alphaBar::Float64
    xiBar::Vector{Float64}
    tauBar::Matrix{Float64}
end

JLD.writeas(pcs::Vector{PsfComponent}) = begin
    [PsfComponentSerial(pc.alphaBar, pc.xiBar, pc.tauBar) for pc in pcs]
end
JLD.readas(pcs::Vector{PsfComponentSerial}) = begin
    [PsfComponent(pc.alphaBar,
                  convert(SVector{2,Float64}, pc.xiBar),
                  convert(SMatrix{2,2,Float64,4}, pc.tauBar)) for pc in pcs]
end


text:
abstract type Strength end
struct Classical{T} <: Strength
    θ::T
end
Classical(;θ = 0.25) = Classical(θ)

function (c::Classical)(At::SparseMatrixCSC{Tv,Ti}) where {Ti,Tv}

    θ = c.θ

    m, n = size(At)
    T = deepcopy(At)

    for i = 1:n
        _m = find_max_off_diag(T, i)
        threshold = θ * _m
        for j in nzrange(T, i)
            row = T.rowval[j]
            val = T.nzval[j]

            if row != i
                if abs(val) >= threshold
                    T.nzval[j] = abs(val)
                else
                    T.nzval[j] = 0
                end
            end

        end
    end

    dropzeros!(T)

    scale_cols_by_largest_entry!(T)

    adjoint(T), T
end

function find_max_off_diag(A, i)
    m = zero(eltype(A))
    for j in nzrange(A, i)
        row = A.rowval[j]
        val = A.nzval[j]
        if row != i
            m = max(m, abs(val))
        end
    end
    m
end

function find_max(A, i)
    m = zero(eltype(A))
    for j in nzrange(A, i)
        row = A.rowval[j]
        val = A.nzval[j]
        m = max(m, val)
    end
    m
end

function scale_cols_by_largest_entry!(A::SparseMatrixCSC)
    n = size(A, 1)
    for i = 1:n
        _m = find_max(A, i)
        for j in nzrange(A, i)
            A.nzval[j] /= _m
        end
    end
    A
end

struct SymmetricStrength{T} <: Strength
    θ::T
end
SymmetricStrength() = SymmetricStrength(0.)

function (s::SymmetricStrength{T})(A, bsr_flag = false) where {T}

    θ = s.θ

    if bsr_flag && θ == 0
        S = SparseMatrixCSC(size(A)...,
                    A.colptr, A.rowval, ones(eltype(A), size(A.rowval)))
        return S, S
    else
        S = deepcopy(A)
    end
    n = size(A, 1)
    diags = Vector{eltype(A)}(undef, n)

    for i = 1:n
        diag = zero(eltype(A))
        for j in nzrange(A, i)
            row = A.rowval[j]
            val = A.nzval[j]
            if row == i
                diag += val
            end
        end
        diags[i] = norm(diag)
    end

    for i = 1:n
        eps_Aii = θ * θ * diags[i]
        for j in nzrange(A, i)
            row = A.rowval[j]
            val = A.nzval[j]
            if row != i
                if val*val < eps_Aii * diags[row]
                    S.nzval[j] = 0
                end
            end
        end
    end

    dropzeros!(S)

    S.nzval .= abs.(S.nzval)
    scale_cols_by_largest_entry!(S)

    S, S
end


text:
# Main file for solving the full information SAA problem

using Distributions, StatsFuns, StatsBase, JuMP, Gurobi

const gurobi_env = Gurobi.Env()

# read these as inputs from the corresponding "instance file"
modRepNum = parse(ARGS[1])
degNum = parse(ARGS[2])
dataRepNum = parse(ARGS[3])
saaRepNum = parse(ARGS[4])

#*========= MODEL INFORMATION ==========
const caseNum = 1
const paramsFile = "params_case" * string(caseNum) * "_fullinf_saa.jl"
#*======================================

#*========= INCLUDE FILES ====================
include("maxParameters.jl")
include("randomSeeds.jl")

include(paramsFile)
# the next two data files are defined in paramsFile
include(modelFile)
include(dataFile)

include("solveModel.jl")
#*======================================

# set seed for reproducibility
srand(startingSeed)

# directory name for storing results
const baseDirName = "case" * string(caseNum) * "_fullinf_saa/" * "mod_" * string(modRepNum) * "/" * "deg_" * string(degNum) * "/"
const subDirName = baseDirName * "rep_" * string(dataRepNum) * "/"
const subDirName2 = subDirName * "saa_" * string(saaRepNum) * "/"
mkpath(subDirName2)

#*========= PRINT OPTIONS ====================
const storeResults = true

const infoFile = "ddsp.txt"
const modelDataFile = "model_data.txt"

const fullinfSAAObjFile = "fullinf_saa_obj.txt"
const fullinfSAATimeFile = "fullinf_saa_time.txt"
const covRealFile = "covariate_obs.txt"
#*======================================

#*========= GENERATE MODEL PARAMETERS ====================

srand(randomSeeds_models[modRepNum])
numCovariates = 3

covariate_mean, covariate_covMat, coeff_true = generateBaseDemandModel(numCovariates)

#*========= STORE RESULTS ====================
if(storeResults)

	# write details to text file, including some key details about the test instance
	details_file = baseDirName * infoFile
	open(details_file, "w") do f
		write(f,"case number: $caseNum \n")
		write(f,"numResources: $numResources \n")
		write(f,"numCustomers: $numCustomers \n")
		write(f,"model replicate number: $modRepNum \n")
		write(f,"degree: $(degree[degNum]) \n")
		write(f,"demand_errors_scaling: $demand_errors_scaling \n")
		write(f,"numMCScenarios: $numMCScenarios \n\n")

		write(f,"randomSeeds_MC: $randomSeeds_MC \n")
		write(f,"randomSeeds_models: $randomSeeds_models \n")
		write(f,"randomSeeds_data: $randomSeeds_data \n")
	end

	mod_data_file = baseDirName * modelDataFile
	open(mod_data_file, "w") do f
		write(f,"covariate_mean = $covariate_mean \n")
		write(f,"covariate_covMat = $covariate_covMat \n")
		write(f,"coeff_true = $coeff_true \n")
	end

end
#*============================================

#*========= GENERATE DATA REPLICATE ====================

srand(randomSeeds_data[dataRepNum])

covariate_obs = generateCovariateReal(numCovariates,covariate_mean,covariate_covMat)

#*========= STORE RESULTS ====================
if(storeResults)

	cov_real_file = subDirName * covRealFile
	open(cov_real_file, "w") do f
		write(f,"covariate_obs = $covariate_obs")
	end

end
#*============================================

#*========= SOLVE THE FULL-INFORMATION SAA MODEL ====================

tic()

# create scenarios from the true conditional distribution
srand(randomSeeds_MC[saaRepNum])

demand_scen_MC = generateTrueCondScenarios(numMCScenarios,covariate_obs,degree[degNum],coeff_true)

# solve SAA problems to estimate true objective value at this conditioning
~, SAAObj = solveSAAModel(demand_scen_MC)

SAATime = toq()

#*========= STORE RESULTS ====================
if(storeResults)

	# write data to text file
	obj_est_file = subDirName2 * fullinfSAAObjFile
	open(obj_est_file, "w") do f
		write(f,"$SAAObj")
	end

	time_obj_file = subDirName2 * fullinfSAATimeFile
	open(time_obj_file, "w") do f
		write(f,"$SAATime")
	end

end
#*============================================

max_stars_repo_name: rohitkannan/DD-SAA

max_issues_repo_name: rohitkannan/DD-SAA

max_forks_repo_name: rohitkannan/DD-SAA


text:
module DanceWebpack

import Dance.Configuration

export setup

"""
    setup(project_path::String=".")

- Copy sample files (`package.json`, `webpack.config.js`, `src` dir) to `static` dir in project
- Obtain `authors`, `name`, `version` from `Project.toml` and set in `package.json` (optional, values can be blank if no `Project.toml` found)
- Overwrite `html_file` in `webpack.config.js`
"""
function setup(project_path::String=".") :: Nothing
    copied_file_path::String = abspath(joinpath(project_path, "static"))

    if !isfile("dance.jl")
        @error "Please add dance.jl to project root"
    elseif !isdir("settings")
        @error "Please add settings dir to project root"
    else
        if Configuration.populate(project_path)
            settings::Dict{Symbol, Any} = Configuration.Settings
            html_base_file_path::String= settings[:html_base_filename]
            routes_path::String = settings[:routes_filename]

            # Copy Webpack files to project root
            sample_files_dir::String = joinpath(@__DIR__, "../files")
            cp(sample_files_dir, copied_file_path)
            if Sys.iswindows()
                run(`icacls.exe $copied_file_path /reset /T /Q`)
                username::String = read(run(`whoami`), String)
                run(`icacls.exe $copied_file_path /grant $username:F /T /Q`)
            else
                run(`chmod -R 755 $copied_file_path`)
            end

            # Obtain `authors`, `name`, `version` values from project Project.toml
            authors::String = ""
            name::String = ""
            version::String = ""
            try
                for line in readlines(joinpath(copied_file_path, "../Project.toml"))
                    for item in [
                        ["authors", authors],
                        ["name", name],
                        ["version", version]
                    ]
                        if occursin(item[1], line)
                            value::String = strip(split(line, '=')[2])
                            if item[1]=="authors"
                                authors = replace(value, "\\\"" => "*")  # temporary flag quotation escapes: \"
                                authors = replace(authors, "\"" => "\\\"")
                                authors = replace(authors, "*" => "\\\"")  # put back flag quotation escapes: \"
                            else
                                item[2] = strip(value, '"')
                                if item[1]=="name"
                                    name = item[2]

                                    # package.json names cannot contain spaces
                                    if occursin(" ", name)
                                        name = replace(name, " " => "_")
                                    end
                                elseif item[1]=="version"
                                    version = item[2]
                                end
                            end
                        end
                    end
                end
            catch e
                nothing
            end

            # Overwrite `author`, `name`, `version` var tags in package.json
            file_string::String = ""
            package_json_file_path::String = joinpath(copied_file_path, "package.json")
            open(package_json_file_path, "r") do io
                file_string = read(io, String)
            end
            for item in [
                ["author", authors],
                ["name", name],
                ["version", version]
            ]
                idx::String = "\$" * item[1]
                file_string = replace(file_string, "$idx" => item[2])
            end
            open(package_json_file_path, "w") do io
                write(io, file_string)
            end

            # Overwrite `$html_file` in webpack.config.js
            webpack_config_file_path::String = joinpath(copied_file_path, "webpack.config.js")
            open(webpack_config_file_path, "r") do io
                file_string = read(io, String)
            end
            idx = "\$html_file"
            file_string = replace(file_string, "$idx" => "\"../" * html_base_file_path * ".html\"")
            open(webpack_config_file_path, "w") do io
                write(io, file_string)
            end

            # Add Webpack static dir output to project routes.jl file
            open(routes_path * ".jl", "a+") do io
                write(io, "static_dir(\"/static\", \"static/dist\")\n")
            end

            @info "Webpack files added to `static` dir in project.\nPlease cd into `static` dir and add dependencies by installing NodeJS and running `npm install`"
        else
            @error "Error setting up project settings. Please see error above."
        end
    end
end

end


text:
using LightGraphs.Experimental.Traversals
const LET = LightGraphs.Experimental.Traversals
struct DummyTraversalState <: LET.AbstractTraversalState end

@testset "Traversals" begin
    d = DummyTraversalState()
    @test LET.initfn!(d, 1) == true
    @test LET.previsitfn!(d, 1) == true
    @test LET.visitfn!(d, 1, 2) == true
    @test LET.newvisitfn!(d, 1, 2) == true
    @test LET.postvisitfn!(d, 1) == true
    @test LET.postlevelfn!(d) == true

    x = [3,1,2]
    @test sort!(x, 3, 1, LET.NOOPSort, Base.Sort.ForwardOrdering()) == x

    @testset "BFS" begin
        g1= smallgraph(:house)
        dg1 = path_digraph(6); add_edge!(dg1, 2, 6)

        b1 = LET.BFS()
        @test b1 == LET.BFS(LET.NOOPSort)

        b2 = LET.BFS(QuickSort)

        @testset "undirected" begin
            for g in testgraphs(g1)
                v1 = @inferred LET.visited_vertices(g, 3, b1)
                v2 = @inferred LET.visited_vertices(g, 3, b2)
                v3 = @inferred LET.visited_vertices(g, [3], b1)
                @test v1 == v2 == v3 == [3, 1, 4, 5, 2]
                p1 = @inferred LET.parents(g, 2, b1)
                p2 = @inferred LET.parents(g, 2, b2)
                @test p1 == p2 == [2, 0, 1, 2, 4]
                d1 = @inferred LET.distances(g, 1, b1)
                d2 = @inferred LET.distances(g, 1 ,b2)
                d3 = @inferred LET.distances(g, [1], b1)
                @test d1 == d2 == d3 == [0, 1, 1, 2, 2]
                t1 = @inferred LET.tree(g, 2, b1)
                t2 = @inferred LET.tree(g, 2, b2)
                t3 = @inferred LET.tree(p1)
                t4 = @inferred LET.tree(p2)
                @test t1 == t2 == t3 == t4

            end
        end
        @testset "directed" begin
            for dg in testgraphs(dg1)
                v1 = @inferred LET.visited_vertices(dg, 1, b1)
                v2 = @inferred LET.visited_vertices(dg, 1, b2)
                v3 = @inferred LET.visited_vertices(dg, [1], b1)
                @test v1 == v2 == v3 == [1, 2, 3, 6, 4, 5]
                @test (@inferred LET.parents(dg, 2, b1)) == LET.parents(dg, 2, b2) == [0, 0, 2, 3, 4, 2]
                p1 = @inferred LET.parents(dg, 2, b1)
                p2 = @inferred LET.parents(dg, 2, b2)
                p3 = @inferred LET.parents(dg, 6, b1, inneighbors)
                p4 = @inferred LET.parents(dg, 6, b2, inneighbors)
                @test p1 == p2 == [0, 0, 2, 3, 4, 2]
                @test p3 == p4 == [2, 6, 4, 5, 6, 0]
                d1 = @inferred LET.distances(dg, 1, b1)
                d2 = @inferred LET.distances(dg, 1, b2)
                d3 = @inferred LET.distances(dg, [1], b1)
                @test d1 == d2 == d3 == [0, 1, 2, 3, 4, 2]
                t1 = @inferred LET.tree(dg, 2, b1)
                t2 = @inferred LET.tree(dg, 2, b2)
                t3 = @inferred LET.tree(p1)
                t4 = @inferred LET.tree(p2)
                @test t1 == t2 == t3 == t4
                t5 = @inferred LET.tree(dg, 6, b1, inneighbors)
                t6 = @inferred LET.tree(dg, 6, b2, inneighbors)
                t7 = @inferred LET.tree(p3)
                t8 = @inferred LET.tree(p4)
                @test t5 == t6 == t7 == t8
            end
        end
    end
    @testset "DFS" begin
        g1 = binary_tree(3)
        dg1 = DiGraph(g1)
        for e in edges(dg1)
            if src(e) > dst(e)
                rem_edge!(dg1, e)
            end
        end

        d = LET.DFS()
        @testset "undirected" begin
            for g in testgraphs(g1)
                v1 = @inferred LET.visited_vertices(g, 1, d)
                v2 = @inferred LET.visited_vertices(g, [1], d)
                @test v1 == v2 == [1, 2, 4, 5, 3, 6, 7]
                p1 = @inferred LET.parents(g, 2, d)
                @test p1 == [2, 0, 1, 2, 2, 3, 3]

                t1 = @inferred LET.tree(g, 2, d)
                t2 = @inferred LET.tree(p1)
                @test t1 == t2
            end
        end
        @testset "directed" begin
            for dg in testgraphs(dg1)
                dg2 = copy(dg)
                add_edge!(dg2, 6, 1)
                v1 = @inferred LET.visited_vertices(dg, 1, d)
                v2 = @inferred LET.visited_vertices(dg, [1], d)
                @test v1 == v2 == [1, 2, 4, 5, 3, 6, 7]

                p1 = @inferred LET.parents(dg, 2, d)
                @test p1 == [0, 0, 0, 2, 2, 0, 0]

                ts1 = @inferred LET.topological_sort(dg)
                @test ts1 == [1, 3, 7, 6, 2, 5, 4]
                @test_throws LET.CycleError topological_sort(dg2)

                t1 = @inferred LET.tree(dg, 2, d)
                t2 = @inferred LET.tree(p1)
                @test t1 == t2

                @test !LET.is_cyclic(dg1)
                @test LET.is_cyclic(dg2)
            end
        end
    end
end


text:
# julia connection to julius
using Sockets, LightXML
# voice decode event
struct VDevent
	vdscore::Float32 # score
	vdconf::Float32 # confidence
	vdstring::String # Julius interp what was said
end
"""
	juliusProcessID

	check if Julius is running by fetching the PID from the system
	knowledge of PID can be used to force shutdown of process later
"""
function juliusProcessID()
	ps = read(`ps ax`,String)
	pss = split(ps,"\n")
	pidid = 0
	for p in pss
		global pidid
		if contains(p,"mod.jconf")
			pid = split(strip(p))[1]
			pidid = parse(Int32,pid)
			return pidid
		end
	end
	return pidid
end
"""
	juliusGetStatus()

	find out what state Julius is in (listening, )
	Input is TCPSocket
	Output is the message
"""
function juliusGetStatus(cs::TCPSocket)
	write(cs,"STATUS\n")
	sleep(0.5)
	flush(cs)
	while true
		msg = readline(cs)
		if contains(msg,"SYSINFO")
			return msg
		end
	end
end
"""
	juliusDecodeVoice()

	Input is the TCPSocket already open
	Output is the complete decode message
	including input status, time, frames, timing,
	and word decode hypotheses
"""
function juliusDecodeVoice(cs::TCPSocket)
	flush(cs)
	msg = ""
	while true
		lin = readline(cs)
		#println(lin)
		msg *= lin
		if contains(lin,"/RECOGOUT")
			break
		end
	end
	msg = replace(msg,">.<"=>">\n<")
	xml = "<REC>"*msg[2:end]*"</REC>"
	#println(xml)
	return parse_xml(xml)
end
"""
	parse_xml()

	inputs a string in XML form and
	extracts the words as a string, the score and
	the mean confidence for the string as a whole
	return as a struct
"""
function parse_xml(xml::String)
	xdoc = parse_string(xml)
	xroot = root(xdoc)
	xscore = xroot["RECOGOUT"][1]["SHYPO"][1]
	score = parse(Float32,attribute(xscore,"SCORE"))
	xwds = xroot["RECOGOUT"][1]["SHYPO"][1]["WHYPO"][2:end-1]
	confs = [parse(Float32,attribute(a,"CM")) for a in xwds]
	nconfs = length(confs)
	conf = sum(confs) ./ nconfs
	heard = join([attribute(a,"WORD") for a in xwds]," ")
	vdevent = VDevent(score,conf,heard)
	return vdevent
end
"""
	juliusRecogToggle()

	tells the Julius server to ignore the speech it hears (off/false)
	or to resume decoding (on/true)
"""
function juliusRecogToggle(cs,onoff::Bool)
	action = onoff ? "RESUME\n" : "PAUSE\n"
	write(cs,action)
end
#
"""
	announce()

	say something through flite
"""
function announce(msg::String)
	cmd = `padsp flite -voice slt -t $msg`
	run(cmd)
end
#
"""
	tidyupandend()

	perform any needed cleanup tasks
	and terminate
"""
function tidyupandend()
	announce("Ending now. Goodbye")
	println("Ending now. Goodbye")
	exit(0)
end
#
# loop
#
function recogLoop(cs::TCPSocket)
	while true
		juliusRecogToggle(cs,false)
		announce("Ready")
		juliusRecogToggle(cs,true)
		vdevent = juliusDecodeVoice(cs)
		println(vdevent.vdstring)
		if vdevent.vdstring == "zulu zulu"
			tidyupandend()
		end
	end
end
#
# main section
# establish the connection to the server and start decoding
#
Jpid = juliusProcessID()
println(Jpid)
if  Jpid > 0
	local cs
	try
		cs = connect(10500)
	catch e
		println(e)
		announce("julius is not running. stopping now")
		error("cannot connect to socket on 10500")
	end
	announce("Julius is ready")
	println(juliusGetStatus(cs))
	recogLoop(cs)
end

max_stars_repo_name: colbec/voiceDM

max_issues_repo_name: colbec/voiceDM

max_forks_repo_name: colbec/voiceDM


text:
using GeometricEquations
using Test

include("functions.jl")
include("initial_conditions.jl")

@testset "$(rpad("Abstract Problem Interface",80))" begin

end


text:
using GeneralizedGenerated

export rand

EmptyNTtype = NamedTuple{(),Tuple{}} where T<:Tuple

@inline function rand(m::JointDistribution)
    return _rand(m.model, m.args)
end

@inline function rand(m::Model)
    return _rand(m, NamedTuple())
end

@gg function _rand(_m::Model, _args)
    type2model(_m) |> sourceRand() |> loadvals(_args, NamedTuple())
end

@gg function _rand(_m::Model, _args::NamedTuple{()})
    type2model(_m) |> sourceRand()
end

export sourceRand
function sourceRand()
    function(m::Model)

        _m = canonical(m)
        proc(_m, st::Assign)  = :($(st.x) = $(st.rhs))
        proc(_m, st::Sample)  = :($(st.x) = rand($(st.rhs)))
        proc(_m, st::Return)  = :(return $(st.rhs))
        proc(_m, st::LineNumber) = nothing

        vals = map(x -> Expr(:(=), x,x),variables(_m))

        wrap(kernel) = @q begin
            $kernel
            $(Expr(:tuple, vals...))
        end

        buildSource(_m, proc, wrap) |> flatten
    end
end


text:
using DrWatson
@quickactivate "OptimalEstimationProject"
using IndirectTrajOpt

# This script attemplts to find an optimal minimum fuel trajectory with initial
# and final conditions ics and fcs respectively. It does so with an indirect
# approaching, first searching for a good selection for the initial co-states
# corresponding to the minimum energy problem. A homotopy continuation is
# then performed to transition PSO initialized co-states from the minumum
# energy problem to that of mimimum fuel.

# Ref: Hecht, G., Botta, E. 2021. "Co-State Initialization with Particle Swarm
# Optimization for Low-Thrust Mimumum-Fuel Trajectory Optimization," AAS/AIAA
# Astrodynamic Specialist Conference.

# IndirectTrajOpt.jl Julia package and its dependancies can be found at:
#   IndirectTrajOpt.jl      - https://github.com/GrantHecht/IndirectTrajOpt.jl
#   IndirectCoStateInit.jl  - https://github.com/GrantHecht/IndirectCoStateInit.jl
#   IndirectShooting.jl     - https://github.com/GrantHecht/IndirectShooting.jl
#   Heuristics.jl           - https://github.com/GrantHecht/Heuristics.jl

function genRefTrajs(swarmSize, numTrials, tf)
    # Set tspan in days
    tspan = (0.0, tf) # [Days]

    # Grab parameters for scaling
    pTemp = IndirectTrajOpt.initCR3BPIndirectParams("Low Thrust 10 CR3BP")
    μ  = pTemp.crp.μ
    LU = pTemp.crp.LU
    TU = pTemp.crp.TU

    # Set initial and final conditions
    ics = [-0.019488511458668, -0.016033479812051, 0.0,
            8.918881923678198, -4.081793688818725, 0.0,
            1.0]
    fcs = [(1.0 - μ) + 6060.483/LU, 19452.284/LU, -34982.968/LU, 0.082677*(TU/LU), 0.006820*(TU/LU), -0.368434*(TU/LU), 0.0]

    # Homotopy continuation
    N = 20
    hVec = [(j^2 - 1.)/(N^2 - 1.) for j in N:-1:1]

    # Instantiate indirect optimization problem
    prob = IndirectOptimizationProblem("Low Thrust 10 CR3BP", ics, fcs, tspan)

    # Create vector of co-state initializers
    itoVec = [IndirectTrajOptimizer(
        prob;
        initOptimizer    = :PSO,
        initCostFunc     = :WSS,
        numParticles     = swarmSize,
        numSwarms        = 8,
        UBs              = [20., 20., 5., 0.1, 0.1, 0.01, 1.0],
        LBs              = [-20., -20., -5., -0.1, -0.1, -0.01, 0.0],
        iUBs             = [1e-6, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6, 1.0],
        iLBs             = [-1e-6, -1e-6, -1e-6, -1e-6, -1e-6, -1e-6, 1e-6],
        weights          = [10, 10, 10, 1, 1, 1, 1],
        MFD              = 1.0,
        displayInterval  = 5,
        funcTol          = 1e-5,
        minNeighborhoodFraction = 0.05,
        maxStallIters    = 50,
        maxIters         = 5000,
        maxTime          = 1800.0,
        homotopyParamVec = hVec,
        dataFolder       = datadir("numParticles"*string(swarmSize)*"_tf"*replace(string(tf), '.' => 'p')),
        writeData        = true) for i in 1:numTrials];

    # Initialize all co-states
    for i in 1:numTrials
        initialize!(itoVec[i])
        solve!(itoVec[i], factor = 3.0, convergenceAttempts = 4, ftol = 1e-11)
    end
end

max_stars_repo_name: GrantHecht/OptimalEstimationProject


text:
include("../src/FluxFederated.jl")
using .FluxFederated
using Sockets
using Flux # Note: this takes awhile to load

# Listen on localhost port 3000
server = listen(3000)

# Initialize model
model = Chain(Dense(10,5,relu),Dense(5,2),softmax)

println("Initial params\n", params(model), "\n")

# Serve model asynchronously
serve_model(server, model)


text:
# Data types declared in `include/blockmat.h`
# Small type names are the C types
# Capitalized are the corresponding Julia types

# Utils

function fptr(x::Vector{T}) where T
    # CSDP starts indexing at 1 so we need to do "- sizeof(T)"
    pointer(x) - sizeof(T)
end

function ptr(x::X) where X
    Base.reinterpret(Base.Ptr{X}, Base.pointer_from_objref(x))
end
export fptr, ptr

function mywrap(X::blockmatrix)
    BlockMatrix(X)
end

function _unsafe_wrap(A, x, n, own::Bool)
    Base.unsafe_wrap(A, x, n, own=own)
end

function mywrap(x::Ptr{T}, len) where T
    # I give false to unsafe_wrap to specify that Julia do not own the array so it should not free it
    # because the pointer it has has an offset
    y = _unsafe_wrap(Array, x + sizeof(T), len, false)
    # fptr takes care of this offset
    #finalizer(s -> Libc.free(fptr(s)), y)
    y
end

# The problem is
# max ⟨C, X⟩
#   ⟨A_i, X⟩ = b_i  ∀ i
#         X  ⪰ 0

# CSDP does not have any C functions manipulate the values of the structures
# so we need to create the structures in Julia.
# Since the arrays are stored in Julia, we need to make sure they are not garbage collected
# so we put both the C structure (named csdp) and the allocated arrays in a wrapper with the same name as the C structure but capitalized

# In CSDP, the matrix C is represented as a blockmatrix
# which is a block diagonal matrix where each block is a symmetric dense (or diagonal) matrix of Cdouble
# * blockdatarec contains a vector containing the values (Cdouble) (1-indexed)
# * blockrec contains the blockdatarec, the category (dense or diagonal), and the block size
# * blockmatrix contains the number of blocks and a vector containing the blocks (blockrec) (1-indexed)

# blockrec
mutable struct BlockRec <: AbstractMatrix{Cdouble}
    _blockdatarec::Vector{Cdouble}
    csdp::blockrec
end
function BlockRec(a::Vector{Cdouble}, n::Int)
    if n > 0
        # /!\ the matrix is 0-indexed -> pointer
        BlockRec(a, blockrec(blockdatarec(pointer(a)), MATRIX, csdpshort(n)))
    else
        # /!\ the diagonal matrix is 1-indexed -> fptr
        BlockRec(a, blockrec(blockdatarec(fptr(a)), DIAG, csdpshort(-n)))
    end
end
BlockRec(a::Vector, n::Int) = BlockRec(Vector{Cdouble}(a), n)
function BlockRec(A::Matrix)
    return BlockRec(reshape(A, length(A)), LinearAlgebra.checksquare(A))
end
function BlockRec(A::Diagonal)
    a = Vector{Cdouble}(diag(A))
    BlockRec(a, -length(a))
end
function blockreczeros(n)
    if n > 0
        BlockRec(zeros(Cdouble, n^2), n)
    else
        BlockRec(zeros(Cdouble, -n), n)
    end
end

function Base.size(A::blockrec)
    n = A.blocksize
    (n, n)
end
Base.size(A::BlockRec) = size(A.csdp)
function Base.getindex(A::BlockRec, i, j)
    n = A.csdp.blocksize
    if A.csdp.blockcategory == MATRIX
        A._blockdatarec[i+(j-1)*n]
    elseif  A.csdp.blockcategory == DIAG
        if i == j
            A._blockdatarec[i]
        else
            0
        end
    else
        error("Invalid category")
    end
end
function Base.setindex!(A::BlockRec, v, i, j)
    n = A.csdp.blocksize
    if A.csdp.blockcategory == MATRIX
        A._blockdatarec[i+(j-1)*n] = v
        A._blockdatarec[j+(i-1)*n] = v
    elseif  A.csdp.blockcategory == DIAG
        if i == j
            A._blockdatarec[i] = v
        else
            error("Cannot set off-diagonal entry of diagonal matrix")
        end
    else
        error("Invalid category")
    end
end

# blockmatrix
mutable struct BlockMatrix <: AbstractBlockMatrix{Cdouble}
    jblocks::Vector{BlockRec}
    blocks::Vector{blockrec}
    csdp::blockmatrix
end

function BlockMatrix(jblocks::AbstractVector{BlockRec})
    blocks = map(block->block.csdp, jblocks)
    csdp = blockmatrix(length(blocks), fptr(blocks))
    BlockMatrix(jblocks, blocks, csdp)
end
BlockMatrix(As::AbstractVector) = BlockMatrix(map(BlockRec, As))
BlockMatrix(As::AbstractMatrix...) = BlockMatrix(collect(As))
blockmatzeros(blkdims) = BlockMatrix(map(blockreczeros, blkdims))

function BlockMatrix(csdp::blockmatrix)
    # I give false so that Julia does not try to free it
    blocks = _unsafe_wrap(Array, csdp.blocks + sizeof(blockrec), csdp.nblocks, false)
    jblocks = map(blocks) do csdp
        let n = csdp.blocksize, c = csdp.blockcategory, d = csdp.data._blockdatarec
            if c == MATRIX
                _blockdatarec = _unsafe_wrap(Array, d, n^2, false)
            elseif c == DIAG
                _blockdatarec = _unsafe_wrap(Array, d + sizeof(Cdouble), n, false)
            else
                error("Unknown block category $(c)")
            end
            BlockRec(_blockdatarec, csdp)
        end
    end
    BlockMatrix(jblocks, blocks, csdp)
end

## Old function for that...
#   function Blockmatrix(X::CSDP.blockmatrix)
#       bs = pointer_to_array(X.blocks + sizeof(CSDP.blockrec), X.nblocks)
#       Bs = map(bs) do b
#           let s = b.blocksize, c = b.blockcategory, d = b.data._blockdatarec
#               if b.blockcategory == CSDP.MATRIX
#                   pointer_to_array(d, (s, s))
#               elseif b.blockcategory == CSDP.DIAG
#                   diagm(pointer_to_array(d + sizeof(Cdouble), s))
#               else
#                   error("Unknown block category $(b.blockcategory)")
#               end
#           end
#       end
#       Blockmatrix(Bs, bs)
#   end

# In CSDP, the matrices A_i are represented as a constraintmatrix
# which is a block diagonal matrix where each block is a symmetric sparse matrix of Cdouble
# * sparseblock contains a sparse description of the entries of a block
# * constraintmatrix contains a linked list of the blocks

mutable struct SparseBlock <: AbstractMatrix{Cdouble}
    i::Vector{csdpshort}
    j::Vector{csdpshort}
    v::Vector{Cdouble}
    n::CSDP_INT
    csdp::sparseblock
    function SparseBlock(i::Vector{csdpshort}, j::Vector{csdpshort}, v::Vector{Cdouble}, n::Integer, csdp)
        new(i, j, v, n, csdp)
    end
end

function SparseBlock(i::Vector{csdpshort}, j::Vector{csdpshort}, v::Vector{Cdouble}, n::Integer)
    @assert length(i) == length(j) == length(v)
    block = sparseblock(C_NULL,    # next
                        C_NULL,    # nextbyblock
                        fptr(v),   # entries
                        fptr(i),   # iindices
                        fptr(j),   # jindices
                        length(i), # numentries
                        0,         # blocknum
                        n,         # blocksize
                        0,         # constraintnum
                        1)         # issparse
    SparseBlock(i, j, v, n, block)
end

SparseBlock(A::SparseBlock) = A
function SparseBlock(A::SparseMatrixCSC{Cdouble})
    n = LinearAlgebra.checksquare(A)
    nn = nnz(A)
    I = csdpshort[]
    J = csdpshort[]
    V = Cdouble[]
    vals = nonzeros(A)
    rows = rowvals(A)
    for col = 1:n
        for j in nzrange(A, col)
            row = rows[j]
            if row <= col
                push!(I, row)
                push!(J, col)
                push!(V, vals[j])
            end
        end
    end
    SparseBlock(I, J, V, n)
end
SparseBlock(A::AbstractMatrix{Cdouble}) = SparseBlock(SparseMatrixCSC{Cdouble, CSDP_INT}(A))
SparseBlock(A::AbstractMatrix) = SparseBlock(map(Cdouble, A))
Base.convert(::Type{SparseBlock}, A::AbstractMatrix) = SparseBlock(A)

function sparseblockzeros(n)
    SparseBlock(csdpshort[], csdpshort[], Cdouble[], abs(n))
end

function Base.size(A::SparseBlock)
    (A.n, A.n)
end
function findindices(A::SparseBlock, i, j)
    if i > A.n || j > A.n || i <= 0 || j <= 0
        error("Invalid indices")
    end
    for k in 1:length(A.i)
        if A.i[k] == i && A.j[k] == j
            return k
        end
    end
    return 0
end
function Base.getindex(A::SparseBlock, i, j)
    k = findindices(A, i, j)
    if k == 0
        0
    else
        A.v[k]
    end
end
function Base.setindex!(A::SparseBlock, v, i, j)
    k = findindices(A, i, j)
    if k == 0
        push!(A.i, i)
        push!(A.j, j)
        push!(A.v, v)
        @assert A.csdp.numentries + 1 == length(A.i)
        A.csdp.numentries = length(A.i)
        # If push! has reallocated it, we need to change the pointer
        A.csdp.entries = fptr(A.v)
        A.csdp.iindices = fptr(A.i)
        A.csdp.jindices = fptr(A.j)
    else
        A.v[k] = v
    end
end

mutable struct ConstraintMatrix <: AbstractBlockMatrix{Cdouble}
    jblocks::Vector{SparseBlock}
    csdp::constraintmatrix
end

function ConstraintMatrix(constr, jblocks::AbstractVector{SparseBlock})
    if isempty(jblocks)
        error("No variable")
    end
    next = C_NULL
    for blocknum in length(jblocks):-1:1
        jblock = jblocks[blocknum]
        jblock.csdp.next = next
        jblock.csdp.blocknum = blocknum
        jblock.csdp.constraintnum = constr
        next = pointer_from_objref(jblock.csdp)
        #next = Base.unsafe_convert(Ptr{sparseblock}, Ref(block))
    end
    csdp = constraintmatrix(Ptr{sparseblock}(next))
    ConstraintMatrix(jblocks, csdp)
end
ConstraintMatrix(i, bs::AbstractMatrix...) = ConstraintMatrix(i, SparseBlock[b for b in bs])
constrmatzeros(i, blkdims) = ConstraintMatrix(i, map(sparseblockzeros, blkdims))

# I need sparseblock to be immutable for this function to work
# but this is kind of annoying since I want to modify its entries in setindex!(::SparseBlock, ...)
function ConstraintMatrix(csdp::constraintmatrix, k::Integer)
    # I take care to free the blocks array when necessary in mywrap since CSDP won't take care of that (see the code of read_prob)
    blocks = _unsafe_wrap(Array, csdp.blocks, k, true) # FIXME this is a linked list, not an array...
    jblocks = map(blocks) do csdp
        ne = csdp.numentries
        i = mywrap(csdp.iindices, ne)
        j = mywrap(csdp.jindices, ne)
        v = mywrap(csdp.entries, ne)
        SparseBlock(i, j, v, csdp.blocksize, csdp)
    end
    ConstraintMatrix(jblocks, csdp)
end

# Needed by MPB_wrapper
function Base.getindex(A::Union{blockmatrix, BlockMatrix, ConstraintMatrix}, i::Integer)
    block(A, i)
end

block(A::blockmatrix, i::Integer) = getblockrec(A, i)
function block(A::Union{BlockMatrix, ConstraintMatrix}, i::Integer)
    A.jblocks[i]
end
nblocks(A::blockmatrix) = A.nblocks
nblocks(A::Union{BlockMatrix, ConstraintMatrix}) = length(A.jblocks)

export BlockMatrix, ConstraintMatrix

"""Solver status"""
mutable struct Csdp
    n::CSDP_INT
    k::CSDP_INT
    X::BlockMatrix
    y::Vector{Cdouble}
    constant_offset::Cdouble
    constraints::Vector{SparseBlock}
    pobj::Cdouble
    dobj::Cdouble
end


text:
import AlphaZero.GI

using StaticArrays
using Crayons

const Cell = UInt8
const EMPTY = 0x00
const ARROW = 0x03

tmp = [0 1 0 0 1 0;
       1 0 0 0 0 1;
       0 0 0 0 0 0;
       0 0 0 0 0 0;
       2 0 0 0 0 2;
       0 2 0 0 2 0]

const NUM_COLS = size(tmp,1)
const NUM_ROWS = size(tmp,2)
const NUM_CELLS = length(tmp)
const Board = SMatrix{NUM_COLS, NUM_ROWS, Cell, NUM_CELLS}

const Player = UInt8
const WHITE = 0x01
const BLACK = 0x02

other(p::Player) = 0x03 - p

const INITIAL_BOARD = SMatrix{NUM_COLS,NUM_ROWS}(Player.(tmp))
const INITIAL_STATE = (board=INITIAL_BOARD, curplayer=WHITE, phase=0x00,selected=0x00)

struct GameSpec <: GI.AbstractGameSpec end

mutable struct GameEnv <: GI.AbstractGameEnv
  board :: Board
  curplayer :: Player
  finished :: Bool
  winner :: Player
  phase :: UInt8
  selected :: UInt8
  #selected works as long as size of board is not too big
end

GI.spec(::GameEnv) = GameSpec()

function GI.init(::GameSpec)
  board = INITIAL_STATE.board
  curplayer = INITIAL_STATE.curplayer
  finished = false
  winner = 0x00
  phase = INITIAL_STATE.phase
  selected = INITIAL_STATE.selected
  return GameEnv(board, curplayer, finished, winner, phase, selected)
end

GI.two_players(::GameSpec) = true

GI.actions(::GameSpec) = collect(1:NUM_CELLS)

flip_cell_color(c::Cell) = (0 < c < 3) ? other(c) : c

flip_colors(board::Board) = flip_cell_color.(board)

function GI.vectorize_state(::GameSpec, state)
  tmp = zeros(Float32,NUM_ROWS,NUM_COLS,6)
  board = state.curplayer == WHITE ? state.board : flip_colors(state.board)
  for i = 1:3
    tmp[:,:,i] = board .== i
  end
  if state.phase == 0
    tmp[:,:,4] = tmp[:,:,1]
  else
    tmp[:,:,state.phase+4] = get_queen_moves(board,state.selected)
  end
  return tmp
end

function GI.set_state!(g::GameEnv, state)
  g.board = state.board
  g.curplayer = state.curplayer
  g.phase = state.phase
  g.selected = state.selected
  if any(GI.actions_mask(g))
    g.winner = 0x00
    g.finished = false
  else
    g.winner = other(g.curplayer)
    g.finished = true
  end
  return
end

GI.current_state(g::GameEnv) =
  (board=g.board, curplayer=g.curplayer,phase=g.phase,selected=g.selected)

# TODO: maybe MCTS should make the copy itself. The performance cost should not be great
# and it would probably avoid people a lot of pain.

GI.game_terminated(g::GameEnv) = g.finished

GI.white_playing(g::GameEnv) = g.curplayer == WHITE

function get_queen_moves(b::Board,loc_linear)
  loc = CartesianIndex(xy_of_pos(loc_linear))
  valid_moves = falses(size(b))
  for dir∈CartesianIndex.([(-1,-1),(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1)])
    for i=1:max(size(b)...)
      tmploc = loc + i * dir
      if checkbounds(Bool,b,tmploc) && b[tmploc]==0
        valid_moves[tmploc] = true
      else
        break
      end
    end
  end
  return valid_moves[:]
end

function GI.actions_mask(g::GameEnv)
  if g.phase == 0
    return (g.board .== g.curplayer)[:]
  else
    return get_queen_moves(g.board,g.selected)
  end
end

function GI.play!(g::GameEnv, a)
  if g.phase == 0
    g.selected = a
    if !any(get_queen_moves(g.board,a))
      g.finished = true
      g.winner = other(g.curplayer)
    end
  elseif g.phase == 1
    g.board = Base.setindex(g.board,0,g.selected)
    g.board = Base.setindex(g.board,g.curplayer,a)
    g.selected = a
  else
    g.board = Base.setindex(g.board,3,a)
    g.curplayer = other(g.curplayer)
    g.selected = 0x00
  end
  g.phase = (g.phase+1)%3
end

function GI.white_reward(g::GameEnv)
  g.winner == WHITE && (return  1.)
  g.winner == BLACK && (return -1.)
  return 0.
end

pos_of_xy((x, y)) = (y - 1) * NUM_COLS + (x - 1) + 1

xy_of_pos(pos) = ((pos - 1) % NUM_COLS + 1, (pos - 1) ÷ NUM_COLS + 1)

function generate_dihedral_symmetries()
  N = NUM_COLS
  rot((x, y)) = (y, N - x + 1) # 90° rotation
  flip((x, y)) = (x, N - y + 1) # flip along vertical axis
  ap(f) = p -> pos_of_xy(f(xy_of_pos(p)))
  sym(f) = map(ap(f), collect(1:NUM_CELLS))
  rot2 = rot ∘ rot
  rot3 = rot2 ∘ rot
  return [
    sym(rot), sym(rot2), sym(rot3),
    sym(flip), sym(flip ∘ rot), sym(flip ∘ rot2), sym(flip ∘ rot3)]
end

const SYMMETRIES = generate_dihedral_symmetries()

function GI.symmetries(::GameSpec, s)
  return [
    ((board=Board(s.board[sym]), curplayer=s.curplayer,phase=s.phase,selected=s.selected), sym)
    for sym in SYMMETRIES]
end

#####
##### Interface for interactive exploratory tools
#####

player_color(p) = p == WHITE ? crayon"light_red" : crayon"light_blue"
player_name(p)  = p == WHITE ? "Red" : "Blue"
function cell_mark(c)
  if c == EMPTY
    return "."
  elseif c == ARROW
    return "x"
  else
    return "o"
  end
end

function cell_color(c;selected=false)
  if selected
    return crayon"light_green"
  elseif c == ARROW || c == EMPTY
    return crayon""
  else
    return player_color(c)
  end
end

function GI.render(g::GameEnv; with_position_names=true, botmargin=true)
  pname = player_name(g.curplayer)
  pcol = player_color(g.curplayer)
  print(pcol, pname)
  if g.phase == 0
    print(" selects:")
  elseif g.phase == 1
    print(" moves:")
  else
    print(" shoots")
  end
  print(crayon"reset", "\n\n")
  # Print legend
  print(" ")
  for col in 1:NUM_COLS
    print(col, " ")
  end
  print("\n")
  # Print board
  for row in 1:NUM_ROWS
    print(row," ")
    for col in 1:NUM_COLS
      x = g.board[row,col]
      selected = CartesianIndex(row,col) == g.selected
      print(cell_color(x;selected=selected), cell_mark(x), crayon"reset", " ")
    end
    print("\n")
  end
  botmargin && print("\n")
end

function GI.action_string(::GameSpec, a)
  return string(a.I[1])*" "*string(a.I[2])
end

function GI.parse_action(::GameSpec, str)
  try
    x,y = parse(Int, str[1]),parse(Int,str[3])
    ((1 <= x,y <= NUM_COLS)|>all) ? CartesianIndex(x,y) : nothing
  catch
    nothing
  end
end

#=
function GI.read_state(::GameSpec)
  board = Array(INITIAL_BOARD)
  try
    for col in 1:NUM_COLS
      input = readline()
      for (row, c) in enumerate(input)
        c = lowercase(c)
        if c ∈ ['o', 'w', '1']
          board[col, row] = WHITE
        elseif c ∈ ['x', 'b', '2']
          board[col, row] = BLACK
        end
      end
    end
    nw = count(==(WHITE), board)
    nb = count(==(BLACK), board)
    if nw == nb
      curplayer = WHITE
    elseif nw == nb + 1
      curplayer = BLACK
    else
      return nothing
    end
    return (board=board, curplayer=curplayer)
  catch e
    return nothing
  end
end

=#


text: true && false && true

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
function show_heading(io::IO, requirer)
    print(io, "INFO: POMDPLinter requirements for ")
    printstyled(io, handle_method(requirer), color=:blue)
    println(io, " and dependencies. ([✔] = implemented correctly; [X] = not implemented; [?] = could not determine)")
end

function show_requirer(io::IO, r::AbstractRequirementSet)
    print(io, "For ")
    printstyled(io, "$(handle_method(r.requirer))", color=:blue)
    if r.parent == nothing
        println(io, ":")
    else
        println(io, " (in $(handle_method(r.parent))):")
    end
end

function show_checked_list(io::IO, cl::AbstractVector{T}) where T <: Tuple
    for item in cl
        if ismissing(first(item))
            printstyled(io, "  [?] $(format_method(item[2], item[3]))", color=:yellow)
            println(io)
        elseif first(item) == true
            printstyled(io, "  [✔] $(format_method(item[2], item[3]))", color=:green)
            println(io)
        else
            @assert first(item) == false
            printstyled(io, "  [X] $(format_method(item[2], item[3]))", color=:red)
            println(io)
        end
    end
end

function show_incomplete(io, r::RequirementSet)
    @assert r.exception != nothing
    extype = typeof(r.exception)
    printstyled(io, "  WARNING: Some requirements may not be shown because a $(extype) was thrown.", color=:yellow)
    println(io)
end

handle_method(str::Any) = string(str)
handle_method(str::Req) = format_method(str...)
short_method(str::Any) = string(str)
short_method(str::Req) = string(first(str))

function format_method(f::Function, argtypes::TupleType; module_names=false, color=nothing)
    fname = f
    typenames = argtypes.parameters
    if !module_names
        fname = typeof(f).name.mt.name
        mless_typenames = []
        for t in argtypes.parameters
            if isa(t, Union)
                str = "Union{"
                for (i, tt) in enumerate(fieldnames(typeof(t)))
                    str = string(str, getfield(t, tt), i<length(fieldnames(typeof(t))) ? ',' : '}')
                end
                push!(mless_typenames, str)
            elseif isa(t, UnionAll)
                push!(mless_typenames, string(t))
            else
                str = string(t.name.name)
                if !isempty(t.parameters)
                    ps = map(t.parameters) do p
                        if p isa Symbol
                            return Meta.quot(p)
                        else
                            return p
                        end
                    end
                    str = string(str, "{$(ps...)}")
                end
                push!(mless_typenames, str)
            end
        end
        typenames = mless_typenames
    end
    str = "$fname("
    len = length(typenames)
    for (i, t) in enumerate(typenames)
        str = string(str, "::$t")
        if i < len
            str = string(str, ", ")
        end
    end
    str = string(str, ")")
end

logger_context(::AbstractLogger) = IOContext()
logger_context(l::ConsoleLogger) = IOContext(l.stream)
logger_context(l::SimpleLogger) = IOContext(l.stream)
logger_context() = logger_context(current_logger())

"""
Return a String with the req checked if it is implemented.
"""
function schecked(req::Req; kwargs...)
    checked_list = [(implemented(req), req...)]
    sprint(show_checked_list, checked_list; kwargs...)
end