text:
export Linear
export linear

struct Linear
    W::Var
    b::Var
end

"""
    Linear(T::Type, insize::Int, outsize::Int, [init_w=Xavier()], [init_b=Fill(0)])

Linear function (a.k.a. affine transformation).

```math
f(x) = w^{T}x + b
```
where ``w`` is a weight matrix and ``b`` is a bias vector.

```julia
T = Float32
x = Var(rand(T,10,5))
f = Linear(T,10,7)
y = f(x)
```
"""
function Linear(::Type{T}, insize::Int, outsize::Int;
    init_W=Xavier(), init_b=Fill(0)) where T

    W = init_W(T, insize, outsize)
    b = init_b(T, outsize)
    Linear(parameter(W), parameter(b))
end
(f::Linear)(x) = linear(x, f.W, f.b)

function linear(x::Var, W::Var, b::Var)
    T = eltype(x)
    if ndims(x) == 1
        ydata = gemv('T', W.data, x.data)
        addto!(ydata, b.data)
    elseif ndims(x) == 2
        ydata = gemm('T', 'N', W.data, x.data)
        broadcast_addto!(ydata, b.data)
    else
        throw("Invalid ndims of x: $(ndims(x))")
    end
    Var(ydata, ∇linear!, (x,W,b))
end
linear(x::Node, W, b) = Node(linear, (x,W,b))

function ∇linear!(y::Var, x::Var, W::Var, b::Var)
    T = eltype(x)
    if ndims(x) == 1
        gy = reshape(y.grad, length(y.grad), 1)
        xdata = reshape(x.data, length(x.data), 1)
        isnothing(W.grad) || gemm!('N', 'T', T(1), xdata, gy, T(1), W.grad)
        isnothing(x.grad) || gemv!('N', T(1), W.data, y.grad, T(1), x.grad)
        isnothing(b.grad) || addto!(b.grad, y.grad)
    elseif ndims(x) == 2
        isnothing(x.grad) || gemm!('N', 'N', T(1), W.data, y.grad, T(1), x.grad)
        isnothing(W.grad) || gemm!('N', 'T', T(1), x.data, y.grad, T(1), W.grad)
        isnothing(b.grad) || addto!(b.grad, sum(y.grad,dims=2))
    else
        throw("Invalid ndims of x: $(ndims(x))")
    end
end


text:
# Busjle et. al. 2009
# ===================

const _MI_MAT_TYPE = NamedArray{Float64,2,PairwiseListMatrix{Float64,false,Vector{Float64}},NTuple{2,OrderedDict{String,Int}}}

function _buslje09(aln, alphabet::A, clusters, lambda, apc) where A
    mi = mapcolpairfreq!(mutual_information, aln,
                         Counts{Float64,2,A}(ContingencyTable(Float64,Val{2},alphabet)),
                         Val{false}, # use diagonal
                         pseudocounts = AdditiveSmoothing{Float64}(lambda),
                         weights = clusters,
                         diagonalvalue = NaN)
    if apc
        APC!(mi)
    end
    mi
end

"""
`buslje09` takes a MSA or a file and a `FileFormat` as first arguments. It calculates a Z score
and a corrected MI/MIp as described on **Busjle et. al. 2009**.

keyword argument, type, default value and descriptions:

```
  - lambda      Float64   0.05    Low count value
  - clustering  Bool      true    Sequence clustering (Hobohm I)
  - threshold             62      Percent identity threshold for clustering
  - maxgap      Float64   0.5     Maximum fraction of gaps in positions included in calculation
  - apc         Bool      true    Use APC correction (MIp)
  - samples     Int       100     Number of samples for Z-score
  - fixedgaps   Bool      true    Fix gaps positions for the random samples
  - alphabet    ResidueAlphabet UngappedAlphabet()  Residue alphabet to be used
```

This function returns:

```
  - Z score
  - MI or MIp
```
"""
function buslje09(aln::AbstractMatrix{Residue};
                  lambda::Float64 = 0.05,
                  clustering::Bool = true,
                  threshold = 62,
                  maxgap::Float64 = 0.5,
                  apc::Bool = true,
                  samples::Int = 100,
                  alphabet::ResidueAlphabet = UngappedAlphabet(),
                  fixedgaps::Bool = true)::NTuple{2,_MI_MAT_TYPE}
    aln = filtercolumns(aln, gapfraction(aln,1) .<= maxgap)
    clusters = clustering ? hobohmI(aln, threshold) : NoClustering()
    mi = _buslje09(aln, alphabet, clusters, lambda, apc)
    if samples > 0
        rand_mi = Array{PairwiseListMatrix{Float64,false,Vector{Float64}}}(undef, samples)
        zmi = copy(mi)
        residuematrix = getresidues(aln)
        for ns in 1:samples
            shuffle!(residuematrix,1,fixedgaps)
            rand_mi[ns] = getarray(_buslje09(aln, alphabet, clusters, lambda, apc))
        end
        PairwiseListMatrices.zscore!(rand_mi, getarray(zmi))
        return(zmi, mi)
    else
        return(fill!(copy(mi),0.0), mi)
    end
end

function buslje09(filename::String, format::Type{T}; kargs...) where T <: FileFormat
    aln = read(filename, T, AnnotatedMultipleSequenceAlignment, generatemapping=true)
    buslje09(aln; kargs...)
end

# MIToS BLMI: Blosum MI
# =====================

function _BLMI(aln, clusters, alpha, beta, apc, lambda::Float64=0.0)
    mi = mapcolpairfreq!(mutual_information, aln,
                         Probabilities{Float64,2,UngappedAlphabet}(
                            ContingencyTable(Float64,Val{2},UngappedAlphabet())
                         ),
                         Val{false}, # use diagonal
                         pseudocounts = AdditiveSmoothing{Float64}(lambda),
                         weights = clusters,
                         diagonalvalue = NaN,
                         pseudofrequencies=BLOSUM_Pseudofrequencies(alpha,beta))
    if apc
        APC!(mi)
    end
    mi
end

"""
`BLMI` takes a MSA or a file and a `FileFormat` as first arguments. It calculates a Z score
(ZBLMI) and a corrected MI/MIp as described on **Busjle et. al. 2009** but using using
BLOSUM62 pseudo frequencies instead of a fixed pseudocount.

Keyword argument, type, default value and descriptions:

```
  - beta        Float64   8.512   β for BLOSUM62 pseudo frequencies
  - lambda      Float64   0.0     Low count value
  - threshold             62      Percent identity threshold for sequence clustering (Hobohm I)
  - maxgap      Float64   0.5     Maximum fraction of gaps in positions included in calculation
  - apc         Bool      true    Use APC correction (MIp)
  - samples     Int       50      Number of samples for Z-score
  - fixedgaps   Bool      true    Fix gaps positions for the random samples
```

This function returns:

```
  - Z score (ZBLMI)
  - MI or MIp using BLOSUM62 pseudo frequencies (BLMI/BLMIp)
```
"""
function BLMI(aln::AbstractMatrix{Residue};
              beta::Float64 = 8.512,
              threshold = 62,
              maxgap::Float64 = 0.5,
              apc::Bool = true,
              samples::Int = 50,
              fixedgaps::Bool = true,
              lambda::Float64 = 0.0)::NTuple{2,_MI_MAT_TYPE}
    aln = filtercolumns(aln, gapfraction(aln,1) .<= maxgap)
    clusters = hobohmI(aln, threshold)
    numbercl = Float64(nclusters(clusters))
    mi = _BLMI(aln, clusters, numbercl, beta, apc, lambda)
    if samples > 0
        rand_mi = Array{PairwiseListMatrix{Float64,false,Vector{Float64}}}(undef, samples)
        zmi = copy(mi)
        residuematrix = getresidues(aln)
        for ns in 1:samples
            shuffle!(residuematrix,1,fixedgaps)
            rand_mi[ns] = getarray(_BLMI(aln, clusters, numbercl, beta, apc, lambda))
        end
        PairwiseListMatrices.zscore(rand_mi, getarray(zmi))
        return(zmi, mi)
    else
        return(fill!(copy(mi),0.0), mi)
    end
end

function BLMI(filename::String, format::Type{T}; kargs...) where T <: FileFormat
    aln = read(filename, T, AnnotatedMultipleSequenceAlignment, generatemapping=true)
    BLMI(aln; kargs...)
end

ext: jl


text:
module MolenaarAlt
using JuMP
using KNITRO

function molenaar_alt(sub_out_con)

    molenaar = Model(
        optimizer_with_attributes(
            KNITRO.Optimizer,
            "ms_enable" => 1,
            "opttol" => 1E-16,
            "opttolabs" => 1e-16,
            "honorbnds" => 1,
            "ms_maxsolves" => 10,
        ),
    )
    JuMP.set_silent(molenaar)

    michaelis_menten(metabolite, protein, kcat, Km) =
        (kcat * protein * metabolite) / (Km + metabolite)
    register(molenaar, :michaelis_menten, 4, michaelis_menten; autodiff = true)

    # Constraints (used the same for now)
    UB = 100
    LB = 1e-8

    @variables molenaar begin
        # growth rate
        LB <= mu <= UB

        # Concentrations
        LB <= sub_out <= UB # substrate outside concentration
        LB <= sub_in <= UB # subtrate inside concentration
        LB <= ribo <= UB # ribosome concentration
        LB <= lipid <= UB # lipid concentration
        LB <= prec <= UB # precursor concentration
        LB <= transp <= UB # transporter concentration
        LB <= lip_syn <= UB # lipid biosynthesis enzyme concentration

        # new
        LB <= catef <= UB # catabolic ineffecient enzyme
        LB <= metef <= UB # metabolic effecient enzyme

        # Rates
        LB <= v_transp <= UB # transporter rate
        LB <= v_lip <= UB # lipid biosynthesis rate
        LB <= v_ribo <= UB # ribosomale rate

        # new
        LB <= v_catef <= UB # catabolic enzyme rate
        LB <= v_metef <= UB # metabolic enzyme rate

        # Cellular constants
        LB <= beta <= UB # surface to area ratio
        0 <= used_ribo[1:5] <= 1  # ribosomes used for synthesis
    end

    @NLparameters molenaar begin
        kcat_ribo == 3
        Km_ribo == 1
        kcat_transp == 7
        Km_transp == 1
        kcat_metab == 5
        Km_metab == 1
        kcat_lip == 5
        Km_lip == 1
        sub_out == sub_out_con

        #new
        kcat_catef == 10
        Km_catef == 1
        kcat_metef == 5
        Km_metef == 1

        gamma == 0.8
    end

    @NLconstraints molenaar begin
        v_ribo == michaelis_menten(prec, ribo, kcat_ribo, Km_ribo)
        v_transp == michaelis_menten(sub_out, transp, kcat_transp, Km_transp)
        v_lip == michaelis_menten(prec, lip_syn, kcat_lip, Km_lip)

        #new
        v_catef == michaelis_menten(sub_in, catef, kcat_catef, Km_catef)
        v_metef == michaelis_menten(sub_in, metef, kcat_metef, Km_metef)

        beta * (lipid + transp) == 1     # membrane proportions
        sum(used_ribo[i] for i = 1:5) == 1     # ribosome proportions
        catef + metef + ribo + lip_syn <= 1.0   # intracellular density
        transp <= lipid # transporter density

        # ribosome activity producing enzymes
        used_ribo[1] * v_ribo - mu * transp == 0
        used_ribo[2] * v_ribo - mu * ribo == 0
        used_ribo[3] * v_ribo - mu * lip_syn == 0
        used_ribo[4] * v_ribo - mu * catef == 0
        used_ribo[5] * v_ribo - mu * metef == 0

        # metabolic enzyme balances
        (v_transp - v_metef - v_catef) - mu * sub_in == 0
        (v_metef + gamma * v_catef - v_lip - v_ribo) - mu * prec == 0
        v_lip - mu * lipid == 0
    end

    @objective(molenaar, Max, mu)
    optimize!(molenaar)

    res =
        Dict("mu" => value(mu), "met_flux" => value(v_metef), "cat_flux" => value(v_catef))
    return res
