text:
# PatMat hinge
function f(μ::Real, α0::AbstractVector, β0::AbstractVector, δ0::Real, C1::Real, C2::Real)
    return C2*δ0 + C2^2*sum(max.(β0 .- μ, 0)) - μ
end

function h(μ::Real, α0::AbstractVector, β0::AbstractVector, δ0::Real, C1::Real, C2::Real)
    λ = f(μ, α0, β0, δ0, C1, C2)
    return sum(min.(max.(α0 .- λ, 0), C1)) - sum(min.(max.(β0 .+ λ, 0), λ + μ))
end

function projection(α0::AbstractVector, β0::AbstractVector, δ0::Real, C1::Real, C2::Real)

    if δ0 <= - C2*sum(max.(β0 .+ maximum(α0), 0))
        return zero(α0), zero(β0), zero(δ0)
    end

    μ_ub1 = minimum(β0)
    μ_ub2 = minimum(C2*δ0 + C2^2*sum(β0) .- α0)/(C2^2*length(β0) + 1)
    μ_ub  = min(μ_ub1, μ_ub2) - 1e8
    μ_lb  = maximum(β0) + C2*max(δ0, 0) + 1e8

    μ = find_root(μ -> h(μ, α0, β0, δ0, C1, C2), (μ_lb, μ_ub))
    λ = f(μ, α0, β0, δ0, C1, C2)

    α = @. min(max(α0 - λ, 0), C1)
    β = @. min(max(β0 + λ, 0), λ + μ)
    δ = (λ + μ)/C2

    return α, β, δ
end

# PatMat quadratic
function f(λ::Real, α0::AbstractVector, β0::AbstractVector, δ0::Real)
    return sum(max.(α0 .- λ, 0)) - sum(max.(β0 .+ λ, 0))
end

function projection(α0::AbstractVector, β0::AbstractVector, δ0::Real)

    if - maximum(β0) > maximum(α0)
        return zero(α0), zero(β0), zero(δ0)
    end

    λ = find_root(λ -> f(λ, α0, β0, δ0), (- maximum(β0), maximum(α0)))

    α = @. max(α0 - λ, 0)
    β = @. max(β0 + λ, 0)
    δ = max(δ0, 1e-4)

    return α, β, δ
end

# TopPushK hinge
function f(μ::Real, s::AbstractVector, K::Integer)
    i, j, d  = 2, 1, 1
    λ, λ_old = s[1], 0
    g, g_old = -K*μ, -K*μ

    while g < 0
        g_old = g
        λ_old = λ

        if s[i] <= s[j] + μ
            g += d*(s[i] - λ)
            λ  = s[i]
            d += 1
            i += 1
        else
            g += d*(s[j] + μ - λ)
            λ  = s[j] + μ
            d -= 1
            j += 1
        end
    end
    return -(λ - λ_old)/(g - g_old)*g_old + λ_old
end

function h(μ::Real, s::AbstractVector, α0::AbstractVector, β0::AbstractVector, C::Real, K::Integer)
    λ = f(μ, s, K)
    return sum(min.(max.(α0 .- λ .+ sum(max.(β0 .+ λ .- μ, 0))/K, 0), C)) - K*μ
end

function projection(α0::AbstractVector, β0::AbstractVector, C::Real, K::Real)

    if K >= length(β0)
        @error "No feasible solution: K < length(β0) needed."
        return α0, β0
    end

    if mean(partialsort(β0, 1:K; rev = true)) + maximum(α0) <= 0
        return zero(α0), zero(β0)
    end

    s    = vcat(.- sort(β0; rev = true), Inf)
    μ_lb = 1e-10
    μ_ub = length(α0)*C/K + 1e-6

    μ = find_root(μ -> h(μ, s, α0, β0, C, K), (μ_lb, μ_ub))
    λ = f(μ, s, K)
    δ = sum(max.(β0 .+ λ .- μ, 0))/K

    α = @. max(min(α0 - λ + δ, C), 0)
    β = @. max(min(β0 + λ, μ), 0)

    return α, β
end

# TopPushK quadratic
function h(μ::Real, s::AbstractVector, α0::AbstractVector, β0::AbstractVector, K::Integer)
    λ = f(μ, s, K)
    return sum(max.(α0 .- λ .+ sum(max.(β0 .+ λ .- μ, 0))/K, 0)) - K*μ
end

function projection(α0::AbstractVector, β0::AbstractVector, K::Integer)

    if K >= length(β0)
        @error "No feasible solution: K < length(β0) needed."
        return α0, β0
    end

    if mean(partialsort(β0, 1:K; rev = true)) + maximum(α0) <= 0
        return zero(α0), zero(β0)
    end

    s    = vcat(.- sort(β0; rev = true), Inf)
    μ_lb = 1e-10
    μ_ub = length(α0)*(maximum(α0) + maximum(β0))/K

    μ = find_root(μ -> h(μ, s, α0, β0, K), (μ_lb, μ_ub))
    λ = f(μ, s, K)
    δ = sum(max.(β0 .+ λ .- μ, 0))/K

    α = @. max(α0 - λ + δ, 0)
    β = @. max(min(β0 + λ, μ), 0)

    return α, β
end

function find_root(f, lims)
    try
        Roots.find_zero(f, lims)
    catch
        Roots.fzero(f, sum(lims)/2)
    end
end


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

include("delays.jl")

using Actors, Test, .Threads, .Delays
import Actors: spawn, info, diag, newLink

Base.:(==)(l1::Link, l2::Link) = hash(l1) == hash(l2)

t1 = Ref{Task}()
t2 = Ref{Task}()
t3 = Ref{Task}()

println("Testing error handling with some failures:")

#
# test peer connections
#
act1 = spawn(connect, taskref=t1)
a1   = diag(act1, :act)
act2 = spawn(connect, taskref=t2)
a2   = diag(act2, :act)
send(act1, act2)
@test @delayed a1.conn[1] isa Actors.Peer
@test a1.conn[1].lk == act2
@test @delayed a2.conn[1] isa Actors.Peer
@test a2.conn[1].lk == act1
@test t1[].state == :runnable
@test t2[].state == :runnable
send(act1, "boom")
@test @delayed t1[].state == :failed
@test t2[].state == :done

act1 = spawn(connect, taskref=t1)
act2 = spawn(connect, taskref=t2)
act3 = spawn(connect, taskref=t3)
send(act1, act2)
send(act3, act2)
sleep(0.5)
send(act1, "boom")
@test @delayed t1[].state == :failed
@test @delayed t2[].state == :done
@test @delayed t3[].state == :done

act1 = spawn(connect, taskref=t1)
act2 = spawn(connect, taskref=t2)
act3 = spawn(connect, taskref=t3)
send(act1, act2)
send(act3, act2)
trapExit(act3)
@test @delayed act3.mode == :sticky
send(act1, "boom")
@test @delayed t1[].state == :failed
@test t2[].state == :done
@test t3[].state == :runnable
@test length(diag(act3, :err)) == 1
@test diag(act3, :err)[1].exception isa MethodError

# disconnect
act1 = spawn(connect, taskref=t1)
act2 = spawn(connect, taskref=t2)
act3 = spawn(connect, taskref=t3)
send(act1, act2)
send(act3, act2)
a1   = diag(act1, :act)
a2   = diag(act2, :act)
a3   = diag(act3, :act)
@test @delayed a1.conn[1].lk == act2
@test @delayed a2.conn[1].lk in (act1, act3)
@test @delayed a2.conn[2].lk in (act1, act3)
@test @delayed a3.conn[1].lk == act2
become!(act1, disconnect)
send(act1, act2)
become!(act2, disconnect)
send(act2, act3)
@test @delayed isempty(a1.conn)
@test @delayed isempty(a2.conn)
@test @delayed isempty(a3.conn)

# connect and disconnect to _ROOT
act1 = spawn(threadid, taskref=t1)
connect(act1)
a1   = diag(act1, :act)
at1  = diag(act1, :task)
rt   = diag(Actors._ROOT, :act)
@test a1.conn[1].lk == Actors._ROOT
@test rt.conn[1].lk == act1
send(act1, "boom")
@test @delayed t1[].state == :failed
@test Actors.diag(Actors._ROOT, :task).state == :runnable
@test isempty(rt.conn)
errors = exec(Actors._ROOT, Actors.errored)
@test errors[end] == at1
@test errors[end].exception isa MethodError
act1 = spawn(threadid, taskref=t1)
connect(act1)
a1   = diag(act1, :act)
@test @delayed a1.conn[1].lk == Actors._ROOT
@test rt.conn[1].lk == act1
disconnect(act1)
@test @delayed isempty(a1.conn)
@test isempty(rt.conn)

#
# test monitors
#
me = newLink()
act1 = spawn(threadid, taskref=t1)
act2 = spawn(monitor, act1, taskref=t2)
send(act2, send, me)
a1 = diag(act1, :act)
a2 = diag(act2, :act)
@test a1.conn[1] isa Actors.Monitor
@test a1.conn[1].lk == act2
@test a2.conn[1] isa Actors.Monitored
@test a2.conn[1].lk == act1
@test a2.conn[1].action.f == send
send(act1, "boom")
f1 = receive(me)
@test f1 isa MethodError
@test Actors.diag(act2, :task).state == :runnable
@test isempty(a2.conn)
act1 = spawn(threadid, taskref=t1)
become!(act2, monitor, act1)
send(act2)
a1 = diag(act1, :act)
@test @delayed a1.conn[1] isa Actors.Monitor
@test a1.conn[1].lk == act2
@test isempty(a2.conn)
send(act1, "boom")
@test @delayed t1[].state == :failed
@test isempty(me.chn)
act1 = spawn(threadid, taskref=t1)
become!(act2, monitor, act1)
send(act2, send, me)
a1 = diag(act1, :act)
@test @delayed a1.conn[1].lk == act2
@test a2.conn[1].lk == act1
become!(act2, demonitor)
send(act2, act1)
@test @delayed isempty(a1.conn)
@test isempty(a2.conn)

# monitor with _ROOT
act1 = spawn(threadid, taskref=t1)
monitor(act1, send, me)
a1 = diag(act1, :act)
rt = diag(Actors._ROOT, :act)
@test @delayed a1.conn[1] isa Actors.Monitor
@test a1.conn[1].lk == Actors._ROOT
@test rt.conn[1] isa Actors.Monitored
@test rt.conn[1].lk == act1
@test rt.conn[1].action.f == send
send(act1, "boom")
f1 = receive(me)
@test f1 isa MethodError
@test Actors.diag(Actors._ROOT, :task).state == :runnable
@test isempty(rt.conn)
act1 = spawn(threadid, taskref=t1)
monitor(act1)
a1 = diag(act1, :act)
@test @delayed a1.conn[1] isa Actors.Monitor
@test a1.conn[1].lk == Actors._ROOT
@test isempty(rt.conn)
send(act1, "boom")
@test @delayed isempty(me.chn)
act1 = spawn(threadid, taskref=t1)
monitor(act1, send, me)
a1 = diag(act1, :act)
rt = diag(Actors._ROOT, :act)
@test a1.conn[1].lk == Actors._ROOT
@test rt.conn[1].lk == act1
demonitor(act1)
a1 = diag(act1, :act)
rt = diag(Actors._ROOT, :act)
@test isempty(a1.conn)
@test isempty(rt.conn)

# monitor a task
struct Stop end
me = newLink()
tc = Channel(1)
function ttask(ch)
    while true
        msg = take!(ch)
        msg isa Stop ? break : error(msg)
    end
end
tt = Threads.@spawn ttask(tc)
act1 = spawn(Bhv(monitor, tt, timeout=2), taskref=t1)
mt = call(act1, send, me)
@test mt[].state == :runnable
a1 = diag(act1, :act)
@test @delayed a1.conn[1] isa Actors.Monitored
@test a1.conn[1].lk == tt
@test receive(me) == :timed_out
@test isempty(a1.conn)
mt = call(act1, send, me)
@test mt[].state == :runnable
put!(tc, Stop())
@test receive(me) == :normal
@test mt[].state == :done
@test tt.state == :done
tt = Threads.@spawn ttask(tc)
become!(act1, Bhv(monitor, tt, timeout=2))
mt = call(act1, send, me)
@test mt[].state == :runnable
put!(tc, "boom")
msg = receive(me)
@test msg isa ErrorException
@test msg.msg == "boom"
@test tt.state == :failed
@test mt[].state == :done
@test isempty(a1.conn)

ext: jl


text:
#___  ____       ___ _____   _   _       _ _
#|  \/  | |     |_  /  __ \ | | | |     (_| |
#| .  . | |       | | /  \/ | | | |_ __  _| |_ ___
#| |\/| | |       | | |     | | | | '_ \| | __/ _ \
#| |  | | |___/\__/ | \__/\ | |_| | | | | | || (_) |
#_______\_____\____/ \____/  _____|_|___|_____\_____ _____ _____
#| ___ \        (_)         | | \ \ / / / __  |  _  / __  |  _  |
#| |_/ _ __ ___  _  ___  ___| |_ \ V /  `' / /| |/' `' / /| |/' |
#|  __| '__/ _ \| |/ _ \/ __| __|/   \    / / |  /| | / / |  /| |
#| |  | | | (_) | |  __| (__| |_/ /^\ \ ./ /__\ |_/ ./ /__\ |_/ /
#\_|  |_|  \___/| |\___|\___|\__\/   \/ \_____/\___/\_____/\___/
#              _/ |
#             |__/
#
# This code is part of the proposal of the team "MLJC UniTo" - University of Turin
# for "ProjectX 2020" Climate Change for AI.
# The code is licensed under MIT 3.0
# Please read readme or comments for credits and further information.

# Compiler: Julia 1.5.3
#= Short description of this file: LOTKA-VOLTERRA EQUATIONS.

