text:
import Base.Random: rand, zeros
import Base: ==, hash, show

export Genotype, contribs, fitness

@doc """A genotype representation.
"""
type Genotype{T <: Landscape}
  alleles::AlleleString
  landscape::T
end

Genotype{T <: Landscape}(alleles::Integer, landscape::T) = Genotype(AlleleString(alleles), landscape)

Genotype{T <: Landscape}(g::Genotype{T}) = Genotype(g.alleles, g.landscape)

@doc """contribs(g::Genotype, update::Function)

Return a vector of contributions where the ith element in
the vector is the contribution made by the ith allele in the
genotype, given the values of the k alleles to which it is
epistatically linked.
"""
function contribs(g::Genotype, update::Function)
  return map(1:g.landscape.n) do i
    linksmask = AlleleMask(g.landscape.links[i])
    contribstring = AlleleString(g.alleles & linksmask)
    return get!(update, g.landscape.contribs[i], contribstring)
  end
end

@doc """contribs(g::Genotype)
"""
contribs(g::Genotype{NKLandscape}) = contribs(g, () -> rand())

@doc """contribs(g::Genotype)
"""
contribs(g::Genotype{NKqLandscape}) = contribs(g, () -> rand(0:(g.landscape.q - 1)))

@doc """contribs(g::Genotype)
"""
function contribs(g::Genotype{NKpLandscape})
  update = () -> begin
    if rand() < g.landscape.p
      return 0.0
    else
      return rand()
    end
  end
  return contribs(g, update)
end

# TODO: Do we include the fake NKp zeros in the sum or just let the fitness range lower?
@doc """fitness(g::Genotype{NKqLandscape})
"""
fitness(g::Genotype{NKqLandscape}) = mean(contribs(g)) / (g.landscape.q - 1)

@doc """fitness(g::Genotype)

Compute the fitness of a particular genotype.
"""
fitness(g::Genotype) = mean(contribs(g))

function rand(::Type{Genotype}, ls::Landscape)
  nmask = (AlleleMask(1) << ls.n) - 1
  Genotype(rand(AlleleString) & nmask, ls)
end

zeros(::Type{Genotype}, ls::Landscape) = Genotype(0, ls)

==(left::Genotype, right::Genotype) = left.alleles == right.alleles && is(left.landscape, right.landscape)

hash(g::Genotype) = hash(Genotype) $ hash(g.landscape) $ hash(g.alleles)

show(io::Base.IO, g::Genotype) = print(io, "Genotype($(bits(g.alleles)[(end - g.landscape.n + 1):end]), $(g.landscape))")


text:
#############################################################################
# Joulia
# A Large-Scale Spatial Power System Model for Julia
# See https://github.com/JuliaEnergy/Joulia.jl
#############################################################################
# This file contains functions and structs for power plants

mutable struct PowerPlants
    id
    node
    fuel
    capacity
    mc

    function PowerPlants(pp::DataFrame; avail::DataFrame,
            prices::DataFrame)

        P = pp[1]
        pp_dict = df_to_dict_with_id(pp)
        avail_dict = df_to_dict(avail)
        prices_dict = df_to_dict(prices)

        return new(P, pp_dict[:Node],
                   pp_dict[:Fuel],
                   get_avail_powerplants(pp_dict, avail_dict),
                   calc_mc(pp_dict, prices_dict),
                   )
    end
end

"""
    calc_mc(pp_dict, prices_dict)

Calculate marginal costs from fuel price, efficiency, CO2 price and variable costs.
"""
function calc_mc(pp_dict, prices_dict)

    mc = Dict{String, Array{Float64, 1}}()

    for p in keys(pp_dict[:Fuel])
        f = Symbol(pp_dict[:Fuel][p])
        price = prices_dict[f]
        eff = pp_dict[:Efficiency][p]
        co_price = prices_dict[:CO2]
        emission = pp_dict[:Emission][p]
        vc = pp_dict[:VariableCost][p]

        mc[p] = price / eff + co_price * emission .+ vc
    end

    return mc
end

"""
    get_avail_powerplants(pp_dict, avail)

Determine power plant availability
"""
function get_avail_powerplants(pp_dict, avail)
    P = keys(pp_dict[:Fuel])
    avail_dict = Dict(p => avail[Symbol(pp_dict[:Fuel][p])] for p in P)
    return Dict(p => avail_dict[p] .* pp_dict[:Capacity][p] for p in P)
end

ext: jl


text:
using Minerva

function run_episodes()
    env = CartPole()
    # render(env)
    scores = []
    for i = 1:1000
        _ = reset!(env)
        total_reward = 0
        while true
            a = rand(env.action_space)
            s′, r, done, info = step!(env, a)
            total_reward += r
            # @info step, s′, r, done, info
            done && break
        end
        push!(scores, total_reward)
    end
    return scores
end

# scores = run_episodes()

function sample_run()
    env = CartPole()
    _ = reset!(env)
    render(env)
    total_reward = 0
    steps = 0
    while true
        a = rand(env.action_space)
        s′, r, done, info = step!(env, a)
        total_reward += r
        steps += 1
        # @info step, s′, r, done, info
        done && break
    end
    @info steps
    return total_reward, steps
end

_ = sample_run()


text:
# Shaded line plots -- to be added after implot v0.4

function PlotShaded(x::AbstractArray{T,1}, y1::AbstractArray{T,1}, y2::AbstractArray{T,1};
                    label::String="", count::Integer=min(length(x), length(y2),
                    length(y2)), offset::Integer=0, stride::Integer=1) where {T}

    if eltype(x) == Float32
        LibCImPlot.PlotShadedFloatPtrFloatPtrFloatPtr(label, x, y1, y2, count, offset,
                                                      stride * sizeof(Float32))
    elseif eltype(x) == Float64
        LibCImPlot.PlotShadeddoublePtrdoublePtrdoublePtr(label, x, y1, y2, count, offset,
                                                         stride * sizeof(Float64))
    else
        x = convert.(Float32, x)
        y1 = convert.(Float32, y1)
        y2 = convert.(Float32, y2)
        LibCImPlot.PlotShadedFloatPtrFloatPtrFloatPtr(label, x, y1, y2, count, offset,
                                                         stride * sizeof(Float32))
    end
end

function PlotShaded(x::AbstractArray{T1,1}, y1::AbstractArray{T2,1}, y2::AbstractArray{T3,1};
                    kwargs...) where {T1,T2,T3}
    (x, y1, y2) = promote(x, y1, y2)
    PlotShaded(x, y1, y2; kwargs...)
end

function PlotShaded(x::AbstractArray{T,1}, y::AbstractArray{T,1}, y_ref::T;
                    count::Integer=min(length(x), length(y)), offset::Integer=0,
                    stride::Integer=1, label_id::String="") where {T}
    if eltype(x) == Float32
        LibCImPlot.PlotShadedFloatPtrFloatPtrInt(label, x, y, count, y_ref, offset,
                                                      stride * sizeof(Float32))
    elseif eltype(x) == Float64
        LibCImPlot.PlotShadeddoublePtrdoublePtrInt(label, x, y, count, y_ref,
                                                         offset, stride * sizeof(Float64))
    else
        x = convert.(Float32, x)
        y = convert.(Float32, y)
        y_ref = Float32(y_ref)
        LibCImPlot.PlotShadedFloatPtrFloatPtrInt(label, x, y, count, y_ref,
                                                         offset, stride * sizeof(Float32))
    end
end

function PlotShaded(x::AbstractArray{T1,1}, y::AbstractArray{T2,1}, y_ref::T3;
                    kwargs...) where {T1,T2,T3}
    (x, y) = promote(x, y)
    y_ref = eltype(x)(y_ref)
    PlotShaded(x, y, y_ref; kwargs...)
end

function PlotShaded(x::UnitRange{<:Integer}, y1::AbstractArray{T,1}, y2::AbstractArray{T,1}; label::String="") where {T}
    count = length(x)
    offset = x.start >= 1 ? x.start - 1 : throw("Offset out of bounds")

    if eltype(y1) == Float32
        x = collect(Float32, x)
        LibCImPlot.PlotShadedFloatPtrFloatPtrFloatPtr(label, x, y1, y2, count, offset,
                                                      sizeof(Float32))
    elseif eltype(y1) == Float64
        x = collect(Float64, x)
        LibCImPlot.PlotShadeddoublePtrdoublePtrdoublePtr(label, x, y1, y2, count, offset,
                                                         sizeof(Float64))
    else
        y1 = convert.(Float32, y1)
        y2 = convert.(Float32, y2)
        x = collect(Float32, x)
        LibCImPlot.PlotShadedFloatPtrFloatPtrFloatPtr(label, x, y1, y2, count, offset,
                                                         sizeof(Float32))
    end
end

function PlotShaded(x::UnitRange{<:Integer}, y::AbstractArray{T,1}, y_ref::T; label::String="") where {T}
    count = length(x)
    offset = x.start >= 1 ? x.start - 1 : throw("Offset out of bounds")

    if eltype(y) == Float32
        x = collect(Float32, x)
        LibCImPlot.PlotShadedFloatPtrFloatPtrInt(label, x, y, count, y_ref, offset,
                                                      sizeof(Float32))
    elseif eltype(y) == Float64
        x = collect(Float64, x)
        LibCImPlot.PlotShadeddoublePtrdoublePtrInt(label, x, y, count, y_ref,
                                                         offset,sizeof(Float64))
    else
        x = collect(Float32, x)
        y = convert.(Float32, y)
        y_ref = Float32(y_ref)
        LibCImPlot.PlotShadeddoublePtrdoublePtrInt(label, x, y, count, y_ref,
                                                         offset, sizeof(Float32))
    end
end

function PlotShaded(x::StepRange, y1::AbstractArray{T,1}, y2::AbstractArray{T,1}; label::String="") where {T}
    x.stop < 1 && throw("Range out of bounds")
    count = length(x)
    offset = x.start >= 1 ? x.start - 1 : throw("Range out of bounds")
    stride = x.step >= 1 ? x.step : throw("Stride must be greater than zero.")
    if eltype(y1) == Float32
        x = collect(Float32, x)
        LibCImPlot.PlotShadedFloatPtrFloatPtrFloatPtr(label, x, y1, y2, count, offset,
                                                      stride * sizeof(Float32))
    elseif eltype(y1) == Float64
        x = collect(Float64, x)
        LibCImPlot.PlotShadeddoublePtrdoublePtrdoublePtr(label, x, y1, y2, count, offset,
                                                         stride * sizeof(Float64))
    else
        y1 = convert.(Float32, y1)
        y2 = convert.(Float32, y2)
        x = collect(Float32, x)
        LibCImPlot.PlotShadedFloatPtrFloatPtrFloatPtr(label, x, y1, y2, count, offset,
                                                         stride * sizeof(Float32))
    end