end

end #module


text:
@testset "Path" begin
    p, q = synthetic_gradient(SynthSphere(), radius = 50)
    Z = convert_gradient(Path(), p, q)
    @test Z[38,38] ≈ 0.0
    @test Z[1,1] ≈ 0.0
    @test Z[76,76] ≈ 42.71216023490449
end


text:
#----------------------------------------------------------#
# Build the objective function for LaplacianOptModel here  #
#----------------------------------------------------------#

function objective_maximize_algebraic_connectivity(lom::LaplacianOptModel)
    JuMP.@objective(lom.model, Max, lom.variables[:γ_var])
end

function objective_maximize_spanning_tree_cost(lom::LaplacianOptModel)
    n = lom.data["num_nodes"]
    edge_weights = lom.data["edge_weights"]

    JuMP.@objective(lom.model, Max, sum(edge_weights[i,j] * lom.variables[:z_var][i,j] for i=1:(n-1), j=(i+1):n))
end


text:
# --------------------- Atom related States

import Base: promote_rule, *, +, real

const AtomState{T} = ACE.State{(:mu, :mu0, :rr), Tuple{AtomicNumber, AtomicNumber, SVector{3, T}}}
const PopAtomState{T} = ACE.State{(:mu, :mu0, :rr, :pop), Tuple{AtomicNumber, AtomicNumber, SVector{3, T}, T}}

+(X::TX, u::SVector{3}) where {TX <: AtomState} = TX( (rr = X.rr + u, mu = X.mu, mu0 = X.mu0) )
+(X::TX, u::SVector{3}) where {TX <: PopAtomState} = TX( (rr = X.rr + u, mu = X.mu, mu0 = X.mu0, pop = X.pop) )
+(u::SVector{3}, X::TX) where {TX <: ACE.DState{(:rr,)}} = u + X.rr

real(X::AtomState{T}) where {T} =
      AtomState{real(T)}( (rr = real.(X.rr), mu = X.mu, mu0 = X.mu0) )


text:
"""
    Fly dipoles
    ===========

    flydi
        trajectories for a range of starting conditions

    flydi_last
        final position for a range of starting conditions

"""

function flydi(fa::FastAdjust, voltages::Union{Function, Vector{Float64}},
               di::DataFrame,
               dipole::Union{Float64, Function}, mass::Float64, dt::Float64, method::Function;
               max_field::Float64=NaN, max_t::Float64=1.0, max_iterations::Int64=Int64(1e6))
    """ Dipole trajectory
    """
    # fly dipoles
    num, cols = size(di)
    result = Array{DataFrame, 1}(undef, num)
    for i in 1:num
        t0, r0, v0 = di[i, :t], Vector(di[i, [:x, :y, :z]]), Vector(di[i, [:vx, :vy, :vz]])
        result[i] = traj(fa, voltages, t0, r0, v0, dipole, mass, dt, method,
                         max_field=max_field, max_t=max_t, max_iterations=max_iterations, df=true)
    end
    return result
end

function flydi_last(fa::FastAdjust, voltages::Union{Function, Vector{Float64}},
                     di::DataFrame,
                     dipole::Union{Float64, Function}, mass::Float64, dt::Float64, method::Function;
                     max_field::Float64=NaN, max_t::Float64=1.0, max_iterations::Int64=Int64(1e6))
    """ Dipole trajectory (final position)
    """
    # fly dipoles
    t_initial = Array(di[:t])
    r_initial = convert(Matrix, di[[:x, :y, :z]])
    v_initial = convert(Matrix, di[[:vx, :vy, :vz]])
    num, dims = size(r_initial)
    result = Array{Float64, 2}(undef, num, 7)
    for i in 1:num
        t0, r0, v0 = t_initial[i], r_initial[i, :], v_initial[i, :]
        result[i, :] = traj_last(fa, voltages, t0, r0, v0, dipole, mass, dt, method,
                                 max_field=max_field, max_t=max_t, max_iterations=max_iterations)
    end
    result = DataFrame(result)
    rename!(result, [:t, :x, :y, :z, :vx, :vy, :vz]);
    return result
end


text:
export
	DynamicsConstraint,
	integration

############################################################################################
#                              DYNAMICS CONSTRAINTS										   #
############################################################################################

abstract type AbstractDynamicsConstraint{W<:Coupled,P} <: AbstractConstraint{Equality,W,P} end
state_dim(con::AbstractDynamicsConstraint) = size(con.model)[1]
control_dim(con::AbstractDynamicsConstraint) = size(con.model)[2]
Base.length(con::AbstractDynamicsConstraint) = size(con.model)[1]

""" $(TYPEDEF)
An equality constraint imposed by the discretized system dynamics. Links adjacent time steps.
Supports both implicit and explicit integration methods. Can store values internally for
more efficient computation of dynamics and dynamics Jacobians over the entire trajectory,
particularly for explicit methods. These constraints are used in Direct solvers, where
the dynamics are explicit stated as constraints in a more general optimization method.

# Constructors
```julia
DynamicsConstraint{Q}(model::AbstractModel, N)
```
where `N` is the number of knot points and `Q<:QuadratureRule` is the integration method.
"""
struct DynamicsConstraint{Q<:QuadratureRule,L<:AbstractModel,T,N,M,NM} <: AbstractDynamicsConstraint{Coupled,N}
	model::L
    fVal::Vector{SVector{N,T}}
    xMid::Vector{SVector{N,T}}
    ∇f::Vector{SizedMatrix{N,NM,T,2}}
	A::Vector{SubArray{T,2,SizedMatrix{N,NM,T,2},Tuple{UnitRange{Int},UnitRange{Int}},false}}
	B::Vector{SubArray{T,2,SizedMatrix{N,NM,T,2},Tuple{UnitRange{Int},UnitRange{Int}},false}}
	grad::Vector{GradientExpansion{T,N,M}}
	∇fMid::Vector{SizedMatrix{N,NM,T,2}}
	Am::Vector{SubArray{T,2,SizedMatrix{N,NM,T,2},Tuple{UnitRange{Int},UnitRange{Int}},false}}
	Bm::Vector{SubArray{T,2,SizedMatrix{N,NM,T,2},Tuple{UnitRange{Int},UnitRange{Int}},false}}
end

function DynamicsConstraint{Q}(model::L, N) where {Q,L}
	T = Float64  # TODO: get this from somewhere
	n,m = size(model)
	fVal = [@SVector zeros(n) for k = 1:N]
	xMid = [@SVector zeros(n) for k = 1:N]
	∇f   = [SizedMatrix{n,n+m}(zeros(n,n+m)) for k = 1:N]
	∇fm  = [SizedMatrix{n,n+m}(zeros(n,n+m)) for k = 1:3]
	ix,iu = 1:n, n .+ (1:m)
	A  = [view(∇f[k], ix,ix) for k = 1:N]
	B  = [view(∇f[k], ix,iu) for k = 1:N]
	Am = [view(∇fm[k],ix,ix) for k = 1:3]
	Bm = [view(∇fm[k],ix,iu) for k = 1:3]
	NM = n+m
	grad  = [GradientExpansion{T}(n,m) for k = 1:3]
	DynamicsConstraint{Q,L,T,n,m,NM}(model, fVal, xMid, ∇f, A, B,
		grad, ∇fm, Am, Bm)
end

@inline DynamicsConstraint(model, N) = DynamicsConstraint{DEFAULT_Q}(model, N)
integration(::DynamicsConstraint{Q}) where Q = Q

width(con::DynamicsConstraint{<:Implicit,L,T,N,M,NM}) where {L,T,N,M,NM} = 2N+M
width(con::DynamicsConstraint{<:Explicit,L,T,N,M,NM}) where {L,T,N,M,NM} = 2NM
####!

# Implicit
function evaluate!(vals::Vector{<:AbstractVector}, con::DynamicsConstraint{Q},
		Z::Traj, inds=1:length(Z)-1) where Q<:Implicit
	for k in inds
		vals[k] = discrete_dynamics(Q, con.model, Z[k]) - state(Z[k+1])
	end
end

function jacobian!(∇c::Vector{<:SizedMatrix}, con::DynamicsConstraint{Q,L,T,N},
		Z::Vector{<:AbstractKnotPoint{T,n,m}}, inds=1:length(Z)-1) where {Q<:Implicit,L,T,N,n,m}
	In = Diagonal(@SVector ones(N))
	zinds = [Z[1]._x; Z[1]._u]
	for k in inds
		∇f = uview(∇c[k].data, 1:n, 1:n+m+1)
		discrete_jacobian!(Q, ∇f, con.model, Z[k])
		∇f2 = uview(∇c[k].data, 1:n, n+m .+ (1:n))
		∇f2 .= -Diagonal(@SVector ones(n))
	end
end

struct DynamicsVals{T,N,A}
    fVal::Vector{SVector{N,T}}
    xMid::Vector{SVector{N,T}}
    ∇f::Vector{A}
end

function DynamicsVals(dyn_con::DynamicsConstraint)
	DynamicsVals(dyn_con.fVal, dyn_con.xMid, dyn_con.∇f)
end


text:
"""
This is a module for studying neutrino interactions and oscillations in a two flavor model according to the work of Guidry and Billings, 2020.

https://arxiv.org/abs/1812.00035

"""
module neutrinos

# Include plotting support
using PlotlyJS

"""
Fermi's weak coupling constant in MeV cm^2
"""
const fermiConst = 8.955E-44;

"""
The difference in neutrino masses of types one and two in MeV^2
"""
const deltaMSquared = 7.5E-17;

"""
The energy of the neutrinos in MeV
"""
const energy = 10.0;

"""
This function computes the matter mixing angle (theta_m) based on the vacuum mixing angle (theta) and the coupling strength (chi).
"""
function getMatterMixingAngle(vacuumMixingAngle, couplingStrength)

    twoTheta = 2.0*vacuumMixingAngle;
    cosTwoTheta = cos(twoTheta);
    matterMixingAngle = 0.5*atan(sin(twoTheta)/(cos(twoTheta) - couplingStrength));
    scalingFactor = getScalingFactor(vacuumMixingAngle,couplingStrength);
    println("2T = $twoTheta, chi = $couplingStrength, cosT2M = $cosTwoTheta, theta_m = $matterMixingAngle, scalingFactor = $scalingFactor");

    matterMixingAngle = 0.5*asin(sin(twoTheta)/scalingFactor);

    return matterMixingAngle;
end