Numerical approach using DifferentialEquations.jl and Neural ODE layer approach
using DiffEqFlux.jl based on Chris Rackauckas' paper " DiffEqFlux.jl -
A Julia Library for Neural Differential Equations"
(https://arxiv.org/abs/1902.02376 )
=#

using Flux, DiffEqFlux, DifferentialEquations, Plots

#= To solve the equations numerically, we first define a problem type by giving
the equation, the initial condition, and timespan to solve over.
=#
function lotka_volterra(du,u,p,t) #Define equations
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = -δ*y + γ*x*y
end
u0 = [1.0,1.0] #initial condition
tspan = (0.0,10.0) #timespan
p = [1.5,1.0,3.0,1.0] #parameters
prob = ODEProblem(lotka_volterra,u0,tspan,p)

#= Solve the ODE problem using 5th order Tsitouras method (Tsit5) and plot the
solution.
=#
sol = solve(prob,Tsit5())
plot(sol)

# Data generation from the above ODE solution
sol = solve(prob,Tsit5(),saveat=0.1)
A = sol[1,:] # length 101 vector
t = 0:0.1:10.0
scatter!(t,A) #plot (t,A) over the ODE's solution

#= Build a neural network with the function as our single layer,
and define a loss function as the squared distance of the above generated data
from 1.
 =#

p = param([2.2, 1.0, 2.0, 0.4]) # Initial Parameter Vector
function predict_rd() # Our 1-layer neural network
  diffeq_rd(p,prob,Tsit5(),saveat=0.1)[1,:]
end
loss_rd() = sum(abs2,x-1 for x in predict_rd()) # loss function

#= 100 epoch neural network training to minimize loss function and thus get the
optimized parameters.
=#

data = Iterators.repeated((), 100)
opt = ADAM(0.1)
cb = function () #callback function to observe training
  display(loss_rd())
  # using `remake` to re-create our `prob` with current parameters `p`
  display(plot(solve(remake(prob,p=Flux.data(p)),Tsit5(),saveat=0.1),ylim=(0,6)))
end
# Display the ODE with the initial parameter values.
cb()

Flux.train!(loss_rd, [p], data, opt, cb = cb)

#= Flux trains the neural network, finding its parameters (p) that minimize the
cost function: the forward pass of the neural network includes solving an ODE.
=#


text:
# __precompile__()

"""
Provides a basic interface for defining and working with probabilistic graphical models
"""
module ProbabilisticGraphicalModels

using Reexport
@reexport using Distributions
@reexport using DataFrames

export
    ProbabilisticGraphicalModel,

    Assignment,                    # variable assignment type, complete or partial, for a Bayesian Network
    NodeName,                      # variable name type
    NodeNames,                     # vector of names
    NodeNameUnion,                 # either a NodeName or NodeNames

    Sampler,
    GraphSearchStrategy,

    # Inference
    InferenceMethod,
    InferenceState,
    infer,

    # Common Functions that are not in Base or required packages
    markov_blanket,
    is_independent,
    consistent

include("nodenames.jl")
include("assignments.jl")

abstract type ProbabilisticGraphicalModel end

Base.get(pgm::ProbabilisticGraphicalModel, i::Int) = error("get(pgm, Int) not implemented for $(typeof(pgm))!")
Base.get(pgm::ProbabilisticGraphicalModel, name::NodeName) = error("get(pgm, NodeName) not implemented for $(typeof(pgm))!")
Base.get(pgm::ProbabilisticGraphicalModel, names::NodeNames) = [get(pgm, name) for name in names]
Base.get(pgm::ProbabilisticGraphicalModel, names::Base.AbstractSet{NodeName}) = get(pgm, collect(names))

"""
    length(PGM)
Returns the number of variables in the probabilistic graphical model
"""
Base.length(pgm::ProbabilisticGraphicalModel) = error("length not implemented for $(typeof(pgm))!")

"""
    names(PGM)
Returns a list of NodeNames
"""
Base.names(pgm::ProbabilisticGraphicalModel) = error("names not implemented for $(typeof(pgm))!")

"""
markov_blanket(PGM)
Returns the list of NodeNames forming the Markov blanket for the PGM
"""
markov_blanket(pgm::ProbabilisticGraphicalModel) = error("markov_blanket not implemented for $(typeof(pgm))!")

"""
is_independent(PGM, x::NodeNames, y::NodeNames, given::NodeNames)
Returns whether the set of node names `x` is d-separated from the set `y` given the set `given`
"""
is_independent(pgm::ProbabilisticGraphicalModel, x::NodeNames, y::NodeNames, given::NodeNames) = error("is_independent not implemented for $(typeof(pgm))!")

"""
The pdf of a given assignment after conditioning on the values
"""
Distributions.pdf(pgm::ProbabilisticGraphicalModel, assignment::Assignment) = error("pdf not implemented for $(typeof(pgm))!")
Distributions.pdf(pgm::ProbabilisticGraphicalModel, pair::Pair{NodeName}...) = pdf(pgm, Assignment(pair))

"""
The logpdf of a given assignment after conditioning on the values
"""
Distributions.logpdf(pgm::ProbabilisticGraphicalModel, assignment::Assignment) = error("logpdf not implemented for $(typeof(pgm))!")
Distributions.logpdf(pgm::ProbabilisticGraphicalModel, pair::Pair{NodeName}...) = logpdf(pgm, Assignment(pair))

"""
The logpdf of a set of assignment after conditioning on the values
"""
function Distributions.logpdf(pgm::ProbabilisticGraphicalModel, df::DataFrame)

    logl = 0.0

    a = Assignment()
    nodenames = names(pgm)
    for i in 1 : nrow(df)

        for name in nodenames
            a[name] = df[i, name]
        end

        logl += logpdf(pgm, a)
    end

    return logl
end

"""
The pdf of a set of assignments after conditioning on the values
"""
Distributions.pdf(pgm::ProbabilisticGraphicalModel, df::DataFrame) = exp(logpdf(pgm, df))

include("sampling.jl")
include("learning.jl")
include("inference.jl")

end # module ProbabilisticGraphicalModels


text:
abstract type AbstractMeasureSpec end

#holds the parameters for forming the weights
#not meant to hold any meaningful data- just parameters
struct MeasureSpec{TFvol<:NSymbol, TFDvol<:NSymbol} <: AbstractMeasureSpec
  Fξs::Vector{Symbol}
  ξs::Vector{Characteristic}
  ξidx::Dict{Symbol, Characteristic}
  Fvol::TFvol
  FDvol::TFDvol

  FW::Vector{Symbol} #these are global controls (use case: intercept, controls by k or n)
  FWbyt::Vector{Symbol} #these are controls grouped by t only (use case: time fixed effects)
  FWbykt::Dict{Symbol, Vector{Symbol}} #these are controls grouped by t AND k

  K::Int
end

Fweights(ms::MeasureSpec, s::Symbol) = (ξ->ξ.X[s]).(ms.ξs)
Fgroupedcontrols(ms::MeasureSpec) = [ms.FWbyt...;(Fξ->ms.FWbykt[Fξ]).(ms.Fξs)...;]
Fcontrols(ms::MeasureSpec) = [ms.FW...; Fgroupedcontrols(ms)...;]
grossleverages(ms::MeasureSpec) = (ξ->ξ.ws.gross).(ms.ξs)

#creates the control names as needed
function getcontrolnames(Fξs::Vector{Symbol};
    controlsprefix::Symbol = PARAM[:controlsprefix],
    controlsglobal::Vector{Symbol} = PARAM[:controlsglobal],
    controlsbykt::Vector{Symbol} = PARAM[:controlsbykt],
    controlsbyt::Vector{Symbol} = PARAM[:controlsbyt])

  FW = (s->Symbol(controlsprefix,s)).(controlsglobal)
  FWbyt = (s->Symbol(controlsprefix,s)).(controlsbyt)
  FWbykt::Dict{Symbol,Vector{Symbol}} = Dict{Symbol,Vector{Symbol}}()
  for Fξ ∈ Fξs
    FWbykt[Fξ] = [Symbol(controlsprefix, s, Fξ) for s ∈ controlsbykt]
  end

  ntout = (;FW=FW, FWbyt=FWbyt, FWbykt=FWbykt)
  #println("ntout: $(ntout), typeof: $(typeof(ntout))")
  return ntout
end

function MeasureSpec(ξstuple::Characteristic...;
  Fvol::NSymbol = PARAM[:Fvol],
  FDvol::NSymbol = nothing, #set the keyword to Symbol(:D, Fvol) to enable
  ξs::Vector{Characteristic} = collect(ξstuple),
  Fξs::Vector{Symbol} = (ξ->ξ.Fξ).(ξs),
  controlnames::NamedTuple=getcontrolnames(Fξs)
  )

  ξidx::Dict{Symbol, Characteristic} = Dict(zip(Fξs, ξs))

  #permnodate::Vector{Tuple{Int, Date}} = (Tuple).(zip(panel.permno, panel.date))
  #println("controlnames: $(controlnames), typeof: $(typeof(controlnames))")

  return MeasureSpec(Fξs, ξs, ξidx, Fvol, FDvol,
    controlnames.FW, controlnames.FWbyt, controlnames.FWbykt, length(Fξs))
end

#set up a get index notation
Base.getindex(ms::MeasureSpec, Fξ::Symbol) = ms.ξidx[Fξ]
Base.getindex(ms::MeasureSpec, Fξs::Symbol...) = (Fξ->ms.idx[Fξ]).(Fξs)
Base.getindex(ms::MeasureSpec, i::Int) = ms.ξiBase.dx[i]
Base.getindex(ms::MeasureSpec, is::Int...) = ms.ξidx[is]

Base.length(ms::MeasureSpec) = length(Fξs)

#make iteration work
function Base.iterate(ms::MeasureSpec, state::Int = 1)
  (state > ms.K) && return nothing
  return (ms[state], state + 1)
end

#some convenience methods
getT(panel::AbstractDataFrame, ms::MeasureSpec) = length(unique(panel.date))
getN(panel::AbstractDataFrame, ms::MeasureSpec) = length(unique(panel.permno))
getdims(panel::AbstractDataFrame, ms::MeasureSpec) = (
  T=getT(panel, ms), N = getN(panel, ms), K = ms.K)

function gettidx(panel::AbstractDataFrame, ms::MeasureSpec)
  date::Vector{Date} = sort(unique(panel.date))
  tidx::Dict{Date,Int} = Dict(date[t]=>t for t ∈ 1:length(date))

  return tidx
end

#=makes sure the MeasureSpec matches the ordering in the dataframe
function verifymeasureintegrity(panel::DataFrame, ms::MeasureSpec)
  Nrows::Int = size(panel, 1)

  permno::Vector{Int} = (t->t[1]).(ms.permnodate)
  date::Vector{Date} = (t->t[2]).(ms.permnodate)

  (sum(permno .== panel.permno) == Nrows) || error("permno integrity check failed")
  (sum(date .== panel.date) == Nrows) || error("date integrity check failed")
end=#

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
import Obj
import Obj.Objective

gd() = GdOptimizer(

    # name
    "gd",

    # initial state
    EmptyState(),

    # update
    function (obj::Objective, k::Int, w::Array{Float64,1}, b::Array{Int,1}, inf, state::State)
        gw = Obj.g(obj, w)
        return Obj.f(obj, w), gw, -gw, state
    end

)

max_stars_repo_name: mirkobunse/lso

max_issues_repo_name: mirkobunse/lso

max_forks_repo_name: mirkobunse/lso


text:
struct PointScalarValues{dim,T<:Real,refshape<:AbstractRefShape} <: CellValues{dim,T,refshape}
    N::Vector{T}
end

function Base.show(io::IO, ::MIME"text/plain", pv::PointScalarValues)
    print(io, "$(typeof(pv)) with $(getnbasefunctions(pv)) shape functions.")
end

function PointScalarValues(quad_rule::QuadratureRule, func_interpol::Interpolation)
    PointScalarValues(Float64, quad_rule, func_interpol)
end

function PointScalarValues(::Type{T}, quad_rule::QuadratureRule{dim,shape}, func_interpol::Interpolation) where {dim,T,shape<:AbstractRefShape}

    length(getweights(quad_rule)) == 1 || error("PointScalarValues supports only a single point.")

    # Function interpolation
    n_func_basefuncs = getnbasefunctions(func_interpol)
    N    = fill(zero(T)          * T(NaN), n_func_basefuncs)

    ξ = quad_rule.points[1]
    for i in 1:n_func_basefuncs
        N[i] = value(func_interpol, i, ξ)
    end

    PointScalarValues{dim,T,shape}(N)
end

# PointScalarValues only have one quadrature point anyways
function PointScalarValues(coord::Vec{dim,T}, ip::Interpolation{dim, refshape}) where {dim,refshape,T}
    qr = QuadratureRule{dim,refshape,T}([one(T)], [coord])
    return PointScalarValues(qr, ip)
end

# allow to use function_value with any
Base.@pure _valuetype(::PointScalarValues{dim}, ::Vector{T}) where {dim, T<:AbstractTensor} = T

# allow on-the-fly updating
function reinit!(pv::PointScalarValues{dim,T,refshape}, coord::Vec{dim,T}, func_interpol::Interpolation{dim,refshape,order}) where {dim,T,refshape,order}
    n_func_basefuncs = getnbasefunctions(func_interpol)
    for i in 1:n_func_basefuncs
        pv.N[i] = value(func_interpol, i, coord)
    end
    return pv
end
# TODO: need a show method for PointScalarValues


text: # import HttpServer


text:
using FileIO
using Printf: @sprintf
using Profile
using ProfileView
using Geometry: TriangleMesh, CameraIntrinsics, to_matrix, I4, Pose6D, SceneGraphNode

using DepthRenderer

width = 64
height = 64
cam = CameraIntrinsics(width=width, height=height, fx=width, fy=height, cx=div(width, 2), cy=div(height, 2), near=0.001, far=5.0)
r = Renderer(cam)

# mug
mug_mesh = TriangleMesh("mug.obj")
register_mesh!(r, mug_mesh)

# scene graph
root = SceneGraphNode(
    children=[
        SceneGraphNode(mesh=mug_mesh, transform=to_matrix(Pose6D(0., 0., -0.5, 0., 0., 0.)))])

model = I4(Float32)
view = I4(Float32)

function do_render_test(n::Int)
    for i=1:n
        #println("i: $i")
        draw!(r, root, model, view)
        depth_image = get_depth_image!(r; show_in_window=false)
        #depth_image = depth_image'[end:-1:1,:]
        #save(@sprintf("imgs/depth-%04d.png", i), depth_image ./ maximum(depth_image))
    end
end

#@time do_render_test(500)
@time do_render_test(10000)
@time do_render_test(10000)

#Profile.clear()
#@profile do_render_test(10000)
#li, lidict = Profile.retrieve()
#using JLD
#@save "profdata.jld" li lidict

destroy!(r)

max_stars_repo_name: probcomp/DepthRenderer

max_issues_repo_name: probcomp/DepthRenderer

max_forks_repo_name: probcomp/DepthRenderer


text:
# Julia Implementation of Game of Life
# Any live cell with fewer than two live neighbors dies, as if caused by under-population.
# Any live cell with two or three live neighbors lives on to the next generation.
# Any live cell with more than three live neighbors dies, as if by over-population..
# Any dead cell with exactly three live neighbors becomes a live cell, as if by reproduction.

"""
function count_alive :- Calculates the number of alive neighbours present for a cell(x,y) on the board
Input:
    board :- 2D Grid/board
    row :- x coordinate of the current cell
    column :- y coordinate of the current cell

Output:
    count :- Total number of alive neighbours of the current cell
"""
function count_alive(board,row,column)
    check = [[0 1], [1 0], [1 1], [0 -1], [-1 0], [1 -1], [-1 1],[-1,-1]]
    count = 0
    for c in check
        if (row + c[1]) <= 0 || (column + c[2]) <= 0 || (row + c[1]) > size(board)[1] || (column + c[2]) > size(board[1])[1]
            continue
        end
        if board[row + c[1]][column + c[2]] == 1
            count += 1
        end
    end

    return count

end

"""
function gameOfLife :- Generates the next cell structure for the grid
Input:
    board :- Input 2D Grid/board
Output:
    board :- Next generated cell structure according to the rules
"""
function gameOfLife(board,generations)
    new_board = deepcopy(board)
    rows = size(board)
    column = size(board[1])
    for generation in 1:generations
        changed = []
        for i in 1:rows[1]
            for j in 1:column[1]
                neighbours = count_alive(new_board, i, j)
                if new_board[i][j] == 1
                    if neighbours < 2 || neighbours>3
                        append!(changed,[[i,j]])
                    end
                else
                    if neighbours==3
                        append!(changed,[[i,j]])
                    end
                end
            end
        end

        for toggle_index in changed
            x = toggle_index[1]
            y = toggle_index[2]
            new_board[x][y]  = 1 - new_board[x][y]
        end

        #print(new_board)

    end
    return new_board
end

function print_as_grid(input_array,output_array,expected,count,generations,failed)
    println("********************************************************")
    println("Test Case:-",count)
    if output_array==expected
        println("Test Case ", count," has passed")
    else
        println("Test Case ", count," has failed")
        append!(failed,count)
    end

    println("Input Array ===> Output Array ===> Expected Output")
    println("This test case has ",generations," generations")
    println(input_array[1]," ===>",output_array[1]," ===>",expected[1])
    println(input_array[2]," ===>",output_array[2]," ===>",expected[2])
    println(input_array[3]," ===>",output_array[3]," ===>",expected[3])
    println(input_array[4]," ===>",output_array[4]," ===>",expected[4])
    println("********************************************************")
end

function run_test_cases()
    input1 = [[0,1,0],[0,0,1],[1,1,1],[0,0,0]]
    expected1 = [[0,0,0],[1,0,1],[0,1,1],[0,1,0]]
    input2 = [[1,1,0,1],[0,0,0,1],[1,0,1,1],[0,1,0,0]]
    expected2 = [[0,0,1,0],[1,0,0,1],[0,1,1,1],[0,1,1,0]]
    input3 = [[1,0,0,1],[1,1,0,1],[1,0,0,1],[0,1,1,0]]
    expected3 = [[1,1,1,0],[1,1,0,1],[1,0,0,1],[0,1,1,0]]
    input4 = [[1,0,0],[0,0,1],[1,0,0],[0,0,1]]
    expected4 = [[0,0,0],[0,1,0],[0,1,0],[0,0,0]]
    input5 = [[1,1,0,0,0],[0,0,0,0,1],[0,0,0,1,1],[0,0,0,0,0]]
    expected5 = [[0,0,0,0,0],[0,0,0,1,1],[0,0,0,1,1],[0,0,0,0,0]]
    input6 = [[0,1,0],[0,1,0],[0,1,0],[0,0,0]]
    expected6 = [[0,1,0],[0,1,0],[0,1,0],[0,0,0]]

    failed = []
    output1 = gameOfLife(input1,1)
    output2 = gameOfLife(input2,1)
    output3 = gameOfLife(input3,1)
    output4 = gameOfLife(input4,1)
    output5 = gameOfLife(input5,1)
    output6 = gameOfLife(input6,2)

    print_as_grid(input1,output1,expected1,1,1,failed)
    print_as_grid(input2,output2,expected2,2,1,failed)
    print_as_grid(input3,output3,expected3,3,1,failed)
    print_as_grid(input4,output4,expected4,4,1,failed)
    print_as_grid(input5,output5,expected5,5,1,failed)
    print_as_grid(input6,output6,expected6,6,2,failed)
    println()
    if size(failed)[1]==0
        println("Congratulations, all test Cases Passed")
    else
        println("Out of 6 test cases the following test cases failed:-")
        println(failed)
    end
end

#testing
run_test_cases()


text:
mutable struct SingleIndexModel

    # model -
    α::Float64          # firm specific unexplained return
    β::Float64          # relationship between the firm and the market
    r::Float64          # risk free rate of return
    ϵ::Distribution     # random shocks

    # constructor -
    SingleIndexModel() = new()
end

mutable struct GeometricBrownianMotionModel

    # model parameters -
    μ::Float64
    σ::Float64

    GeometricBrownianMotionModel() = new()
end

# --- ASSETS ------------------------------------------------- #
abstract type AbstractAsset end
mutable struct CallContract <: AbstractAsset

    # data -
    ticker::String
    expiration::Date
    K::Float64
    C::Float64
    sense::Symbol

    # constructor
    CallContract() = new()
end

mutable struct PutContract <: AbstractAsset

    # data -
    ticker::String
    expiration::Date
    K::Float64
    C::Float64
    sense::Symbol

    # constructor
    PutContract() = new()
end

mutable struct Equity <: AbstractAsset

    # data -
    ticker::String
    C::Float64
    sense::Symbol

    # constructor
    Equity() = new()
end

# --- Polygon API model --------------------------------------------- #
abstract type AbstractPolygonEndpointModel end
mutable struct PolygonAggregatesEndpointModel <: AbstractPolygonEndpointModel

    # data -
    adjusted::Bool
    limit::Int64
    sortdirection::String
    apikey::String
    ticker::String
    to::Date
    from::Date
    multiplier::Int64
    timespan::String

    # constructor -
    PolygonAggregatesEndpointModel() = new()
end

# Setup call for options data -
@enum PolygonOptionSortKey ticker=1 underlying_ticker expiration_date strike_price
@enum PolygonOptionContractType call=1 put
mutable struct PolygonOptionsContractReferenceEndpoint <: AbstractPolygonEndpointModel

    # data -
    ticker::Union{Nothing,String}
    underlying_ticker::Union{Nothing, String}
    contract_type::PolygonOptionContractType
    expiration_date::Date
    limit::Int64
    order::String
    sort::PolygonOptionSortKey
    apikey::String

    # constructor -
    PolygonOptionsContractReferenceEndpoint() = new()
end

max_stars_repo_name: Pooksoft/Serenity

max_issues_repo_name: Alba-Intelligence/Serenity

max_forks_repo_name: Alba-Intelligence/Serenity


text:
using EasyStream
using DataFrames
using Test

include("connector.jl")
include("stream.jl")
include("modifiers.jl")


text:
"""
    load(fpath::String; gamma::Union{Nothing,Float64}=nothing)

Read a `.png` image from file at `fpath`. `gamma` can be used to override the automatic gamma
correction, a value of 1.0 means no gamma correction.
Returns a matrix.

The result will be an 8 bit (N0f8) image if the source bit depth is <= 8 bits, 16 bit (N0f16)
otherwise.
The number of channels (and transparency) of the source determines the color type of the output:
    1 channel  -> Gray
    2 channels -> GrayA
    3 channels -> RGB
    4 channels -> RGBA
"""
function load(fpath::String; gamma::Union{Nothing,Float64}=nothing)
    fp = open_png(fpath)
    png_ptr = create_read_struct()
    info_ptr = create_info_struct(png_ptr)
    png_init_io(png_ptr, fp)
    png_set_sig_bytes(png_ptr, PNG_BYTES_TO_CHECK)
    png_read_info(png_ptr, info_ptr)

    width = png_get_image_width(png_ptr, info_ptr)
    height = png_get_image_height(png_ptr, info_ptr)
    color_type_orig = png_get_color_type(png_ptr, info_ptr)
    color_type = color_type_orig
    bit_depth_orig = png_get_bit_depth(png_ptr, info_ptr)
    bit_depth = bit_depth_orig
    num_channels = png_get_channels(png_ptr, info_ptr)
    interlace_type = png_get_interlace_type(png_ptr, info_ptr)

    # TODO: verify this is needed
    backgroundp = png_color_16p()
    if png_get_bKGD(png_ptr, info_ptr, Ptr{png_color_16p}(backgroundp)) != 0
        png_set_background(png_ptr, backgroundp, PNG_BACKGROUND_GAMMA_FILE, 1, 1.0)
    end

    screen_gamma = PNG_DEFAULT_sRGB
    image_gamma = Ref{Cdouble}(-1.0)
    intent = Ref{Cint}(-1)
    if isnothing(gamma)
        if png_get_valid(png_ptr, info_ptr, PNG_INFO_sRGB) != 0
            if png_get_sRGB(png_ptr, info_ptr, intent) != 0
                png_set_gamma(png_ptr, screen_gamma, PNG_DEFAULT_sRGB);
            else
                if png_get_gAMA(png_ptr, info_ptr, image_gamma) != 0
                    png_set_gamma(png_ptr, screen_gamma, image_gamma[])
                else
                    image_gamma[] = 0.45455
                    png_set_gamma(png_ptr, screen_gamma, image_gamma[])
                end
            end
        elseif png_get_valid(png_ptr, info_ptr, PNG_INFO_gAMA) != 0
            if png_get_gAMA(png_ptr, info_ptr, image_gamma) != 0
                png_set_gamma(png_ptr, screen_gamma, image_gamma[])
            else
                image_gamma[] = 0.45455
                png_set_gamma(png_ptr, screen_gamma, image_gamma[])
            end
        end
    elseif gamma != 1
        image_gamma[] = gamma
        png_set_gamma(png_ptr, screen_gamma, image_gamma[])
    end

    if color_type == PNG_COLOR_TYPE_PALETTE
        png_set_palette_to_rgb(png_ptr)
        color_type = PNG_COLOR_TYPE_RGB
    end

    if color_type == PNG_COLOR_TYPE_GRAY && bit_depth < 8
        png_set_expand_gray_1_2_4_to_8(png_ptr)
        png_set_packing(png_ptr)
        bit_depth = UInt8(8)
    end

    if png_get_valid(png_ptr, info_ptr, PNG_INFO_tRNS) != 0
        png_set_tRNS_to_alpha(png_ptr)
        if color_type == PNG_COLOR_TYPE_GRAY || color_type == PNG_COLOR_TYPE_RGB
            color_type |= PNG_COLOR_MASK_ALPHA
        end
    end

    buffer_eltype = _buffer_color_type(color_type, bit_depth)
    n_passes = png_set_interlace_handling(png_ptr)
    bit_depth == 16 && png_set_swap(png_ptr)
    png_read_update_info(png_ptr, info_ptr)

    # We transpose to work around libpng expecting row-major arrays
    buffer = Array{buffer_eltype}(undef, width, height)

    @debug(
        "Read PNG info:",
        fpath,
        height,
        width,
        color_type_orig,
        color_type,
        bit_depth_orig,
        bit_depth,
        num_channels,
        interlace_type,
        gamma,
        image_gamma[],
        screen_gamma,
        intent[],
        n_passes,
        buffer_eltype,
        PNG_HEADER_VERSION_STRING
    )

    png_read_image(png_ptr, map(pointer, eachcol(buffer)))
    png_read_end(png_ptr, info_ptr)
    png_destroy_read_struct(Ref{Ptr{Cvoid}}(png_ptr), Ref{Ptr{Cvoid}}(info_ptr), C_NULL)
    close_png(fp)
    return transpose(buffer)
end

function _buffer_color_type(color_type, bit_depth)
    bit_depth = Int(bit_depth)
    if color_type == PNG_COLOR_TYPE_GRAY
        colors_type = Gray{bit_depth > 8 ? Normed{UInt16,bit_depth} : Normed{UInt8,bit_depth}}
    elseif color_type == PNG_COLOR_TYPE_PALETTE
        colors_type = RGB{bit_depth == 16 ? N0f16 : N0f8}
    elseif color_type == PNG_COLOR_TYPE_RGB
        colors_type = RGB{bit_depth == 16 ? N0f16 : N0f8}
    elseif color_type == PNG_COLOR_TYPE_RGB_ALPHA
        colors_type = RGBA{bit_depth == 16 ? N0f16 : N0f8}
    elseif color_type == PNG_COLOR_TYPE_GRAY_ALPHA
        colors_type = GrayA{bit_depth > 8 ? Normed{UInt16,bit_depth} : Normed{UInt8,bit_depth}}
    else
        throw(error("Unknown color type: $color_type"))
    end
    return colors_type
end

### Write ##########################################################################################
"""
    save(fpath::String, image::AbstractArray;
         compression_level::Integer=0, compression_strategy::Integer=3, filters::Integer=4)

Writes `image` as a png to file at `fpath`.

## Arguments
`compression_level`: 0 (Z_NO_COMPRESSION), 1 (Z_BEST_SPEED), ..., 9 (Z_BEST_COMPRESSION)
`compression_strategy`: 0 (Z_DEFAULT_STRATEGY), 1 (Z_FILTERED), 2 (Z_HUFFMAN_ONLY), 3 (Z_RLE), 4 (Z_FIXED)
`filters`: 0 (None), 1 (Sub), 2 (Up), 3 (Average), 4 (Paeth).

The saved image will have 16 bits of depth if the `image` has eltype that is based on `UInt16`,
8 bits otherwise.
The number of channels and element type of the `image` determines the color type of the
output:
    0/1 channel Float / Integer / Normed or Gray eltype        -> PNG_COLOR_TYPE_GRAY
    2 channels Float  / Integer / Normed or GrayA eltype       -> PNG_COLOR_TYPE_GRAY_ALPHA
    3 channels Float  / Integer / Normed or RGB / BGR eltype   -> PNG_COLOR_TYPE_RGB
    4 channels Float  / Integer / Normed or ARGB / ABGR eltype -> PNG_COLOR_TYPE_RGB_ALPHA
"""
function save(
        fpath,
        image::S;
        compression_level::Integer=Z_NO_COMPRESSION,
        compression_strategy::Integer=Z_RLE,
        filters::Integer=Int(PNG_FILTER_PAETH)
    ) where {
        T,
        S<:Union{AbstractMatrix,AbstractArray{T,3}}
    }
    @assert Z_DEFAULT_STRATEGY <= compression_strategy <= Z_FIXED
    @assert Z_NO_COMPRESSION <= compression_level <= Z_BEST_COMPRESSION
    @assert 2 <= ndims(image) <= 3
    @assert size(image, 3) <= 4

    fp = ccall(:fopen, Ptr{Cvoid}, (Cstring, Cstring), fpath, "wb")
    fp == C_NULL && error("Could not open $(fpath) for writing")

    png_ptr = create_write_struct(png_error_fn, png_warn_fn)
    info_ptr = create_info_struct(png_ptr)
    png_init_io(png_ptr, fp)
    png_set_filter(png_ptr, PNG_FILTER_TYPE_BASE, UInt32(filters))
    png_set_compression_level(png_ptr, compression_level)
    png_set_compression_strategy(png_ptr, compression_strategy)

    height, width = size(image)[1:2]
    bit_depth = _get_bit_depth(image)
    color_type = _get_color_type(image)
    interlace = PNG_INTERLACE_NONE
    compression_type = PNG_COMPRESSION_TYPE_BASE
    filter_type = PNG_FILTER_TYPE_BASE

    elt = eltype(image)
    if (elt <: BGR || elt <: BGRA || elt <: ABGR)
       png_set_bgr(png_ptr)
    end

    if (elt <: ABGR || elt <: ARGB)
        png_set_swap_alpha(png_ptr)
    end

    if color_type == PNG_COLOR_TYPE_GRAY && bit_depth < 8
        png_set_packing(png_ptr)
        bit_depth = 8  # TODO: support 1, 2, 4 bit-depth gray images
    end

    # this gAMA and cHRM should be added for compatibility with older systems
    png_set_sRGB_gAMA_and_cHRM(png_ptr, info_ptr, 0)

    @debug(
        "Write PNG info:",
        fpath,
        height,
        width,
        bit_depth,
        color_type,
        interlace,
        compression_type,
        filter_type,
        filters,
        compression_level,
        compression_strategy,
        typeof(image)
    )

    # TODO: on error this throws an abort signal because we currently don't handle `longjmp`
    png_set_IHDR(
        png_ptr,
        info_ptr,
        width,
        height,
        bit_depth,
        color_type,
        interlace,
        compression_type,
        filter_type,
    )
    png_write_info(png_ptr, info_ptr)
    bit_depth == 16 && png_set_swap(png_ptr) # Handles endianness for 16 bit

    # We transpose to work around libpng expecting row-major arrays
    _write_image(transpose(_prepare_buffer(image)), png_ptr, info_ptr)

    png_destroy_write_struct(Ref(png_ptr), Ref(info_ptr))
    close_png(fp)
end

function _write_image(buf::AbstractArray{T,2}, png_ptr::Ptr{Cvoid}, info_ptr::Ptr{Cvoid}) where {T}
    ccall(
        (:png_write_image, libpng),
        Cvoid,
        (Ptr{Cvoid}, Ptr{Ptr{T}}),
        png_ptr,
        map(pointer, eachcol(buf)),
    )
    png_write_end(png_ptr, info_ptr)
end

_prepare_buffer(x::BitArray) = _prepare_buffer(collect(x))
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:Colorant{<:Normed}} = x
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:UInt8} = reinterpret(Gray{N0f8}, x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:UInt16} = reinterpret(Gray{N0f16}, x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:Normed} = reinterpret(Gray{T}, x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:Gray{<:AbstractFloat}} =
    convert(Array{Gray{N0f8}}, x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:GrayA{<:AbstractFloat}} =
    convert(Array{GrayA{N0f8}}, x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:RGB{<:AbstractFloat}} =
    convert(Array{RGB{N0f8}}, x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:RGBA{<:AbstractFloat}} =
    convert(Array{RGBA{N0f8}}, x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:BGR{<:AbstractFloat}} =
    convert(Array{BGR{N0f8}}, x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:BGRA{<:AbstractFloat}} =
    convert(Array{BGRA{N0f8}}, x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:ARGB{<:AbstractFloat}} =
    convert(Array{ARGB{N0f8}}, x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:ABGR{<:AbstractFloat}} =
    convert(Array{ABGR{N0f8}}, x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:Union{AbstractFloat,Bool}} =
    reinterpret(Gray{N0f8}, convert(Array{N0f8}, x))
