text:
#
# impl.jl -
#
# Implementation of fast discrete correlation methods.
#
module Impl

import Base.OneTo

abstract type Operation end
struct Direct  <: Operation end # apply direct transform
struct Adjoint <: Operation end # apply adjoint transform

# Options.
#
# - Option `JUMP` specifies that, if the dimension of operation is the
#   first one (that is `d = 1`), then the innermost loop is performed along
#   the second dimension of the arrays.
#
const JUMP        =  1 # alow non-unit jump
const LOCALVALUES =  2 # store values (instead of indices) in local variables
const SPLIT       =  4 # split cases
const SIMD        =  8 # vectorization
const REFERENCE   = 16 # reference implementation

#const Indices{N} = NTuple{N,AbstractUnitRange{Int}}
const Dims{N} = NTuple{N,Int}

# If there is no version managing a specific combination of options and
# dimensions, the following will dispatch to methods using CartesianIndices
# for the indices before and after the dimension of operation.
function correlate!(dir::Operation,
                    opt::Int,
                    dst::AbstractArray{T,N},
                    src::AbstractArray{T,N},
                    off::Int,
                    wgt::NTuple{S,T},
                    d::Int,
                    dims::Dims{N},
                    m::Int,
                    n::Int) where {T<:AbstractFloat,N,S}
    correlate!(dir, opt, dst, src, off, wgt, CartesianIndices(dims[1:d-1]),
               CartesianIndices(dims[d+1:N]), m, n)
end

# Direct 1D correlation.
@generated function correlate!(::Direct,
                               opt::Int,
                               dst::AbstractArray{T,1},
                               src::AbstractArray{T,1},
                               off::Int,
                               wgt::NTuple{S,T},
                               d::Int,
                               dims::Dims{1},
                               m::Int,
                               n::Int) where {T<:AbstractFloat,S}
    @assert S ≥ 2
    W = ntuple(k -> Symbol(:w,k), Val(S)) # weights
    X = ntuple(k -> Symbol(:x,k), Val(S)) # neighbor values
    J = ntuple(k -> Symbol(:j,k), Val(S)) # neighbor indices
    Jlast = :(clamp((S-off)+i,1,n)) # index of last neighbor
    init_inds = ntuple(k -> :($(J[k+1]) = clamp($(k+1)-off,1,n)), Val(S-1))
    update_inds = assign_expr(J, (J[2:S]..., Jlast))
    Xload = ntuple(k -> :(src[$(J[k])]), Val(S))
    quote
        @assert n ≥ 1
        $(tuple_expr(W)) = wgt
        $(init_inds...)
        @inbounds @simd for i in 1:m
            $(update_inds)
            dst[i] = $(dot_expr(Xload, W))
        end
    end
end

# Adjoint 1D correlation.
@generated function correlate!(::Adjoint,
                               opt::Int,
                               dst::AbstractArray{T,1},
                               src::AbstractArray{T,1},
                               off::Int,
                               wgt::NTuple{S,T},
                               d::Int,
                               dims::Dims{1},
                               m::Int,
                               n::Int) where {T<:AbstractFloat,S}
    @assert S ≥ 2
    W = ntuple(k -> Symbol(:w,k), Val(S)) # weights
    J = ntuple(k -> Symbol(:j,k), Val(S)) # neighbor indices
    Jlast = :(clamp((S-off)+i,1,n)) # index of last neighbor
    init_inds = ntuple(k -> :($(J[k+1]) = clamp($(k+1)-off,1,n)), Val(S-1))
    update_inds = assign_expr(J, (J[2:S]..., Jlast))
    incr_dest = (:(tmp = src[i]),
                 ntuple(k -> :(dst[$(J[k])] += tmp*$(W[k])), Val(S))...)
    quote
        @assert n ≥ 1
        $(tuple_expr(W)) = wgt
        $(init_inds...)
        @inbounds for i in 1:m
            $(update_inds)
            $(incr_dest...)
        end
    end
end

# Direct 2D correlation.
@generated function correlate!(::Direct,
                               opt::Int,
                               dst::AbstractArray{T,2},
                               src::AbstractArray{T,2},
                               off::Int,
                               wgt::NTuple{S,T},
                               d::Int,
                               dims::Dims{2},
                               m::Int,
                               n::Int) where {T<:AbstractFloat,S}
    @assert S ≥ 2
    W = ntuple(k -> Symbol(:w,k), Val(S)) # weights
    X = ntuple(k -> Symbol(:x,k), Val(S)) # neighbor values
    J = ntuple(k -> Symbol(:j,k), Val(S)) # neighbor indices
    P = ntuple(k -> Symbol(:p,k), Val(S)) # pre-initialized indices
    Jlast = :(clamp((S-off)+i,1,n)) # index of last neighbor
    init_inds = assign_expr(J[2:S], P[2:S])
    update_inds = assign_expr(J, (J[2:S]..., Jlast))
    Xload1 = ntuple(k -> :(src[$(J[k]),t]), Val(S))
    Xload2 = ntuple(k -> :(src[t,$(J[k])]), Val(S))
    quote
        @assert n ≥ 1
        @inbounds o = dims[3-d]
        $(tuple_expr(W)) = wgt
        $(ntuple(k -> :($(P[k+1]) = clamp($(k+1)-off,1,n)), Val(S-1))...)
        if d > 1
            $(init_inds)
            @inbounds for i in 1:m
                $(update_inds)
                for t in 1:o # @simd does not help here
                    dst[t,i] = $(dot_expr(Xload2, W))
                end
            end
        elseif (opt & JUMP) != 0
            $(init_inds)
            @inbounds for i in 1:m
                $(update_inds)
                for t in 1:o # @simd does not help here
                    dst[i,t] = $(dot_expr(Xload1, W))
                end
            end
        else
            @inbounds for t in 1:o
                $(init_inds)
                @simd for i in 1:m # @simd does not help here
                    $(update_inds)
                    dst[i,t] = $(dot_expr(Xload1, W))
                end
            end
        end
    end
end

# Adjoint 2D correlation.
@generated function correlate!(::Adjoint,
                               opt::Int,
                               dst::AbstractArray{T,2},
                               src::AbstractArray{T,2},
                               off::Int,
                               wgt::NTuple{S,T},
                               d::Int,
                               dims::Dims{2},
                               m::Int,
                               n::Int) where {T<:AbstractFloat,S}
    @assert S ≥ 2
    W = ntuple(k -> Symbol(:w,k), Val(S)) # weights
    P = ntuple(k -> Symbol(:p,k), Val(S)) # initial indices
    J = ntuple(k -> Symbol(:j,k), Val(S)) # neighbor indices
    Jlast = :(clamp((S-off)+i,1,n)) # index of last neighbor
    init_inds = assign_expr(J[2:S], P[2:S])
    update_inds = assign_expr(J, (J[2:S]..., Jlast))
    incr_dest1 = (:(tmp = src[i,t]),
                  ntuple(k -> :(dst[$(J[k]),t] += tmp*$(W[k])), Val(S))...)
    incr_dest2 = (:(tmp = src[t,i]),
                  ntuple(k -> :(dst[t,$(J[k])] += tmp*$(W[k])), Val(S))...)
    quote
        @assert n ≥ 1
        @inbounds o = dims[3-d]
        $(tuple_expr(W)) = wgt
        $(ntuple(k -> :($(P[k+1]) = clamp($(k+1)-off,1,n)), Val(S-1))...)
        if d > 1
            $(init_inds)
            @inbounds for i in 1:m
                $(update_inds)
                for t in 1:o # @simd does not help here
                    $(incr_dest2...)
                end
            end
        elseif (opt & JUMP) != 0
            $(init_inds)
            @inbounds for i in 1:m
                $(update_inds)
                for t in 1:o # @simd does not help here
                    $(incr_dest1...)
                end
            end
        else
            @inbounds for t in 1:o
                $(init_inds)
                @simd for i in 1:m # @simd does not help here
                    $(update_inds)
                    $(incr_dest1...)
                end
            end
        end
    end
end

# Default version for other number of dimensions.
@generated function correlate!(::Direct,
                               opt::Int,
                               dst::AbstractArray{T,N},
                               src::AbstractArray{T,N},
                               off::Int,
                               wgt::NTuple{S,T},
                               Ipre::CartesianIndices{Npre},
                               Ipost::CartesianIndices,
                               m::Int,
                               n::Int) where {T,S,N,Npre}
    # Symbols for local variables.
    W = ntuple(k -> Symbol(:w,k), Val(S)) # weights
    J = ntuple(k -> Symbol(:j,k), Val(S)) # neighbor indices
    P = ntuple(k -> Symbol(:p,k), Val(S)) # pre-initialized indices

    # Extract weights and pre-initialize indices..
    preamble = (assign_expr(W, :wgt),
                ntuple(k -> :($(P[k+1]) = clamp($(k+1)-off,1,n)), S-1)...)

    # Index of last neighbor.
    Jlast = :(clamp((S-off)+i,1,n))

    # Initialize indices before the loop along the active dimension.
    init_inds = assign_expr(J[2:S], P[2:S])

    # Update indices into the loop along the active dimension.
    update_inds = assign_expr(J, (J[2:S]..., Jlast))

    if Npre > 0
        # Operate along another dimension than the 1st one.
        Xload = ntuple(k -> :(src[ipre,$(J[k]),ipost]), Val(S))
        return quote
            $(preamble...)
            @inbounds for ipost in Ipost
                $(init_inds)
                for i in 1:m
                    $(update_inds)
                    for ipre in Ipre
                        dst[ipre,i,ipost] = $(dot_expr(Xload, W))
                    end
                end
            end
        end
    else
        # Operate along the 1st dimension.
        X = ntuple(k -> Symbol(:x,k), Val(S)) # neighbor values
        Xload = ntuple(k -> :(src[$(J[k]),ipost]), Val(S))
        init_vals = assign_expr(X[2:S], (k -> :(src[$(P[k+1]),ipost]), S-1))
        update_vals = assign_expr(X, (X[2:S]..., :(src[$(Jlast),ipost])))
        return quote
            $(preamble...)
            if (opt & LOCALVALUES) != 0
                @inbounds for ipost in Ipost
                    $(init_vals)
                    for i in 1:m
                        $(update_vals)
                        dst[i,ipost] = $(dot_expr(X, W))
                    end
                end
            else
                @inbounds for ipost in Ipost
                    $(init_inds)
                    for i in 1:m
                        $(update_inds)
                        dst[i,ipost] = $(dot_expr(Xload, W))
                    end
                end
            end
        end
    end
end

@generated function correlate!(::Adjoint,
                               opt::Int,
                               dst::AbstractArray{T,N},
                               src::AbstractArray{T,N},
                               off::Int,
                               wgt::NTuple{S,T},
                               Ipre::CartesianIndices{Npre},
                               Ipost::CartesianIndices,
                               m::Int,
                               n::Int) where {T,S,N,Npre}
    # Symbols for local variables.
    W = ntuple(k -> Symbol(:w,k), Val(S)) # weights
    J = ntuple(k -> Symbol(:j,k), Val(S)) # neighbor indices
    P = ntuple(k -> Symbol(:p,k), Val(S)) # pre-initialized indices

    # Extract weights and pre-initialize indices..
    preamble = (assign_expr(W, :wgt),
                ntuple(k -> :($(P[k+1]) = clamp($(k+1)-off,1,n)), S-1)...)

    # Initialize indices before the loop along the active dimension.
    init_inds = assign_expr(J[2:S], P[2:S])

    # Update indices into the loop along the active dimension.
    update_inds = assign_expr(J, (J[2:S]..., :(clamp((S-off)+i,1,n))))

    if Npre > 0
        # Operate along another dimension than the 1st one.
        return quote
            $(preamble...)
            @inbounds for ipost in Ipost
                $(init_inds)
                for i in 1:m
                    $(update_inds)
                    for ipre in Ipre
                        tmp = src[ipre,i,ipost]
                        $(ntuple(k -> :(dst[ipre,$(J[k]),ipost]
                                        += tmp*$(W[k])), S)...)
                    end
                end
            end
        end
    else
        # Operate along the 1st dimension.
        return quote
            $(preamble...)
            @inbounds for ipost in Ipost
                $(init_inds)
                for i in 1:m
                    $(update_inds)
                    tmp = src[i,ipost]
                    $(ntuple(k -> :(dst[$(J[k]),ipost] += tmp*$(W[k])), S)...)
                end
            end
        end
    end
end

#------------------------------------------------------------------------------
# METAPROGRAMMING UTILITIES

# Types for expression operands.
const Operand = Union{Expr,Symbol}

# Tuple of operands.
const Operands = Tuple{Vararg{Operand}}

# Paris that can be converted to a tuple of operands.
const Listable = Union{Tuple{Function,Integer},
                       Tuple{Function,Val}}

"""
    assign_expr(lval, rval)

yields an expression to assign `rval` to `lval`.

"""
assign_expr(lval::Operand, rval::Operand) = Expr(:(=), lval, rval)
assign_expr(lval::Operands, rval) = assign_expr(tuple_expr(lval), rval)
assign_expr(lval::Listable, rval) = assign_expr(tuple_expr(lval...), rval)
assign_expr(lval::Operand, rval::Operands) =
    assign_expr(lval, tuple_expr(rval))
assign_expr(lval::Operand, rval::Listable) =
    assign_expr(lval, tuple_expr(rval...))

"""
    sum_expr(args)

yields an expression which is the sum of the expressions in `args`.

"""
sum_expr(args::Tuple{Vararg{Operand}}) = Expr(:call, :(+), args...)
sum_expr(args::Operand...) = sum_expr(args)

