text:
_vega_model(m::GeomBMModel) = GeomBMModel(m.core, Dual(m.volatility, 1))

_delta_model(m::CoreModel) = CoreModel(Dual(m.startprice, 1), m.yieldcurve, m.carrycurve)
_delta_model(m::GeomBMModel) = GeomBMModel(_delta_model(m.core), m.volatility)

_gamma_model(m::CoreModel) = CoreModel(Dual(Dual(m.startprice, 1), 1), m.yieldcurve, m.carrycurve)
_gamma_model(m::GeomBMModel) = GeomBMModel(_gamma_model(m.core), m.volatility)

_rho_model(y::ConstantContinuousYieldCurve) =
    ConstantContinuousYieldCurve(y.dc, Dual(y.rate,1), y.reference_date)
_rho_model(m::CoreModel) =
    CoreModel(m.startprice, _rho_model(m.yieldcurve), m.carrycurve)
_rho_model(m::GeomBMModel) = GeomBMModel(_rho_model(m.core), m.volatility)

"""
    vega(m::GeomBMModel, c::Contract)
    vega(m::GeomBMModel, c::Contract, MonteCarloModel, npaths::Integer)

Calculate the vega Greek: the first derivative of the option price
with respect to volatility for the contract `c` under the model `m`

Inputs
======

* `m::GeomBMModel` - The model under which the contract is being valued
* `c::Contract` - The contract whose Greek is desired.
* `npaths::Integer` - Number of Monte Carlo paths to use.
"""
function vega(m::GeomBMModel, c::Contract, args...; autodiff=true)
    val = value(_vega_model(m), c, args...)
    extract_derivative(val)
end

"""
    delta(m::GeomBMModel, c::Contract)
    delta(m::GeomBMModel, c::Contract, MonteCarloModel, npaths::Integer)

Calculate the delta Greek: the first derivative of the option price with respect to the
spot price of the underlying for the contract `c` under the model `m`

Inputs
======

* `m::GeomBMModel` - The model under which the contract is being valued
* `c::Contract` - The contract whose greek is desired.
* `npaths::Integer` - Number of Monte Carlo paths to use.
"""
function delta(m::GeomBMModel, c::Contract, args...; autodiff=true)
    val = value(_delta_model(m), c, args...)
    deriv = extract_derivative(val)
    currval(deriv)
end

"""
    rho(m::GeomBMModel, c::Contract)
    rho(m::GeomBMModel, c::Contract, MonteCarloModel, npaths::Integer)

Calculate the rho Greek: the first derivative of the option price with respect to the
interest rate, for the contract `c` under the model `m`

Inputs
======

* `m::GeomBMModel` - The model under which the contract is being valued
* `c::Contract` - The contract whose greek is desired.
* `npaths::Integer` - Number of Monte Carlo paths to use.
"""
function rho(m::GeomBMModel, c::Contract, args...; autodiff=true)
    val = value(_rho_model(m), c, args...)
    deriv = extract_derivative(val)
end

"""
    gamma(m::GeomBMModel, c::Contract)
    gamma(m::GeomBMModel, c::Contract, MonteCarloModel, npaths::Integer)

Calculate the gamma Greek: the second derivative of the option price
with respect to the price of the underlying for the contract `c` under the model `m`

Inputs
======

* `m::GeomBMModel` - The model under which the contract is being valued
* `c::Contract` - The contract whose greek is desired.
* `npaths::Integer` - Number of Monte Carlo paths to use.
"""
function gamma(m::GeomBMModel, c::Contract, args...; autodiff=true)
    val = value(_gamma_model(m), c, args...)
    deriv = extract_derivative(extract_derivative(val))
    currval(deriv)
end

"""
    greek(m::GeomBMModel, c::Contract; metric = :vol, n = 1)

Calculate the nth derivative of the option price with respect to different metrics.

Inputs
======

* `m::GeomBMModel` - The model under which the contract is being valued
* `c::Contract` - The contract whose greek is desired.
* `metric::Symbol` (optional) - The metric with respect to which the derivative
is evaluated. (default = :vol, other options - :strike, :interest)
* `n::Int` - (optional) The derivative desired. (default = 1)

"""
function greek(m::GeomBMModel, c::Contract; metric = :vol, n = 1)
    if metric == :vol
        vol = ndual(m.volatility, n)
        gbmm = GeomBMModel(m.core, vol)
        val = value(gbmm, c)
        return nderiv(val, n)
    elseif metric == :strike
        s = ndual(m.core.startprice, n)
        vol = m.volatility
        gbmm = GeomBMModel(m.core.carrycurve.reference_date, s,
                            m.core.yieldcurve.rate, m.core.carrycurve.rate, vol)
        val = value(gbmm, c)
        return currval(nderiv(val, n))
    elseif metric == :interest
        r = ndual(m.core.yieldcurve.rate, 5)
        vol = m.volatility
        gbmm = GeomBMModel(m.core.carrycurve.reference_date, m.core.startprice,
                            r, m.core.carrycurve.rate, vol)
        val = value(gbmm, c)
        return nderiv(val,n)
    end
end

"""
    `greek(m::MonteCarloModel, c::Contract; metric = :vol, n = 1)`

    Calculate the nth derivative of the option price with respect to different
    metrics.

Inputs
======

* `m::MonteCarloModel` - The model under which the contract is being valued
* `c::Contract` - The contract whose greek is desired.
* `metric::Symbol` (optional) - The metric with respect to which the derivative
is evaluated. (default = :vol, other options - :strike, :interest)
* `n::Int` - (optional) The derivative desired. (default = 1)

"""
function greek(m::MonteCarloModel, c::Contract; metric = :vol, n = 1)
    if metric == :vol
        vol = ndual(m.volatility, n)
        gbmm = GeomBMModel(m.core, vol)
        mc = montecarlo(gbmm, m.dates, size(m.paths,1))
        val = value(mc, c)
        return nderiv(val, n)
    elseif metric == :strike
        s = ndual(m.core.startprice, n)
        vol = m.volatility
        gbmm = GeomBMModel(m.core.carrycurve.reference_date, s,
                            m.core.yieldcurve.rate, m.core.carrycurve.rate, vol)
        mc = montecarlo(gbmm, m.dates, size(m.paths,1))
        val = value(mc, c)
        return currval(nderiv(val, n))
    elseif metric == :interest
        r = ndual(m.core.yieldcurve.rate, 5)
        vol = m.volatility
        gbmm = GeomBMModel(m.core.carrycurve.reference_date, m.core.startprice,
                            r, m.core.carrycurve.rate, vol)
        mc = montecarlo(gbmm, m.dates, size(m.paths,1))
        val = value(mc, c)
        return nderiv(val, n)
    end
end

function ndual(val, n::Int)
    if n == 1
        return Dual(val, 1)
    else
        return Dual(ndual(val, n-1), 1)
    end
end

function nderiv(x, n::Int)
    if n == 1
        return extract_derivative(x)
    else
        return extract_derivative(nderiv(x, n-1))
    end
end