end

function PlotShaded(x::StepRange, y::AbstractArray{T,1}, y_ref::T; label::String="") where {T}
    x.stop < 1 && throw("Range out of bounds")
    count = length(x)
    offset = x.start >= 1 ? x.start - 1 : throw("Range out of bounds")
    stride = x.step >= 1 ? x.step : throw("Stride must be greater than zero.")

    if eltype(y) == Float32
        x = collect(Float32, x)
        LibCImPlot.PlotShadedFloatPtrFloatPtrInt(label, x, y, count, y_ref, offset,
                                                      stride * sizeof(Float32))
    elseif eltype(y) == Float64
        x = collect(Float64, x)
        LibCImPlot.PlotShadeddoublePtrdoublePtrInt(label, x, y, count, y_ref,
                                                         offset, stirde * sizeof(Float64))
    else
    x = collect(Float32, x)
        y = convert.(Float32, y)
        y_ref = Float32(y_ref)
        LibCImPlot.PlotShadedFloatPtrFloatPtrInt(label, x, y, count, y_ref,
                                                         offset, stride * sizeof(Float32))
    end
end

function PlotShaded(x::OrdinalRange, y1::AbstractArray{T1,1}, y2::AbstractArray{T2,1};
                    kwargs...) where {T1,T2}
    (y1, y2) = promote(y1, y2)
    PlotShaded(x, y1, y2; kwargs...)
end

function PlotShaded(x::OrdinalRange, y::AbstractArray{T1,1}, y_ref::T2;
                    kwargs...) where {T1,T2}
    y_ref = eltype(y)(y_ref)
    PlotShaded(x, y, y_ref; kwargs...)
end


text:
struct Plsrda
    fm
    lev::AbstractVector
    ni::AbstractVector
end

struct KplsrDa
    fm
    lev::AbstractVector
    ni::AbstractVector
end

"""
    plsrda(X, y, weights = ones(size(X, 1)); nlv)
Discrimination (DA) based on partial least squares regression (PLSR).
* `X` : X-data.
* `y` : y-data (class membership).
* `weights` : Weights of the observations.
* `nlv` : Nb. latent variables (LVs) to compute.

This is the common "PLSDA".
The training variable y (univariate class membership) is transformed
to a dummy table (Ydummy) containing nlev columns, where nlev is the number
of classes present in `y`. Each column of Ydummy is a dummy variable (0/1).
Then, a PLS2 is implemented on the X-data and Ydummy,
returning `nlv` latent variables (LVs).
Finally, a multiple linear regression (MLR) is run between the LVs and each
column of Ydummy, returning predictions of the dummy variables
(= object `posterior` returned by function `predict`).
These predictions can be considered as unbounded (i.e. eventuall outside of [0, 1])
estimates of the class membership probabilities.
For a given observation, the final prediction is the class corresponding
to the dummy variable for which the prediction is the highest.
"""
function plsrda(X, y, weights = ones(size(X, 1)); nlv)
    z = dummy(y)
    fm = plskern(X, z.Y, weights; nlv = nlv)
    Plsrda(fm, z.lev, z.ni)
end

"""
    kplsrda(X, y, weights = ones(size(X, 1)); nlv, kern = "krbf", kwargs...)
Discrimination (DA) based on kernel partial least squares regression (KPLSR).
* `X` : X-data.
* `y` : Univariate class membership.
* `weights` : Weights of the observations.
* `nlv` : Nb. latent variables (LVs) to compute.

This is the common "PLSDA".
The training variable y (univariate class membership) is transformed
to a dummy table (Ydummy) containing nlev columns, where nlev is the number
of classes present in `y`. Each column of Ydummy is a dummy variable (0/1).
Then, a kernel PLS2 is implemented on the X-data and Ydummy,
returning `nlv` latent variables (LVs).
Finally, a multiple linear regression (MLR) is run between the LVs and each
column of Ydummy, returning predictions of the dummy variables
(= object `posterior` returned by function `predict`).
These predictions can be considered as unbounded (i.e. eventuall outside of [0, 1])
estimates of the class membership probabilities.
For a given observation, the final prediction is the class corresponding
to the dummy variable for which the prediction is the highest.
"""
function kplsrda(X, y, weights = ones(size(X, 1)); nlv, kern = "krbf", kwargs...)
    z = dummy(y)
    fm = kplsr(X, z.Y, weights; nlv = nlv, kern = kern, kwargs...)
    KplsrDa(fm, z.lev, z.ni)
end

"""
    predict(object::Union{Plsrda, KplsrDa}, X; nlv = nothing)
Compute Y-predictions from a fitted model.
* `object` : The fitted model.
* `X` : X-data for which predictions are computed.
* `nlv` : Nb. LVs, or collection of nb. LVs, to consider.
    If nothing, it is the maximum nb. LVs.
"""
function predict(object::Union{Plsrda, KplsrDa}, X; nlv = nothing)
    X = ensure_mat(X)
    m = size(X, 1)
    a = size(object.fm.T, 2)
    isnothing(nlv) ? nlv = a : nlv = (max(minimum(nlv), 0):min(maximum(nlv), a))
    le_nlv = length(nlv)
    pred = list(le_nlv, Matrix{Float64})
    posterior = list(le_nlv, Matrix{Float64})
    @inbounds for i = 1:le_nlv
        zp = predict(object.fm, X; nlv = nlv[i]).pred
        #if softmax
        #    @inbounds for j = 1:m
        #        zp[j, :] .= mweights(exp.(zp[j, :]))
        #   end
        #end
        z =  mapslices(argmax, zp; dims = 2)  # if equal, argmax takes the first
        pred[i] = reshape(replacebylev(z, object.lev), m, 1)
        posterior[i] = zp
    end
    if le_nlv == 1
        pred = pred[1]
        posterior = posterior[1]
    end
    (pred = pred, posterior = posterior)
end


text:
using Test
using MianalyzerBackend
const RandomForestClassifier = MianalyzerBackend.RandomForestClassifier

@testset "RandomForestClassifier" begin
  @testset ".feature_importances" begin
    actual = RandomForestClassifier.feature_importances([1 1 3 -3; 1 1.1 12 -12; 1 -1 100 -100], [3; 12; 100], ["l1", "l2", "l3", "l4"])
    @test actual.columns == ["l1", "l2", "l3", "l4"]
    @test size(actual.feature_importances) == (4,)
  end
end


text:
module BayesianNetworks

using SparseArrays

"""
    AbstractNode

Supertype for all types of `Node`
"""
abstract type AbstractNode end

"""
    DiscreteNode <: AbstractNode

# Fields
- `name::String`: the node name
- `parents::Vector{DiscreteNode}`: the parents of the node
- `categories::Vector{String}`: the categories of the node
- `parameters::SparseMatrixCSC`: the parameters of the node

# Constructors
- `DiscreteNode(name::String, categories::Vector{String})`
"""
mutable struct DiscreteNode <: AbstractNode
    name::String
    parents::Vector{DiscreteNode}

    categories::Vector{String}
    parameters::SparseMatrixCSC
end
DiscreteNode(name::String, categories::Vector) =
    DiscreteNode(name, [], categories, sparse([[];], [[];], [[];]))

"""
    name(node::DiscreteNode)

Returns the name of the node
"""
name(node::DiscreteNode) = node.name

"""
    parents(node::DiscreteNode)

Returns the parents of the node
"""
parents(node::DiscreteNode) = node.parents

"""
    categories(node::DiscreteNode)

Returns the categories of the node
"""
categories(node::DiscreteNode) = node.categories

"""
    parameters(node::DiscreteNode)

Returns the CPT of the node
"""
parameters(node::DiscreteNode) = node.parameters

"""
    add_parent!(node::DiscreteNode, parent::DiscreteNode)

Add a parent to the node
"""
add_parent!(node::DiscreteNode, parent::DiscreteNode) =
    append!(node.parents, [parent])

"""
    set_parameters(node::DiscreteNode, parameter::SparseMatrixCSC)

Add a parameter to the node
"""
set_parameters!(node::DiscreteNode, parameters::SparseMatrixCSC) =
    node.parameters = parameters

"""
    BayesianNetwork{T}

# Fields
- `V::Vector{T}`: The array of nodes
- `E::Vector{Tuple{T}}`: The array of edges

# Reference
- Neapolitan, R. Learning Bayesian Networks, Ch. 1, p. 40
"""
mutable struct BayesianNetwork{T}
    V::Vector{T}
    E::Vector{Tuple{T, T}}
end

export
    AbstractNode,
    DiscreteNode,
    name,
    parents,
    categories,
    parameters,
    add_parent!,
    set_parameters!,

    BayesianNetwork

end # module


text:
#
# Constraint Template Definitions
#
# Constraint templates help simplify data wrangling across multiple Power
# Flow formulations by providing an abstraction layer between the network data
# and network constraint definitions.  The constraint template's job is to
# extract the required parameters from a given network data structure and
# pass the data as named arguments to the Power Flow formulations.
#
# Constraint templates should always be defined over "GenericPowerModel"
# and should never refer to model variables
#

### Voltage Constraints ###

constraint_voltage(pm::GenericPowerModel) = constraint_voltage(pm, pm.cnw)
# no data, so no further templating is needed, constraint goes directly to the formulations

constraint_voltage_on_off(pm::GenericPowerModel) = constraint_voltage_on_off(pm, pm.cnw)
# no data, so no further templating is needed, constraint goes directly to the formulations

constraint_voltage_ne(pm::GenericPowerModel) = constraint_voltage_ne(pm, pm.cnw)
# no data, so no further templating is needed, constraint goes directly to the formulations

### Generator Constraints ###

""
function constraint_active_gen_setpoint(pm::GenericPowerModel, n::Int, i::Int)
    gen = ref(pm, n, :gen, i)
    constraint_active_gen_setpoint(pm, n, gen["index"], gen["pg"])
end
constraint_active_gen_setpoint(pm::GenericPowerModel, i::Int) = constraint_active_gen_setpoint(pm, pm.cnw, i::Int)

""
function constraint_reactive_gen_setpoint(pm::GenericPowerModel, n::Int, i::Int)
    gen = ref(pm, n, :gen, i)
    constraint_reactive_gen_setpoint(pm, gen["index"], gen["qg"])
end
constraint_reactive_gen_setpoint(pm::GenericPowerModel, i::Int) = constraint_reactive_gen_setpoint(pm, pm.cnw, i::Int)