sum_expr(fn::Function, n::Integer) = sum_expr(fn, Val(Int(n)))
sum_expr(fn::Function, val::Val) = sum_expr(ntuple(fn, val))

"""
    dot_expr(X,Y)

yields an expression which is the sum of the elementwise product of
`X` and `Y`, 2 tuples of expressions and/or symbols.

"""
dot_expr(X::NTuple{N,Operand}, Y::NTuple{N,Operand}) where {N} =
    sum_expr(k -> :($(X[k])*$(Y[k])), Val(N))

# Make X into an operand.
dot_expr(X::Tuple{Function,Integer}, Y) = dot_expr((X[1], Val(Int(X[2]))), Y)
dot_expr(X::Tuple{Function,Val}, Y) = dot_expr(ntuple(X...), Y)

# Make Y into an operand.
dot_expr(X::NTuple{N,Operand}, Y::Tuple{Function,Integer}) where {N} =
    dot_expr(X, (Y[1], Val(Int(Y[2]))))
dot_expr(X::NTuple{N,Operand}, Y::Tuple{Function,Val}) where {N} =
    dot_expr(X, ntuple(Y...))

"""
    tuple_expr(args)

yields an expression which is a tuple of the expressions in `args`.

"""
tuple_expr(args::Tuple{Vararg{Operand}}) = Expr(:tuple, args...)
tuple_expr(args::Operand...) = tuple_expr(args)

tuple_expr(fn::Function, n::Integer) = tuple_expr(fn, Val(Int(n)))
tuple_expr(fn::Function, val::Val) = tuple_expr(ntuple(fn, val))

end # module


text:
"""
$(TYPEDEF)
"""
mutable struct Point3Point3{D <: IIF.SamplableBelief} <: AbstractRelativeRoots
    Zij::D
    # empty constructor
    Point3Point3{T}() where T = new{T}()
    # regular constructor
    Point3Point3{T}(x::T) where {T <: IIF.SamplableBelief} = new{T}(x)
end
# convenience and default object helper
Point3Point3(x::T=MvNormal(zeros(3),LinearAlgebra.diagm([0.1;0.1;0.1]))) where {T <: IIF.SamplableBelief} = Point3Point3{T}(x)

function getSample(cfo::CalcFactor{<:Point3Point3})
  return rand(cfo.factor.Zij)
end
function (cf::CalcFactor{<:Point3Point3})(meas,
                                            xi,
                                            xj )
  #
  return meas[1:3] .- (xj[1:3] .- xi[1:3])
end

"""
$(TYPEDEF)

Serialization type for `Point3Point3`.
"""
mutable struct PackedPoint3Point3 <: IncrementalInference.PackedInferenceType
    str::String
    PackedPoint3Point3() = new()
    PackedPoint3Point3(s::AS) where {AS <: AbstractString} = new(s)
end
function convert(::Type{Point3Point3}, d::PackedPoint3Point3)
  return Point3Point3( convert(SamplableBelief, d.str) )
end
function convert(::Type{PackedPoint3Point3}, d::Point3Point3)
  return PackedPoint3Point3( convert(PackedSamplableBelief, d.Zij) )
end

ext: jl


text:
#Defines the control system type, including its discrete time behaviour
#include("DiscreteTimeSystem.jl")
#= ControlFunction should input the "measured state" of the vehicle obtained from the
    Navigation system, the attitude/rate commands from the guidance system, and time.
=#

#Define ControlSystem as a subtype of the DiscreteTimeSystem
abstract type AbstractControlSystem <: DiscreteTimeSystem end

########### Template function ###########################
mutable struct ControlSystem{T <: Number} <: AbstractControlSystem
    Active::Bool #Activation Boolean
    Ts::T #Sample Time
    ControlRelation::Function #Encodes
    Tprev::T # Storage for previous activation values
end
ControlSystem() = ControlSystem(false, 0.0, x -> nothing, 0.0)
ControlSystem(Ts::Number) = ControlSystem(false, Float64(Ts), x -> nothing, 0.0)
ControlSystem(Ts::Number, csfunc::Function) = ControlSystem(true, Float64(Ts), csfunc, 0.0)
ControlSystem(Active::Bool, Ts::Number, csfunc::Function) = ControlSystem(Active, Float64(Ts), csfunc, 0.0)

#Set the callback affect condition of the ControlSystem
#placeholder!
function affect!(int,dts::Type{<:AbstractControlSystem})
    for cu in user_cache(int)
        cu.u = [dts.ControlRelation(int.u)...]
    end
  dts.Tprev = int.t
end

#=cb = ContinuousCallback((t, u, int) -> condition(det,t,u,int),
                          (int) -> affect!(int,det))
=#


text:
#Grid_vs_Bayes_XOR.jl

#Here we will optimize the learning rate and sigmoid parameter for the XOR neural network

include("XOR_MD.jl")
include("Kernals.jl")
include("gaussian_process.jl")

#Set the random seed:
srand(1234) #Seed for stalzer srand(1234)

"""
We will compare the effect of randomly selecting a learning rate and sigmoid
hyperparemeter vs the use of Bayesian Optimization for finding the optimial LR on the MSE.
Suppose we have limited computing time of 100000 epochs and that we haveju N tries to
minimise the MSE. Let us say that the learning rate is between a and b

"""

#Initialise Layers and params ==========================================

Layer_1=uniform(0,1,2,2)

Layer_2=uniform(0,1,2,1)

epoch_vec=linspace(100,1000,10)   # For the final report up this to 100 and leave for 10 minute to get smoothest graph
epoch_bayes_result=zeros(epoch_vec)
epoch_random_result=zeros(epoch_vec)

a=0.001  #Change to 0.001
b=1

c=0.001 #Was the same as above
d=1 #Change back to what it was before which was same as a, b

N=20

#Curry the sigmoid functions:

function hyper_curry(h)
    return (x->sigmoid(x,h))
end

function hyper_curry_deriv(h)
    return (x->sigmoid_deriv(x,h))
end

Random_Learning_Rates=uniform(a,b,N,1)
Random_Hyperparameters=uniform(c,d,N,1)
Random_Mat=cat(2,Random_Learning_Rates,Random_Hyperparameters)
Random_MSE=zeros(N)

timingsrand=zeros(length(epoch_vec))
timingsbays=zeros(length(epoch_vec))
convergence=zeros(length(epoch_vec))

for p=1:length(epoch_vec)

    epochs=epoch_vec[p]

    # Random Learning Rates Examples ========================================

    #Random_Mat conjoins Random_Learning_Rates and Random_Hyperparameters
    # Random_Mat is a Nx2 matrix where Random_Mat[1,:] is the first entry
    #with LR_1 and hyperparemeter 1.

    tic()
    for i=1:length(Random_Learning_Rates)
        node_function=hyper_curry(Random_Mat[i,2])
        node_deriv=hyper_curry_deriv(Random_Mat[i,2])
        learning_rate=Random_Mat[i,1]
        Random_MSE[i]=Train_Neural_Net_Loop(epochs,Layer_1,Layer_2,learning_rate,node_function,node_deriv)[3]
        println("Epoch Complete")
    end

    println("Random Learning Rates Training Completed")
    q=toc()
    timingsrand[p]=q

    #Bayesian Optimization Examples===================================================

    #Here are the points we can pick from in the Optimization

    LR_Test=linspace(a,b,40)
    HP_Test=linspace(c,d,40)

    #Here is the carteisan product of these written as a vector
    Test=gen_points([LR_Test,HP_Test])[1]

    #We first have to pick a random point to begin bayesian optimization:

    #currently starts with the midpoint, possibly randomise this:
    Bayesian_Points=[Test[Int(round(length(Test)/2))]]

    #Bayesian_Points is an vector of arrays where in each array first entry is LR second entry is Hyper-Parameters:

    #Define hyperparemeter functions:
    node_function=hyper_curry(Bayesian_Points[1][2])
    node_deriv=hyper_curry_deriv(Bayesian_Points[1][2])

    #Define Learning Rate:
    learning_rate=Bayesian_Points[1][1]

    #Run first train before Bayesian Optimization:
    Bayesian_MSE=[Train_Neural_Net_Loop(epochs,Layer_1,Layer_2,learning_rate,node_function,node_deriv)[3]]

    #Begin Bayesian Optimization:
    tic()
    for k=2:N
        D=[(Bayesian_Points[i],Bayesian_MSE[i]) for i=1:length(Bayesian_Points)]
        mu, sigma, D=gaussian_process_chol(std_exp_square_ker,D,1e-6,Test)
        println("Gaussian Process Complete","\r")
        mu=reshape(mu,length(mu));
        sigma=reshape(sigma,length(sigma))

        new_point=findmin(mu-sigma)[2]

        #Here we will need to change the number 2 to k
        Bayesian_Points=cat(1,Bayesian_Points,[Test[new_point]])

        learning_rate=Bayesian_Points[k][1]

        node_function=hyper_curry(Bayesian_Points[k][2])

        node_deriv=hyper_curry_deriv(Bayesian_Points[k][2])

        value_to_be_appended=Train_Neural_Net_Loop(epochs,Layer_1,Layer_2,learning_rate,node_function,node_deriv)[3]

        if value_to_be_appended !=Bayesian_MSE[k-1]
            Bayesian_MSE=cat(1,Bayesian_MSE,[value_to_be_appended])
            println("Epoch Complete")
            if k==N
                convergence[p]=N
            end

        else
            println("Found Optimum on the ", k-1, " iteration of ", N, " iterations")
            Bayesian_Points=Bayesian_Points[1:length(Bayesian_Points)-1]
            convergence[p]=k-1

            break
        end

    end
    q=toc()
    timingsbays[p]=q

    # Bayesian Plotting =========================================================

    println(" The optimium is located at ",Bayesian_Points[end])

    epoch_random_result[p]=minimum(Random_MSE)
    epoch_bayes_result[p]=minimum(Bayesian_MSE)

    #Move this to the bottom
    println("Bayesian_Learning_Rates Training Complete")
    println("The minimum MSE by Bayesian Optimization was", minimum(Bayesian_MSE))
    println("The mininmum MSE by Random Selection was", minimum(Random_MSE))

    println("completed cycle ",p, " out of overall cycle", length(epoch_vec))
end

using PyPlot
# fig = figure("pyplot_subplot_mixed",figsize=(7,7))
# ax=axes()
plot(epoch_vec,epoch_bayes_result,label="Bayesian Optimization")

plot(epoch_vec,epoch_random_result,label="Random Grid Search",alpha=0.7)
title("MSE Plot for different epochs")
xlabel("Epochs")
ylabel("MSE")
legend(loc="upper right",fancybox="true")
grid("on")
show()

using PyPlot

plot(epoch_vec,timingsrand,label="Random Timings")
plot(epoch_vec,timingsbays, label="Bayesian Op Timings")
title("Timings")
xlabel("Epochs")
ylabel("Time (s)")
legend(loc="upper right",fancybox="true")
grid("on")
show()

using PyPlot
plot(epoch_vec,convergence)
title("convergance")
xlabel("Epochs")
ylabel("convergance val")
grid("on")
show()


text:
function calculate_transcriptional_control_array(t::Float64,x::Array{Float64,1}, problem::Dict{String,Any})::Float64

    # initialize -
    u_variable = 0.0

    # alias elements of the species vector -
    mRNA = x[1]
    G = x[2]
    σ70 = x[3]

    # get stuff from the problem dictionary -
    E1 = problem["E1"]
    E2 = problem["E2"]
    R = problem["ideal_gas_constant_R"]
    T_K = problem["temperature_K"]
    KD = problem["inducer_dissociation_constant"]
    n = problem["inducer_cooperativity_parameter"]

    # TODO: write u-varible function here
    f_I = (σ70^n)/(KD^n + σ70^n)
    W1 = exp(-E1/(R*T_K))
    W2 = exp(-E2/(R*T_K))
    Z = 1 + W1 + f_I*W2
    u_variable = (W1 + f_I*W2)/Z

    # return -
    return u_variable
end


text:
"""
Appends a tracking process to determine the time of the callback to be used in
the reverse pass. The rationale is explain in:

https://github.com/SciML/DiffEqSensitivity.jl/issues/4
"""
track_callbacks(cb,t,u,p,sensealg) = track_callbacks(CallbackSet(cb),t,u,p,sensealg)
track_callbacks(cb::CallbackSet,t,u,p,sensealg) = CallbackSet(
   map(cb->_track_callback(cb,t,u,p,sensealg), cb.continuous_callbacks),
   map(cb->_track_callback(cb,t,u,p,sensealg), cb.discrete_callbacks))

mutable struct ImplicitCorrection{T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,RefType}
  gt_val::T1
  gu_val::T2
  gt::T3
  gu::T4
  gt_conf::T5
  gu_conf::T6
  condition::T7
  Lu_left::T8
  Lu_right::T9
  dy_left::T10
  dy_right::T11
  cur_time::RefType # initialized as "dummy" Ref that gets overwritten by Ref of loss
  terminated::Bool
end

ImplicitCorrection(cb::DiscreteCallback,t,u,p,sensealg) = nothing
function ImplicitCorrection(cb,t,u,p,sensealg)
  condition = cb.condition

  gt_val = similar(u,1)
  gu_val = similar(u)

  fakeinteg = FakeIntegrator(u,p,t,t)
  gt, gu = build_condition_wrappers(cb,condition,u,t,fakeinteg)

  gt_conf = build_deriv_config(sensealg,gt,gt_val,t)
  gu_conf = build_grad_config(sensealg,gu,u,p)

  dy_left = similar(u)
  dy_right = similar(u)

  Lu_left = similar(u)
  Lu_right = similar(u)

  cur_time = Ref(1) # initialize the Ref, set to Ref of loss below
  terminated = false

  ImplicitCorrection(gt_val,gu_val,gt,gu,gt_conf,gu_conf,condition,Lu_left,Lu_right,dy_left,dy_right,cur_time,terminated)