"""
This function computes the coupling strength (chi) from the neutrino energy (E) and the electron number density (n_e). The energy should be in units of erg.
"""
function getCouplingStrength(energy, e_number_density)

    k = ((2.0*sqrt(2.0)*energy*fermiConst)/deltaMSquared);
    println("$k");
    coupling_strength = k*e_number_density;

    return coupling_strength;
end

"""
This function computes the electron number density as a function of the radial fraction, R/R_sun.
"""
function getElectronNumberDensity(radiusFraction)
    return 1.475E26*exp(-10.54*radiusFraction);
end

"""
This function computes the probability of a neutrino to be in the electron neutrino flavor as a function of its position in the sun (given by the r/r_sun fraction, radiusFraction), and the vacuum mixing angle, (theta).
"""
function getElectronNeutrinoProbability(radiusFraction, theta)
    # Compute the number density
    electronNumberDensity = getElectronNumberDensity(radiusFraction);
    # Compute the coupling strength
    couplingStrength = getCouplingStrength(energy, electronNumberDensity);
    # Compute the matter mixing angle
    matterMixingAngle = getMatterMixingAngle(theta, couplingStrength);
    println("$electronNumberDensity, $couplingStrength, $matterMixingAngle");

    # Compute the probability
    probability = 0.5*(1.0 + (cos(2.0*theta) * cos(2.0*matterMixingAngle)));

    return probability;
end

"""
This function computes the probability of a neutrino to be in the muon neutrino flavor as a function of its position in the sun (given by the r/r_sun fraction, radiusFraction), and the mixing angle, (theta).
"""
function getMuonNeutrinoProbability(radiusFraction, mixingAngle)
    # This is just the opposite of the electron neutrino probability.
    return 1.0 - getElectronNeutrinoProbability(radiusFraction, mixingAngle);
end

"""
This function returns the scaling factor f(theta,chi) for oscillations in matter where theta is the vacuum mixing angle and chi is the coupling strength.
"""
function getScalingFactor(theta, couplingStrength)
    return sqrt(1.0-2.0*couplingStrength*cos(2.0*theta) + couplingStrength*couplingStrength);
end

"""
This function plots the electron number density as a function of the the radial fraction, R/R_sun.
"""
function plotElectronNumberDensity()
    n = 100;
    xValues = 0.0:0.01:n;
    yValues = Array{Float32}(undef,n);
    for i = 1:n
        yValues[i] = log10(getElectronNumberDensity(xValues[i]));
    end
    layout = Layout(;title="Electron Number Density");
    plot(xValues,yValues,layout);

end

"""
This function plots the scaling factor f(theta,chi).
"""
function plotScalingFactor(theta)
    n = 400;
    xValues = 0.0:0.01:n;
    yValues = Array{Float32}(undef,n);
    for i = 1:n
        yValues[i] = getMatterMixingAngle(theta, xValues[i]);
    end
    layout = Layout(;title="Scaling Factor");
    plot(xValues,yValues,layout);
end

"""
This function plots the matter mixing angle as function of vacuum mixing angle and coupling strength.
"""
function plotMatterMixingAngle(theta)
    n = 400;
    xValues = 0.0:0.01:n;
    yValues = Array{Float32}(undef,n);
    for i = 1:n
        yValues[i] = (xValues[i] > 1.0)  ? getMatterMixingAngle(theta, xValues[i]) : pi/2.0 - getMatterMixingAngle(theta, xValues[i]);
    end
    layout = Layout(;title="Matter Mixing Angle");
    plot(xValues,yValues,layout);
end

function plotCouplingStrength()

    values[100];
    for i = 1:100
        #values[i] = #getCouplingStrength(i*neutrinos.energy,)
    end

end

end # module


text:
file_dir = Pkg.dir("D4M")*"/examples/2Apps/1EntityAnalysis/Entity.jld";

using JLD

E = load(file_dir)["E"]

Es = E;
E = logical(E);

Ae = sqIn(E);
Ae
using PyPlot
figure;
spy(Ae);


text:
struct Parameters{T<:AbstractFloat}
    W::Array{Array{T,2},1}
    b::Array{Array{T,2},1}
    layer_sizes::Array{Int,1}
    activation_functions::Array{Function,1}
    loss_function::Function

    function Parameters{T}(params_dict::Dict) where {T <: AbstractFloat}
        layer_sizes = params_dict[:layer_sizes]
        input_layer_size = params_dict[:input_layer_size]
        rng = params_dict[:rng]

        L = length(layer_sizes)
        W = Array{Array{T,2},1}(undef, L)
        b = Array{Array{T,2},1}(undef, L)
        for l = 1:L
            if l ≠ 1
                W[l] =
                    randn(rng, layer_sizes[l], layer_sizes[l-1]) .*
                    sqrt(2 / layer_sizes[l-1])
            else
                W[l] =
                    randn(rng, layer_sizes[l], input_layer_size) .*
                    sqrt(2 / input_layer_size)
            end

            b[l] = zeros(layer_sizes[l], 1)
        end

        params_dict[:W] = W
        params_dict[:b] = b
        params_dict[:activation_functions] = get(
            params_dict,
            :activation_functions,
            push!(Function[relu for l = 1:L-1], logistic),
        )
        params_dict[:loss_function] = get(params_dict, :loss_function, cross_entropy_binary)

        return new(
            params_dict[:W],
            params_dict[:b],
            params_dict[:layer_sizes],
            params_dict[:activation_functions],
            params_dict[:loss_function],
        )
    end
end

function initialize_zeros(W_dims, object_types)
    if object_types == :matrices
        return [zeros(Float32, n, m) for (n, m) in W_dims]
    elseif object_types == :vectors
        return [zeros(Float32, n, 1) for (n, m) in W_dims]
    end
end

mutable struct Cache{T<:AbstractFloat}
    A::Array{Array{T,2},1}
    Z::Array{Array{T,2},1}
    dA::Array{Array{T,2},1}
    dZ::Array{Array{T,2},1}
    dW::Array{Array{T,2},1}
    db::Array{Array{T,2},1}
    v_dW::Array{Array{T,2},1}
    v_db::Array{Array{T,2},1}
    s_dW::Array{Array{T,2},1}
    s_db::Array{Array{T,2},1}

    function Cache{T}(params_dict::Dict) where {T <: AbstractFloat}
        # Create array containing the dimensions of the W matrices
        W_dims = pushfirst!(copy(params_dict[:layer_sizes]), params_dict[:input_layer_size])
        W_dims = [(W_dims[l+1], W_dims[l]) for l = 1:length(W_dims)-1]

        # Initialize the cache variables
        A = initialize_zeros(W_dims, :vectors)
        Z = initialize_zeros(W_dims, :vectors)
        dA = initialize_zeros(W_dims, :vectors)
        dZ = initialize_zeros(W_dims, :vectors)
        dW = initialize_zeros(W_dims, :matrices)
        db = initialize_zeros(W_dims, :vectors)

        # v and s are for the Adam and momentum optimization cache
        v_dW = initialize_zeros(W_dims, :matrices)
        v_db = initialize_zeros(W_dims, :vectors)
        s_dW = initialize_zeros(W_dims, :matrices)
        s_db = initialize_zeros(W_dims, :vectors)

        return new(A, Z, dA, dZ, dW, db, v_dW, v_db, s_dW, s_db)
    end
end

mutable struct HyperParameters{T<:AbstractFloat}
    learning_rate::T
    mini_batch_size::Int
    optimization::Symbol
    λ::T
    # Momentum optimization parameters
    β::T
    # Adam optimization parameters
    β₁::T
    β₂::T
    ϵ::T
    t::Int

    function HyperParameters{T}(hparams_dict::Dict) where {T <: AbstractFloat}

        return new(
            convert(T, get(hparams_dict, :learning_rate, 0.01)),
            get(hparams_dict, :mini_batch_size, 64),
            get(hparams_dict, :optimization, :gd),
            convert(T, get(hparams_dict, :λ, 0.1)),
            convert(T, get(hparams_dict, :β, 0.9)),
            convert(T, get(hparams_dict, :β₁, 0.9)),
            convert(T, get(hparams_dict, :β₂, 0.999)),
            convert(T, get(hparams_dict, :ϵ, 1e-8)),
            get(hparams_dict, :t, 0),
        )
    end
end

struct MiniBatch{T<:AbstractFloat}
    X::Array{T,2}
    Y::Array{T,2}
end

mutable struct NeuralNetwork{T<:AbstractFloat}
    X::Array{T,2}
    Y::Array{T,2}
    m::T
    L::Int
    params::Parameters
    hparams::HyperParameters
    cache::Cache
    rng::MersenneTwister
    J::Array{T,1}
    epoch::Int

    function NeuralNetwork{T}(
        X::Array,
        Y::Array;
        params_dict::Dict = Dict(),
        hparams_dict::Dict = Dict(),
    ) where {T <: AbstractFloat}

        # Explicitly typecast to T to ensure all calcs use the same precision
        X = convert(Matrix{T}, X)
        Y = convert(Matrix{T}, Y)

        # Explicitly typecast to Any dicts to avoid type errors
        params_dict = convert(Dict{Any,Any}, params_dict)
        hparams_dict = convert(Dict{Any,Any}, hparams_dict)

        # Set values for useful variables
        params_dict[:input_layer_size] = size(X, 1)
        params_dict[:layer_sizes] = get(params_dict, :layer_sizes, [3, 1])

        # Create the rng and seed it with the value from params_dict
        rng = MersenneTwister(get(params_dict, :seed, 8))
        params_dict[:rng] = rng

        # Initialize instances of custom types
        params = Parameters{T}(params_dict)
        hparams = HyperParameters{T}(hparams_dict)
        cache = Cache{T}(params_dict)

        # Set starting values of variables
        m = convert(T, size(X, 2))
        L = length(params_dict[:layer_sizes])
        J = []
        epoch = 0

        NN = new{T}(X, Y, m, L, params, hparams, cache, rng, J, epoch)
    end
end


text:
#=
[eegplotter.jl]
Version =  0.022
Author = "William Herrera"
Copyright = "Copyright 2018 William Herrera"
Created = "12 Jan 2018"
Purpose = "EEG file routines viewer example"
=#

using EDFPlus
using DSP
using Plots
import FileIO
pyplot()
using PyPlot

linspace(start, stop, len) = LinRange{Float64}(start, stop, len)

function averagereference(edfh, channels=edfh.mapped_signals)
    data = EDFPlus.signaldata(edfh)
    recs = edfh.datarecords
    siglen = maximum(x->x.smp_per_record, edfh.signalparam)
    avgref = zeros(siglen*recs)
    spans = map(chan->EDFPlus.signalindices(edfh, chan), channels)
    chancount = length(channels)
    for recnum in 1:size(data)[1], channum in 1:chancount
        span = spans[channum]
        addlen = span[2] - span[1]
        addstart = (recnum-1)*siglen+1
        avgref[addstart:addstart+addlen] = data[recnum, span[1]:span[2]][:] ./ chancount
    end
    avgref
end

