text:
struct SpectralConv{T, W, I, O} <: AbstractSpectralOperator
    transform::T
    weights::W
    in_channels::I
    out_channels::O
end

"""
    SpectralConv(transform, channels; init = glorot_uniform)

Constructs a spectral convolution operator with spectral transform `trafo`, channels `channels`,
and initialization function `init`.

# Example
```
julia> trafo = FourierTransform(modes = (12, 3, 4, 12))
julia> channels = 1 => 4
julia> conv = SpectralConv(trafo, channels)
```
"""
function SpectralConv(
    transform::T,
    ch::Pair{S, S};
    init = Flux.glorot_uniform,
) where {T, S <: Int}
    # initialize optimizable weights
    FT = eltype(transform)
    in, out = ch
    scale = one(FT) / (in * out)

    # depending on the transform, the size of the coefficient
    # matrix may be different
    space_dims = size(transform)

    # depending on the transform, we may need to initialize
    # with real or complex numbers
    if FT <: Complex
        weights =
            init(in, space_dims..., out) + init(in, space_dims..., out) * im
    else
        weights = init(in, space_dims..., out)
    end
    weights = scale * weights

    return SpectralConv(transform, weights, in, out)
end

Flux.@functor SpectralConv

function (conv::SpectralConv)(x::AbstractArray)
    # TODO: (1) perhaps dispatch directly instead of calling ndims(conv)
    #       (2) perhaps something to infer size(x) at compile time
    c = OperatorFlux.forward(conv.transform, x)
    ct = OperatorFlux.truncate_modes(conv.transform, c)
    wc = OperatorFlux.tensor_contraction(conv.weights, ct, Val(ndims(conv)))
    wcp = OperatorFlux.pad_modes(conv.transform, wc, size(c)[2:(end - 1)])
    y = OperatorFlux.inverse(conv.transform, wcp)

    return y
end

struct SpectralCovariance{T, U, V, W, I, O} <: AbstractSpectralOperator
    transform::T
    weights_μ::U
    weights_d::V
    weights_w::W
    in_channels::I
    out_channels::O
end

"""
    SpectralCovariance(transform, rank, channels; init = glorot_uniform)

Constructs a spectral covariance operator with spectral transform `trafo`, low-rank matrix rank `rank`, channels `channels`,
and initialization function `init`.

# Example
```
julia> trafo = FourierTransform(modes = (12, 3, 4, 12))
julia> channels = 1 => 4
julia> rank = 3
julia> conv = SpectralCovariance(trafo, rank, channels)
```
"""
function SpectralCovariance(
    transform::T,
    rank::R,
    ch::Pair{S, S};
    init = Flux.glorot_uniform,
) where {T, R, S <: Int}
    FT = eltype(transform)
    in, out = ch
    scale = one(FT) / (in * out)

    # depending on the transform, we may need to initialize
    # with real or complex numbers
    if FT <: Complex
        weights_μ = init(in, out) + init(in, out) * im
        weights_d = init(in, out) + init(in, out) * im
        weights_w = init(in, rank, out) + init(in, rank, out) * im
    else
        weights_μ = init(in, out)
        weights_d = init(in, out)
        weights_w = init(in, rank, out)
    end
    weights_μ = scale * weights_μ
    weights_d = scale * weights_d
    weights_w = scale * weights_w

    return SpectralCovariance(
        transform,
        weights_μ,
        weights_d,
        weights_w,
        in,
        out,
    )
end

Flux.@functor SpectralCovariance

function (so::SpectralCovariance)(x::AbstractArray)
    N_space = length(size(x)[2:(end - 1)])

    # Bring input into truncated spectral representation
    c = OperatorFlux.forward(so.transform, x)
    ct = OperatorFlux.truncate_modes(so.transform, c)

    # paramterize the mean for the low-rank approximation
    μ = OperatorFlux.sparse_mean(so.weights_μ, ct, Val(N_space))
    μ = OperatorFlux.pad_modes(so.transform, μ, size(c)[2:(end - 1)])
    μ = OperatorFlux.inverse(so.transform, μ)

    # parameterize the diagonal matrix for the low-rank approximation
    D = OperatorFlux.sparse_mean(so.weights_d, ct, Val(N_space))
    D = OperatorFlux.pad_modes(so.transform, D, size(c)[2:(end - 1)])
    D = OperatorFlux.inverse(so.transform, D)
    D = @. log(1 + exp(D)) # keeping a stable positive attitude

    # parameterize the low-rank matrix for the low-rank approximation
    rank = size(so.weights_w)[end - 1]
    V = OperatorFlux.sparse_covariance(so.weights_w, ct, Val(N_space))
    V = OperatorFlux.pad_modes(so.transform, V, (size(c)[2:(end - 1)]..., rank))
    # we need to bring the array into a shape that is compatible with the
    # transform API
    store_v_size = size(V)
    V = reshape(
        V,
        (store_v_size[1:(end - 2)]..., prod(store_v_size[(end - 1):end]))...,
    )
    V = OperatorFlux.inverse(so.transform, V)
    V = reshape(V, (size(V)[1:(end - 1)]..., store_v_size[(end - 1):end]...)...)

    return μ, D, V
end

struct SpectralKernelOperator{L, C, A} <: AbstractSpectralOperator
    linear::L
    conv::C
    σ::A
end

"""
    SpectralKernelOperator(transform, channels, σ = identity; init = glorot_uniform)

Constructs a spectral kernel operator with spectral transform `trafo`, channels `channels`,
activation function `σ`, and initialization function `init`.

References:
Fourier Neural Operator for Parametric Partial Differential Equations - https://arxiv.org/abs/2010.08895

# Example
```
julia> trafo = FourierTransform(modes = (12, 3, 4, 12))
julia> channels = 1 => 4
julia> σ = relu
julia> conv = SpectralKernelOperator(trafo, channels, σ)
```
"""
function SpectralKernelOperator(
    transform::AbstractTransform,
    channels::Pair{S, S},
    σ = identity;
    init = Flux.glorot_uniform,
) where {S <: Int, N}
    linear = Dense(channels.first, channels.second)
    conv = SpectralConv(transform, channels, init = init)

    return SpectralKernelOperator(linear, conv, σ)
end

Flux.@functor SpectralKernelOperator

function (so::SpectralKernelOperator)(x)
    # The definition of the Spectral Kernel Operators include a linear transform as well.
    return so.σ.(so.linear(x) + so.conv(x))
end

# Base extensions
Base.ndims(so::AbstractSpectralOperator) = ndims(so.transform)

function Base.show(io::IO, so::SpectralConv)
    print(
        io,
        "SpectralConv($(so.in_channels) => $(so.out_channels), $(so.transform))",
    )
end

function Base.show(io::IO, so::SpectralCovariance)
    print(
        io,
        "SpectralConv($(so.in_channels) => $(so.out_channels), $(so.transform)))",
    )
end

function Base.show(io::IO, so::SpectralKernelOperator)
    print(
        io,
        "SpectralKernelOperator(" *
        "$(sco.conv.in_channel) => $(sco.conv.out_channel), " *
        "$(sco.conv.transform), " *
        "σ=$(string(sco.σ)) " *
        ")",
    )
end


text:
# get test data from EconDatasets
intRates = []
try
    intRates = dataset("Treasuries")
catch
    getDataset("Treasuries")
    intRates = dataset("Treasuries")
end

# create subset
dats = Date(1991, 1, 2):Date(2013, 12, 31)
intData = intRates[dats, 4]
names!(intData.vals, [:Y1])

# save subset as test data to make it available to R
fname = joinpath(Pkg.dir("Econometrics"), "test/data/no_git_intRates.csv")
writeTimedata(fname, intData)

ext: jl


text:
"""
    A = JopFilter(spc, responsetype, designmethod)

where `A` is a filter applied to a signal in `spc::JotSpace`, and built
using `responsetype` and `designmethod`.  The `responsetype` and
`designmethod` are described in `https://github.com/JuliaDSP/DSP.jl`.
The filter is applied along the fast dimension of the space.

# Examples

## 1D
```
using JetPackDSP, Jets
A = JopFilter(JetSpace(Float64,512), Highpass(10.0, fs=125), Butterworth(4))
d = A*rand(domain(A))
```

## 2D
```
using JetPackDSP, Jets
A = JopFilter(JetSpace(Float64,512,10), Highpass(10.0, fs=125), Butterworth(4))
d = A*rand(domain(A))
```
"""
function JopFilter(spc::JetSpace{T}, responsetype::FilterType, designmethod::FilterCoefficients) where {T}
	n = size(spc,1)
	tmp1 = zeros(T,n)
	tmp2 = zeros(T,n)
	tmp1[div(n,2)] = 1;
	tmp2 .= filtfilt(digitalfilter(responsetype, designmethod), tmp1)
	ztmp2 = fft(tmp2)
	filter = real.(abs.(ztmp2))
	JopLn(dom = spc, rng = spc, df! = JopFilter_df!, df′! = JopFilter_df′!, s = (filter=filter,))
end

export JopFilter

function JopFilter_df!(d::AbstractArray{T}, m::AbstractArray{T}; filter, kwargs...) where {T}
	n = size(m, 1)
    _d = reshape(d, n, :)
    _m = reshape(m, n, :)

    ztmp = zeros(Complex{T}, n)
    for k = 1:size(_d, 2)
        @inbounds begin
            ztmp .= fft(_m[:,k])
            ztmp .*= filter
            _d[:,k] .= real.(ifft(ztmp))
        end
    end
    d
end

function JopFilter_df′!(m::AbstractArray{T}, d::AbstractArray{T}; filter, kwargs...) where {T}
    n = size(m, 1)
    _d = reshape(d, n, :)
    _m = reshape(m, n, :)

    ztmp = zeros(Complex{T}, n)
    for k = 1:size(_d,2)
        @inbounds begin
            ztmp .= fft(_d[:,k])
            ztmp .*= filter
            _m[:,k] .= real.(ifft(ztmp))
        end
    end
    m
end


text:
"""
spin18.jl - sampling robustness
"""

WDIR = joinpath(@__DIR__, "../../")
include(joinpath(WDIR, "src", "spin", "spin.jl"))

using Altro
using HDF5
using LinearAlgebra
using RobotDynamics
using StaticArrays
using TrajectoryOptimization
const RD = RobotDynamics
const TO = TrajectoryOptimization

# paths
const EXPERIMENT_META = "spin"
const EXPERIMENT_NAME = "spin18"
const SAVE_PATH = joinpath(WDIR, "out", EXPERIMENT_META, EXPERIMENT_NAME)

# problem
const CONTROL_COUNT = 1
const STATE_COUNT = 2
const SAMPLE_COUNT = 8
const ASTATE_SIZE_BASE = STATE_COUNT * HDIM_ISO + 3 * CONTROL_COUNT
const ASTATE_SIZE = ASTATE_SIZE_BASE + SAMPLE_COUNT * HDIM_ISO
const ACONTROL_SIZE = CONTROL_COUNT
const INITIAL_STATE1 = [1., 0, 0, 0]
const INITIAL_STATE2 = [0., 1, 0, 0]
const INITIAL_STATE3 = [1., 0, 0, 1] ./ sqrt(2)
const INITIAL_STATE4 = [1., -1, 0, 0] ./ sqrt(2)
# state indices
const STATE1_IDX = 1:HDIM_ISO
const STATE2_IDX = STATE1_IDX[end] + 1:STATE1_IDX[end] + HDIM_ISO
const INTCONTROLS_IDX = STATE2_IDX[end] + 1:STATE2_IDX[end] + CONTROL_COUNT
const CONTROLS_IDX = INTCONTROLS_IDX[end] + 1:INTCONTROLS_IDX[end] + CONTROL_COUNT
const DCONTROLS_IDX = CONTROLS_IDX[end] + 1:CONTROLS_IDX[end] + CONTROL_COUNT
const S1_IDX = DCONTROLS_IDX[end] + 1:DCONTROLS_IDX[end] + HDIM_ISO
const S2_IDX = S1_IDX[end] + 1:S1_IDX[end] + HDIM_ISO
const S3_IDX = S2_IDX[end] + 1:S2_IDX[end] + HDIM_ISO
const S4_IDX = S3_IDX[end] + 1:S3_IDX[end] + HDIM_ISO
const S5_IDX = S4_IDX[end] + 1:S4_IDX[end] + HDIM_ISO
const S6_IDX = S5_IDX[end] + 1:S5_IDX[end] + HDIM_ISO
const S7_IDX = S6_IDX[end] + 1:S6_IDX[end] + HDIM_ISO
const S8_IDX = S7_IDX[end] + 1:S7_IDX[end] + HDIM_ISO
# control indices
const D2CONTROLS_IDX = 1:CONTROL_COUNT

# model
struct Model <: AbstractModel
    namp::Float64
end
@inline RD.state_dim(::Model) = ASTATE_SIZE
@inline RD.control_dim(::Model) = ACONTROL_SIZE

# This cost puts a gate error cost on
# the sample states and a LQR cost on the other terms.
# The hessian w.r.t the state and controls is constant.
struct Cost{N,M,T} <: TO.CostFunction
    Q::Diagonal{T, SVector{N,T}}
    R::Diagonal{T, SVector{M,T}}
    q::SVector{N, T}
    c::T
    hess_astate::Symmetric{T, SMatrix{N,N,T}}
    target_states::Array{SVector{HDIM_ISO, T}, 1}
    q_ss1::T
    q_ss2::T
    q_ss3::T
    q_ss4::T
end

function Cost(Q::Diagonal{T,SVector{N,T}}, R::Diagonal{T,SVector{M,T}},
              xf::SVector{N,T}, target_states::Array{SVector{HDIM_ISO}, 1},
              q_ss1::T, q_ss2::T, q_ss3::T, q_ss4::T) where {N,M,T}
    q = -Q * xf
    c = 0.5 * xf' * Q * xf
    hess_astate = zeros(N, N)
    # For reasons unknown to the author, throwing a -1 in front
    # of the gate error Hessian makes the cost function work.
    # This is strange, because the gate error Hessian has been
    # checked against autodiff.
    hess_state1 = -1 * q_ss1 * hessian_gate_error_iso2(target_states[1])
    hess_state2 = -1 * q_ss2 * hessian_gate_error_iso2(target_states[2])
    hess_state3 = -1 * q_ss3 * hessian_gate_error_iso2(target_states[3])
    hess_state4 = -1 * q_ss4 * hessian_gate_error_iso2(target_states[4])
    hess_astate[S1_IDX, S1_IDX] = hess_state1
    hess_astate[S2_IDX, S2_IDX] = hess_state2
    hess_astate[S3_IDX, S3_IDX] = hess_state3
    hess_astate[S4_IDX, S4_IDX] = hess_state4
    hess_astate[S5_IDX, S5_IDX] = hess_state1
    hess_astate[S6_IDX, S6_IDX] = hess_state2
    hess_astate[S7_IDX, S7_IDX] = hess_state3
    hess_astate[S8_IDX, S8_IDX] = hess_state4
    hess_astate += Q
    hess_astate = Symmetric(SMatrix{N, N}(hess_astate))
    return Cost{N,M,T}(Q, R, q, c, hess_astate, target_states, q_ss1, q_ss2, q_ss3, q_ss4)
end

@inline TO.state_dim(cost::Cost{N,M,T}) where {N,M,T} = N
@inline TO.control_dim(cost::Cost{N,M,T}) where {N,M,T} = M
@inline Base.copy(cost::Cost{N,M,T}) where {N,M,T} = Cost{N,M,T}(
    cost.Q, cost.R, cost.q, cost.c, cost.hess_astate,
    cost.target_states, cost.q_ss1, cost.q_ss2, cost.q_ss3, cost.q_ss4
)