_prepare_buffer(x::AbstractArray{T,3}) where {T<:Union{AbstractFloat,Bool}} =
    __prepare_buffer(convert(Array{N0f8}, x))
_prepare_buffer(x::AbstractArray{T,3}) where {T<:Union{UInt8,Int8}} =
    __prepare_buffer(reinterpret(N0f8, x))
_prepare_buffer(x::AbstractArray{T,3}) where {T<:Union{UInt16,Int16}} =
    __prepare_buffer(reinterpret(N0f16, x))
_prepare_buffer(x::AbstractArray{T,3}) where {T<:Normed} = __prepare_buffer(x)

function __prepare_buffer(x::AbstractArray{T,3}) where {T}
    nchannels = size(x, 3)
    if nchannels == 1
        ifelse(ndims(x) == 3, _prepare_buffer(dropdims(x, dims=3)), x)
    elseif nchannels == 2
        convert(Array{GrayA}, colorview(GrayA, view(x, :, :, 1), view(x, :, :, 2)))
    elseif nchannels == 3
        convert(Array{RGB}, colorview(RGB, view(x, :, :, 1), view(x, :, :, 2), view(x, :, :, 3)))
    elseif nchannels == 4
        convert(
            Array{RGBA},
            colorview(RGBA, view(x, :, :, 1), view(x, :, :, 2), view(x, :, :, 3), view(x, :, :, 4)),
        )
    else
        error("Unsupported number of channels $(nchannels) in input. Only <= 4 is expected.")
    end