end

struct TrackedAffect{T,T2,T3,T4,T5,T6}
    event_times::Vector{T}
    tprev::Vector{T}
    uleft::Vector{T2}
    pleft::Vector{T3}
    affect!::T4
    correction::T5
    event_idx::Vector{T6}
end

TrackedAffect(t::Number,u,p,affect!::Nothing,correction) = nothing
TrackedAffect(t::Number,u,p,affect!,correction) = TrackedAffect(Vector{typeof(t)}(undef,0),Vector{typeof(t)}(undef,0),
                                                                Vector{typeof(u)}(undef,0),Vector{typeof(p)}(undef,0),affect!,correction,
                                                                Vector{Int}(undef,0))

function (f::TrackedAffect)(integrator,event_idx=nothing)
    uleft = deepcopy(integrator.u)
    pleft = deepcopy(integrator.p)
    if event_idx===nothing
      f.affect!(integrator)
    else
      f.affect!(integrator,event_idx)
    end
    if integrator.u_modified
        if isempty(f.event_times)
          push!(f.event_times,integrator.t)
          push!(f.tprev,integrator.tprev)
          push!(f.uleft,uleft)
          push!(f.pleft,pleft)
          if event_idx !== nothing
            push!(f.event_idx,event_idx)
          end
        else
          if !maximum(.≈(integrator.t, f.event_times, rtol=0.0, atol=1e-14))
            push!(f.event_times,integrator.t)
            push!(f.tprev,integrator.tprev)
            push!(f.uleft,uleft)
            push!(f.pleft,pleft)
            if event_idx !== nothing
              push!(f.event_idx, event_idx)
            end
          end
        end
    end
end

function _track_callback(cb::DiscreteCallback,t,u,p,sensealg)
    correction = ImplicitCorrection(cb,t,u,p,sensealg)
    DiscreteCallback(cb.condition,
                     TrackedAffect(t,u,p,cb.affect!,correction),
                     cb.initialize,
                     cb.finalize,
                     cb.save_positions)
end

function _track_callback(cb::ContinuousCallback,t,u,p,sensealg)
    correction = ImplicitCorrection(cb,t,u,p,sensealg)
    ContinuousCallback(
        cb.condition,
        TrackedAffect(t,u,p,cb.affect!,correction),
        TrackedAffect(t,u,p,cb.affect_neg!,correction),
        cb.initialize,
        cb.finalize,
        cb.idxs,
        cb.rootfind,cb.interp_points,
        cb.save_positions,
        cb.dtrelax,cb.abstol,cb.reltol,cb.repeat_nudge)
end

function _track_callback(cb::VectorContinuousCallback,t,u,p,sensealg)
    correction = ImplicitCorrection(cb,t,u,p,sensealg)
    VectorContinuousCallback(
               cb.condition,
               TrackedAffect(t,u,p,cb.affect!,correction),
               TrackedAffect(t,u,p,cb.affect_neg!,correction),
               cb.len,cb.initialize,cb.finalize,cb.idxs,
               cb.rootfind,cb.interp_points,
               collect(cb.save_positions),
               cb.dtrelax,cb.abstol,cb.reltol,cb.repeat_nudge)
end

struct FakeIntegrator{uType,P,tType,tprevType}
    u::uType
    p::P
    t::tType
    tprev::tprevType
end

struct CallbackSensitivityFunction{fType,Alg<:DiffEqBase.AbstractSensitivityAlgorithm,C<:AdjointDiffCache,pType} <: SensitivityFunction
    f::fType
    sensealg::Alg
    diffcache::C
    prob::pType
end

getprob(S::CallbackSensitivityFunction) = S.prob
inplace_sensitivity(S::CallbackSensitivityFunction) = true

"""
Sets up callbacks for the adjoint pass. This is a version that has an effect
at each event of the forward pass and defines the reverse pass values via the
vjps as described in https://arxiv.org/pdf/1905.10403.pdf Equation 13.

For more information, see https://github.com/SciML/DiffEqSensitivity.jl/issues/4
"""
setup_reverse_callbacks(cb,sensealg,g,cur_time,terminated) = setup_reverse_callbacks(CallbackSet(cb),sensealg,g,cur_time,terminated)
function setup_reverse_callbacks(cb::CallbackSet,sensealg,g,cur_time,terminated)
    cb = CallbackSet(_setup_reverse_callbacks.(cb.continuous_callbacks,(sensealg,),(g,),(cur_time,),(terminated,))...,
                     reverse(_setup_reverse_callbacks.(cb.discrete_callbacks,(sensealg,),(g,),(cur_time,),(terminated,)))...)
    return cb
end

function _setup_reverse_callbacks(cb::Union{ContinuousCallback,DiscreteCallback,VectorContinuousCallback},sensealg,g,loss_ref,terminated)

    if cb isa Union{ContinuousCallback,VectorContinuousCallback} && cb.affect! !== nothing
      cb.affect!.correction.cur_time = loss_ref # set cur_time
      cb.affect!.correction.terminated = terminated # flag if time evolution was terminated by callback
    end

    # ReverseLossCallback adds gradients before and after the callback if save_positions is (true, true).
    # This, however, means that we must check the save_positions setting within the callback.
    # if save_positions = [1,1] is true the loss gradient is accumulated correctly before and after callback.
    # if save_positions = [0,0] no extra gradient is added.
    # if save_positions = [0,1] the gradient contribution is added before the callback but no additional gradient is added afterwards.
    # if save_positions = [1,0] the gradient contribution is added before, and in principle we would need to correct the adjoint state again. Thefore,

    cb.save_positions == [1,0] && error("save_positions=[1,0] is currently not supported.")

    function affect!(integrator)

        indx, pos_neg = get_indx(cb,integrator.t)
        tprev = get_tprev(cb,indx,pos_neg)
        event_idx = cb isa VectorContinuousCallback ? get_event_idx(cb,indx,pos_neg) : nothing

        w = let tprev=tprev, pos_neg=pos_neg, event_idx=event_idx
              function (du,u,p,t)
                _affect! = get_affect!(cb,pos_neg)
                fakeinteg = FakeIntegrator([x for x in u],[x for x in p],t,tprev)
                if cb isa VectorContinuousCallback
                  _affect!(fakeinteg,event_idx)
                else
                  _affect!(fakeinteg)
                end
                du .= fakeinteg.u
              end
        end
        S = integrator.f.f # get the sensitivity function

        # Create a fake sensitivity function to do the vjps
        fakeS = CallbackSensitivityFunction(w,sensealg,S.diffcache,integrator.sol.prob)

        du = first(get_tmp_cache(integrator))
        λ,grad,y,dλ,dgrad,dy = split_states(du,integrator.u,integrator.t,S)

        # if save_positions[2] = false, then the right limit is not saved. Thus, for
        # the QuadratureAdjoint we would need to lift y from the left to the right limit.
        # However, one also needs to update dgrad later on.
        if (sensealg isa QuadratureAdjoint && !cb.save_positions[2]) # || (sensealg isa InterpolatingAdjoint && ischeckpointing(sensealg))
          # lifting for InterpolatingAdjoint is not needed anymore. Callback is already applied.
          w(y,y,integrator.p,integrator.t)
        end

        if cb isa Union{ContinuousCallback,VectorContinuousCallback}
          # correction of the loss function sensitivity for continuous callbacks
          # wrt dependence of event time t on parameters and initial state.
          # Must be handled here because otherwise it is unclear if continuous or
          # discrete callback was triggered.
          @unpack correction = cb.affect!
          @unpack dy_right, Lu_right = correction
          # compute #f(xτ_right,p_right,τ(x₀,p))
          compute_f!(dy_right,S,y,integrator)
          # if callback did not terminate the time evolution, we have to compute one more correction term.
          if cb.save_positions[2] && !correction.terminated
            loss_indx = correction.cur_time[] + 1
            loss_correction!(Lu_right,y,integrator,g,loss_indx)
          else
            Lu_right .*= false
          end
        end

        update_p = copy_to_integrator!(cb,y,integrator.p,integrator.t,indx,pos_neg)
        # reshape u and du (y and dy) to match forward pass (e.g., for matrices as initial conditions). Only needed for BacksolveAdjoint
        if sensealg isa BacksolveAdjoint
          _size = pos_neg ? size(cb.affect!.uleft[indx]) : size(cb.affect_neg!.uleft[indx])
          y = reshape(y, _size)
          dy = reshape(dy, _size)
        end

        if cb isa Union{ContinuousCallback,VectorContinuousCallback}
          # compute the correction of the right limit (with left state limit inserted into dgdt)
          @unpack dy_left, cur_time = correction
          compute_f!(dy_left,S,y,integrator)
          dgdt(dy_left,correction,sensealg,y,integrator,tprev,event_idx)
          if !correction.terminated
            implicit_correction!(Lu_right,dλ,λ,dy_right,correction)
            correction.terminated = false # additional callbacks might have happened which didn't terminate the time evolution
          end
        end

        if update_p
          # changes in parameters
          if !(sensealg isa QuadratureAdjoint)

            wp = let tprev=tprev, pos_neg=pos_neg, event_idx=event_idx
              function (dp,p,u,t)
                _affect! = get_affect!(cb,pos_neg)
                fakeinteg = FakeIntegrator([x for x in u],[x for x in p],t,tprev)
                if cb isa VectorContinuousCallback
                  _affect!(fakeinteg, event_idx)
                else
                  _affect!(fakeinteg)
                end
                dp .= fakeinteg.p
              end
            end
            fakeSp = CallbackSensitivityFunction(wp,sensealg,S.diffcache,integrator.sol.prob)
            #vjp with Jacobin given by dw/dp before event and vector given by grad
            vecjacobian!(dgrad, integrator.p, grad, y, integrator.t, fakeSp;
                                  dgrad=nothing, dy=nothing)
            grad .= dgrad
          end
        end

        vecjacobian!(dλ, y, λ, integrator.p, integrator.t, fakeS;
                              dgrad=dgrad, dy=dy)

        if cb isa Union{ContinuousCallback,VectorContinuousCallback}
          # second correction to correct for left limit
          @unpack Lu_left = correction
          implicit_correction!(Lu_left,dλ,dy_left,correction)
          dλ .+= Lu_left - Lu_right

          if cb.save_positions[1] == true
            # if the callback saved the first position, we need to implicitly correct this value as well
            loss_indx = correction.cur_time[]
            implicit_correction!(Lu_left,dy_left,correction,y,integrator,g,loss_indx)
            dλ .+= Lu_left
          end
        end

        λ .= dλ

        if !(sensealg isa QuadratureAdjoint)
          grad .-= dgrad
        end

    end

    times = if typeof(cb) <: DiscreteCallback
        cb.affect!.event_times
    else
        [cb.affect!.event_times;cb.affect_neg!.event_times]
    end

    PresetTimeCallback(times,
                       affect!,
                       save_positions = (false,false))
end

get_indx(cb::DiscreteCallback,t) = (searchsortedfirst(cb.affect!.event_times,t), true)
function get_indx(cb::Union{ContinuousCallback,VectorContinuousCallback}, t)
  if !isempty(cb.affect!.event_times) || !isempty(cb.affect_neg!.event_times)
    indx = searchsortedfirst(cb.affect!.event_times,t)
    indx_neg = searchsortedfirst(cb.affect_neg!.event_times,t)
    if !isempty(cb.affect!.event_times) && cb.affect!.event_times[min(indx,length(cb.affect!.event_times))]==t
      return indx, true
    elseif !isempty(cb.affect_neg!.event_times) && cb.affect_neg!.event_times[min(indx_neg,length(cb.affect_neg!.event_times))]==t
      return indx_neg, false
    else
      error("Event was triggered but no corresponding event in ContinuousCallback was found. Please report this error.")
    end
  else
    error("No event was recorded. Please report this error.")
  end
end

get_tprev(cb::DiscreteCallback,indx,bool) = cb.affect!.tprev[indx]
function get_tprev(cb::Union{ContinuousCallback,VectorContinuousCallback}, indx, bool)
  if bool
    return cb.affect!.tprev[indx]
  else
    return cb.affect_neg!.tprev[indx]
  end
end

function get_event_idx(cb::VectorContinuousCallback, indx, bool)
  if bool
    return cb.affect!.event_idx[indx]
  else
    return cb.affect_neg!.event_idx[indx]
  end
end

function copy_to_integrator!(cb::DiscreteCallback, y, p, t, indx, bool)
  copyto!(y, cb.affect!.uleft[indx])
  update_p = (p != cb.affect!.pleft[indx])
  update_p && copyto!(p, cb.affect!.pleft[indx])
  update_p
end

function copy_to_integrator!(cb::Union{ContinuousCallback,VectorContinuousCallback}, y, p, t, indx, bool)
  if bool
    copyto!(y, cb.affect!.uleft[indx])
    update_p = (p != cb.affect!.pleft[indx])
    update_p && copyto!(p, cb.affect!.pleft[indx])
  else
    copyto!(y, cb.affect_neg!.uleft[indx])
    update_p = (p != cb.affect_neg!.pleft[indx])
    update_p && copyto!(p, cb.affect_neg!.pleft[indx])
  end
  update_p
end

function compute_f!(dy,S,y,integrator)
  p, t = integrator.p, integrator.t

  if inplace_sensitivity(S)
    S.f(dy,y,p,t)
  else
    dy[:] .= S.f(y,p,t)
  end
  return nothing
end