@inline TO.stage_cost(cost::Cost{N,M,T}, astate::SVector{N}) where {N,M,T} = (
    0.5 * astate' * cost.Q * astate + cost.q'astate + cost.c
    + cost.q_ss1 * gate_error_iso2(astate, cost.target_states[1], S1_IDX[1] - 1)
    + cost.q_ss2 * gate_error_iso2(astate, cost.target_states[2], S2_IDX[1] - 1)
    + cost.q_ss3 * gate_error_iso2(astate, cost.target_states[3], S3_IDX[1] - 1)
    + cost.q_ss4 * gate_error_iso2(astate, cost.target_states[4], S4_IDX[1] - 1)
    + cost.q_ss1 * gate_error_iso2(astate, cost.target_states[1], S5_IDX[1] - 1)
    + cost.q_ss2 * gate_error_iso2(astate, cost.target_states[2], S6_IDX[1] - 1)
    + cost.q_ss3 * gate_error_iso2(astate, cost.target_states[3], S7_IDX[1] - 1)
    + cost.q_ss4 * gate_error_iso2(astate, cost.target_states[4], S8_IDX[1] - 1)
)

@inline TO.stage_cost(cost::Cost{N,M,T}, astate::SVector{N}, acontrol::SVector{M}) where {N,M,T} = (
    TO.stage_cost(cost, astate) + 0.5 * acontrol' * cost.R * acontrol
)

function TO.gradient!(E::TO.QuadraticCostFunction, cost::Cost{N,M,T}, astate::SVector{N,T}) where {N,M,T}
    E.q = (cost.Q * astate + cost.q + [
        @SVector zeros(ASTATE_SIZE_BASE);
        cost.q_ss1 * jacobian_gate_error_iso2(astate, cost.target_states[1], S1_IDX[1] - 1);
        cost.q_ss2 * jacobian_gate_error_iso2(astate, cost.target_states[2], S2_IDX[1] - 1);
        cost.q_ss3 * jacobian_gate_error_iso2(astate, cost.target_states[3], S3_IDX[1] - 1);
        cost.q_ss4 * jacobian_gate_error_iso2(astate, cost.target_states[4], S4_IDX[1] - 1);
        cost.q_ss1 * jacobian_gate_error_iso2(astate, cost.target_states[1], S5_IDX[1] - 1);
        cost.q_ss2 * jacobian_gate_error_iso2(astate, cost.target_states[2], S6_IDX[1] - 1);
        cost.q_ss3 * jacobian_gate_error_iso2(astate, cost.target_states[3], S7_IDX[1] - 1);
        cost.q_ss4 * jacobian_gate_error_iso2(astate, cost.target_states[4], S8_IDX[1] - 1);
    ])
    return false
end

function TO.gradient!(E::TO.QuadraticCostFunction, cost::Cost{N,M,T}, astate::SVector{N,T},
                      acontrol::SVector{M,T}) where {N,M,T}
    TO.gradient!(E, cost, astate)
    E.r = cost.R * acontrol
    E.c = 0
    return false
end

function TO.hessian!(E::TO.QuadraticCostFunction, cost::Cost{N,M,T}, astate::SVector{N,T}) where {N,M,T}
    E.Q = cost.hess_astate
    return true
end

function TO.hessian!(E::TO.QuadraticCostFunction, cost::Cost{N,M,T}, astate::SVector{N,T},
                     acontrol::SVector{M,T}) where {N,M,T}
    TO.hessian!(E, cost, astate)
    E.R = cost.R
    E.H .= 0
    return true
end

# dynamics
function RD.discrete_dynamics(::Type{RD.RK3}, model::Model, astate::StaticVector,
                              acontrols::StaticVector, time::Real, dt::Real) where {SC}
    camp = astate[CONTROLS_IDX[1]]
    negi_hc = camp * NEGI_H1_ISO
    h_prop = exp((FQ_NEGI_H0_ISO + negi_hc) * dt)
    state1 = h_prop * astate[STATE1_IDX]
    state2 = h_prop * astate[STATE2_IDX]
    intcontrols = astate[INTCONTROLS_IDX[1]] + dt * astate[CONTROLS_IDX[1]]
    controls = astate[CONTROLS_IDX[1]] + dt * astate[DCONTROLS_IDX[1]]
    dcontrols = astate[DCONTROLS_IDX[1]] + dt * acontrols[D2CONTROLS_IDX[1]]

    hp_prop = exp((FQ_NEGI_H0_ISO + (camp + model.namp) * NEGI_H1_ISO) * dt)
    hn_prop = exp((FQ_NEGI_H0_ISO + (camp - model.namp) * NEGI_H1_ISO) * dt)
    s1 = hp_prop * astate[S1_IDX]
    s2 = hp_prop * astate[S2_IDX]
    s3 = hp_prop * astate[S3_IDX]
    s4 = hp_prop * astate[S4_IDX]
    s5 = hn_prop * astate[S5_IDX]
    s6 = hn_prop * astate[S6_IDX]
    s7 = hn_prop * astate[S7_IDX]
    s8 = hn_prop * astate[S8_IDX]

    astate_ = [
        state1; state2; intcontrols; controls; dcontrols;
        s1; s2; s3; s4; s5; s6; s7; s8;
    ]

    return astate_
end

# main
function run_traj(;gate_type=xpiby2, evolution_time=56.8, solver_type=altro,
                  sqrtbp=false, integrator_type=rk3, qs=[1e0, 1e0, 1e0, 1e-1, 1e0, 1e0, 1e0, 1e0, 1e-1],
                  dt_inv=Int64(1e1), smoke_test=false, constraint_tol=1e-8, al_tol=1e-4,
                  pn_steps=2, max_penalty=1e11, verbose=true, save=true, max_iterations=Int64(2e5),
                  namp=NAMP_PREFACTOR)
    # model configuration
    model = Model(namp)
    n = state_dim(model)
    m = control_dim(model)
    t0 = 0.
    tf = evolution_time

    # initial state
    x0 = SVector{n}([
        INITIAL_STATE1;
        INITIAL_STATE2;
        zeros(3 * CONTROL_COUNT);
        repeat([INITIAL_STATE1; INITIAL_STATE2; INITIAL_STATE3; INITIAL_STATE4], 2);
    ])
    # target state
    gate_unitary = GT_GATE[gate_type]
    target_states = Array{SVector{HDIM_ISO}, 1}(undef, 4)
    target_states[1] = gate_unitary * INITIAL_STATE1
    target_states[2] = gate_unitary * INITIAL_STATE2
    target_states[3] = gate_unitary * INITIAL_STATE3
    target_states[4] = gate_unitary * INITIAL_STATE4
    xf = SVector{n}([
        target_states[1];
        target_states[2];
        zeros(3 * CONTROL_COUNT);
        repeat([target_states[1]; target_states[2];
                target_states[3]; target_states[4]], 2);
    ])

    # control amplitude constraint
    x_max = SVector{n}([
        fill(Inf, STATE_COUNT * HDIM_ISO);
        fill(Inf, CONTROL_COUNT);
        fill(MAX_CONTROL_NORM_0, 1); # control
        fill(Inf, CONTROL_COUNT);
        fill(Inf, SAMPLE_COUNT * HDIM_ISO);
    ])
    x_min = SVector{n}([
        fill(-Inf, STATE_COUNT * HDIM_ISO);
        fill(-Inf, CONTROL_COUNT);
        fill(-MAX_CONTROL_NORM_0, 1); # control
        fill(-Inf, CONTROL_COUNT);
        fill(-Inf, SAMPLE_COUNT * HDIM_ISO);
    ])
    # control amplitude constraint at boundary
    x_max_boundary = SVector{n}([
        fill(Inf, STATE_COUNT * HDIM_ISO);
        fill(Inf, CONTROL_COUNT);
        fill(0, 1); # control
        fill(Inf, CONTROL_COUNT);
        fill(Inf, SAMPLE_COUNT * HDIM_ISO);
    ])
    x_min_boundary = SVector{n}([
        fill(-Inf, STATE_COUNT * HDIM_ISO);
        fill(-Inf, CONTROL_COUNT);
        fill(0, 1); # control
        fill(-Inf, CONTROL_COUNT);
        fill(-Inf, SAMPLE_COUNT * HDIM_ISO);
    ])

    # initial trajectory
    dt = dt_inv^(-1)
    N = Int(floor(evolution_time * dt_inv)) + 1
    U0 = [SVector{m}([
        fill(1e-4, CONTROL_COUNT);
    ]) for k = 1:N-1]
    X0 = [SVector{n}([
        fill(NaN, n);
    ]) for k = 1:N]
    Z = Traj(X0, U0, dt * ones(N))

    # cost function
    Q = Diagonal(SVector{n}([
        fill(qs[1], STATE_COUNT * HDIM_ISO); # ψ1, ψ2
        fill(qs[2], 1); # ∫a
        fill(qs[3], 1); # a
        fill(qs[4], 1); # ∂a
        fill(0, SAMPLE_COUNT * HDIM_ISO);
    ]))
    Qf = Q * N
    R = Diagonal(SVector{m}([
        fill(qs[9], CONTROL_COUNT);
    ]))
    cost_k = Cost(Q, R, xf, target_states, qs[5], qs[6], qs[7], qs[8])
    cost_f = Cost(Qf, R, xf, target_states, N * qs[5], N * qs[6], N * qs[7], N * qs[8])
    objective = TO.Objective(cost_k, cost_f, N)

    # must satisfy control amplitude bound
    control_bnd = BoundConstraint(n, m, x_max=x_max, x_min=x_min)
    # must statisfy conrols start and end at 0
    control_bnd_boundary = BoundConstraint(n, m, x_max=x_max_boundary, x_min=x_min_boundary)
    # must reach target state, must have zero net flux
    target_astate_constraint = GoalConstraint(xf, [STATE1_IDX; STATE2_IDX; INTCONTROLS_IDX])
    # must obey unit norm
    norm_constraints = [NormConstraint(n, m, 1, TO.Equality(), idxs) for idxs in (
        STATE1_IDX, STATE2_IDX, S1_IDX, S2_IDX, S3_IDX, S4_IDX, S5_IDX, S6_IDX, S7_IDX, S8_IDX,
    )]
    constraints = ConstraintList(n, m, N)
    add_constraint!(constraints, control_bnd, 2:N-2)
    add_constraint!(constraints, control_bnd_boundary, N-1:N-1)
    add_constraint!(constraints, target_astate_constraint, N:N);
    for norm_constraint in norm_constraints
        add_constraint!(constraints, norm_constraint, 2:N-1)
    end

    # solve problem
    prob = Problem{IT_RDI[integrator_type]}(model, objective, constraints, x0, xf, Z, N, t0, evolution_time)
    solver = ALTROSolver(prob)
    verbose_pn = verbose ? true : false
    verbose_ = verbose ? 2 : 0
    projected_newton = solver_type == altro ? true : false
    constraint_tolerance = solver_type == altro ? constraint_tol : al_tol
    iterations_inner = smoke_test ? 1 : 300
    iterations_outer = smoke_test ? 1 : 30
    n_steps = smoke_test ? 1 : pn_steps
    set_options!(solver, square_root=sqrtbp, constraint_tolerance=constraint_tolerance,
                 projected_newton_tolerance=al_tol, n_steps=n_steps,
                 penalty_max=max_penalty, verbose_pn=verbose_pn, verbose=verbose_,
                 projected_newton=projected_newton, iterations_inner=iterations_inner,
                 iterations_outer=iterations_outer, iterations=max_iterations)
    Altro.solve!(solver)

    # post-process
    acontrols_raw = TO.controls(solver)
    acontrols_arr = permutedims(reduce(hcat, map(Array, acontrols_raw)), [2, 1])
    astates_raw = TO.states(solver)
    astates_arr = permutedims(reduce(hcat, map(Array, astates_raw)), [2, 1])
    Q_raw = Array(Q)
    Q_arr = [Q_raw[i, i] for i in 1:size(Q_raw)[1]]
    Qf_raw = Array(Qf)
    Qf_arr = [Qf_raw[i, i] for i in 1:size(Qf_raw)[1]]
    R_raw = Array(R)
    R_arr = [R_raw[i, i] for i in 1:size(R_raw)[1]]
    cidx_arr = Array(CONTROLS_IDX)
    d2cidx_arr = Array(D2CONTROLS_IDX)
    cmax = TO.max_violation(solver)
    cmax_info = TO.findmax_violation(TO.get_constraints(solver))
    iterations_ = Altro.iterations(solver)

    result = Dict(
        "acontrols" => acontrols_arr,
        "controls_idx" => cidx_arr,
        "d2controls_dt2_idx" => d2cidx_arr,
        "evolution_time" => evolution_time,
        "astates" => astates_arr,
        "Q" => Q_arr,
        "Qf" => Qf_arr,
        "R" => R_arr,
        "cmax" => cmax,
        "cmax_info" => cmax_info,
        "dt" => dt,
        "sample_count" => SAMPLE_COUNT,
        "solver_type" => Integer(solver_type),
        "sqrtbp" => Integer(sqrtbp),
        "max_penalty" => max_penalty,
        "constraint_tol" => constraint_tol,
        "al_tol" => al_tol,
        "gate_type" => Integer(gate_type),
        "save_type" => Integer(jl),
        "integrator_type" => Integer(integrator_type),
        "iterations" => iterations_,
        "max_iterations" => max_iterations,
    )

    # save
    if save
        save_file_path = generate_file_path("h5", EXPERIMENT_NAME, SAVE_PATH)
        println("Saving this optimization to $(save_file_path)")
        h5open(save_file_path, "cw") do save_file
            for key in keys(result)
                write(save_file, key, result[key])
            end
        end
        result["save_file_path"] = save_file_path
    end

    return result
end

max_stars_repo_name: SchusterLab/rbqoc

max_issues_repo_name: SchusterLab/rbqoc

max_forks_repo_name: SchusterLab/rbqoc


text:
mutable struct POMDPEnv{T,Ts,Ta}
    model::T
    state::Ts
    actions::Ta
    actionspace::DiscreteSpace
end
POMDPEnv(model) = POMDPEnv(model, initialstate(model, ENV_RNG), actions(model),
                           DiscreteSpace(n_actions(model), 1))

mutable struct MDPEnv{T,Ts,Ta}
    model::T
    state::Ts
    actions::Ta
    actionspace::DiscreteSpace
end
MDPEnv(model) = MDPEnv(model, initialstate(model, ENV_RNG), actions(model),
                       DiscreteSpace(n_actions(model), 1))

actionspace(env::Union{MDPEnv, POMDPEnv}) = env.actionspace
observationindex(env, o) = Int64(o) + 1

function interact!(env::POMDPEnv, action)
    s, o, r = generate_sor(env.model, env.state, env.actions[action], ENV_RNG)
    env.state = s
    (observation = observationindex(env.model, o),
     reward = r,
     isdone = isterminal(env.model, s))
end
function reset!(env::Union{POMDPEnv, MDPEnv})
    env.state = initialstate(env.model, ENV_RNG)
    (observation = env.state,)
end
function getstate(env::POMDPEnv)
    (observation = observationindex(env.model, generate_o(env.model, env.state, ENV_RNG)),
     isdone = isterminal(env.model, env.state))
end

function interact!(env::MDPEnv, action)
    s = rand(ENV_RNG, transition(env.model, env.state, env.actions[action]))
    r = POMDPs.reward(env.model, env.state, env.actions[action])
    env.state = s
    (observation = stateindex(env.model, s),
     reward = r,
     isdone = isterminal(env.model, s))
end
function getstate(env::MDPEnv)
    (observation = stateindex(env.model, env.state),
     isdone = isterminal(env.model, env.state))
end


text:
### A Pluto.jl notebook ###
# v0.12.4

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 39b0d07a-1a70-11eb-2379-f7db05b0309f
begin
	using Plots, PlutoUI
	import CSV, DataFrames, Dates
end

# ╔═╡ 71c4521a-1a6f-11eb-248f-2b3c6a1fd63c
md"""
# Covid Data Exploration
"""

# ╔═╡ 7729bea0-1a71-11eb-1bcd-35dbb24fad87
begin
	url = "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv";
	download(url, "covid_global_cases.csv")
	data = CSV.read("covid_global_cases.csv") |>
		df -> DataFrames.rename(df, 1 => "province", 2 => "country")
	rm("covid_global_cases.csv")