### Bus - Setpoint Constraints ###

""
constraint_theta_ref(pm::GenericPowerModel, i::Int) = constraint_theta_ref(pm, pm.cnw, i)
# no data, so no further templating is needed, constraint goes directly to the formulations

""
function constraint_voltage_magnitude_setpoint(pm::GenericPowerModel, n::Int, i::Int)
    bus = ref(pm, n, :bus, i)
    constraint_voltage_magnitude_setpoint(pm, n, bus["index"], bus["vm"])
end
constraint_voltage_magnitude_setpoint(pm::GenericPowerModel, i::Int) = constraint_voltage_magnitude_setpoint(pm, pm.cnw, i::Int)

### Bus - KCL Constraints ###

""
function constraint_kcl_shunt(pm::GenericPowerModel, n::Int, i::Int)
    if !haskey(pm.con[:nw][n], :kcl_p)
        pm.con[:nw][n][:kcl_p] = Dict{Int,ConstraintRef}()
    end
    if !haskey(pm.con[:nw][n], :kcl_q)
        pm.con[:nw][n][:kcl_q] = Dict{Int,ConstraintRef}()
    end

    bus = ref(pm, n, :bus, i)
    bus_arcs = ref(pm, n, :bus_arcs, i)
    bus_arcs_dc = ref(pm, n, :bus_arcs_dc, i)
    bus_gens = ref(pm, n, :bus_gens, i)

    constraint_kcl_shunt(pm, n, i, bus_arcs, bus_arcs_dc, bus_gens, bus["pd"], bus["qd"], bus["gs"], bus["bs"])
end
constraint_kcl_shunt(pm::GenericPowerModel, i::Int) = constraint_kcl_shunt(pm, pm.cnw, i::Int)

""
function constraint_kcl_shunt_ne(pm::GenericPowerModel, n::Int, i::Int)
    bus = ref(pm, n, :bus, i)
    bus_arcs = ref(pm, n, :bus_arcs, i)
    bus_arcs_dc = ref(pm, n, :bus_arcs_dc, i)
    bus_arcs_ne = ref(pm, n, :ne_bus_arcs, i)
    bus_gens = ref(pm, n, :bus_gens, i)

    constraint_kcl_shunt_ne(pm, n, i, bus_arcs, bus_arcs_dc, bus_arcs_ne, bus_gens, bus["pd"], bus["qd"], bus["gs"], bus["bs"])
end
constraint_kcl_shunt_ne(pm::GenericPowerModel, i::Int) = constraint_kcl_shunt_ne(pm, pm.cnw, i::Int)

### Branch - Ohm's Law Constraints ###