function dgdt(dy,correction,sensealg,y,integrator,tprev,event_idx)
  # dy refers to f evaluated on left limit
  @unpack gt_val, gu_val, gt, gu, gt_conf, gu_conf, condition = correction

  p, t = integrator.p, integrator.t

  fakeinteg = FakeIntegrator([x for x in y],p,t,tprev)

  # derivative and gradient of condition with respect to time and state, respectively
  gt.u = y
  gt.integrator = fakeinteg

  gu.t = t
  gu.integrator = fakeinteg

  # for VectorContinuousCallback we also need to set the event_idx.
  if gt isa VectorConditionTimeWrapper
    gt.event_idx = event_idx
    gu.event_idx = event_idx

    # safety check: evaluate condition to check if several conditions were true.
    # This is currently not supported
    condition(gt.out_cache,y,t,integrator)
    gt.out_cache .= abs.(gt.out_cache) .< 1000*eps(eltype(gt.out_cache))
    (sum(gt.out_cache)!=1 || gt.out_cache[event_idx]!=1) && error("Either several events were triggered or `event_idx` was falsely identified. Output of conditions $(gt.out_cache)")
  end

  derivative!(gt_val, gt, t, sensealg, gt_conf)
  gradient!(gu_val, gu, y, sensealg, gu_conf)

  gt_val .+= dot(gu_val,dy)
  @. gt_val = inv(gt_val) # allocates?

  @. gu_val *= -gt_val
  return nothing
end

function loss_correction!(Lu,y,integrator,g,indx)
  # ∂L∂t correction should be added if L depends explicitly on time.
  p, t = integrator.p, integrator.t
  g(Lu,y,p,t,indx)
  return nothing
end

function implicit_correction!(Lu,dλ,λ,dy,correction)
  @unpack gu_val = correction

  # remove gradients from adjoint state to compute correction factor
  @. dλ = λ - Lu
  Lu .= dot(dλ,dy)*gu_val

  return nothing
end

function implicit_correction!(Lu,λ,dy,correction)
  @unpack gu_val = correction

  Lu .= dot(λ,dy)*gu_val

  return nothing
end

function implicit_correction!(Lu,dy,correction,y,integrator,g,indx)
  @unpack gu_val = correction

  p, t = integrator.p, integrator.t

  # loss function gradient (not condition!)
  # ∂L∂t correction should be added, also ∂L∂p is missing.
  # correct adjoint
  g(Lu,y,p,t,indx)

  Lu .= dot(Lu,dy)*gu_val

  # note that we don't add the gradient Lu here again to the correction because it will be added by the ReverseLossCallback.
  return nothing
end

# ConditionTimeWrapper: Wrapper for implicit correction for ContinuousCallback
# VectorConditionTimeWrapper: Wrapper for implicit correction for VectorContinuousCallback
function build_condition_wrappers(cb::ContinuousCallback,condition,u,t,fakeinteg)
  gt = ConditionTimeWrapper(condition,u,fakeinteg)
  gu = ConditionUWrapper(condition,t,fakeinteg)
  return gt, gu
end
function build_condition_wrappers(cb::VectorContinuousCallback,condition,u,t,fakeinteg)
  out = similar(u, cb.len) # create a cache for condition function (out,u,t,integrator)
  gt = VectorConditionTimeWrapper(condition,u,fakeinteg,1,out)
  gu = VectorConditionUWrapper(condition,t,fakeinteg,1,out)
  return gt, gu
end
mutable struct ConditionTimeWrapper{F,uType,Integrator} <: Function
  f::F
  u::uType
  integrator::Integrator
end
(ff::ConditionTimeWrapper)(t) = [ff.f(ff.u,t,ff.integrator)]
mutable struct ConditionUWrapper{F,tType,Integrator} <: Function
  f::F
  t::tType
  integrator::Integrator
end
(ff::ConditionUWrapper)(u) = ff.f(u,ff.t,ff.integrator)
mutable struct VectorConditionTimeWrapper{F,uType,Integrator,outType} <: Function
  f::F
  u::uType
  integrator::Integrator
  event_idx::Int
  out_cache::outType
end
(ff::VectorConditionTimeWrapper)(t) = (ff.f(ff.out_cache,ff.u,t,ff.integrator); [ff.out_cache[ff.event_idx]])

mutable struct VectorConditionUWrapper{F,tType,Integrator,outType} <: Function
  f::F
  t::tType
  integrator::Integrator
  event_idx::Int
  out_cache::outType
end
(ff::VectorConditionUWrapper)(u) = (out = similar(u,length(ff.out_cache)); ff.f(out,u,ff.t,ff.integrator); out[ff.event_idx])

DiffEqBase.terminate!(i::FakeIntegrator) = nothing

# get the affect function of the callback. For example, allows us to get the `f` in PeriodicCallback without the integrator.tstops handling.
get_affect!(cb::DiscreteCallback,bool) = get_affect!(cb.affect!)
get_affect!(cb::Union{ContinuousCallback,VectorContinuousCallback},bool) = bool ? get_affect!(cb.affect!) : get_affect!(cb.affect_neg!)
get_affect!(affect!::TrackedAffect) = get_affect!(affect!.affect!)
get_affect!(affect!) = affect!
get_affect!(affect!::DiffEqCallbacks.PeriodicCallbackAffect) = affect!.affect!


text:
# # Plankton mixing and blooming
#
# In this example, we simulate the mixing of phytoplankton by convection
# that decreases in time and eventually shuts off, thereby precipitating a
# phytoplankton bloom. A similar scenario was simulated by
# [Taylor and Ferrari (2011)](https://aslopubs.onlinelibrary.wiley.com/doi/abs/10.4319/lo.2011.56.6.2293),
# providing evidence that the
# ["critical turbulence hypothesis"](https://en.wikipedia.org/wiki/Critical_depth#Critical_Turbulence_Hypothesis)
# explains the explosive bloom of oceanic phytoplankton
# observed in spring.
#
# The phytoplankton in our model are advected, diffuse, grow, and die according to
#
# ```math
# ∂_t P + \boldsymbol{u ⋅ ∇} P - κ ∇²P = (μ₀ \exp(z / λ) - m) \, P \, ,
# ```
#
# where ``\boldsymbol{v}`` is the turbulent velocity field, ``κ`` is an isotropic diffusivity,
#  ``μ₀`` is the phytoplankton growth rate at the surface, ``λ`` is the scale over
# which sunlight attenuates away from the surface, and ``m`` is the mortality rate
# of phytoplankton due to viruses and grazing by zooplankton. We use Oceananigans'
#  `Forcing` abstraction to implement the phytoplankton dynamics described by the right
# side of the phytoplankton equation above.
#
# This example demonstrates
#
#   * How to use a user-defined forcing function to
#     simulate the dynamics of phytoplankton growth in sunlight
#     and grazing by zooplankton.
#   * How to set time-dependent boundary conditions.
#   * How to use the `TimeStepWizard` to adapt the simulation time-step.
#   * How to use `AveragedField` to diagnose spatial averages of model fields.
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add Oceananigans, Plots, JLD2, Measures"
# ```

# ## The grid
#
# We use a two-dimensional grid with 64² points and 1 m grid spacing and assign `Flat`
# to the `y` direction:

using Oceananigans
using Oceananigans.Units: minutes, hour, hours, day

grid = RegularRectilinearGrid(size=(64, 64), extent=(64, 64), topology=(Periodic, Flat, Bounded))

# ## Boundary conditions
#
# We impose a surface buoyancy flux that's initially constant and then decays to zero,

buoyancy_flux(x, y, t, p) = p.initial_buoyancy_flux * exp(-t^4 / (24 * p.shut_off_time^4))

buoyancy_flux_parameters = (initial_buoyancy_flux = 1e-8, # m² s⁻³
                                    shut_off_time = 2hours)

buoyancy_flux_bc = FluxBoundaryCondition(buoyancy_flux, parameters = buoyancy_flux_parameters)

# The fourth power in the argument of `exp` above helps keep the buoyancy flux relatively
# constant during the first phase of the simulation. We produce a plot of this time-dependent
# buoyancy flux for the visually-oriented,

using Plots, Measures

time = range(0, 12hours, length=100)

flux_plot = plot(time ./ hour, [buoyancy_flux(0, 0, t, buoyancy_flux_parameters) for t in time],
                 linewidth = 2, xlabel = "Time (hours)", ylabel = "Surface buoyancy flux (m² s⁻³)",
                 size = (800, 300), margin = 5mm, label = nothing)

# The buoyancy flux effectively shuts off after 6 hours of simulation time.
#
# !!! info "The flux convention in Oceananigans.jl"
#     Fluxes are defined by the direction a quantity is carried: _positive_ velocities
#     produce _positive_ fluxes, while _negative_ velocities produce _negative_ fluxes.
#     Diffusive fluxes are defined with the same convention. A positive flux at the _top_
#     boundary transports buoyancy _upwards, out of the domain_. This means that a positive
#     flux of buoyancy at the top boundary reduces the buoyancy of near-surface fluid,
#     causing convection.
#
# The initial condition and bottom boundary condition impose the constant buoyancy gradient

N² = 1e-4 # s⁻²

buoyancy_gradient_bc = GradientBoundaryCondition(N²)

# In summary, the buoyancy boundary conditions impose a destabilizing flux
# at the top and a stable buoyancy gradient at the bottom:

buoyancy_bcs = TracerBoundaryConditions(grid, top = buoyancy_flux_bc, bottom = buoyancy_gradient_bc)

# ## Phytoplankton dynamics: light-dependent growth and uniform mortality
#
# We use a simple model for the growth of phytoplankton in sunlight and decay
# due to viruses and grazing by zooplankton,

growing_and_grazing(x, y, z, t, P, p) = (p.μ₀ * exp(z / p.λ) - p.m) * P
nothing # hide

# with parameters

plankton_dynamics_parameters = (μ₀ = 1/day,   # surface growth rate
                                 λ = 5,       # sunlight attenuation length scale (m)
                                 m = 0.1/day) # mortality rate due to virus and zooplankton grazing

# We tell `Forcing` that our plankton model depends
# on the plankton concentration `P` and the chosen parameters,

plankton_dynamics = Forcing(growing_and_grazing, field_dependencies = :P,
                            parameters = plankton_dynamics_parameters)

# ## The model
#
# The name "`P`" for phytoplankton is specified in the
# constructor for `IncompressibleModel`. We additionally specify a fifth-order
# advection scheme, third-order Runge-Kutta time-stepping, isotropic viscosity and diffusivities,
# and Coriolis forces appropriate for planktonic convection at mid-latitudes on Earth.

model = IncompressibleModel(
                   grid = grid,
              advection = UpwindBiasedFifthOrder(),
            timestepper = :RungeKutta3,
                closure = IsotropicDiffusivity(ν=1e-4, κ=1e-4),
               coriolis = FPlane(f=1e-4),
                tracers = (:b, :P), # P for Plankton
               buoyancy = BuoyancyTracer(),
                forcing = (P=plankton_dynamics,),
    boundary_conditions = (b=buoyancy_bcs,)
)

# ## Initial condition
#
# We set the initial phytoplankton at ``P = 1 \, \rm{μM}``.
# For buoyancy, we use a stratification that's mixed near the surface and
# linearly stratified below, superposed with surface-concentrated random noise.

mixed_layer_depth = 32 # m

stratification(z) = z < -mixed_layer_depth ? N² * z : - N² * mixed_layer_depth

noise(z) = 1e-4 * N² * grid.Lz * randn() * exp(z / 4)

initial_buoyancy(x, y, z) = stratification(z) + noise(z)

set!(model, b=initial_buoyancy, P=1)

# ## Adaptive time-stepping, logging, output and simulation setup
#
# We use a `TimeStepWizard` that limits the
# time-step to 2 minutes, and adapts the time-step such that CFL
# (Courant-Freidrichs-Lewy) number hovers around `1.0`,

wizard = TimeStepWizard(cfl=1.0, Δt=2minutes, max_change=1.1, max_Δt=2minutes)

# We also write a function that prints the progress of the simulation

using Printf

progress(sim) = @printf("Iteration: %d, time: %s, Δt: %s\n",
                        sim.model.clock.iteration,
                        prettytime(sim.model.clock.time),
                        prettytime(sim.Δt.Δt))

simulation = Simulation(model, Δt=wizard, stop_time=24hours,
                        iteration_interval=20, progress=progress)

# We add a basic `JLD2OutputWriter` that writes velocities and both
# the two-dimensional and horizontally-averaged plankton concentration,

averaged_plankton = AveragedField(model.tracers.P, dims=(1, 2))

outputs = (w = model.velocities.w,
           plankton = model.tracers.P,
           averaged_plankton = averaged_plankton)

simulation.output_writers[:simple_output] =
    JLD2OutputWriter(model, outputs,
                     schedule = TimeInterval(20minutes),
                     prefix = "convecting_plankton",
                     force = true)

# !!! info "Using multiple output writers"
#     Because each output writer is associated with a single output `schedule`,
#     it often makes sense to use _different_ output writers for different types of output.
#     For example, reduced fields like `AveragedField` usually consume less disk space than
#     two- or three-dimensional fields, and can thus be output more frequently without
#     blowing up your hard drive. An arbitrary number of output writers may be added to
#     `simulation.output_writers`.
#
# The simulation is set up. Let there be plankton:

run!(simulation)

# Notice how the time-step is reduced at early times, when turbulence is strong,
# and increases again towards the end of the simulation when turbulence fades.

# ## Visualizing the solution
#
# We'd like to a make a plankton movie. First we load the output file
# and build a time-series of the buoyancy flux,

using JLD2

file = jldopen(simulation.output_writers[:simple_output].filepath)

iterations = parse.(Int, keys(file["timeseries/t"]))

times = [file["timeseries/t/$iter"] for iter in iterations]