function eegpages(edfh; channels=collect(3:8), secsperpage=15.0)
    epages = Array{Array{Array{Float64,1},1},1}([])
    fs = samplerate(edfh, channels[1])
    starts = linspace(0.0, edfh.file_duration, div(edfh.file_duration,secsperpage))
    for t1 in starts[1:end-1]
        t2 = t1 + secsperpage
        epage = multichanneltimesegment(edfh, channels, t1, t2, true)
        for i in 1:size(epage)[1]
            epage[i] = highpassfilter(epage[i], fs, 1.5)
            epage[i] = lowpassfilter(epage[i], fs, 24.0)
        end
        push!(epages, epage)
    end
    avgref = averagereference(edfh, collect(1:8))
    spectpage = spectrogram(avgref)
    spectsegment = Plots.contourf(spectpage.time, spectpage.freq, log.(spectpage.power),
                        xaxis=false, yaxis=false, colorbar=false)
    imgname = "tmp.png"
    Plots.savefig(spectsegment, imgname)
    timepoints = linspace(0.0, secsperpage, length(epages[1][1]))
    channelnames = [trim(edfh.signalparam[chan].label) for chan in channels]
    epages, timepoints, channelnames, imgname, edfh.path
end

function plotpage(timepoints, page, ylabels, imgname, title)
    nchannels = size(page)[1]
    plt = Plots.plot(timepoints, page, layout=(nchannels+1,1),
                           xticks=collect(timepoints[1]:1:timepoints[end]),
                           yticks=false, legend=false, linetype=:line)
    for i in 1:nchannels
        Plots.plot!(yaxis=true, xaxis=false, subplot=i, tight_layout=true, ylabel = ylabels[i])
    end
    plt[nchannels][:xaxis][:showaxis] = true
    Plots.plot!(title = title, subplot=1, size=dims)
    img = FileIO.load(imgname)
    plot!(img, aspect_ratio="auto",subplot=nchannels+1)
    plt
end

function vieweeg(filename)
    currentpage = 1
    maxpage = 1
    channelcount = 8
    needupdate = false
    function onclick(event)
        println(event)
        println(event[:x], "  ", event[:y], "   ", event[:canvas][:get_width_height]())
        x = event[:x]
        y = event[:y]
        (width, height) = event[:canvas][:get_width_height]()
        global dims = (width, height)
        yfrac = div(height, channelcount+2)
        ybottom = div(yfrac, 2)
        ytop = yfrac + ybottom
        if y < ytop && y > ybottom && x < width && x > 0
            tmp = Int(round(x * maxpage / width))
            oldpage = currentpage
            currentpage = (tmp < 1) ? 1 : (tmp > maxpage ? maxpage : tmp)
            if oldpage != currentpage
                needupdate = true
            end
        end
    end
    function onkeypress(event)
        println(event)
        println(event[:key])
        (width, height) = event[:canvas][:get_width_height]()
        global dims = (width, height)
        tmp = currentpage
        if event[:key] in ("up", "right", "pageup", "+")
            tmp = currentpage + 1
        elseif event[:key] in ("down", "left", "pagedown", "-")
            tmp = currentpage - 1
        else
            return
        end
        oldpage = currentpage
        currentpage = (tmp < 1) ? 1 : (tmp > maxpage ? maxpage : tmp)
        if oldpage != currentpage
            needupdate = true
        end
    end
    edfh = loadfile(filename)
    epages, timepoints, channelnames, imgname, title = eegpages(edfh)
    maxpage = length(epages)
    plt = plotpage(timepoints, epages[2], channelnames, imgname, title)
    display(plt)
    fig = PyPlot.gcf()
    fig[:canvas][:mpl_connect]("button_press_event", onclick)
    fig[:canvas][:mpl_connect]("key_press_event", onkeypress)
    println("entering loop...")
    while true
        yield()
        if needupdate
            needupdate = false
            newpage = epages[currentpage]
            xdelta = (currentpage - 1) * timepoints[end]
            mins, secs = Int.(floor.(divrem(xdelta, 60)))
            newtitle = replace(title, r".+\/([^\/]+)$" => s"\1")
            newtitle *= " (Page $currentpage of $maxpage: at $mins minutes, $secs seconds)"
            display(plotpage(timepoints .+ xdelta, newpage, channelnames, imgname, newtitle))
            println("updated")
            yield()
        end
        sleep(0.1)
    end
end

global dims = (1200, 400)
const filename = "../test/EDFPlusTestFile.edf"
vieweeg(filename)


text:
@testset "Alternates" begin
   @test collect(Alternates(0)) == [0]
   @test collect(Alternates(1)) == [0, 1, -1 ]
   @test collect(Alternates(2)) == [0, 1, -1, 2, -2 ]
   @test collect(Alternates(3)) == [0, 1, -1, 2, -2, 3, -3 ]
end


text:
"""
Forward Learner (ie Value iteration) with Bayes-factor mediated updates
(via particle filtering)
core

-- TSAS is initialized to Uniform (1/nr_states)
-- Q value at goal state is set to 0.
"""
mutable struct ForwardLearnerBF{TTEstimate}
    nr_states::Int
    nr_actions::Int
    γ::Float64
    initvalue::Float64
    R::Array{Float64, 1}
    QSA::Array{Float64, 2}
    TSAS::Array{Float64, 3}
    Testimate::TTEstimate
    do_debug_print::Bool
end
function ForwardLearnerBF(; nr_states = 10, nr_actions = 4,
    γ = .9, init_R = 0.,
    initvalue = init_R / (1. - γ),
    R = zeros(Float64, nr_states) .+ init_R,
    QSA = zeros(Float64, nr_states, nr_actions) .+ initvalue,
    TSAS = (1.0 / nr_states)*ones(nr_states, nr_actions, nr_states),
    do_debug_print = false,
    nparticles = 3, surpriseprobability = .01, stochasticity = .01, seedlearner = 3,)

    Testimate = TParticleFilterJump(nr_states = nr_states,
                            nr_actions = nr_actions,
                            nparticles = nparticles,
                            surpriseprobability = surpriseprobability,
                            stochasticity = stochasticity,
                            seed = seedlearner)

    ForwardLearnerBF(nr_states, nr_actions, γ, initvalue, R, QSA, TSAS,
        Testimate, do_debug_print)
end
export ForwardLearnerBF

""" Update: main function """
function update!(learner::ForwardLearnerBF, state, action, next_state, reward)

    updatet!(learner.Testimate, state, action, next_state, reward>0)

    for sPrime in 1:learner.nr_states
        learner.TSAS[state, action, sPrime] = copy(learner.Testimate.Ps1a0s0[sPrime][(action, state)])
    end

    learner.R[next_state] = reward #rewards are observed. R is deterministic
    dovalueiteration!(learner)
    if abs(sum(learner.TSAS[state, action, :]) - 1.) > 0.001
        error("Tsas does not sum to 1 in row but to $(sum(TSAS[state,action,:]))")
    end
    if learner.do_debug_print
        do_debug_print_final(learner)
    end
end
export update!


text:
# Julia wrapper for header: /Volumes/BOOTCAMP/programs/compa_libs/iup/include/iup_scintilla.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function IupScintillaOpen()
  ccall( (:IupScintillaOpen, iup), Cvoid, (), )
end

function IupScintilla()
  ccall( (:IupScintilla, iup), Ptr{Cint}, (), )
end


text:
# ---
# title: Figure tweaking
# cover: assets/figure_tweaking.png
# description: Setting figure attributes.
# author: "[Pietro Vertechi](https://github.com/piever)"
# ---

using AlgebraOfGraphics, CairoMakie
set_aog_theme!() #src

# To tweak figure attributes, simply use the `figure` keyword when plotting. For example

df = (x=rand(100), y=rand(100), z=rand(100), c=rand(["a", "b"], 100))
xyc = data(df) * mapping(:x, :y, layout=:c)
layers = linear() + mapping(color=:z)
plt = xyc * layers
fg = draw(
    plt,
    axis=(aspect=1,),
    figure=(figure_padding=10, backgroundcolor=:gray80, resolution=(800, 400))
)

# save cover image #src
mkpath("assets") #src
save("assets/figure_tweaking.png", fg) #src


text:
using PatModules
include("utils.jl")
@once module two
    import ..utils

    function say_bye()
        utils.say("bye")
    end
end


text:
"""
    fetch(t::Vector)

Lorem ipsum dolor sit amet.
"""
function fetch(t::Vector)
  map(Base.fetch, t) #Vectorize fetch
end


text:
@testsuite "vectors" AT->begin
    a = Float32[]
    x = AT(a)
    @test length(x) == 0
    push!(x, 12f0)
    @test length(x) == 1
    @test x[1] == 12f0

    a = Float32[0]
    x = AT(a)
    @test length(x) == 1
    @test length(GPUArrays.buffer(x)) == 1
    push!(x, 12)
    @test length(GPUArrays.buffer(x)) == GPUArrays.grow_dimensions(0, 1, 1)
    resize!(x, 5)
    @test length(x) == 5
    @test length(GPUArrays.buffer(x)) == 5

    resize!(x, 3)
    @test length(x) == 3
    # we don't shrink buffers yet... TODO shrink them... or should we?
    @test length(GPUArrays.buffer(x)) == 5

    x = AT(Array{Float32}(undef, 16))
    reshape!(x, (2, 2, 2, 2))
    @test size(x) == (2, 2, 2, 2)
    x = AT(Array{Float32}(undef, 16))
    y = AT(rand(Float32, 32))
    update!(x, y)
    @test size(x) == (32,)
end


text:
#Constant mean function

"""
    MeanVec <: Mean

Mean function for a 2D GP,
with constant mean along one dimension
```math
m(x) = [β_1, ..., β_n]
```
with constants ``β_1``, ..., ``β_n``.
"""
mutable struct MeanVec <: Mean
    "Constant vector"
    β::Array{Float64,1}
    "x points corresponding to `β`"
    x0::Array{Float64,1}
    "Priors for mean parameters"
    priors::Array
    """
        MeanVec(β::Array{Float64,1})

    Create `MeanVec` with constant vector `β`.
    """
    MeanVec(β::Array{Float64,1}, x0::Array{Float64,1}) = new(β, x0, [])
end

mean(mVec::MeanVec, x::AbstractVector) = mVec.β[ [findfirst(isequal(_x[2]),mVec.x0) for _x in x] ]
function mean(mVec::MeanVec, X::AbstractMatrix)
    mean_vector = zeros(size(X,2))
    for i in 1:size(X,2)
        mean_vector[i] = mVec.β[findfirst(isequal(X[2,i]),mVec.x0)]
    end
    return mean_vector
end

get_params(mVec::MeanVec) = mVec.β
get_param_names(::MeanVec) = [:β]
num_params(mVec::MeanVec) = length(mVec.β)
function set_params!(mVec::MeanVec, hyp::AbstractVector)
    length(hyp) == length(mVec.β) || throw(ArgumentError("ConstVec mean function only has $(mVec.dim) parameters"))
    copyto!(mVec.β, hyp)
    copyto!(mVec.x0, mVec.x0)
end
function grad_mean(mVec::MeanVec, x::AbstractVector)
    idx = findfirst(isequal(x[2]),mVec.x0)
    dM_theta = zeros(length(mVec.β))
    dM_theta[idx] = 1.0
    return dM_theta
end


text:
import Measurements: Measurement
import ..UncertainValues: UncertainValue
import Distributions: Normal

resample(m::Measurement{T}) where T = resample(UncertainValue(Normal, m.val, m.err))

function resample(m::Measurement{T}, n::Int) where T
    uval = UncertainValue(Normal, m.val, m.err)

    [resample(uval) for i = 1:n]
end