end

# ╔═╡ f3c46c64-1a6f-11eb-046f-ef04da3cc23c
begin
	date_labels = names(data)[5:end]
	date_format = Dates.DateFormat("m/d/y")
	dates = Dates.Date.(date_labels, date_format) .+ Dates.Year(2000)
end

# ╔═╡ f3a333fa-1a6f-11eb-38cf-63eb06a4e8a0
countries = unique(data[:, :country]);

# ╔═╡ f38336e0-1a6f-11eb-35d1-cd48853a70d7
data_by_country = data |>
	df -> DataFrames.groupby(df, :country) |>
	df -> DataFrames.combine(df, (date_labels .=> sum .=> date_labels));

# ╔═╡ f366cfc8-1a6f-11eb-2e5e-b77251d64915
data_by_country;

# ╔═╡ c6bc690a-1a75-11eb-0348-8f17b8ac5298
md"""
## Per Country Helper Functions
"""

# ╔═╡ f32a56b0-1a6f-11eb-3d9b-018e87833269
function make_first_plot()
	base_plot = plot()
	xlabel!(base_plot, "Dates")
	ylabel!(base_plot, "Confirmed Cases")
	title!(base_plot, "Confirmed Covid Cases")

	base_plot
end

# ╔═╡ 47aefd04-1a77-11eb-242b-91ca5d5e1c43
function get_country_data(country)
	data_by_country |>
		df -> filter(:country => val -> val == country, df) |>
		df -> df[1, 2:end] |>
		df -> convert(Vector, df)
end

# ╔═╡ f30b96e2-1a6f-11eb-35eb-bbbc58be9c7b
function add_plot!(target_plot, country, index)
	country_data = get_country_data(country)
	plot!(target_plot,
	 	dates[1:index], country_data[1:index],
	 	xticks    = dates[1:12:end],
	 	xrotation = 45,
	 	legend = :topleft,
	 	label = country,
	 	lw = 3)
end

# ╔═╡ 2f89793e-1a79-11eb-0324-0d0e6365b50e
@bind selected_countries MultiSelect(
	[ctry => ctry for ctry in countries],
	default = ["US", "Brazil", "China", "Italy", "India"])

# ╔═╡ 7ca306b8-1a79-11eb-3f41-df9344c279f3
@bind dates_index Slider(1:length(dates), default = 40)

# ╔═╡ f2ef88a2-1a6f-11eb-1d8f-efd4f136e8b7
begin
	output_plot = make_first_plot()
	for country in selected_countries
		add_plot!(output_plot, country, dates_index)
	end
	output_plot
end

# ╔═╡ f2d07866-1a6f-11eb-1e71-fbf0edb3e3b9
data_by_country |>
	df -> df[df.country .== "US", 2:end] |>
	df -> convert(Vector, df[1, :])

# ╔═╡ f2969d30-1a6f-11eb-377c-5db9eed4a9c6

# ╔═╡ f2790d7e-1a6f-11eb-013a-ddf620c7d40a

# ╔═╡ f25dd522-1a6f-11eb-3e31-79412281c1cb

# ╔═╡ f24049f6-1a6f-11eb-232a-cd95846ccc90

# ╔═╡ f2043556-1a6f-11eb-2ff8-3945b6356950

# ╔═╡ Cell order:
# ╟─71c4521a-1a6f-11eb-248f-2b3c6a1fd63c
# ╠═39b0d07a-1a70-11eb-2379-f7db05b0309f
# ╠═7729bea0-1a71-11eb-1bcd-35dbb24fad87
# ╠═f3c46c64-1a6f-11eb-046f-ef04da3cc23c
# ╠═f3a333fa-1a6f-11eb-38cf-63eb06a4e8a0
# ╠═f38336e0-1a6f-11eb-35d1-cd48853a70d7
# ╠═f366cfc8-1a6f-11eb-2e5e-b77251d64915
# ╟─c6bc690a-1a75-11eb-0348-8f17b8ac5298
# ╠═f32a56b0-1a6f-11eb-3d9b-018e87833269
# ╠═47aefd04-1a77-11eb-242b-91ca5d5e1c43
# ╠═f30b96e2-1a6f-11eb-35eb-bbbc58be9c7b
# ╟─2f89793e-1a79-11eb-0324-0d0e6365b50e
# ╠═7ca306b8-1a79-11eb-3f41-df9344c279f3
# ╠═f2ef88a2-1a6f-11eb-1d8f-efd4f136e8b7
# ╠═f2d07866-1a6f-11eb-1e71-fbf0edb3e3b9
# ╠═f2969d30-1a6f-11eb-377c-5db9eed4a9c6
# ╠═f2790d7e-1a6f-11eb-013a-ddf620c7d40a
# ╠═f25dd522-1a6f-11eb-3e31-79412281c1cb
# ╠═f24049f6-1a6f-11eb-232a-cd95846ccc90
# ╠═f2043556-1a6f-11eb-2ff8-3945b6356950

max_stars_repo_name: martin-fabbri/colab-notebooks

max_issues_repo_name: martin-fabbri/colab-notebooks

max_forks_repo_name: martin-fabbri/colab-notebooks


text:
"""
Performs a simple single-threaded Simulated Annealing optimization.

  Args:
    f (Function): Objective function.
    x (Array): Initial design variable values.
    t (Function): A function representing the annealing schedule (temperature).
    k_max (Integer): Maximum iterations.

  Returns:
    Tuple(Array, Float): A tuple representing the minimizers and minimum for
      the given objective function.
"""
function simulated_annealing(f, x, t, k_max)
    y = f(x)
    x_best, y_best = x, y
    for k in 1 : k_max
        x′ = x + randn(length(x))
        y′ = f(x′)
        Δy = y′ - y
        if Δy ≤ 0 || rand() < exp(-Δy/t(k))
            x, y = x′, y′
        end
        if y′ < y_best
            x_best, y_best = x′, y′
        end
    end
    return x_best, y_best
end

"""
Performs a simple multithreaded Simulated Annealing optimization.

  Args:
    f (Function): Objective function.
    x (Array): Initial design variable values.
    t (Function): A function representing the annealing schedule (temperature).
    k_max (Integer): Maximum iterations.

  Returns:
    Tuple(Array, Float): A tuple representing the minimizers and minimum for
      the given objective function.
"""
function simulated_annealing_p(f, x, t, k_max)
    y = f(x)
    x_best, y_best = x, y
    # create a lock to use in reading/writing current best solution
    lk = ReentrantLock()

    # split max iterations evenly over threads
    Threads.@threads for k in 1 : k_max
        x′ = x + randn(length(x))
        y′ = f(x′)
        Δy = y′ - y
        if Δy ≤ 0 || rand() < exp(-Δy/t(k))
            x, y = x′, y′
        end

        # acquire lock to potentially update current best while avoiding race conditions
        lock(lk) do
            if y′ < y_best
                x_best, y_best = x′, y′
            end
        end
    end
    return x_best, y_best
end

max_stars_repo_name: camirmas/distributed-ga

max_issues_repo_name: camirmas/distributed-ga

max_forks_repo_name: camirmas/distributed-ga


text:
abstract type PersistentMap{K, V} <: AbstractDict{K, V} end

struct NotFound end

struct PersistentArrayMap{K, V} <: PersistentMap{K, V}
    kvs::Vector{Pair{K, V}}
end
PersistentArrayMap{K, V}() where {K, V} =
    PersistentArrayMap{K, V}(Pair{K, V}[])
PersistentArrayMap(kvs::(Union{Tuple{K, V}, Pair{K, V}})...) where {K, V} =
    PersistentArrayMap{K, V}(Pair{K, V}[Pair(k, v) for (k, v) in kvs])
PersistentArrayMap(; kwargs...) = PersistentArrayMap(kwargs...)

Base.isequal(m1::PersistentArrayMap, m2::PersistentArrayMap) =
    isequal(Set(m1.kvs), Set(m2.kvs))
==(m1::PersistentArrayMap, m2::PersistentArrayMap) =
    Set(m1.kvs) == Set(m2.kvs)

Base.length(m::PersistentArrayMap)  = length(m.kvs)
Base.isempty(m::PersistentArrayMap) = length(m) == 0

findkeyidx(m::PersistentArrayMap, k) = findfirst(kv -> kv[1] == k, m.kvs)

function _get(m::PersistentArrayMap, k, default, hasdefault::Bool)
    for kv in m.kvs
        kv[1] == k && return kv[2]
    end
    hasdefault ? default : default()
end

Base.get(m::PersistentArrayMap, k) =
    _get(m, k, ()->error("key not found: $k"), false)
Base.get(m::PersistentArrayMap, k, default) =
    _get(m, k, default, true)
Base.getindex(m::PersistentArrayMap, k) = get(m, k)

Base.haskey(m::PersistentArrayMap, k) = get(m, k, NotFound()) != NotFound()

function assoc(m::PersistentArrayMap{K, V}, k, v) where {K, V}
    idx = findkeyidx(m, k)
    idx === nothing && return PersistentArrayMap{K, V}(push!(m.kvs[1:end], Pair{K,V}(k, v)))

    kvs = m.kvs[1:end]
    kvs[idx] = Pair{K,V}(k, v)
    PersistentArrayMap{K, V}(kvs)
end

function dissoc(m::PersistentArrayMap{K, V}, k) where {K, V}
    idx = findkeyidx(m, k)
    idx === nothing && return m

    kvs = m.kvs[1:end]
    splice!(kvs, idx)
    PersistentArrayMap{K, V}(kvs)
end

function Base.iterate(m::PersistentArrayMap, i = 1)
    if i > length(m)
        return nothing
    else
        return (m.kvs[i], i + 1)
    end
end

Base.map(f::( Union{DataType, Function}), m::PersistentArrayMap) =
    PersistentArrayMap([f(kv) for kv in m]...)

Base.show(io::IO, ::MIME"text/plain", m::PersistentArrayMap{K, V}) where {K, V} =
    print(io, "Persistent{$K, $V}$(m.kvs)")

# Persistent Hash Maps
# ====================

struct PersistentHashMap{K, V} <: PersistentMap{K, V}
    trie::SparseBitmappedTrie{PersistentArrayMap{K, V}}
    length::Int
end
PersistentHashMap{K, V}() where {K, V} =
    PersistentHashMap{K, V}(SparseNode(PersistentArrayMap{K, V}), 0)

function PersistentHashMap(itr)
    if length(itr) == 0
        return PersistentHashMap()
    end
    K, V = typejoin(map(typeof, itr)...).types
    m = PersistentHashMap{K, V}()
    for (k, v) in itr
        m = assoc(m, k, v)
    end
    m
end

function PersistentHashMap(kvs::(Tuple{Any, Any})...)
    PersistentHashMap([kvs...])
end

function PersistentHashMap(kvs::(Pair)...)
    PersistentHashMap([kvs...])
end

function PersistentHashMap(; kwargs...)
    isempty(kwargs) ?
    PersistentHashMap{Any, Any}() :
    PersistentHashMap(kwargs...)
end

Base.length(m::PersistentHashMap) = m.length
Base.isempty(m::PersistentHashMap) = length(m) == 0

zipd(x,y) = map(p -> p[1] => p[2], zip(x,y))
Base.isequal(m1::PersistentHashMap, m2::PersistentHashMap) =
    length(m1) == length(m2) && all(x -> isequal(x...), zipd(m1, m2))

tup_eq(x) = x[1] == x[2]
==(m1::PersistentHashMap, m2::PersistentHashMap) =
    length(m1) == length(m2) && all(x -> x[1] == x[2], zipd(m1, m2))

function _update(f::Function, m::PersistentHashMap{K, V}, key) where {K, V}
    keyhash = reinterpret(Int, hash(key))
    arraymap = get(m.trie, keyhash, PersistentArrayMap{K, V}())
    newmap = f(arraymap)
    newtrie, _ = update(m.trie, keyhash, newmap)
    PersistentHashMap{K, V}(newtrie,
                            m.length + (length(newmap) < length(arraymap) ? -1 :
                                        length(newmap) > length(arraymap) ? 1 :
                                        0))
end

function assoc(m::PersistentHashMap{K, V}, key, value) where {K, V}
    _update(m, key) do arraymap
        assoc(arraymap, key, value)
    end
end

function dissoc(m::PersistentHashMap, key)
    _update(m, key) do arraymap
        dissoc(arraymap, key)
    end
end

function Base.getindex(m::PersistentHashMap, key)
    val = get(m.trie, reinterpret(Int, hash(key)), NotFound())
    (val === NotFound()) && error("key not found")
    val[key]
end

Base.get(m::PersistentHashMap, key) = m[key]
function Base.get(m::PersistentHashMap, key, default)
    val = get(m.trie, reinterpret(Int, hash(key)), NotFound())
    (val === NotFound()) && return default
    val[key]
end

function Base.haskey(m::PersistentHashMap, key)
    val = get(m.trie, reinterpret(Int, hash(key)), NotFound())
    (val != NotFound()) && haskey(val, key)
end

function Base.iterate(m::PersistentHashMap)
    trie_iter_result = iterate(m.trie)
    if trie_iter_result === nothing
        return nothing
    else
        arrmap, triestate = trie_iter_result
        kvs = arrmap.kvs
        return iterate(m, (kvs, triestate))
    end
end

function Base.iterate(m::PersistentHashMap, (kvs, triestate))
    if isempty(kvs) && isempty(triestate)
        return nothing
    else
        if isempty(kvs)
            arrmap, triestate = iterate(m.trie, triestate)
            return iterate(m, (arrmap.kvs, triestate))
        else
            return (kvs[1], (kvs[2:end], triestate))
        end
    end
end

function Base.map(f::( Union{Function, DataType}), m::PersistentHashMap)
    PersistentHashMap([f(kv) for kv in m]...)
end

function Base.filter(f::Function, m::PersistentHashMap{K, V}) where {K, V}
    arr = Array{Pair{K, V},1}()
    for el in m
        f(el) && push!(arr, el)
    end
    isempty(arr) ? PersistentHashMap{K, V}() : PersistentHashMap(arr...)
end

# Suppress ambiguity warning while allowing merging with array
function _merge(d::PersistentHashMap, others...)
    acc = d
    for other in others
        for (k, v) in other
            acc = assoc(acc, k, v)
        end
    end
    acc
end

# This definition suppresses ambiguity warning
Base.merge(d::PersistentHashMap, others::AbstractDict...) =
    _merge(d, others...)
Base.merge(d::PersistentHashMap, others...) =
    _merge(d, others...)

 function Base.show(io::IO, ::MIME"text/plain", m::PersistentHashMap{K, V}) where {K, V}
    print(io, "Persistent{$K, $V}[")
    print(io, join(["$k => $v" for (k, v) in m], ", "))
    print(io, "]")
end


text:
open("titanic.jl", "w") do io
    println(io, "DataFrame(")
    for n in names(df)
        println(io, "    ", repr(n), " => ", repr(df[:, n]), ",")
    end
    println(io, ")")
end


text:
include("triangle.jl");
include("trianglev2.jl");

adj_mtx_file = "../../../data/A_adj.tsv"
inc_mtx_file = "../../../data/A_inc.tsv"

T1 = @elapsed N1 = triangle(adj_mtx_file)
T2 = @elapsed N2 = trianglev2(adj_mtx_file)

println("\nimplementation I  : time = ", T1, " triangle count = ", N1);
println("implementation II : time = ", T2, " triangle count = ", N2);

#######################################################
# Graph Challenge Benchmark
# Developer : Dr. Siddharth Samsi (ssamsi@mit.edu)
#
# MIT
#########################################################
# (c) <2017> Massachusetts Institute of Technology
########################################################

max_stars_repo_name: bkmgit/GraphChallenge

max_issues_repo_name: bkmgit/GraphChallenge

max_forks_repo_name: bkmgit/GraphChallenge


text:
using Paxos
using Paxos.Transports.Common
using Paxos.Transports.Memory
using Paxos.Transports.TCP
using Paxos.Utils

using Logging

using Test