buoyancy_flux_time_series = [buoyancy_flux(0, 0, t, buoyancy_flux_parameters) for t in times]
nothing # hide

# and then we construct the ``x, z`` grid,

xw, yw, zw = nodes(model.velocities.w)
xp, yp, zp = nodes(model.tracers.P)
nothing # hide

# Finally, we animate plankton mixing and blooming,

using Plots

@info "Making a movie about plankton..."

w_lim = 0   # the maximum(abs(w)) across the whole timeseries

for (i, iteration) in enumerate(iterations)
    w = file["timeseries/w/$iteration"][:, 1, :]

    global w_lim = maximum([w_lim, maximum(abs.(w))])
end

anim = @animate for (i, iteration) in enumerate(iterations)

    @info "Plotting frame $i from iteration $iteration..."

    t = file["timeseries/t/$iteration"]
    w = file["timeseries/w/$iteration"][:, 1, :]
    P = file["timeseries/plankton/$iteration"][:, 1, :]
    averaged_P = file["timeseries/averaged_plankton/$iteration"][1, 1, :]

    P_min = minimum(P) - 1e-9
    P_max = maximum(P) + 1e-9
    P_lims = (0.95, 1.1)

    w_levels = range(-w_lim, stop=w_lim, length=20)

    P_levels = collect(range(P_lims[1], stop=P_lims[2], length=20))
    P_lims[1] > P_min && pushfirst!(P_levels, P_min)
    P_lims[2] < P_max && push!(P_levels, P_max)

    kwargs = (xlabel="x (m)", ylabel="y (m)", aspectratio=1, linewidth=0, colorbar=true,
              xlims=(0, model.grid.Lx), ylims=(-model.grid.Lz, 0))

    w_contours = contourf(xw, zw, w';
                          color = :balance,
                          levels = w_levels,
                          clims = (-w_lim, w_lim),
                          kwargs...)

    P_contours = contourf(xp, zp, clamp.(P, P_lims[1], P_lims[2])';
                          color = :matter,
                          levels = P_levels,
                          clims = P_lims,
                          kwargs...)

    P_profile = plot(averaged_P, zp,
                     linewidth = 2,
                     label = nothing,
                     xlims = (0.9, 1.3),
                     ylabel = "z (m)",
                     xlabel = "Plankton concentration (μM)")

    flux_plot = plot(times ./ hour, buoyancy_flux_time_series,
                     linewidth = 1,
                     label = "Buoyancy flux time series",
                     color = :black,
                     alpha = 0.4,
                     legend = :topright,
                     xlabel = "Time (hours)",
                     ylabel = "Buoyancy flux (m² s⁻³)",
                     ylims = (0.0, 1.1 * buoyancy_flux_parameters.initial_buoyancy_flux))

    plot!(flux_plot, times[1:i] ./ hour, buoyancy_flux_time_series[1:i],
          color = :steelblue,
          linewidth = 6,
          label = nothing)

    scatter!(flux_plot, times[i:i] / hour, buoyancy_flux_time_series[i:i],
             markershape = :circle,
             color = :steelblue,
             markerstrokewidth = 0,
             markersize = 15,
             label = "Current buoyancy flux")

    layout = Plots.grid(2, 2, widths=(0.7, 0.3))

    w_title = @sprintf("Vertical velocity (m s⁻¹) at %s", prettytime(t))
    P_title = @sprintf("Plankton concentration (μM) at %s", prettytime(t))

    plot(w_contours, flux_plot, P_contours, P_profile,
         title=[w_title "" P_title ""],
         layout=layout, size=(1000.5, 1000.5))
end

mp4(anim, "convecting_plankton.mp4", fps = 8) # hide


text:
# Autogenerated wrapper script for lrslib_jll for x86_64-w64-mingw32
export liblrs, liblrsnash, lrs, lrsnash