""
function constraint_ohms_yt_from(pm::GenericPowerModel, n::Int, i::Int)
    branch = ref(pm, n, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    g, b = calc_branch_y(branch)
    tr, ti = calc_branch_t(branch)
    c = branch["br_b"]
    tm = branch["tap"]

    constraint_ohms_yt_from(pm, n, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm)
end
constraint_ohms_yt_from(pm::GenericPowerModel, i::Int) = constraint_ohms_yt_from(pm, pm.cnw, i)

""
function constraint_ohms_yt_to(pm::GenericPowerModel, n::Int, i::Int)
    branch = ref(pm, n, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    g, b = calc_branch_y(branch)
    tr, ti = calc_branch_t(branch)
    c = branch["br_b"]
    tm = branch["tap"]

    constraint_ohms_yt_to(pm, n, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm)
end
constraint_ohms_yt_to(pm::GenericPowerModel, i::Int) = constraint_ohms_yt_to(pm, pm.cnw, i)

""
function constraint_ohms_y_from(pm::GenericPowerModel, n::Int, i::Int)
    branch = ref(pm, n, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    g, b = calc_branch_y(branch)
    c = branch["br_b"]
    tm = branch["tap"]
    ta = branch["shift"]

    constraint_ohms_y_from(pm, n, f_bus, t_bus, f_idx, t_idx, g, b, c, tm, ta)
end
constraint_ohms_y_from(pm::GenericPowerModel, i::Int) = constraint_ohms_y_from(pm, pm.cnw, i)

""
function constraint_ohms_y_to(pm::GenericPowerModel, n::Int, i::Int)
    branch = ref(pm, n, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    g, b = calc_branch_y(branch)
    c = branch["br_b"]
    tm = branch["tap"]
    ta = branch["shift"]

    constraint_ohms_y_to(pm, n, f_bus, t_bus, f_idx, t_idx, g, b, c, tm, ta)
end
constraint_ohms_y_to(pm::GenericPowerModel, i::Int) = constraint_ohms_y_to(pm, pm.cnw, i)

### DC LINES ###

""
function constraint_dcline(pm::GenericPowerModel, n::Int, i::Int)
    dcline = ref(pm, n, :dcline, i)
    f_bus = dcline["f_bus"]
    t_bus = dcline["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)
    loss0 = dcline["loss0"]
    loss1 = dcline["loss1"]

    constraint_dcline(pm, n, f_bus, t_bus, f_idx, t_idx, loss0, loss1)
end
constraint_dcline(pm::GenericPowerModel, i::Int) = constraint_dcline(pm, pm.cnw, i)

function constraint_active_dcline_setpoint(pm::GenericPowerModel, n::Int, i::Int)
    dcline = ref(pm, n, :dcline, i)
    f_bus = dcline["f_bus"]
    t_bus = dcline["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)
    pf = dcline["pf"]
    pt = dcline["pt"]

    constraint_active_dcline_setpoint(pm, n, f_idx, t_idx, pf, pt)
end
constraint_active_dcline_setpoint(pm::GenericPowerModel, i::Int) = constraint_active_dcline_setpoint(pm, pm.cnw, i)

### Branch - On/Off Ohm's Law Constraints ###

""
function constraint_ohms_yt_from_on_off(pm::GenericPowerModel, n::Int, i::Int)
    branch = ref(pm, n, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    g, b = calc_branch_y(branch)
    tr, ti = calc_branch_t(branch)
    c = branch["br_b"]
    tm = branch["tap"]

    vad_min = ref(pm, n, :off_angmin)
    vad_max = ref(pm, n, :off_angmax)

    constraint_ohms_yt_from_on_off(pm, n, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, vad_min, vad_max)
end
constraint_ohms_yt_from_on_off(pm::GenericPowerModel, i::Int) = constraint_ohms_yt_from_on_off(pm, pm.cnw, i)

""
function constraint_ohms_yt_to_on_off(pm::GenericPowerModel, n::Int, i::Int)
    branch = ref(pm, n, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    g, b = calc_branch_y(branch)
    tr, ti = calc_branch_t(branch)
    c = branch["br_b"]
    tm = branch["tap"]

    vad_min = ref(pm, n, :off_angmin)
    vad_max = ref(pm, n, :off_angmax)

    constraint_ohms_yt_to_on_off(pm, n, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, vad_min, vad_max)
end
constraint_ohms_yt_to_on_off(pm::GenericPowerModel, i::Int) = constraint_ohms_yt_to_on_off(pm, pm.cnw, i)

""
function constraint_ohms_yt_from_ne(pm::GenericPowerModel, n::Int, i::Int)
    branch = ref(pm, n, :ne_branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    g, b = calc_branch_y(branch)
    tr, ti = calc_branch_t(branch)
    c = branch["br_b"]
    tm = branch["tap"]

    vad_min = ref(pm, n, :off_angmin)
    vad_max = ref(pm, n, :off_angmax)

    constraint_ohms_yt_from_ne(pm, n, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, vad_min, vad_max)
end
constraint_ohms_yt_from_ne(pm::GenericPowerModel, i::Int) = constraint_ohms_yt_from_ne(pm, pm.cnw, i)

""
function constraint_ohms_yt_to_ne(pm::GenericPowerModel, n::Int, i::Int)
    branch = ref(pm, n, :ne_branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    g, b = calc_branch_y(branch)
    tr, ti = calc_branch_t(branch)
    c = branch["br_b"]
    tm = branch["tap"]

    vad_min = ref(pm, n, :off_angmin)
    vad_max = ref(pm, n, :off_angmax)

    constraint_ohms_yt_to_ne(pm, n, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, vad_min, vad_max)
end
constraint_ohms_yt_to_ne(pm::GenericPowerModel, i::Int) = constraint_ohms_yt_to_ne(pm, pm.cnw, i)

### Branch - Current ###

""
function constraint_power_magnitude_sqr(pm::GenericPowerModel, n::Int, i::Int)
    branch = ref(pm, n, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    arc_from = (i, f_bus, t_bus)

    tm = branch["tap"]

    constraint_power_magnitude_sqr(pm, n, f_bus, t_bus, arc_from, tm)
end
constraint_power_magnitude_sqr(pm::GenericPowerModel, i::Int) = constraint_power_magnitude_sqr(pm, pm.cnw, i)

""
function constraint_power_magnitude_sqr_on_off(pm::GenericPowerModel, n::Int, i::Int)
    branch = ref(pm, n, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    arc_from = (i, f_bus, t_bus)

    tm = branch["tap"]

    constraint_power_magnitude_sqr_on_off(pm, n, i, f_bus, arc_from, tm)
end
constraint_power_magnitude_sqr_on_off(pm::GenericPowerModel, i::Int) = constraint_power_magnitude_sqr_on_off(pm, pm.cnw, i)

""
function constraint_power_magnitude_link(pm::GenericPowerModel, n::Int, i::Int)
    branch = ref(pm, n, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    arc_from = (i, f_bus, t_bus)

    g, b = calc_branch_y(branch)
    tr, ti = calc_branch_t(branch)
    c = branch["br_b"]
    tm = branch["tap"]

    constraint_power_magnitude_link(pm, n, f_bus, t_bus, arc_from, g, b, c, tr, ti, tm)
end
constraint_power_magnitude_link(pm::GenericPowerModel, i::Int) = constraint_power_magnitude_link(pm, pm.cnw, i)

""
function constraint_power_magnitude_link_on_off(pm::GenericPowerModel, n::Int, i::Int)
    branch = ref(pm, n, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    arc_from = (i, f_bus, t_bus)

    g, b = calc_branch_y(branch)
    tr, ti = calc_branch_t(branch)
    c = branch["br_b"]
    tm = branch["tap"]

    constraint_power_magnitude_link_on_off(pm, n, i, arc_from, g, b, c, tr, ti, tm)
end
constraint_power_magnitude_link_on_off(pm::GenericPowerModel, i::Int) = constraint_power_magnitude_link_on_off(pm, pm.cnw, i)

### Branch - Thermal Limit Constraints ###

"""

    constraint_thermal_limit_from(pm::GenericPowerModel, n::Int, i::Int)

Adds the (upper and lower) thermal limit constraints for the desired branch to the PowerModel.

"""
function constraint_thermal_limit_from(pm::GenericPowerModel, n::Int, i::Int)
    if !haskey(pm.con[:nw][n], :sm_fr)
        pm.con[:nw][n][:sm_fr] = Dict{Int,Any}() # note this can be a constraint or variable
    end

    branch = ref(pm, n, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)

    constraint_thermal_limit_from(pm, n, f_idx, branch["rate_a"])
end
constraint_thermal_limit_from(pm::GenericPowerModel, i::Int) = constraint_thermal_limit_from(pm, pm.cnw, i)

""
function constraint_thermal_limit_to(pm::GenericPowerModel, n::Int, i::Int)
    if !haskey(pm.con[:nw][n], :sm_to)
        pm.con[:nw][n][:sm_to] = Dict{Int,Any}() # note this can be a constraint or variable
    end

    branch = ref(pm, n, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    t_idx = (i, t_bus, f_bus)

    constraint_thermal_limit_to(pm, n, t_idx, branch["rate_a"])
end
constraint_thermal_limit_to(pm::GenericPowerModel, i::Int) = constraint_thermal_limit_to(pm, pm.cnw, i)

""
function constraint_thermal_limit_from_on_off(pm::GenericPowerModel, n::Int, i::Int)
    branch = ref(pm, n, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)

    constraint_thermal_limit_from_on_off(pm, n, i, f_idx, branch["rate_a"])
end
constraint_thermal_limit_from_on_off(pm::GenericPowerModel, i::Int) = constraint_thermal_limit_from_on_off(pm, pm.cnw, i)

""
function constraint_thermal_limit_to_on_off(pm::GenericPowerModel, n::Int, i::Int)
    branch = ref(pm, n, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    t_idx = (i, t_bus, f_bus)

    constraint_thermal_limit_to_on_off(pm, n, i, t_idx, branch["rate_a"])
end
constraint_thermal_limit_to_on_off(pm::GenericPowerModel, i::Int) = constraint_thermal_limit_to_on_off(pm, pm.cnw, i)

""
function constraint_thermal_limit_from_ne(pm::GenericPowerModel, n::Int, i::Int)
    branch = ref(pm, n, :ne_branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)

    constraint_thermal_limit_from_ne(pm, n, i, f_idx, branch["rate_a"])
end
constraint_thermal_limit_from_ne(pm::GenericPowerModel, i::Int) = constraint_thermal_limit_from_ne(pm, pm.cnw, i)

""
function constraint_thermal_limit_to_ne(pm::GenericPowerModel, n::Int, i::Int)
    branch = ref(pm, n, :ne_branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    t_idx = (i, t_bus, f_bus)

    constraint_thermal_limit_to_ne(pm, n, i, t_idx, branch["rate_a"])
end
constraint_thermal_limit_to_ne(pm::GenericPowerModel, i::Int) = constraint_thermal_limit_to_ne(pm, pm.cnw, i)

### Branch - Phase Angle Difference Constraints ###

""
function constraint_voltage_angle_difference(pm::GenericPowerModel, n::Int, i::Int)
    branch = ref(pm, n, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    pair = (f_bus, t_bus)
    buspair = ref(pm, n, :buspairs, pair)

    if buspair["branch"] == i
        constraint_voltage_angle_difference(pm, n, f_bus, t_bus, buspair["angmin"], buspair["angmax"])
    end
end
constraint_voltage_angle_difference(pm::GenericPowerModel, i::Int) = constraint_voltage_angle_difference(pm, pm.cnw, i)

""
function constraint_voltage_angle_difference_on_off(pm::GenericPowerModel, n::Int, i::Int)
    branch = ref(pm, n, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]

    vad_min = ref(pm, n, :off_angmin)
    vad_max = ref(pm, n, :off_angmax)

    constraint_voltage_angle_difference_on_off(pm, n, i, f_bus, t_bus, branch["angmin"], branch["angmax"], vad_min, vad_max)
end
constraint_voltage_angle_difference_on_off(pm::GenericPowerModel, i::Int) = constraint_voltage_angle_difference_on_off(pm, pm.cnw, i)

""
function constraint_voltage_angle_difference_ne(pm::GenericPowerModel, n::Int, i::Int)
    branch = ref(pm, n, :ne_branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]

    vad_min = ref(pm, n, :off_angmin)
    vad_max = ref(pm, n, :off_angmax)

    constraint_voltage_angle_difference_ne(pm, n, i, f_bus, t_bus, branch["angmin"], branch["angmax"], vad_min, vad_max)
end
constraint_voltage_angle_difference_ne(pm::GenericPowerModel, i::Int) = constraint_voltage_angle_difference_ne(pm, pm.cnw, i)

### Branch - Loss Constraints ###

""
function constraint_loss_lb(pm::GenericPowerModel, n::Int, i::Int)
    branch = ref(pm, n, :branch, i)
    @assert branch["br_r"] >= 0
    @assert branch["br_x"] >= 0
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    c = branch["br_b"]
    tr = branch["tr"]

    constraint_loss_lb(pm, n, f_bus, t_bus, f_idx, t_idx, c, tr)
end
constraint_loss_lb(pm::GenericPowerModel, i::Int) = constraint_loss_lb(pm, pm.cnw, i)


text:
using Revise
using POMDPs
using GridGoalMDPs
using DiscreteValueIteration

m = TwoAgentMDP(size=(25,25))
@show n_states(m)
@show n_actions(m)

solver = SparseValueIterationSolver(verbose=true)

p = solve(solver, m)

@show sizeof(p.qmat)


text:
using BinaryBuilder

# These are the platforms built inside the wizard
platforms = [
    BinaryProvider.Linux(:x86_64, :glibc)
]

# If the user passed in a platform (or a few, comma-separated) on the
# command-line, use that instead of our default platforms
if length(ARGS) > 0
    platforms = platform_key.(split(ARGS[1], ","))
end
info("Building for $(join(triplet.(platforms), ", "))")

# Collection of sources required to build linux
sources = [
    "https://github.com/torvalds/linux/archive/v4.15-rc2.tar.gz" =>
    "98d1d188424d280a3e07b3dbff589422e44124145f7ed2003f38c81be57ce9c4",
]

script = raw"""
cd $WORKSPACE/srcdir
export PATH=/usr/bin:$PATH
cd linux-4.15-rc2/
curl -OL https://raw.githubusercontent.com/JuliaPackaging/BinaryBuilder.jl/e6ac598dc48023f57f470d2a2fb6d10d450dbced/deps/linuxkernel.config
cp linuxkernel.config arch/x86/configs/binarybuilder_defconfig
apk add libelf-dev openssl-dev libelf-dev gcc linux-headers musl-dev bc
make binarybuilder_defconfig
make -j40
cp vmlinux $DESTDIR
cp arch/x86/boot/bzImage $DESTDIR

"""

products = prefix -> [
    ExecutableProduct(prefix,"vmlinux")
]

# Build the given platforms using the given sources
autobuild(pwd(), "linux", platforms, sources, script, products)


text:
using Documenter, Polynomials

DocMeta.setdocmeta!(Polynomials, :DocTestSetup, :(using Polynomials); recursive=true)

makedocs(modules = [Polynomials],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "Polynomials.jl",
    authors = "Jameson Nash, Keno Fischer, and other contributors",
    pages = [
        "Manual" => "index.md",
    ],
)

deploydocs(repo = "github.com/JuliaMath/Polynomials.jl.git")


text:
## Problems

function build_max_import(pm::_PM.AbstractBFModel)
    _FP.post_flex_tnep(pm; objective = false, intertemporal_constraints = false)
    constraint_same_investments(pm)
    objective_max_import(pm)
end

function build_max_export(pm::_PM.AbstractBFModel)
    _FP.post_flex_tnep(pm; objective = false, intertemporal_constraints = false)
    for n in _PM.nw_ids(pm)
        constraint_load_pcurt_zero(pm, n)
    end
    constraint_same_investments(pm)
    objective_max_export(pm)
end

function build_min_cost_with_same_investments_in_all_sub_nws(pm::_PM.AbstractBFModel)
    _FP.post_flex_tnep(pm; objective = false, intertemporal_constraints = false)
    for n in _PM.nw_ids(pm)
        constraint_load_pcurt_zero(pm, n)
    end
    constraint_same_investments(pm)
    objective_min_investment_cost(pm)
end

function build_max_import_with_current_investments(pm::_PM.AbstractBFModel)
    _FP.post_flex_tnep(pm; objective = false, intertemporal_constraints = false)
    for n in _PM.nw_ids(pm)
        constraint_ne_branch_indicator_fix(pm, n)
        constraint_ne_storage_indicator_fix(pm, n)
    end
    objective_max_import(pm)
end

function build_max_export_with_current_investments(pm::_PM.AbstractBFModel)
    _FP.post_flex_tnep(pm; objective = false, intertemporal_constraints = false)
    for n in _PM.nw_ids(pm)
        constraint_ne_branch_indicator_fix(pm, n)
        constraint_ne_storage_indicator_fix(pm, n)
        constraint_load_pcurt_zero(pm, n)
    end
    objective_max_export(pm)
end

function build_min_cost_with_fixed_investments(pm::_PM.AbstractBFModel)
    _FP.post_flex_tnep(pm; objective = true, intertemporal_constraints = true)
    for n in _PM.nw_ids(pm)
        constraint_ne_branch_indicator_fix(pm, n)
        constraint_ne_storage_indicator_fix(pm, n)
    end
end

function build_max_flex_band_with_bounded_cost(pm::_PM.AbstractBFModel)
    _FP.post_flex_tnep(pm; objective = false, intertemporal_constraints = false)
    for n in _PM.nw_ids(pm)
        constraint_load_pcurt_zero(pm, n)
    end
    constraint_same_investments(pm)
    constraint_investment_cost_max(pm; sub_nw = 1)
    objective_max_flex_band_2sub(pm)
end

## Constraints

"Disable involuntary curtailment of loads"
function constraint_load_pcurt_zero(pm::_PM.AbstractPowerModel, n::Int)
    for i in _PM.ids(pm, n, :load)
        pcurt = _PM.var(pm, n, :pcurt, i)
        JuMP.@constraint(pm.model, pcurt == 0.0)
    end
end

"Fix investment decisions on candidate branches according to values in data structure"
function constraint_ne_branch_indicator_fix(pm::_PM.AbstractPowerModel, n::Int)
    for i in _PM.ids(pm, n, :ne_branch)
        indicator = _PM.var(pm, n, :branch_ne, i)
        value = _PM.ref(pm, n, :ne_branch, i, "built")
        JuMP.@constraint(pm.model, indicator == value)
    end
end

"Fix investment decisions on candidate storage according to values in data structure"
function constraint_ne_storage_indicator_fix(pm::_PM.AbstractPowerModel, n::Int)
    for i in _PM.ids(pm, n, :ne_storage)
        indicator = _PM.var(pm, n, :z_strg_ne, i)
        value = _PM.ref(pm, n, :ne_storage, i, "isbuilt")
        JuMP.@constraint(pm.model, indicator == value)
    end
end

"Put an upper bound on investment cost related to a single subnetwork"
function constraint_investment_cost_max(pm::_PM.AbstractPowerModel; sub_nw = 1)
    inv_cost = sum(
            _FP.calc_ne_branch_cost(pm, n)
            + _FP.calc_ne_storage_cost(pm, n)
            + _FP.calc_load_investment_cost(pm, n)
        for n in _FP.nw_ids(pm; hour=1, scenario=1, sub_nw))
    max_cost = pm.ref[:it][_PM.pm_it_sym][:max_cost]

    JuMP.@constraint(pm.model, inv_cost <= max_cost)
end

"Ensure that investment decisions are the same in each subnetwork"
function constraint_same_investments(pm::_PM.AbstractPowerModel)
    for y in 1:_FP.dim_length(pm, :year)
        n_1, rest = Iterators.peel(_FP.nw_ids(pm; hour=1, scenario=1, year=y))
        for n_2 in rest
            for i in _PM.ids(pm, :ne_branch, nw = n_2)
                _FP.constraint_ne_branch_investment_same(pm, n_1, n_2, i)
            end
            for i in _PM.ids(pm, :ne_storage, nw = n_2)
                _FP.constraint_ne_storage_investment_same(pm, n_1, n_2, i)
            end
            n_1 = n_2
        end
    end
end

## Objectives

function objective_max_import(pm::_PM.AbstractPowerModel)
    return JuMP.@objective(pm.model, Max,
        sum( calc_td_coupling_power_active(pm, n) for (n, nw_ref) in _PM.nws(pm) )
    )
end

function objective_max_export(pm::_PM.AbstractPowerModel)
    return JuMP.@objective(pm.model, Min,
        sum( calc_td_coupling_power_active(pm, n) for (n, nw_ref) in _PM.nws(pm) )
    )
end

function calc_td_coupling_power_active(pm::_PM.AbstractPowerModel, n::Int)
    pcc_gen = _FP.dim_prop(pm, n, :sub_nw, "d_gen")
    p = _PM.var(pm, n, :pg, pcc_gen)
    return p
end

function objective_max_flex_band_2sub(pm::_PM.AbstractPowerModel)
    return JuMP.@objective(pm.model, Max,
        sum( calc_td_coupling_power_active(pm, n) for n in _FP.nw_ids(pm; hour=1, scenario=1, sub_nw=1) )
        - sum( calc_td_coupling_power_active(pm, n) for n in _FP.nw_ids(pm; hour=1, scenario=1, sub_nw=2) )
    )
end

function objective_min_investment_cost(pm::_PM.AbstractPowerModel)
    investment = sum(
        _FP.calc_ne_branch_cost(pm, n)
        + _FP.calc_ne_storage_cost(pm, n)
        + _FP.calc_load_investment_cost(pm, n)
        for n in _FP.nw_ids(pm; hour=1, scenario=1)
    )
    JuMP.@objective(pm.model, Min, investment)
end


text:
module Layouts

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element

export layout, pageContainer, page

function __init__()
  Genie.Renderer.Html.register_normal_element("q__layout", context = Genie.Renderer.Html)
  Genie.Renderer.Html.register_normal_element("q__page__container", context = Genie.Renderer.Html)
  Genie.Renderer.Html.register_normal_element("q__page", context = Genie.Renderer.Html)
end

function layout(args...;
                wrap::Function = StippleUI.DEFAULT_WRAPPER,
                kwargs...)
     wrap() do
      Genie.Renderer.Html.q__layout(args...; attributes([kwargs...], StippleUI.API.ATTRIBUTES_MAPPINGS)...)
     end
end

function page_container(args...;
                wrap::Function = StippleUI.DEFAULT_WRAPPER,
                kwargs...)
      wrap() do
        Genie.Renderer.Html.q__page__container(args...; attributes([kwargs...], StippleUI.API.ATTRIBUTES_MAPPINGS)...)
      end
end

function page(args...;
              wrap::Function = StippleUI.DEFAULT_WRAPPER,
              kwargs...)
      wrap() do
        Genie.Renderer.Html.q__page(args...; attributes([kwargs...], StippleUI.API.ATTRIBUTES_MAPPINGS)...)
      end
end

end


text:
function simulate_single_echo_mag(mag, TEs, TE_SE=mean(TEs))
    weighting = get_single_echo_weighting(TEs, TE_SE)
    weighting = to_dim(weighting, 4)
    return dropdims(sum(mag .* weighting; dims=4); dims=4)
end

function get_single_echo_weighting(TEs, TE_SE)
    if TE_SE < TEs[1] || TE_SE > TEs[end]
        error("Not possible to simulate TE=$TE_SE from $TEs !")
    end
    ΔTE = TEs[2] - TEs[1]
    echoend_ME = TEs[end] + ΔTE / 2
    echostart_ME = TEs[1] - ΔTE / 2
    echowidth_sim = 2min(echoend_ME - TE_SE, TE_SE - echostart_ME)

    echostart_sim = TE_SE - echowidth_sim / 2
    echoend_sim = TE_SE + echowidth_sim / 2

    return get_single_echo_weighting(TEs, echostart_sim, echoend_sim)
end

function get_single_echo_weighting(TEs, echostart_sim, echoend_sim)
    ΔTE = TEs[2] - TEs[1]
    echoend_ME = TEs[end] + ΔTE / 2
    echostart_ME = TEs[1] - ΔTE / 2
    if (echostart_sim + 1e-5) < echostart_ME || (echoend_sim - 1e-5) > echoend_ME || (echoend_sim - echostart_sim + 1e-5) < ΔTE
        error("Not possible to simulate [$(echostart_sim);$(echoend_sim)] from $TEs !")
    end
    weighting = ones(length(TEs))
    if echostart_sim > echostart_ME
        lowechoborder = argmin(abs.(TEs .- echostart_sim))
        weighting[1:lowechoborder] .= 0
        weighting[lowechoborder] = (TEs[lowechoborder] + ΔTE/2 - echostart_sim) / ΔTE
    elseif echoend_sim < echoend_ME
        highechoborder = argmin(abs.(TEs .- echoend_sim))
        weighting[highechoborder:end] .= 0
        weighting[highechoborder] = (echoend_sim - (TEs[highechoborder] - ΔTE/2)) / ΔTE
    end
    return weighting
end


text:
#  Outside a domain recognized by hirlam.org create a ~/.netrc file (chmod 600)
#
#  machine hirlam.org
#  login <username>
#  password <password>

function fetch_wikipages()
    # wikiurl = "https://hirlam.org/trac/wiki/"

    files = "$(@__DIR__)/hsdfiles"

    @sync for file in readlines(files)
       #  @async run(`curl -n -sS --create-dirs $wikiurl$file\?format=txt -o $WIKIDIR/$file.wiki`)
       println("Fetching $wikiurl$file?format=txt")
       run(`curl -n -sS --create-dirs $wikiurl$file\?format=txt -o $WIKIDIR/$file.wiki`)
    end
end


text:
__precompile__()

module Sundials

using Compat

const depsfile = joinpath(dirname(dirname(@__FILE__)),"deps","deps.jl")
if isfile(depsfile)
    include(depsfile)
else
    error("Sundials is not properly installed. Please run Pkg.build(\"Sundials\")")
end

##################################################################
# Deprecations
##################################################################

@deprecate nvlength length
@deprecate asarray convert
@deprecate nvector NVector

# some definitions from the system C headers wrapped into the types_and_consts.jl
const DBL_MAX = prevfloat(Inf)
const DBL_MIN = nextfloat(-Inf)
const DBL_EPSILON = eps(Cdouble)

typealias FILE Void
typealias __builtin_va_list Ptr{Void}

include("types_and_consts.jl")

include("handle.jl")
include("nvector_wrapper.jl")

include("nvector.jl")
include("libsundials.jl")
if isdefined(:libsundials_cvodes)
    include("cvodes.jl")
else
    include("cvode.jl")
end
if isdefined(:libsundials_idas)
    include("idas.jl")
else
    include("ida.jl")
end
include("kinsol.jl")

include("simple.jl")

end # module


text:
# Test whether PolyhedralCone projects points onto same locations as equivalent FiniteCones

using Distributions
using Random

Random.seed!(1)

print("\tTest 1 - dim == 2...")
Z = Normal()
μ = rand(Z, 2)
A = rand(Z, 2, 2)
Σ = A'A
dist = MvNormal(μ, Σ)

K_f = FiniteCone([1 0; 0 1])
K_p = PolyhedralCone([1 0; 0 1])

for i in 1:1000
    x = rand(dist)
    @test project(K_f, x) ≈ project(K_p, x) atol=1e-4
end

println("success!")


text:
export StochasticEnv

using StatsBase: sample, Weights

struct StochasticEnv{E<:AbstractEnv,R} <: AbstractEnvWrapper
    env::E
    rng::R
end

function StochasticEnv(env; rng = Random.GLOBAL_RNG)
    ChanceStyle(env) === EXPLICIT_STOCHASTIC ||
        throw(ArgumentError("only environments of EXPLICIT_STOCHASTIC style is supported"))
    env = StochasticEnv(env, rng)
    reset!(env)
    env
end

function RLBase.reset!(env::StochasticEnv)
    reset!(env.env)
    while current_player(env.env) == chance_player(env.env)
        p = prob(env.env)
        A = action_space(env.env)
        x = A[sample(env.rng, Weights(p, 1.0))]
        env.env(x)
    end
end

function (env::StochasticEnv)(a)
    env.env(a)
    while current_player(env.env) == chance_player(env.env)
        p = prob(env.env)
        A = action_space(env.env)
        x = A[sample(env.rng, Weights(p, 1.0))]
        env.env(x)
    end
end

RLBase.ChanceStyle(::StochasticEnv) = STOCHASTIC
RLBase.players(env::StochasticEnv) =
    [p for p in players(env.env) if p != chance_player(env.env)]
Random.seed!(env::StochasticEnv, s) = Random.seed!(env.rng, s)


text:
using NLPModels
using SolverTools

abstract type TR1DModel{T} end

abstract type QuadraticModel{T} <: TR1DModel{T} end

"Quadratic model solution"
function solveModel(M:: QuadraticModel{T}, t :: T) where T
    return -M.g/M.H
end

"Quadratic model solution in region  Δn ≤ t ≤ Δp "
function solveModel(M:: QuadraticModel{T}, t :: T, Δp :: T, Δn :: T) where T
    dN = -M.g/M.H

    if (Δm(M, Δp, Δn) < 0.0) #| (Δn==0.0)
        d=Δp
        @debug "dp"
    else
        d=Δn
        @debug "dn"
    end

    if  (dN < Δp) & (dN > Δn) & (Δm(M, d, dN) > 0.0)
        d=dN
        @debug "dN = ", d
    end

    return d
end

"Quadratic model evaluation"
function m(M:: QuadraticModel{T}, t :: T) where T
    return M.f + M.g * t + T(0.5)*M.H * t^2
end

"Quadratic model predicted reduction"
function Δm(M:: QuadraticModel{T}, t :: T) where T
    return M.g * t  +  T(0.5) * M.H * t^2
end

"Quadratic model predicted reduction given two points"
function Δm(M:: QuadraticModel{T}, t1 :: T, t2 :: T) where T
    #q1 = M.g * t1  +  T(0.5) * M.H * t1^2
    #q2 = M.g * t2  +  T(0.5) * M.H * t2^2
    #diff = q1 - q2
    diff = (M.g + T(0.5) * M.H * (t1+t2)) * (t1 - t2)
    return diff
end

"Taylor quadratic polynomial model"
mutable struct Taylor2{T} <: QuadraticModel{T}
    f :: T
    g :: T
    H :: T
end
Taylor2{T}() where T = Taylor2{T}(0,0,0)

function initModel(M:: Taylor2{T}, ϕ :: OneDModel, t :: T, ϕt :: T, dϕt :: T) where T
    M.f = ϕt
    M.g = dϕt
    M.H = hess(ϕ, t)
    return M
end

function updateModel!(M:: Taylor2{T}, ϕ :: OneDModel, t :: T, ϕt :: T, dϕt :: T) where T
    M.f = ϕt
    M.g  = dϕt
    M.H = hess(ϕ, t)
    return M
end

function initH_N(ϕ :: LineModel, t :: Float64, bidon)
    return hess(ϕ, t)
end

function updateH_N(ϕ :: LineModel, t :: Float64, dϕ :: Float64)
    return hess(ϕ, t)
end

"Secant quadratic model"
mutable struct secant{T} <: QuadraticModel{T}
    f :: T
    g :: T
    H :: T
end
secant{T}() where T = secant{T}(0,0,0)

function initModel(M:: secant{T}, ϕ :: OneDModel, t :: T, ϕt :: T, dϕt :: T) where T
    M.f = ϕt
    M.g = dϕt
    M.H = initH_s!(ϕ, t, dϕt)
    return M
end

function updateModel!(M:: secant{T}, ϕ :: OneDModel, t :: T, ϕt :: T, dϕt :: T) where T
    M.f = ϕt
    M.g = dϕt
    M.H = updateH_s!(ϕ, t, dϕt)
    return M
end

let (tpred, dϕpred) = (NaN, NaN)

    global function initH_s!(ϕ :: OneDModel, t :: Float64, dϕ :: Float64)
        (tpred, dϕpred) = (t, dϕ)
        return 1.0
    end

    global function updateH_s!(ϕ :: OneDModel, t :: Float64, dϕt :: Float64)
        s = t - tpred
        y = dϕt - dϕpred
        H = y/s

        tpred, dϕpred = t, dϕt

        return H
    end
end


text:
include("ceil_floor.jl")

"""
    median(nums)

Finds median of a vector of numbers

# Example

```julia
median([2,1,3,4])                   # returns 2.5
median([2, 70, 6, 50, 20, 8, 4])    # returns 8
median([0])                         # returns 0
```

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee)
"""
function median(nums)
    sorted = sort(nums)
    len = length(sorted)
    mid_index = Int(floor_val(len / 2))
    return len % 2 == 0 ? (sorted[mid_index+1] + sorted[mid_index]) / 2 : sorted[mid_index+1]
end

max_stars_repo_name: arubhardwaj/Julia

max_issues_repo_name: arubhardwaj/Julia

max_forks_repo_name: arubhardwaj/Julia


text:
# TODO make sure that setting nlayers=1 works

module MultilayerQG

export
  fwdtransform!,
  invtransform!,
  streamfunctionfrompv!,
  pvfromstreamfunction!,
  updatevars!,

  set_q!,
  set_psi!,
  energies,
  fluxes

using
  FFTW,
  LinearAlgebra,
  Reexport

@reexport using FourierFlows

using LinearAlgebra: mul!, ldiv!
using FFTW: rfft, irfft
using FourierFlows: getfieldspecs, varsexpression, parsevalsum, parsevalsum2, superzeros

nothingfunction(args...) = nothing

"""
    Problem(; parameters...)

Construct a multi-layer QG problem.
"""
function Problem(;
    # Numerical parameters
          nx = 128,
          Lx = 2π,
          ny = nx,
          Ly = Lx,
          dt = 0.01,
    # Physical parameters
     nlayers = 2,                       # number of fluid layers
          f0 = 1.0,                     # Coriolis parameter
           β = 0.0,                     # y-gradient of Coriolis parameter
           g = 1.0,                     # gravitational constant
           U = zeros(ny, nlayers),      # imposed zonal flow U(y) in each layer
           H = [0.2, 0.8],              # rest fluid height of each layer
           ρ = [4.0, 5.0],              # density of each layer
         eta = zeros(nx, ny),           # topographic PV
    # Bottom Drag and/or (hyper)-viscosity
           μ = 0.0,
           ν = 0.0,
          nν = 1,
    # Timestepper and eqn options
     stepper = "RK4",
      calcFq = nothingfunction,
      linear = false,
           T = Float64)

   grid = TwoDGrid(nx, Lx, ny, Ly; T=T)
   params = Params(nlayers, g, f0, β, Array{T}(ρ), Array{T}(H), Array{T}(U), Array{T}(eta), μ, ν, nν, grid, calcFq=calcFq)
   vars = calcFq == nothingfunction ? Vars(grid, params) : ForcedVars(grid, params)
   eqn = linear ? LinearEquation(params, grid) : Equation(params, grid)

  FourierFlows.Problem(eqn, stepper, dt, grid, vars, params)
end

abstract type BarotropicParams <: AbstractParams end

struct Params{T} <: AbstractParams
  # prescribed params
   nlayers :: Int            # Number of fluid layers
         g :: T              # Gravitational constant
        f0 :: T              # Constant planetary vorticity
         β :: T              # Planetary vorticity y-gradient
         ρ :: Array{T,3}     # Array with density of each fluid layer
         H :: Array{T,3}     # Array with rest height of each fluid layer
         U :: Array{T,3}     # Array with imposed constant zonal flow U(y) in each fluid layer
       eta :: Array{T,2}     # Array containing topographic PV
         μ :: T              # Linear bottom drag
         ν :: T              # Viscosity coefficient
        nν :: Int            # Hyperviscous order (nν=1 is plain old viscosity)
   calcFq! :: Function       # Function that calculates the forcing on QGPV q

  # derived params
        g′ :: Array{T,1}     # Array with the reduced gravity constants for each fluid interface
        Qx :: Array{T,3}     # Array containing x-gradient of PV due to eta in each fluid layer
        Qy :: Array{T,3}     # Array containing y-gradient of PV due to β, U, and eta in each fluid layer
         S :: Array{T,4}     # Array containing coeffients for getting PV from  streamfunction
      invS :: Array{T,4}     # Array containing coeffients for inverting PV to streamfunction
  rfftplan :: FFTW.rFFTWPlan{T,-1,false,3}  # rfft plan for FFTs
end

struct SingleLayerParams{T} <: BarotropicParams
  # prescribed params
         g :: T              # Gravitational constant
        f0 :: T              # Constant planetary vorticity
         β :: T              # Planetary vorticity y-gradient
         ρ :: T              # Fluid layer density
         H :: T              # Fluid layer rest height
         U :: Array{T,1}     # Imposed constant zonal flow U(y)
       eta :: Array{T,2}     # Array containing topographic PV
         μ :: T              # Linear bottom drag
         ν :: T              # Viscosity coefficient
        nν :: Int            # Hyperviscous order (nν=1 is plain old viscosity)
   calcFq! :: Function       # Function that calculates the forcing on QGPV q

  # derived params
        Qx :: Array{T,2}     # Array containing x-gradient of PV due to eta in each fluid layer
        Qy :: Array{T,2}     # Array containing meridional PV gradient due to β, U, and eta in each fluid layer
  rfftplan :: FFTW.rFFTWPlan{T,-1,false,2}  # rfft plan for FFTs
end

function Params(nlayers, g, f0, β, ρ, H, U::Array{T,2}, eta, μ, ν, nν, grid::AbstractGrid{T}; calcFq=nothingfunction, effort=FFTW.MEASURE) where T

   ny, nx = grid.ny , grid.nx
  nkr, nl = grid.nkr, grid.nl
   kr, l  = grid.kr , grid.l

  U = reshape(U, (1, ny, nlayers))

  # g′ = g*(ρ[2:nlayers]-ρ[1:nlayers-1]) ./ ρ[1:nlayers-1] # definition match PYQG
  g′ = g*(ρ[2:nlayers]-ρ[1:nlayers-1]) ./ ρ[2:nlayers] # correct definition

  Fm = @. f0^2 / ( g′*H[2:nlayers  ] )
  Fp = @. f0^2 / ( g′*H[1:nlayers-1] )

   ρ = reshape(ρ, (1,  1, nlayers))
   H = reshape(H, (1,  1, nlayers))

  Uyy = repeat(irfft( -l.^2 .* rfft(U, [1, 2]),  1, [1, 2]), outer=(1, 1, 1))
  Uyy = repeat(Uyy, outer=(nx, 1, 1))

  etah = rfft(eta)
  etax = irfft(im*kr.*etah, nx)
  etay = irfft(im*l .*etah, nx)

  Qx = zeros(nx, ny, nlayers)
  @views @. Qx[:, :, nlayers] += etax

  Qy = zeros(nx, ny, nlayers)
  Qy[:, :, 1] = @. β - Uyy[:, :, 1] - Fp[1]*( U[:, :, 2] - U[:, :, 1] )
  for j = 2:nlayers-1
    Qy[:, :, j] = @. β - Uyy[:, :, j] - Fp[j]*( U[:, :, j+1] - U[:, :, j] ) - Fm[j-1]*( U[:, :, j-1] - U[:, :, j] )
  end
  @views Qy[:, :, nlayers] = @. β - Uyy[:, :, nlayers] - Fm[nlayers-1]*( U[:, :, nlayers-1] - U[:, :, nlayers] )
  @views @. Qy[:, :, nlayers] += etay

  S = Array{T}(undef, (nkr, nl, nlayers, nlayers))
  calcS!(S, Fp, Fm, grid)

  invS = Array{T}(undef, (nkr, nl, nlayers, nlayers))
  calcinvS!(invS, Fp, Fm, grid)

  rfftplanlayered = plan_rfft(Array{T,3}(undef, grid.nx, grid.ny, nlayers), [1, 2]; flags=effort)

  if nlayers == 1
    SingleLayerParams{T}(g, f0, β, ρ, H, U, eta, μ, ν, nν, calcFq, Qx, Qy, grid.rfftplan)
  else
    Params{T}(nlayers, g, f0, β, ρ, H, U, eta, μ, ν, nν, calcFq, g′, Qx, Qy, S, invS, rfftplanlayered)
  end
end

function Params(nlayers, g, f0, β, ρ, H, U::Array{T,1}, eta, μ, ν, nν, grid::AbstractGrid{T}; calcFq=nothingfunction, effort=FFTW.MEASURE) where T
  U = reshape(U, (1, nlayers))
  U = repeat(U, outer=(grid.ny, 1))
  Params(nlayers, g, f0, β, ρ, H, U, eta, μ, ν, nν, grid; calcFq=calcFq, effort=effort)
end

# ---------
# Equations
# ---------

function hyperdissipation(ν, nν, Krsq, nkr, nl, nlayers, T)
  L = Array{Complex{T}}(undef, (nkr, nl, nlayers))
  @. L = -ν*Krsq^nν
  @views @. L[1, 1, :] = 0
  L
end

function LinearEquation(p, gr::AbstractGrid{T}) where T
  L = hyperdissipation(p.ν, p.nν, gr.Krsq, gr.nkr, gr.nl, p.nlayers, T)
  FourierFlows.Equation(L, calcNlinear!, gr)
end

function Equation(p, gr::AbstractGrid{T}) where T
  L = hyperdissipation(p.ν, p.nν, gr.Krsq, gr.nkr, gr.nl, p.nlayers, T)
  FourierFlows.Equation(L, calcN!, gr)
end

function Equation(p::SingleLayerParams, gr::AbstractGrid{T}) where T
  L = @. -p.μ - p.ν*gr.Krsq^p.nν + im*p.β*gr.kr*gr.invKrsq
  L[1, 1] = 0
  FourierFlows.Equation(L, calcN!, gr)
end

# ----
# Vars
# ----

abstract type BarotropicVars <: AbstractVars end

const physicalvars = [:q, :psi, :u, :v]
const fouriervars = [ Symbol(var, :h) for var in physicalvars ]
const forcedfouriervars = cat(fouriervars, [:Fqh], dims=1)

varspecs = cat(
  getfieldspecs(physicalvars, :(Array{T,3})),
  getfieldspecs(fouriervars, :(Array{Complex{T},3})),
  dims=1)

forcedvarspecs = cat(
  getfieldspecs(physicalvars, :(Array{T,3})),
  getfieldspecs(forcedfouriervars, :(Array{Complex{T},3})),
  dims=1)

singlelayervarsspecs = cat(
  getfieldspecs(physicalvars, :(Array{T,2})),
  getfieldspecs(fouriervars, :(Array{Complex{T},2})),
  dims=1)

# Construct Vars types
eval(varsexpression(:Vars, physicalvars, fouriervars))
eval(varsexpression(:ForcedVars, physicalvars, forcedfouriervars))
eval(varsexpression(:SingleLayerVars, singlelayervarsspecs; parent=:BarotropicVars, typeparams=:T))

"""
    Vars(g)

Returns the vars for unforced multi-layer QG problem with grid gr.
"""
function Vars(gr::AbstractGrid{T}, p) where T
  @zeros T (gr.nx, gr.ny, p.nlayers) q psi u v
  @zeros Complex{T} (gr.nkr, gr.nl, p.nlayers) qh psih uh vh
  Vars(q, psi, u, v, qh, psih, uh, vh)
end

"""
    ForcedVars(g)

Returns the vars for forced multi-layer QG problem with grid gr.
"""
function ForcedVars(gr::AbstractGrid{T}, p) where T
  v = Vars(gr, p)
  Fqh = zeros(Complex{T}, (gr.nkr, gr.nl, p.nlayers))
  ForcedVars(getfield.(Ref(v), fieldnames(typeof(v)))..., Fqh)
end

function SingleLayerVars(gr::AbstractGrid{T}) where T
  @zeros T (gr.nx, gr.ny) q psi u v
  @zeros Complex{T} (gr.nkr, gr.nl) qh psih uh vh
  SingleLayerParams(q, psi, u, v, qh, psih, uh, vh)
end

fwdtransform!(varh, var, p::AbstractParams) = mul!(varh, p.rfftplan, var)
invtransform!(var, varh, p::AbstractParams) = ldiv!(var, p.rfftplan, varh)

function streamfunctionfrompv!(psih, qh, invS, gr)
  for j=1:gr.nl, i=1:gr.nkr
    @views psih[i, j, :] .= invS[i, j, :, :]*qh[i, j, :]
  end
end

function pvfromstreamfunction!(qh, psih, S, gr)
  for j=1:gr.nl, i=1:gr.nkr
    @views qh[i, j, :] .= S[i, j, :, :]*psih[i, j, :]
  end
end

"""
    calcS!(S, Fp, Fm, gr)

Constructs the stretching matrix S that connects q and ψ: q_{k,l} = S * ψ_{k,l}.
"""
function calcS!(S, Fp, Fm, gr)
  F = Matrix(Tridiagonal(Fm, -([Fp; 0] + [0; Fm]), Fp))
  for n=1:gr.nl, m=1:gr.nkr
     k2 = gr.Krsq[m, n]
    Skl = -k2*I + F
    @views S[m, n, :, :] .= Skl
  end
  nothing
end

"""
    calcinvS!(S, Fp, Fm, g)

Constructs the inverse of the stretching matrix S that connects q and ψ:
ψ_{k,l} = invS * q_{k,l}.
"""
function calcinvS!(invS, Fp, Fm, gr)
  F = Matrix(Tridiagonal(Fm, -([Fp; 0] + [0; Fm]), Fp))
  for n=1:gr.nl, m=1:gr.nkr
    k2 = gr.Krsq[m, n]
    if k2 == 0
      k2 = 1
    end
    Skl = -k2*I + F
    @views invS[m, n, :, :] .= I/Skl
  end
  @views invS[1, 1, :, :] .= 0
  nothing
end

# -------
# Solvers
# -------

function calcN!(N, sol, t, cl, v, p, gr)
  calcN_advection!(N, sol, v, p, gr)
  @views @. N[:, :, p.nlayers] += p.μ*gr.Krsq*v.psih[:, :, p.nlayers]   # bottom linear drag
  addforcing!(N, sol, t, cl, v, p, gr)
  nothing
end

function calcNlinear!(N, sol, t, cl, v, p, gr)
  calcN_linearadvection!(N, sol, v, p, gr)
  @views @. N[:, :, p.nlayers] += p.μ*gr.Krsq*v.psih[:, :, p.nlayers]   # bottom linear drag
  addforcing!(N, sol, t, cl, v, p, gr)
  nothing
end

"""
    calcN_advection!(N, sol, v, p, g)

Calculates the advection term.
"""
function calcN_advection!(N, sol, v, p, gr)
  @. v.qh = sol

  streamfunctionfrompv!(v.psih, v.qh, p.invS, gr)

  @. v.uh = -im*gr.l *v.psih
  @. v.vh =  im*gr.kr*v.psih

  invtransform!(v.u, v.uh, p)
  @. v.u += p.U                         # add the imposed zonal flow U
  @. v.q = v.u*p.Qx
  fwdtransform!(v.uh, v.q, p)
  @. N = -v.uh                          # -(U+u)*∂Q/∂x

  invtransform!(v.v, v.vh, p)
  @. v.q = v.v*p.Qy
  fwdtransform!(v.vh, v.q, p)
  @. N -= v.vh                          # -v*∂Q/∂y

  invtransform!(v.q, v.qh, p)

  @. v.u *= v.q # u*q
  @. v.v *= v.q # v*q

  fwdtransform!(v.uh, v.u, p)
  fwdtransform!(v.vh, v.v, p)

  @. N -= im*gr.kr*v.uh + im*gr.l*v.vh    # -∂[(U+u)q]/∂x-∂[vq]/∂y

  nothing
end

"""
    calcN_linearadvection!(N, sol, v, p, g)

Calculates the advection term of the linearized equations.
"""
function calcN_linearadvection!(N, sol, v, p, gr)
  @. v.qh = sol

  streamfunctionfrompv!(v.psih, v.qh, p.invS, gr)

  @. v.uh = -im*gr.l *v.psih
  @. v.vh =  im*gr.kr*v.psih

  invtransform!(v.u, v.uh, p)
  @. v.u += p.U                        # add the imposed zonal flow U
  @. v.q = v.u*p.Qx
  fwdtransform!(v.uh, v.q, p)
  @. N = -v.uh                         # -(U+u)*∂Q/∂x

  invtransform!(v.v, v.vh, p)
  @. v.q = v.v*p.Qy
  fwdtransform!(v.vh, v.q, p)
  @. N -= v.vh                         # -v*∂Q/∂y

  invtransform!(v.q, v.qh, p)
  @. v.u = p.U
  @. v.u *= v.q # u*q

  fwdtransform!(v.uh, v.u, p)

  @. N -= im*gr.kr*v.uh                 # -∂[U*q]/∂x

  nothing
end

addforcing!(N, sol, t, cl, v::Vars, p, gr) = nothing

function addforcing!(N, sol, t, cl, v::ForcedVars, p, gr)
  p.calcFq!(v.Fqh, sol, t, cl, v, p, gr)
  @. N += v.Fqh
  nothing
end

# ----------------
# Helper functions
# ----------------

"""
    updatevars!(prob)

Update `prob.vars` using `prob.sol`.
"""
function updatevars!(prob)
  p, v, gr, sol = prob.params, prob.vars, prob.grid, prob.sol

  @. v.qh = sol
  streamfunctionfrompv!(v.psih, v.qh, p.invS, gr)
  @. v.uh = -im*gr.l*v.psih
  @. v.vh =  im*gr.kr*v.psih

  invtransform!(v.q, deepcopy(v.qh), p)
  invtransform!(v.psi, deepcopy(v.psih), p)
  invtransform!(v.u, deepcopy(v.uh), p)
  invtransform!(v.v, deepcopy(v.vh), p)
  nothing
end

"""
    set_q!(prob)

Set the solution `prob.sol` as the transform of `q` and updates variables.
"""
function set_q!(prob, q)
  p, v, sol = prob.params, prob.vars, prob.sol

  fwdtransform!(v.qh, q, p)
  @. v.qh[1, 1, :] = 0
  @. sol = v.qh

  updatevars!(prob)
  nothing
end

"""
    set_psi!(prob)

Set the solution `prob.sol` to correspond to a streamfunction `psi` and
updates variables.
"""
function set_psi!(prob, psi)
  p, v, gr = prob.params, prob.vars, prob.grid

  fwdtransform!(v.psih, psi, p)
  pvfromstreamfunction!(v.qh, v.psih, p.S, gr)
  invtransform!(v.q, v.qh, p)
  set_q!(prob, v.q)

  nothing
end

"""
    energies(prob)

Returns the kinetic energy of each fluid layer KE_1,...,KE_nlayers, and the
potential energy of each fluid interface PE_{3/2},...,PE_{nlayers-1/2}.
"""
function energies(prob)
  v, p, gr, sol = prob.vars, prob.params, prob.grid, prob.sol
  KE, PE = zeros(p.nlayers), zeros(p.nlayers-1)

  @. v.qh = sol
  streamfunctionfrompv!(v.psih, v.qh, p.invS, gr)

  @. v.uh = gr.Krsq * abs2(v.psih)
  for j=1:p.nlayers
    KE[j] = 1/(2*gr.Lx*gr.Ly)*parsevalsum(v.uh[:, :, j], gr)*p.H[j]/sum(p.H)
  end

  for j=1:p.nlayers-1
    PE[j] = 1/(2*gr.Lx*gr.Ly)*p.f0^2/p.g′[j]*parsevalsum(abs2.(v.psih[:, :, j+1].-v.psih[:, :, j]), gr)
  end

  KE, PE
end

"""
    fluxes(prob)

Returns the lateral eddy fluxes within each fluid layer
lateralfluxes_1,...,lateralfluxes_nlayers and also the vertical eddy fluxes for
each fluid interface verticalfluxes_{3/2},...,verticalfluxes_{nlayers-1/2}
"""
function fluxes(prob)
  v, p, gr, sol = prob.vars, prob.params, prob.grid, prob.sol
  lateralfluxes, verticalfluxes = zeros(p.nlayers), zeros(p.nlayers-1)

  updatevars!(prob)

  @. v.uh = -im*gr.l*v.uh
  invtransform!(v.u, v.uh, p)

  lateralfluxes = (sum(@. p.H*p.U*v.v*v.u; dims=(1,2)))[1, 1, :]
  lateralfluxes *= gr.dx*gr.dy/(gr.Lx*gr.Ly*sum(p.H))

  for j=1:p.nlayers-1
    verticalfluxes[j] = sum( @views @. p.f0^2/p.g′[j] * (p.U[: ,:, j] - p.U[:, :, j+1])*v.v[:, :, j+1]*v.psi[:, :, j] ; dims=(1,2) )[1]
    verticalfluxes[j] *= gr.dx*gr.dy/(gr.Lx*gr.Ly*sum(p.H))
  end

  lateralfluxes, verticalfluxes
end

end # module


text:
"""
Constant stomatal conductance for CO₂ struct.

# Arguments

- `g0`: intercept.
- `gs`: stomatal conductance.
- `gs_min = 0.001`: residual conductance. We consider the residual conductance being different
 from `g0` because in practice `g0` can be negative when fitting real-world data.

Then used as follows:
Gs = ConstantGs(0.0,0.1)
Gₛ = Gs.g0 + Gs.gs * A
"""
struct ConstantGs{T} <: AbstractGsModel
    g0::T
    gs::T
    gs_min::T
end

function ConstantGs(g0, gs, gs_min)
    ConstantGs(promote(g0, gs, gs_min)...)
end

ConstantGs(g0,gs) = ConstantGs(g0, gs, oftype(gs, 0.001))

function inputs(::ConstantGs)
    (:Gₛ,)
end

function outputs(::ConstantGs)
    (:Gₛ,)
end

Base.eltype(x::ConstantGs) = typeof(x).parameters[1]

"""
Constant stomatal closure

# Note

`meteo` is just declared here for compatibility with other formats of calls.
"""
function gs_closure(leaf::LeafModels{I,E,A,Gs,S}, meteo = missing) where {I,E,A,Gs <: ConstantGs,S}
    leaf.stomatal_conductance.gs - leaf.stomatal_conductance.g0
end

"""
Constant stomatal conductance for CO₂ (mol m-2 s-1).

# Note

`meteo` or `gs_mod` are just declared here for compatibility with the call from
photosynthesis (need a constant way of calling the functions).
"""
function gs(leaf::LeafModels{I,E,A,Gs,S}, gs_mod) where {I,E,A,Gs <: ConstantGs,S}
    leaf.stomatal_conductance.gs
end

function gs(leaf::LeafModels{I,E,A,Gs,S}, meteo::M) where {I,E,A,Gs <: ConstantGs,S,M <: AbstractAtmosphere}
    leaf.stomatal_conductance.gs
end

function gs!(leaf::LeafModels{I,E,A,Gs,S}, gs_mod) where {I,E,A,Gs <: ConstantGs,S}
    leaf.status.Gₛ = leaf.stomatal_conductance.gs
end

function gs!(leaf::LeafModels{I,E,A,Gs,S}, meteo::M) where {I,E,A,Gs <: ConstantGs,S,M <: AbstractAtmosphere}
    leaf.status.Gₛ = leaf.stomatal_conductance.gs
end


text:
function fibo(n::Int)
    if n < 2
        return n;
    end
    fibo(n-1) + fibo(n-2);
end

println(fibo(45))


text:
# Handy helpers

function include_everywhere(filepath)
    include(filepath) # Load on Node 1 first, triggering any precompile
    if nprocs() > 1
        fullpath = joinpath(@__DIR__, filepath)
        @sync for p in workers()
            @async remotecall_wait(include, p, fullpath)
        end
    end
end

function zip_word_sets(U, V)
    zipped = ""
    for (u, v) in zip(U, V)
        zipped *= u * " " * v * " "
    end
    return zipped
end

# Make a bunch of test "documents" filled with WEAT words
function make_test_docs(weat_word_set)
    SA = zip_word_sets(weat_word_set.S, weat_word_set.A)
    SB = zip_word_sets(weat_word_set.S, weat_word_set.B)
    TA = zip_word_sets(weat_word_set.T, weat_word_set.A)
    TB = zip_word_sets(weat_word_set.T, weat_word_set.B)
    return (SA=SA, SB=SB, TA=TA, TB=TB)
end

# Regex Helper
function extract(string, regex; trim=(0,0), cast=nothing)
    m = match(regex, string)
    tmp = nothing
    if (m != nothing)
        tmp = m.match[1+trim[1]:end-trim[2]]
        if (cast != nothing)
            tmp = parse(cast, tmp)
        end
    end
    return tmp
end

# Get file info from naming convention
function fileinfo(filepath)
    corpus = extract(filepath, r"C[0-9]+", trim=(1,0), cast=Int64)
    min_vocab = extract(filepath, r"V[0-9]+", trim=(1,0), cast=Int64)
    window = extract(filepath, r"W[0-9]+", trim=(1,0), cast=Int64)
    dimension = extract(filepath, r"D[0-9]+", trim=(1,0), cast=Int64)
    eta = extract(filepath, r"R[0-9]+.[0-9]+", trim=(1,0), cast=Float64)
    max_iters = extract(filepath, r"E[0-9]+", trim=(1,0), cast=Int64)
    seed = extract(filepath, r"S[0-9]+", trim=(1,0), cast=Int64)
    tmp = extract(filepath, r".[0-9]{3}.bin$", trim=(1, 4), cast=Int64)
    iters = tmp == nothing ? max_iters : tmp
    return (corpus=corpus, min_vocab=min_vocab, window=window, dimension=dimension,
            eta=eta, max_iters=max_iters, seed=seed, iters=iters)
end


text:
using Test
using TestSetExtensions
using FluxTraining
using Animations
using Colors
using FluxTraining: EpochEnd, LearningRate, protect, Protected,
    ProtectedException, Read, Write, Events, Phases, runstep, runepoch, epoch!, step!
using FluxTraining: getcallback, setcallbacks!, replacecallback!, removecallback!, addcallback!

import FluxTraining.Phases: Phase
using Flux: trainable
using Suppressor
using FluxTraining: CHECKS, SanityCheckException
using Zygote

include("./testdata.jl");


text:
function NegateDetector(arg0::EventDetector)
    return NegateDetector((EventDetector,), arg0)
end

function g(obj::NegateDetector, arg0::SpacecraftState)
    return jcall(obj, "g", jdouble, (SpacecraftState,), arg0)
end

function get_original(obj::NegateDetector)
    return jcall(obj, "getOriginal", EventDetector, ())
end

function init(obj::NegateDetector, arg0::SpacecraftState, arg1::AbsoluteDate)
    return jcall(obj, "init", void, (SpacecraftState, AbsoluteDate), arg0, arg1)
end


text:
module Consumers
    include("../../utils/Utils.jl")
    import .Utils

    using AMQPClient

    export base_consumer

    """ Функция базового консьюмера для подписки на очередь.

    channel: канал связи с сервисом очередей,
    message: сообщение от сервиса очередей,
    user_consumer_reaction: функция пользовательской реакции
    """
    function  base_consumer(channel::AMQPClient.MessageChannel, message::AMQPClient.Message, user_consumer_reaction::Function)::Bool
        # Логирование для режима отладки
        @debug "Параметры вызова функции:"
        @debug "> канал соединения с серверов RabbitMQ " channel
        @debug "> сообщение из очереди " message
        @debug "> функция консьюмера пользователя " user_consumer_reaction
        try
            @debug "Запускается пользовательская операция обработки реакции " user_consumer_reaction
            user_consumer_reaction(channel, message)
            AMQPClient.basic_ack(channel, message.delivery_tag)
            return true
        catch error
            @error Utils.Chronometer.message_with_time("Не удалось отреагировать") error
            throw(error)
        end
    end
end