function greetingResponder(greeting, response, messenger)
    for message in receivedMessages(messenger)
        sendMessage(messenger, "$response: $message")
    end
end

function testRoundtrip(transport, address, greeting, response)
    bounded(5) do
        worked = false
        @debug "Creating listener"
        listenerTask = listener(transport, address) do messenger
            greetingResponder(greeting, response, messenger)
        end
        @debug "Pausing for listener to be alive"
        sleep(2)
        try
            @debug "Listener created"
            @debug "Connecting"
            connection(transport, address) do messenger
                @debug "Connected"
                @debug "Sending test message"
                sendMessage(messenger, greeting)
                @debug "Sent"
                answer = take!(receivedMessages(messenger))
                worked = (answer == "$response: $greeting")
            end
        catch ex
            worked = false
            @error "Error during roundtrip" exception = (ex, stacktrace(catch_backtrace()))
        finally
            close(listenerTask)
        end
        worked
    end
end

function calls(fn, transport, addressesAndtimeouts, overallTimeout, greeting)
    bounded(overallTimeout) do
        listenerTasks = []
        messengers = Array{Messenger, 1}()
        nextResponse = 1
        try
            for addressAndTimeout in addressesAndtimeouts
                address, timeout = addressAndTimeout
                response = nextResponse
                nextResponse += 1
                listenerTask = listener(transport, address) do messenger
                    sleep(timeout)
                    greetingResponder(greeting, response, messenger)
                end
                push!(listenerTasks, listenerTask)
                sleep(0.5)
                messenger = messengerTo(transport, address)
                push!(messengers, messenger)
            end
            responses = gcall(messengers, overallTimeout, greeting)
            fn(responses)
        catch ex
            @error "Exception caught making calls" exception=(ex, stacktrace(catch_backtrace()))
            rethrow()
        finally
            closeAll(listenerTasks)
            closeAll(messengers)
        end
    end
end

function testGoodCalls(transport)
    worked = false
    addressesAndTimeouts = [
        (("localhost", 8001),1),
        (("localhost", 8002),1),
        (("localhost", 8003),1)
    ]
    calls(transport, addressesAndTimeouts, 3, "hello") do responses
        expected = [
            "1: hello",
            "2: hello",
            "3: hello"
        ]
        worked = (sort(responses) == expected)
    end
    worked
end

function testBadCalls(transport)
    worked = false
    addressesAndTimeouts = [
        (("localhost", 8004),0.1),
        (("localhost", 8005),2),
        (("localhost", 8006),0.1)
    ]
    calls(transport, addressesAndTimeouts, 2, "hello") do responses
        expected = [
            "1: hello",
            "3: hello"
        ]
        worked = (sort(responses) == expected)
    end
    worked
end


text:
# ---
# title: 924. Minimize Malware Spread
# id: problem924
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Depth-first Search, Union Find
# link: <https://leetcode.com/problems/minimize-malware-spread/description/>
# hidden: true
# ---
#
# In a network of nodes, each node `i` is directly connected to another node `j`
# if and only if `graph[i][j] = 1`.
#
# Some nodes `initial` are initially infected by malware.  Whenever two nodes
# are directly connected and at least one of those two nodes is infected by
# malware, both nodes will be infected by malware.  This spread of malware will
# continue until no more nodes can be infected in this manner.
#
# Suppose `M(initial)` is the final number of nodes infected with malware in the
# entire network, after the spread of malware stops.
#
# We will remove one node from the initial list.  Return the node that if
# removed, would minimize `M(initial)`.  If multiple nodes could be removed to
# minimize `M(initial)`, return such a node with the smallest index.
#
# Note that if a node was removed from the `initial` list of infected nodes, it
# may still be infected later as a result of the malware spread.
#
#
#
# **Example 1:**
#
#
#
#     Input: graph = [[1,1,0],[1,1,0],[0,0,1]], initial = [0,1]
#     Output: 0
#
#
# **Example 2:**
#
#
#
#     Input: graph = [[1,0,0],[0,1,0],[0,0,1]], initial = [0,2]
#     Output: 0
#
#
# **Example 3:**
#
#
#
#     Input: graph = [[1,1,1],[1,1,1],[1,1,1]], initial = [1,2]
#     Output: 1
#
#
#
#
# **Note:**
#
#   * `1 < graph.length = graph[0].length <= 300`
#   * `0 <= graph[i][j] == graph[j][i] <= 1`
#   * `graph[i][i] == 1`
#   * `1 <= initial.length <= graph.length`
#   * `0 <= initial[i] < graph.length`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using CyclicArrays, OrdinaryDiffEq, IndividualDisplacements