end

_get_bit_depth(img::BitArray) = 8  # TODO: write 1 bit-depth images
_get_bit_depth(img::AbstractArray{C}) where {C<:Colorant} = __get_bit_depth(eltype(C))
_get_bit_depth(img::AbstractArray{T}) where {T<:Normed} = __get_bit_depth(T)
# __get_bit_depth(::Type{Normed{T,1}}) where T = 1  # TODO: write 1 bit-depth images
# __get_bit_depth(::Type{Normed{T,2}}) where T = 2  # TODO: write 2 bit-depth images
# __get_bit_depth(::Type{Normed{T,4}}) where T = 4  # TODO: write 4 bit-depth images
__get_bit_depth(::Type{Normed{T,8}}) where T = 8
__get_bit_depth(::Type{Normed{T,16}}) where T = 16
__get_bit_depth(::Type{Normed{T,N}}) where {T,N} = ifelse(N <= 8, 8, 16)
__get_bit_depth(::Type{<:AbstractFloat}) = 8
_get_bit_depth(img::AbstractArray{T}) where {T<:AbstractFloat} = 8
_get_bit_depth(img::AbstractArray{<:Bool}) = 8  # TODO: write 1 bit-depth images
_get_bit_depth(img::AbstractArray{<:UInt8}) = 8
_get_bit_depth(img::AbstractArray{<:UInt16}) = 16

_get_color_type(x::AbstractArray{<:Gray{T}}) where {T} = PNG_COLOR_TYPE_GRAY
_get_color_type(x::AbstractArray{<:GrayA{T}}) where {T} = PNG_COLOR_TYPE_GRAY_ALPHA
_get_color_type(x::AbstractArray{<:RGB{T}}) where {T} = PNG_COLOR_TYPE_RGB
_get_color_type(x::AbstractArray{<:RGBA{T}}) where {T} = PNG_COLOR_TYPE_RGBA
_get_color_type(x::AbstractArray{<:BGR{T}}) where {T} = PNG_COLOR_TYPE_RGB
_get_color_type(x::AbstractArray{<:BGRA{T}}) where {T} = PNG_COLOR_TYPE_RGBA
_get_color_type(x::AbstractArray{<:ARGB{T}}) where {T} = PNG_COLOR_TYPE_RGBA
_get_color_type(x::AbstractArray{<:ABGR{T}}) where {T} = PNG_COLOR_TYPE_RGBA
function _get_color_type(
        x::AbstractArray{T, N}
    ) where {
        T<:Union{Normed,Unsigned,Bool,AbstractFloat},
        N
    }
    if N == 2
        return PNG_COLOR_TYPE_GRAY
    elseif N == 3
        d = size(x, 3)
        d == 1 && (return PNG_COLOR_TYPE_GRAY)
        d == 2 && (return PNG_COLOR_TYPE_GRAY_ALPHA)
        d == 3 && (return PNG_COLOR_TYPE_RGB)
        d == 4 && (return PNG_COLOR_TYPE_RGBA)
    end
    error("Number of dimensions $(N) in image not supported (only 2D or 3D Arrays are expected).")
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

module MolecularGraph

    export
        MolecularGraphUtil,
        MolecularGraphGeometry,
        MolecularGraphModel

    module MolecularGraphUtil
        include("./util/meta.jl")
        include("./util/iterator.jl")
    end

    module MolecularGraphModel
        using ..MolecularGraphUtil

        include("./graph/interface.jl")

        include("./graph/plaingraph.jl")
        include("./graph/plaindigraph.jl")
        include("./graph/plainhypergraph.jl")

        include("./graph/generator.jl")
        include("./graph/multigraph.jl")
        include("./graph/dag.jl")
        include("./graph/inducedsubgraph.jl")
        include("./graph/product.jl")
        include("./graph/disjointunion.jl")
        include("./graph/linegraph.jl")

        include("./graph/dfs.jl")
        include("./graph/shortestpath.jl")
        include("./graph/bipartite.jl")
        include("./graph/triangle.jl")
        include("./graph/clique.jl")
        include("./graph/connectivity.jl")
        include("./graph/cycle.jl")
        include("./graph/planarity.jl")

        include("./graph/isomorphism/edgeinduced.jl")
        include("./graph/isomorphism/vf2.jl")
        include("./graph/isomorphism/cliquemcs.jl")
    end

    module MolecularGraphGeometry
        using LinearAlgebra
        using ..MolecularGraphUtil
        using ..MolecularGraphModel

        include("./geometry/interface.jl")

        include("./geometry/coords2d.jl")
        include("./geometry/coords3d.jl")
        include("./geometry/coordsinternal.jl")
    end

    using LinearAlgebra
    using Statistics
    using YAML
    using ..MolecularGraphUtil
    using ..MolecularGraphGeometry
    using ..MolecularGraphModel

    include("./model/interface.jl")
    include("./draw/interface.jl")

    include("./model/atom.jl")
    include("./model/bond.jl")
    include("./model/molgraph.jl")

    include("substructure.jl")
    include("mcs.jl")
    include("preprocess.jl")
    include("properties.jl")
    include("wclogp.jl")
    include("funcgroup.jl")

    include("./smarts/base.jl")
    include("./smarts/atom.jl")
    include("./smarts/bond.jl")
    include("./smarts/logicaloperator.jl")
    include("./smarts/molecule.jl")

    include("./draw/draw2d.jl")
    include("./draw/svg.jl")

    include("download.jl")
    include("sdfilereader.jl")
    include("sdfilewriter.jl")

end


text:
function max_pooling2d_fwd!{T}(x::AbstractArray{T,4}, y::AbstractArray{T,4},
                               width::Int, height::Int, channels::Int, num::Int, pooled_width::Int,
                               pooled_height::Int, kernel_w::Int, kernel_h::Int, pad_w::Int, pad_h::Int,
                               stride_w::Int, stride_h::Int)
  for n = 1:num, c = 1:channels, ph = 1:pooled_height, pw = 1:pooled_width
    hstart = (ph - 1)*stride_h - pad_h
    wstart = (pw - 1)*stride_w - pad_w
    hend   = min(hstart + kernel_h, height)
    wend   = min(wstart + kernel_w, width)

    hstart = max(hstart, 0) + 1
    wstart = max(wstart, 0) + 1

    y[pw, ph, c, n] = maximum(x[wstart:wend, hstart:hend, c, n])
  end
end

function maxpool2d!(y::AbstractArray{T,4}, x::AbstractArray{T,4};
                    window::Dims{2}=(2,2), padding::Dims{2}=(0,0),
                    stride::Dims{2}=window) where T
    Wx,Hx,Cx,Nx = size(x)
    Wy,Hy,Cy,Ny = size(y)
    (w1,w2) = window
    (p1,p2) = padding
    (s1,s2) = stride
    max_pooling2d_fwd!(x,y,Wx,Hx,Cx,Nx,Wy,Hy,w1,w2,p1,p2,s1,s2)
    return y
end