using GMP_jll
JLLWrappers.@generate_wrapper_header("lrslib")
JLLWrappers.@declare_library_product(liblrs, "liblrs-0.dll")
JLLWrappers.@declare_library_product(liblrsnash, "liblrsnash.dll")
JLLWrappers.@declare_executable_product(lrs)
JLLWrappers.@declare_executable_product(lrsnash)
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll)
    JLLWrappers.@init_library_product(
        liblrs,
        "bin\\liblrs-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        liblrsnash,
        "bin\\liblrsnash.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        lrs,
        "bin\\lrs.exe",
    )

    JLLWrappers.@init_executable_product(
        lrsnash,
        "bin\\lrsnash.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Conda
using BinDeps
using Base.Test
using Compat

env = Conda.Environment(:test_conda_jl)
@test Conda.exists("curl", env)
Conda.add("curl", env)

if is_unix()
    curl_path = joinpath(Conda.prefix(env), "bin", "curl-config")
end
if is_windows()
    curl_path = joinpath(Conda.lib_dir(env), "curl.exe")
end

@test isfile(curl_path)

@test isfile(joinpath(Conda.bin_dir(env), basename(curl_path)))

Conda.rm("curl", env)
if is_unix()
    @test !isfile(curl_path)
end

@test isfile(Conda.conda_bin(env))
Conda.add("python", env)
@test isfile(joinpath(Conda.python_dir(env), "python" * (is_windows() ? ".exe": "")))

channels = Conda.channels()
@test channels == ["defaults"]

Conda.add_channel("foo", env)
@test Conda.channels(env) == ["foo", "defaults"]
# Testing that calling the function twice do not fail
Conda.add_channel("foo", env)

Conda.rm_channel("foo", env)
channels = Conda.channels(env)
@test channels == ["defaults"]


text:
# -*- coding: utf-8 -*-
# AES tests from:
# https://git.lysator.liu.se/nettle/nettle/blob/master/testsuite/aes-test.c

# AES 128
for (key,text,encrypted) in [
    (
        "00010203050607080A0B0C0D0F101112",
        "506812A45F08C889B97F5980038B8359",
        "D8F532538289EF7D06B506A4FD5BE9C9"
    ),(
        "14151617191A1B1C1E1F202123242526",
        "5C6D71CA30DE8B8B00549984D2EC7D4B",
        "59AB30F4D4EE6E4FF9907EF65B1FB68C"
    ),(
        "28292A2B2D2E2F30323334353738393A",
        "53F3F4C64F8616E4E7C56199F48F21F6",
        "BF1ED2FCB2AF3FD41443B56D85025CB1",
    ),(
        "A0A1A2A3A5A6A7A8AAABACADAFB0B1B2",
        "F5F4F7F684878689A6A7A0A1D2CDCCCF",
        "CE52AF650D088CA559425223F4D32694",
    ),(
        "2b7e151628aed2a6abf7158809cf4f3c",
        "6bc1bee22e409f96e93d7e117393172aae2d8a571e03ac9c9eb76fac45af8e5130c81c46a35ce411e5fbc1191a0a52eff69f2445df4f9b17ad2b417be66c3710",
        "3ad77bb40d7a3660a89ecaf32466ef97f5d3d58503b9699de785895a96fdbaaf43b1cd7f598ece23881b00e3ed0306887b0c785e27e8ad3f8223207104725dd4",
    )
]
    @test encrypt("aes128", hex2bytes(key), hex2bytes(text)) == hex2bytes(encrypted)
    @test decrypt("aes128", hex2bytes(key), hex2bytes(encrypted)) == hex2bytes(text)
end

# AES192
for (key,text,encrypted) in [
    (
        "00010203050607080A0B0C0D0F10111214151617191A1B1C",
        "2D33EEF2C0430A8A9EBF45E809C40BB6",
        "DFF4945E0336DF4C1C56BC700EFF837F",
    ),(
        "8e73b0f7da0e6452c810f32b809079e562f8ead2522c6b7b",
        "6bc1bee22e409f96e93d7e117393172aae2d8a571e03ac9c9eb76fac45af8e5130c81c46a35ce411e5fbc1191a0a52eff69f2445df4f9b17ad2b417be66c3710",
        "bd334f1d6e45f25ff712a214571fa5cc974104846d0ad3ad7734ecb3ecee4eefef7afd2270e2e60adce0ba2face6444e9a4b41ba738d6c72fb16691603c18e0e",
    )
]
    @test encrypt("aes192", hex2bytes(key), hex2bytes(text)) == hex2bytes(encrypted)
    @test decrypt("aes192", hex2bytes(key), hex2bytes(encrypted)) == hex2bytes(text)
end

# AES256
for (key,text,encrypted) in [
    (
        "00010203050607080A0B0C0D0F10111214151617191A1B1C1E1F202123242526",
        "834EADFCCAC7E1B30664B1ABA44815AB",
        "1946DABF6A03A2A2C3D0B05080AED6FC",
    ),(
        "603deb1015ca71be2b73aef0857d77811f352c073b6108d72d9810a30914dff4",
        "6bc1bee22e409f96e93d7e117393172aae2d8a571e03ac9c9eb76fac45af8e5130c81c46a35ce411e5fbc1191a0a52eff69f2445df4f9b17ad2b417be66c3710",
        "f3eed1bdb5d2a03c064b5a7e3db181f8591ccb10d410ed26dc5ba74a31362870b6ed21b99ca6f4f9f153e7b1beafed1d23304b7a39f9f3ff067d8d8f9e24ecc7",
    )
]
    @test encrypt("aes256", hex2bytes(key), hex2bytes(text)) == hex2bytes(encrypted)
    @test decrypt("aes256", hex2bytes(key), hex2bytes(encrypted)) == hex2bytes(text)
end

# Test lower-level API
key = "this key's exactly 32 bytes long"
enc = Encryptor("AES256", key)
plaintext = "this is 16 chars"
ciphertext = encrypt(enc, Vector{UInt8}(plaintext))

dec = Decryptor("AES256", key)
deciphertext = decrypt(dec, ciphertext)
@test Vector{UInt8}(plaintext) == deciphertext # no bytestring

willcauseassertion = "this is 16 (∀).." # case of length(::String) == 16
@test length(willcauseassertion) == 16
@test sizeof(willcauseassertion) == 18
# @test_throws AssertionError Vector{UInt8}(willcauseassertion) == decrypt(dec, encrypt(enc, Vector{UInt8}(willcauseassertion))) # can not catch this c assertion
@test_throws ArgumentError Vector{UInt8}(willcauseassertion) == decrypt(dec, encrypt(enc, Vector{UInt8}(willcauseassertion)))
# @test_throws AssertionError Vector{UInt8}(willcauseassertion) == decrypt(dec, encrypt(enc, Vector{UInt8}(willcauseassertion))) # can not catch this c assertion
@test_throws ArgumentError Vector{UInt8}(willcauseassertion) == decrypt(dec, encrypt(enc, Vector{UInt8}(willcauseassertion)))

willbebroken = "this is 16 (∀)" # case of length(::String) != 16
@test length(willbebroken) == 14
@test sizeof(willbebroken) == 16
@test Vector{UInt8}(willbebroken) == decrypt(dec, encrypt(enc, Vector{UInt8}(willbebroken)))
@test willbebroken == String(decrypt(dec, encrypt(enc, Vector{UInt8}(willbebroken))))
@test Vector{UInt8}(willbebroken) == decrypt(dec, encrypt(enc, Vector{UInt8}(willbebroken)))
@test willbebroken == String(decrypt(dec, encrypt(enc, Vector{UInt8}(willbebroken))))

criticalbytes = hex2bytes("6e6f74555446382855aa552de2888029")
@test length(criticalbytes) == 16
@test sizeof(criticalbytes) == 16
@test criticalbytes == decrypt(dec, encrypt(enc, criticalbytes))

# This one will pass, but may be caught UnicodeError exception when evaluate it by julia ide.
dummy = String(decrypt(dec, encrypt(enc, criticalbytes)))
@test isa(dummy, AbstractString)
if !isdefined(Core, :String) || !isdefined(Core, :AbstractString)
    @test !isa(dummy, ASCIIString)
end
@test isa(dummy, String) # gray zone
@test sizeof(dummy) == 16
@test length(Vector{UInt8}(dummy)) == 16
@test length(dummy) != 16

# Test errors
@test_throws ArgumentError Encryptor("this is not a cipher", key)
@test_throws ArgumentError Decryptor("this is not a cipher either", key)

@test_throws ArgumentError Encryptor("AES256", "bad key length")
@test_throws ArgumentError Decryptor("AES256", "bad key length")

enc = Encryptor("AES256", key)
dec = Decryptor("AES256", key)
result = Vector{UInt8}(undef, sizeof(plaintext) - 1)
@test_throws ArgumentError encrypt!(enc, result, plaintext)
@test_throws ArgumentError decrypt!(dec, result, plaintext)

# Test show methods
println("Testing cipher show methods:")
println(get_cipher_types()["AES256"])
println(Encryptor("AES256", key))
println(Decryptor("AES256", key))

println("Testing cipher AES256CBC:")
# AES256CBC
for (pw,salt,iv,key,text,encrypted) in [
    (
        Vector(b"Secret Passphrase"),
        "a3e550e89e70996c",
        "7c7ed9434ddb9c2d1e1fcc38b4bf4667",
        "e299ff9d8e4831f07e5323913c53e5f0fec3a040a211d6562fa47607244d0051",
        "4d657373616765",
        "da8aab1b904205a7e49c1ecc7118a8f4",
    ),(
        Vector(b"Secret Passphrase"),
        "a3e550e89e70996c",
        "7c7ed9434ddb9c2d1e1fcc38b4bf4667",
        "e299ff9d8e4831f07e5323913c53e5f0fec3a040a211d6562fa47607244d0051",
        "4d657373616765090909090909090909",
        "da8aab1b904205a7e49c1ecc7118a8f4804bef7be79216196739de7845da182d",
    )
]
    (key32, iv16) = (hex2bytes(key), hex2bytes(iv))
    @test gen_key32_iv16(pw, hex2bytes(salt)) == (key32, iv16)
    @test encrypt("aes256", :CBC, iv16, key32, add_padding_PKCS5(hex2bytes(text), 16)) == hex2bytes(encrypted)
    @test trim_padding_PKCS5(decrypt("aes256", :CBC, iv16, key32, hex2bytes(encrypted))) == hex2bytes(text)
end

# Test errors
badkey = "this key's exactly 32(∪∩∀ДЯ)...."
@test length(badkey) == 32
@test sizeof(badkey) == 40
@test_throws ArgumentError Encryptor("AES256", badkey)
@test_throws ArgumentError Decryptor("AES256", badkey)

iv = Vector(b"this is 16 chars")
key = Vector(b"this key's exactly 32 bytes long")
enc = Encryptor("AES256", key)
dec = Decryptor("AES256", key)
shortresult = Vector{UInt8}(undef, sizeof(plaintext) - 1)
@test_throws ArgumentError encrypt!(enc, :CBC, iv, shortresult, plaintext)
@test_throws ArgumentError decrypt!(dec, :CBC, iv, shortresult, plaintext)
result = Vector{UInt8}(undef, sizeof(plaintext))
shorttext = Vector{UInt8}(undef, sizeof(plaintext) - 1)
@test_throws ArgumentError encrypt!(enc, :CBC, iv, result, shorttext)
@test_throws ArgumentError decrypt!(dec, :CBC, iv, result, shorttext)
longresult = Vector{UInt8}(undef, sizeof(plaintext) + 1)
longtext = Vector{UInt8}(undef, sizeof(plaintext) + 1)
@test_throws ArgumentError encrypt!(enc, :CBC, iv, longresult, longtext)
@test_throws ArgumentError decrypt!(dec, :CBC, iv, longresult, longtext)
shortiv = Vector{UInt8}(undef, sizeof(iv) - 1)
@test_throws ArgumentError encrypt!(enc, :CBC, shortiv, result, plaintext)
@test_throws ArgumentError decrypt!(dec, :CBC, shortiv, result, plaintext)
longiv = Vector{UInt8}(undef, sizeof(iv) + 1)
@test_throws ArgumentError encrypt!(enc, :CBC, longiv, result, plaintext)
@test_throws ArgumentError decrypt!(dec, :CBC, longiv, result, plaintext)

# encrypt!(enc, :GCM, iv, result, plaintext)
# decrypt!(dec, :CCM, iv, result, plaintext)
@test_throws ArgumentError encrypt!(enc, :UNKNOWN, iv, result, plaintext)
@test_throws ArgumentError decrypt!(dec, :UNKNOWN, iv, result, plaintext)

println("Cipher AES256CBC OK.")


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMQuad.jl/blob/master/LICENSE

# Tensorial quadrature rules in 1, 2, 3 dimensions

function tensor_product(w::Tuple, p::Tuple, dim::Int)
    @assert length(w) == length(p)

    N = length(w)
    weights = Float64[]
    points = NTuple{dim, Float64}[]
    for i in CartesianIndices(ntuple(i -> 1:N, dim))
        push!(weights, prod(w[k] for k in Tuple(i)))
        push!(points, ntuple(k -> p[i[k]], dim))
    end
    return zip(Tuple(weights), Tuple(points))
end

names = [:GLSEG, :GLQUAD, :GLHEX]
names2 = ["segment", "quadrilateral", "hexahedron"]

for n in 1:length(QUAD_DATA)
    points, weights = QUAD_DATA[n]
    order = 2(n-1)+1
    for dim in (1, 2, 3)
        n_points = length(points)^dim
        quadname = QuoteNode(Symbol(string(names[dim], n_points)))
        z = tensor_product(weights, points, dim)
        @eval begin
            @doc """
                get_quadrature_points(::Type{Val{:$($(quadname))})

            Gauss-Legendre quadrature, $($(n_points)) point rule on $($(names2[dim]))."""
            get_quadrature_points(::Type{Val{$(quadname)}}) = $z
        end
        @eval get_order(::Type{Val{$(quadname)}}) = $order
    end
end


text:
# This file is not an active part of the package. This is the code
# that uses the Clang.jl package to wrap Hypre using the headers.

# Find all headers
# To run the script from Julia console:
# include(joinpath(Pkg.dir("Hypre"), "src", "wrap_hypre.jl"));
push!(Libdl.DL_LOAD_PATH, "/usr/local/Cellar/llvm/4.0.0_1/lib")
using Clang.wrap_c
const wdir = dirname(@__FILE__)
const pkgbasedir = joinpath(wdir, "..")

# `outpath` specifies, where the julian wrappers would be generated.
# If the generated .jl files are ok, they have to be copied to the "src" folder
# overwriting the old ones
const outpath = joinpath(pkgbasedir, "new")
rm(outpath, recursive = true, force = true)
mkpath(outpath)
mkpath(outpath)

const incpath = joinpath(pkgbasedir, "deps", "usr", "include")
if !isdir(incpath)
  error("Run Pkg.build(\"Hypre\") before trying to wrap C headers.")
end

info("Scanning Hypre headers in $incpath...")
const Hypre_header_files = ["sstruct_ls", "struct_ls", "parcsr_ls", "DistributedMatrixPilutSolver_types", "DistributedMatrixPilutSolver_protos", "lobpcg"]
const Hypre_headers =[joinpath(incpath, y) for y in (string("HYPRE_", x, ".h") for x in Hypre_header_files)]
# map(x->joinpath(incpath, x), readdir(incpath))

const clang_path = "/usr/local/Cellar/llvm/4.0.0_1/lib/clang/4.0.0/" # change to your clang location
const includes = [
    joinpath(clang_path, "include"),
    incpath,
    "/usr/local/Cellar/open-mpi/2.1.0/include/"
]

function find_outfile(s)
    joinpath(outpath, string(first(splitext(basename(s))), ".jl"))
end
find_libfile(s) = "libHypre"

const context = wrap_c.init(
    headers = Hypre_headers,
    common_file = joinpath(outpath, "types_and_consts.jl"),
    clang_args = [
        "-D", "__STDC_LIMIT_MACROS",
        "-D", "__STDC_CONSTANT_MACROS",
        # "-v"
    ],
    # clang_diagnostics = true,
    header_library = find_libfile,
    header_outputfile = find_outfile,
    clang_includes = includes,
)

info("Generating .jl wrappers for Hypre in $outpath...")
run(context)
info("Done generating .jl wrappers for Hypre in $outpath")


text:
## Script which prints out the hostnames of the worker processes
# Used in conjunction with a job script to ensure multi-node parallelism

hosts = @parallel for i=1:120
       run(`hostname`) end
println(hosts)


text:
using CpuId

register_size = simdbytes()
register_count = cpufeature(CpuId.AVX512F) ? 32 : 16
register_capacity = register_size * register_count
cache_size = CpuId.cachesize()
num_cores = CpuId.cpucores()

cpu_info_string = """
const REGISTER_SIZE = $register_size
const REGISTER_COUNT = $register_count
const REGISTER_CAPACITY = $register_capacity
const FP256 = $(cpufeature(CpuId.FP256)) # Is AVX2 fast?
const CACHELINE_SIZE = $(cachelinesize())
const CACHE_SIZE = $cache_size
const NUM_CORES = $num_cores
const FMA3 = $(cpufeature(CpuId.FMA3))
const AVX2 = $(cpufeature(CpuId.AVX2))
const AVX512F = $(cpufeature(CpuId.AVX512F))
const AVX512ER = $(cpufeature(CpuId.AVX512ER))
const AVX512PF = $(cpufeature(CpuId.AVX512PF))
const AVX512VL = $(cpufeature(CpuId.AVX512VL))
const AVX512BW = $(cpufeature(CpuId.AVX512BW))
const AVX512DQ = $(cpufeature(CpuId.AVX512DQ))
const AVX512CD = $(cpufeature(CpuId.AVX512CD))

"""

open(joinpath(@__DIR__, "..", "src", "cpu_info.jl"), "w") do f
    write(f, cpu_info_string)
end


text:
module Legolas

using Tables, Arrow

"""
    lift(f, x)

Return `f(x)` unless `x isa Union{Nothing,Missing}`, in which case return `missing`.

This is particularly useful when handling values from `Arrow.Table`, whose null values
may present as either `missing` or `nothing` depending on how the table itself was
originally constructed.
"""
lift(::Any, ::Union{Nothing,Missing}) = missing
lift(f, x) = f(x)

"""
    lift(f)

Returns a curried function, `x -> lift(f,x)`
"""
lift(f) = Base.Fix1(lift, f)

const LEGOLAS_SCHEMA_QUALIFIED_METADATA_KEY = "legolas_schema_qualified"

include("rows.jl")
include("tables.jl")

end # module


text:
using GLVisualize, FileIO, ImageIO
robj = visualize(file"nasa.jpg")
println(robj[:preferred_camera])
view(robj)
renderloop()


text:
doc"""
`Atom(Z; R=(0.,0.,0.), params...)`

creates an atom at `R` with nuclear charge `Z`,
with any additional parameters `params`.
"""
type Atom <: Associative{Symbol,Any}
    Z::Float64              # atomic number
    R::Point3f              # position
    d::Dict{Symbol,Any}     # additional properties
end

Base.delete!(at::Atom, key)                                    = delete!(at.d, key)
Base.done(at::Atom, i)                                         = done(at.d, i)
Base.empty!(at::Atom)                                          = empty!(at.d)
Base.copy(at::Atom)                                            = Atom(at.Z,at.R,copy(at.d))
Base.get(at::Atom, key, default)                               = get(at.d, key, default)
Base.get(default::Union{DataType,Function}, at::Atom, key)     = get(default, at.d, key)
Base.get!(at::Atom, key0, default)                             = get!(at.d, key0, default)
Base.get!(default::Union{DataType,Function}, at::Atom, key0)   = get!(default, at.d, key0)
Base.getindex(at::Atom, key)                                   = getindex(at.d, key)
Base.getkey(at::Atom, key, default)                            = getkey(at.d, key, default)
Base.haskey(at::Atom, key)                                     = haskey(at.d, key)
Base.isempty(at::Atom)                                         = isempty(at.d)
Base.length(at::Atom)                                          = length(at.d)
Base.next(at::Atom, i)                                         = next(at.d, i)
Base.pop!(at::Atom, key)                                       = pop!(at.d, key)
Base.pop!(at::Atom, key, default)                              = pop!(at.d, key, default)
Base.setindex!(at::Atom, v0, key0)                             = setindex!(at.d, v0, key0)
Base.similar(at::Atom)                                         = similar(at.d)
Base.sizehint!(at::Atom, newsz)                                = sizehint!(at.d, newsz)
Base.start(at::Atom)                                           = start(at.d)
Base.merge!(at::Atom, d::Dict)                                 = begin merge!(at.d, d); d end
Base.merge(at::Atom, d::Dict)                                  = Atom(merge(at.d, d))
Base.writemime(io::IO, mime::MIME"text/plain", at::Atom)       = show_symbol_dict(io,at)
Base.summary(at::Atom)                                         = string(typeof(at),"(Z=$(at.Z),R=$(at.R))")
Base.show(io::IO, at::Atom)                                    = print(io,summary(at))

function show_symbol_dict(io::IO, d::Associative{Symbol,Any}; cols=100)
    shown_set = get(task_local_storage(), :SHOWNSET, nothing)
    if shown_set === nothing
        shown_set = ObjectIdDict()
        task_local_storage(:SHOWNSET, shown_set)
    end
    d in keys(shown_set) && (print(io, "#= circular reference =#"); return)

    try
        println(io, summary(d))

        shown_set[d] = true
        ks = map(string,keys(d))
        keylen = length(ks) > 0 ? maximum(map(length,ks)) : 0

        for (i,k) in enumerate(keys(d))
            print(io, Base.rpad(ks[i], keylen))
            print(io, " = ")
            val = Base.with_output_limit(()->sprint(show, d[k]))
            val = Base._truncate_at_width_or_chars(val, cols - keylen, "\r\n")
            print(io, val)
            print(io, '\n')
        end
    finally
        delete!(shown_set, d)
    end
end

function Base.call(::Type{Atom};
        Z         ::Real         = 1.,
        R         ::Point3{Real} = (0.,0.,0.),
        makebasis ::Bool         = false,
        makegrid  ::Bool         = false,
        etacut    ::Float64      = 1e-4,
        a         ::Float64      = 10.,
        w         ::Float64      = 2.,
        r0max     ::Float64      = 20.,
        npts      ::Int          = 100,
        nr        ::Int          = 80,
        na        ::Int          = 38,
        params...)

    at = Atom(convert(Float64,Z), convert(Point3f,R), Dict{Symbol,Any}(params))

    if makebasis
        dft!(at)
        splineradial!(at;npts=npts)
        naobasis!(at;etacut=etacut,a=a,w=w,r0max=r0max,npts=npts)
    end

    if makegrid
        makegrid!(at;nr=nr,na=na)
    end

    at
end

typealias Atoms Vector{Atom}

doc"""
`setposition(at,R)`

sets the position of atom `at` to `R`
"""
function setposition(at::Atom, R::Point3f)
    at.R = R
    for k in keys(at)
        setposition(at[k],R)
    end
end

function setposition(a::AbstractArray, R::Point3f)
    for i in 1:length(a)
        setposition(a[i],R)
    end
end

setposition(a, R::Point3f) = return

doc"""
`deepcopy(at,R)`

returns a copy of `at` at the position `R`
"""
function Base.deepcopy(at::Atom, R::Point3f)
    at_ = deepcopy(at)
    setposition(at_,R)
    at_
end

doc"""
`splineradial!(at; npts=100)`

splines the electronic density `rho`, hartree potential `vh`,
and Kohn-Sham potential `vs` and creates radial functions.
Additionally, the Coulomb potential `v` is created.

assumes the `rho`, `vh` and `vs` have been calculated and are
given in `at` discretized on the radial grid `r`.

*sets*

`at[:v]` (see `CoulombPotential`)

`at[:rho]`, `at[:vh]`, `at[:vs]` (see `RadialFunc`).

`npts` - the number of evenly distributed spline points to use

e.g.
```julia
at = dft(Z=1.)
splineradial!(at)
plot(at[:rho].u,0,5,label="rho")
plot!(at[:vh].u,0,5,label="vh")
plot!(at[:vs].u,0,5,label="vs")
ylims!(-2,1)
```
"""
function splineradial!(at::Atom; npts::Int=100)
    r = at[:r]
    Y00 = 1/sqrt(4π)
    q = 4π*multipolemoment(r,at[:rho],0)/Y00
    at[:v]   = CoulombPotential(at.Z,R=at.R)
    at[:rho] = RadialFunc(r,at[:rho],R=at.R,rmax=r[end],npts=npts)
    at[:vh]  = MultipoleFunc(r,at[:vh]./Y00,l=0,m=0,R=at.R,q=q,rcut=r[end],npts=npts)
    at[:vs]  = RadialFunc(r,at[:vs],R=at.R,rmax=r[end],npts=npts)
    return
end

function splineradial(at::Atom; npts::Int=100)
    let at = copy(at)
        splineradial!(at;npts=npts)
        return at
    end
end

function makegrid!(at::Atom; nr::Int=80, na::Int=38, rmax::Float64=Inf)
    if isinf(rmax)
        @assert :basis in keys(at) "create an atomic basis before making a grid"
        rmax = maximum([φ.rmax for φ in at[:basis]])
    end
    r = atomicradgrid(nr,rmax)
    Ω = lebgrid(na)
    at[:g] = atomicgrid(r,Ω,R=at.R)
    return
end

function makegrid(at::Atom; nr::Int=80, na::Int=38, rmax::Float64=Inf)
    let at = copy(at)
        makegrid!(at;nr=nr,na=na,rmax=rmax)
        return at
    end
end

function makegrid!(atoms::Atoms; nr::Int=80, na::Int=38, rmax::Float64=Inf)
    for at in atoms
        makegrid!(at,nr=nr,na=na,rmax=rmax)
    end
end

function makegrid(atoms::Atoms; nr::Int=80, na::Int=38, rmax::Float64=Inf)
    let atoms = copy(atoms)
        makegrid!(atoms;nr=nr,na=na,rmax=rmax)
        return atoms
    end
end

hasbasis(at::Atom) = :basis in keys(at)
hasgrid(at::Atom) = :g in keys(at)

print_stats(at::Atom) = print_stats(at[:stats])

max_stars_repo_name: johncsnyder/Element


text:
module StochasticGene

using Distributions
using StatsBase
using Statistics
using DelimitedFiles
using Dates
using Distributed
using LinearAlgebra
using Plots
using SparseArrays
using DifferentialEquations
using LSODA
using DataFrames
using FFTW
using Downloads
using CSV
using MultivariateStats

export
    rna_setup,
    makeswarm,
    fit_rna,
    make_dataframes,
    write_dataframes,
    write_dataframes_only,
    write_winners,
    write_augmented,
    fix_filenames,
    fix,
    large_deviance,
    large_rhat,
    assemble_all,
    write_histograms,
    plot_histogram,
    simulatorGM

### Source files

# Type system
include("common.jl")

# commonly used functions
include("utilities.jl")

# Metropolis Hastings MCMC for computing posterior distributions of model parameters
include("metropolis_hastings.jl")

# Input output functions
include("io.jl")

# Chemical master equation solutions of stochastic models for likelihood functions in fitting algorithms
include("chemical_master.jl")

# Transition rate matrices of stochastic models defining master equations
include("transition_rate_matrices.jl")

# Probability distributions by direct simulation of stochastic models using Gillespie and Gibson-Bruck algorithms
include("simulator.jl")
include("telegraphsplice.jl")

# functions specific for Gene Trap experiments of Wan et al.
include("genetrap.jl")

# functions for scRNA and FISH experiments
include("rna.jl")

# functions for post fit analysis and plotting
include("analysis.jl")

# functions for use on NIH cluster Biowulf
include("biowulf.jl")

# include("/Users/carsonc/Dropbox/Larson/GeneTrap_analysis/code/GillespieSimulatorExInt.jl")
# include("scRNA.jl")
# include("tff1.jl")

"""
A Julia module for simulation and Bayesian inference of parameters of stochastic models of gene transcription.

API Overview:

"""

end #Module

max_stars_repo_name: nih-niddk-mbs/Transcription

max_issues_repo_name: nih-niddk-mbs/Transcription

max_forks_repo_name: nih-niddk-mbs/Transcription


text:
import StatsBase: RealVector, RealMatrix
using KernelEstimator, Optim
import KernelEstimator.lcv, KernelEstimator.bwlcv, KernelEstimator.midrange

function dirichletkernel_oneval(x::RealVector, xdata::RealVector, λ::Real; log = false)
    if(log)
        logpdf(Dirichlet(1 .+ xdata/λ), x)
    else
        pdf(Dirichlet(1 .+ xdata/λ), x)
    end
end

function dirichletkernel(x::RealVector, xdata::RealMatrix, λ::Real, w::Vector, n::Int; log = false)
    # w .= mapslices(xdatavec -> dirichletkernel_oneval(x, xdatavec, λ), xdata, 2) |> vec
    for i in eachindex(w)
        w[i] = dirichletkernel_oneval(x, xdata[i,:], λ; log = log)
    end
    nothing
end

function dirichletkernel(x::RealVector, xdata::Array{Array{Float64,1},1}, λ::Real, w::Vector, n::Int; log = false)
    # w .= mapslices(xdatavec -> dirichletkernel_oneval(x, xdatavec, λ), xdata, 2) |> vec
    for i in eachindex(w)
        w[i] = dirichletkernel_oneval(x, xdata[i], λ; log = log)
    end
    nothing
end

function midrange(x::RealMatrix)
    mapslices(xvec -> quantile(xvec, [.25, .75]), x, dims = 1) |> x -> x[2,:] - x[1,:] |> maximum
#     lq, uq = quantile(x, [.25, .75])
#     uq - lq
end

# midrange(x::Array{Array{Float64,1},1}) = midrange.(x) |> maximum
# midrange(x::Array{Array{Int64,1},1}) = midrange.(x) |> maximum

ff(x)::Float64 = max(midrange([xi[1] for xi in x]), midrange([xi[2] for xi in x]), midrange([xi[3] for xi in x]))

midrange(x::Array{Array{Float64,1},1})::Float64 = ff(x)

midrange(x::Array{Array{Int64,1},1})::Float64 = ff(x)

function lcv(xdata::RealMatrix, kernel::Function, h::Real, w::Vector, n::Int)
#     -mean(kerneldensity(xdata,xdata,kernel,h)) + mean(map(kernel, xdata, xdata, h))
    ind = 1
    ind_end = 1+n
    ll = 0.0
    @inbounds while ind < ind_end
        # kernel(xdata[ind,:], xdata, h, w, n)
        kernel(view(xdata, ind, :), xdata, h, w, n)
        w[ind] = 0.0
        ll += log(mean(w))
        ind += 1
    end
    -ll
end

function lcv(xdata::Array{Array{Float64,1},1}, kernel::Function, h::Real, w::Vector, n::Int)
    #     -mean(kerneldensity(xdata,xdata,kernel,h)) + mean(map(kernel, xdata, xdata, h))
    ind = 1
    ind_end = 1+n
    ll = 0.0
    @inbounds while ind < ind_end
        kernel(xdata[ind], xdata, h, w, n)
        w[ind] = 0.0
        ll += log(mean(w))
        ind += 1
    end
    -ll
end

function bwlcv(xdata::RealMatrix, kernel::Function)
    #This function seems to hit against the higher bound
    n = size(xdata,1)
    w = zeros(n)
    h0=midrange(xdata)
    if h0==0 #Algorithm returns 0 when h0=0. This happens when there are many ties, mixdrange can return 0
        h0 = median(xdata)
    end
    hlb = h0/n^2
    hub = h0
    if kernel==betakernel
        hub = 0.25
    end
    return Optim.minimizer(Optim.optimize(h->lcv(xdata,kernel,h,w,n), hlb, hub, iterations=200,abs_tol=h0/n^2))
end

function bwlcv(xdata::Array{Array{Float64,1},1}, kernel::Function)
    n = length(xdata)
    w = zeros(n)
    h0=midrange(xdata)
    if h0==0 #Algorithm returns 0 when h0=0. This happens when there are many ties, mixdrange can return 0
        h0 = median(xdata |> Base.Iterators.flatten)
    end
    hlb = h0/n^2
    hub = h0
    if kernel==betakernel
        hub = 0.25
    end
    return Optim.minimizer(Optim.optimize(h->lcv(xdata,kernel,h,w,n), hlb, hub, iterations=200,abs_tol=h0/n^2))
end

function bwlcv_large_bounds(xdata::RealMatrix, kernel::Function)
    n = size(xdata,1)
    w = zeros(n)
    h0=midrange(xdata)
    if h0==0 #Algorithm returns 0 when h0=0. This happens when there are many ties, mixdrange can return 0
        h0 = median(xdata)
    end
    hlb = h0/n^2
    hub = h0*10
    if kernel==betakernel
        hub = 0.25
    end
    res = Optim.minimizer(Optim.optimize(h->lcv(xdata,kernel,h,w,n), hlb, hub, iterations=200,abs_tol=h0/n^2))
    i = 1
    while (res==hub)&(i<10)
        res = Optim.minimizer(Optim.optimize(h->lcv(xdata,kernel,h,w,n), hlb, 10*i*hub, iterations=200,abs_tol=h0/n^2))
        i = i+1
    end
    return res
end

function minus_log_leaveoneout(xdata::RealMatrix, kernel::Function, h::Real, w::Vector, n::Int)

    ll = 0.0
    @inbounds for ind in 1:n
        kernel(view(xdata, ind, :), xdata, h, w, n; log = false)
        w[ind] = 0.0
        ll += log(sum(w))
    end
    - ll + log(n-1)
end

function bwloo_large_bounds(xdata::RealMatrix, kernel::Function)
    n = size(xdata,1)
    w = zeros(n)
    h0=midrange(xdata)
    if h0==0 #Algorithm returns 0 when h0=0. This happens when there are many ties, mixdrange can return 0
        h0 = median(xdata)
    end
    hlb = h0/n^2
    hub = h0#*10
    if kernel==betakernel
        hub = 0.25
    end
#     println(hlb)
#     println(hub)
#     res = Optim.minimizer(Optim.optimize(h->lcv(xdata,kernel,h,w,n), hlb, hub, iterations=200,abs_tol=h0/n^2))
    res = Optim.minimizer(Optim.optimize(h->minus_log_leaveoneout(xdata,kernel,h,w,n), hlb, hub, iterations=200))
    i = 1
    while (isapprox(res,hub, atol = h0/n^2))&(i<10)
        res = Optim.minimizer(Optim.optimize(h->minus_log_leaveoneout(xdata,kernel,h,w,n), hlb, 10^i*hub, iterations=200,abs_tol=h0/n^2))
        i = i+1
    end
    return res
end

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
#=
  Test set 2: Construct a spine from all 3 tables and link all records.
=#

println("\nTEST SET 2")

println("Link all health service events")
outdir2 = run_linkage(joinpath("config", "link_all_health_service_events.toml"))

println("Reporting")
outfile = joinpath(outdir, "linkage_report.csv")
summarise_linkage_run(outdir2, outfile)
result  = table_to_set_of_dicts(outfile)

@test length(result) == 9
@test in(Dict(:tablename => "LINKAGE RUNS",            :status => outdir2,                     :nrecords => missing), result)
@test in(Dict(:tablename => "spine",                   :status => "existent",                  :nrecords => 6), result)
@test in(Dict(:tablename => "hospital_admissions",     :status => "linked with criteria ID 1", :nrecords => 5), result)
@test in(Dict(:tablename => "emergency_presentations", :status => "linked with criteria ID 2", :nrecords => 4), result)
@test in(Dict(:tablename => "emergency_presentations", :status => "linked with criteria ID 3", :nrecords => 1), result)
@test in(Dict(:tablename => "influenza_cases",         :status => "linked with criteria ID 4", :nrecords => 3), result)
@test in(Dict(:tablename => "influenza_cases",         :status => "linked with criteria ID 5", :nrecords => 2), result)
@test in(Dict(:tablename => "influenza_cases",         :status => "linked with criteria ID 6", :nrecords => 1), result)
@test in(Dict(:tablename => "influenza_cases",         :status => "linked with criteria ID 7", :nrecords => 2), result)


text:
"""
    Projector(ϕs, orbitals)

Represents the projector on the subspace spanned by the radial
orbitals `ϕs` (corresponding to `orbitals`).
"""
struct Projector{T,B<:AbstractQuasiMatrix,RO<:RadialOrbital{T,B},O} <: NBodyOperator{1}
    ϕs::Vector{RO}
    orbitals::Vector{O}
end

Base.iszero(me::OrbitalMatrixElement{1,A,<:Projector,B}) where {A<:SpinOrbital,B<:SpinOrbital} =
    me.a[1] ∉ me.o.orbitals || me.b[1] ∉ me.o.orbitals

function Base.show(io::IO, projector::Projector)
    write(io, "P(")
    write(io, join(string.(projector.orbitals), " "))
    write(io, ")")
end

projectout!(y::RO, ::Nothing) where RO = y

"""
    projectout!(y, projector)

Project out all components of `y` parallel to the radial orbitals
`projector.ϕs`.
"""
function projectout!(y::RO, projector::Proj) where {RO,Proj<:Projector}
    yc = y.args[2]

    for ϕ in projector.ϕs
        c = materialize(applied(*, ϕ', y))
        yc .-= c*ϕ.args[2]
        # y -= c*ϕ
    end
    y
end


text:
######################################################################
# points.jl: projective points on elliptic curves
######################################################################

export EllipticPoint, Point

export coordinates, base_curve, isinfinity, normalize!, normalized, samefields

"""
Concrete type for projective points on elliptic curves given by a projective planar equation.
"""

mutable struct EllipticPoint{T<:Nemo.RingElem} <: ProjectivePoint{T}
	X::T
	Y::T
	Z::T
	curve::EllipticCurve{T}
end

coordinates(P::EllipticPoint) = (P.X, P.Y, P.Z)

base_curve(P::EllipticPoint) = P.curve

######################################################################
# Calling with the 'Point' constructor
######################################################################

function Point(x::T, y::T, curve::EllipticCurve{T}) where T
	R = base_ring(x)
	return EllipticPoint(x, y, R(1), curve)
end

function Point(X::T, Y::T, Z::T, curve::EllipticCurve{T}) where T
	return EllipticPoint(X, Y, Z, curve)
end

######################################################################
# Basic methods for projective points
######################################################################

base_ring(P::EllipticPoint) = parent(P.X)

"""
Decide whether a given elliptic point is valid.
"""
function isvalid(P::EllipticPoint)
	Eq = projective_equation(P.curve)
	x, y, z = coordinates(P)
	return Eq(x, y, z) == 0 && (x, y, z) != (0, 0, 0)
end

"""
Decide whether two projective points are given by the exact same coordinates.
"""
samefields(P::EllipticPoint, Q::EllipticPoint) = (P.curve == Q.curve) & (P.X == Q.X) & (P.Y == Q.Y) & (P.Z == Q.Z)

"""
Describe a projective point giving its X, Y and Z coordinates, and the curve it lives on.
"""
show(io::IO, P::EllipticPoint) = print("Point (", P.X, ":", P.Y, ":", P.Z, ") on ", P.curve)

"""
Decide whether a projective point is at infinity.
"""
isinfinity(P::EllipticPoint) = iszero(P.Z)

"""
Get a new normalized projective point from any projective point.
"""
function normalized(P::EllipticPoint)
    K = base_ring(P)
    if isinfinity(P)
        return EllipticPoint(zero(K),
                               one(K),
                               zero(K),
                               P.curve)
    else
        return EllipticPoint(P.X // P.Z,
                               P.Y // P.Z,
                               one(K),
                               P.curve)
    end
end

"""
Normalize a projective point in place.
"""
function normalize!(P::EllipticPoint{T}) where T<:Nemo.FieldElem
    K = base_ring(P)
    if isinfinity(P)
        P.X = zero(K)
        P.Y = one(K)
        P.Z = zero(K)
    else
        P.X = P.X // P.Z
        P.Y = P.Y // P.Z
        P.Z = one(K)
    end
    return
end

"""
Decides whether two projective points are the same.
"""

==(P::EllipticPoint, Q::EllipticPoint) = samefields(normalized(P), normalized(Q))


text:
using Test
using SumOfSquares
using DynamicPolynomials

# Like `term_test` but with a variable `y` fixed to 1.
function term_fixed_test(
    optimizer, config::MOIT.Config, cone::SumOfSquares.PolyJuMP.PolynomialSet)
    atol = config.atol
    rtol = config.rtol

    model = _model(optimizer)

    @variable(model, α)

    @polyvar x y
    set = @set y == 1
    cref = @constraint(model, (α - 1) * x^2 * y in cone, domain = set)

    @objective(model, Min, α)
    optimize!(model)

    @test termination_status(model) == MOI.OPTIMAL
    @test objective_value(model) ≈ 1.0 atol=atol rtol=rtol

    @test primal_status(model) == MOI.FEASIBLE_POINT
    @test value(α) ≈ 1.0 atol=atol rtol=rtol

    test_constraint_primal(cref, 0.0)

    p = gram_matrix(cref)
    @test getmat(p) ≈ zeros(1, 1) atol=atol rtol=rtol
    @test p.basis.monomials == [x]

    @test dual_status(model) == MOI.FEASIBLE_POINT
    for (m, μ) in [(x^2 * y, dual(cref)), (x^2, moments(cref))]
        @test μ isa AbstractMeasure{Float64}
        @test length(moments(μ)) == 1
        @test moment_value(moments(μ)[1]) ≈ 1.0 atol=atol rtol=rtol
        @test monomial(moments(μ)[1]) == m
    end

    ν = moment_matrix(cref)
    @test getmat(ν) ≈ ones(1, 1) atol=atol rtol=rtol
    @test ν.basis.monomials == [x]

    S = SumOfSquares.SOSPolynomialSet{
        typeof(set), Monomial{true}, MonomialVector{true}, SumOfSquares.Certificate.Remainder{SumOfSquares.Certificate.Newton{typeof(cone), MonomialBasis, Tuple{}}}
    }
    @test list_of_constraint_types(model) == [(Vector{JuMP.AffExpr}, S)]
    test_delete_bridge(
        model, cref, 1,
        ((MOI.VectorOfVariables, MOI.Nonnegatives, 0),
         (MOI.VectorAffineFunction{Float64},
          SumOfSquares.PolyJuMP.ZeroPolynomialSet{
              typeof(set), MonomialBasis,
              Monomial{true}, MonomialVector{true}},
          0)))
end
sos_term_fixed_test(optimizer, config)   = term_fixed_test(optimizer, config, SOSCone())
sd_tests["sos_term_fixed"] = sos_term_fixed_test
sdsos_term_fixed_test(optimizer, config) = term_fixed_test(optimizer, config, SDSOSCone())
soc_tests["sdsos_term_fixed"] = sdsos_term_fixed_test
dsos_term_fixed_test(optimizer, config)  = term_fixed_test(optimizer, config, DSOSCone())
linear_tests["dsos_term_fixed"] = dsos_term_fixed_test

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
module PhysicsCodesLabXD

# The structure of the package is the following:

# 1. Define tensor category objects: simple objects, dual, braiding structure, tensor product

end


text:
module SequencerJ

using Reexport
@reexport using Distances
@reexport using Logging
using Printf: @sprintf
using SparseArrays: sparse, spzeros
using LightGraphs
using SimpleWeightedGraphs
using StatsBase: mean, ecdf, ECDF, sample
using LinearAlgebra: issymmetric, Symmetric, I

import Distances: Euclidean, KLDivergence, PeriodicEuclidean
import LightGraphs: bfs_parents, _bfs_parents
import Base: show

include("bfs.jl")
include("utils.jl")
include("distancemetrics.jl")
include("sequencer.jl")

# functions
export sequence, autoscale, D, mst, elong, elongation, order, prettyp, loss, ensuregrid!
export emd, energy

# types
export SequencerResult, EMD, Energy

# constants
export L2, WASS1D, KLD, ENERGY, ALL_METRICS

end


text:
# Standard deps.
using InstantiateFromURL, Test, Pkg

# TOML deps.
import Pkg.TOML: parsefile

@testset "Main Tests" begin include("maintest.jl") end
@testset "TOML Tests" begin include("tomltest.jl") end
@testset "Directory Names Tests" begin include("namestest.jl") end # same as maintest, but in a different directory
@testset "Test activate_github_path" begin include("pathtest.jl") end


text:
using Test
using DataFrames
using CSV
using ..ChemistryFeaturization.FeatureDescriptor

const cf = ChemistryFeaturization

@testset "ElementFeatureDescriptor" begin
    # test for errors
    @test_throws AssertionError ElementFeatureDescriptor("heffalump")

    # construct a few from built-in features...
    local fnames = ["X", "Block", "Atomic mass"]
    local lengths = [default_nbins, 4, default_nbins]
    local cats = [false, true, false]
    local conts = [false, false, false]

    for i = 1:3
        f = ElementFeatureDescriptor(fnames[i])
        @test output_shape(f) == lengths[i]
        @test f.categorical == cats[i]
    end

    # # encode/decode some stuff...
    He_mol = AtomGraph(Float32.([1 0; 0 1]), ["He", "He"])
    triangle_C = AtomGraph(Float32.([0 1 1; 1 0 1; 1 1 0]), ["C", "C", "C"])
    X, block, amass = ElementFeatureDescriptor.(fnames)

    @testset "Encode" begin
        @test_throws AssertionError X(He_mol)

        @test block(He_mol) == Float64.([0 0; 1 1; 0 0; 0 0])
        @test amass(He_mol)[3, :] == ones(2)

        # now let's try some options (and make sure they're ignored when appropriate...)...
        X = ElementFeatureDescriptor("X", nbins = 8, logspaced = true)
        @test X(triangle_C)[6, :] == ones(3)

        block = ElementFeatureDescriptor("Block", nbins = 5)
        @test block(triangle_C)[2, :] == ones(3)

        @testset "Explicitly specify Codec" begin
            ohoc_codec = amass.encoder_decoder
            amass_v2 = ElementFeatureDescriptor("Atomic mass", ohoc_codec)
            @test amass_v2(He_mol) == amass(He_mol)
        end
    end

    @testset "Decode" begin
        @test decode(block, block(He_mol)) == ["p", "p"]

        true_He_amass = atom_data_df[2, Symbol("Atomic mass")]
        He_amass_min, He_amass_max = decode(amass, amass(He_mol)[:, 1])
        @test He_amass_min < true_He_amass < He_amass_max

        true_X = atom_data_df.X[6]
        decoded_X_range = decode(X, X(triangle_C)[:, 1])
        @test decoded_X_range[1] < true_X < decoded_X_range[2]
        @test decode(block, block(triangle_C)) == ["p", "p", "p"]
    end

    # and make a custom lookup table...
    @testset "Custom Lookup Table" begin
        df = CSV.read(abspath(@__DIR__, "..", "test_data", "lookup_table.csv"), DataFrame)
        meaning = ElementFeatureDescriptor("MeaningOfLife", df)
        @test meaning(triangle_C)[10, :] == ones(3)

        @testset "Encodable Elements" begin
            @test encodable_elements(meaning) == ["C", "As", "Tc"]
            @test encodable_elements("MeaningOfLife", df) == ["C", "As", "Tc"]
        end
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    mean(solution)

Mean of simulation `solution`.
"""
function mean(solution::SimulationSolution)
  data = DataFrame([variable => mean(reals) for (variable, reals) in solution.realizations])
  georef(data, solution.domain)
end

"""
    var(solution)

Variance of simulation `solution`.
"""
function var(solution::SimulationSolution)
  data = DataFrame([variable => var(reals) for (variable, reals) in solution.realizations])
  georef(data, solution.domain)
end

"""
    quantile(solution, p)

p-quantile of simulation `solution`.
"""
function quantile(solution::SimulationSolution, p::Number)
  cols = []
  for (variable, reals) in solution.realizations
    quantiles = map(1:nelms(solution.domain)) do location
      slice = getindex.(reals, location)
      quantile(slice, p)
    end
    push!(cols, variable => quantiles)
  end
  georef(DataFrame(cols), solution.domain)
end

quantile(solution::SimulationSolution, ps::AbstractVector) = [quantile(solution, p) for p in ps]


text:
# Copyright (c) 2019 Guido Kraemer
# Copyright (c) 2019 Simon Castano
# Distributed under the MIT software license, see the accompanying
# file COPYING or http://www.opensource.org/licenses/mit-license.php.

"""
    Witness

Witness data type has a single `data` field, stored as `Vector{Vector{UInt8}}`
in which each inner `Vector{UInt8}` represent a stack item.

The `Witness` is a serialization of all witness data of the transaction.
Each `TxIn` is associated with a witness field. A witness field starts with a
`CompactSizeUInt` to indicate the number of stack items for the `TxIn`.
It is followed by stack items, with each item starts with a `CompactSizeUInt`
to indicate the length.
Witness data is NOT script.

A non-witness program `TxIn` MUST be associated with an empty witness field,
represented by a `0x00`. If all `TxIn`s are not witness program, a
transaction's `wtxid` is equal to its `txid`.
"""
struct Witness
    data :: Vector{Vector{UInt8}}
end

Witness() = Witness([[0x00]])

"""
    Witness(io::IO) -> Witness

Parse `Witness` from an IOBuffer
"""
function Witness(io::IO)
    n_items = CompactSizeUInt(io).value
    data = Vector{Vector{UInt8}}(undef, n_items)

    for i in 1:n_items
        l = CompactSizeUInt(io).value
        data[i] = read!(io, Array{UInt8, 1}(undef, l))
    end

    Witness(data)
end

function Base.show(io::IO, program::Witness)
    println("Witness Program:")
    for item in program.data
        println(io, bytes2hex(item))
    end
end

"""
    serialize(program::Witness) -> Vector{UInt8}

Serialize a `Witness` data type starting with a `CompactSizeUInt` to indicate
the number of stack items for the `TxIn`.
It is followed by stack items, with each item starts with a `CompactSizeUInt`
to indicate the length.
"""
function serialize(program::Witness)
    l = CompactSizeUInt(length(program.data))
    result = serialize(l)
    for item in program.data
        l = CompactSizeUInt(length(item))
        append!(result, serialize(l))
        append!(result, item)
    end
    return result
end