function cyclicarray_example()
    𝑃=cyclicarray_setup()

    x0=collect(0:20:490).%100 .+ 170
    y0=reshape(repeat(collect(0:20:80),1,5)',25) .+ 160
    uInit=Float64.([x0';y0'])

    N=1500
    tspan = (0.0,N*𝑃.dt)

    prob = ODEProblem(dxy_dt_CyclicArray,uInit,tspan,𝑃)
    #sol = solve(prob,Tsit5(),reltol=1e-5,abstol=1e-5)
    sol = solve(prob,Euler(),dt=𝑃.dt)

    x=[sol[1,i,end] for i in 1:length(x0)]
    y=[sol[2,i,end] for i in 1:length(x0)]

    return x,y
end

function cyclicarray_setup()
    nx=ny=100
    nfaces=1;

    faces=zeros(nfaces,2,2,4);
    faces[1,1,1,:]=[1,1,2,0];
    faces[1,1,2,:]=[1,1,1,0];
    faces[1,2,1,:]=[1,2,2,0];
    faces[1,2,2,:]=[1,2,1,0];

    g=CyclicArray(faces);

    ind=range(1, 360, length = 360)
    x = y = ind/180*pi
    f(y,x) = sin(x) + sin(y)

    X = repeat(reshape(x, 1, :), length(y), 1)
    Y = repeat(y, 1, length(x))

    F=f.(Y,X)

    phi=CyclicArray(F,g)
    xg=CyclicArray(repeat(ind, 1, length(x))',g)
    yg=CyclicArray(repeat(ind, 1, length(x)),g)
    u=diff(phi,dims=2)
    v=-diff(phi,dims=1)

    (u = u, v = v, xg = xg, yg = yg, dt = 100)
end


text:
@testset "4.5.4.1 (a+b sec)^m (A+B sec+C sec^2)" begin
    (A, B, C, b, c, d, e, f, m, x, ) = @variables A B C b c d e f m x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(b*sec(e+f*x))^m*(A+B*sec(e+f*x)+C*sec(e+f*x)^2)=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(b*sec(e+f*x))^m*(A+C*sec(e+f*x)^2)=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*sec(e+f*x)^m*(A+C*sec(e+f*x)^2)=#

    @test_int [sec(c + d*x)^6*(A + C*sec(c + d*x)^2), x, 3, ((7*A + 6*C)*tan(c + d*x))/(7*d) + (C*sec(c + d*x)^6*tan(c + d*x))/(7*d) + (2*(7*A + 6*C)*tan(c + d*x)^3)/(21*d) + ((7*A + 6*C)*tan(c + d*x)^5)/(35*d)]
    @test_int [sec(c + d*x)^5*(A + C*sec(c + d*x)^2), x, 4, ((6*A + 5*C)*atanh(sin(c + d*x)))/(16*d) + ((6*A + 5*C)*sec(c + d*x)*tan(c + d*x))/(16*d) + ((6*A + 5*C)*sec(c + d*x)^3*tan(c + d*x))/(24*d) + (C*sec(c + d*x)^5*tan(c + d*x))/(6*d)]
    @test_int [sec(c + d*x)^4*(A + C*sec(c + d*x)^2), x, 3, ((5*A + 4*C)*tan(c + d*x))/(5*d) + (C*sec(c + d*x)^4*tan(c + d*x))/(5*d) + ((5*A + 4*C)*tan(c + d*x)^3)/(15*d)]
    @test_int [sec(c + d*x)^3*(A + C*sec(c + d*x)^2), x, 3, ((4*A + 3*C)*atanh(sin(c + d*x)))/(8*d) + ((4*A + 3*C)*sec(c + d*x)*tan(c + d*x))/(8*d) + (C*sec(c + d*x)^3*tan(c + d*x))/(4*d)]
    @test_int [sec(c + d*x)^2*(A + C*sec(c + d*x)^2), x, 3, ((3*A + 2*C)*tan(c + d*x))/(3*d) + (C*sec(c + d*x)^2*tan(c + d*x))/(3*d)]
    @test_int [sec(c + d*x)^1*(A + C*sec(c + d*x)^2), x, 2, ((2*A + C)*atanh(sin(c + d*x)))/(2*d) + (C*sec(c + d*x)*tan(c + d*x))/(2*d)]
    @test_int [sec(c + d*x)^0*(A + C*sec(c + d*x)^2), x, 3, A*x + (C*tan(c + d*x))/d]
    @test_int [cos(c + d*x)^1*(A + C*sec(c + d*x)^2), x, 2, (C*atanh(sin(c + d*x)))/d + (A*sin(c + d*x))/d]
    @test_int [cos(c + d*x)^2*(A + C*sec(c + d*x)^2), x, 2, (1/2)*(A + 2*C)*x + (A*cos(c + d*x)*sin(c + d*x))/(2*d)]
    @test_int [cos(c + d*x)^3*(A + C*sec(c + d*x)^2), x, 3, ((A + C)*sin(c + d*x))/d - (A*sin(c + d*x)^3)/(3*d)]
    @test_int [cos(c + d*x)^4*(A + C*sec(c + d*x)^2), x, 3, (1/8)*(3*A + 4*C)*x + ((3*A + 4*C)*cos(c + d*x)*sin(c + d*x))/(8*d) + (A*cos(c + d*x)^3*sin(c + d*x))/(4*d)]
    @test_int [cos(c + d*x)^5*(A + C*sec(c + d*x)^2), x, 4, ((A + C)*sin(c + d*x))/d - ((2*A + C)*sin(c + d*x)^3)/(3*d) + (A*sin(c + d*x)^5)/(5*d)]
    @test_int [cos(c + d*x)^6*(A + C*sec(c + d*x)^2), x, 4, (1/16)*(5*A + 6*C)*x + ((5*A + 6*C)*cos(c + d*x)*sin(c + d*x))/(16*d) + ((5*A + 6*C)*cos(c + d*x)^3*sin(c + d*x))/(24*d) + (A*cos(c + d*x)^5*sin(c + d*x))/(6*d)]

    @test_int [sec(c + d*x)^m*(-C*m/(m + 1) + C*sec(c + d*x)^2), x, 1, (C*sec(c + d*x)^(1 + m)*sin(c + d*x))/(d*(1 + m))]
    @test_int [sec(c + d*x)^m*(A - A*(m + 1)/m*sec(c + d*x)^2), x, 1, -((A*sec(c + d*x)^(1 + m)*sin(c + d*x))/(d*m))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(b*sec(e+f*x))^(m/2)*(A+C*sec(e+f*x)^2)=#

    @test_int [(A + C*sec(c + d*x)^2)*(b*sec(c + d*x))^(5/2), x, 4, (2*b^2*(7*A + 5*C)*sqrt(cos(c + d*x))*Elliptic.F((1/2)*(c + d*x), 2)*sqrt(b*sec(c + d*x)))/(21*d) + (2*b*(7*A + 5*C)*(b*sec(c + d*x))^(3/2)*sin(c + d*x))/(21*d) + (2*C*(b*sec(c + d*x))^(5/2)*tan(c + d*x))/(7*d)]
    @test_int [(A + C*sec(c + d*x)^2)*(b*sec(c + d*x))^(3/2), x, 4, -((2*b^2*(5*A + 3*C)*Elliptic.E((1/2)*(c + d*x), 2))/(5*d*sqrt(cos(c + d*x))*sqrt(b*sec(c + d*x)))) + (2*b*(5*A + 3*C)*sqrt(b*sec(c + d*x))*sin(c + d*x))/(5*d) + (2*C*(b*sec(c + d*x))^(3/2)*tan(c + d*x))/(5*d)]
    @test_int [(A + C*sec(c + d*x)^2)*(b*sec(c + d*x))^(1/2), x, 3, (2*(3*A + C)*sqrt(cos(c + d*x))*Elliptic.F((1/2)*(c + d*x), 2)*sqrt(b*sec(c + d*x)))/(3*d) + (2*C*sqrt(b*sec(c + d*x))*tan(c + d*x))/(3*d)]
    @test_int [(A + C*sec(c + d*x)^2)/(b*sec(c + d*x))^(1/2), x, 3, (2*(A - C)*Elliptic.E((1/2)*(c + d*x), 2))/(d*sqrt(cos(c + d*x))*sqrt(b*sec(c + d*x))) + (2*C*tan(c + d*x))/(d*sqrt(b*sec(c + d*x)))]
    @test_int [(A + C*sec(c + d*x)^2)/(b*sec(c + d*x))^(3/2), x, 3, (2*(A + 3*C)*sqrt(cos(c + d*x))*Elliptic.F((1/2)*(c + d*x), 2)*sqrt(b*sec(c + d*x)))/(3*b^2*d) + (2*A*tan(c + d*x))/(3*d*(b*sec(c + d*x))^(3/2))]
    @test_int [(A + C*sec(c + d*x)^2)/(b*sec(c + d*x))^(5/2), x, 3, (2*(3*A + 5*C)*Elliptic.E((1/2)*(c + d*x), 2))/(5*b^2*d*sqrt(cos(c + d*x))*sqrt(b*sec(c + d*x))) + (2*A*tan(c + d*x))/(5*d*(b*sec(c + d*x))^(5/2))]
    @test_int [(A + C*sec(c + d*x)^2)/(b*sec(c + d*x))^(7/2), x, 4, (2*(5*A + 7*C)*sqrt(cos(c + d*x))*Elliptic.F((1/2)*(c + d*x), 2)*sqrt(b*sec(c + d*x)))/(21*b^4*d) + (2*(5*A + 7*C)*sin(c + d*x))/(21*b^3*d*sqrt(b*sec(c + d*x))) + (2*A*tan(c + d*x))/(7*d*(b*sec(c + d*x))^(7/2))]
    @test_int [(A + C*sec(c + d*x)^2)/(b*sec(c + d*x))^(9/2), x, 4, (2*(7*A + 9*C)*Elliptic.E((1/2)*(c + d*x), 2))/(15*b^4*d*sqrt(cos(c + d*x))*sqrt(b*sec(c + d*x))) + (2*(7*A + 9*C)*sin(c + d*x))/(45*b^3*d*(b*sec(c + d*x))^(3/2)) + (2*A*tan(c + d*x))/(9*d*(b*sec(c + d*x))^(9/2))]

    @test_int [(3 + 3*sec(c + d*x)^2)/sec(c + d*x)^(1/2), x, 1, (6*sqrt(sec(c + d*x))*sin(c + d*x))/d]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*sec(e+f*x)^m*(A+C*sec(e+f*x)^2)*when*A*(m+1)+C*m=0=#

    @test_int [sec(e + f*x)^m*(m - (m + 1)*sec(e + f*x)^2), x, 1, -((sec(e + f*x)^(1 + m)*sin(e + f*x))/f)]

    @test_int [sec(e + f*x)^5*(5 - 6*sec(e + f*x)^2), x, 1, -((sec(e + f*x)^5*tan(e + f*x))/f)]
    @test_int [sec(e + f*x)^4*(4 - 5*sec(e + f*x)^2), x, 1, -((sec(e + f*x)^4*tan(e + f*x))/f)]
    @test_int [sec(e + f*x)^3*(3 - 4*sec(e + f*x)^2), x, 1, -((sec(e + f*x)^3*tan(e + f*x))/f)]
    @test_int [sec(e + f*x)^2*(2 - 3*sec(e + f*x)^2), x, 1, -((sec(e + f*x)^2*tan(e + f*x))/f)]
    @test_int [sec(e + f*x)^1*(1 - 2*sec(e + f*x)^2), x, 1, -((sec(e + f*x)*tan(e + f*x))/f)]
    @test_int [sec(e + f*x)^0*(0 - 1*sec(e + f*x)^2), x, 2, -(tan(e + f*x)/f)]
    @test_int [cos(e + f*x)^1*(-1 - 0*sec(e + f*x)^2), x, 1, -(sin(e + f*x)/f)]
    @test_int [cos(e + f*x)^2*(-2 + 1*sec(e + f*x)^2), x, 1, -((cos(e + f*x)*sin(e + f*x))/f)]
    @test_int [cos(e + f*x)^3*(-3 + 2*sec(e + f*x)^2), x, 1, -((cos(e + f*x)^2*sin(e + f*x))/f)]
    @test_int [cos(e + f*x)^4*(-4 + 3*sec(e + f*x)^2), x, 1, -((cos(e + f*x)^3*sin(e + f*x))/f)]
    @test_int [cos(e + f*x)^5*(-5 + 4*sec(e + f*x)^2), x, 1, -((cos(e + f*x)^4*sin(e + f*x))/f)]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(b*sec(e+f*x))^m*(B*sec(e+f*x)+C*sec(e+f*x)^2)=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*sec(e+f*x)^m*(B*sec(e+f*x)+C*sec(e+f*x)^2)=#

    @test_int [sec(c + d*x)^3*(B*sec(c + d*x) + C*sec(c + d*x)^2), x, 7, (3*C*atanh(sin(c + d*x)))/(8*d) + (B*tan(c + d*x))/d + (3*C*sec(c + d*x)*tan(c + d*x))/(8*d) + (C*sec(c + d*x)^3*tan(c + d*x))/(4*d) + (B*tan(c + d*x)^3)/(3*d)]
    @test_int [sec(c + d*x)^2*(B*sec(c + d*x) + C*sec(c + d*x)^2), x, 6, (B*atanh(sin(c + d*x)))/(2*d) + (C*tan(c + d*x))/d + (B*sec(c + d*x)*tan(c + d*x))/(2*d) + (C*tan(c + d*x)^3)/(3*d)]
    @test_int [sec(c + d*x)^1*(B*sec(c + d*x) + C*sec(c + d*x)^2), x, 6, (C*atanh(sin(c + d*x)))/(2*d) + (B*tan(c + d*x))/d + (C*sec(c + d*x)*tan(c + d*x))/(2*d)]
    @test_int [sec(c + d*x)^0*(B*sec(c + d*x) + C*sec(c + d*x)^2), x, 4, (B*atanh(sin(c + d*x)))/d + (C*tan(c + d*x))/d]
    @test_int [cos(c + d*x)^1*(B*sec(c + d*x) + C*sec(c + d*x)^2), x, 4, B*x + (C*atanh(sin(c + d*x)))/d]
    @test_int [cos(c + d*x)^2*(B*sec(c + d*x) + C*sec(c + d*x)^2), x, 3, C*x + (B*sin(c + d*x))/d]
    @test_int [cos(c + d*x)^3*(B*sec(c + d*x) + C*sec(c + d*x)^2), x, 5, (B*x)/2 + (C*sin(c + d*x))/d + (B*cos(c + d*x)*sin(c + d*x))/(2*d)]
    @test_int [cos(c + d*x)^4*(B*sec(c + d*x) + C*sec(c + d*x)^2), x, 6, (C*x)/2 + (B*sin(c + d*x))/d + (C*cos(c + d*x)*sin(c + d*x))/(2*d) - (B*sin(c + d*x)^3)/(3*d)]
    @test_int [cos(c + d*x)^5*(B*sec(c + d*x) + C*sec(c + d*x)^2), x, 7, (3*B*x)/8 + (C*sin(c + d*x))/d + (3*B*cos(c + d*x)*sin(c + d*x))/(8*d) + (B*cos(c + d*x)^3*sin(c + d*x))/(4*d) - (C*sin(c + d*x)^3)/(3*d)]
    @test_int [cos(c + d*x)^6*(B*sec(c + d*x) + C*sec(c + d*x)^2), x, 7, (3*C*x)/8 + (B*sin(c + d*x))/d + (3*C*cos(c + d*x)*sin(c + d*x))/(8*d) + (C*cos(c + d*x)^3*sin(c + d*x))/(4*d) - (2*B*sin(c + d*x)^3)/(3*d) + (B*sin(c + d*x)^5)/(5*d)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(b*sec(e+f*x))^(m/2)*(B*sec(e+f*x)+C*sec(e+f*x)^2)=#

    @test_int [(B*sec(c + d*x) + C*sec(c + d*x)^2)*(b*sec(c + d*x))^(3/2), x, 10, -((6*b^2*C*Elliptic.E((1/2)*(c + d*x), 2))/(5*d*sqrt(cos(c + d*x))*sqrt(b*sec(c + d*x)))) + (2*b*B*sqrt(cos(c + d*x))*Elliptic.F((1/2)*(c + d*x), 2)*sqrt(b*sec(c + d*x)))/(3*d) + (6*b*C*sqrt(b*sec(c + d*x))*sin(c + d*x))/(5*d) + (2*B*(b*sec(c + d*x))^(3/2)*sin(c + d*x))/(3*d) + (2*C*(b*sec(c + d*x))^(5/2)*sin(c + d*x))/(5*b*d)]
    @test_int [(B*sec(c + d*x) + C*sec(c + d*x)^2)*(b*sec(c + d*x))^(1/2), x, 9, -((2*b*B*Elliptic.E((1/2)*(c + d*x), 2))/(d*sqrt(cos(c + d*x))*sqrt(b*sec(c + d*x)))) + (2*C*sqrt(cos(c + d*x))*Elliptic.F((1/2)*(c + d*x), 2)*sqrt(b*sec(c + d*x)))/(3*d) + (2*B*sqrt(b*sec(c + d*x))*sin(c + d*x))/d + (2*C*(b*sec(c + d*x))^(3/2)*sin(c + d*x))/(3*b*d)]
    @test_int [(B*sec(c + d*x) + C*sec(c + d*x)^2)/(b*sec(c + d*x))^(1/2), x, 8, -((2*C*Elliptic.E((1/2)*(c + d*x), 2))/(d*sqrt(cos(c + d*x))*sqrt(b*sec(c + d*x)))) + (2*B*sqrt(cos(c + d*x))*Elliptic.F((1/2)*(c + d*x), 2)*sqrt(b*sec(c + d*x)))/(b*d) + (2*C*sqrt(b*sec(c + d*x))*sin(c + d*x))/(b*d)]
    @test_int [(B*sec(c + d*x) + C*sec(c + d*x)^2)/(b*sec(c + d*x))^(3/2), x, 7, (2*B*Elliptic.E((1/2)*(c + d*x), 2))/(b*d*sqrt(cos(c + d*x))*sqrt(b*sec(c + d*x))) + (2*C*sqrt(cos(c + d*x))*Elliptic.F((1/2)*(c + d*x), 2)*sqrt(b*sec(c + d*x)))/(b^2*d)]
    @test_int [(B*sec(c + d*x) + C*sec(c + d*x)^2)/(b*sec(c + d*x))^(5/2), x, 8, (2*C*Elliptic.E((1/2)*(c + d*x), 2))/(b^2*d*sqrt(cos(c + d*x))*sqrt(b*sec(c + d*x))) + (2*B*sqrt(cos(c + d*x))*Elliptic.F((1/2)*(c + d*x), 2)*sqrt(b*sec(c + d*x)))/(3*b^3*d) + (2*B*sin(c + d*x))/(3*b^2*d*sqrt(b*sec(c + d*x)))]
    @test_int [(B*sec(c + d*x) + C*sec(c + d*x)^2)/(b*sec(c + d*x))^(7/2), x, 9, (6*B*Elliptic.E((1/2)*(c + d*x), 2))/(5*b^3*d*sqrt(cos(c + d*x))*sqrt(b*sec(c + d*x))) + (2*C*sqrt(cos(c + d*x))*Elliptic.F((1/2)*(c + d*x), 2)*sqrt(b*sec(c + d*x)))/(3*b^4*d) + (2*B*sin(c + d*x))/(5*b^2*d*(b*sec(c + d*x))^(3/2)) + (2*C*sin(c + d*x))/(3*b^3*d*sqrt(b*sec(c + d*x)))]
    @test_int [(B*sec(c + d*x) + C*sec(c + d*x)^2)/(b*sec(c + d*x))^(9/2), x, 10, (6*C*Elliptic.E((1/2)*(c + d*x), 2))/(5*b^4*d*sqrt(cos(c + d*x))*sqrt(b*sec(c + d*x))) + (10*B*sqrt(cos(c + d*x))*Elliptic.F((1/2)*(c + d*x), 2)*sqrt(b*sec(c + d*x)))/(21*b^5*d) + (2*B*sin(c + d*x))/(7*b^2*d*(b*sec(c + d*x))^(5/2)) + (2*C*sin(c + d*x))/(5*b^3*d*(b*sec(c + d*x))^(3/2)) + (10*B*sin(c + d*x))/(21*b^4*d*sqrt(b*sec(c + d*x)))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(b*sec(e+f*x))^m*(A+B*sec(e+f*x)+C*sec(e+f*x)^2)=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*sec(e+f*x)^m*(A+B*sec(e+f*x)+C*sec(e+f*x)^2)=#

    @test_int [sec(c + d*x)^4*(A + B*sec(c + d*x) + C*sec(c + d*x)^2), x, 7, (3*B*atanh(sin(c + d*x)))/(8*d) + ((5*A + 4*C)*tan(c + d*x))/(5*d) + (3*B*sec(c + d*x)*tan(c + d*x))/(8*d) + (B*sec(c + d*x)^3*tan(c + d*x))/(4*d) + (C*sec(c + d*x)^4*tan(c + d*x))/(5*d) + ((5*A + 4*C)*tan(c + d*x)^3)/(15*d)]
    @test_int [sec(c + d*x)^3*(A + B*sec(c + d*x) + C*sec(c + d*x)^2), x, 6, ((4*A + 3*C)*atanh(sin(c + d*x)))/(8*d) + (B*tan(c + d*x))/d + ((4*A + 3*C)*sec(c + d*x)*tan(c + d*x))/(8*d) + (C*sec(c + d*x)^3*tan(c + d*x))/(4*d) + (B*tan(c + d*x)^3)/(3*d)]
    @test_int [sec(c + d*x)^2*(A + B*sec(c + d*x) + C*sec(c + d*x)^2), x, 6, (B*atanh(sin(c + d*x)))/(2*d) + ((3*A + 2*C)*tan(c + d*x))/(3*d) + (B*sec(c + d*x)*tan(c + d*x))/(2*d) + (C*sec(c + d*x)^2*tan(c + d*x))/(3*d)]
    @test_int [sec(c + d*x)^1*(A + B*sec(c + d*x) + C*sec(c + d*x)^2), x, 5, ((2*A + C)*atanh(sin(c + d*x)))/(2*d) + (B*tan(c + d*x))/d + (C*sec(c + d*x)*tan(c + d*x))/(2*d)]
    @test_int [sec(c + d*x)^0*(A + B*sec(c + d*x) + C*sec(c + d*x)^2), x, 4, A*x + (B*atanh(sin(c + d*x)))/d + (C*tan(c + d*x))/d]
    @test_int [cos(c + d*x)^1*(A + B*sec(c + d*x) + C*sec(c + d*x)^2), x, 4, B*x + (C*atanh(sin(c + d*x)))/d + (A*sin(c + d*x))/d]
    @test_int [cos(c + d*x)^2*(A + B*sec(c + d*x) + C*sec(c + d*x)^2), x, 4, (1/2)*(A + 2*C)*x + (B*sin(c + d*x))/d + (A*cos(c + d*x)*sin(c + d*x))/(2*d)]
    @test_int [cos(c + d*x)^3*(A + B*sec(c + d*x) + C*sec(c + d*x)^2), x, 6, (B*x)/2 + ((A + C)*sin(c + d*x))/d + (B*cos(c + d*x)*sin(c + d*x))/(2*d) - (A*sin(c + d*x)^3)/(3*d)]
    @test_int [cos(c + d*x)^4*(A + B*sec(c + d*x) + C*sec(c + d*x)^2), x, 6, (1/8)*(3*A + 4*C)*x + (B*sin(c + d*x))/d + ((3*A + 4*C)*cos(c + d*x)*sin(c + d*x))/(8*d) + (A*cos(c + d*x)^3*sin(c + d*x))/(4*d) - (B*sin(c + d*x)^3)/(3*d)]
    @test_int [cos(c + d*x)^5*(A + B*sec(c + d*x) + C*sec(c + d*x)^2), x, 8, (3*B*x)/8 + ((A + C)*sin(c + d*x))/d + (3*B*cos(c + d*x)*sin(c + d*x))/(8*d) + (B*cos(c + d*x)^3*sin(c + d*x))/(4*d) - ((2*A + C)*sin(c + d*x)^3)/(3*d) + (A*sin(c + d*x)^5)/(5*d)]
    @test_int [cos(c + d*x)^6*(A + B*sec(c + d*x) + C*sec(c + d*x)^2), x, 7, (1/16)*(5*A + 6*C)*x + (B*sin(c + d*x))/d + ((5*A + 6*C)*cos(c + d*x)*sin(c + d*x))/(16*d) + ((5*A + 6*C)*cos(c + d*x)^3*sin(c + d*x))/(24*d) + (A*cos(c + d*x)^5*sin(c + d*x))/(6*d) - (2*B*sin(c + d*x)^3)/(3*d) + (B*sin(c + d*x)^5)/(5*d)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(b*sec(e+f*x))^(m/2)*(A+B*sec(e+f*x)+C*sec(e+f*x)^2)=#

    @test_int [(A + B*sec(c + d*x) + C*sec(c + d*x)^2)*(b*sec(c + d*x))^(3/2), x, 8, -((2*b^2*(5*A + 3*C)*Elliptic.E((1/2)*(c + d*x), 2))/(5*d*sqrt(cos(c + d*x))*sqrt(b*sec(c + d*x)))) + (2*b*B*sqrt(cos(c + d*x))*Elliptic.F((1/2)*(c + d*x), 2)*sqrt(b*sec(c + d*x)))/(3*d) + (2*b*(5*A + 3*C)*sqrt(b*sec(c + d*x))*sin(c + d*x))/(5*d) + (2*B*(b*sec(c + d*x))^(3/2)*sin(c + d*x))/(3*d) + (2*C*(b*sec(c + d*x))^(3/2)*tan(c + d*x))/(5*d)]
    @test_int [(A + B*sec(c + d*x) + C*sec(c + d*x)^2)*(b*sec(c + d*x))^(1/2), x, 7, -((2*b*B*Elliptic.E((1/2)*(c + d*x), 2))/(d*sqrt(cos(c + d*x))*sqrt(b*sec(c + d*x)))) + (2*(3*A + C)*sqrt(cos(c + d*x))*Elliptic.F((1/2)*(c + d*x), 2)*sqrt(b*sec(c + d*x)))/(3*d) + (2*B*sqrt(b*sec(c + d*x))*sin(c + d*x))/d + (2*C*sqrt(b*sec(c + d*x))*tan(c + d*x))/(3*d)]
    @test_int [(A + B*sec(c + d*x) + C*sec(c + d*x)^2)/(b*sec(c + d*x))^(1/2), x, 6, (2*(A - C)*Elliptic.E((1/2)*(c + d*x), 2))/(d*sqrt(cos(c + d*x))*sqrt(b*sec(c + d*x))) + (2*B*sqrt(cos(c + d*x))*Elliptic.F((1/2)*(c + d*x), 2)*sqrt(b*sec(c + d*x)))/(b*d) + (2*C*tan(c + d*x))/(d*sqrt(b*sec(c + d*x)))]
    @test_int [(A + B*sec(c + d*x) + C*sec(c + d*x)^2)/(b*sec(c + d*x))^(3/2), x, 6, (2*B*Elliptic.E((1/2)*(c + d*x), 2))/(b*d*sqrt(cos(c + d*x))*sqrt(b*sec(c + d*x))) + (2*(A + 3*C)*sqrt(cos(c + d*x))*Elliptic.F((1/2)*(c + d*x), 2)*sqrt(b*sec(c + d*x)))/(3*b^2*d) + (2*A*tan(c + d*x))/(3*d*(b*sec(c + d*x))^(3/2))]
    @test_int [(A + B*sec(c + d*x) + C*sec(c + d*x)^2)/(b*sec(c + d*x))^(5/2), x, 7, (2*(3*A + 5*C)*Elliptic.E((1/2)*(c + d*x), 2))/(5*b^2*d*sqrt(cos(c + d*x))*sqrt(b*sec(c + d*x))) + (2*B*sqrt(cos(c + d*x))*Elliptic.F((1/2)*(c + d*x), 2)*sqrt(b*sec(c + d*x)))/(3*b^3*d) + (2*B*sin(c + d*x))/(3*b^2*d*sqrt(b*sec(c + d*x))) + (2*A*tan(c + d*x))/(5*d*(b*sec(c + d*x))^(5/2))]
    @test_int [(A + B*sec(c + d*x) + C*sec(c + d*x)^2)/(b*sec(c + d*x))^(7/2), x, 8, (6*B*Elliptic.E((1/2)*(c + d*x), 2))/(5*b^3*d*sqrt(cos(c + d*x))*sqrt(b*sec(c + d*x))) + (2*(5*A + 7*C)*sqrt(cos(c + d*x))*Elliptic.F((1/2)*(c + d*x), 2)*sqrt(b*sec(c + d*x)))/(21*b^4*d) + (2*B*sin(c + d*x))/(5*b^2*d*(b*sec(c + d*x))^(3/2)) + (2*(5*A + 7*C)*sin(c + d*x))/(21*b^3*d*sqrt(b*sec(c + d*x))) + (2*A*tan(c + d*x))/(7*d*(b*sec(c + d*x))^(7/2))]

    #= ::Title:: =#
    #=Integrands*of*the*form*(a+b*sec(e+f*x))^m*(A+B*sec(e+f*x)+C*sec(e+f*x)^2)=#
    end

max_stars_repo_path: test/4 Trig functions/4.5 Secant/4.5.4.1 (a+b sec)^m (A+B sec+C sec^2).jl


text:
function prime_field(K::NumField)
  return FlintRationalField()
end

function prime_field(F::FqNmodFiniteField; cached::Bool = true)
  return GF(Int(characteristic(F)), cached = cached)
end

function prime_field(F::FqFiniteField; cached::Bool = true)
  return GF(characteristic(F), cached = cached)
end

function prime_field(F::T; cached::Bool = true) where T <: Union{GaloisField, GaloisFmpzField}
  return F
end

function Nemo.hom(F::FinField, K::FinField, a::FinFieldElem; check::Bool = true)
  @assert parent(a) == K

  if check
    @assert iszero(defining_polynomial(F)(a))
  end

  #We need a preimage function
  p = characteristic(K)
  Kp = prime_field(K, cached = false)
  Kpx = PolynomialRing(Kp, "x", cached = false)[1]
  Fp = prime_field(F, cached = false)
  Fpx = PolynomialRing(Fp, "x", cached = false)[1]
  M = zero_matrix(Kp, degree(F), degree(K))
  el = one(K)
  M[1, 1] = one(Kp)
  for i = 2:nrows(M)
    el = mul!(el, el, a)
    for j = 1:ncols(M)
      M[i, j] = coeff(el, j-1)
    end
  end

  aux = zero_matrix(Kp, 1, degree(F))
  aux1 = zero_matrix(Kp, 1, degree(K))
  function img(x::FinFieldElem)
    @assert parent(x) == F
    for i = 1:degree(F)
      aux[1, i] = Kp(coeff(x, i-1))
    end
    mul!(aux1, aux, M)
    pol = Kpx(elem_type(Kp)[aux1[1, j] for j = 1:degree(K)])
    return K(pol)
  end

  function preimg(x::FinFieldElem)
    @assert parent(x) == K
    for i = 1:degree(K)
      aux1[1, i] = Kp(coeff(x, i-1))
    end
    fl, y = can_solve_with_solution(M, aux1, side = :left)
    if !fl
      error("The element is not in the image!")
    end
    polF = Fpx(elem_type(Fp)[Fp(lift(y[1, j])) for j = 1:degree(F)])
    return F(polF)
  end

  return Nemo.FinFieldMorphism(F, K, img, preimg)
end

#Function that finds an embedding of k into K
function Nemo.embed_any(k::FinField, K::FinField)
  f = defining_polynomial(k)
  rt = roots(f, K)[1]
  return hom(k, K, rt)
end

function tr(a::FinFieldElem, mp::Nemo.FinFieldMorphism)
  @assert parent(a) == codomain(mp)
  q = order(domain(mp))
  d = divexact(degree(codomain(mp)), degree(domain(mp)))
  el = a
  aq = a
  for i = 1:d-1
    aq = aq^q
    el += aq
  end
  return mp\(el)
end

function norm(a::FinFieldElem, mp::Nemo.FinFieldMorphism)
  @assert parent(a) == codomain(mp)
  q = order(domain(mp))
  d = divexact(degree(codomain(mp)), degree(domain(mp)))
  el = a
  aq = a
  for i = 1:d-1
    aq = aq^q
    el *= aq
  end
  return mp\(el)
end

#TODO: compute via linear algebra.
function minpoly(a::FinFieldElem, mp::Nemo.FinFieldMorphism)
  conjs = Set{typeof(a)}()
  q = order(domain(mp))
  push!(conjs, a)
  el = a^q
  while !(el in conjs)
    push!(conjs, el)
    el = el^q
  end
  t = PolynomialRing(codomain(mp), "x", cached = false)[2]
  pol = prod(t-x for x in conjs)
  return map_coeffs(preimage_map(mp), pol)
end

function Nemo.generator_minimal_polynomial(mp::Nemo.FinFieldMorphism{FqFiniteField, FqFiniteField})
  return minpoly(gen(codomain(mp)), mp)
end

function Nemo.any_root(f::PolyElem{T}) where T <: FinFieldElem
  return roots(f)[1]
end


text:
import .CUDA
import Adapt
"""
    struct TextureEx end

Executes the computation using CUDA Texture arrays.
Computations are inaccurate, but extremly fast.
"""
struct TextureEx end

struct AxesProduct{T,N,A} <: AbstractArray{T,N}
    axes::A
    function AxesProduct(axes::NTuple{N,AbstractVector}) where {N}
        T = Tuple{map(eltype, axes)...}
        A = typeof(axes)
        return new{T,N,A}(axes)
    end
end
Base.size(o::AxesProduct) = map(length, o.axes)
Adapt.@adapt_structure AxesProduct

Base.@propagate_inbounds function Base.getindex(o::AxesProduct, I::Integer...)
    map(getindex, o.axes, Tuple(I))
end

Base.IndexStyle(::Type{<:AxesProduct}) = IndexCartesian()

resolve_address_mode(::LI.Replicate) = CUDA.ADDRESS_MODE_CLAMP
resolve_address_mode(::LI.Reflect)   = CUDA.ADDRESS_MODE_MIRROR

@noinline function check_data_to(data_to::NamedTuple)
    for ax in Tuple(data_to.axes)
        @argcheck ax isa AbstractRange
        @argcheck eltype(ax) === Float32
        @argcheck length(ax) >= 2
        @argcheck issorted(ax)
    end
end

function _pullback_namedtuple(f, axes_from, data_to, ::TextureEx; kw...)
    out_vals = CUDA.CuArray{Float32}(undef, length.(Tuple(axes_from)))
    out = (axes=axes_from, values=out_vals)
    _pullback_namedtuple!(f, out, data_to, TextureEx(); kw...)
end

function _pullback_namedtuple!(f, out, data_to, ::TextureEx;
                               extrapolate=LI.Replicate())
    check_data_to(data_to)
    @argcheck extrapolate === LI.Replicate() # others are not handled properly by CUDA.jl
    AC.check_consistency(out)
    address_mode = resolve_address_mode(extrapolate)
    vals_tex = CUDA.CuTexture(CUDA.CuTextureArray(data_to.values);
        address_mode,
        interpolation=CUDA.LinearInterpolation()
    )
    pts = AxesProduct(Tuple(out.axes))
    axes_to = Tuple(data_to.axes)
    broadcast!(out.values, pts, Ref(vals_tex)) do pt, tex
        pt = Tuple(f(pt))
        pt = indexcoords(pt, axes_to)
        tex[pt...]
    end
    out
end

indexcoords(pt, axes) = map(indexcoord, pt, axes)
function indexcoord(x, r)
    (x - first(r)) / step(r) + firstindex(r)
end


text:
using HDF5

# first, plot acquisition benchmark with input-noise-corrected data
path = "SARA/NatCom2020/outer/results/"
f = h5open(path * "outer_acquisition_benchmark.h5")
nrandom = read(f, "nrandom")
r2 = read(f, "r2")
close(f)

f = h5open(path * "outer_acquisition_benchmark_input_noise_0.h5")
r2_no_input_noise = read(f, "r2")
close(f)

using Statistics
function random_statistics(x::AbstractMatrix, nrandom::Int)
    random_statistics(@view(x[:, 1:nrandom]))
end
function random_statistics(x::AbstractMatrix)
    m = vec(mean(x, dims = 2))
    s = vec(std(x, dims = 2))
    return m, s
end

function plot_results(r2::AbstractMatrix, nrandom::Int; linestyle = :solid)
    # decide how many iterations to plot
    nsample = 64
    r2 = @view r2[1:nsample, :]
    r2_random, r2_random_std = random_statistics(r2, nrandom)

    r2 = @view r2[:, nrandom+1:end]
    r2_sigma = r2[:, 1]
    r2_stripe_sigma = r2[:, 2]
    nucb = (size(r2, 2)-3) ÷ 2
    r2_ucb = r2[:, 3:2+nucb]
    r2_stripe_ucb = r2[:, 3+nucb:end]

    # choose the ucb parameter that is best at iteration ind
    ind = nsample ÷ 2
    k = argmax(r2_ucb[ind, :])
    r2_ucb = r2_ucb[:, k]
    k = argmax(r2_stripe_ucb[ind, :])
    r2_stripe_ucb = r2_stripe_ucb[:, k]

    p!(y, label, color) = plot!(y, label = label, linestyle = linestyle, linecolor = color)
    data = [r2_random, r2_sigma, r2_stripe_sigma, r2_ucb, r2_stripe_ucb]
    labels = ["random", "sigma", "stripe sigma", "ucb", "stripe ucb"]
    colors = palette(:seaborn_muted)[1:length(data)]
    p!.(data, labels, colors)
end

using Plots
plotly()
plot(legend = :bottomright, linewidth = 1.,
     xlabel = "number of samples", ylabel = "R²ₛ score",
     framealpha = 0, foreground_color_legend = nothing, background_color_legend = nothing,  grid=true, framestyle = :box)
ylims!((0, .8)) # .8 is approximately the converged r2 value
plot_results(r2, nrandom)
plot_results(r2_no_input_noise, nrandom, linestyle = :dash)
gui()


text:
module EconModels

using Base.Dates, DataFrames

const VERBOSITY    = Dict(:none => 0, :low => 1, :high => 2)
const DATE_FORMAT  = "yymmdd"
const DATASERIES_DELIM = "__"

export

    # interface.jl
    ModelInterface, description, spec, subspec,
    dataroot, saveroot, inpath, workpath, rawpath, tablespath, figurespath, logpath, filestring,

    # settings.jl
    Setting, get_setting, default_settings!, default_test_settings!, data_vintage,

    # data_utils.jl
    prev_quarter, next_quarter, get_quarter_ends, stringstodates, quartertodate, subtract_quarters,
    format_dates!, iterate_quarters, na2nan!, iterate_quarters,

    # constants
    VERBOSITY, DATE_FORMAT, DATASERIES_DELIM

include("interface.jl")
include("settings.jl")
include("data_utils.jl")

end # module


text:
module Zoned

include("create_zonation_funcs.jl")
include("zonation_forward.jl")
include("zonation_objective.jl")

# zonation Crank-Nicolson
export smooth_zon_objective, create_constraint_zon, create_jac_constraint, grad_smooth_zon_objective # main zonation funcs
export create_zon_loglike_objective, create_grad_loglike_objective
export create_jac_constraint_monotonic_cooling, create_constraint_monotonic_cooling
export decrease_constraint, jac_decrease_constraint

# objective
export decreasing_zon_objective, grad_decreasing_zon_objective
export grad_null_objective, null_objective
export create_decreasing_objective, create_grad_decreasing_objective
export zonation_forward # zonation func for testing
export ind_dPb, ind_dPb_vec

end


text:
using Test, LinearAlgebra

using XPORTA: write_ieq, IEQ

using BellScenario

@testset "X-2-4 scenario" begin

@time for X in 5:9
    filepath = "./data/quick_adjacency_decomposition/"
    filename = "$X-2-4"

    println("$filename scenario")

    PM = LocalSignaling(X,4,2)
    BG_seed = BellGame(cat([1 0 0 0;1 0 0 0;0 1 0 0;0 0 1 0],zeros(Int64,(4,X-4)), dims=2),2)

    vertices = LocalPolytope.vertices(PM)
    canonical_facets = LocalPolytope.adjacency_decomposition(vertices, BG_seed, PM, max_vertices=200)
    games = collect(keys(canonical_facets))

    pretty_print_txt(games, filepath*filename)
    write_ieq(filepath*"ieq/"*filename, convert(IEQ, games))
end

end


text:
using ElementarySymmetricFunctions
using Test

function naive_esf(x::AbstractVector{T}) where T <: Real
    n = sum(x .!= 0)
    xidx = x .!= 0
    S = zeros(T, n+1)
    states = hcat(reverse.(digits.(0:2^n-1,base=2,pad=n))...)'

    r_states = vec(mapslices(sum, states, dims=2))

    for r in 0:n
        idx = findall(r_states .== r)
        S[r+1] = sum(mapslices(x->prod(x[x .!= 0]), states[idx, :] .* x[xidx]', dims=2))
    end
    return S
end

function naive_pb(p::AbstractVector{T}) where T <: Real
    x = p ./ (1 .- p)
    naive_esf(x) * prod(1 ./ (1 .+ x))
end

@testset "ElementarySymmetricFunctions.jl" begin
    x = [3.5118, .6219, .2905, .8450, 1.8648]
    n = length(x)

    naive_sol = naive_esf(x)
    naive_sol_reg = naive_sol ./ binomial.(n,0:n)

    #test recursive and FFT algos
    @test esf_sum(x) ≈ naive_sol
    @test esf_sum_reg(x) ≈ naive_sol_reg
    @test esf_sum_log(x) ≈ log.(naive_sol)

    @test esf_dc_fft(x) ≈ naive_sol
    @test esf_dc_fft_reg(x) ≈ naive_sol_reg

    #test derivatives and difference algos
    G = zeros(n,n)
    [G[:,i] .= x for i in 1:n]
    [G[i,i] = 0 for i in 1:n]
    naive_P = mapslices(naive_esf, G, dims=1)'
    naive_P_reg = mapslices(x->naive_esf(x) ./ binomial.(n, 1:n), G, dims=1)'

    S,P = esf_diff(x)
    @test S ≈ naive_sol
    @test P ≈ naive_P

    S,P,err = esf_diff_updown(x)
    @test S ≈ naive_sol
    @test P ≈ naive_P

    S,P = esf_diff_reg(x)
    @test S ≈ naive_sol_reg
    @test P ≈ naive_P_reg

    y = zeros(n)
    H_naive = zeros(n,n,n+1)
    for i in 1:n
        for j in 1:n
            y .= x; y[i]=0; y[j]=0;
            H_naive[i,j,1:(sum(y .!=0)+1)] .= naive_esf(y)
        end
    end

    @test esf_sum_dervs_1(x)[2] ≈ vcat([H_naive[i,i,:] for i in 1:n]'...)
    @test esf_sum_dervs_2(x)[2] ≈ H_naive

    y = zeros(n)
    H_naive_reg = zeros(n,n,n+1)
    for i in 1:n
        for j in 1:n
            y .= x; y[i]=0; y[j]=0;
            H_naive_reg[i,j,1:(sum(y .!=0)+1)] .= naive_esf(y) ./ binomial.(n, (0+sum(y .== 0)):n)
        end
    end
    @test esf_sum_dervs_1_reg(x)[2] ≈ vcat([H_naive_reg[i,i,:] for i in 1:n]'...)
    @test esf_sum_dervs_2_reg(x)[2] ≈ H_naive_reg

    #test poisson-binomial algos
    p = x ./ (1 .+ x)
    naive_sol = naive_pb(p)

    @test poisbin_sum_taub(p) ≈ naive_sol
    @test poisbin_fft(p) ≈ naive_sol
    @test poisbin_fft_cf(p) ≈ naive_sol
    @test poisbin_chen(p) ≈ naive_sol
    @test poisbin_dc_fft(p) ≈ naive_sol
    @test poisbin_sum_taub_log(p) ≈ log.(naive_sol)

end


text:
using Random
using Distributions

Random.seed!(1234)
description = "Hummert six-state model without conserved charge."

#F/RT
ξ = (96.485332)/(298*8.314462)

#index after which all states are open
oi = 3
# indices of activated states
ai = [1, 2, 4, 5]

"""
Notes on how performance can be improved:
* Modify the objective function so that `db_pars` is called only once per evoluation of the objective function.
    i.e. `tm(...)` should not have to call `db_pars`
* Modify `tm(...)` to update a pre-allocated array in-place
* Compute voltage dependent rates `a1, b1, ...` in `tm(...)` in one line by directly multiplying the coefficients `a0, b0, ...`
* Use a static vector `SVector` for parameters, e.g. the output of `db_pars`
* Use views of the transition matrix to define the ODE, i.e. argument `Q` in `rhs!(...)`

Similar performance tips are available in `main_v4.jl`.
The above apply equally to `four_3s_ramentol.jl`.
"""

function db_pars(pars::AbstractVector{}, db; fmt="x")::Vector{}
    """
    Uses OR -> OA (a2; c0, sc) and OA-OH (a3; e0, se) transitions to set detailed balance
        (minimum spanning tree method described by Colquhoun)
    pars = parameters
    DB_CONDITION = uses above method to set detailed balance, unless `x`
    """
    if db == "x"
        a10, b10, s1, a20, b20, s2, a30, b30, s3, a40, b40, s4, g1, g2, h1, h2 = exp.(pars)
    elseif db == "h2"
        a10, b10, s1, a20, b20, s2, a30, b30, s3, a40, b40, s4, g1, g2, h1  = exp.(pars)
        h2 = (g2*b40*b30*h1*a10*a20)/(a40*a30*g1*b10*b20)
    elseif db == "a3"
        a10, b10, s1, a20, b20, s2, b30, a40, b40, s4, g1, g2, h1, h2  = exp.(pars)
        a30 = (a10*a20*g2*b40*b30*h1)/(a40*h2*b20*b10*h1)
        s3 = 1/((1/s1) + (1/s2) - (1/s4))
    else
        println("   `db` not understood. Current options are 'h2' and 'a3'")
    end

    if fmt == "log"
        return log.([a10, b10, s1, a20, b20, s2, a30, b30, s3, a40, b40, s4, g1, g2, h1, h2])
    else
        return [a10, b10, s1, a20, b20, s2, a30, b30, s3, a40, b40, s4, g1, g2, h1, h2]
    end
end

function get_vdep_pars(p::AbstractVector{}, v)
    a1, a2, a3, a4 = [p[i]*exp(-v/p[i+2]) for i = 1:3:10]
    b1, b2, b3, b4 = [p[i]*exp(v/p[i+1]) for i = 2:3:11]

    return [a1, b1, a2, b2, a3, b3, a4, b4]
end

function tm(v, p::AbstractVector{}, db)
    a10, b10, s1, a20, b20, s2, a30, b30, s3, a40, b40, s4, g1, g2, h1, h2 = db_pars(p, db)

    a1 = a10*exp(-v/s1)
    b1 = b10*exp(v/s1)
    a3 = a30*exp(-v/s2)
    b3 = b30*exp(v/s2)

    a2 = a20*exp(-v/s3)
    b2 = b20*exp(v/s3)
    a4 = a40*exp(-v/s4)
    b4 = b40*exp(v/s4)

    Q = [
        a1 -(b1+a2+a1) (b2-a1) -a1 -a1 -a1;
        0 a2 -(b2+g2) 0 0 h2;
        g1 -g1 -g1 -(h1+a3+g1) (b3-g1) -g1;
        0 0 0 a3 -(b3+a4) b4;
        0 0 g2 0 a4 -(h2+b4);
    ]
    return Q
end

function rhs!(du, u, Q, t)
    du .= Q[:,2:end] * u + Q[:,1]
end
function prepare_BoundsAndParams(pars::Array{Float64,1}, bounds, db)
    # remove entries from parameter and bound arrays given detailed balance
    if db == true || db == "x"
        return pars, bounds
    else
        if DB_CONDITION == "h2"
            deleteat!(pars, 16)
            deleteat!(bounds, 16)
        elseif DB_CONDITION == "a3"
            deleteat!(pars, [7, 9])
            deleteat!(bounds, [7, 9])
        end
    end
    return pars, bounds
end

bds = [(1e-10, 1e3), (1e-10, 1e3), (6, 1e3), #a1/b1
        (1e-10, 1e3), (1e-10, 1e3), (6, 1e3),
        (1e-10, 1e3), (1e-10, 1e3), (6, 1e3), #a3/b3
        (1e-10, 1e3), (1e-10, 1e3), (6, 1e3),
        (1e-5, 10),
        (1e-8, 10), #g1
        (1e-4, 10),
        (1e-8, 1e-1), #h1
]

ln_par_bds = [log.(x) for x in bds]

p0 = [-1.366795614121635, -4.731276419445293, 4.146181643448856,
        -3.0831331501017187, 2.3286672592590754, 2.9630705986068224,
        -3.698762226164891, 1.5730159834565232, 4.15,
        1.998606200149607, 0.8069046570541617, 2.96,
        -2.097864956067902, -6.917709307586828, -5.092722946312469, -4.884952245097995]

max_stars_repo_name: haganenoneko/FunnyChannelGatingReview

max_issues_repo_name: haganenoneko/FunnyChannelGatingReview

max_forks_repo_name: haganenoneko/FunnyChannelGatingReview


text:
using OnionSockets

using CryptoGroups
using DiffieHellman
using Random
using Sockets
import Multiplexers: Multiplexer, route, forward

G = CryptoGroups.MODP160Group()
chash(envelopeA,envelopeB,key) = hash("$envelopeA $envelopeB $key")
dh = DH(data->(data,nothing),envelope->envelope,G,chash,() -> rand(1:BigInt(1)<<100))

ballotmember = dh
memberballot = dh
ballotgate = dh
gateballot = dh
membergate = dh
gatemember = dh

N = 2

@sync begin
    @async let
        routers = listen(2001)
        try
            @show "Router"
            serversocket = accept(routers)
            secureserversocket = accept(serversocket,nothing,ballotgate)

            mux = Multiplexer(secureserversocket,N) # Perhaps I could put route process with in mux.
            task = @async route(mux)

            susersockets = []
            for i in 1:N
                securesocket = accept(mux.lines[i],ballotmember)
                push!(susersockets,securesocket)
            end

            for i in 1:N
                serialize(susersockets[i],"A secure message from the router")
                @show deserialize(susersockets[i])
            end

            serialize(secureserversocket,:Terminate)
            wait(task)
        finally
            close(routers)
        end
    end

    @async let
        server = listen(2000)
        try
            @show "Server"

            routersocket = connect(2001)
            secureballotbox = connect(routersocket,nothing,gateballot)

            usersockets = IO[]

            while length(usersockets)<N
                secureusersocket = accept(server,nothing,gatemember)
                push!(usersockets,secureusersocket)
            end

            forward(usersockets,secureballotbox)
        finally
            close(server)
        end
    end

    @async let
        @show "User 1"
        usersocket = connect(2000)
        securesocket = connect(usersocket,nothing,membergate)
        sroutersocket = connect(securesocket,nothing,memberballot)

        @show deserialize(sroutersocket)
        serialize(sroutersocket,"A scuere msg from user 1")
    end

    @async let
        @show "User 2"
        usersocket = connect(2000)
        securesocket = connect(usersocket,nothing,membergate)
        sroutersocket = connect(securesocket,nothing,memberballot)

        @show deserialize(sroutersocket)
        serialize(sroutersocket,"A scuere msg from user 2")
    end
end


text:
function dist(a,b)
    ret = 0
    for i = 1:length(a)
        c, v = a[i], b[i]
        ret += c!=v ? 1 : 0
    end
    ret
end

a=readline()
b=readline()

sol=999999999
for i = 1:length(a)-length(b)+1
    global sol
    sa = a[i:i+length(b)-1]
    sol=min(dist(sa,b), sol)
end
println(sol)

max_stars_repo_name: sugitanishi/competitive-programming

max_issues_repo_name: sugitanishi/competitive-programming

max_forks_repo_name: sugitanishi/competitive-programming


text:
###############################################################################
#
# Cache of values to speed the canopy_fluxes!
#
###############################################################################
"""
    mutable struct CFCache{FT}

Cache to speed [`canopy_fluxes!`](@ref) by pre-allocating arrays

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct CFCache{FT}
    "absorbed energy from wave lengths"
    abs_wave   ::Array{FT,1}
    "absfs' * lidf [nAzi]"
    absfs_lidf ::Array{FT,1}
    "wave length energy [same as dWL]"
    E_all      ::Array{FT,1}
    "wave length energy [same as iPAR]"
    E_iPAR     ::Array{FT,1}
    "lPs [nLayer]"
    lPs        ::Array{FT,1}
    "kChlrel [same as iPAR]"
    kChlrel    ::Array{FT,1}
    "diffusive PAR [same as iPAR]"
    PAR_diff   ::Array{FT,1}
    "diffusive PAR for photosynthesis [same as iPAR]"
    PAR_diffCab::Array{FT,1}
    "direct PAR [same as iPAR]"
    PAR_dir    ::Array{FT,1}
    "diffusive PAR for photosynthesis [same as iPAR]"
    PAR_dirCab ::Array{FT,1}
end

###############################################################################
#
# Cache of values to speed the canopy_geometry!
#
###############################################################################
"""
    mutable struct CGCache{FT}

Cache to speed [`canopy_geometry!`](@ref) by pre-allocating arrays

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct CGCache{FT}
    # 1D arrays
    "cos_ttli .* cos(vza) dim: nIncl"
    _Co ::Array{FT,1}
    "cos_ttli .* cos(sza) dim: nIncl"
    _Cs ::Array{FT,1}
    "sin_ttli .* sin(vza) dim: nIncl"
    _So ::Array{FT,1}
    "sin_ttli .* sin(sza) dim: nIncl"
    _Ss ::Array{FT,1}

    # 2D arrays
    "maxtrix filled with 1 dim: (1, nAzi)"
    _1s ::Array{FT,2}
    "2D array to speed up _cds and _cdo dim: (nIncl, nAzi)"
    _2d ::Array{FT,2}
    "_Co * _1s .+ _So * cos_philo' dim: (nIncl, nAzi)"
    _cdo::Array{FT,2}
    "_Cs * _1s .+ _Ss * cos_ttlo' dim: (nIncl, nAzi)"
    _cds::Array{FT,2}
end

###############################################################################
#
# Cache of values to speed the short_wave!
#
###############################################################################
"""
    mutable struct SFCache{FT}

Cache to speed [`SIF_fluxes!`](@ref) by pre-allocating arrays

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct SFCache{FT}
    # 1D array
    M⁻_sun        ::Array{FT,1}
    M⁺_sun        ::Array{FT,1}
    wfEs          ::Array{FT,1}
    sfEs          ::Array{FT,1}
    sbEs          ::Array{FT,1}
    M⁺⁻           ::Array{FT,1}
    M⁺⁺           ::Array{FT,1}
    M⁻⁺           ::Array{FT,1}
    M⁻⁻           ::Array{FT,1}
    sun_dwl_iWlE  ::Array{FT,1}
    tmp_dwl_iWlE  ::Array{FT,1}
    ϕ_cosΘ_lidf   ::Array{FT,1}
    vfEplu_shade  ::Array{FT,1}
    vbEmin_shade  ::Array{FT,1}
    vfEplu_sun    ::Array{FT,1}
    vbEmin_sun    ::Array{FT,1}
    sigfEmin_shade::Array{FT,1}
    sigbEmin_shade::Array{FT,1}
    sigfEmin_sun  ::Array{FT,1}
    sigbEmin_sun  ::Array{FT,1}
    sigfEplu_shade::Array{FT,1}
    sigbEplu_shade::Array{FT,1}
    sigfEplu_sun  ::Array{FT,1}
    sigbEplu_sun  ::Array{FT,1}
    zeroB         ::Array{FT,1}
    tmp_1d_nWlF   ::Array{FT,1}
    tmp_1d_nLayer ::Array{FT,1}
    dnorm         ::Array{FT,1}

    # 2D array
    "transmission of diffusive light?"
    τ_dd                ::Array{FT,2}
    "extinction of diffuse light?"
    ρ_dd                ::Array{FT,2}
    Xdd                 ::Array{FT,2}
    Rdd                 ::Array{FT,2}
    Y                   ::Array{FT,2}
    U                   ::Array{FT,2}
    S⁻                  ::Array{FT,2}
    S⁺                  ::Array{FT,2}
    piLs                ::Array{FT,2}
    piLd                ::Array{FT,2}
    Fsmin               ::Array{FT,2}
    Fsplu               ::Array{FT,2}
    Fdmin               ::Array{FT,2}
    Fdplu               ::Array{FT,2}
    Femo                ::Array{FT,2}
    M⁺                  ::Array{FT,2}
    M⁻                  ::Array{FT,2}
    ϕ_cosΘ              ::Array{FT,2}
    F⁻                  ::Array{FT,2}
    F⁺                  ::Array{FT,2}
    net_diffuse         ::Array{FT,2}
    tmp_2d_nWlF_nLayer  ::Array{FT,2}
    tmp_2d_nWlF_nLayer_2::Array{FT,2}
end

###############################################################################
#
# Cache of values to speed the short_wave!
#
###############################################################################
"""
    mutable struct SWCache{FT}

Cache to speed [`short_wave!`](@ref) by pre-allocating arrays

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct SWCache{FT}
    # 1D arrays
    "dnorm?"
    dnorm ::Array{FT,1}
    "pi * Lo"
    piLo  ::Array{FT,1}
    "pi * Lo from canopy"
    piLoc ::Array{FT,1}
    "pi * Lo from soil"
    piLos ::Array{FT,1}

    # 2D arrays
    "pi * Lo from canopy 2D matrix"
    piLoc2::Array{FT,2}
    "extinction of diffuse light?"
    ρ_dd  ::Array{FT,2}
    "extinction of direct light?"
    ρ_sd  ::Array{FT,2}
    "transmission of diffusive light?"
    τ_dd  ::Array{FT,2}
    "transmission of direct light?"
    τ_sd  ::Array{FT,2}
end

###############################################################################
#
# Cache of values to speed up the calculations
#
###############################################################################
"""
    mutable struct RTCache{FT}

Collection of caches to speed up RT module

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct RTCache{FT}
    "[`CFCache`](@ref) type cache"
    cf_con::CFCache{FT}
    "[`CGCache`](@ref) type cache"
    cg_con::CGCache{FT}
    "[`SFCache`](@ref) type cache"
    sf_con::SFCache{FT}
    "[`SWCache`](@ref) type cache"
    sw_con::SWCache{FT}
end

max_stars_repo_name: climate-machine/Land


text:
function runcode(app)

display("run code...");
mpirun = app.mpirun;
if app.platform == "cpu"
    if app.mpiprocs==1
        str = "./app/serial" * app.appname * " datain/ dataout/out";
    else
        str = mpirun * " -np " * string(app.mpiprocs) * " ./app/mpi" * app.appname * " datain/ dataout/out";
    end
    run(string2cmd(str));
elseif app.platform == "gpu"
    if app.mpiprocs==1
        str = "./app/gpu" * app.appname * " datain/ dataout/out";
    else
        str = mpirun * " -np " * string(app.mpiprocs) * " ./app/gpumpi" * app.appname * " datain/ dataout/out";
    end
    run(string2cmd(str));
else
    error("app.platform must be either cpu or gpu");
end

return str

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
# Copyright 2017-2019 Eric S. Tellez
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import StatsBase: predict, fit
export NearNeighborRegression, optimize!, predict, predict_proba
using Statistics
using SimilaritySearch:
    Sequential, KnnResult, empty!

mutable struct NearNeighborRegression{IndexType,DataType}
    dist::Function
    X::IndexType
    y::Vector{DataType}
    k::Int
    summarize::Function

    function NearNeighborRegression(dist::Function, X::AbstractVector{ItemType}, y::AbstractVector{DataType}; summarize=mean, k::Int=1) where {ItemType, DataType}
        index = fit(Sequential, X)
        new{typeof(index), DataType}(dist, index, y, k, summarize)
    end
end

function predict(nnc::NearNeighborRegression{IndexType,DataType}, vector) where {IndexType,DataType}
    [predict_one(nnc, item) for item in vector]
end

function predict_one(nnc::NearNeighborRegression{IndexType,DataType}, item) where {IndexType,DataType}
    res = KnnResult(nnc.k)
    search(nnc.X, nnc.dist, item, res)
    DataType[nnc.y[p.objID] for p in res] |> nnc.summarize
end

function _train_create_table_reg(dist::Function, train_X, train_y, test_X, k::Int)
    index = fit(Sequential, train_X)
    res = KnnResult(k)
    function f(x)
        empty!(res)  # this is thread unsafe
        search(index, dist, x, res)
        [train_y[p.objID] for p in res]
    end

    f.(test_X)
end

function _train_predict(nnc::NearNeighborRegression{IndexType,DataType}, table, test_X, k) where {IndexType,DataType}
    A = Vector{DataType}(undef, length(test_X))
    for i in 1:length(test_X)
        row = table[i]
        A[i] = nnc.summarize(row[1:k])
    end

    A
end

function gmean(X)
    prod(X)^(1/length(X))
end

function hmean(X)
    d = 0.0
    for x in X
        d += 1.0 / x
    end
    length(X) / d
end

function optimize!(nnr::NearNeighborRegression, scorefun::Function; summarize_list=[mean, median, gmean, hmean], runs=3, trainratio=0.5, testratio=0.5, folds=0, shufflefolds=true)
    mem = Dict{Tuple,Float64}()

    function f(train_X, train_y, test_X, test_y)
        _nnr = NearNeighborRegression(nnr.dist, train_X, train_y)
        kmax = sqrt(length(train_y)) |> round |> Int
        table = _train_create_table_reg(nnr.dist, train_X, train_y, test_X, kmax)
        k = 2
        while k <= kmax
            _nnr.k = k - 1
            for summarize in summarize_list
                _nnr.summarize = summarize
                pred_y = _train_predict(_nnr, table, test_X, _nnr.k)
                score = scorefun(test_y, pred_y)
                key = (k - 1, summarize)
                mem[key] = get(mem, key, 0.0) + score
            end
            k += k
        end
        0
    end

    if folds > 1
        kfolds(f, nnr.X.db, nnr.y, folds=folds, shuffle=shufflefolds)
        bestlist = [(score/folds, conf) for (conf, score) in mem]
    else
        montecarlo(f, nnr.X.db, nnr.y, runs=runs, trainratio=trainratio, testratio=testratio)
        bestlist = [(score/runs, conf) for (conf, score) in mem]
    end

    sort!(bestlist, by=x -> (-x[1], x[2][1]))
    best = bestlist[1]
    nnr.k = best[2][1]
    nnr.summarize = best[2][2]

    bestlist
end


text:
@testset "Buffer" begin
  w = rand(2,3)
  b = rand(2)
  c = rand(2)

  buff = Zygote.Buffer([])
  @test length(buff) == 0

  push!(buff, w)
  push!(buff, b)
  push!(buff, c)
  @test length(buff) == 3
  @test buff[1] === w
  @test buff[2] === b
  @test buff[3] === c

  deleteat!(buff, 2)
  @test length(buff) == 2
  @test buff[1] === w
  @test buff[2] === c

  deleteat!(buff, 1)
  @test length(buff) == 1
  @test buff[1] === c
end

@testset "isderiving" begin

  function f(x)
    if Zygote.isderiving(x)
      x^2
    else
      2x^2
    end
  end

  # Test higher order derivatives
  gs = gradient(4) do x
    gradient(x) do y
      f(y)
    end[1]
  end

  @test gs == (2,)

  struct Tester
    cpu_offload::Float64
  end

  function Tester(p)
      @show Zygote.isderiving(p)
      cpu_offload = Zygote.isderiving(p) ? 0.0 : 0.2
      Tester(cpu_offload)
  end

  function f(p)
    sum(Tester(p).cpu_offload .* p)
  end

  p = [1.0]
  gs = gradient(f, p)
  @test gs[1] == [0.]

end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
abstract type AbstractProblem end

struct CellProblemAdvecTemp{MeshType,TbfType,TypeK,TypeS,VType,BcondType,LaplacianStruct,AdvectionStruct} <: AbstractProblem
    Tc::Vector{Float64}
    k::TypeK
    s::TypeS
    ρC::Float64
    u::VType
    bcond::BcondType
    mesh::MeshType
    laplacian!::LaplacianStruct
    advection!::AdvectionStruct
    ∇Tc::Vector{Vec2D{Float64}}
    Tbf::TbfType
end

function CellProblemAdvecTemp(Tc,mesh,d)
    bcond = boundary_conditions(d)

    k = ConstVec{d[:conductivity]}()
    s = ConstVec{d[:source]}()
    ρC = d[Symbol("rho*C")]
    uvec = ConstVec{d[:velocity]}()
    ∇Tc = zeros(Vec2D{Float64},length(Tc))
    Tbf = FieldAtBoundary(Tc,mesh,bcond)

    advection = get_advection_method(d)
    laplacian = get_laplacian_method(Tc,bcond,k,mesh,d)
    types = typeof.((mesh,Tbf,k,s,uvec,bcond,laplacian,advection))
    return CellProblemAdvecTemp{types...}(Tc,k,s,ρC,uvec,bcond,mesh,laplacian,advection,∇Tc,Tbf)
end

@inline needs_gradient_calculation(p) = (typeof(p.advection!) === UpWind2ndOrder) && (typeof(p.laplacian!) === MethodB) && (typeof(p.u) <: ConstVec && p.u[1] == zero(eltype(p.u)))

function calculate_rhs!(rhs,p::CellProblemAdvecTemp)
    fill!(rhs,zero(eltype(rhs)))

    is_implicit(p.laplacian!) || p.laplacian!(rhs,p)

    needs_gradient_calculation(p) && gradient!(p.∇Tc,FaceSimpleInterpolation(p.Tc,p.mesh),p.Tbf,p.mesh.f2cloops)

    if !(typeof(p.u) <: ConstVec && p.u[1] == zero(eltype(p.u)))
        p.advection!(rhs,p)
    end
    if !(typeof(p.s) <: ConstVec && p.s[1] == zero(eltype(p.s)))
        add_source!(rhs,p)
    end

    is_implicit(p.laplacian!) && p.laplacian!(rhs,p)

    return nothing
end

function add_source!(rhs,p)
    s = p.s
    @inbounds @simd for i=1:length(rhs)
        rhs[i] += s[i]
    end
end

struct StokesProblem{UfType,UbfType,MeshType,uBcondType,LaplacianStruct,Ma,PMa,PMap} <: AbstractProblem
    u::Vec2DArray{Float64}
    u_old::Vec2DArray{Float64}
    uf::UfType
    ubf::UbfType
    δu::Vec2DArray{Float64}
    ru::Vec2DArray{Float64}
    p::Vector{Float64}
    δp::Vector{Float64}
    ν::ConstVec{Float64}
    ρ::ConstVec{Float64}
    bcond::uBcondType
    mesh::MeshType
    laplacian!::LaplacianStruct
    ∇u::Vector{FiniteVolumeMesh.Ten2D{Float64}}
    A::Ma
    pcgA::PMa
    s::Vector{Float64}
    Ap::PoissonP{MeshType}
    pcgAp::PMap
    dt::Float64
end

function StokesProblem(u,mesh,d)
    u_old = similar(u)
    uf = FaceSimpleInterpolation(u,mesh)
    bcond = uboundary_conditions(d)
    ubf = FieldAtBoundary(u,mesh,bcond)
    δu = similar(u)
    ru = similar(u)
    ν = ConstVec(d[:viscosity])
    ρ = ConstVec(d[:density])
    ∇u = zeros(Ten2D{Float64},length(u))
    p = zeros(length(u))
    δp = zeros(length(u))

    laplacian = get_laplacian_method(u,bcond,ν,mesh,d)
    A = aIpDbG(d,mesh)
    pcgA = PCG(δu)
    pcgAp = PCG(δp)
    types = typeof.((uf,ubf,mesh,bcond,laplacian,A,pcgA,pcgAp))

    s = zeros(length(mesh.cells))
    Ap = PoissonP(d,mesh)
    dt = d[:dt]
    return StokesProblem{types...}(u,u_old,uf,ubf,δu,ru,p,δp,ν,ρ,bcond,mesh,laplacian,∇u,A,pcgA,s,Ap,pcgAp,dt)
end

struct NSProblem{UfType,UbfType,MeshType,uBcondType,LaplacianStruct,AdvecStruct,Ma,PMa,PMap,TBcondType,TLaplacianStruct,TAdvectionStruct,TbfType,MaT,PMaT} <: AbstractProblem
    u::Vec2DArray{Float64}
    uf::UfType
    ubf::UbfType
    δu::Vec2DArray{Float64}
    ru::Vec2DArray{Float64}
    p::Vector{Float64}
    δp::Vector{Float64}
    ν::ConstVec{Float64}
    ρ::ConstVec{Float64}
    bcond::uBcondType
    mesh::MeshType
    laplacian!::LaplacianStruct
    uadvection!::AdvecStruct
    ∇u::Vector{FiniteVolumeMesh.Ten2D{Float64}}
    A::Ma
    pcgA::PMa
    sp::Vector{Float64}
    Ap::PoissonP{MeshType}
    pcgAp::PMap
    dt::Float64
    Tc::Vector{Float64}
    ∇Tc::Vec2DArray{Float64}
    rT::Vector{Float64}
    δT::Vector{Float64}
    k::ConstVec{Float64}
    s::ConstVec{Float64}
    ρC::Float64
    Tbcond::TBcondType
    Tlaplacian!::TLaplacianStruct
    Tadvection!::TAdvectionStruct
    Tbf::TbfType
    AT::MaT
    pcgAT::PMaT
end

function NSProblem(u,mesh,d)
    uf = FaceSimpleInterpolation(u,mesh)
    bcond = uboundary_conditions(d)
    ubf = FieldAtBoundary(u,mesh,bcond)
    δu = Vec2DArray{Float64}(length(u))
    ru = Vec2DArray{Float64}(length(u))
    ν = ConstVec(d[:viscosity])
    ρ = ConstVec(d[:density])
    ∇u = zeros(Ten2D{Float64},length(u))
    p = zeros(length(u))
    δp = zeros(length(u))

    laplacian = get_laplacian_method(u,bcond,ConstVec{Float64}(0.5*ν[1]),mesh,d)
    A = aIpDbG(ConstVec(d[:density]/d[:dt]),-0.5*ν[1],d,mesh)
    pcgA = PCG(δu,KrylovCtrl("vPCG.set"))
    pcgAp = PCG(δp,KrylovCtrl("pPCG.set"))
    uadvection = UpWindAdvection(u,ubf,u,ubf,mesh)

    sp = zeros(length(mesh.cells))
    Ap = PoissonP(d,mesh)
    dt = d[:dt]

    Tc = zeros(Float64,length(mesh.cells))
    ∇Tc = Vec2DArray{Float64}(length(u))
    rT = zeros(Float64,length(mesh.cells))
    dT = zeros(Float64,length(mesh.cells))
    Tbcond = boundary_conditions(d)

    k = ConstVec{Float64}(d[:conductivity])
    s = ConstVec{Float64}(d[:source])
    ρC = d[Symbol("rho*C")]
    Tbf = FieldAtBoundary(Tc,mesh,Tbcond)

    Tadvection = UpWindAdvection(Tc,Tbf,u,ubf,mesh)
    Tlaplacian = get_laplacian_method(Tc,Tbcond,ConstVec{Float64}(0.5*k[1]),mesh,d)
    AT = aIpDbG(ConstVec{Float64}(ρC/dt),-0.5*k[1],d,mesh)
    pcgAT = PCG(dT,KrylovCtrl("tPCG.set"))

    types = typeof.((uf,ubf,mesh,bcond,laplacian,uadvection,A,pcgA,pcgAp,Tbcond,Tlaplacian,Tadvection,Tbf,AT,pcgAT))
    return NSProblem{types...}(u,uf,ubf,δu,ru,p,δp,ν,ρ,bcond,mesh,laplacian,uadvection,∇u,A,pcgA,sp,Ap,pcgAp,dt,Tc,∇Tc,rT,dT,k,s,ρC,Tbcond,Tlaplacian,Tadvection,Tbf,AT,pcgAT)
end


text:
#   Debug.Meta:
# ===============
# Metaprogramming tools used throughout the Debug package

module Meta
using Compat, Debug.AST
export Ex, quot, is_expr
export isblocknode, is_function, is_scope_node, is_in_type, introduces_scope
export untyped_comprehensions, typed_comprehensions, comprehensions
export headof, argsof, argof, nargsof

typealias Ex @compat(Union{Expr, ExNode})

quot(ex) = QuoteNode(ex)

is_expr(ex::Ex, head)          = headof(ex) === head
is_expr(ex::Ex, heads::Set)    = headof(ex) in heads
is_expr(ex::Ex, heads::Vector) = headof(ex) in heads
is_expr(ex,     head)          = false
is_expr(ex,     head, n::Int)  = is_expr(ex, head) && nargsof(ex) == n

isblocknode(node) = is_expr(node, :block)

# for both kinds of AST:s
is_function(node)     = false
is_function(node::Ex) = is_expr(node, [:function, :->], 2) ||
    (is_expr(node, :(=), 2) && is_expr(argof(node,1), :call))

const untyped_comprehensions = [:comprehension, :dict_comprehension]
const typed_comprehensions   = [:typed_comprehension,
                                :typed_dict_comprehension]
const comprehensions = [untyped_comprehensions; typed_comprehensions]

const scope_heads = Set([:while; :try; :for; :let; comprehensions...])
is_scope_node(ex) = is_expr(ex, scope_heads) || is_function(ex)

# only for Node/Nothing
if VERSION >= v"0.4.0-dev"
    is_in_type(::Void) = false
else
    is_in_type(::Nothing) = false
end
function is_in_type(node::Node)
    if isa(node.state, Rhs)
        isa(envof(node), LocalEnv) && is_expr(envof(node).source, :type)
    else
        is_in_type(parentof(node))
    end
end

# only for Node
introduces_scope(node::Node) = node.introduces_scope

## Accessors that work on both Expr:s and ExNode:s ##

headof(ex::Expr)   = ex.head
headof(ex::ExNode) = valueof(ex).head

argsof(ex::Expr)   = ex.args
argsof(ex::ExNode) = valueof(ex).args

nargsof(ex)  = length(argsof(ex))
argof(ex, k) = argsof(ex)[k]

end # module