function max_pooling2d_bwd!{T}(x::AbstractArray{T,4}, y::AbstractArray{T,4},
  grad_output::AbstractArray{T,4}, grad_input::AbstractArray{T,4}, width::Int, height::Int,
  channels::Int, num::Int, pooled_width::Int, pooled_height::Int, kernel_w::Int,
  kernel_h::Int, pad_w::Int, pad_h::Int, stride_w::Int, stride_h::Int)

  grad_input[:, :, :, :] = 0
  #pragma omp parallel for
  for n = 1:num, c = 1:channels, ph = 1:pooled_height, pw = 1:pooled_width
    hstart = (ph - 1) * stride_h - pad_h
    wstart = (pw - 1) * stride_w - pad_w
    hend   = min(hstart + kernel_h, height)
    wend   = min(wstart + kernel_w, width)
    hstart = max(hstart, 0) + 1
    wstart = max(wstart, 0) + 1
    maxval = y[pw, ph, c, n]
    d_maxval = grad_output[pw, ph, c, n]
    for h = hstart:hend, w = wstart:wend
      if x[w, h, c, n] == maxval
        grad_input[w, h, c, n] += d_maxval
      end
    end
  end
end

function maxpool2d_grad!(dx::AbstractArray{T,4}, dy::AbstractArray{T,4}, y::AbstractArray{T,4}, x::AbstractArray{T,4};
                         window::Dims{2}=(2,2), padding::Dims{2}=(0,0),
                         stride::Dims{2}=window) where T
    Wx,Hx,Cx,Nx = size(x)
    Wy,Hy,Cy,Ny = size(y)
    (w1,w2) = window
    (p1,p2) = padding
    (s1,s2) = stride
    max_pooling2d_bwd!(x,y,dy,dx,Wx,Hx,Cx,Nx,Wy,Hy,w1,w2,p1,p2,s1,s2)
    return dx
end

function mean_pooling2d_fwd!{T}(x::AbstractArray{T,4}, y::AbstractArray{T,4},
  width::Int, height::Int, channels::Int, num::Int, pooled_width::Int,
  pooled_height::Int, kernel_w::Int, kernel_h::Int,pad_w::Int, pad_h::Int,
  stride_w::Int, stride_h::Int)
  kernel_size = kernel_w * kernel_h
  for n = 1:num, c = 1:channels, ph = 1:pooled_height, pw = 1:pooled_width
    hstart = (ph - 1) * stride_h - pad_h
    wstart = (pw - 1) * stride_w - pad_w
    hend   = min(hstart + kernel_h, height)
    wend   = min(wstart + kernel_w, width)

    hstart = max(hstart, 0) + 1
    wstart = max(wstart, 0) + 1

    y[pw, ph, c, n] = sum(x[wstart:wend, hstart:hend, c, n]) / kernel_size
  end
end

function meanpool2d!(y::AbstractArray{T,4}, x::AbstractArray{T,4};
                     window::Dims{2}=(2,2), padding::Dims{2}=(0,0),
                     stride::Dims{2}=window) where T
    Wx,Hx,Cx,Nx = size(x)
    Wy,Hy,Cy,Ny = size(y)
    (w1,w2) = window
    (p1,p2) = padding
    (s1,s2) = stride
    mean_pooling2d_fwd!(x,y,Wx,Hx,Cx,Nx,Wy,Hy,w1,w2,p1,p2,s1,s2)
    return y
end

function mean_pooling2d_bwd!{T}(x::AbstractArray{T,4}, y::AbstractArray{T,4},
  width::Int, height::Int, channels::Int, num::Int, pooled_width::Int,
  pooled_height::Int, kernel_w::Int, kernel_h::Int, pad_w::Int, pad_h::Int,
  stride_w::Int, stride_h::Int)

  x[:, :, :, :] = 0
  kernel_size = kernel_w * kernel_h

  #pragma omp parallel for
  for n = 1:num, c = 1:channels, ph = 1:pooled_height, pw = 1:pooled_width
    hstart = (ph - 1) * stride_h - pad_h
    wstart = (pw - 1) * stride_w - pad_w
    hend   = min(hstart + kernel_h, height)
    wend   = min(wstart + kernel_w, width)
    hstart = max(hstart, 0) + 1
    wstart = max(wstart, 0) + 1

    oval = y[pw, ph, c, n] / kernel_size
    x[wstart:wend, hstart:hend, c, n] += oval
  end
end

function meanpool2d_grad!(dx::AbstractArray{T,4}, dy::AbstractArray{T,4}, y::AbstractArray{T,4}, x::AbstractArray{T,4};
                         window::Dims{2}=(2,2), padding::Dims{2}=(0,0),
                         stride::Dims{2}=window) where T
    Wx,Hx,Cx,Nx = size(x)
    Wy,Hy,Cy,Ny = size(y)
    (w1,w2) = window
    (p1,p2) = padding
    (s1,s2) = stride
    mean_pooling2d_bwd!(dx,dy,Wx,Hx,Cx,Nx,Wy,Hy,w1,w2,p1,p2,s1,s2)
    return dx
end

function max_pooling3d_fwd!{T}(x::AbstractArray{T,5}, y::AbstractArray{T,5},
  width::Int, height::Int, depth::Int, channels::Int, num::Int, pooled_width::Int,
  pooled_height::Int, pooled_depth::Int, kernel_w::Int, kernel_h::Int, kernel_d::Int,
  pad_w::Int, pad_h::Int, pad_d::Int, stride_w::Int, stride_h::Int, stride_d::Int)
  for n = 1:num, c = 1:channels, pd = 1:pooled_depth, ph = 1:pooled_height, pw = 1:pooled_width
    dstart = (pd - 1)* stride_d - pad_d
    hstart = (ph - 1)* stride_h - pad_h
    wstart = (pw - 1)* stride_w - pad_w

    dend   = min(dstart + kernel_d, depth)
    hend   = min(hstart + kernel_h, height)
    wend   = min(wstart + kernel_w, width)

    dstart = max(dstart, 0) + 1
    hstart = max(hstart, 0) + 1
    wstart = max(wstart, 0) + 1

    y[pw, ph, pd, c, n] =
    maximum(x[wstart:wend, hstart:hend, dstart:dend, c, n])
  end
end

function maxpool3d!(y::AbstractArray{T,5}, x::AbstractArray{T,5};
                    window::Dims{3}=(2,2,2), padding::Dims{3}=(0,0,0),
                    stride::Dims{3}=window) where T
    Wx,Hx,Dx,Cx,Nx = size(x)
    Wy,Hy,Dy,Cy,Ny = size(y)
    (w1,w2,w3) = psize(window, x)
    (p1,p2,p3) = psize(padding, x)
    (s1,s2,s3) = psize(stride, x)
    max_pooling3d_fwd!(x,y,Wx,Hx,Dx,Cx,Nx,Wy,Hy,Dy,w1,w2,w3,p1,p2,p3,s1,s2,s3)
    return y
end

function max_pooling3d_bwd!{T}(x::AbstractArray{T,5}, y::AbstractArray{T,5},
  grad_output::AbstractArray{T,5}, grad_input::AbstractArray{T,5}, width::Int, height::Int, depth::Int,
  channels::Int, num::Int, pooled_width::Int, pooled_height::Int, pooled_depth::Int,
  kernel_w::Int, kernel_h::Int, kernel_d::Int, pad_w::Int, pad_h::Int, pad_d::Int,
  stride_w::Int, stride_h::Int, stride_d::Int)

  grad_input[:, :, :, :, :] = 0

  #pragma omp parallel for
  for n = 1:num, c = 1:channels, pd = 1:pooled_depth, ph = 1:pooled_height, pw = 1:pooled_width
    dstart = (pd - 1) * stride_h - pad_h
    hstart = (ph - 1) * stride_h - pad_h
    wstart = (pw - 1) * stride_w - pad_w

    dend   = min(dstart + kernel_d, depth)
    hend   = min(hstart + kernel_h, height)
    wend   = min(wstart + kernel_w, width)

    dstart = max(dstart, 0) + 1
    hstart = max(hstart, 0) + 1
    wstart = max(wstart, 0) + 1

    maxval = y[pw, ph, pd, c, n]
    d_maxval = grad_output[pw, ph, pd, c, n]
    for d = dstart:dend, h = hstart:hend, w = wstart:wend
      if x[w, h, d, c, n] == maxval
        grad_input[w, h, d, c, n] += d_maxval
      end
    end
  end
end

function maxpool3d_grad!(dx::AbstractArray{T,5}, dy::AbstractArray{T,5}, y::AbstractArray{T,5}, x::AbstractArray{T,5};
                         window::Dims{3}=(2,2,2), padding::Dims{3}=(0,0,0),
                         stride::Dims{3}=window) where T
    Wx,Hx,Dx,Cx,Nx = size(x)
    Wy,Hy,Dy,Cy,Ny = size(y)
    (w1,w2,w3) = psize(window, x)
    (p1,p2,p3) = psize(padding, x)
    (s1,s2,s3) = psize(stride, x)
    max_pooling3d_bwd!(x,y,dy,dx,Wx,Hx,Dx,Cx,Nx,Wy,Hy,Dy,w1,w2,w3,p1,p2,p3,s1,s2,s3)
    return dx
end

function mean_pooling3d_fwd!{T}(x::AbstractArray{T,5}, y::AbstractArray{T,5},
  width::Int, height::Int, depth::Int, channels::Int, num::Int, pooled_width::Int,
  pooled_height::Int, pooled_depth::Int, kernel_w::Int, kernel_h::Int, kernel_d::Int,
  pad_w::Int, pad_h::Int, pad_d::Int, stride_w::Int, stride_h::Int, stride_d::Int)

  kernel_size = kernel_w * kernel_h * kernel_d
  #pragma omp parallel for
  for n = 1:num, c = 1:channels, pd = 1:pooled_depth, ph = 1:pooled_height, pw = 1:pooled_width
    dstart = (pd - 1) * stride_d - pad_d
    hstart = (ph - 1) * stride_h - pad_h
    wstart = (pw - 1) * stride_w - pad_w

    dend   = min(dstart + kernel_d, depth)
    hend   = min(hstart + kernel_h, height)
    wend   = min(wstart + kernel_w, width)

    dstart = max(dstart, 0) + 1
    hstart = max(hstart, 0) + 1
    wstart = max(wstart, 0) + 1

    y[pw, ph, pd, c, n] =
    sum(x[wstart:wend, hstart:hend, dstart:dend, c, n]) / kernel_size
  end
end

function meanpool3d!(y::AbstractArray{T,5}, x::AbstractArray{T,5};
                     window::Dims{3}=(2,2), padding::Dims{3}=(0,0),
                     stride::Dims{3}=window) where T
    Wx,Hx,Dx,Cx,Nx = size(x)
    Wy,Hy,Dy,Cy,Ny = size(y)
    (w1,w2,w3) = psize(window, x)
    (p1,p2,p3) = psize(padding, x)
    (s1,s2,s3) = psize(stride, x)
    mean_pooling3d_fwd!(x,y,Wx,Hx,Dx,Cx,Nx,Wy,Hy,Dy,w1,w2,w3,p1,p2,p3,s1,s2,s3)
    return y
end

function mean_pooling3d_bwd!{T}(grad_input::AbstractArray{T,5}, grad_output::AbstractArray{T,5},
  width::Int, height::Int, depth::Int, channels::Int, num::Int, pooled_width::Int,
  pooled_height::Int, pooled_depth::Int, kernel_w::Int, kernel_h::Int, kernel_d::Int,
  pad_w::Int, pad_h::Int, pad_d::Int, stride_w::Int, stride_h::Int, stride_d::Int)

  kernel_size = kernel_w * kernel_h * kernel_d
  fill!(grad_input, 0.0)

  #pragma omp parallel for
  for n = 1:num, c = 1:channels, pd = 1:pooled_depth, ph = 1:pooled_height, pw = 1:pooled_width
    dstart = (pd - 1) * stride_d - pad_d
    hstart = (ph - 1) * stride_h - pad_h
    wstart = (pw - 1) * stride_w - pad_w
    dend   = min(dstart + kernel_d, depth)
    hend   = min(hstart + kernel_h, height)
    wend   = min(wstart + kernel_w, width)
    dstart = max(dstart, 0) + 1
    hstart = max(hstart, 0) + 1
    wstart = max(wstart, 0) + 1

    grad_input[wstart:wend, hstart:hend, dstart:dend, c, n] += grad_output[pw, ph, pd, c, n] / kernel_size
  end
end

function meanpool3d_grad!(dx::AbstractArray{T,5}, dy::AbstractArray{T,5}, y::AbstractArray{T,5}, x::AbstractArray{T,5};
                         window::Dims{3}=(2,2,2), padding::Dims{3}=(0,0,0),
                         stride::Dims{3}=window) where T
    Wx,Hx,Dx,Cx,Nx = size(x)
    Wy,Hy,Dy,Cy,Ny = size(y)
    (w1,w2,w3) = psize(window, x)
    (p1,p2,p3) = psize(padding, x)
    (s1,s2,s3) = psize(stride, x)
    mean_pooling3d_bwd!(dx,dy,Wx,Hx,Dx,Cx,Nx,Wy,Hy,Dy,w1,w2,w3,p1,p2,p3,s1,s2,s3)
    return dx
end


text:
include("../Packages/CatalecticantMatrices.jl")

### Create the model
Σ = Catalecticant(2, 3)

### Compute the ML-degree of Σ
ml_degree_witness(Σ; dual=true)
# MLDegreeWitness:
#  • ML degree → 36
#  • model dimension → 15
#  • dual → true

max_stars_repo_name: LauraBMo/InvertingCats

max_issues_repo_name: LauraBMo/InvertingCats

max_forks_repo_name: LauraBMo/InvertingCats


text:
using Test;

struct Rule
    Low::Int
    High::Int
end

function parse_rule(r)
    s = split(r, "-")
    low = parse(Int, s[1])
    high = parse(Int, s[2])

    return Rule(low, high)
end

function valid_for_some_rule(num, rule_lookup)
    for k in keys(rule_lookup)
        for r in rule_lookup[k]
            if r.Low <= num && num <= r.High
                return true
            end
        end
    end

    return false
end