text:
"""
```julia
lorenz(u0=[0.0, 10.0, 0.0]; σ = 10.0, ρ = 28.0, β = 8/3) -> ds
```
```math
\\begin{aligned}
\\dot{X} &= \\sigma(Y-X) \\\\
\\dot{Y} &= -XZ + \\rho X -Y \\\\
\\dot{Z} &= XY - \\beta Z
\\end{aligned}
```
The famous three dimensional system due to Lorenz [1], shown to exhibit
so-called "deterministic nonperiodic flow". It was originally invented to study a
simplified form of atmospheric convection.

Currently, it is most famous for its strange attractor (occuring at the default
parameters), which resembles a butterfly. For the same reason it is
also associated with the term "butterfly effect" (a term which Lorenz himself disliked)
even though the effect applies generally to dynamical systems.
Default values are the ones used in the original paper.

The parameter container has the parameters in the same order as stated in this
function's documentation string.

[1] : E. N. Lorenz, J. atmos. Sci. **20**, pp 130 (1963)
"""
function lorenz(u0=[0.0, 10.0, 0.0]; σ = 10.0, ρ = 28.0, β = 8/3)
    return CDS(loop, u0, [σ, ρ, β], loop_jac)
end
function loop(u, p, t)
    @inbounds begin
        σ = p[1]; ρ = p[2]; β = p[3]
        du1 = σ*(u[2]-u[1])
        du2 = u[1]*(ρ-u[3]) - u[2]
        du3 = u[1]*u[2] - β*u[3]
        return SVector{3}(du1, du2, du3)
    end
end
function loop_jac(u, p, t)
    @inbounds begin
        σ, ρ, β = p
        J = @SMatrix [-σ  σ  0;
        ρ - u[3]  (-1)  (-u[1]);
        u[2]   u[1]  -β]
        return J
    end
end

function lorenz_iip(u0=[0.0, 10.0, 0.0]; σ = 10.0, ρ = 28.0, β = 8/3)
    return CDS(liip, u0, [σ, ρ, β], liip_jac)
end
function liip(du, u, p, t)
    @inbounds begin
        σ = p[1]; ρ = p[2]; β = p[3]
        du[1] = σ*(u[2]-u[1])
        du[2] = u[1]*(ρ-u[3]) - u[2]
        du[3] = u[1]*u[2] - β*u[3]
        return nothing
    end
end
function liip_jac(J, u, p, t)
    @inbounds begin
    σ, ρ, β = p
    J[1,1] = -σ; J[1, 2] = σ; J[1,3] = 0
    J[2,1] = ρ - u[3]; J[2,2] = -1; J[2,3] = -u[1]
    J[3,1] = u[2]; J[3,2] = u[1]; J[3,3] = -β
    return nothing
    end
end

"""
```julia
roessler(u0=rand(3); a = 0.2, b = 0.2, c = 5.7)
```
```math
\\begin{aligned}
\\dot{x} &= -y-z \\\\
\\dot{y} &= x+ay \\\\
\\dot{z} &= b + z(x-c)
\\end{aligned}
```
This three-dimensional continuous system is due to Rössler [1].
It is a system that by design behaves similarly
to the `lorenz` system and displays a (fractal)
strange attractor. However, it is easier to analyze qualitatively, as for example
the attractor is composed of a single manifold.
Default values are the same as the original paper.

The parameter container has the parameters in the same order as stated in this
function's documentation string.

[1] : O. E. Rössler, Phys. Lett. **57A**, pp 397 (1976)
"""
function roessler(u0=rand(3); a = 0.2, b = 0.2, c = 5.7)
    return CDS(roessler_eom, u0, [a, b, c], roessler_jacob)
end
function roessler_eom(u, p, t)
    @inbounds begin
    a, b, c = p
    du1 = -u[2]-u[3]
    du2 = u[1] + a*u[2]
    du3 = b + u[3]*(u[1] - c)
    return SVector{3, Float64}(du1, du2, du3)
    end
end
function roessler_jacob(u, p, t)
    a, b, c = p
    return @SMatrix [0.0 (-1.0) (-1.0);
                     1.0 a 0.0;
                     u[3] 0.0 (u[1]-c)]
end

"""
    double_pendulum(u0 = [π/2, 0, 0, rand()];
                    G=10.0, L1 = 1.0, L2 = 1.0, M1 = 1.0, M2 = 1.0)
Famous chaotic double pendulum system (also used for our logo!). Keywords
are gravity (G), lengths of each rod and mass of each ball (all assumed SI units).

The variables order is [θ1, dθ1/dt, θ2, dθ2/dt].

Jacobian is created automatically (thus methods that use the Jacobian will be slower)!

(please contribute the Jacobian and the e.o.m. in LaTeX :smile:)

The parameter container has the parameters in the same order as stated in this
function's documentation string.
"""
function double_pendulum(u0=[π/2, 0, 0, rand()]; G=10.0, L1 = 1.0, L2 = 1.0, M1 = 1.0, M2 = 1.0)
    return CDS(doublependulum_eom, u0, [G, L1, L2, M1, M2])
end
@inbounds function doublependulum_eom(state, p, t)
    G, L1, L2, M1, M2 = p

    du1 = state[2]
    del_ = state[3] - state[1]
    den1 = (M1 + M2)*L1 - M2*L1*cos(del_)*cos(del_)
    du2 = (M2*L1*state[2]*state[2]*sin(del_)*cos(del_) +
               M2*G*sin(state[3])*cos(del_) +
               M2*L2*state[4]*state[4]*sin(del_) -
               (M1 + M2)*G*sin(state[1]))/den1

    du3 = state[4]

    den2 = (L2/L1)*den1
    du4 = (-M2*L2*state[4]*state[4]*sin(del_)*cos(del_) +
               (M1 + M2)*G*sin(state[1])*cos(del_) -
               (M1 + M2)*L1*state[2]*state[2]*sin(del_) -
               (M1 + M2)*G*sin(state[3]))/den2
    return SVector{4}(du1, du2, du3, du4)
end