function delta(m::MonteCarloModel{C,D,T}, args...) where {C<:CoreForwardModel,D,T}
    x = [forwardprice(m.core, SingleStock(), dt) for dt in m.dates]
    ref = m.paths ./ x'
    ForwardDiff.gradient(x) do x
        mm = MonteCarloModel(m.core, m.dates, ref .* x')
        value(mm, args...)
    end
end


text:
"""
A container for displaying widgets in a grid.
"""
mutable struct Grid <: Container{GtkGrid}
    widget::GtkWidget
    function Grid(; homogeneous=false, spacing=0, props...)
        self = new(@widget GtkGrid())
        self["column-homogeneous"] = homogeneous
        self["row-homogeneous"] = homogeneous
        if spacing isa Real
            self["row-spacing"] = spacing
            self["column-spacing"] = spacing
        elseif spacing isa Tuple{Real, Real}
            self["row-spacing"] = spacing[1]
            self["column-spacing"] = spacing[2]
        end
        return self
    end
end

function attach!(grid::Grid, child::Widget, x::Int, y::Int, w::Int, h::Int)
    grid[x: x + (w - 1), y:y + (h - 1)] = child
end

function attach!(grid::Grid, sibling::Widget, child::Widget, side::Symbol, w::Int, h::Int)
    gtk_grid_attach_next_to(
        gwidget(grid),
        gwidget(sibling),
        gwidget(child),
        side,
        w, h
    )
end

"""
Represents a grid element that may span multiple rows and/or columns.

You can use the [`span`](@ref) function to create a empty `GridElement`.

Calling a empty `GridElement` as a function and passing a `Widget` as
argument will transform that widget into a `GridElement`, making it
span multiple cells in a [`Grid`](@ref) container.

## Usage:

```julia
grid = Grid([
    # The Button will span 3 columns
    [Button("Example 1") |> span(cols=3)],
    # The Button will span 2 rows
    [Button("Example 2") |> span(rows=2)],
    # The Button will span an 2 rows and 3 columns
    [Button("Example 3") |> span(2, 3)],
])
```
"""
struct GridElement <: Widget{GtkWidget}
    widget::Widget
    w::Int
    h::Int
    GridElement() = new(NullContainer(), 1, 1)
    GridElement(w::Int, h::Int) = new(NullContainer(), w, h)
    GridElement(widget::Widget, w::Int, h::Int) = new(widget, w, h)
    GridElement(e::GridElement) = new(e.widget, e.w, e.h)
    GridElement(widget::Widget) = new(widget, 1, 1)
end

GridElement(::typeof(-)) = GridElement()

gwidget(this::GridElement) = gwidget(this.widget)

(cell::GridElement)(widget::Widget) = GridElement(widget, cell.w, cell.h)
"""
        span(; rows = 1, cols = 1) -> GridElement

Creates a empty [`GridElement`](@ref) that spans multiple `rows` and `cols`.
"""
span(; rows::Int = 1, cols::Int = 1) = GridElement(cols, rows)
"""
        span(rows, cols) -> GridElement

Creates a empty [`GridElement`](@ref) that spans multiple `rows` and `cols`.
"""
span(cols::Int, rows::Int) = GridElement(cols, rows)
"""
        rspan(rows) -> GridElement

Creates a empty [`GridElement`](@ref) that spans multiple `rows`
"""
rspan(rows::Int) = span(1, rows)
"""
        cspan(columns) -> GridElement

Creates a empty [`GridElement`](@ref) that spans multiple `columns`
"""
cspan(columns::Int) = span(columns, 1)

Base.isempty(e::GridElement) = e.widget isa NullContainer

const GridCell = Union{<:Widget, Tuple{}}

"""
        Grid(rows::Vector; props...) -> Grid

Creates a Grid and adds to it each element in `rows`.

The `rows` argument must be a list of `Vector{<:Widget}`.
The position of each element in `rows` determines its
position within the grid.

You can also use [`GridElement`](@ref) to make widgets span
multiple cells within the grid.

## Example:

```julia
grid = Grid(
    Label("A"), Label("B"),     |,
    Label("C"), Label("D"),     |,
    Label("E") |> span(cols=2), |,
)
```
"""
function Grid(cells::Union{<:Widget, typeof(-), typeof(|)}...; props...)
    self = Grid(; props...)
    rows = splitv(x -> x == |, collect(cells))
    for (j, row) in enumerate(rows)
        for (i, cell) in enumerate(row)
            if cell isa GridElement && !isempty(cell)
                w, h = cell.w, cell.h
                x, y = i:i + w - 1, j: j + h - 1
                self[x, y] = cell.widget
            elseif cell isa Widget
                self[i, j] = cell
            else
                # ignore
            end
        end
    end
    return self
end

ext: jl


text:
include("tracer_basic.jl")
include("tracer_provider.jl")


text:
using LinearAlgebra, Statistics, PosDefManifold, Test,
      Diagonalizations

# for more tests see the `Examples` folder

const err=1e-6

# compare two matrices
function compare(a::AbstractArray, b::AbstractArray)
     @test norm(a-b)/length(a) < err
end

# check that cov(F'XX'F)=D arranging terms for the two possible values of `dims`
function compare(d::Int, X::AbstractArray, F::AbstractArray, D::AbstractArray)
    d==1 ? Y=X*F : Y=F'*X
    d==1 ? C=((Y'*Y)/size(Y, 1))[1:size(D, 1), 1:size(D, 2)] :
           C=((Y*Y')/size(Y, 2))[1:size(D, 1), 1:size(D, 2)]
    compare(C, D)
end

# check that D is approximately diagonal
isApproxD(D)=norm(D-Diagonal(D))<err

####  Create data for testing the case k=1, m>1
# `t` is the number of samples, e.g.,
    # in Tangent space transfer learning is the number of trials,
# `m` is the number of datasets, e.g., subjects
# `n` is the number of variables, e.g.,
    # in Tangent space transfer learning is the number of elements of the tangent vectors
# `noise` must be smaller than 1.0. The smaller the noise, the more data are correlated
function getData(t, m, n, noise)
    # create m identical data matrices and rotate them by different
    # random orthogonal matrices V_1,...,V_m
    𝐕=[randU(n) for i=1:m] # random orthogonal matrices
    X=randn(n, t)  # data common to all subjects
    # each subject has this common part plus a random part
    𝐗=[𝐕[i]'*((1-noise)*X + noise*randn(n, t)) for i=1:m]
    return 𝐗, 𝐕
end

# As before but output k vectors of data as per the other method
# used for the case m>1, k>1
function getData(t, m, k, n, noise)
    # create m identical data matrices and rotate them by different
    # random orthogonal matrices V_1,...,V_m
    𝐕=[randU(n) for i=1:m] # random orthogonal matrices
    # variables common to all subjects with unique variance profile across k
    X=[(abs2.(randn(n))).*randn(n, t) for s=1:k]
    # each subject has this common part plus a random part
    𝐗=[[𝐕[i]*((1-noise)*X[s] + noise*randn(n, t)) for i=1:m] for s=1:k]
    return 𝐗, 𝐕
end

# Generate data for tests
n, t=3, 10
X=[ 0.0161263  -0.228679  -0.610958;
    0.884571    1.53687    0.569085;
   -1.26195    -1.35456   -0.303478;
   -0.275838    0.775055  -0.461387;
    1.35221    -0.890299  -0.224943;
   -0.204674    0.318228   0.266376;
   -0.203296   -1.83068   -0.869954;
   -0.726922   -0.529861  -0.609168;
   -0.0809636  -0.107088   0.0315276;
    0.177739    0.124636   0.495547]

Y=[ 2.16758   -2.29723   -1.92301;
    3.48593   -0.357984  -0.457085;
    3.02924   -1.5643    -1.38311;
    1.53344   -1.79264   -1.75126;
    0.990636  -0.618866  -0.643236;
    1.37825    0.704042   0.847845;
    0.324096   0.43795    0.142093;
   -3.68015    0.539658   0.436929;
    0.141023  -0.303311  -0.384726;
   -0.903474   1.20838    0.944415]

Cx=Symmetric((X'*X)/t)
Cy=Symmetric((Y'*Y)/t)
Cxy=(X'*Y)/t

𝐗 = [genDataMatrix(t, n) for i = 1:5]
𝐘 = [genDataMatrix(t, n) for i = 1:5]
𝐂x=Vector{Hermitian}([Hermitian((X'*X)/t) for X in 𝐗])
𝐂y=Vector{Hermitian}([Hermitian((Y'*Y)/t) for Y in 𝐘])
𝐂xy=Vector{Matrix}([(X'*Y)/size(X, 1) for (X, Y) in zip(𝐗, 𝐘)])

X_=[X, Matrix(X')]

############  Begin Test  ###################################

@testset "PCA" begin
    ## method (1)
    pC=pca(Cx; simple=true)
    @test Cx≈pC.F*pC.D*pC.iF

    ## method (2)
    for d=1:2
        p=pca(X_[d]; dims=d)
        compare(d, X_[d], p.F, p.D)
        p=pca(X_[d]; dims=d, eVar=3)
        compare(d, X_[d], p.F, p.D)
        p=pca(X_[d]; dims=d, eVar=0.5)
        compare(d, X_[d], p.F, p.D)
        p=pca(X_[d]; dims=d, covEst=LShrLW, eVar=0.5)
        # p=pca(X_[d]; covEst=NShrLW, eVar=0.5)
    end

    p=pca(X; simple=true)
    @test p==pC

    ## method (3)
    p=pca(𝐗; covEst=SCM, eVar=0.5, metric=logEuclidean)
    G=PosDefManifold.mean(logEuclidean, 𝐂x)
    D=p.F'*G*p.F
    @test isApproxD(D)
    @test Diagonal(D)≈Diagonal(reverse(eigvals(G))[1: size(D, 1)])
end

@testset "Whitening" begin
    ## method (1)
    wC=whitening(Cx; simple=true)
    @test wC.F'*Cx*wC.F≈I

    ## method (2)
    for d=1:2
        w=whitening(X_[d]; dims=d)
        compare(d, X_[d], w.F, Diagonal(ones(eltype(w.F), size(w.F, 2))))
        w=whitening(X_[d]; dims=d, eVar=3, eVarMeth=searchsortedlast)
        compare(d, X_[d], w.F, Diagonal(ones(eltype(w.F), size(w.F, 2))))
        w=whitening(X_[d]; dims=d, eVar=0.5, eVarMeth=searchsortedlast)
        compare(d, X_[d], w.F, Diagonal(ones(eltype(w.F), size(w.F, 2))))
        w=whitening(X_[d]; covEst=LShrLW, eVar=0.5)
        # w=whitening(X_[d]; covEst=NShrLW, eVar=0.5)
    end
    w=whitening(X; simple=true)
    @test w==wC

    ## method (3)
    w=whitening(𝐗; covEst=SCM, simple=true, metric=logEuclidean)
    G=PosDefManifold.mean(logEuclidean, 𝐂x)
    @test w.F'*G*w.F≈I
end

@testset "MCA" begin
    ## method (1)
    mC=mca(Cxy, simple=true)
    @test Cxy≈mC.F[1]*mC.D*mC.F[2]'

    ## method (2)
    m=mca(X, Y, simple=true)
    @test isApproxD(m.F[1]'*Cxy*m.F[2])
    @test m==mC

    ## method (3)
    m=mca(𝐗, 𝐘; eVar=0.99)
    m𝐂=mca(mean(𝐂xy); eVar=0.99)
    @test m==m𝐂
end

@testset "CCA" begin
    c=cca(X, Y)
    @test c.F[1]'*Cx*c.F[1]≈I
    @test c.F[2]'*Cy*c.F[2]≈I
    @test isApproxD(c.F[1]'*Cxy*c.F[2])

    # CCA of two vectors of data matrices
    c=cca(𝐗, 𝐘; metric=Euclidean, eVar=0.99)
    c𝐂=cca(Symmetric(mean(𝐂x)), Symmetric(mean(𝐂y)), mean(𝐂xy); eVar=0.99)
    @test c==c𝐂
end

@testset "CSP" begin
    t, n=20, 10
    X=genDataMatrix(t, n)
    Y=genDataMatrix(t, n)
    Cx=Symmetric((X'*X)/t)
    Cy=Symmetric((Y'*Y)/t)

    c=csp(X, Y)
    c=csp(Cx, Cy)
    @test isApproxD(c.F'*Cx*c.F)
    @test isApproxD(c.F'*Cy*c.F)
    @test c.F'*(Cx+Cy)*c.F≈I

    d=csp(Cx, Cy; eVarC=0.)
    @test isApproxD(d.F'*Cx*d.F)
    @test isApproxD(d.F'*Cy*d.F)
    @test d.F'*(Cx+Cy)*d.F≈I

    # the two solutions are identical out of a possible sign ambiguity on the evectors
    @test abs.(c.iF*d.F)≈I
end

@testset "CSTP" begin
    t, n=10, 20
    𝐗 = [genDataMatrix(t, n) for i = 1:10]
    Xfixed=randn(t, n)./1
    for i=1:length(𝐗) 𝐗[i]+=Xfixed end
    C1=Hermitian( mean((X'*X)/t for X∈𝐗) )
    C2=Hermitian( mean((X*X')/n for X∈𝐗) )
    c=cstp(Xfixed, C1, C2; simple=true)
    @test c.F[1]'*C2*c.F[1]≈I
    @test c.F[2]'*C1*c.F[2]≈I
    Z=c.F[1]'*Xfixed*c.F[2]
    n=minimum(size(Z))
    @test isApproxD(Z[1:n, 1:n])
end

@testset "AJD" begin
    t, n, k=50, 10, 50
    A=randn(n, n) # mixing matrix
    𝐗 = [genDataMatrix(t, n) for i = 1:k]
    Xfixed=randn(t, n)./1
    for i=1:length(𝐗) 𝐗[i]+=Xfixed end
    𝐂 = ℍVector([ℍ((𝐗[s]'*𝐗[s])/t) for s=1:k])
    aX=ajd(𝐗; simple=true)
    aC=ajd(𝐂; simple=true)
    @test aX≈aC

    # # # Test orthogonal AJD algorithms
    # create 20 random commuting matrices
    # they all have the same eigenvectors
    𝐂=randP(3, 20; eigvalsSNR=Inf, commuting=true)
    # estimate the approximate joint diagonalizer (ajd) using orthogonal solvers
    # the ajd must be equivalent to the eigenvector matrix of any of the matrices in 𝐂
    a=ajd(𝐂; algorithm=:OJoB)
    @test norm([spForm(a.F'*eigvecs(C)) for C ∈ 𝐂])/20<√err
    a=ajd(𝐂; algorithm=:JADE)
    @test norm([spForm(a.F'*eigvecs(C)) for C ∈ 𝐂])/20<√err
    a=ajd(𝐂; algorithm=:JADEmax)
    @test norm([spForm(a.F'*eigvecs(C)) for C ∈ 𝐂])/20<√err

    # # # Test non-orthogonal AJD algorithms
    # generate positive definite matrices with model A*D_κ*D, where A is an
    # invertible mixing matrix and D_κ, for all κ=1:k, are diagonal matrices.
    # The estimated AJD matrix must be the inverse of A
    # and all transformed matrices bust be diagonal
    n, k=3, 10
    Dest=PosDefManifold.randΛ(eigvalsSNR=Inf, n, k)
    # make the problem identifiable
    for i=1:k Dest[k][1, 1]*=i/(k/2) end
    for i=1:k Dest[k][3, 3]/=i/(k/2) end
    A=randn(n, n) # non-singular mixing matrix
    Cset3=Vector{Hermitian}([Hermitian(A*D*A') for D ∈ Dest])
    a=ajd(Cset3; algorithm=:NoJoB, eVarC=n)
    @test spForm(a.F'*A)<√err
    @test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<err
    a=ajd(Cset3; algorithm=:LogLike, eVarC=n)
    @test spForm(a.F'*A)<√err
    @test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<err
    a=ajd(Cset3; algorithm=:LogLikeR, eVarC=n)
    @test spForm(a.F'*A)<√err
    @test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<err
    a=ajd(Cset3; algorithm=:GAJD, eVarC=n)
    @test spForm(a.F'*A)<√err*10 # GAJD has problems sometimes
    @test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<err
    a=ajd(Cset3; algorithm=:QNLogLike, eVarC=n)
    @test spForm(a.F'*A)<√err
    @test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<err
end

@testset "gMCA" begin
    t, m, n, noise = 200, 2, 6, 0.1
    𝐗, 𝐕 = getData(t, m, n, noise)
    Cx=(𝐗[1]*𝐗[1]')/t
    Cy=(𝐗[2]*𝐗[2]')/t
    Cxy=(𝐗[1]*𝐗[2]')/t

    # check that for the case m=2 GMCA gives the same result as MCA
    gm=gmca(𝐗; simple=true)

    # do MCA
    m=mca(Cxy; simple=true)

    @test (m.F[1]'*Cxy*m.F[2]) ≈ (gm.F[1]'*Cxy*gm.F[2])
    # the following must be the identity matrix out of a possible sign ambiguity
    @test abs.(m.F[1]'*gm.F[1]) ≈ I
    @test abs.(m.F[2]'*gm.F[2]) ≈ I
end

@testset "gCCA" begin
    t, m, n, noise = 200, 2, 6, 0.1
    𝐗, 𝐕 = getData(t, m, n, noise)
    Cx=(𝐗[1]*𝐗[1]')/t
    Cy=(𝐗[2]*𝐗[2]')/t
    Cxy=(𝐗[1]*𝐗[2]')/t

    # check that for the case m=2 GCCA gives the same result as CCA
    gc=gcca(𝐗; simple=true)
    gm=gmca(𝐗; simple=true)
    # do CCA
    c=cca(𝐗[1], 𝐗[2]; simple=true)

    @test isApproxD(gc.F[1]'*Cxy*gc.F[2])

    @test (c.F[1]'*Cxy*c.F[2]) ≈ (gc.F[1]'*Cxy*gc.F[2])
    @test (gc.F[1]'*Cx*gc.F[1]) ≈ I
    @test (gc.F[2]'*Cy*gc.F[2]) ≈ I

    # the following must be the identity matrix out of a possible sign ambiguity
    @test abs.(c.iF[1]*gc.F[1]) ≈ I
    @test abs.(c.iF[2]*gc.F[2]) ≈ I

    # test that gCCA is equivalent to GMCA on pre-whitened data
    w1=whitening(𝐗[1])
    w2=whitening(𝐗[2])
    Y=[w1.F'*𝐗[1], w2.F'*𝐗[2]]
    C1=(Y[1]*Y[1]')/t
    C2=(Y[2]*Y[2]')/t
    C12=(Y[1]*Y[2]')/t
    gm2=gmca(Y; simple=true, algorithm=:OJoB)

    gm2=gcca(Y; simple=true, algorithm=:OJoB)
    @test gm2.F[1]'*C1*gm2.F[1] ≈ I
    @test gm2.F[2]'*C2*gm2.F[2] ≈ I
    @test isApproxD(gm2.F[1]'*C12*gm2.F[2])
    @test gm2.F[1]'*w1.F'*Cx*w1.F*gm2.F[1] ≈ I
    @test gm2.F[2]'*w2.F'Cy*w2.F*gm2.F[2] ≈ I
    @test isApproxD(gm2.F[1]'*w1.F'*Cxy*w2.F*gm2.F[2])
end

@testset "mAJD" begin
    # do joint blind source separation of non-stationary data
    t, m, n, k, noise = 200, 10, 3, 20, 0.01
    𝐗, 𝐕=getData(t, m, k, n, noise)
    aX=majd(𝐗; fullModel=true, algorithm=:OJoB)
    # the spForm index of the estimated demixing matrices times the true
    # mixing matrix must be low
    @test mean(spForm(aX.F[i]'*𝐕[i]) for i=1:m)<0.05
end


text:
## header to provide surface and context
using Cairo
using Colors

c = CairoRGBSurface(256,256)
cr = CairoContext(c)

save(cr)

z = Array{RGB24}(undef, 2, 2)
c1 = convert(RGB24,colorant"grey20")
c2 = convert(RGB24,colorant"grey80")
z[1,1] = c1
z[1,2] = c2
z[2,1] = c2
z[2,2] = c1

img = CairoImageSurface(z)
pattern = CairoPattern(img)
pattern_set_extend(pattern, Cairo.EXTEND_REPEAT)
pattern_set_filter(pattern, Cairo.FILTER_BILINEAR)

m = CairoMatrix(1/8.0,0,0,1/8.0,0,0)

set_matrix(pattern, m)
set_source(cr, pattern)
paint(cr)
restore(cr)

save(cr)

# 5 uses of set_source

# color
save(cr)
rectangle(cr,16,32,224,32)
set_source(cr, colorant"red4")
fill(cr)
restore(cr)

# color with alpha
save(cr)
rectangle(cr,16,72,224,32)
set_source(cr, alphacolor(colorant"blue",0.5))
fill(cr)
restore(cr)

# image from surface
save(cr)
rectangle(cr,16,112,224,32)
clip(cr)
new_path(cr)
s = read_from_png("data/mulberry.png")
set_source(cr,s)
paint_with_alpha(cr,0.6)
restore(cr)

# image from context (from surface)
save(cr)
rectangle(cr,16,152,224,32)
clip(cr)
new_path(cr)
set_source(cr,creategc(s))
paint(cr)
restore(cr)

# image for surface with offset
save(cr)
rectangle(cr,16,192,224,32)
clip(cr)
set_source(cr,s,0.0,40.0)
paint(cr)
restore(cr)

## mark picture with current date
restore(cr)
move_to(cr,0.0,12.0)
set_source_rgb(cr, 0,0,0)
show_text(cr,Libc.strftime(time()))
write_to_png(c,"sample_alpha_paint.png")


text:
"""
    function Element1D(matrix::Array{Int, 1}, i::AbstractIntVar, x::AbstractIntVar, trailer::Trailer)

Create a constraint stating that `matrix[i] == x`
"""

function Element1D(matrix::Array{Int, 1}, i::AbstractIntVar, x::AbstractIntVar, trailer::Trailer)
    return Element2D(matrix[:, :], i, IntVar(1, 1, "SEAPEARL_one", trailer), x, trailer)
end


text:
#using Morbit
#using Test

@testset "batch_objectives" begin
global counter_1 = 0;
global counter_2 = 0;

f = x -> sum( x.^2 )

# this function does NOT handle array of vector inputs and uses standard broadcasting
function objf1( x_vec )
    global counter_1
    counter_1 += 1;
    f(x_vec)
end

# this function handles batch evaluation all by itself,
# and broadcasting should be overwritten internally
function objf2( x_vec_or_array )
    global counter_2
    counter_2 += 1;
    if isa( x_vec_or_array , Vector{<:Real} )
        return f(x_vec_or_array)
    elseif isa(x_vec_or_array , Vector{<:Vector{<:Real}} )
        return [ f(x) for x ∈ x_vec_or_array ]
    end
end

n_vars = 3
sites = [ rand(n_vars) for i = 1 : 5 ];

mop = MixedMOP(n_vars)
add_objective!(mop, objf1, :expensive, 1, false)
add_objective!(mop, objf2, :expensive, 1, true)

Y = Morbit.eval_all_objectives.(mop, sites)

@test counter_1 == 5
@test counter_2 == 1
@test all( y[1] == y[2] for y in Y )
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: mediatailor
using AWS.Compat
using AWS.UUIDs

"""
    DeletePlaybackConfiguration()

Deletes the playback configuration for the specified name.

# Required Parameters
- `Name`: The identifier for the playback configuration.

"""
delete_playback_configuration(Name; aws_config::AWSConfig=global_aws_config()) = mediatailor("DELETE", "/playbackConfiguration/$(Name)"; aws_config=aws_config)
delete_playback_configuration(Name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediatailor("DELETE", "/playbackConfiguration/$(Name)", args; aws_config=aws_config)

"""
    GetPlaybackConfiguration()

Returns the playback configuration for the specified name.

# Required Parameters
- `Name`: The identifier for the playback configuration.

"""
get_playback_configuration(Name; aws_config::AWSConfig=global_aws_config()) = mediatailor("GET", "/playbackConfiguration/$(Name)"; aws_config=aws_config)
get_playback_configuration(Name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediatailor("GET", "/playbackConfiguration/$(Name)", args; aws_config=aws_config)

"""
    ListPlaybackConfigurations()

Returns a list of the playback configurations defined in AWS Elemental MediaTailor. You can specify a maximum number of configurations to return at a time. The default maximum is 50. Results are returned in pagefuls. If MediaTailor has more configurations than the specified maximum, it provides parameters in the response that you can use to retrieve the next pageful.

# Optional Parameters
- `MaxResults`: Maximum number of records to return.
- `NextToken`: Pagination token returned by the GET list request when results exceed the maximum allowed. Use the token to fetch the next page of results.
"""
list_playback_configurations(; aws_config::AWSConfig=global_aws_config()) = mediatailor("GET", "/playbackConfigurations"; aws_config=aws_config)
list_playback_configurations(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = mediatailor("GET", "/playbackConfigurations", args; aws_config=aws_config)

"""
    ListTagsForResource()

Returns a list of the tags assigned to the specified playback configuration resource.

# Required Parameters
- `ResourceArn`: The Amazon Resource Name (ARN) for the playback configuration. You can get this from the response to any playback configuration request.

"""
list_tags_for_resource(ResourceArn; aws_config::AWSConfig=global_aws_config()) = mediatailor("GET", "/tags/$(ResourceArn)"; aws_config=aws_config)
list_tags_for_resource(ResourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediatailor("GET", "/tags/$(ResourceArn)", args; aws_config=aws_config)

"""
    PutPlaybackConfiguration()

Adds a new playback configuration to AWS Elemental MediaTailor.

# Optional Parameters
- `AdDecisionServerUrl`: The URL for the ad decision server (ADS). This includes the specification of static parameters and placeholders for dynamic parameters. AWS Elemental MediaTailor substitutes player-specific and session-specific parameters as needed when calling the ADS. Alternately, for testing you can provide a static VAST URL. The maximum length is 25,000 characters.
- `AvailSuppression`: The configuration for Avail Suppression. Ad suppression can be used to turn off ad personalization in a long manifest, or if a viewer joins mid-break.
- `Bumper`: The configuration for bumpers. Bumpers are short audio or video clips that play at the start or before the end of an ad break.
- `CdnConfiguration`: The configuration for using a content delivery network (CDN), like Amazon CloudFront, for content and ad segment management.
- `DashConfiguration`: The configuration for DASH content.
- `LivePreRollConfiguration`: The configuration for pre-roll ad insertion.
- `ManifestProcessingRules`: The configuration for manifest processing rules. Manifest processing rules enable customization of the personalized manifests created by MediaTailor.
- `Name`: The identifier for the playback configuration.
- `PersonalizationThresholdSeconds`: The maximum duration of underfilled ad time (in seconds) allowed in an ad break.
- `SlateAdUrl`: The URL for a high-quality video asset to transcode and use to fill in time that's not used by ads. AWS Elemental MediaTailor shows the slate to fill in gaps in media content. Configuring the slate is optional for non-VPAID configurations. For VPAID, the slate is required because MediaTailor provides it in the slots that are designated for dynamic ad content. The slate must be a high-quality asset that contains both audio and video.
- `TranscodeProfileName`: The name that is used to associate this playback configuration with a custom transcode profile. This overrides the dynamic transcoding defaults of MediaTailor. Use this only if you have already set up custom profiles with the help of AWS Support.
- `VideoContentSourceUrl`: The URL prefix for the master playlist for the stream, minus the asset ID. The maximum length is 512 characters.
- `tags`: The tags to assign to the playback configuration.
"""
put_playback_configuration(; aws_config::AWSConfig=global_aws_config()) = mediatailor("PUT", "/playbackConfiguration"; aws_config=aws_config)
put_playback_configuration(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = mediatailor("PUT", "/playbackConfiguration", args; aws_config=aws_config)

"""
    TagResource()

Adds tags to the specified playback configuration resource. You can specify one or more tags to add.

# Required Parameters
- `ResourceArn`: The Amazon Resource Name (ARN) for the playback configuration. You can get this from the response to any playback configuration request.
- `tags`: A comma-separated list of tag key:value pairs. For example:
 {
 \"Key1\": \"Value1\",
 \"Key2\": \"Value2\"
 }

"""
tag_resource(ResourceArn, tags; aws_config::AWSConfig=global_aws_config()) = mediatailor("POST", "/tags/$(ResourceArn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(ResourceArn, tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediatailor("POST", "/tags/$(ResourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), args)); aws_config=aws_config)

"""
    UntagResource()

Removes tags from the specified playback configuration resource. You can specify one or more tags to remove.

# Required Parameters
- `ResourceArn`: The Amazon Resource Name (ARN) for the playback configuration. You can get this from the response to any playback configuration request.
- `tagKeys`: A comma-separated list of the tag keys to remove from the playback configuration.

"""
untag_resource(ResourceArn, tagKeys; aws_config::AWSConfig=global_aws_config()) = mediatailor("DELETE", "/tags/$(ResourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(ResourceArn, tagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediatailor("DELETE", "/tags/$(ResourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), args)); aws_config=aws_config)


text:
####
#### Graph
####

"""
    Graph

Represents a simple, undirected, unweighted graph.
"""
mutable struct Graph
    num_edges::Int
    adj_list::Vector{Vector{Int}}
end

"""
    Graph(n=0)

Create an empty (no edges) `Graph` with `n` vertices.
"""
Graph(n::Integer) = add_vertex!(Graph(), n)
Graph() = Graph(0, Vector{Int}[])

## Make Graph behave as a scalar in broadcasting expressions. That is, when
## automatically broadcasting a scalar function over arrays, make Graph behave
## like a scalar would.
Base.Broadcast.broadcastable(g::Graph) = Ref(g)

# Print an informative, but non-machine-readable representation.
function Base.show(io::IO, g::Graph)
    print(io, "Graph(nv=$(nv(g)), ne=$(ne(g)))")
    return nothing
end

"""
    adj_list(g::Graph)

Return `g` as an adjacency list. This does *not* return a copy.
"""
adj_list(g::Graph) = g.adj_list

"""
    adj_list(g::Graph, v::Integer)

Return a list of vertices adjacent to `v`.
"""
adj_list(g::Graph, v::Integer) = adj_list(g)[v]

Base.length(g::Graph) = length(adj_list(g))

"""
    ne(g::Graph)

Return the number of edges in `g`.
"""
ne(g::Graph) = g.num_edges

"""
    nv(g::Graph)

Return the number of vertices in `g`.
"""
nv(g::Graph) = length(g)

# True if there are no edges.
Base.isempty(g::Graph) = iszero(nv(g))

## TODO: Is returning essentially 1:0 a good idea for an empty graph? It allows
## iterating correctly over an empty graph, but could lead to bugs if one asks
## for `first(vertices(g))`
"""
    vertices(g::Graph)

Return an iterator over the vertices in `g`.
"""
function vertices(g::Graph)
    return Base.OneTo(nv(g))
end

"""
    add_vertex!(g::Graph, n::Integer=1)

Add `n` vertices to `g`. No edges are added.
"""
function add_vertex!(g::Graph, n::Integer=1)
    for i in 1:n
        push!(adj_list(g), Int[])
    end
    return g
end

"""
    add_edge!(adj_list_entry::Vector, v)

Add vertex `v` to the list of edges `adj_list_entry`, which represents all edges
terminating in a single vertex.  Information about this single vertex is not
neccessary nor obtainable from within this function.
"""
function add_edge!(adj_list::Vector, v)
    i = searchsortedfirst(adj_list, v)
    len = length(adj_list)
    if len == 0
        insert!(adj_list, 1, v)
    else
#        i > len || adj_list[i] == v && error("Edge already present")
        i > len || adj_list[i] == v && return nothing
        insert!(adj_list, i, v)
    end
    return adj_list
end

"""
    add_edge!(g::Graph, v1, v2)

Add the edge `(v1, v2)` to `g`.
"""
function add_edge!(g::Graph, v1, v2)
    result = add_edge!(g.adj_list[v1], v2)
    result == nothing && return nothing
    add_edge!(g.adj_list[v2], v1)
    g.num_edges += 1
    return g
end

"""
    del_edge!(adj_list_entry::Vector, v::Integer)

Delete the edge terminating in`v` from `adj_list_entry`.
The other vertex is neither needed nor referenced in this
function.
"""
function del_edge!(adj_list::Vector, v::Integer)
    i = searchsortedfirst(adj_list, v)
    i > length(adj_list) && error("Edge not present")
    deleteat!(adj_list, i)
    return adj_list
end

"""
    del_edge!(g::Graph, v1, v2)

Delete the edge `(v1, v2)` from `g`.
"""
function del_edge!(g::Graph, v1, v2)
    del_edge!(g.adj_list[v1], v2)
    del_edge!(g.adj_list[v2], v1)
    g.num_edges -= 1
    return g
end

####
#### Shortest path length
####

## Dijkstra's Algorithm.
## See https://en.wikipedia.org/wiki/Dijkstra%27s_algorithm
"""
    shortest_path_length(g::Graph, v0::Integer)

Return a list of distances from all vertices in `g` to vertex `v0` "distance"
means the length of a shortest path.
"""
function shortest_path_length(g::Graph, v0::Integer)
    q = PriorityQueue()
    dist = zeros(Int, nv(g))
    dist[v0] = 0 # Re-zero for insurance against future change.

    for v in vertices(g)
        if v != v0 # The distance from `v0` to `v0` is always zero
            dist[v] = typemax(eltype(dist)) # all distances are infinity initially
        end
        insert!(q, dist[v], v) # Put vertex `v` in the queue with priority `dist[v]`.
    end

    while ! isempty(q)
        u = pop!(q) # Get the vertex (provisionally) closest to `v0`.
        for v in adj_list(g, u) # Iterate over all edges `(v, u)` terminating in `u`.
            ! in(v, q) && continue # Skip vertices already examined.
            alt = dist[u] + 1 # Add the distance between two vertices sharing an edge: 1
            if alt < dist[v] # Did we find a shorter path ?
                dist[v] = alt
                decrease_priority!(q, alt, v)
            end
        end
    end
    for i in eachindex(dist)
        if dist[i] < 0 || dist[i] == typemax(Int)
            dist[i] = -1 # replace representations of infinity by -1 for clarity.
        end
    end
    return dist
end

## FIXME: It would be more efficient to terminate Dijkstra's algorthm after the
## target vertex has been examined. Instead, we are computing all distances.
"""
    shortest_path_length(g::Graph, v0::Integer, v1::Integer)

Return a length of the shortest path between `v0` and `v1`.
"""
function shortest_path_length(g::Graph, v0::Integer, v1::Integer)
    dist = shortest_path_length(g, v0)
    return dist[v1]
end


text:
mkdir(resultsdir)
mkdir(resultsdir*"/tags")
mkdir(resultsdir*"/images")

fid = open(resultsdir*"/readme.txt", "w")
println(fid, currdirtime)
println(fid, datafolder)
println(fid, camidxs)
println(fid, ARGS)
close(fid)

fid = open(resultsparentdir*"/racecar.log", "a")
println(fid, "$(currdirtime), $datafolder, $(camidxs), $(ARGS)")
close(fid)


text:
module DeprecatedTest

using Base.Test
using ForwardDiff, DiffResults

using ForwardDiff: AbstractConfig, GradientConfig,
                   JacobianConfig, HessianConfig,
                   MultithreadConfig

include(joinpath(dirname(@__FILE__), "utils.jl"))

function similar_duals(a::AbstractArray, b::AbstractArray)
    return typeof(a) == typeof(b) && size(a) == size(b)
end

similar_duals(a::Tuple, b::Tuple) = all(similar_duals.(a, b))

function similar_config(a::AbstractConfig, b::AbstractConfig)
    return a.seeds == b.seeds && similar_duals(a.duals, b.duals)
end

function similar_config(a::HessianConfig, b::HessianConfig)
    return (similar_config(a.gradient_config, b.gradient_config) &&
            similar_config(a.jacobian_config, b.jacobian_config))
end

x = rand(3)
y = rand(3)
out = DiffResults.HessianResult(x)
N = 1
chunk = ForwardDiff.Chunk{N}()

info("The following tests print lots of deprecation warnings on purpose.")

@test similar_config(GradientConfig{N}(x), GradientConfig(nothing, x, chunk))
@test similar_config(JacobianConfig{N}(x), JacobianConfig(nothing, x, chunk))
@test similar_config(JacobianConfig{N}(y, x), JacobianConfig(nothing, y, x, chunk))
@test similar_config(HessianConfig{N}(x), HessianConfig(nothing, x, chunk))
@test similar_config(HessianConfig{N}(out, x), HessianConfig(nothing, out, x, chunk))
@test similar_config(MultithreadConfig(GradientConfig(nothing, x, chunk)), GradientConfig(nothing, x, chunk))

info("Deprecation testing is now complete, so any further deprecation warnings are real.")

end # module


text:
include(joinpath("structures", "cluster.jl"))
include(joinpath("structures", "dataitem.jl"))
include(joinpath("structures", "eva.jl"))
include(joinpath("structures", "fittedeva.jl"))


text:
module PMTK4BasicModels

    include("gamma/gamma.jl")
    export normpdf
    include("norm/norm.jl")
    export contingencyTableUnpairedDiffPostMc
    include("sub/contingencyTableUnpairedDiffPostMc.jl")

end


text:
function _default_agent_update_logic!(agent::VLAgentModel, order::VLOrderModel, price_array::Array{Float64,2}, iteration_index::Int64)

    try

        # get stuff from the order -
        order_quantity = order.quantity                   # what is the order quantity?
        order_action = order.action                       # what is the order action?
        asset_index = order.asset_index                   # what asset was associated with this order

        # get stuff the agent -
        new_agent_position_array = copy(agent.agent_position_array)

        # update the position array -
        current_position = agent.agent_position_array[iteration_index - 1, asset_index]
        𝒪 = 0
        if (order_action == 1)
            𝒪 = current_position + order_quantity
        elseif (order_action == -1)
            𝒪 = current_position + order_quantity
        else
            𝒪 = current_position
        end

        if (𝒪 >= 0)
            new_agent_position_array[iteration_index, asset_index] = 𝒪
        end

        # update -
        agent.agent_position_array = new_agent_position_array

    catch error
        throw(error)
    end
end

function _default_agent_trade_logic()

    try

    catch error
        throw(error)
    end
end


text:
using CSFML.LibCSFML
using ModernGL
include("ModernGL_utils.jl")

# https://www.tomdalling.com/blog/modern-opengl/01-getting-started-in-xcode-and-visual-cpp/

### Init

mode = sfVideoMode(700, 720, 32)
# settings = SF::ContextSettings.new(
#   depth: 24, stencil: 8, antialiasing: 4,
#   major: 3, minor: 0
# )
window = sfRenderWindow_create(mode, "SFML window", sfResize | sfClose, C_NULL)
@assert window != C_NULL
sfWindow_setVerticalSyncEnabled(window,true)
sfWindow_setActive(window, sfTrue)
println(createcontextinfo()) # used by ModernGL_utils;
glEnable(GL_DEPTH_TEST)
glDepthFunc(GL_LESS)
glClearColor(0.4 ,0.2, 0.3, 1.0)

### shaders

vsh = """
    $(get_glsl_version_string())
    //layout(location = 0) in vec3 v_position;
    //layout(location = 1) in vec3 v_color;

    in vec3 v_position;
    //in vec3 v_color;
    //out vec3 initial_Color;
    void main() {
        //initial_Color = v_color;
        gl_Position = vec4(v_position, 1.0);
    }
    """

fsh = """
    $(get_glsl_version_string())
    // in vec3 initial_Color;
    out vec4 outColor;
    void main() {
        outColor = vec4(5.0, 2.0, 3.0, 1.0);
        //outColor = vec4(frag_colour, 1.0 );
    }
    """

vertexShader = createShader(vsh, GL_VERTEX_SHADER)
fragmentShader = createShader(fsh, GL_FRAGMENT_SHADER)

function bindAttrib(shader_program)
    #glBindAttribLocation(shader_program, 0, "v_position");
    #glBindAttribLocation(shader_program, 1, "v_color");
end

program = createShaderProgram( bindAttrib, vertexShader, fragmentShader)
#glUseProgram(program)
#println("program : $program")
print("shader ---> ")
glCheckError()

### vertex, colors and buffers

vao = glGenVertexArray()
glBindVertexArray(vao)

points_vbo = glGenBuffer()
glBindBuffer(GL_ARRAY_BUFFER, points_vbo)

vertexData = GLfloat[
        #  X     Y     Z       U     V
         0.0, 0.8, 0.0,
        -0.8,-0.8, 0.0,
         0.8,-0.8, 0.0
]

glBufferData(GL_ARRAY_BUFFER, sizeof(vertexData), vertexData, GL_STATIC_DRAW)
println("points_vbo : $points_vbo")
glCheckError()

glEnableVertexAttribArray(0);
glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 0 , C_NULL)

glBindBuffer(GL_ARRAY_BUFFER, 0)
glBindVertexArray(0)

###

event_ref = Ref{sfEvent}()
clock = sfClock_create()
# Let's express time in seconds.
function get_time(clock::Ptr{sfClock})
    sfTime_asSeconds(sfClock_getElapsedTime(clock))
end
function restart(clock::Ptr{sfClock})
     sfTime_asSeconds(sfClock_restart(clock))
end

target_frequency = 24.0
running = true
t0 = 9999999
while (running)

    while Bool(sfRenderWindow_isOpen(window)) && running
        frame_timestamp = get_time(clock)
        frame_time = frame_timestamp - t0
        global t0 = frame_timestamp
        actual_frequency = 1 / frame_time

        # process events
        while Bool(sfRenderWindow_pollEvent(window, event_ref))
            # close window : exit
            event_ref[].type == sfEvtClosed && sfRenderWindow_close(window)
            event_ref[].type == sfEvtResized && println("Trigger sfEvtResized.")
            event_ref[].type == sfEvtLostFocus && println("Trigger sfEvtLostFocus.")
            event_ref[].type == sfEvtGainedFocus && println("Trigger sfEvtGainedFocus.")
            event_ref[].type == sfEvtTextEntered && println("Trigger sfEvtTextEntered: $(event_ref[].text.unicode)")
            event_ref[].type == sfEvtKeyPressed && println("Trigger sfEvtKeyPressed: $(event_ref[].key.code)")
            event_ref[].type == sfEvtKeyReleased && println("Trigger sfEvtKeyReleased: $(event_ref[].key.code)")
            event_ref[].type == sfEvtMouseWheelMoved && println("Trigger sfEvtMouseWheelMoved: $(event_ref[].mouseWheel.x), $(event_ref[].mouseWheel.y)")
            event_ref[].type == sfEvtMouseWheelScrolled && println("Trigger sfEvtMouseWheelScrolled: $(event_ref[].mouseWheelScroll.wheel)")
            event_ref[].type == sfEvtMouseButtonPressed && println("Trigger sfEvtMouseButtonPressed: $(event_ref[].mouseButton.button)")
            event_ref[].type == sfEvtMouseButtonReleased && println("Trigger sfEvtMouseButtonReleased: $(event_ref[].mouseButton.x), $(event_ref[].mouseButton.y)")
            event_ref[].type == sfEvtMouseMoved && println("Trigger sfEvtMouseMoved: $(event_ref[].mouseMove.x), $(event_ref[].mouseMove.y)")

            if event_ref[].type == sfEvtClosed
                global running = false
                println("Render window closed.")
            elseif event_ref[].type == sfEvtResized
                glViewport(9, 9, event_ref[].size.width, event_ref[].size.height)
                println("Viewport updated.")
            elseif event_ref[].type == sfEvtKeyPressed
                if event_ref[].key.code == sfKeyQ
                    global running = false
                    println("Quiting application.")
                end
            end
        end
        time_event_processing = round(get_time(clock) - t0, digits=9)

        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

        glUseProgram(program)
        glBindVertexArray(vao)
        glDrawArrays(GL_TRIANGLES, 0 , 3)
        glBindVertexArray(0)
        glUseProgram(0)

        time_render = round(get_time(clock) - time_event_processing - t0, digits=9)
        sfRenderWindow_display(window)
        time_to_sleep = 1 / target_frequency - time_render - time_event_processing

        if time_to_sleep > 0.005
            sleep(time_to_sleep-0.005)
        end
    end
end

sfRenderWindow_destroy(window)

max_stars_repo_path: GL_playground/tomdalling_openGL 0x01.jl

max_stars_repo_name: leroy-walton/playgrounds

max_issues_repo_path: GL_playground/tomdalling_openGL 0x01.jl

max_issues_repo_name: leroy-walton/playgrounds

max_forks_repo_path: GL_playground/tomdalling_openGL 0x01.jl

max_forks_repo_name: leroy-walton/playgrounds


text:
module FHColoredLightbeam

using ..Ahorn, Maple
using Random

@mapdef Entity "FrostHelper/ColoredLightbeam" ColoredLightbeam(x::Integer, y::Integer, width::Integer=32, height::Integer=24, flag::String="", rotation::Integer=0, color::String="ccffff")

const colors = sort(collect(keys(Ahorn.XNAColors.colors)))

function lightbeamFinalizer(entity::ColoredLightbeam)
    x, y = Ahorn.position(entity)
    nx, ny = Int.(entity.nodes[1])

    theta = atan(y - ny, x - nx) + pi / 2
    magnitude = sqrt((y - ny)^2 + (x - nx)^2)

    entity.rotation = round(Int, rad2deg(theta))
    entity.height = round(Int, max(24, magnitude))
end

Ahorn.editingOptions(entity::ColoredLightbeam) = Dict{String, Any}(
    "color" => colors
)

const placements = Ahorn.PlacementDict(
    "Colored Lightbeam (Frost Helper)" => Ahorn.EntityPlacement(
        ColoredLightbeam,
        "line",
        Dict{String, Any}(),
        lightbeamFinalizer
    )
)

Ahorn.minimumSize(entity::ColoredLightbeam) = 32, 24
Ahorn.resizable(entity::ColoredLightbeam) = true, true

function rotate(point::Tuple{Number, Number}, theta::Number)
    res = [
        cos(theta)  -sin(theta);
        sin(theta) cos(theta)
    ] * [point[1]; point[2]]

    return (res[1], res[2])
end

function getSelectionRect(entity::ColoredLightbeam)
    x, y = Ahorn.position(entity)

    theta = deg2rad(get(entity.data, "rotation", 0))
    width = round(Int, get(entity.data, "width", 32))
    height = round(Int, get(entity.data, "height", 24))

    points = Tuple{Number, Number}[
        (-width / 2, 0),
        (width / 2, 0),
        (-width / 2, height - 4),
        (width / 2, height - 4)
    ]

    rotated = rotate.(points, theta)

    tlx, tly = minimum(p -> p[1], rotated), minimum(p -> p[2], rotated)
    brx, bry = maximum(p -> p[1], rotated), maximum(p -> p[2], rotated)

    return Ahorn.Rectangle(x + tlx, y + tly, brx - tlx, bry - tly)
end

function Ahorn.selection(entity::ColoredLightbeam)
    return getSelectionRect(entity)
end

texture = "util/lightbeam"

function renderLightbeam(ctx::Ahorn.Cairo.CairoContext, entity::ColoredLightbeam)
    x, y = Ahorn.position(entity)

    lightbeam = Ahorn.getSprite(texture, "Gameplay")

    theta = deg2rad(get(entity.data, "rotation", 0))
    width = round(Int, get(entity.data, "width", 32))
    height = round(Int, get(entity.data, "height", 24))

    rng = Ahorn.getSimpleEntityRng(entity)

    Ahorn.Cairo.save(ctx)

    Ahorn.translate(ctx, x, y)
    Ahorn.rotate(ctx, theta)
    Ahorn.translate(ctx, width / 2, 0)
    Ahorn.rotate(ctx, pi / 2)
    Ahorn.scale(ctx, (height - 4) / lightbeam.width, 1)
    Ahorn.set_antialias(ctx, 1)

    Ahorn.Cairo.save(ctx)

    for i in 0:width - 1
        Ahorn.drawImage(ctx, texture, 0, i, alpha=0.3)
    end

    for i in 0:4:width - 1
        num = i * 0.6
        lineWidth = 4.0 + sin(num * 0.5 + 1.2) * 4.0
        length = height + sin(num * 0.25) * 8
        alpha = 0.6 + sin(num + 0.8) * 0.3
        offset = sin((num + i * 32) * 0.1 + sin(num * 0.05 + i * 0.1) * 0.25) * (width / 2.0 - lineWidth / 2.0)

        # Not part of official calculations, but makes the render a bit less boring
        offsetMultiplier = (rand(rng) - 0.5) * 2

        for j in 0:3
            Ahorn.drawImage(ctx, texture, 0, offset * offsetMultiplier + width / 2 + j, alpha=alpha / 1.5)
        end
    end

    Ahorn.Cairo.restore(ctx)
end

Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::ColoredLightbeam, room::Maple.Room) = renderLightbeam(ctx, entity)

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
module SetCover

using SparseArrays
using ..Diderot
using ..Diderot: Arc, Node, Layer

struct Instance
    values::Vector{Float64}         # should be negative (is maximized)
    sets::SparseMatrixCSC{Bool,Int} # rows == sets to cover

    function Instance(values, sets)
        @assert length(values) == size(sets, 2)
        @assert all(sum.(eachrow(sets)) .> 0) # feasible
        return new(values, sets)
    end
end

Base.length(instance::Instance) = length(instance.values)
Diderot.domain_type(instance::Instance) = Bool
Diderot.value_type(instance::Instance) = Float64
Diderot.initial_state(instance::Instance) = BitSet(1:size(instance.sets, 1))

function Diderot.transitions(instance::Instance, state, variable)
    results = Dict{Arc{BitSet, Bool, Float64}, BitSet}()

    # Sets that are covered by this element.
    covered = BitSet(findall(instance.sets[:, variable]))

    # true: select element, remove all sets that are covered.
    results[Arc(state, true, instance.values[variable])] =
        setdiff(state, covered)

    # false: don't select element, state is unchanged.
    #
    # Need to check if this could be infeasible. For now, assume that variable
    # order is static, then we can check whether this variable was the "last
    # chance" for any set to cover.
    for set in intersect(covered, state)
        if variable == last(findall(instance.sets[set, :]))
            # Must add this element, so we don't add "false" transition.
            return results
        end
    end

    # Don't need this element, so we add "false" transition.
    results[Arc(state, false, 0.0)] = state

    return results
end

# For relaxation:
# - Keep states with smallest sets-to-cover.
# - Merge states by intersection.
sort_by(tuple) = length(tuple.first)
merge = intersect
relax(max_width) = Diderot.RelaxAllInOne(max_width, merge, sort_by)

end


text:
# Part 1
checksum = 0;
for line in eachline(STDIN)
  input = parse.(Int, split(line, ' '))
  checksum = checksum + (maximum(input) - minimum(input))
end

println("Part 1: $checksum")


text:
# Translating
# https://proplot.readthedocs.io/en/latest/basics.html#Changing-rc-settings
# to Julia

# Import proplot using PyCall
using PyCall
proplot = pyimport("proplot")

# Update global settings in several different ways
proplot.rc.cycle = "colorblind"
proplot.rc.color = "gray6"
proplot.rc.update(Dict("fontname" => "Noto Sans"))
proplot.rc["figure.facecolor"] = "gray3"
proplot.rc.axesfacecolor = "gray4"

# Apply settings to figure with context()
f, axs = proplot.subplots(ncols=2, aspect=1, width=6, span=false, sharey=2)
[ax.format(suptitlesize = 11, toplabelcolor="gray6", linewidth=1.5) for ax in axs]

# Plot lines
N, M = 100, 6
values = 1:M+1
for ax in axs
    data = cumsum(rand(N, M) .- 0.5, dims=1)
    ax.plot(data, linewidth=3, cycle="Grays")
end

# Apply settings to axes with format()
[ax.format(
    grid=false, xlabel="x label", ylabel="y label",
    collabels=["Column label 1", "Column label 2"],
    suptitle="Rc settings demo",
    suptitlecolor="gray7",
    abc=true, abcloc="l", abcstyle="A)",
    title="Title", titleloc="r", titlecolor="gray7"
) for ax in axs]
ay = axs[end].twinx()
ay.format(ycolor="red", linewidth=1.5, ylabel="secondary axis")
ay.plot(cumsum(rand(100) .- 0.2), color="r", lw=3)

# Reset persistent modifications from head of cell
proplot.rc.reset()
f.savefig("figs/Changing-rc-settings.svg", transparent=false)

max_stars_repo_name: natgeo-wong/JuliaProPlotExamples

max_issues_repo_name: natgeo-wong/JuliaProPlotExamples

max_forks_repo_name: natgeo-wong/JuliaProPlotExamples


text:
IDmm = DataFrame(NPS = [0.125,
	0.25,
	0.375,
	0.5,
	0.75,
	1,
	1.25,
	1.5,
	2,
	2.5,
	3,
	3.5,
	4,
	5,
	6,
	8,
	10,
	12,
	14,
	16,
	18,
	20,
	22,
	24,
	26,
	28,
	30,
	32,
	36,
	42,
	48,
	0.125,
	0.25,
	0.375,
	0.5,
	0.75,
	1,
	1.25,
	1.5,
	2,
	2.5,
	3,
	3.5,
	4,
	5,
	6,
	8,
	10,
	12,
	14,
	16,
	18,
	20,
	22,
	24,
	26,
	28,
	30,
	32,
	36,
	42,
	48,
	0.125,
	0.25,
	0.375,
	0.5,
	0.75,
	1,
	1.25,
	1.5,
	2,
	2.5,
	3,
	3.5,
	4,
	5,
	6,
	8,
	10,
	12,
	14,
	16,
	18,
	20,
	22,
	24,
	26,
	28,
	30,
	32,
	36,
	42,
	48,
	0.5,
	0.75,
	1,
	1.25,
	1.5,
	2,
	2.5,
	3,
	4,
	5,
	6,
	8,
	10,
	12,
	14,
	16,
	18,
	20,
	22,
	24,
	26,
	28,
	30,
	32,
	36,
	8,
	10,
	12,
	14,
	16,
	18,
	20,
	22,
	24,
	26,
	28,
	30,
	32,
	36,
	8,
	10,
	12,
	14,
	16,
	18,
	20,
	22,
	24,
	26,
	28,
	30,
	32,
	36,
	0.125,
	0.25,
	0.375,
	0.5,
	0.75,
	1,
	1.25,
	1.5,
	2,
	2.5,
	3,
	3.5,
	4,
	5,
	6,
	8,
	10,
	12,
	14,
	16,
	18,
	20,
	24,
	32,
	36,
	8,
	10,
	12,
	14,
	16,
	18,
	20,
	22,
	24,
	0.125,
	0.25,
	0.375,
	0.5,
	0.75,
	1,
	1.25,
	1.5,
	2,
	2.5,
	3,
	3.5,
	4,
	5,
	6,
	8,
	10,
	12,
	14,
	16,
	18,
	20,
	22,
	24,
	8,
	10,
	12,
	14,
	16,
	18,
	20,
	24,
	4,
	5,
	6,
	8,
	10,
	12,
	14,
	16,
	18,
	20,
	24,
	8,
	10,
	12,
	14,
	16,
	18,
	20,
	24,
	0.5,
	0.75,
	1,
	1.25,
	1.5,
	2,
	2.5,
	3,
	4,
	5,
	6,
	8,
	10,
	12,
	14,
	16,
	18,
	20,
	24],
Schedule = ["OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"OD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"STD",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XS",
"XXS",
"XXS",
"XXS",
"XXS",
"XXS",
"XXS",
"XXS",
"XXS",
"XXS",
"XXS",
"XXS",
"XXS",
"XXS",
"XXS",
"SCH10",
"SCH10",
"SCH10",
"SCH10",
"SCH10",
"SCH10",
"SCH10",
"SCH10",
"SCH10",
"SCH10",
"SCH10",
"SCH20",
"SCH20",
"SCH20",
"SCH20",
"SCH20",
"SCH20",
"SCH20",
"SCH20",
"SCH20",
"SCH20",
"SCH20",
"SCH20",
"SCH20",
"SCH20",
"SCH30",
"SCH30",
"SCH30",
"SCH30",
"SCH30",
"SCH30",
"SCH30",
"SCH30",
"SCH30",
"SCH30",
"SCH30",
"SCH30",
"SCH30",
"SCH30",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH40",
"SCH60",
"SCH60",
"SCH60",
"SCH60",
"SCH60",
"SCH60",
"SCH60",
"SCH60",
"SCH60",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH80",
"SCH100",
"SCH100",
"SCH100",
"SCH100",
"SCH100",
"SCH100",
"SCH100",
"SCH100",
"SCH120",
"SCH120",
"SCH120",
"SCH120",
"SCH120",
"SCH120",
"SCH120",
"SCH120",
"SCH120",
"SCH120",
"SCH120",
"SCH140",
"SCH140",
"SCH140",
"SCH140",
"SCH140",
"SCH140",
"SCH140",
"SCH140",
"SCH160",
"SCH160",
"SCH160",
"SCH160",
"SCH160",
"SCH160",
"SCH160",
"SCH160",
"SCH160",
"SCH160",
"SCH160",
"SCH160",
"SCH160",
"SCH160",
"SCH160",
"SCH160",
"SCH160",
"SCH160",
"SCH160"],
Idmm = [10.287,
	13.716,
	17.145,
	21.336,
	26.67,
	33.401,
	42.164,
	48.26,
	60.325,
	73.025,
	88.9,
	101.6,
	114.3,
	141.3002,
	168.275,
	219.075,
	273.05,
	323.85,
	355.6,
	406.4,
	457.2,
	508,
	558.8,
	609.6,
	660.4,
	711.2,
	762,
	812.8,
	914.4,
	1066.8,
	1219.2,
	6.84,
	9.22,
	11.56,
	15.7,
	20.9,
	26.6,
	35,
	40.9,
	52.5,
	62.6,
	77.9,
	90.2,
	102.3,
	128.1,
	154.1,
	202.7,
	254.5,
	304.9,
	336.54,
	387.34,
	438.14,
	488.94,
	539.74,
	590.54,
	641.34,
	692.14,
	742.94,
	793.74,
	895.34,
	1047.74,
	1200.14,
	5.48,
	7.66,
	10.7,
	13.9,
	18.9,
	24.4,
	32.4,
	38.1,
	49.3,
	59,
	73.7,
	85.4,
	97.1,
	122.3,
	146.3,
	193.7,
	247.7,
	298.5,
	330.2,
	381,
	431.8,
	482.6,
	533.4,
	584.2,
	635,
	685.8,
	736.6,
	787.4,
	889,
	1041.4,
	1193.8,
	6.3,
	11.1,
	15.2,
	22.8,
	27.9,
	38.1,
	45,
	58.5,
	80.1,
	103.1,
	124.5,
	174.7,
	222.3,
	273.1,
	342.9,
	393.7,
	444.5,
	495.3,
	546.1,
	596.9,
	644.56,
	695.36,
	746.16,
	796.96,
	898.56,
	206.3,
	260.3,
	311.1,
	339.76,
	390.56,
	441.36,
	488.94,
	539.74,
	590.54,
	635,
	685.8,
	736.6,
	787.4,
	889,
	205.1,
	257.5,
	307.1,
	336.54,
	387.34,
	434.94,
	482.6,
	533.4,
	581.06,
	628.64,
	679.44,
	730.24,
	781.04,
	882.64,
	6.84,
	9.22,
	11.56,
	15.7,
	20.9,
	26.6,
	35,
	40.9,
	52.5,
	62.6,
	77.9,
	90.2,
	102.3,
	128.1,
	154.1,
	202.7,
	254.5,
	303.3,
	333.34,
	381,
	428.66,
	477.82,
	574.64,
	777.84,
	876.3,
	198.5,
	247.7,
	295.3,
	325.42,
	373.08,
	419.1,
	466.76,
	514.34,
	560.38,
	5.48,
	7.66,
	10.7,
	13.9,
	18.9,
	24.4,
	32.4,
	38.1,
	49.3,
	59,
	73.7,
	85.4,
	97.1,
	122.3,
	146.3,
	193.7,
	242.9,
	288.9,
	317.5,
	363.52,
	409.54,
	455.62,
	501.6,
	547.68,
	188.9,
	236.5,
	281.1,
	307.94,
	354.02,
	398.48,
	442.92,
	531.82,
	92.1,
	115.9,
	139.7,
	182.5,
	230.3,
	273.1,
	300.02,
	344.48,
	387.36,
	431.8,
	517.56,
	177.9,
	222.3,
	266.7,
	292.1,
	333.36,
	377.82,
	419.1,
	504.82,
	11.7,
	15.5,
	20.6,
	29.4,
	34.1,
	42.9,
	54,
	66.7,
	87.3,
	109.5,
	131.7,
	173.1,
	215.9,
	257.3,
	284.18,
	325.42,
	366.72,
	407.98,
	490.52] )


text:
using .RationalFunctionFields: RationalFunctionField, compute_groebner!,
                    contains_using_groebner, contains_randomized

logger = Logging.SimpleLogger(stderr, Logging.Debug)
Logging.global_logger(logger)

@testset "Basics tests" begin

end

@testset "Contains tests" begin
    FF = Nemo.QQ
    R, (x1, x2) = Nemo.PolynomialRing(FF, ["x1", "x2"])
    # set = [ (a^2 + b^2) // 1, a*b // 1 ]
    # set = [ a*b // (a + b), (a + b) // 1, (a^3 - b) // 1 ]
    set = [
          ( x1^2 + x2^2 ) // R(1),
          ( x1^3 + x2^3 ) // R(1),
          ( x1^4 + x2^4 ) // R(1)
    ]

    FF = RationalFunctionField(set)
    compute_groebner!(FF)

#    @test contains_using_groebner( FF, x1*x2 // R(1) )
 #   @test contains_using_groebner( FF, (x1 + x2) // R(1) )

    # note that may fail occasionaly
    @test contains_randomized( FF, (x1*x2) // R(1) )
    @test contains_randomized( FF, (x1 + x2 + x1*x2) // R(1) )

end


text:
"""
Diffusion2D provides 2D diffusion codes (pure 2D and axisymmetric)
"""
module Diffusion2D
export diffusion2D_AxiSymm_step!, diffusion2D_step!, bc2D_x!, bc2D_z!, bc2D_z_bottom!;

using ParallelStencil
using ParallelStencil.FiniteDifferences2D

@init_parallel_stencil(Threads, Float64, 2);    # initialize parallel stencil in 2D

#------------------------------------------------------------------------------------------
# Solve one diffusion timestep in axisymmetric geometry, including latent heat, with spatially variable Rho, Cp and K
@parallel function diffusion2D_AxiSymm_step!(Tnew, T, R, Rc, qr, qz, K, Kr, Kz, Rho, Cp, dt, dr, dz, L, dPhi_dt)
    @all(Kr)    =  @av_xa(K);                                       # average K in r direction
    @all(Kz)    =  @av_ya(K);                                       # average K in z direction
    @all(qr)    =  @all(Rc).*@all(Kr).*@d_xa(T)./dr;                # heatflux in r
    @all(qz)    =            @all(Kz).*@d_ya(T)./dz;                # heatflux in z
    @inn(Tnew)  =  @inn(T) + dt./(@inn(Rho)*@inn(Cp)).*
                                (1.0./@inn(R).*@d_xi(qr)./dr  +     # 2nd derivative in r
                                               @d_yi(qz)./dz) +     # 2nd derivative in z
                                @inn(Rho)*L*@inn(dPhi_dt);          # latent heat

    return
end

# Solve one diffusion timestep in 2D geometry, including latent heat, with spatially variable Rho, Cp and K
@parallel function diffusion2D_step!(Tnew, T, qx, qz, K, Kx, Kz, Rho, Cp, dt, dx, dz, L, dPhi_dt)
    @all(Kx)    =  @av_xa(K);                                       # average K in x direction
    @all(Kz)    =  @av_ya(K);                                       # average K in z direction
    @all(qx)    =  @all(Kx).*@d_xa(T)./dx;                          # heatflux in x
    @all(qz)    =  @all(Kz).*@d_ya(T)./dz;                          # heatflux in z
    @inn(Tnew)  =  @inn(T) + dt./(@inn(Rho)*@inn(Cp)).*
                                 (@d_xi(qx)./dx + @d_yi(qz)./dz) +  # 2nd derivative
                                 @inn(Rho)*L*@inn(dPhi_dt);         # latent heat

    return
end

# Set x- boundary conditions to be zero-flux
@parallel_indices (iy) function bc2D_x!(T::Data.Array) # apply zero flux BC's at the x-side boundaries
    T[1  , iy] = T[2    , iy]
    T[end, iy] = T[end-1, iy]
    return
end

# Set z- boundary conditions to be zero-flux
@parallel_indices (ix) function bc2D_z!(T::Data.Array) # apply zero flux BC's at side boundaries
    T[ix,1 ]    = T[ix, 2    ]
    T[ix, end]  = T[ix, end-1]
    return
end

# Set z- boundary conditions @ bottom to be zero-flux
@parallel_indices (ix) function bc2D_z_bottom!(T::Data.Array) # apply zero flux BC's at side boundaries
    T[ix,1 ]    = T[ix, 2    ]
    return
end

end

"""
Diffusion3D provides 2D diffusion codes (pure 2D and axisymmetric)
"""
module Diffusion3D

# load required julia packages
using ParallelStencil
using ParallelStencil.FiniteDifferences3D

export diffusion3D_step_varK!, bc3D_x!, bc3D_y!

ParallelStencil.@reset_parallel_stencil()       # reset, as we initialized parallel_stencil already above (if we don't do this here, we dont seem to be able to define the functions below)

@init_parallel_stencil(Threads, Float64, 3);    # initialize parallel stencil in 2D

# Solve one diffusion timestep in 3D geometry, including latent heat, with spatially variable Rho, Cp and K
#  Note: needs the 3D stencil routines; hence part is commented
@parallel function diffusion3D_step_varK!(Tnew, T, qx, qy, qz, K, Kx, Ky, Kz, Rho, Cp, dt, dx, dy, dz, L, dPhi_dt)
    @all(Kx)    =  @av_xa(K);                                       # average K in x direction
    @all(Ky)    =  @av_ya(K);                                       # average K in y direction
    @all(Kz)    =  @av_za(K);                                       # average K in z direction
    @all(qx)    =  @all(Kx).*@d_xa(T)./dx;                          # heatflux in x
    @all(qy)    =  @all(Ky).*@d_ya(T)./dy;                          # heatflux in y
    @all(qz)    =  @all(Kz).*@d_za(T)./dz;                          # heatflux in z

    @inn(Tnew)  =  @inn(T) + dt./(@inn(Rho)*@inn(Cp)).*
                 (@d_xi(qx)./dx + @d_yi(qy)./dy + @d_zi(qz)./dz) +  # 2nd derivative
                               @inn(Rho)*L*@inn(dPhi_dt);           # latent heat

    return
end

# Set x- boundary conditions to be zero-flux
@parallel_indices (iy,iz) function bc3D_x!(T) # apply zero flux BC's at the x-side boundaries
    T[1  , iy,iz] = T[2    , iy,iz]
    T[end, iy,iz] = T[end-1, iy,iz]
    return
end

# Set y- boundary conditions to be zero-flux
@parallel_indices (ix,iz) function bc3D_y!(T) # apply zero flux BC's at the y-side boundaries
    T[ix  , 1,iz] = T[ix, 2,iz]
    T[ix, end,iz] = T[ix, end-1,iz]
    return
end

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualNetworkPeeringPropertiesFormat <: SwaggerModel
    allowVirtualNetworkAccess::Any # spec type: Union{ Nothing, Bool } # spec name: allowVirtualNetworkAccess
    allowForwardedTraffic::Any # spec type: Union{ Nothing, Bool } # spec name: allowForwardedTraffic
    allowGatewayTransit::Any # spec type: Union{ Nothing, Bool } # spec name: allowGatewayTransit
    useRemoteGateways::Any # spec type: Union{ Nothing, Bool } # spec name: useRemoteGateways
    remoteVirtualNetwork::Any # spec type: Union{ Nothing, SubResource } # spec name: remoteVirtualNetwork
    remoteAddressSpace::Any # spec type: Union{ Nothing, AddressSpace } # spec name: remoteAddressSpace
    peeringState::Any # spec type: Union{ Nothing, String } # spec name: peeringState
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function VirtualNetworkPeeringPropertiesFormat(;allowVirtualNetworkAccess=nothing, allowForwardedTraffic=nothing, allowGatewayTransit=nothing, useRemoteGateways=nothing, remoteVirtualNetwork=nothing, remoteAddressSpace=nothing, peeringState=nothing, provisioningState=nothing)
        o = new()
        validate_property(VirtualNetworkPeeringPropertiesFormat, Symbol("allowVirtualNetworkAccess"), allowVirtualNetworkAccess)
        setfield!(o, Symbol("allowVirtualNetworkAccess"), allowVirtualNetworkAccess)
        validate_property(VirtualNetworkPeeringPropertiesFormat, Symbol("allowForwardedTraffic"), allowForwardedTraffic)
        setfield!(o, Symbol("allowForwardedTraffic"), allowForwardedTraffic)
        validate_property(VirtualNetworkPeeringPropertiesFormat, Symbol("allowGatewayTransit"), allowGatewayTransit)
        setfield!(o, Symbol("allowGatewayTransit"), allowGatewayTransit)
        validate_property(VirtualNetworkPeeringPropertiesFormat, Symbol("useRemoteGateways"), useRemoteGateways)
        setfield!(o, Symbol("useRemoteGateways"), useRemoteGateways)
        validate_property(VirtualNetworkPeeringPropertiesFormat, Symbol("remoteVirtualNetwork"), remoteVirtualNetwork)
        setfield!(o, Symbol("remoteVirtualNetwork"), remoteVirtualNetwork)
        validate_property(VirtualNetworkPeeringPropertiesFormat, Symbol("remoteAddressSpace"), remoteAddressSpace)
        setfield!(o, Symbol("remoteAddressSpace"), remoteAddressSpace)
        validate_property(VirtualNetworkPeeringPropertiesFormat, Symbol("peeringState"), peeringState)
        setfield!(o, Symbol("peeringState"), peeringState)
        validate_property(VirtualNetworkPeeringPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type VirtualNetworkPeeringPropertiesFormat

const _property_map_VirtualNetworkPeeringPropertiesFormat = Dict{Symbol,Symbol}(Symbol("allowVirtualNetworkAccess")=>Symbol("allowVirtualNetworkAccess"), Symbol("allowForwardedTraffic")=>Symbol("allowForwardedTraffic"), Symbol("allowGatewayTransit")=>Symbol("allowGatewayTransit"), Symbol("useRemoteGateways")=>Symbol("useRemoteGateways"), Symbol("remoteVirtualNetwork")=>Symbol("remoteVirtualNetwork"), Symbol("remoteAddressSpace")=>Symbol("remoteAddressSpace"), Symbol("peeringState")=>Symbol("peeringState"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_VirtualNetworkPeeringPropertiesFormat = Dict{Symbol,String}(Symbol("allowVirtualNetworkAccess")=>"Bool", Symbol("allowForwardedTraffic")=>"Bool", Symbol("allowGatewayTransit")=>"Bool", Symbol("useRemoteGateways")=>"Bool", Symbol("remoteVirtualNetwork")=>"SubResource", Symbol("remoteAddressSpace")=>"AddressSpace", Symbol("peeringState")=>"String", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ VirtualNetworkPeeringPropertiesFormat }) = collect(keys(_property_map_VirtualNetworkPeeringPropertiesFormat))
Swagger.property_type(::Type{ VirtualNetworkPeeringPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualNetworkPeeringPropertiesFormat[name]))}
Swagger.field_name(::Type{ VirtualNetworkPeeringPropertiesFormat }, property_name::Symbol) =  _property_map_VirtualNetworkPeeringPropertiesFormat[property_name]

const _allowed_VirtualNetworkPeeringPropertiesFormat_peeringState = ["Initiated", "Connected", "Disconnected"]

function check_required(o::VirtualNetworkPeeringPropertiesFormat)
    true
end

function validate_property(::Type{ VirtualNetworkPeeringPropertiesFormat }, name::Symbol, val)
    if name === Symbol("peeringState")
        Swagger.validate_param(name, "VirtualNetworkPeeringPropertiesFormat", :enum, val, _allowed_VirtualNetworkPeeringPropertiesFormat_peeringState)
    end
end


text:
# # [Cathepsin Enzyme Reactions](@id enzyme_example)
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/examples/enzymes/enzyme_reactions.ipynb)

using AlgebraicPetri
using Catlab.Programs
using Catlab.Graphics
using Catlab.WiringDiagrams
using Catlab.CategoricalAlgebra

using DifferentialEquations
using Plots

display_uwd(ex, prog="neato") = to_graphviz(ex, box_labels=:name, junction_labels=:variable, graph_attrs=Dict(:overlap => "false"), prog=prog)
ode(x, t) = ODEProblem(vectorfield(x), concentrations(x), t, rates(x));

# ## Define objects and initial conditions

ob(x) = codom(Open([first(x)], LabelledReactionNet{Number,Int}(x), [first(x)])).ob;
K = :K=>33000;
S = :S=>33000;
L = :L=>33000;
Kinact = :Kinact=>0;
Sinact = :Sinact=>0;
Linact = :Linact=>0;
E = :E=>700000;
G = :G=>1300000;

# ## Helper functions for generating primitive Petri net operations

inact(in,on::Number) = begin
  inact = Symbol(first(in), :inact)
  Open(LabelledReactionNet{Number,Int}(unique((in, inact=>0)), ((Symbol(:inact_,first(in)),on),first(in)=>inact)))
end;

bind(in1, in2, on::Number, off::Number) = begin
  out = Symbol(first(in1),first(in2))
  Open(LabelledReactionNet{Number,Int}(unique((in1, in2,out=>0)), ((Symbol(:bind_,first(in1),first(in2)),on),(first(in1),first(in2))=>out),
                                                           ((Symbol(:unbind_,out),off),out=>(first(in1),first(in2)))))
end;

deg(prod1,prod2,on::Number) = begin
  in = Symbol(first(prod1),first(prod2))
  prod2str = String(first(prod2))
  degprod2 = Symbol(endswith(prod2str, "inact") ? first(prod2str) : prod2str, :deg)
  Open(LabelledReactionNet{Number,Int}(unique((in=>0, prod1,degprod2=>0)), ((Symbol(:deg_,in),on),in=>(first(prod1),degprod2))));
end;

# ## Cathepsin *X* reacting with itself

catX = @relation (X, Xinact, Xdeg) where (X, Xinact, Xdeg, XX, XXinact) begin
  inactX(X, Xinact)
  bindXX(X, XX)
  degXX(XX, X, Xdeg)
  bindXXinact(X, Xinact, XXinact)
  degXXinact(XXinact, X, Xdeg)
end
display_uwd(catX)

# ## Cathepsin *X* reacting with Substrate *Y*

catXsubY = @relation (X, Xinact, Xdeg, Y, Ydeg) where (X, Xinact, Xdeg, Y, XY, Ydeg) begin
  bindXY(X, Y, XY)
  degXY(XY, X, Ydeg)
end
display_uwd(catXsubY)

# ## Cathepsin *X* reacting with Cathepsin *Y*

catXY = @relation (X, Xinact, Xdeg, Y, Yinact, Ydeg) where (X, Xinact, Xdeg, Y, Yinact, Ydeg, XY, XYinact) begin
  bindXY(X, Y, XY)
  degXY(XY, X, Ydeg)
  bindXYinact(X, Yinact, XYinact)
  degXYinact(XYinact, X, Ydeg)
end
display_uwd(catXY)

# ## Define all enzyme reactions

rxns = Dict(
  :K => [inact(K, 7.494e-10)
         bind(K, K, 7.814e-4, 3.867e-3)
         deg(K, K, 2.265e-1)
         bind(K, Kinact, 7.814e-4, 3.867e-3)
         deg(K, Kinact, 2.265e-1)],
  :S => [inact(S, 7.494e-10)
         bind(S, S, 7.814e-4, 3.867e-3)
         deg(S, S, 2.265e-1)
         bind(S, Sinact, 7.814e-4, 3.867e-3)
         deg(S, Sinact, 2.265e-1)],
  :L => [inact(L, 7.494e-10)
         bind(L, L, 7.814e-4, 3.867e-3)
         deg(L, L, 2.265e-1)
         bind(L, Linact, 7.814e-4, 3.867e-3)
         deg(L, Linact, 2.265e-1)],
  :KE => [bind(K, E, 9.668e-6, 1e-2)
          deg(K, E, 1.728e0)],
  :KG => [bind(K, G, 2.764e-6, 8.78e-1)
          deg(K, G, 1.502)],
  :SE => [bind(S, E, 4.197e-7, 1.06e-3)
          deg(S, E, 1.384e4)],
  :SG => [bind(S, G, 5.152e-8, 3.894e-3)
          deg(S, G, 8.755e-1)],
  :LE => [bind(L, E, 1.977e-8, 1e-2)
          deg(L, E, 1.066e2)],
  :LG => [bind(L, G, 3.394e-8, 2.365e1)
          deg(L, G, 4.352)],
  :KS => [bind(K, S, 8.822e-4, 4.114e5)
          deg(K, S, 9e-10)
          bind(K, Sinact, 8.822e-4, 4.114e5)
          deg(K, Sinact, 9e-10)],
  :KL => [bind(K, L, 1.756e-4, 3.729e4)
          deg(K, L, 6.505e6)
          bind(K, Linact, 1.756e-4, 3.729e4)
          deg(K, Linact, 6.505e6)],
  :SK => [bind(S, K, 8.822e-4, 4.114e5)
          deg(S, K, 9e-10)
          bind(S, Kinact, 8.822e-4, 4.114e5)
          deg(S, Kinact, 9e-10)],
  :SL => [bind(S, L, 1e-3, 5e2)
          deg(S, L, 1e-7)
          bind(S, Linact, 1e-3, 5e2)
          deg(S, Linact, 1e-7)],
  :LK => [bind(L, K, 1e-3, 4.118e3)
          deg(L, K, 3.234e1)
          bind(L, Kinact, 1e-3, 4.118e3)
          deg(L, Kinact, 3.234e1)],
  :LS => [bind(L, S, 1.056e-12, 5e2)
          deg(L, S, 5e-1)
          bind(L, Sinact, 1.056e-12, 5e2)
          deg(L, Sinact, 5e-1)]
);

# ## Helper functions to generate oapply calls

cat(cat) = begin
  catsym = first(cat)
  out = oapply(catX, Dict([:inactX, :bindXX, :degXX, :bindXXinact, :degXXinact] .=> rxns[catsym]), Dict(
    :X=>ob(cat),
    :Xinact=>ob(Symbol(catsym,:inact)=>0),
    :Xdeg=>ob(Symbol(catsym,:deg)=>0),
    :XX=>ob(Symbol(catsym,catsym)=>0),
    :XXinact=>ob(Symbol(catsym,catsym,:inact)=>0)))
  bundle_legs(out, [[1,2,3]])
end

cat_sub(cat1, sub) = begin
  catsym = first(cat1)
  subsym = first(sub)
  catsub = Symbol(catsym, subsym)
  out = oapply(catXsubY, Dict([:bindXY, :degXY] .=> rxns[catsub]), Dict(
    :X=>ob(cat1),
    :Xinact=>ob(Symbol(catsym,:inact)=>0),
    :Xdeg=>ob(Symbol(catsym,:deg)=>0),
    :Y=>ob(sub),
    :XY=>ob(Symbol(catsym,subsym)=>0),
    :Ydeg=>ob(Symbol(subsym,:deg)=>0)))
  bundle_legs(out, [[1,2,3], [4,5]])
end

cat_cat(cat1, cat2) = begin
  cat1sym = first(cat1)
  cat2sym = first(cat2)
  catcat = Symbol(cat1sym, cat2sym)
  out = oapply(catXY, Dict([:bindXY, :degXY, :bindXYinact, :degXYinact] .=> rxns[catcat]), Dict(
    :X=>ob(cat1),
    :Xinact=>ob(Symbol(cat1sym,:inact)=>0),
    :Xdeg=>ob(Symbol(cat1sym,:deg)=>0),
    :Y=>ob(cat2),
    :Yinact=>ob(Symbol(cat2sym,:inact)=>0),
    :Ydeg=>ob(Symbol(cat2sym,:deg)=>0),
    :XY=>ob(catcat=>0),
    :XYinact=>ob(Symbol(catcat,:inact)=>0)))
  bundle_legs(out, [[1,2,3], [4,5,6]])
end

functor(x) = oapply(x, Dict(
  :catK=>cat(K),
  :catS=>cat(S),
  :catL=>cat(L),
  :catKcatS=>cat_cat(K,S),
  :catKcatL=>cat_cat(K,L),
  :catScatK=>cat_cat(S,K),
  :catScatL=>cat_cat(S,L),
  :catLcatK=>cat_cat(L,K),
  :catLcatS=>cat_cat(L,S),
  :catKsubE=>cat_sub(K,E),
  :catSsubE=>cat_sub(S,E),
  :catLsubE=>cat_sub(L,E),
  :catKsubG=>cat_sub(K,G),
  :catSsubG=>cat_sub(S,G),
  :catLsubG=>cat_sub(L,G)));

function enzyme_uwd(enzymes::Array{Symbol}, substrates::Array{Symbol})
  rel = RelationDiagram{Symbol}(0)

  chemicals = vcat(substrates, enzymes)

  subs = add_junctions!(rel, length(substrates), variable=substrates)
  enzs = add_junctions!(rel, length(enzymes), variable=enzymes)
  nsubs = length(subs)
  nenzs = length(enzs)

  catx = add_parts!(rel, :Box, nenzs, name=[Symbol("cat$i") for i in enzymes])
  add_parts!(rel, :Port, nenzs, junction=enzs, box=catx)

  for x in 1:nenzs
    for y in 1:nenzs
      if y != x
        catxy = add_part!(rel, :Box, name=Symbol("cat$(enzymes[x])cat$(enzymes[y])"))
        add_parts!(rel, :Port, 2, junction=[enzs[x], enzs[y]], box=catxy)
      end
    end
  end

  for x in 1:nenzs
    for y in 1:nsubs
      catxy = add_part!(rel, :Box, name=Symbol("cat$(enzymes[x])sub$(substrates[y])"))
      add_parts!(rel, :Port, 2, junction=[enzs[x], subs[y]], box=catxy)
    end
  end
  add_parts!(rel, :OuterPort, length(chemicals), outer_junction = vcat(subs, enzs))
  rel
end

enzyme_reaction(args...) = enzyme_uwd(args...) |> functor |> apex

# # Manually Defining Models

KSE = @relation (K, S, E) begin
  catK(K)
  catS(S)
  catKcatS(K, S)
  catScatK(S, K)
  catKsubE(K, E)
  catSsubE(S, E)
end
display_uwd(KSE)

# # Defining Models using API

# ## CatK, CatS, Elastin

KSE = enzyme_uwd([:K, :S], [:E])
display_uwd(KSE)
#-

KSE_petri = apex(functor(KSE))
ode_prob = ode(KSE_petri, (0.0, 120.0))
sol = solve(ode_prob)

plot(sol)

#-

plot(sol, lw = 2, ylims = (0, 15000), xlims = (0, 30))

# ## CatK, CatS, CatL, Elastin

KSLE = enzyme_reaction([:K, :S, :L], [:E])
ode_prob = ode(KSLE, (0.0,120.0))
sol = solve(ode_prob)
plot(sol, lw = 1, size = (1066, 600))

#-

plot(sol, ylims = (0, 20000), xlims = (0, 30), lw = 2, size = (1066, 600))

# ## CatK, CatS, CatL, Elastin, Gelatin

KSLEG = enzyme_reaction([:K, :S, :L], [:E, :G])
ode_prob = ode(KSLEG, (0.0,120.0))
sol = solve(ode_prob)
plot(sol, lw = 1, size = (1066, 600))

#-

plot(sol, ylims = (0, 30000), xlims = (0, 50), lw = 2, size = (1066, 600))


text:
using Flux, MLDataPattern, Mill, JsonGrinder, JSON, IterTools, Statistics, BenchmarkTools, ThreadTools
import JsonGrinder: suggestextractor, ExtractCategorical, ExtractDict, ExtractString, MultipleRepresentation, extractscalar
import Mill: mapdata, sparsify, reflectinmodel

###############################################################
# start by loading all samples
###############################################################
samples = open("data/recipes.json","r") do fid
	Vector{Dict}(JSON.parse(read(fid, String)))
end
JSON.print(samples[1],2)

###############################################################
# create schema of the JSON
###############################################################
sch = JsonGrinder.schema(samples)

###############################################################
# create extractor and split it into one for loading targets and
# one for loading data, using custom function to set conditions for using n-gram representation
###############################################################
delete!(sch.childs,:id)

extractor = suggestextractor(sch)

extract_data = ExtractDict(deepcopy(extractor.dict))
extract_target = ExtractDict(deepcopy(extractor.dict))
delete!(extract_target.dict, :ingredients)
delete!(extract_data.dict, :cuisine)
extract_target.dict[:cuisine] = JsonGrinder.ExtractCategorical(keys(sch[:cuisine]))

extract_data(JsonGrinder.sample_synthetic(sch))
###############################################################
# we convert JSONs to Datasets
###############################################################
# advised to use all the samples, this is just speedup to demonstrate functionality
data = tmap(extract_data, samples[1:5_000])
data = reduce(catobs, data)
target = tmap(extract_target, samples[1:5_000])
target = reduce(catobs, target)[:cuisine].data

###############################################################
# 	create the model according to the data
###############################################################
m = reflectinmodel(sch, extract_data,
	k -> Dense(k,20,relu),
	d -> SegmentedMeanMax(d),
	fsm = Dict("" => k -> Dense(k, size(target, 1))),
)

###############################################################
#  train
###############################################################
opt = Flux.Optimise.ADAM()
loss = (x,y) -> Flux.logitcrossentropy(m(x).data, y)
valdata = data[1:1_000],target[:,1:1_000]
data, target = data[1_001:5_000], target[:,1001:5_000]
# for less recourse-chungry training, we use only part of data for trainng, but it is advised to used all, as i following line:
# data, target = data[1001:nobs(data)], target[:,1001:size(target,2)]

cb = () -> println("accuracy = ",mean(Flux.onecold(m(valdata[1]).data) .== Flux.onecold(valdata[2])))
#iterations = 1000
iterations = 100
ps = Flux.params(m)
mean(Flux.onecold(m(data).data) .== Flux.onecold(target))

@info "testing the gradient"
loss(data, target)
gs = gradient(() -> loss(data, target), ps)
Flux.Optimise.update!(opt, ps, gs)

Flux.Optimise.train!(loss, ps, repeatedly(() -> (data, target), 20), opt, cb = Flux.throttle(cb, 2))
# feel free to train for longer period of time, this example is learns only 20 iterations, so it runs fast
# Flux.Optimise.train!(loss, ps, repeatedly(() -> (data, target), 500), opt, cb = Flux.throttle(cb, 10))

#calculate the accuracy
mean(Flux.onecold(m(data).data) .== Flux.onecold(target))

# samples = open("recipes_test.json","r") do fid
# 	Array{Dict}(JSON.parse(readstring(fid)))
# end
# ids = map(s -> s["id"],samples)
# tstdata = @>> samples map(extract_data);
# tstdata = cat(tstdata...);
# tstdata = mapdata(sentence2ngrams,tstdata)
# tstdata = mapdata(i -> sparsify(Float32.(i),0.05),tstdata)
# names = extract_target.other["cuisine"].items
# y = map(i -> names[i],Flux.argmax(m(tstdata)));

# using DataFrames
# CSV.write("cuisine.csv",DataFrame(id = ids,cuisine = y ),delim=',')

# schema can also be created in parallel for better performance, compare:

# single threaded
# @btime JsonGrinder.schema(samples)
# multi threaded
# @btime merge(tmap(x->JsonGrinder.schema(collect(x)), Iterators.partition(samples, 5_000))...)


text:
"""
    perm_pvals(fun, data, nPerms; permFun, funArgs...)

Obtains permutation p-values.

# Arguments

- fun = function that returns a test statistic
- data = RawData object
- nPerms = number of permutations. Defaults to `1000`.

# Keyword arguments

- permFun = function used to permute `Y`. Defaults to `shuffle_rows`
  (shuffles rows of `Y`).
- funArgs = variable keyword arguments to be passed into `fun`

# Value

Tuple
- testStats: 2d array of floats; t-statistics
- pvals: 2d array of floats; permutation p-values

# Some notes

Permutations are computed in parallel when possible.

"""
function perm_pvals(fun::Function, data::RawData, nPerms::Int64=1000;
                    permFun::Function=shuffle_rows, funArgs...)

    # Calculate test statistics for data
    testStats = fun(data; funArgs...)
    # Take the absolute value of the test statistics
    absTestStats= abs.(testStats)

    # Initialze array to store permutation p-values
    pvals = convert(SharedArrays.SharedArray{Float64,2},
                    zeros(size(testStats)))
    # Initialize array to store test statistics based on permuted data
    absPermTestStats = Array{Float64}(undef, size(testStats))
    # Initialize RawData object to store permuted Y
    dataPerm = RawData(Response(zeros(size(get_Y(data)))), data.predictors)

    # For every permutation
    @sync @distributed for i = 1:nPerms
        # Permute Y
        dataPerm.response.Y[:,:] = permFun(get_Y(data))
        # Calculate test statistics based on permuted data
        absPermTestStats[:,:] = abs.(fun(dataPerm; funArgs...))

        # Increment p-values based on permuted test statistics
        for k in 1:length(pvals)
            pvals[k] += (absTestStats[k]<=absPermTestStats[k])/nPerms
        end
    end

    # Return test statistics and permutation p-values
    return testStats, convert(Array{Float64,2}, pvals)
end


text:
using Geodesy, DataFrames, JLD, Base.Dates

MANHATTAN_CENTER = LLA(40.782, -73.9706)

PTIME,DTIME,PLONG,PLAT,DLONG,DLAT,FARE,TIP = 1,2,3,4,5,6,7,8
for j in 9:12
  df = DataFrame(pTime= DateTime[], dTime= DateTime[],pX=Float32[], pY=Float32[], dX=Float32[], dY=Float32[], fare=Float32[], tip=Float32[])
  println("Opening data/reduced_trip_$j.csv")
  f = open("data/reduced_trip_$j.csv")
  for (i, ln) in enumerate(eachline(f))
    if i == 1
      continue
    end
    if i%1000000 == 0
      println(i)
    end
    s = split(ln,",")
    pENU = ENU(LLA(parse(Float32,s[PLAT]),parse(Float32,s[PLONG])),MANHATTAN_CENTER)
    dENU = ENU(LLA(parse(Float32,s[DLAT]),parse(Float32,s[DLONG])),MANHATTAN_CENTER)

    pT = DateTime(s[PTIME], "y-m-d H:M:S")
    dT = DateTime(s[DTIME], "y-m-d H:M:S")

    push!(df,[pT,dT,pENU.east, pENU.north, dENU.east, dENU.north, parse(Float32, s[FARE]), parse(Float32, s[TIP])])
  end
  save("data/reduced_trip_$j.jld","df",df)
end

max_stars_repo_name: semartin/TaxiDataExtraction

max_issues_repo_name: ericytung/TaxiDataTools

max_forks_repo_name: ericytung/TaxiDataTools


text:
"""
Determine the null and alternate hypotheses under z-tests given:

* `ns`: Per-group sample size(s)
* `δ`: Deviation from the null
* `σs`: Per-group standard deviations
"""
function hypotheses(::Type{OneSampleZTest}, n::Real, δ::Real, σ::Real)
    se = sqrt(1 / n) * σ
    null = Normal(0, se)
    alt = Normal(δ, se)
    return null, alt
end

"""
Determine the null and alternate hypotheses under t-tests given:

* `ns`: Per-group sample size(s)
* `δ`: Deviation from the null
* `σs`: Per-group standard deviations
"""
function hypotheses(::Type{EqualVarianceZTest}, n::Real, δ::Real, σ::Real)
    se = sqrt(2 / n) * σ
    null = Normal(0, se)
    alt = Normal(δ, se)
    return null, alt
end

# @doc """
# Determine the values of z that lead to overestimates of the magnitude of δ.
# """ ->
thresholds(alt::Normal) = -mean(alt), mean(alt)