function run(fname)
    f = open(fname, "r")
    lines = readlines(f)
    close(f)

    section = "rules"
    rule_lookup = Dict()
    total = 0
    for l in lines
        # check if we changed section in the doc..
        if l == ""
            continue
        elseif l == "your ticket:"
            section = ""
            section = "your_ticket"
            continue
        elseif l == "nearby tickets:"
            section = "nearby_tickets"
            continue
        end

        # in a section currently.. handle the line
        if section == "rules"
            parts = split(l, ": ")
            rule_name = parts[1]
            ranges = split(parts[2], " or ")
            rules = map(parse_rule, ranges)
            rule_lookup[rule_name] = rules
        elseif section == "nearby_tickets"
            nums = map((x) -> parse(Int, x), split(l, ","))
            for n in nums
                if ! valid_for_some_rule(n, rule_lookup)
                    println("NOT VALID: ", n)
                    total += n
                end
            end
        end
    end

    # println("Rules:", rule_lookup)
    println("[Part 1] File = ", fname, " ... result = ", total)
    return total
end

function compute_rule_positions(tickets, rule_lookup)
    # to begin, each position could be any of the rules
    possibles = []
    for i in 1:length(tickets[1])
        p = Set()
        for k in keys(rule_lookup)
            push!(p, k)
        end
        push!(possibles, p)
    end

    for t in tickets
        for (idx, val) in enumerate(t)
            # check if it fails to satisfy any of the remaining possible rules
            to_delete = []
            for rule_name in possibles[idx]
                rules = rule_lookup[rule_name]
                invalid_count = 0
                for r in rules
                    if !(r.Low <= val && val <= r.High)
                        invalid_count += 1
                    end
                end
                if invalid_count == 2
                    push!(to_delete, rule_name)
                end
            end
            for rule_name in to_delete
                delete!(possibles[idx], rule_name)
            end
        end
    end

    println("Possibles", possibles)

    # Simplify based on things which have only one option
    progress = true
    done = Set()
    while progress
        progress = false

        for (idx, p) in enumerate(possibles)
            if length(p)  == 1 && ! in(idx, done)
                # record the value in this set
                val = Nothing
                for v in p
                    val = v
                end

                # remove this item from other sets
                for i in 1:length(possibles)
                    if i != idx
                        delete!(possibles[i], val)
                    end
                end

                progress = true
                push!(done, idx)
                break
            end
        end
    end

    out = []
    for p in possibles
        val = Nothing
        for v in p
            push!(out, v)
            # assume just one, else we messed up
            break
        end
    end

    return out
end

function run2(fname)
    f = open(fname, "r")
    lines = readlines(f)
    close(f)

    section = "rules"
    rule_lookup = Dict()
    valid_tickets = []
    your_ticket = Nothing # not defined yet
    for l in lines
        # check if we changed section in the doc..
        if l == ""
            continue
        elseif l == "your ticket:"
            section = ""
            section = "your_ticket"
            continue
        elseif l == "nearby tickets:"
            section = "nearby_tickets"
            continue
        end

        # in a section currently.. handle the line
        if section == "rules"
            parts = split(l, ": ")
            rule_name = parts[1]
            ranges = split(parts[2], " or ")
            rules = map(parse_rule, ranges)
            rule_lookup[rule_name] = rules
        elseif section == "your_ticket"
            nums = map((x) -> parse(Int, x), split(l, ","))
            your_ticket = nums
        elseif section == "nearby_tickets"
            nums = map((x) -> parse(Int, x), split(l, ","))
            valid = true
            for n in nums
                if ! valid_for_some_rule(n, rule_lookup)
                    valid = false
                    break
                end
            end
            if valid
                push!(valid_tickets, nums)
            end
        end
    end

    println("Rules:", rule_lookup)
    println("Your ticket:", your_ticket)
    println("Valid tickets:", valid_tickets)

    rule_list = compute_rule_positions(valid_tickets, rule_lookup)
    println("rule_list = ", rule_list)

    result = 1
    for (idx, r) in enumerate(rule_list)
        if startswith(r, "departure")
            val = your_ticket[idx]
            println(r, " ", val)
            result *=  val
        end
    end

    println("[Part 2] File = ", fname, " ... result = ", result)
    return result
end

# @test run("16ex.txt") == 71
# run("16.txt")

run2("16ex.part2.txt")

run2("16.txt")


text:
function ueberpruefeStartpunkt(SPX1,SPX2,PX1,PX2)

if round(SPX1;digits=2) == round(PX1[1,2];digits = 2)
    if round(SPX2;digits=2) == round(PX2[1,2];digits = 2)
        println("Ihr habt die korrekten Koordinaten für den Startpunkt ausgewählt!")
    else
        println("Überprüft erneut eure Eingabe für die X2 Koordinate des Startpunkts.")
    end
elseif round(SPX1;digits = 2) == round(PX1[2,1];digits = 2)
    if round(SPX2;digits = 2) == round(PX2[2,1];digits = 2)
        println("Eure Eingaben stimmen noch nicht ganz. Die erste Zahl in der Klammer steht für die Zeile und die zweite für die Spalte.")
    else
        println("Eure Eingabe für X1 stimmt noch nicht ganz. Die erste Zahl in der Klammer steht für die Zeile und die zweite für die Spalte.")
    end
else
    println("Überprüft erneut eure Eingaben für beide Koordinaten des Startpunkts.")
end
end


text:
import DataFrames
import PredictMD
import Test

a = DataFrames.DataFrame()
a[:x] = ["foo", "foo", "foo", "bar", "bar", "bar"]
contrasts = PredictMD.DataFrameFeatureContrasts(a, [:x])
transformer = PredictMD.MutableDataFrame2DecisionTreeTransformer([:x], :y)
PredictMD.set_feature_contrasts!(transformer, contrasts)

Test.@test length(PredictMD.transform(transformer, a)) == 6

b = DataFrames.DataFrame()
b[:x] = ["bar", "foo"]
Test.@test length(PredictMD.transform(transformer, b)) == 2

c = DataFrames.DataFrame()
c[:x] = ["foo", "bar", "baz"]
Test.@test_throws KeyError PredictMD.transform(transformer, c)

d = DataFrames.DataFrame()
Test.@test_throws ErrorException PredictMD.DataFrameFeatureContrasts(d, [:x, :x])


text:
using Paxos.Ballots
using Paxos.Ledgers
using Paxos.Nodes
using Paxos.Transports.TCP
using Paxos.Transports.Memory

using UUIDs

function ledgerAddEntry(ledger=Ledger())
  req = Request(RequestID(uuid4(),NodeID(),SequenceNumber(1)), Operation(:inc))
  leaderID = NodeID()
  addEntry(ledger, leaderID, req)
  ledger
end


text:
abstract Node

type NodeGene <: Node
    #= A node gene encodes the basic artificial neuron model.
    nodetype should be "INPUT", "HIDDEN", or "OUTPUT" =#
    id::Int64
    ntype::Symbol
    bias::Float64
    response::Float64
    activation::Symbol
    timeConstant::Float64
    function NodeGene(id::Int64, nodetype::Symbol, bias::Float64=0., response::Float64=1., # 4.924273,
                      activation::Symbol=:sigm, timeConstant::Float64=1.0)
        @assert  activation in [:none, :sigm, :tanh, :relu]
        new(id, nodetype, bias, response, activation, timeConstant)
    end
end

function Base.show(io::IO, ng::NodeGene)
    @printf(io, "Node %2d %6s, bias=%+2.10s, response=%+2.10s, %s(), time constant=%+2.5s\n",
            ng.id, ng.ntype, ng.bias, ng.response,ng.activation,ng.timeConstant)
    return
end

function get_child(ng::NodeGene, other::NodeGene)
    # Creates a new NodeGene ramdonly inheriting its attributes from parents
    assert(ng.id == other.id)
    ng = NodeGene(ng.id, ng.ntype,
                    randbool()? ng.bias : other.bias,
                    randbool()? ng.response : other.response, ng.activation,
                    randbool()? ng.timeConstant : other.timeConstant)
    return ng
end

function mutate_bias!(ng::NodeGene, cf::Config)
    ng.bias += randn() * cf.bias_mutation_power
    if (ng.bias > cf.max_weight)
        ng.bias = cf.max_weight
    elseif (ng.bias < cf.min_weight)
        ng.bias = cf.min_weight
    end
end

function mutate_response!(ng::NodeGene, cf::Config)
    #  Mutates the neuron's average firing response.
    ng.response += randn() * cf.bias_mutation_power
end

function mutate_time_constant(ng::NodeGene, cf::Config)
    # Warning: pertubing the time constant (tau) may result in numerical instability
    ng.timeConstant += randn() * .001
    if ng.timeConstant > cf.max_weight
        ng.timeConstant = cf.max_weight
    elseif ng.timeConstant < cf.min_weight
        ng.timeConstant = cf.min_weight
    end
end

function mutate!(ng::NodeGene, cf::Config)
    if rand() < cf.prob_mutatebias mutate_bias!(ng, cf) end
    if rand() < cf.prob_mutatebias mutate_response!(ng, cf) end
#     if rand() < 0.1 ng.mutate_time_constant() end
end

get_new_innov_number(g::Global) =  g.innov_number += 1

type ConnectionGene
    inId::Int64
    outId::Int64
    weight::Float64
    enable::Bool
    key::(Int64,Int64)
    innovNumber::Int64
    function ConnectionGene(g::Global,inId::Int64, outId::Int64, weight::Float64, enable::Bool=true, innov::Int64=0)
        key = (inId, outId)
        if innov == 0
            if haskey(g.innovations,key)
                innovNumber = g.innovations[key]
            else
                innovNumber = get_new_innov_number(g)
                g.innovations[key] = innovNumber
            end
        else
            innovNumber = innov
        end
#         println("innovNumber = $innovNumber")
        new(inId, outId, weight, enable, key, innovNumber)
    end
end

function mutate!(cg::ConnectionGene, cf::Config)
    if rand() < cf.prob_mutate_weight
        mutate_weight!(cg, cf)
    end
    if rand() < cf.prob_togglelink
        cg.enable = true
    end
end

function mutate_weight!(cg::ConnectionGene, cf::Config)
#     cg.weight += (rand() * 2 -1) * cf.weight_mutation_power
    cg.weight += randn() * cf.weight_mutation_power
    if cg.weight > cf.max_weight
        cg.weight = cf.max_weight
    elseif cg.weight < cf.min_weight
        cg.weight = cf.min_weight
    end
end

function weight_replaced!(cg::ConnectionGene, cf::Config)
        # cg.weight = random.uniform(-Config.random_range, Config.random_range)
        cg.weight = randn() * cf.weight_stdev
end

function split(g::Global,cg::ConnectionGene, node_id::Int64)
    # Splits a connection, creating two new connections and disabling this one """
    cg.enable = false
    new_conn1 = ConnectionGene(g, cg.inId, node_id, 1.0, true)
    new_conn2 = ConnectionGene(g, node_id, cg.outId, cg.weight, true)
    return new_conn1, new_conn2
end

is_same_innov(self::ConnectionGene, other::ConnectionGene) = self.innovNumber == other.innovNumber

get_child(self::ConnectionGene, other::ConnectionGene) = randbool() ? self : other

function maxInnov(cgs::Dict{(Int64,Int64), ConnectionGene})
    @assert length(cgs) > 0

    cgsKeys = collect(keys(cgs))
    maxCg = cgs[cgsKeys[1]]
    for key in cgsKeys
        maxCg = cgs[key].innovNumber > maxCg.innovNumber? cgs[key]: maxCg
    end
    return maxCg
end

function Base.show(io::IO, cg::ConnectionGene)
    @printf(io, "In: %2d, Out: %2d, Weight: %+3.5f, %6s, InnovID: %d",
            cg.inId, cg.outId, cg.weight,(cg.enable? "Enabled":"Disabled"), cg.innovNumber)
    return
end

max_stars_repo_name: Andy-P/NEAT

max_issues_repo_name: Andy-P/NEAT

max_forks_repo_name: Andy-P/NEAT


text:
using Test
using DebuggerFramework

# This file will define a minimalistic programming language and implement a
# debugger for this language. The language is not supposed to be useful, it
# is a proving ground for the DebuggerFramework APIs.

# Language semantics:
# Every line has to start with its line number relative to the current program,
# followed by a colon. The colon is followed by a whitespace separated list of
# integers, indicating either a line to call, or 0 which indicates a return.
# This language has no conditional branches, so all loops are infinite. The
# computed value is the number of calls that were made

module BrainStack
    using DebuggerFramework
    struct BrainStackAST
        stmts::Vector{Vector{Int}}
    end

    function Base.parse(::Type{BrainStackAST}, data)
        BrainStackAST(map(enumerate(split(data, '\n', keepempty=false))) do args
            lineno, line = args
            parts = split(line, r"\s", keepempty=false)
            @assert parts[1] == "$lineno:"
            map(x->parse(Int,x), parts[2:end])
        end)
    end

    mutable struct BrainStackInterpreterState
        ast::BrainStackAST
        stack::Vector{Tuple{Int, Int}}
        idx::Int
        pos::Int
        ncalls::Int
    end

    # Helpfully split up the interpreter to be re-usable from the debugger below
    enter(ast::BrainStackAST) = BrainStackInterpreterState(ast, [(0,0)], 1, 1, 0)
    function step_one!(state::BrainStackInterpreterState)
        if state.idx == 0
            return false
        end
        op = state.ast.stmts[state.idx][state.pos]
        if op == 0
            state.idx, state.pos = pop!(state.stack)
            state.pos += 1
        else
            push!(state.stack, (state.idx, state.pos))
            state.idx, state.pos = op, 1
            state.ncalls += 1
        end
        return true
    end
    function interpret(ast::BrainStackAST)
        state = enter(ast)
        while step_one!(state); end
        state.ncalls
    end

    # Debugger support
    struct BrainStackStackRef <: DebuggerFramework.StackFrame
        state::BrainStackInterpreterState
        idx::Int
    end

    function DebuggerFramework.debug(ast::BrainStackAST, args...)
        state = enter(ast)
        stack = [BrainStackStackRef(state, 0)]
        DebuggerFramework.RunDebugger(stack, args...)
    end

    function idxpos(frame)
        if frame.idx == 0
            idx, pos = frame.state.idx, frame.state.pos
        else
            idx, pos = frame.state.stack[end-(frame.idx-1)]
        end
        idx, pos
    end

    function DebuggerFramework.print_status_synthtic(io::IO, state, frame::BrainStackStackRef, lines_before, total_lines)
        ast = frame.state.ast
        idx, pos = idxpos(frame)
        print(io, "$idx: ")
        for (opno, op) in enumerate(ast.stmts[idx])
            if opno == pos
                printstyled(io, string(op); bold = true, color=:yellow)
            else
                print(io, op)
            end
            (opno != length(ast.stmts[idx])) && print(io, ' ')
       end
       return 1
    end

    DebuggerFramework.locdesc(frame::BrainStackStackRef) = "statement $(idxpos(frame)[1])"

    function update_stack!(ds, state, stack)
        ds.stack = [BrainStackStackRef(state, i) for i = 0:length(stack)-1]
        if ds.level > length(ds.stack)
            ds.level = length(ds.stack)
        end
    end

    function DebuggerFramework.execute_command(state, frame::BrainStackStackRef, cmd::Val{:si}, command)
        step_one!(frame.state)
        update_stack!(state, frame.state, frame.state.stack)
        return true
    end