"""
    henonheiles(u0=[0, -0.25, 0.42081,0])
```math
\\begin{aligned}
\\dot{x} &= p_x \\\\
\\dot{y} &= p_y \\\\
\\dot{p}_x &= -x -2 xy \\\\
\\dot{p}_y &= -y - (x^2 - y^2)
\\end{aligned}
```

The Hénon–Heiles system [1] was introduced as a simplification of the motion
of a star around a galactic center. It was originally intended to study the
existence of a "third integral of motion" (which would make this 4D system integrable).
In that search, the authors encountered chaos, as the third integral existed
for only but a few initial conditions.

The default initial condition is a typical chaotic orbit.

[1] : Hénon, M. & Heiles, C., The Astronomical Journal **69**, pp 73–79 (1964)
"""
function henonheiles(u0=[0, -0.25, 0.42081, 0]#=; conserveE::Bool = true=#)

    i = one(eltype(u0))
    o = zero(eltype(u0))
    J = zeros(eltype(u0), 4, 4)

    # Vhh(q1, q2) = 1//2 * (q1^2 + q2^2 + 2q1^2 * q2 - 2//3 * q2^3)
    # Thh(p1, p2) = 1//2 * (p1^2 + p2^2)
    # Hhh(q1, q2, p1, p2) = Thh(p1, p2) + Vhh(q1, q2)
    # Hhh(u::AbstractVector) = Hhh(u...)
    #
    # E = Hhh(u0)
    #
    # ghh! = (resid, u) -> begin
    #     resid[1] = Hhh(u[1],u[2],u[3],u[4]) - E
    #     resid[2:4] .= 0
    # end

    # if conserveE
    #     cb = ManifoldProjection(ghh!, nlopts=Dict(:ftol=>1e-13), save = false)
    #     prob = ODEProblem(hheom!, u0, (0., 100.0),  callback=cb)
    # else
        # prob = ODEProblem(hheom!, u0, (0., 100.0))
    # end
    return CDS(hheom!, u0, nothing, hhjacob!, J)
end
function hheom!(du, u, p, t)
    @inbounds begin
        du[1] = u[3]
        du[2] = u[4]
        du[3] = -u[1] - 2u[1]*u[2]
        du[4] = -u[2] - (u[1]^2 - u[2]^2)
        return nothing
    end
end
function hhjacob!(J, u, p, t)
    @inbounds begin
    o = 0.0; i = 1.0
    J[1,:] .= (o,    o,     i,    o)
    J[2,:] .= (o,    o,     o,    i)
    J[3,:] .= (-i - 2*u[2],   -2*u[1],   o,   o)
    J[4,:] .= (-2*u[1],  -1 + 2*u[2],  o,   o)
    return nothing
    end
end

"""
    lorenz96(N::Int, u0 = rand(M); F=0.01)

```math
\\frac{dx_i}{dt} = (x_{i+1}-x_{i-2})x_{i-1} - x_i + F
```

`N` is the chain length, `F` the forcing. Jacobian is created automatically.
(parameter container only contains `F`)
"""
function lorenz96(N::Int, u0 = rand(N); F=0.01)
    @assert N ≥ 4 "`N` must be at least 4"
    lor96 = Lorenz96{N}() # create struct
    return CDS(lor96, u0, [F])
end
struct Lorenz96{N} end # Structure for size type
function (obj::Lorenz96{N})(dx, x, p, t) where {N}
    F = p[1]
    # 3 edge cases
    @inbounds dx[1] = (x[2] - x[N - 1]) * x[N] - x[1] + F
    @inbounds dx[2] = (x[3] - x[N]) * x[1] - x[2] + F
    @inbounds dx[N] = (x[1] - x[N - 2]) * x[N - 1] - x[N] + F
    # then the general case
    for n in 3:(N - 1)
      @inbounds dx[n] = (x[n + 1] - x[n - 2]) * x[n - 1] - x[n] + F
    end
    return nothing
end

"""
    duffing(u0 = [rand(), rand(), 0]; ω = 2.2, f = 27.0, d = 0.2, β = 1)
The (forced) duffing oscillator, that satisfies the equation
```math
\\ddot{x} + d\\cdot\\dot{x} + β*x + x^3 = f\\cos(\\omega t)
```
with `f, ω` the forcing strength and frequency and `d` the dampening.

The parameter container has the parameters in the same order as stated in this
function's documentation string.
"""
function duffing(u0 = [rand(), rand()]; ω = 2.2, f = 27.0, d = 0.2, β = 1)

    J = zeros(eltype(u0), 2, 2)
    J[1,2] = 1
    return CDS(duffing_eom, u0, [ω, f, d, β], duffing_jacob)
end
@inbounds function duffing_eom(x, p, t)
    ω, f, d, β = p
    dx1 = x[2]
    dx2 = f*cos(ω*t) - β*x[1] - x[1]^3 - d * x[2]
    return SVector{2, Float64}(dx1, dx2)
end
@inbounds function duffing_jacob(u, p, t)
    ω, f, d, β = p
    return @SMatrix [0 1 ;
    (-β - 3u[1]^2) -d]
end

"""
    shinriki(u0 = [-2, 0, 0.2]; R1 = 22.0)
Shinriki oscillator with all other parameters (besides `R1`) set to constants.
*This is a stiff problem, be careful when choosing solvers and tolerances*.
"""
function shinriki(u0 = [-2, 0, 0.2]; R1 = 22.0)
    # # Jacobian caller for Shinriki:
    # shinriki_eom(::Type{Val{:jac}}, J, u, p, t) = (shi::Shinriki)(t, u, J)
    return CDS(shinriki_eom, u0, [R1])
end
shinriki_voltage(V) = 2.295e-5*(exp(3.0038*V) - exp(-3.0038*V))
function shinriki_eom(u, p, t)
    R1 = p[1]

    du1 = (1/0.01)*(
    u[1]*(1/6.9 - 1/R1) - shinriki_voltage(u[1] - u[2]) - (u[1] - u[2])/14.5
    )

    du2 = (1/0.1)*(
    shinriki_voltage(u[1] - u[2]) + (u[1] - u[2])/14.5 - u[3]
    )

    du3 = (1/0.32)*(-u[3]*0.1 + u[2])
    return SVector{3}(du1, du2, du3)
end

"""
```julia
gissinger(u0 = 3rand(3); μ = 0.119, ν = 0.1, Γ = 0.9)
```
```math
\\begin{aligned}
\\dot{Q} &= \\mu Q - VD \\\\
\\dot{D} &= -\\nu D + VQ \\\\
\\dot{V} &= \\Gamma -V + QD
\\end{aligned}
```
A continuous system that models chaotic reversals due to Gissinger [1], applied
to study the reversals of the magnetic field of the Earth.

The parameter container has the parameters in the same order as stated in this
function's documentation string.

[1] : C. Gissinger, Eur. Phys. J. B **85**, 4, pp 1-12 (2012)
"""
function gissinger(u0 = 3rand(3); μ = 0.119, ν = 0.1, Γ = 0.9)
    return CDS(gissinger_eom, u0, [μ, ν, Γ])
end
function gissinger_eom(u, p, t)
    μ, ν, Γ = p
    du1 = μ*u[1] - u[2]*u[3]
    du2 = -ν*u[2] + u[1]*u[3]
    du3 = Γ - u[3] + u[1]*u[2]
    return SVector{3}(du1, du2, du3)
end

"""
```julia
rikitake(u0 = [1, 0, 0.6]; μ = 1.0, α = 1.0)
```
```math
\\begin{aligned}
\\dot{x} &= -\\mu x +yz \\\\
\\dot{y} &= -\\mu y +x(z-\\alpha) \\\\
\\dot{V} &= 1 - xz
\\end{aligned}
```
Rikitake's dynamo is a system that tries to model the magnetic reversal events
by means of a double-disk dynamo system.

[1] : T. Rikitake Math. Proc. Camb. Phil. Soc. **54**, pp 89–105, (1958)
"""
function rikitake(u0 = [1, 0, 0.6]; μ = 1.0, α = 1.0)
    return CDS(rikitake_eom, u0, [μ, α])
end
function rikitake_eom(u, p, t)
    μ, α = p
    x,y,z = u
    xdot = -μ*x + y*z
    ydot = -μ*y + x*(z - α)
    zdot = 1 - x*y
    return SVector{3}(xdot, ydot, zdot)
end

"""
```julia
nosehoover(u0 = [0, 0.1, 0])
```
```math
\\begin{aligned}
\\dot{x} &= y \\\\
\\dot{y} &= yz - x \\\\
\\dot{V} &= 1 - y^2
\\end{aligned}
```
Three dimensional conservative continuous system, taken from the book
"Elegant Chaos" by J. C. Sprott.
"""
nosehoover(u0 = [0, 0.1, 0]) = CDS(nosehoover_eom, u0, nothing)
function nosehoover_eom(u, p, t)
    x,y,z = u
    xdot = y
    ydot = y*z - x
    zdot  = 1.0 - y*y
    return SVector{3}(xdot, ydot, zdot)
end

"""
```julia
labyrinth(u0 = [1.0, 0, 0])
```
```math
\\begin{aligned}
\\dot{x} &= \\sin(y) \\\\
\\dot{y} &= \\sin(z) \\\\
\\dot{V} &= \\sin(x)
\\end{aligned}
```
Three dimensional conservative continuous system, whose evolution in 3D space looks
like a speudo-random walk, the orbit moving around like in a labyrinth.
Taken from the book "Elegant Chaos" by J. C. Sprott.
"""
labyrinth(u0 = [1.0, 0, 0]) = CDS(labyrinth_eom, u0, nothing)
function labyrinth_eom(u, p, t)
    x,y,z = u
    xdot = sin(y)
    ydot = sin(z)
    zdot = sin(x)
    return SVector{3}(xdot, ydot, zdot)
end


text:
"""
`AUC(scores::PairwiseListMatrix, msacontacts::PairwiseListMatrix)`

Returns the Area Under a ROC (Receiver Operating Characteristic) Curve (AUC) of the
`scores` for `msacontact` prediction. `score` and `msacontact` lists are vinculated
(inner join) by their labels (i.e. column number in the file). `msacontact` should have
1.0 for true contacts and 0.0 for not contacts (NaN or other numbers for missing
values). You need to do `using ROCAnalysis` before using this function.
"""
function ROCAnalysis.AUC(scores::NamedArray{L,2,PairwiseListMatrix{L,false,VL},NL},
    msacontacts::NamedArray{R,2,PairwiseListMatrix{L,false,VR},NR}) where {L <: AbstractFloat,
        R <: AbstractFloat,VL,VR,NL,NR}
    sco, con = join(scores, msacontacts, kind=:inner)
    true_contacts, false_contacts = getcontactmasks(con)
    scores_list = getlist(getarray(sco))
    ROCAnalysis.AUC(
        ROCAnalysis.roc(
            scores_list[true_contacts  .& .!isnan.(scores_list)],
            scores_list[false_contacts .& .!isnan.(scores_list)]))
end

export AUC


text:
using LinearAlgebra
using StaticArrays
using ParametricSurfaces
using Random
using Nystrom
Random.seed!(1)

# parameters
k   = 4π
λ   = 2π/k
ppw = 16
dx  = λ/ppw

# select pde
pdes = [Laplace(dim=3),
        Helmholtz(dim=3,k=k),
        Elastostatic(dim=3,μ=1,λ=2),
        Maxwell(dim=3,k=k)]
pde = pdes[2]

# select operator
ops = [SingleLayerOperator,
       DoubleLayerOperator]
op = ops[1]

# geometry
clear_entities!()
geo = ParametricSurfaces.Sphere(;radius=1)
Ω   = Domain(geo)
Γ   = boundary(Ω)
np  = ceil(2/dx)
M   = meshgen(Γ,(np,np))
msh = NystromMesh(M,Γ;order=1)
nx = length(msh.dofs);
iop = op(pde,msh)

# full product estimation
T = Nystrom.default_density_eltype(pde)
I = rand(1:nx,1000)
B = rand(T,nx) # CHECK: randn or rand?
Xpts = msh.dofs
Ypts = msh.dofs
K = (x,y) -> iop.kernel(x,y)*Nystrom.weight(y)
tfull = @elapsed exa = [sum(K(Xpts[i],Ypts[j])*B[j] for j in 1:nx) for i in I]
@info "Estimated time for full product: $(tfull*nx/1000)"

# IFGF
nmax = 100   # max points per leaf box
p = (3,5,5)  # interpolation points per dimension
compress = Nystrom.IFGFCompressor(;p,nmax,_profile=true);
A = compress(iop);
C = zeros(T,nx)
Nystrom.IFGF.@hprofile C = A*B
er = norm(C[I]-exa,2) / norm(exa,2)
@info "" typeof(pde) op (er,nx)

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
module Main

export main

function main()

  input = """1163751742
  1381373672
  2136511328
  3694931569
  7463417111
  1319128137
  1359912421
  3125421639
  1293138521
  2311944581"""
  ints = reshape(map(x -> parse(Int, x), collect(replace(input, '\n' => ""))), (10, 10))
  println(ints)

  function next_pos(curr_pos::Vector{Int}, height::Int, width::Int)::Vector{Vector{Int}}
    offsets = [[0, 1], [1, 0], [1, 1], [-1, -1]]
    neighbors = map(.+, offsets, curr_pos)
    is_within_bounds(x) = all(x .> 0) && x[1] < height && x[2] < width
    directions = findall(is_within_bounds, neighbors)
    offsets[directions], directions
  end

  # x = next_pos([1, 1], 10, 10)
  # println(x)

  function dfs()
    visited = Dict()
    deque = []
    push!(deque, ([1, 1], "down", ints[[1, 1]]))
    push!(deque, ([1, 1], "right", ints[[1, 1]]))
    while true
      curr_pos, direction, count = pop!(deque)
      if (curr_pos, direction) in keys(visited)
        if count < visited((curr_pos, direction))
          visited((curr_pos, direction)) = count
        end
      end
      visited[(curr_pos, direction)] = count
      deque = [next_pos; deque]
    end
  end

end

end

main()

max_stars_repo_name: acjensen/advent-of-code

max_issues_repo_name: acjensen/advent-of-code

max_forks_repo_name: acjensen/advent-of-code


text:
## mapreduce implementation that skips missing

function skipmissing_init(f, op, na::BitArray, data::Array, ifirst::Int, ilast::Int)
    # Get first non-missing element
    ifirst = Base.findnextnot(na, ifirst)
    @inbounds d1 = data[ifirst]

    # Get next non-missing element
    ifirst = Base.findnextnot(na, ifirst+1)
    @inbounds d2 = data[ifirst]

    # Reduce first two elements
    (op(f(d1), f(d2)), ifirst)
end

function mapreduce_seq_impl_skipmissing(f, op, T, A::DataArray, ifirst::Int, ilast::Int)
    data = A.data
    na = A.na
    chunks = na.chunks

    v, i = skipmissing_init(f, op, na, data, ifirst, ilast)

    while i < ilast
        i += 1
        @inbounds na = Base.unsafe_bitgetindex(chunks, i)
        na && continue
        @inbounds d = data[i]
        v = op(v, f(d))
    end
    v
end

# Pairwise map-reduce
function mapreduce_pairwise_impl_skipmissing(f, op, A::DataArray{T}, bytefirst::Int, bytelast::Int, n_notna::Int, blksize::Int) where T
    if n_notna <= blksize
        ifirst = 64*(bytefirst-1)+1
        ilast = min(64*bytelast, length(A))
        # Fall back to Base implementation if no missings in block
        return ilast - ifirst + 1 == n_notna ? Base.mapreduce_seq_impl(f, op, A.data, ifirst, ilast) :
                                               mapreduce_seq_impl_skipmissing(f, op, T, A, ifirst, ilast)
    end

    # Find byte in the middle of range
    # The block size is restricted so that there will always be at
    # least two non-missing elements in the returned range
    chunks = A.na.chunks
    nmid = 0
    imid = bytefirst-1
    while nmid < (n_notna >> 1)
        imid += 1
        @inbounds nmid += count_zeros(chunks[imid])
    end

    v1 = mapreduce_pairwise_impl_skipmissing(f, op, A, bytefirst, imid, nmid, blksize)
    v2 = mapreduce_pairwise_impl_skipmissing(f, op, A, imid+1, bytelast, n_notna-nmid, blksize)
    op(v1, v2)
end

if isdefined(Base, :pairwise_blocksize)
    sum_pairwise_blocksize(f) = Base.pairwise_blocksize(f, +)
else
    const sum_pairwise_blocksize = Base.sum_pairwise_blocksize
end

mapreduce_impl_skipmissing(f, op, A::DataArray{T}) where {T} =
    mapreduce_seq_impl_skipmissing(f, op, T, A, 1, length(A.data))
mapreduce_impl_skipmissing(f, op::typeof(+), A::DataArray) =
    mapreduce_pairwise_impl_skipmissing(f, op, A, 1, length(A.na.chunks),
                                     length(A.na)-countnz(A.na),
                                     max(128, sum_pairwise_blocksize(f)))

## general mapreduce interface

function _mapreduce_skipmissing(f, op, A::DataArray{T}) where T
    n = length(A)
    na = A.na

    nna = countnz(na)
    nna == n && return Base.mr_empty(f, op, T)
    nna == n-1 && return Base.r_promote(op, f(A.data[Base.findnextnot(na, 1)]))
    nna == 0 && return Base.mapreduce_impl(f, op, A.data, 1, n)

    mapreduce_impl_skipmissing(f, op, A)
end

# This is only safe when we can guarantee that if a function is passed
# missing, it returns missing. Otherwise we will fall back to the implementation
# in Base, which is slow because it's type-unstable, but guarantees the
# correct semantics
const SafeMapFuns = Union{typeof(identity), typeof(abs), typeof(abs2),
                            typeof(exp), typeof(log), typeof(Base.centralizedabs2fun)}
const SafeReduceFuns = Union{typeof(+), typeof(*), typeof(max), typeof(min)}
function Base._mapreduce(f::SafeMapFuns, op::SafeReduceFuns, A::DataArray)
    any(A.na) && return missing
    Base._mapreduce(f, op, A.data)
end

function Base.mapreduce(f, op::Function, A::DataArray;
                        skipmissing::Bool=false, skipna::Bool=false)
    if skipna && !skipmissing
        Base.depwarn("skipna=$skipna is deprecated, use skipmissing=$skipna instead", :mapreduce)
        skipmissing = true
    end
    (op === +) ? (skipmissing ? _mapreduce_skipmissing(f, +, A) : Base._mapreduce(f, +, A)) :
    (op === *) ? (skipmissing ? _mapreduce_skipmissing(f, *, A) : Base._mapreduce(f, *, A)) :
    (op === &) ? (skipmissing ? _mapreduce_skipmissing(f, &, A) : Base._mapreduce(f, &, A)) :
    (op === |) ? (skipmissing ? _mapreduce_skipmissing(f, |, A) : Base._mapreduce(f, |, A)) :
    skipmissing ? _mapreduce_skipmissing(f, op, A) : Base._mapreduce(f, op, A)
end

# To silence deprecations, but could be more efficient
function Base.mapreduce(f, op::Union{typeof(|), typeof(&)}, A::DataArray;
                        skipmissing::Bool=false, skipna::Bool=false)
    if skipna && !skipmissing
        Base.depwarn("skipna=$skipna is deprecated, use skipmissing=$skipna instead", :mapreduce)
        skipmissing = true
    end
    skipmissing ? _mapreduce_skipmissing(f, op, A) : Base._mapreduce(f, op, A)
end

function Base.mapreduce(f, op, A::DataArray;
                        skipmissing::Bool=false, skipna::Bool=false)
    if skipna && !skipmissing
        Base.depwarn("skipna=$skipna is deprecated, use skipmissing=$skipna instead", :mapreduce)
        skipmissing = true
    end
    skipmissing ? _mapreduce_skipmissing(f, op, A) : Base._mapreduce(f, op, A)
end

Base.reduce(op, A::DataArray;
            skipmissing::Bool=false, skipna::Bool=false) =
    mapreduce(identity, op, A; skipmissing=skipmissing, skipna=skipna)

## usual reductions

for (fn, op) in ((:(Base.sum), +),
                 (:(Base.prod), *),
                 (:(Base.minimum), min),
                 (:(Base.maximum), max))
    @eval begin
        $fn(f::Union{Function,$(supertype(typeof(abs)))}, a::DataArray;
            skipmissing::Bool=false, skipna::Bool=false) =
            mapreduce(f, $op, a; skipmissing=skipmissing, skipna=skipna)
        $fn(a::DataArray; skipmissing::Bool=false, skipna::Bool=false) =
            mapreduce(identity, $op, a; skipmissing=skipmissing, skipna=skipna)
    end
end

## mean

Base.mean(a::DataArray; skipmissing::Bool=false, skipna::Bool=false) =
    sum(a; skipmissing=skipmissing, skipna=skipna) / (length(a.na)-(skipna || skipmissing ? countnz(a.na) : 0))

## variance

function Base.varm(A::DataArray{T}, m::Number;
                   corrected::Bool=true, skipmissing::Bool=false, skipna::Bool=false) where T
    if skipna || skipmissing
        if skipna && !skipmissing
            Base.depwarn("skipna=$skipna is deprecated, use skipmissing=$skipna instead", :mapreduce)
        end

        n = length(A)
        na = A.na

        nna = countnz(na)
        nna == n && return convert(Base.momenttype(T), NaN)
        nna == n-1 && return convert(Base.momenttype(T),
                                     abs2(A.data[Base.findnextnot(na, 1)] - m)/(1 - corrected))

        /(nna == 0 ? Base.centralize_sumabs2(A.data, m, 1, n) :
                     mapreduce_impl_skipmissing(Base.centralizedabs2fun(m), +, A),
          n - nna - corrected)
    else
        any(A.na) && return missing
        Base.varm(A.data, m; corrected=corrected)
    end
end
Base.varm(A::DataArray{T}, m::Missing;
          corrected::Bool=true, skipmissing::Bool=false, skipna::Bool=false) where {T} = missing

function Base.var(A::DataArray;
                  corrected::Bool=true, mean=nothing, skipmissing::Bool=false, skipna::Bool=false)
    mean == 0 ? Base.varm(A, 0; corrected=corrected, skipmissing=skipmissing, skipna=skipna) :
    mean == nothing ? varm(A, Base.mean(A; skipmissing=skipmissing, skipna=skipna);
                           corrected=corrected, skipmissing=skipmissing, skipna=skipna) :
    isa(mean, Union{Number,Missing}) ?
        varm(A, mean; corrected=corrected, skipmissing=skipmissing, skipna=skipna) :
        throw(ErrorException("Invalid value of mean."))
end

Base.stdm(A::DataArray, m::Number;
          corrected::Bool=true, skipmissing::Bool=false, skipna::Bool=false) =
    sqrt(varm(A, m; corrected=corrected, skipmissing=skipmissing, skipna=skipna))

Base.std(A::DataArray;
         corrected::Bool=true, mean=nothing, skipmissing::Bool=false, skipna::Bool=false) =
    sqrt(var(A; corrected=corrected, mean=mean, skipmissing=skipmissing, skipna=skipna))

## weighted mean

function Base.mean(a::DataArray, w::Weights; skipmissing::Bool=false, skipna::Bool=false)
    if skipna || skipmissing
        v = a .* w.values
        sum(v; skipmissing=true) / sum(DataArray(w.values, v.na); skipmissing=true)
    else
        any(ismissing, a) ? missing : mean(a.data, w)
    end
end

function Base.mean(a::DataArray, w::Weights{W,V};
                   skipmissing::Bool=false, skipna::Bool=false) where {W,V<:DataArray}
    if skipna || skipmissing
        v = a .* w.values
        sum(v; skipmissing=true) / sum(DataArray(w.values.data, v.na); skipmissing=true)
    else
        any(ismissing, a) || any(ismissing, w.values) ? missing : wsum(a.data, w.values.data) / w.sum
    end
end


text:
using Oulu
using Test

@testset "Oulu.jl" begin
    # Write your tests here.
    @test Oulu.addNumbers(1, 2) == 3
    @test Oulu.addNumbers(1, 2) != 4
    @code_warntype Oulu.addNumbers(1, 2)
end


text:
function launch()

    @static if Sys.isapple()
        # OpenGL 3.2 + GLSL 150
        glsl_version = 150
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
        GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
        GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac
    else
        # OpenGL 3.0 + GLSL 130
        glsl_version = 130
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
        # GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
        # GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # 3.0+ only
    end

    # setup GLFW error callback
    error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
    GLFW.SetErrorCallback(error_callback)

    # create window
    window = GLFW.CreateWindow(1280, 720, "ElectroDermal Activity Analysis")
    @assert window != C_NULL
    GLFW.MakeContextCurrent(window)
    GLFW.SwapInterval(1)  # enable vsync

    # setup Dear ImGui context
    ctx = CImGui.CreateContext()

    # setup Dear ImGui style
    # CImGui.StyleColorsDark()
    # CImGui.StyleColorsClassic()
    CImGui.StyleColorsLight()

    # setup Platform/Renderer bindings
    ImGui_ImplGlfw_InitForOpenGL(window, true)
    ImGui_ImplOpenGL3_Init(glsl_version)
    should_show_dialog = true
    clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]

    open_file_dialog = OpenFileDialog(pwd(),"", pwd(),"", false, false)
    open_file_dialog2 = OpenFileDialog(pwd(),"",pwd(),"",false,false)
    open_files = false
    open_files2 = false
    data = DataFrame()
    dt = Float32[]
    while !GLFW.WindowShouldClose(window)

        GLFW.PollEvents()
        # start the Dear ImGui frame
        ImGui_ImplOpenGL3_NewFrame()
        ImGui_ImplGlfw_NewFrame()
        CImGui.NewFrame()

        if CImGui.BeginMainMenuBar()
            if CImGui.BeginMenu("File")
                populate_file_menu!(open_file_dialog) #visiblity = true
                CImGui.EndMenu()
            end
            CImGui.EndMainMenuBar()
        end

        if isvisible(open_file_dialog)
            display_dialog!(open_file_dialog)
            if has_pending_action(open_file_dialog)
                data = perform_dialog_action(open_file_dialog)
                consume_action!(open_file_dialog)
                dt = Cfloat.(data[3:end,1])
                Base.display(data)
                open_files = true
            end
        end
        if open_files
            begin
                CImGui.Begin("Plot")
                CImGui.Text("Plot")
                st = data[1,1]
                freq = data[2,1]
                CImGui.PlotLines("Result", dt, length(dt))

                CImGui.Text("Primitives")
                    sz, thickness, col = @cstatic sz=Cfloat(36.0) thickness=Cfloat(4.0) col=Cfloat[1.0,0.0,0.4,0.2] begin
                        @c CImGui.DragFloat("Size", &sz, 0.2, 2.0, 72.0, "%.0f")
                        @c CImGui.DragFloat("Thickness", &thickness, 0.05, 1.0, 8.0, "%.02f")
                        CImGui.ColorEdit4("Color", col)
                    end
                    p = CImGui.GetCursorScreenPos()
                    col32 = CImGui.ColorConvertFloat4ToU32(ImVec4(col...))
                    begin
                        width = 1000
                        height = 200
                        CImGui.PlotLines("EDA Measurements", dt, length(dt), 0 , "EDA", 0, 1.0, (width,height))
                        draw_list = CImGui.GetWindowDrawList()
                        #x::Cfloat = p.x + 4.0
                        #y::Cfloat = p.y + 4.0
                        x::Cfloat = p.x
                        y::Cfloat = p.y
                        spacing = 8.0
                        # Draws (almost transparent) horizontal bars
                        for yₙ in range(y, step = 40, stop = y + height - 40)
                            CImGui.AddRectFilled(draw_list, ImVec2(x, yₙ), ImVec2(x+width, yₙ+20), col32);
                        end
                    end
                    #draw the x axis
                    begin
                        p = CImGui.GetCursorScreenPos()
                        width = 1000
                        col = Cfloat[0.0,0.0,0.0,1.0]
                        col32 = CImGui.ColorConvertFloat4ToU32(ImVec4(col...))
                        draw_list = CImGui.GetWindowDrawList()
                        x = p.x
                        y = p.y
                        #CImGui.display(typeof(p))

                        time = Dates.unix2datetime(st)
                        # Base.display(time)
                        #min = minute(time)
                        #sec = second(time)
                        f = Cfloat(1 / freq)
                        CImGui.AddLine(draw_list, ImVec2(x, y), ImVec2(x+width, y), col32, Cfloat(1));
                        for xₙ in range(x, step = 40, stop = x + width)
                            CImGui.AddLine(draw_list, ImVec2(xₙ, y), ImVec2(xₙ, y-5), col32, Cfloat(1));
                            #p.x = xₙ;
                            min = minute(time);
                            sec = second(time);
                            mil = millisecond(time);
                            if mil ==0
                                CImGui.AddText(draw_list, ImVec2(xₙ, y), col32, string(string(min), ":", string(sec)));
                            end
                            next_time = time + Dates.Millisecond(250);
                            time = next_time;
                        end
                    end

                end
                CImGui.End()
            end

        # rendering
        CImGui.Render()
        GLFW.MakeContextCurrent(window)
        display_w, display_h = GLFW.GetFramebufferSize(window)
        glViewport(0, 0, display_w, display_h)
        glClearColor(clear_color...)
        glClear(GL_COLOR_BUFFER_BIT)
        ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())

        GLFW.MakeContextCurrent(window)
        GLFW.SwapBuffers(window)
    end

    # cleanup
    ImGui_ImplOpenGL3_Shutdown()
    ImGui_ImplGlfw_Shutdown()
    CImGui.DestroyContext(ctx)

    GLFW.DestroyWindow(window)