end


text:
@extend struct VanillaOptions <: Options
    e₀::Float64 #Industrial emissions 2010 (GtCO2 per year)
    μ₀::Float64 #Initial emissions control rate for base case 2010
    tnopol::Float64 #Period before which no emissions controls base
    cprice₀::Float64 #Initial base carbon price (2005$ per tCO2)
    gcprice::Float64 #Growth rate of base carbon price per year
    periodfullpart::Float64 #Period at which have full participation
    partfract2010::Float64 #Fraction of emissions under control in 2010
    partfractfull::Float64 #Fraction of emissions under control at full time
end

function options(version::V2013R{VanillaFlavour};
    N::Int = 60, #Number of years to calculate (from 2010 onwards)
    tstep::Int = 5, #Years per Period
    α::Float64 = 1.45, #Elasticity of marginal utility of consumption
    ρ::Float64 = 0.015, #Initial rate of social time preference per year ρ
    γₑ::Float64 = 0.3, #Capital elasticity in production function
    pop₀::Int = 6838, #Initial world population (millions)
    popadj::Float64 = 0.134, #Growth rate to calibrate to 2050 pop projection
    popasym::Int = 10500, #Asymptotic population (millions)
    δk::Float64 = 0.1, #Depreciation rate on capital (per year)
    q₀::Float64 = 63.69, #Initial world gross output (trill 2005 USD)
    k₀::Float64 = 135.0, #Initial capital value (trill 2005 USD)
    a₀::Float64 = 3.8, #Initial level of total factor productivity
    ga₀::Float64 = 0.079, #Initial growth rate for TFP per 5 years
    δₐ::Float64 = 0.006, #Decline rate of TFP per 5 years
    gσ₁::Float64 = -0.01, #Initial growth of sigma (continuous per year)
    δσ::Float64 = -0.001, #Decline rate of decarbonization per period
    eland₀::Float64 = 3.3, #Carbon emissions from land 2010 (GtCO2 per year)
    deland::Float64 = 0.2, #Decline rate of land emissions (per period)
    e₀::Float64 = 33.61, #Industrial emissions 2010 (GtCO2 per year)
    μ₀::Float64 = 0.039, #Initial emissions control rate for base case 2010
    mat₀::Float64 = 830.4, #Initial Concentration in atmosphere 2010 (GtC)
    mu₀::Float64 = 1527.0, #Initial Concentration in upper strata 2010 (GtC)
    ml₀::Float64 = 10010.0, #Initial Concentration in lower strata 2010 (GtC)
    mateq::Float64 = 588.0, #Equilibrium concentration atmosphere  (GtC)
    mueq::Float64 = 1350.0, #Equilibrium concentration in upper strata (GtC)
    mleq::Float64 = 10000.0, #Equilibrium concentration in lower strata (GtC)
    ϕ₁₂::Float64 = 0.088, #Carbon cycle transition matrix coefficient
    ϕ₂₃::Float64 = 0.0025, #Carbon cycle transition matrix coefficient
    t2xco2::Float64 = 2.9, #Equilibrium temp impact (oC per doubling CO2)
    fₑₓ0::Float64 = 0.25, #2010 forcings of non-CO2 GHG (Wm-2)
    fₑₓ1::Float64 = 0.7, #2100 forcings of non-CO2 GHG (Wm-2)
    tocean₀::Float64 = 0.0068, #Initial lower stratum temp change (C from 1900)
    tatm₀::Float64 = 0.8, #Initial atmospheric temp change (C from 1900)
    tatm_up::Float64 = 40.0, #Maximum atmospheric temp change (C from 1900)
    ξ₁::Float64 = 0.098, #Climate equation coefficient for upper level
    ξ₃::Float64 = 0.088, #Transfer coefficient upper to lower stratum
    ξ₄::Float64 = 0.025, #Transfer coefficient for lower level
    η::Float64 = 3.8, #Forcings of equilibrium CO2 doubling (Wm-2)
    ψ₁::Float64 = 0.0, #Damage intercept
    ψ₂::Float64 = 0.00267, #Damage quadratic term
    ψ₃::Float64 = 2.0, #Damage exponent
    θ₂::Float64 = 2.8, #Exponent of control cost function
    pback::Float64 = 344.0, #Cost of backstop 2005$ per tCO2 2010
    gback::Float64 = 0.025, #Initial cost decline backstop cost per period
    limμ1::Float64 = 1.0, #Upper limit on control rate before 2150
    limμ2::Float64 = 1.2, #Upper limit on control rate after 2150
    tnopol::Float64 = 45.0, #Period before which no emissions controls base
    cprice₀::Float64 = 1.0, #Initial base carbon price (2005$ per tCO2)
    gcprice::Float64 = 0.02, #Growth rate of base carbon price per year
    fosslim::Float64 = 6000.0, #Maximum cumulative extraction fossil fuels (GtC)
    scale1::Float64 = 0.016408662, #Multiplicative scaling coefficient
    scale2::Float64 = -3855.106895, #Additive scaling coefficient
    periodfullpart::Float64 = 21.0, #Period at which have full participation
    partfract2010::Float64 = 1.0, #Fraction of emissions under control in 2010
    partfractfull::Float64 = 1.0) #Fraction of emissions under control at full time
    VanillaOptions(N,tstep,α,ρ,γₑ,pop₀,popadj,popasym,δk,q₀,k₀,a₀,ga₀,δₐ,gσ₁,δσ,eland₀,deland,mat₀,mu₀,ml₀,mateq,mueq,mleq,ϕ₁₂,ϕ₂₃,t2xco2,fₑₓ0,fₑₓ1,tocean₀,tatm₀,tatm_up,ξ₁,ξ₃,ξ₄,η,ψ₁,ψ₂,ψ₃,θ₂,pback,gback,limμ1,limµ2,fosslim,scale1,scale2,e₀,μ₀,tnopol,cprice₀,gcprice,periodfullpart,partfract2010,partfractfull)
end

function Base.show(io::IO, ::MIME"text/plain", opt::VanillaOptions)
    println(io, "Options for Vanilla 2013R version");
    println(io, "Time step");
    println(io, "N: $(opt.N), tstep: $(opt.tstep)");
    println(io, "Preferences");
    println(io, "α: $(opt.α), ρ: $(opt.ρ)");
    println(io, "Population and Technology");
    println(io, "γₑ: $(opt.γₑ), pop₀: $(opt.pop₀), popadj: $(opt.popadj), popasym: $(opt.popasym), δk: $(opt.δk)");
    println(io, "q₀: $(opt.q₀), k₀: $(opt.k₀), a₀: $(opt.a₀), ga₀: $(opt.ga₀), δₐ: $(opt.δₐ)");
    println(io, "Emissions Parameters");
    println(io, "gσ₁: $(opt.gσ₁), δσ: $(opt.δσ), eland₀: $(opt.eland₀)");
    println(io, "deland: $(opt.deland), e₀: $(opt.e₀), μ₀: $(opt.μ₀)");
    println(io, "Carbon Cycle");
    println(io, "mat₀: $(opt.mat₀), mu₀: $(opt.mu₀), ml₀: $(opt.ml₀)");
    println(io, "mateq: $(opt.mateq), mueq: $(opt.mueq), mleq: $(opt.mleq)");
    println(io, "Flow Parameters");
    println(io, "ϕ₁₂: $(opt.ϕ₁₂), ϕ₂₃: $(opt.ϕ₂₃)");
    println(io, "Climate Model Parameters");
    println(io, "t2xco2: $(opt.t2xco2), fₑₓ0: $(opt.fₑₓ0), fₑₓ1: $(opt.fₑₓ1)");
    println(io, "tocean₀: $(opt.tocean₀), tatm₀: $(opt.tatm₀), tatm_up: $(opt.tatm_up), ξ₁: $(opt.ξ₁)");
    println(io, "ξ₃: $(opt.ξ₃), ξ₄: $(opt.ξ₄), η: $(opt.η)");
    println(io, "Climate Damage Parameters");
    println(io, "ψ₁: $(opt.ψ₁), ψ₂: $(opt.ψ₂), ψ₃: $(opt.ψ₃)");
    println(io, "Abatement Cost");
    println(io, "θ₂: $(opt.θ₂), pback: $(opt.pback), gback: $(opt.gback), limμ1: $(opt.limμ1), limμ2: $(opt.limμ2)");
    println(io, "tnopol: $(opt.tnopol), cprice₀: $(opt.cprice₀), gcprice: $(opt.gcprice)");
    println(io, "Participation Parameters");
    println(io, "periodfullpart: $(opt.periodfullpart), partfract2010: $(opt.partfract2010), partfractfull: $(opt.partfractfull)");
    println(io, "Fossil Fuel Availability");
    println(io, "fosslim: $(opt.fosslim)");
    println(io, "Scaling Parameters");
    print(io, "scale1: $(opt.scale1), scale2: $(opt.scale2)");
end

@extend struct VanillaParameters <: Parameters
    pbacktime::Array{Float64,1} # Backstop price
    cpricebase::Array{Float64,1} # Carbon price in base case
    rr::Array{Float64,1} # Average utility social discount rate
    optlrsav::Float64 # Optimal savings rate
    partfract::Array{Float64,1} # Fraction of emissions in control regime
end

function generate_parameters(c::VanillaOptions)
    optlrsav::Float64 = (c.δk + .004)/(c.δk + .004c.α + c.ρ)*c.γₑ; # Optimal savings rate
    ϕ₁₁::Float64 = 1 - c.ϕ₁₂; # Carbon cycle transition matrix coefficient
    ϕ₂₁::Float64 = c.ϕ₁₂*c.mateq/c.mueq; # Carbon cycle transition matrix coefficient
    ϕ₂₂::Float64 = 1 - ϕ₂₁ - c.ϕ₂₃; # Carbon cycle transition matrix coefficient
    ϕ₃₂::Float64 = c.ϕ₂₃*c.mueq/c.mleq; # Carbon cycle transition matrix coefficient
    ϕ₃₃::Float64 = 1 - ϕ₃₂; # Carbon cycle transition matrix coefficient
    σ₀::Float64 = c.e₀/(c.q₀*(1-c.μ₀)); # Carbon intensity 2010 (kgCO2 per output 2005 USD 2010)
    λ::Float64 = c.η/c.t2xco2; # Climate model parameter

    # Backstop price
    pbacktime = Array{Float64}(undef, c.N);
    # Growth rate of productivity from 0 to N
    gₐ = Array{Float64}(undef, c.N);
    # Emissions from deforestation
    Etree = Array{Float64}(undef, c.N);
    # Average utility social discount rate
    rr = Array{Float64}(undef, c.N);
    # Carbon price in base case
    cpricebase = Array{Float64}(undef, c.N);

    for i in 1:c.N
        pbacktime[i] = c.pback*(1-c.gback)^(i-1);
        gₐ[i] = c.ga₀*exp.(-c.δₐ*5*(i-1));
        Etree[i] = c.eland₀*(1-c.deland)^(i-1);
        rr[i] = 1/((1+c.ρ)^(c.tstep*(i-1)));
        cpricebase[i] = c.cprice₀*(1+c.gcprice)^(5*(i-1));
    end

    # Initial conditions and offset required
    # Level of population and labor
    L = Array{Float64}(undef, c.N);
    L[1] = c.pop₀;
    # Level of total factor productivity
    A = Array{Float64}(undef, c.N);
    A[1] = c.a₀;
    # Change in sigma (cumulative improvement of energy efficiency)
    gσ = Array{Float64}(undef, c.N);
    gσ[1] = c.gσ₁;
    # CO2-equivalent-emissions output ratio
    σ = Array{Float64}(undef, c.N);
    σ[1] = σ₀;

    for i in 1:c.N-1
        L[i+1] = L[i]*(c.popasym/L[i])^c.popadj;
        A[i+1] = A[i]/(1-gₐ[i]);
        gσ[i+1] = gσ[i]*((1+c.δσ)^c.tstep);
        σ[i+1] = σ[i]*exp.(gσ[i]*c.tstep);
    end

    # Adjusted cost for backstop
    θ₁ = Array{Float64}(undef, c.N);
    # Exogenous forcing for other greenhouse gases
    fₑₓ = Array{Float64}(undef, c.N);
    # Fraction of emissions in control regime
    partfract = Array{Float64}(undef, c.N);

    for i in 1:c.N
        θ₁[i] = pbacktime[i]*σ[i]/c.θ₂/1000.0;
        fₑₓ[i] = if i <= 19
                         c.fₑₓ0+(1/18)*(c.fₑₓ1-c.fₑₓ0)*(i-1)
                     else
                         c.fₑₓ0+(c.fₑₓ1-c.fₑₓ0)
                     end;
        partfract[i] = if i <= c.periodfullpart
                            c.partfract2010+(c.partfractfull-c.partfract2010)*(i-1)/c.periodfullpart
                       else
                            c.partfractfull
                       end;
    end
    partfract[1] = c.partfract2010;
    VanillaParameters(ϕ₁₁,ϕ₂₁,ϕ₂₂,ϕ₃₂,ϕ₃₃,σ₀,λ,gₐ,Etree,L,A,gσ,σ,θ₁,fₑₓ,pbacktime,cpricebase,rr,optlrsav,partfract)
end

function Base.show(io::IO, opt::VanillaParameters)
    ctx = IOContext(io, :limit => true, :compact => true);
    println(io, "Calculated Parameters for Vanilla 2013R");
    println(io, "Optimal savings rate: $(opt.optlrsav)");
    println(io, "Carbon cycle transition matrix coefficients");
    println(io, "ϕ₁₁: $(opt.ϕ₁₁), ϕ₂₁: $(opt.ϕ₂₁), ϕ₂₂: $(opt.ϕ₂₂), ϕ₃₂: $(opt.ϕ₃₂), ϕ₃₃: $(opt.ϕ₃₃)");
    println(io, "2010 Carbon intensity: $(opt.σ₀)");
    println(io, "Climate model parameter: $(opt.λ)");
    println(io, "Backstop price: ");
    show(ctx, opt.pbacktime);
    println(io, "");
    print(io, "Growth rate of productivity: ");
    show(ctx, opt.gₐ);
    println(io, "");
    print(io, "Emissions from deforestation: ");
    show(ctx, opt.Etree);
    println(io, "");
    print(io, "Avg utility social discout rate: ");
    show(ctx, opt.rr);
    println(io, "");
    print(io, "Base case carbon price: ");
    show(ctx, opt.cpricebase);
    println(io, "");
    print(io, "Population and labour: ");
    show(ctx, opt.L);
    println(io, "");
    print(io, "Total factor productivity: ");
    show(ctx, opt.A);
    println(io, "");
    print(io, "Δσ: ");
    show(ctx, opt.gσ);
    println(io, "");
    print(io, "σ: ");
    show(ctx, opt.σ);
    println(io, "");
    print(io, "θ₁: ");
    show(ctx, opt.θ₁);
    println(io, "");
    print(io, "Exogenious forcing: ");
    show(ctx, opt.fₑₓ);
    println(io, "");
    print(io, "Fraction of emissions in control regieme: ");
    show(ctx, opt.partfract);
end

@extend struct VanillaEquations <: Equations
    cc::Array{ConstraintRef{Model,C,Shape} where Shape<:JuMP.AbstractShape where C,1} # Output Consumption
end

function model_eqs(scenario::Scenario, model::Model, config::VanillaOptions, params::VanillaParameters, vars::Variables, isMumps::Bool)
    N = config.N;
    # Equations #
    # Emissions Equation
    eeq = @constraint(model, [i=1:N], vars.E[i] == vars.Eind[i] + params.Etree[i]);
    # Industrial Emissions
    @constraint(model, [i=1:N], vars.Eind[i] == params.σ[i] * vars.YGROSS[i] * (1-vars.μ[i]));
    # Radiative forcing equation
    @NLconstraint(model, [i=1:N], vars.FORC[i] == config.η * (log(vars.Mₐₜ[i]/588.0)/log(2)) + params.fₑₓ[i]);
    # Equation for damage fraction
    @NLconstraint(model, [i=1:N], vars.Ω[i] == config.ψ₁*vars.Tₐₜ[i]+config.ψ₂*vars.Tₐₜ[i]^config.ψ₃);
    # Damage equation
    @constraint(model, [i=1:N], vars.DAMAGES[i] == vars.YGROSS[i]*vars.Ω[i]);
    # Cost of exissions reductions equation
    @NLconstraint(model, [i=1:N], vars.Λ[i] == vars.YGROSS[i] * params.θ₁[i] * vars.μ[i]^config.θ₂ * params.partfract[i]^(1-config.θ₂));
    # Equation for MC abatement
    @NLconstraint(model, [i=1:N], vars.MCABATE[i] == params.pbacktime[i] * vars.μ[i]^(config.θ₂-1));
    # Carbon price equation from abatement
    @NLconstraint(model, [i=1:N], vars.CPRICE[i] == params.pbacktime[i] * (vars.μ[i]/params.partfract[i])^(config.θ₂-1));
    # Output gross equation
    @NLconstraint(model, [i=1:N], vars.YGROSS[i] == params.A[i]*(params.L[i]/1000.0)^(1-config.γₑ)*vars.K[i]^config.γₑ);
    # Output net of damages equation
    @constraint(model, [i=1:N], vars.YNET[i] == vars.YGROSS[i]*(1-vars.Ω[i]));
    # Output net equation
    @constraint(model, [i=1:N], vars.Y[i] == vars.YNET[i] - vars.Λ[i]);
    # Consumption equation
    cc = @constraint(model, [i=1:N], vars.C[i] == vars.Y[i] - vars.I[i]);
    # Per capita consumption definition
    @constraint(model, [i=1:N], vars.CPC[i] == 1000.0 * vars.C[i] / params.L[i]);
    # Savings rate equation
    @constraint(model, [i=1:N], vars.I[i] == vars.S[i] * vars.Y[i]);
    # Period utility
    @constraint(model, [i=1:N], vars.CEMUTOTPER[i] == vars.PERIODU[i] * params.L[i] * params.rr[i]);
    # Instantaneous utility function equation
    @NLconstraint(model, [i=1:N], vars.PERIODU[i] == ((vars.C[i]*1000.0/params.L[i])^(1-config.α)-1)/(1-config.α)-1);

    # Equations (offset) #
    # Cumulative carbon emissions
    @constraint(model, [i=1:N-1], vars.CCA[i+1] == vars.CCA[i] + vars.Eind[i]*5/3.666);
    # Atmospheric concentration equation
    @constraint(model, [i=1:N-1], vars.Mₐₜ[i+1] == vars.Mₐₜ[i]*params.ϕ₁₁ + vars.Mᵤₚ[i]*params.ϕ₂₁ + vars.E[i]*(5/3.666));
    # Lower ocean concentration
    @constraint(model, [i=1:N-1], vars.Mₗₒ[i+1] == vars.Mₗₒ[i]*params.ϕ₃₃ + vars.Mᵤₚ[i]*config.ϕ₂₃);
    # Shallow ocean concentration
    @constraint(model, [i=1:N-1], vars.Mᵤₚ[i+1] == vars.Mₐₜ[i]*config.ϕ₁₂ + vars.Mᵤₚ[i]*params.ϕ₂₂ + vars.Mₗₒ[i]*params.ϕ₃₂);
    # Temperature-climate equation for atmosphere
    @constraint(model, [i=1:N-1], vars.Tₐₜ[i+1] == vars.Tₐₜ[i] + config.ξ₁*((vars.FORC[i+1]-params.λ*vars.Tₐₜ[i])-(config.ξ₃*(vars.Tₐₜ[i]-vars.Tₗₒ[i]))));
    # Temperature-climate equation for lower oceans
    @constraint(model, [i=1:N-1], vars.Tₗₒ[i+1] == vars.Tₗₒ[i] + config.ξ₄*(vars.Tₐₜ[i]-vars.Tₗₒ[i]));
    # Capital balance equation
    @NLconstraint(model, [i=1:N-1], vars.K[i+1] <= (1-config.δk)^config.tstep * vars.K[i] + config.tstep*vars.I[i]);
    # Interest rate equation
    @NLconstraint(model, [i=1:N-1], vars.RI[i] == (1+config.ρ)*(vars.CPC[i+1]/vars.CPC[i])^(config.α/config.tstep)-1);

    # Savings rate for asympotic equilibrium
    for i=51:N
        JuMP.fix(vars.S[i], params.optlrsav; force=true);
    end
    # Initial conditions
    JuMP.fix(vars.CCA[1], 90.0; force=true);
    JuMP.fix(vars.Mₐₜ[1], config.mat₀; force=true);
    JuMP.fix(vars.Mᵤₚ[1], config.mu₀; force=true);
    JuMP.fix(vars.Mₗₒ[1], config.ml₀; force=true);
    JuMP.fix(vars.Tₗₒ[1], config.tocean₀; force=true);

    if !isMumps || typeof(scenario) <: OptimalPriceScenario
        JuMP.fix(vars.K[1], config.k₀; force=true);
        JuMP.fix(vars.Tₐₜ[1], config.tatm₀; force=true);
    elseif typeof(scenario) <: BasePriceScenario
        # We can't fix these, the solution becomes concave.
        # This is something buggy in JuMP I think. Haven't been able to pin it down.
        @NLconstraint(model, vars.K[1] == config.k₀);
        @constraint(model, vars.Tₐₜ[1] == config.tatm₀);
    end

    @constraint(model, vars.UTILITY == config.tstep * config.scale1 * sum(vars.CEMUTOTPER[i] for i=1:N) + config.scale2);

    # Objective function
    @objective(model, Max, vars.UTILITY);
    VanillaEquations(eeq,cc)
end

include("ScenariosVanilla.jl")

function solve(scenario::Scenario, version::V2013R{VanillaFlavour};
    config::VanillaOptions = options(version),
    linear_solver::ipoptLinearSolver=ma97,
    optimizer::MOI.OptimizerWithAttributes = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 5, "max_iter" => 99900, "print_frequency_iter" => 250, "sb" => "yes", "linear_solver" => selectLinearSolver(linear_solver)))

    model = Model(optimizer);

    # Generate a solver test to implement DICE.jl#35 hacks.
    isMumps = JuMP.get_optimizer_attribute(model, "linear_solver") == "mumps";

    params = generate_parameters(config);

    # Rate limit
    μ_ubound = [if t < 30 config.limμ1 else config.limμ2*params.partfract[t] end for t in 1:config.N];

    cprice_ubound = assign_scenario(scenario, config, params);
    cprice_ubound[1] = params.cpricebase[1];
    # Warning: If parameters are changed, the next equation might make base case infeasible.
    # If so, reduce tnopol so that we don't run out of resources.
    for i in 2:config.N
        if i > config.tnopol
            cprice_ubound[i] = 1000.0;
        end
    end

    variables = model_vars(version, model, config.N, config.fosslim, μ_ubound, config.tatm_up, cprice_ubound);

    equations = model_eqs(scenario, model, config, params, variables, isMumps);

    optimize!(model);

    results = model_results(model, config, params, variables, equations);

    DICENarrative(config,params,model,scenario,version,variables,equations,results)
end


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

function stress_driven_uniaxial_increment!(material, dstress11,
        dt; dstrain=[dstress11/200e3, -0.3*dstress11/200e3, -0.3*dstress11/200e3, 0.0, 0.0, 0.0], max_iter=50, norm_acc=1e-9)
    converged = false
    stress0 = tovoigt(material.variables.stress)
    for i=1:max_iter
        material.ddrivers.time = dt
        material.ddrivers.strain = fromvoigt(SymmetricTensor{2,3,Float64}, dstrain; offdiagscale=2.0)
        integrate_material!(material)
        stress = tovoigt(material.variables_new.stress)
        dstress = stress - stress0
        r = dstress
        r[1] -= dstress11
        D = tovoigt(material.variables_new.jacobian)
        dstr = -D\r
        dstrain .+= dstr
        norm(dstr) < norm_acc && (converged = true; break)
    end
    converged || error("No convergence in strain increment")
    return nothing
end


text:
"""
```
write_forecast_centric_packet(m, input_type, cond_type,
    output_vars = [:forecastobs, :forecastpseudo, :shockdecobs, :shockdecpseudo];
    sections = [:estimation, :forecast]
    forecast_string = "", outdir = joinpath(saveroot(m), \"Packets\", spec(m), subspec(m)))
```

Write standard estimation and forecast result packet to `outdir`, with an emphasis
on the forecasting results. The main difference with the standard packet is that the
prior posterior plots at placed at the very end.
"""
function write_forecast_centric_packet(m::AbstractDSGEModel, input_type::Symbol, cond_type::Symbol,
                               output_vars::Vector{Symbol} = [:forecastobs, :forecastpseudo,
                                                              :shockdecobs, :shockdecpseudo];
                               sections::Vector{Symbol} = [:estimation, :forecast],
                               forecast_string::String = "",
                               outdir::String = joinpath(saveroot(m), "Packets", spec(m), subspec(m)),
                               purpose::String = "")
    @assert issubset(sections, [:estimation, :forecast, :irf]) "Section specified in `section` kwarg is not supported. Must be a subset of [:estimation, :forecast, :irf]."

    # Title and authors
    title = "Estimation and Forecasting Results \\\\ " * DSGE.description(m)
    authors = "User"

    # Compute file name
    cdvt_str = cond_type == :none ? "" : "_cdvt=" * cond_vintage(m)
    fn = joinpath(outdir, "results_cond=" * string(cond_type) * cdvt_str * "_fcid=" * forecast_string * "_vint=" * data_vintage(m) * ".tex")
    isdir(dirname(fn)) || mkpath(dirname(fn))

    # Write packet
    open(fn, "w") do fid
        write_preamble(fid, title, authors)
        write_spec_section(fid, m, purpose = purpose)
        if :estimation in sections
            write_estimation_table_section(fid, m)
        end
        if :forecast in sections
            write_forecast_section(fid, m, input_type, cond_type,
                                   setdiff(output_vars, [:irfstates, :irfobs, :irfpseudo]),
                                   forecast_string = forecast_string)
        end
        if :irf in sections
            write_irf_section(fid, m , input_type, cond_type, output_vars,
                              forecast_string = forecast_string)
        end
        if :estimation in sections
            write_estimation_plot_section(fid, m)
        end
        write_postamble(fid)
    end
    println("Wrote " * fn)

end

"""
```
write_estimation_table_section(fid, m; plotroot = "")
```

Write parameter moment tables.
"""
function write_estimation_table_section(fid::IOStream, m::AbstractDSGEModel;
                                  plotroot::String = "")
    @printf fid "\n\n"
    @printf fid "\\clearpage\n"
    @printf fid "\\section{Estimation Moments}\n"
    @printf fid "\n"

    @printf fid "\\input{%s}\n" tablespath(m, "estimate", "moments.tex")

end

"""
```
write_estimation_plot_section(fid, m; plotroot = "")
```

Write prior/posterior plots.
"""
function write_estimation_plot_section(fid::IOStream, m::AbstractDSGEModel;
                                  plotroot::String = "")

    @printf fid "\\clearpage\n"
    @printf fid "\\section{Estimation Histograms}\n"
    write_estimation_plots(fid, m, plotroot = plotroot)
    @printf fid "\n"

end