end

function perform_dialog_action(dialog::OpenFileDialog)
    directory = get_directory(dialog, ConfirmedStatus())
    file_name = get_file(dialog, ConfirmedStatus())
    path = joinpath(directory, file_name)
    data = CSV.File(path ; header= ["EDA"]) |> DataFrame
end

function populate_file_menu!(dialog::AbstractDialog)
    if CImGui.MenuItem("Open", "Ctrl+O")
        set_visibility!(dialog, true)
    end
    if CImGui.MenuItem("Save", "Ctrl+s")

    end
    if CImGui.MenuItem("Save as")

    end
    if CImGui.MenuItem("Quit", "Alt+F4")

    end
end

# function display_file_dialog(should_show_dialog::Bool, path::String)
#     if should_show_dialog
#         @c CImGui.Begin("Open Dialog", &should_show_dialog)
#         path_directories = splitpath(path)
#         selected_directory = Cint(length(path_directories))
#         # Draw a button for each directory that constitutes the current path.
#         for (index, d) in enumerate(path_directories)
#             CImGui.Button(d) && (selected_directory = Cint(index);)
#             CImGui.SameLine()
#         end
#         # If a button is clicked then we truncate the path up to and including the directory that was clicked.
#         path = selected_directory == 1 ? joinpath(first(path_directories)) : joinpath(path_directories[1:selected_directory]...)
#         CImGui.NewLine()
#         CImGui.BeginChild("Directory and File Listing", CImGui.ImVec2(CImGui.GetWindowWidth() * 0.98, -CImGui.GetWindowHeight() * 0.2))
#         CImGui.Columns(1)
#         # Make a list of directories that are visibile from the current directory.
#         visible_directories = filter(p->isdir(joinpath(path, p)), readdir(path))
#         for (n, folder_name) in enumerate(visible_directories)
#             # When the user clicks on a directory then change directory by appending it to the current path.
#             if CImGui.Selectable("[Dir] " * "$folder_name")
#                 @info "Trigger Item $n | find me here: $(@__FILE__) at line $(@__LINE__)"
#                 path = joinpath(path, folder_name)
#             end
#         end
#         # Make a list of files that are visible in the current directory.
#         visible_files = filter(p->isfile(joinpath(path, p)), readdir(path))
#         selected_file = Cint(0)
#         for (n, file_name) in enumerate(visible_files)
#             if CImGui.Selectable("[File] " * "$file_name")
#                 @info "Trigger Item $n | find me here: $(@__FILE__) at line $(@__LINE__)"
#             end
#         end
#         CImGui.EndChild()
#         CImGui.Text("File Name:")
#         CImGui.SameLine()
#         @show
#         #file_name₀ =  "\0"*"\0"^(255)
#         # Reserve 255 characters to receive input text. In the extreme
#         # case where a file is selected that exceeds 255 characters we will append
#         # the filename with a C_NULL character.
#         N = selected_file == 0 ? 255 : max(255 - length(visible_files[selected_file]), 1)
#         file_name₀ =  selected_file == 0 ?  "\0"*"\0"^(N)  : visible_files[selected_file]*"\0"^(N)
#         buf = Cstring(pointer(file_name₀))
#         CImGui.InputText("",buf, length(file_name₀))
#         #@show file_name₀
#         @show selected_file
#
#         #@cstatic  str0="Hello, world! " * "\0"^115 begin
#         # str0="Hello, world! "
#         # @show typeof(str0), str0
#         # ImGui.InputText("input text", str0, length(str0))
#
#
#         #CImGui.InputText("input text", str0, length(str0))
#         # Print the name of the selected file in the textbox.
#         CImGui.Button("Cancel") && (should_show_dialog = false;)
#         CImGui.SameLine()
#         CImGui.Button("Open") && (should_show_dialog = false;)
#
#         #@show current_directory
#         #CImGui.Button("Close Me") && (should_show_dialog = false;)
#         CImGui.End()
#     end
#
#     should_show_dialog, path
# end


text:
mutable struct Course
    name::AbstractString            # Name of the course
    id::Int                         # Unique id for course
    term::Int                       # The term the course belongs to

    credits::Int                    # Course credit hours
    delay::Int                      # The course's delay factor - the value of the longest path the course is on
    blocking::Int                   # The course's blocking factor - the number of courses blocked by this course
    cruciality::Int                 # The sum of the delay and blocking factor
    centrality::Int                 # Course centrality
    reachability::Int               # Course reachability

    prereqs::Array{Course}          # Array of courses that are a prerequisite to this course
    coreqs::Array{Course}           # Array of courses that must be taken in the same term or before
    postreqs::Array{Course}         # Array of courses that has this course as a prerequisite
    termReq::Int                    # Number of terms that must be completed before enrolling

    students::Array{Student}        # Array of students enrolled in this course in a given term

    model::Dict                     # Dictionary to hold information for grade predictions
    passrate::Float64               # Percentage of students that pass the course
    failures::Int                   # Total number of students who do not pass the course
    grades::Array{Float64}          # Array of all grades made in the course
    enrolled::Int                   # Total number of students that are enrolled in the course
    termenrollment::Array{Int}      # An array of enrollment by term
    termpassed::Array{Int}          # An array of the number of students who pass each term

    # Minimum Information Required
    function Course(name::AbstractString, credits::Number, prereqs::Array{Course}, coreqs::Array{Course})
        this = new()

        this.name = name
        this.credits = credits
        this.prereqs = prereqs
        this.postreqs = Course[]
        this.coreqs = Course[]
        this.termReq = 0

        # Assign current course as postreq for prereqs
        for c in this.prereqs
            push!(c.postreqs, this)
        end

        this.passrate = 0
        this.failures = 0
        this.enrolled = 0
        this.grades = Float64[]
        this.students = Student[]

        return this
    end

    # With Passrate
    function Course(name::AbstractString, credits::Number, passrate::Number, prereqs::Array{Course}, coreqs::Array{Course})
        this = Course(name, credits, prereqs, coreqs)
        this.passrate = passrate
        return this
    end

    # With termReq
    function Course(name::AbstractString, credits::Number, prereqs::Array{Course}, coreqs::Array{Course}, termReq::Number)
        this = Course(name, credits, prereqs, coreqs)
        this.termReq = termReq
        return this
    end

    # With termReq and Passrate
    function Course(name::AbstractString, credits::Number, termReq::Number, passrate::Number, prereqs::Array{Course}, coreqs::Array{Course})
        this = Course(name, credits, prereqs, coreqs)
        this.termReq = termReq
        this.passrate = passrate
        return this
    end
end

# Helper Functions

# Compute the cruciality, delay, and blocking factor of a course
function cruciality(course::Course)
    b = blocking(course)
    d = delay(course)

    course.blocking = b
    course.delay = d
    course.cruciality = d + b
end

# Compute the blocking factor of a course
function blocking(course, b=-1, courses=[])
    for c in course.postreqs
        b += blocking(c, 0, courses)
    end

    if !in(course.name, courses)
        b += 1
        push!(courses, course.name)
    end

    return b
end

# Compute the delay factor of a course
function delay(course)
    # Traverse forward from course
    f = []
    forward(course, f)

    # Traverse backwards from course
    b = []
    backward(course, b)

    return maximum(f) + maximum(b) + 1
end

# Determines all path length (l) from course and stores them in path
function forward(course, paths, l=0)
    # Keep traversing forward
    for c in course.postreqs
        forward(c, paths, l+1)
    end

    # If no postreqs, record path length
    if length(course.postreqs) == 0
        push!(paths, l)
    end
end

# Determines all path length (l) behind course and stores them in path
function backward(course, paths, l=0)
    # Keep traversing backward
    for c in [course.prereqs; course.coreqs]
        backward(c, paths, l+1)
    end

    # If no prereqs, record path length
    if length(course.prereqs) == 0
        push!(paths, l)
    end
end

max_stars_repo_name: M-S-H/CASL

max_issues_repo_name: M-S-H/CASL

max_forks_repo_name: M-S-H/CASL


text:
using vOptGeneric, JuMP
using GLPK

M = 50
p1 =   [M 4 5 M M M;
        M M 2 1 2 7;
        M M M 5 2 M;
        M M 5 M M 3;
        M M M M M 4;
        M M M M M M]
p2 =   [M 3 1 M M M;
        M M 1 4 2 2;
        M M M 1 7 M;
        M M 1 M M 2;
        M M M M M 2;
        M M M M M M]

m = vModel(GLPK.Optimizer) ; JuMP.set_silent(m)

@variable(m, x[1:6,1:6], Bin)

@addobjective(m, Min , sum(x[i,j]*p1[i,j] for i=1:6, j=1:6))
@addobjective(m, Min , sum(x[i,j]*p2[i,j] for i=1:6, j=1:6))
@constraint(m, eps, sum(x[1,j] for j = 1:6) == 1)
@constraint(m, sum(x[i,6] for i = 1:6) == 1)
@constraint(m, cstr[i=2:5], sum(x[i,j] for j = 1:6) - sum(x[j,i] for j =1:6) == 0)

vSolve(m, method=:epsilon, step=0.5, verbose=true)

#get the results and print them
Y_N = getY_N(m)
for n = 1:length(Y_N)
    X = value.(x, n)
    for ind in findall(val -> val ≈ 1, X)
        i,j = ind.I
        print("$i->$j ")
    end
    println("| z = ",Y_N[n])
end