text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    NormSpectralBridge{T}

The `NormSpectralCone` is representable with a PSD constraint, since
``t \\ge \\sigma_1(X)`` if and only if ``[tI X^\\top; X tI] \\succ 0``.
"""
struct NormSpectralBridge{T,F,G} <: AbstractBridge
    row_dim::Int # row dimension of X
    column_dim::Int # column dimension of X
    psd_index::MOI.ConstraintIndex{F,MOI.PositiveSemidefiniteConeTriangle}
end

const NormSpectral{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{NormSpectralBridge{T},OT}

function bridge_constraint(
    ::Type{NormSpectralBridge{T,F,G}},
    model::MOI.ModelLike,
    f::G,
    s::MOI.NormSpectralCone,
) where {T,F,G}
    f_scalars = MOI.Utilities.eachscalar(f)
    t = f_scalars[1]
    row_dim = s.row_dim
    column_dim = s.column_dim
    side_dim = row_dim + column_dim
    psd_set = MOI.PositiveSemidefiniteConeTriangle(side_dim)
    psd_func =
        MOI.Utilities.zero_with_output_dimension(F, MOI.dimension(psd_set))
    for i in 1:side_dim
        MOI.Utilities.operate_output_index!(
            +,
            T,
            MOI.Utilities.trimap(i, i),
            psd_func,
            t,
        )
    end
    X_idx = 2
    for j in 1:column_dim, i in (column_dim+1):side_dim
        MOI.Utilities.operate_output_index!(
            +,
            T,
            MOI.Utilities.trimap(i, j),
            psd_func,
            f_scalars[X_idx],
        )
        X_idx += 1
    end
    psd_index = MOI.add_constraint(model, psd_func, psd_set)
    return NormSpectralBridge{T,F,G}(row_dim, column_dim, psd_index)
end

function MOI.supports_constraint(
    ::Type{NormSpectralBridge{T}},
    ::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.NormSpectralCone},
) where {T}
    return true
end

function MOI.Bridges.added_constrained_variable_types(
    ::Type{<:NormSpectralBridge},
)
    return Tuple{Type}[]
end

function MOI.Bridges.added_constraint_types(
    ::Type{NormSpectralBridge{T,F,G}},
) where {T,F,G}
    return Tuple{Type,Type}[(F, MOI.PositiveSemidefiniteConeTriangle)]
end

function concrete_bridge_type(
    ::Type{<:NormSpectralBridge{T}},
    G::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.NormSpectralCone},
) where {T}
    F = MOI.Utilities.promote_operation(
        vcat,
        T,
        MOI.Utilities.scalar_type(G),
        T,
    )
    return NormSpectralBridge{T,F,G}
end

function MOI.get(
    ::NormSpectralBridge{T,F,G},
    ::MOI.NumberOfConstraints{F,MOI.PositiveSemidefiniteConeTriangle},
)::Int64 where {T,F,G}
    return 1
end

function MOI.get(
    bridge::NormSpectralBridge{T,F,G},
    ::MOI.ListOfConstraintIndices{F,MOI.PositiveSemidefiniteConeTriangle},
) where {T,F,G}
    return [bridge.psd_index]
end

function MOI.delete(model::MOI.ModelLike, bridge::NormSpectralBridge)
    MOI.delete(model, bridge.psd_index)
    return
end

function MOI.get(
    model::MOI.ModelLike,
    ::MOI.ConstraintFunction,
    bridge::NormSpectralBridge{T,F,G},
) where {T,F,G}
    psd_func = MOI.Utilities.eachscalar(
        MOI.get(model, MOI.ConstraintFunction(), bridge.psd_index),
    )
    t = psd_func[1]
    side_dim = bridge.row_dim + bridge.column_dim
    X = psd_func[[
        MOI.Utilities.trimap(i, j) for j in 1:bridge.column_dim for
        i in (bridge.column_dim+1):side_dim
    ]]
    return MOI.Utilities.convert_approx(G, MOI.Utilities.operate(vcat, T, t, X))
end

function MOI.get(
    ::MOI.ModelLike,
    ::MOI.ConstraintSet,
    bridge::NormSpectralBridge,
)
    return MOI.NormSpectralCone(bridge.row_dim, bridge.column_dim)
end

function MOI.supports(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimalStart,
    ::Type{NormSpectralBridge{T,F,G}},
) where {T,F,G}
    return MOI.supports(
        model,
        attr,
        MOI.ConstraintIndex{F,MOI.PositiveSemidefiniteConeTriangle},
    )
end

function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimal,MOI.ConstraintPrimalStart},
    bridge::NormSpectralBridge,
)
    primal = MOI.get(model, attr, bridge.psd_index)
    t = primal[1]
    side_dim = bridge.row_dim + bridge.column_dim
    X = primal[[
        MOI.Utilities.trimap(i, j) for j in 1:bridge.column_dim for
        i in (bridge.column_dim+1):side_dim
    ]]
    return vcat(t, X)
end

function MOI.set(
    model::MOI.ModelLike,
    ::MOI.ConstraintPrimalStart,
    bridge::NormSpectralBridge{T},
    value,
) where {T}
    column_dim = bridge.column_dim
    side_dim = bridge.row_dim + column_dim
    primal = zeros(T, div(side_dim * (side_dim + 1), 2))
    for i in 1:side_dim
        primal[MOI.Utilities.trimap(i, i)] = value[1]
    end
    X_idx = 2
    for j in 1:column_dim, i in (column_dim+1):side_dim
        primal[MOI.Utilities.trimap(i, j)] = value[X_idx]
        X_idx += 1
    end
    MOI.set(model, MOI.ConstraintPrimalStart(), bridge.psd_index, primal)
    return
end

# Given [U X'; X V] is dual on PSD constraint, the dual on NormSpectralCone
# constraint is (tr(U) + tr(V), 2X) in NormNuclearCone.
function MOI.get(
    model::MOI.ModelLike,
    ::MOI.ConstraintDual,
    bridge::NormSpectralBridge,
)
    dual = MOI.get(model, MOI.ConstraintDual(), bridge.psd_index)
    column_dim = bridge.column_dim
    side_dim = bridge.row_dim + column_dim
    t = sum(dual[MOI.Utilities.trimap(i, i)] for i in 1:side_dim)
    X =
        2 * dual[[
            MOI.Utilities.trimap(i, j) for j in 1:column_dim for
            i in (column_dim+1):side_dim
        ]]
    return vcat(t, X)
end

"""
    NormNuclearBridge{T}

The `NormNuclearCone` is representable with an SDP constraint and extra
variables, since ``t \\ge \\sum_i \\sigma_i (X)`` if and only if there exists
symmetric matrices ``U, V`` such that ``[U X^\\top; X V] \\succ 0`` and
``t \\ge (tr(U) + tr(V)) / 2``.
"""
struct NormNuclearBridge{T,F,G,H} <: AbstractBridge
    row_dim::Int # row dimension of X
    column_dim::Int # column dimension of X
    U::Vector{MOI.VariableIndex}
    V::Vector{MOI.VariableIndex}
    ge_index::MOI.ConstraintIndex{F,MOI.GreaterThan{T}}
    psd_index::MOI.ConstraintIndex{G,MOI.PositiveSemidefiniteConeTriangle}
end

const NormNuclear{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{NormNuclearBridge{T},OT}

function bridge_constraint(
    ::Type{NormNuclearBridge{T,F,G,H}},
    model::MOI.ModelLike,
    f::H,
    s::MOI.NormNuclearCone,
) where {T,F,G,H}
    f_scalars = MOI.Utilities.eachscalar(f)
    row_dim = s.row_dim
    column_dim = s.column_dim
    side_dim = row_dim + column_dim
    U_dim = div(column_dim * (column_dim + 1), 2)
    V_dim = div(row_dim * (row_dim + 1), 2)
    U = MOI.add_variables(model, U_dim)
    V = MOI.add_variables(model, V_dim)
    diag_vars = vcat(
        [U[MOI.Utilities.trimap(i, i)] for i in 1:column_dim],
        [V[MOI.Utilities.trimap(i, i)] for i in 1:row_dim],
    )
    rhs = MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.(one(T), diag_vars),
        zero(T),
    )
    ge_index = MOI.Utilities.normalize_and_add_constraint(
        model,
        MOI.Utilities.operate(
            -,
            T,
            f_scalars[1],
            MOI.Utilities.operate!(/, T, rhs, T(2)),
        ),
        MOI.GreaterThan(zero(T)),
        allow_modify_function = true,
    )
    psd_set = MOI.PositiveSemidefiniteConeTriangle(side_dim)
    psd_func = MOI.VectorOfVariables(U)
    nuc_dim = 1 + row_dim * column_dim
    for i in 1:row_dim
        row_i = (1+i):row_dim:nuc_dim
        psd_func = MOI.Utilities.operate(vcat, T, psd_func, f_scalars[row_i])
        psd_func = MOI.Utilities.operate(
            vcat,
            T,
            psd_func,
            MOI.VectorOfVariables(V[[MOI.Utilities.trimap(i, j) for j in 1:i]]),
        )
    end
    psd_index = MOI.add_constraint(model, psd_func, psd_set)
    return NormNuclearBridge{T,F,G,H}(
        row_dim,
        column_dim,
        U,
        V,
        ge_index,
        psd_index,
    )
end

function MOI.supports_constraint(
    ::Type{NormNuclearBridge{T}},
    ::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.NormNuclearCone},
) where {T}
    return true
end

function MOI.Bridges.added_constrained_variable_types(
    ::Type{<:NormNuclearBridge},
)
    return Tuple{Type}[]
end

function MOI.Bridges.added_constraint_types(
    ::Type{NormNuclearBridge{T,F,G,H}},
) where {T,F,G,H}
    return Tuple{Type,Type}[
        (F, MOI.GreaterThan{T}),
        (G, MOI.PositiveSemidefiniteConeTriangle),
    ]
end

function concrete_bridge_type(
    ::Type{<:NormNuclearBridge{T}},
    H::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.NormNuclearCone},
) where {T}
    S = MOI.Utilities.scalar_type(H)
    F = MOI.Utilities.promote_operation(
        -,
        T,
        S,
        MOI.Utilities.promote_operation(
            /,
            T,
            MOI.Utilities.promote_operation(+, T, T, T),
            T,
        ),
    )
    G = MOI.Utilities.promote_operation(vcat, T, MOI.VectorOfVariables, H)
    return NormNuclearBridge{T,F,G,H}
end

function MOI.get(bridge::NormNuclearBridge, ::MOI.NumberOfVariables)::Int64
    return length(bridge.U) + length(bridge.V)
end

function MOI.get(bridge::NormNuclearBridge, ::MOI.ListOfVariableIndices)
    return vcat(bridge.U, bridge.V)
end

function MOI.get(
    ::NormNuclearBridge{T,F,G,H},
    ::MOI.NumberOfConstraints{F,MOI.GreaterThan{T}},
)::Int64 where {T,F,G,H}
    return 1
end

function MOI.get(
    ::NormNuclearBridge{T,F,G,H},
    ::MOI.NumberOfConstraints{G,MOI.PositiveSemidefiniteConeTriangle},
)::Int64 where {T,F,G,H}
    return 1
end

function MOI.get(
    bridge::NormNuclearBridge{T,F,G,H},
    ::MOI.ListOfConstraintIndices{F,MOI.GreaterThan{T}},
) where {T,F,G,H}
    return [bridge.ge_index]
end

function MOI.get(
    bridge::NormNuclearBridge{T,F,G,H},
    ::MOI.ListOfConstraintIndices{G,MOI.PositiveSemidefiniteConeTriangle},
) where {T,F,G,H}
    return [bridge.psd_index]
end

function MOI.delete(model::MOI.ModelLike, bridge::NormNuclearBridge)
    MOI.delete(model, bridge.ge_index)
    MOI.delete(model, bridge.psd_index)
    MOI.delete(model, bridge.U)
    MOI.delete(model, bridge.V)
    return
end

function MOI.get(
    model::MOI.ModelLike,
    ::MOI.ConstraintFunction,
    bridge::NormNuclearBridge{T,F,G,H},
) where {T,F,G,H}
    ge_func = MOI.get(model, MOI.ConstraintFunction(), bridge.ge_index)
    psd_func = MOI.Utilities.eachscalar(
        MOI.get(model, MOI.ConstraintFunction(), bridge.psd_index),
    )
    column_dim = bridge.column_dim
    side_dim = bridge.row_dim + column_dim
    t = MOI.Utilities.operate(
        +,
        T,
        ge_func,
        MOI.Utilities.operate(
            /,
            T,
            MOI.Utilities.operate(
                +,
                T,
                [psd_func[MOI.Utilities.trimap(i, i)] for i in 1:side_dim]...,
            ),
            T(2),
        ),
    )
    t = MOI.Utilities.remove_variable(
        MOI.Utilities.remove_variable(t, bridge.U),
        bridge.V,
    )
    X = psd_func[[
        MOI.Utilities.trimap(i, j) for j in 1:bridge.column_dim for
        i in (bridge.column_dim+1):side_dim
    ]]
    return MOI.Utilities.convert_approx(H, MOI.Utilities.operate(vcat, T, t, X))
end

function MOI.get(
    ::MOI.ModelLike,
    ::MOI.ConstraintSet,
    bridge::NormNuclearBridge,
)
    return MOI.NormNuclearCone(bridge.row_dim, bridge.column_dim)
end

function MOI.supports(
    model::MOI.ModelLike,
    attr::MOI.ConstraintDualStart,
    ::Type{NormNuclearBridge{T,F,G,H}},
) where {T,F,G,H}
    ci_1 = MOI.ConstraintIndex{F,MOI.GreaterThan{T}}
    ci_2 = MOI.ConstraintIndex{G,MOI.PositiveSemidefiniteConeTriangle}
    return MOI.supports(model, attr, ci_1) && MOI.supports(model, attr, ci_2)
end

function MOI.get(
    model::MOI.ModelLike,
    ::MOI.ConstraintPrimal,
    bridge::NormNuclearBridge,
)
    ge_primal = MOI.get(model, MOI.ConstraintPrimal(), bridge.ge_index)
    psd_primal = MOI.get(model, MOI.ConstraintPrimal(), bridge.psd_index)
    side_dim = bridge.row_dim + bridge.column_dim
    t =
        ge_primal +
        sum(psd_primal[MOI.Utilities.trimap(i, i)] for i in 1:side_dim) / 2
    X = psd_primal[[
        MOI.Utilities.trimap(i, j) for j in 1:bridge.column_dim for
        i in (bridge.column_dim+1):side_dim
    ]]
    return vcat(t, X)
end

# Given t is dual on GreaterThan constraint and [U X'; X V] is dual on PSD
# constraint, the dual on NormNuclearCone constraint is
# (t, 2X) in NormNuclearCone.
function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintDual,MOI.ConstraintDualStart},
    bridge::NormNuclearBridge,
)
    t = MOI.get(model, attr, bridge.ge_index)
    psd_dual = MOI.get(model, attr, bridge.psd_index)
    side_dim = bridge.row_dim + bridge.column_dim
    X =
        2 * psd_dual[[
            MOI.Utilities.trimap(i, j) for j in 1:bridge.column_dim for
            i in (bridge.column_dim+1):side_dim
        ]]
    return vcat(t, X)
end

function MOI.set(
    model::MOI.ModelLike,
    ::MOI.ConstraintDualStart,
    bridge::NormNuclearBridge{T},
    value,
) where {T}
    MOI.set(model, MOI.ConstraintDualStart(), bridge.ge_index, value[1])
    column_dim = bridge.column_dim
    side_dim = bridge.row_dim + column_dim
    dual = zeros(T, div(side_dim * (side_dim + 1), 2))
    for i in 1:side_dim
        dual[MOI.Utilities.trimap(i, i)] = value[1]
    end
    X_idx = 2
    for j in 1:column_dim, i in (column_dim+1):side_dim
        dual[MOI.Utilities.trimap(i, j)] = value[X_idx] / 2
        X_idx += 1
    end
    MOI.set(model, MOI.ConstraintDualStart(), bridge.psd_index, dual)
    return
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
module Reader
    using GZip
    using BinDeps

    export

        readtrainimages,
        readtestimages,
        readtrainimages,
        readtestlabels,

        download_helper

    # Constants

    const IMAGEOFFSET = 16
    const LABELOFFSET = 8

    const TRAINIMAGES = "train-images-idx3-ubyte.gz"
    const TRAINLABELS = "train-labels-idx1-ubyte.gz"
    const TESTIMAGES  = "t10k-images-idx3-ubyte.gz"
    const TESTLABELS  = "t10k-labels-idx1-ubyte.gz"

    # Includes

    include("readheader.jl")
    include("readimages.jl")
    include("readlabels.jl")
    include("download.jl")
end

ext: jl


text:
push!(LOAD_PATH, pwd())
using AtomicOpt
using Printf
using LinearAlgebra
import Random: seed!, randperm, shuffle
using Distributed
using BSON: @save, @load

function randop(n::Integer, family=:orth)
    if family == :orth
        Q = Matrix(qr(randn(n,n)).Q)
    end
    return Q
end

function randnsparse(n, nnz)
    p = randperm(n)[1:nnz]
    x = zeros(n)
    x[p] = rand(nnz)
    return x
end

struct Signal{opT, atomT}
    n::Int64
    nnz::Int64
    x::Vector{Float64}
    Q::opT
    A::atomT
end

function Signal(n, nnz; family=:orth)
    x = randnsparse(n, nnz)
    Q = randop(n, family)
    λ = norm(x, 1)
    A = λ*(Q*OneBall(n; maxrank = nnz))
    return Signal(n, nnz, x, Q, A)
end

rotate(s::Signal) = s.Q*s.x

function RecordError(m, n, k, nnz, t)
    signals = ntuple(d->Signal(n, nnz), k)
    M = rand(m, n)
    b = M*sum(rotate(s) for s in signals)
    A = mapreduce(s->s.A, ×, signals)
    x, normr = level_set(M, b, A,  α = 0.0, tol = 1e-3, maxIts=10*length(b), logger=0)
    err = 0
    for (i,(x,s)) in enumerate(zip(x, signals))
        err1 = norm(x-s.Q*s.x);
        if err1 > err
            err = err1
        end
    end
    @printf("  m = %3d, n = %3d, k = %1d, nnz = %1d, t = %2d, error = %8.2e\n", m, n, k, nnz, t, err)
    flush(stdout)
    return err
end

function RecordError2(m, n, k, nnz, α, t)
    signals = ntuple(d->Signal(n, nnz), k)
    M = randn(m, n)
    b = M*sum(rotate(s) for s in signals)
    u = randn(m); u .*= α / norm(u)
    b = b + u
    A = mapreduce(s->s.A, ×, signals)
    x, normr = level_set(M, b, A,  α = α^2/2, tol = 1e-3, maxIts=10*length(b), logger=0)
    err = 0
    for (i,(x,s)) in enumerate(zip(x, signals))
        err1 = norm(x-s.Q*s.x);
        if err1 > err
            err = err1
        end
    end
    @printf("  α = %8.2e, t = %3d\n", α, t)
    flush(stdout)
    return err
end

# function AvgError(m, n, k, nnz)
#     out = @distributed (+) for i = 1:100
#         RecordError(m, n, k, nnz)
#     end
#     @printf("  m = %7d, n = %7d, k = %7d, nnz = %7d, error = %8.2e\n", m, n, k, nnz, out/100)
#     flush(stdout)
#     return out/100
# end

# @show RecordError(200, 200, 4, 5, 1)

# m = 125
# n = 200
# k = 3
# nnz = 5
# J = CartesianIndices((31, 31, 50))
# @printf("Start E2!\n")
# E2 = pmap(i->RecordError(500 - (i[1] - 1)*15, 50+(i[2]-1)*15, 2, nnz, i[3]), J)
# @printf("Start E3!\n")
# E3 = pmap(i->RecordError(500 - (i[1] - 1)*15, 50+(i[2]-1)*15, 3, nnz, i[3]), J)
# @printf("Start E4!\n")
# E4 = pmap(i->RecordError(500 - (i[1] - 1)*15, 50+(i[2]-1)*15, 4, nnz, i[3]), J)
# J = CartesianIndices((200, 50))
# E = pmap(i->RecordError2(m, n, k, nnz, 0.01*i[1], i[2]), J)

# # save
# @save "./examples/PhaseTransition/NoisyMeasurement.bson" E
# @save "./examples/PhaseTransition/NoiselessMeasurement.bson" E2 E3 E4

nnz = 3
n = 1000
J = CartesianIndices((19, 9, 10))
E = pmap(i->RecordError(1000 - (i[1] - 1)*50, n, i[2]+1, nnz, i[3]), J)
@save "./examples/PhaseTransition/NoiselessMeasurement2.bson" E


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("Parallel_Computing")

# # GPUs
#
# The graphics processor in your computer is _itself_ like a mini-computer highly
# tailored for massively and embarassingly parallel operations (like computing how light will bounce
# off of every point on a 3D mesh of triangles).
#
# Of course, recently their utility in other applications has become more clear
# and thus the GPGPU was born.
#
# Just like how we needed to send data to other processes, we need to send our
# data to the GPU to do computations there.

#-

# ## How is a GPU different from a CPU?
#
# This is what a typical consumer CPU looks like:
#
# ![](https://raw.githubusercontent.com/JuliaComputing/JuliaAcademyData.jl/master/courses/Parallel_Computing/images/i7.jpg)
#
# And this is what a GPU looks like:
#
# ![](https://raw.githubusercontent.com/JuliaComputing/JuliaAcademyData.jl/master/courses/Parallel_Computing/images/GK110.jpg)
#
# Each SMX isn't just one "core", each is a _streaming multiprocessor_ capable of running hundreds of threads simultaneously itself.  There are so many threads, in fact, that you reason about them in groups of 32 — called a "warp."  No, no [that warp](https://www.google.com/search?tbm=isch&q=warp&tbs=imgo:1), [this one](https://www.google.com/search?tbm=isch&q=warp%20weaving&tbs=imgo:1).
#
# The card above supports up to 6 warps per multiprocessor, with 32 threads each, times 15 multiprocessors... 2880 threads at a time!
#
# Also note the memory interfaces.
#
# --------------
#
# Each thread is relatively limited — and a warp is almost like a SIMD unit that supports branching. Except it's still only executing one instruction even after a branch:
#
# ![](https://raw.githubusercontent.com/JuliaComputing/JuliaAcademyData.jl/master/courses/Parallel_Computing/images/warp-branch.png)

#-

# You can inspect the installed GPUs with nvidia-smi:

#nb ;nvidia-smi
#jl run(`nvidia-smi`)

# ## Example
#
# The deep learning MNIST example: https://fluxml.ai/experiments/mnist/
#
# This is how it looks on the CPU:

using Flux, Flux.Data.MNIST, Statistics
using Flux: onehotbatch, onecold, crossentropy, throttle
using Base.Iterators: repeated, partition

imgs = MNIST.images()
labels = onehotbatch(MNIST.labels(), 0:9)

## Partition into batches of size 32
train = [(cat(float.(imgs[i])..., dims = 4), labels[:,i])
         for i in partition(1:60_000, 32)]
## Prepare test set (first 1,000 images)
tX = cat(float.(MNIST.images(:test)[1:1000])..., dims = 4)
tY = onehotbatch(MNIST.labels(:test)[1:1000], 0:9)

m = Chain(
        Conv((3, 3), 1=>32, relu),
        Conv((3, 3), 32=>32, relu),
        x -> maxpool(x, (2,2)),
        Conv((3, 3), 32=>16, relu),
        x -> maxpool(x, (2,2)),
        Conv((3, 3), 16=>10, relu),
        x -> reshape(x, :, size(x, 4)),
        Dense(90, 10), softmax)

loss(x, y) = crossentropy(m(x), y)
accuracy(x, y) = mean(onecold(m(x)) .== onecold(y))
## opt = ADAM() # <-- Move Flux.params(m) here!
## @time Flux.train!(loss, Flux.params(m), train[1:10], opt, cb = () -> @show(accuracy(tX, tY)))
opt = ADAM(Flux.params(m), ) # <-- Move Flux.params(m) here!
Flux.train!(loss, train[1:1], opt, cb = () -> @show(accuracy(tX, tY)))
@time Flux.train!(loss, train[1:10], opt, cb = () -> @show(accuracy(tX, tY)))

# Now let's re-do it on a GPU. "All" it takes is moving the data there with `gpu`!

include(datapath("scripts/fixupCUDNN.jl")) # JuliaBox uses an old version of CuArrays; this backports a fix for it
gputrain = gpu.(train[1:10])
gpum = gpu(m)
gputX = gpu(tX)
gputY = gpu(tY)
gpuloss(x, y) = crossentropy(gpum(x), y)
gpuaccuracy(x, y) = mean(onecold(gpum(x)) .== onecold(y))
gpuopt = ADAM(Flux.params(gpum), )
Flux.train!(gpuloss, gpu.(train[1:1]), gpuopt, cb = () -> @show(gpuaccuracy(gputX, gputY)))
@time Flux.train!(gpuloss, gputrain, gpuopt, cb = () -> @show(gpuaccuracy(gputX, gputY)))

# ## Defining your own GPU kernels
#
# So that's leveraging Flux's ability to work with GPU arrays — which is magical
# and awesome — but you don't always have a library to lean on like that.
# How might you define your own GPU kernel?
#
# Recall the monte carlo pi example:

function serialpi(n)
    inside = 0
    for i in 1:n
        x, y = rand(), rand()
        inside += (x^2 + y^2 <= 1)
    end
    return 4 * inside / n
end

# How could we express this on the GPU?

using CuArrays.CURAND
function findpi_gpu(n)
    4 * sum(curand(Float64, n).^2 .+ curand(Float64, n).^2 .<= 1) / n
end
findpi_gpu(10_000_000)

#-

using BenchmarkTools
@btime findpi_gpu(10_000_000)
@btime serialpi(10_000_000)

# That leans on broadcast to build the GPU kernel — and is creating three arrays
# in the process — but it's still much faster than our serial pi from before.

#-

# In general, using CuArrays and broadcast is one of the best ways to just
# get everything to work. If you really want to get your hands dirty, you
# can use [CUDAnative.jl](https://github.com/JuliaGPU/CUDAnative.jl) to manually specify exactly how everything works,
# but be forewarned, it's not for the [faint at heart](https://github.com/JuliaGPU/CUDAnative.jl/blob/master/examples/reduce/reduce.jl)! (If you've done CUDA
# programming in C or C++, it's very similar.)

max_stars_repo_path: Courses/Parallel_Computing/080 GPUs.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_path: Courses/Parallel_Computing/080 GPUs.jl

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_path: Courses/Parallel_Computing/080 GPUs.jl

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
#Test if tf listner works correctly
using RobotOS
#
# # Initialize node
# init_node("listner", anonymous=true) # delete

path = @__DIR__

# Need to run the listener and broadcaster in parralllel
script = (`julia --project=@. $path/tf_listener.jl` &
          `julia --project=@. $path/tf_broadcaster.jl` &
          `julia --project=@. $path/tf_applied.jl`);
Base.run(script);


text:
type SigmoidLayer <: Layer
    inputActivation::Vol
    outputActivation::Vol

    function SigmoidLayer(dimensions)
        new(Vol(dimensions, 0.0), Vol(dimensions, 0.0))
    end
end

function forward(self::SigmoidLayer, inputActivation::Vol)
    @assert inputSize(self) == outputSize(self)
    self.inputActivation = inputActivation
    self.outputActivation.w = logistic(inputActivation.w)
    return self.outputActivation
end

function backward(self::SigmoidLayer)
    @assert inputSize(self) == outputSize(self)
    self.inputActivation.dw = self.outputActivation.w .* (1.0 .- self.outputActivation.w) .* self.outputActivation.dw
end

gradientStates(self::SigmoidLayer) = []
weights(self::SigmoidLayer) = []

type RectifiedLinearLayer
    inputActivation::Vol
    outputActivation::Vol

    function RectifiedLinearLayer(dimensions)
        new(Vol(dimensions, 0.0), Vol(dimensions, 0.0))
    end
end

function forward(self::RectifiedLinearLayer, inputActivation::Vol)
    @assert size(inputActivation) == size(self.inputActivation)
    self.inputActivation = inputActivation
    self.outputActivation = max(self.inputActivation, 0)
    return self.outputActivation
end

function backward(self::RectifiedLinearLayer)
    @assert inputSize(self) == outputSize(self)
    fill!(self.inputActivation.dw, 0.0)
    self.inputActivation.dw = (0 .< self.outputActivation.w) .* self.outputActivation.dw
end

gradientStates(self::RectifiedLinearLayer) = []


text:
@testset "open() for DataSet" begin
    proj = DataSets.load_project("Data.toml")

    blob_ds = dataset(proj, "embedded_blob")
    @test open(blob_ds) isa Blob
    @test open(String, blob_ds) == "\0\0\0\0\0\0E@"
    @test read(open(blob_ds), Float64) === 42.0

    @test open(IO, blob_ds) do io
        read(io, String)
    end == "\0\0\0\0\0\0E@"

    @context begin
        @test @!(open(String, blob_ds)) == "\0\0\0\0\0\0E@"

        blob = @! open(blob_ds)
        @test blob isa Blob
        @test @!(open(String, blob)) == "\0\0\0\0\0\0E@"

        @test read(blob, Float64) === 42.0
        @test read(blob) == UInt8[0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x45, 0x40]
    end

    tree_ds = dataset(proj, "embedded_tree")
    @test open(tree_ds) isa BlobTree
    @test open(String, open(tree_ds)[path"d01/a.txt"]) == "1 a content"
    @test open(String, open(tree_ds)[path"d02/b.txt"]) == "2 b content"
    @context begin
        tree = @! open(tree_ds)
        @test tree isa BlobTree

        @test isdir(tree)
        @test !isfile(tree)

        @test readdir(tree) == ["d01", "d02", "d03", "d04"]
        @test readdir(tree["d01"]) == ["a.txt", "b.txt"]

        @test !isdir(tree[path"d01/a.txt"])
        @test isfile(tree[path"d01/a.txt"])

        @test_throws ErrorException tree[path"nonexistent/a/b"]
        @test_throws ErrorException tree["nonexistent"]

        @test @!(open(String, tree[path"d01/a.txt"])) == "1 a content"
        @test @!(open(String, tree[path"d02/b.txt"])) == "2 b content"
    end
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

using Base.Test

@test @__LINE__ == 5

include("test_sourcepath.jl")
thefname = "the fname!//\\&\1*"
@test include_string("include_string_test() = @__FILE__", thefname)() == Base.source_path()
@test include_string("Base.source_path()", thefname) == Base.source_path()
@test basename(@__FILE__) == "loading.jl"
@test isabspath(@__FILE__)

@test isdir(@__DIR__)
@test @__DIR__() == dirname(@__FILE__)

# Issue #5789 and PR #13542:
mktempdir() do dir
    cd(dir) do
        let true_filename = "cAsEtEsT.jl", lowered_filename="casetest.jl"
            touch(true_filename)
            @test Base.isfile_casesensitive(true_filename)
            @test !Base.isfile_casesensitive(lowered_filename)

            # check that case-sensitivity only applies to basename of a path:
            if isfile(lowered_filename) # case-insensitive filesystem
                mkdir("cAsEtEsT")
                touch(joinpath("cAsEtEsT", true_filename))
                @test Base.isfile_casesensitive(joinpath("casetest", true_filename))
                @test !Base.isfile_casesensitive(joinpath("casetest", lowered_filename))
            end
        end

        # Test Unicode normalization; pertinent for OS X
        let nfc_name = "\U00F4.jl"
            touch(nfc_name)
            @test Base.isfile_casesensitive(nfc_name)
        end
    end
end

let paddedname = "Ztest_sourcepath.jl"
    filename = SubString(paddedname, 2, length(paddedname))
    @test Base.find_in_path(filename) == abspath(paddedname[2:end])
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
module CLIBucket

using ...SimpleTasks.Types

import SimpleTasks.Services.Bucket

export CLIBucketService

type Provider
    prefix::Cmd
    command::Cmd
end

type CLIBucketService <: BucketService
    provider::Provider
    name::AbstractString

    CLIBucketService(provider::Provider, name::AbstractString) =
        check_reachable(provider, name) && new(provider, name)
end

const MEGABYTE = 1024 * 1024
const WRITE_BYTES = 50 * MEGABYTE
const READ_BYTES = 50 * MEGABYTE

function check_reachable(provider::Provider, bucket_name::AbstractString)
    ls_cmd = `$(provider.command) ls $(provider.prefix)/$bucket_name`
    try
        # pipeline into DevNull to squelch stdout
        run(pipeline(ls_cmd, stdout=DevNull))
    catch
        throw(ArgumentError("Unable to access bucket \"$bucket_name\" " *
            "with $ls_cmd"))
    end
    return true
end

function Bucket.download(bucket::CLIBucketService,
    remote_file::AbstractString,
    local_file::Union{AbstractString, IO, Void}=nothing)

    download_cmd = `$(bucket.provider.command) cp
        "$(bucket.provider.prefix)/$(bucket.name)/$remote_file" -`

    s3_output = Pipe()
    # open the cmd in write mode. this automatically takes the 2nd arg
    # (stdio) and uses it as redirection of STDOUT.
    (s3_input, process) = open(download_cmd, "w", s3_output)
    close(s3_output.in)

    # Stop gap measure to check if we couldn't locate the file
    timedwait(() -> process_exited(process), 1.0)
    if !process_running(process) && !success(process)
        error("Error downloading $remote_file using command $download_cmd")
    end

    if local_file != nothing
        if isa(local_file, ASCIIString)
            local_file = open(local_file, "w")
        end
        while !eof(s3_output)
            write(local_file, readbytes(s3_output, WRITE_BYTES))
        end
        if !success(process)
            error("Download did not complete cleanly")
        end
    else
        local_file = s3_output
    end

    return local_file
end

function Bucket.upload(bucket::CLIBucketService,
        local_file::Union{AbstractString, IO}, remote_file::AbstractString)

    upload_cmd = `$(bucket.provider.command) cp -
        "$(bucket.provider.prefix)/$(bucket.name)/$remote_file"`

    s3_input = Pipe()
    # open the cmd in read mode. this automatically takes the 2nd arg
    # (stdio) and uses it as redirection of STDIN.
    (s3_output, process) = open(upload_cmd, "r", s3_input)

    # Stop gap measure to check if we couldn't upload the file
    timedwait(() -> process_exited(process), 1.0)
    if !process_running(process) && !success(process)
        error("Error uploading $remote_file using command $upload_cmd")
    end

    if isa(local_file, AbstractString)
        local_file = open(local_file, "r")
    end

    while !eof(local_file)
        write(s3_input, readbytes(local_file, READ_BYTES))
    end

    close(s3_input.in)

    if !success(process)
        error("Upload did not complete cleanly")
    end
end

function Bucket.delete(bucket::CLIBucketService,
    remote_file::AbstractString)
    remove_cmd = `$(bucket.provider.command) rm
        "$(bucket.provider.prefix)/$(bucket.name)/$remote_file"`
    try
        # pipeline into DevNull to squelch stdout
        run(remove_cmd)
    catch
        throw(ArgumentError("Unable to delete $remote_file from bucket " *
            "\"$(bucket.name)\" with command $remove_cmd"))
    end
    return true
end

end # module CLIBucket

max_forks_repo_name: seung-lab/SimpleTasks


text:
"""
This is borrowed from [DataFrames.jl]().
Return between 0 and 2 names from `colnames` closest to `name`.
`colnames` : some iterable collection of symbols
"""
function fuzzymatch(colnames, name::Symbol)
    ucname = uppercase(string(name))
    dist = [(levenshtein(uppercase(string(x)), ucname), x) for x in colnames]
    sort!(dist)
    c = [count(x -> x[1] <= i, dist) for i in 0:2]
    maxd = max(0, searchsortedlast(c, 8) - 1)
    return [s for (d, s) in dist if d <= maxd]
end

"""
Return a nice-ish error message if the Symbol `name` isn't a column name in `table`, otherwise a zero-length string.
"""
function checkcol(table, name::Symbol)
    i = Tables.columnindex(table, name)
    if i == 0 # if no such column
        names = Tables.columnnames(table)
        nearestnames = join(fuzzymatch(names, name),", " )
        return "There isn't a variable called '$name' in your data; the nearest names appear to be: $nearestnames"
    end

    nrows = length(Tables.getcolumn(table, name))
    if nrows == 0
        return "Column '$name' is empty."
    end

    return ""
end

"""
Check that each name in the given model `f` exists in the data source `t` and return a message if not. Return a zero string otherwise.
`t` is something that implements the `Tables` interface.
"""
function checknamesexist(f::FormulaTerm, t)
    if ! Tables.istable(t)
        throw(ArgumentError( "$(typeof(t)) isn't a valid Table type" ))
    end
    for n in StatsModels.termvars(f)
        msg = checkcol(t, n)
        if msg != ""
            return msg
        end
    end
    return ""
end


text:
module ScientificComputing

using Pluto

# Just for the sake of it
start_slides() = begin
    Pluto.run(notebook=joinpath(@__DIR__, "..", "slides","notebook.jl"))
end

export start_slides
end  # module ScientificComputing


text:
@testset "Conv Layers" begin
    in_channel = 3
    out_channel = 5
    N = 4
    T = Float32

    adj1 =  [0 1 0 1
            1 0 1 0
            0 1 0 1
            1 0 1 0]

    g1 = GNNGraph(adj1,
            ndata=rand(T, in_channel, N),
            graph_type=GRAPH_T)

    adj_single_vertex =  [0 0 0 1
                          0 0 0 0
                          0 0 0 1
                          1 0 1 0]

    g_single_vertex = GNNGraph(adj_single_vertex,
                                ndata=rand(T, in_channel, N),
                                graph_type=GRAPH_T)

    test_graphs = [g1, g_single_vertex]

    @testset "GCNConv" begin
        l = GCNConv(in_channel => out_channel)
        for g in test_graphs
            test_layer(l, g, rtol=1e-5, outsize=(out_channel, g.num_nodes))
        end

        l = GCNConv(in_channel => out_channel, tanh, bias=false)
        for g in test_graphs
            test_layer(l, g, rtol=1e-5, outsize=(out_channel, g.num_nodes))
        end

        l = GCNConv(in_channel => out_channel, add_self_loops=false)
        test_layer(l, g1, rtol=1e-5, outsize=(out_channel, g1.num_nodes))

        @testset "edge weights" begin
            s = [2,3,1,3,1,2]
            t = [1,1,2,2,3,3]
            w = [1,2,3,4,5,6]
            g = GNNGraph((s, t, w), graph_type=GRAPH_T)
            x = ones(1, g.num_nodes)
            l = GCNConv(1 => 1, add_self_loops=false, use_edge_weight=true)
            l.weight .= 1
            d = degree(g, dir=:in)
            y = l(g, x)
            @test y[1,1] ≈ w[1] / √(d[1]*d[2]) + w[2] / √(d[1]*d[3])
            @test y[1,2] ≈ w[3] / √(d[2]*d[1]) + w[4] / √(d[2]*d[3])
            if GRAPH_T == :coo
                @test y ≈ l(g, x, w)
            end
        end
    end

    @testset "ChebConv" begin
        k = 2
        l = ChebConv(in_channel => out_channel, k)
        @test size(l.weight) == (out_channel, in_channel, k)
        @test size(l.bias) == (out_channel,)
        @test l.k == k
        for g in test_graphs
            g = add_self_loops(g)
            test_layer(l, g, rtol=1e-5, test_gpu=false, outsize=(out_channel, g.num_nodes))
            if TEST_GPU
                @test_broken test_layer(l, g, rtol=1e-5, test_gpu=true, outsize=(out_channel, g.num_nodes))
            end
        end

        @testset "bias=false" begin
            @test length(Flux.params(ChebConv(2=>3, 3))) == 2
            @test length(Flux.params(ChebConv(2=>3, 3, bias=false))) == 1
        end
    end

    @testset "GraphConv" begin
        l = GraphConv(in_channel => out_channel)
        for g in test_graphs
            test_layer(l, g, rtol=1e-5, outsize=(out_channel, g.num_nodes))
        end

        l = GraphConv(in_channel => out_channel, relu, bias=false, aggr=mean)
        for g in test_graphs
            test_layer(l, g, rtol=1e-5, outsize=(out_channel, g.num_nodes))
        end

        @testset "bias=false" begin
            @test length(Flux.params(GraphConv(2=>3))) == 3
            @test length(Flux.params(GraphConv(2=>3, bias=false))) == 2
        end
    end

    @testset "GATConv" begin

        for heads in (1, 2), concat in (true, false)
            l = GATConv(in_channel => out_channel; heads, concat)
            for g in test_graphs
                test_layer(l, g, rtol=1e-4,
                    outsize=(concat ? heads*out_channel : out_channel, g.num_nodes))
            end
        end

        @testset "bias=false" begin
            @test length(Flux.params(GATConv(2=>3))) == 3
            @test length(Flux.params(GATConv(2=>3, bias=false))) == 2
        end
    end

    @testset "GatedGraphConv" begin
        num_layers = 3
        l = GatedGraphConv(out_channel, num_layers)
        @test size(l.weight) == (out_channel, out_channel, num_layers)

        for g in test_graphs
            test_layer(l, g, rtol=1e-5, outsize=(out_channel, g.num_nodes))
        end
    end

    @testset "EdgeConv" begin
        l = EdgeConv(Dense(2*in_channel, out_channel), aggr=+)
        for g in test_graphs
            test_layer(l, g, rtol=1e-5, outsize=(out_channel, g.num_nodes))
        end
    end

    @testset "GINConv" begin
        nn = Dense(in_channel, out_channel)

        l = GINConv(nn, 0.01f0, aggr=mean)
        for g in test_graphs
            test_layer(l, g, rtol=1e-5, outsize=(out_channel, g.num_nodes))
        end

        @test !in(:eps, Flux.trainable(l))
    end

    @testset "NNConv" begin
        edim = 10
        nn = Dense(edim, out_channel * in_channel)

        l = NNConv(in_channel => out_channel, nn, tanh, bias=true, aggr=+)
        for g in test_graphs
            g = GNNGraph(g, edata=rand(T, edim, g.num_edges))
            test_layer(l, g, rtol=1e-5, outsize=(out_channel, g.num_nodes))
        end
    end

    @testset "SAGEConv" begin
        l = SAGEConv(in_channel => out_channel)
        @test l.aggr == mean

        l = SAGEConv(in_channel => out_channel, tanh, bias=false, aggr=+)
        for g in test_graphs
            test_layer(l, g, rtol=1e-5, outsize=(out_channel, g.num_nodes))
        end
    end

    @testset "ResGatedGraphConv" begin
        l = ResGatedGraphConv(in_channel => out_channel, tanh, bias=true)
        for g in test_graphs
            test_layer(l, g, rtol=1e-5, outsize=(out_channel, g.num_nodes))
        end
    end

    @testset "CGConv" begin
        edim = 10
        l = CGConv((in_channel, edim) => out_channel, tanh, residual=false, bias=true)
        for g in test_graphs
            g = GNNGraph(g, edata=rand(T, edim, g.num_edges))
            test_layer(l, g, rtol=1e-5, outsize=(out_channel, g.num_nodes))
        end
    end

    @testset "AGNNConv" begin
        l = AGNNConv()
        l.β == [1f0]
        for g in test_graphs
            test_layer(l, g, rtol=1e-5, outsize=(in_channel, g.num_nodes))
        end
    end
end


text:
"""
Implementation of Genetic Algorithm

The constructor takes following keyword arguments:

- `populationSize`: The size of the population
- `crossoverRate`: The fraction of the population at the next generation, not including elite children, that is created by the crossover function.
- `mutationRate`: Probability of chromosome to be mutated
- `ɛ`/`epsilon`: Positive integer specifies how many individuals in the current generation are guaranteed to survive to the next generation. Floating number specifies fraction of population.
- `selection`: [Selection](@ref) function (default: [`tournament`](@ref))
- `crossover`: [Crossover](@ref) function (default: [`genop`](@ref))
- `mutation`: [Mutation](@ref) function (default: [`genop`](@ref))
- `metrics` is a collection of convergence metrics.
"""
struct GA{T1,T2,T3} <: AbstractOptimizer
    populationSize::Int
    crossoverRate::Float64
    mutationRate::Float64
    ɛ::Real
    selection::T1
    crossover::T2
    mutation::T3
    metrics::ConvergenceMetrics

    GA(; populationSize::Int=50, crossoverRate::Float64=0.8, mutationRate::Float64=0.1,
        ɛ::Real=0, epsilon::Real=ɛ,
        selection::T1=tournament(2),
        crossover::T2=genop,
        mutation::T3=genop,
        metrics = ConvergenceMetric[AbsDiff(1e-12)]) where {T1, T2, T3} =
        new{T1,T2,T3}(populationSize, crossoverRate, mutationRate, epsilon, selection, crossover, mutation, metrics)
end
population_size(method::GA) = method.populationSize
default_options(method::GA) = (iterations=1000,)
summary(m::GA) = "GA[P=$(m.populationSize),x=$(m.crossoverRate),μ=$(m.mutationRate),ɛ=$(m.ɛ)]"
show(io::IO,m::GA) = print(io, summary(m))

mutable struct GAState{T,IT} <: AbstractOptimizerState
    N::Int
    eliteSize::Int
    fitness::T
    fitpop::Vector{T}
    fittest::IT
end
value(s::GAState) = s.fitness
minimizer(s::GAState) = s.fittest

"""Initialization of GA algorithm state"""
function initial_state(method::GA, options, objfun, population)
    T = typeof(value(objfun))
    N = length(first(population))
    fitness = zeros(T, method.populationSize)

    # setup state values
    eliteSize = isa(method.ɛ, Int) ? method.ɛ : round(Int, method.ɛ * method.populationSize)

    # Evaluate population fitness
    fitness = map(i -> value(objfun, i), population)
    minfit, fitidx = findmin(fitness)

    # setup initial state
    return GAState(N, eliteSize, minfit, fitness, copy(population[fitidx]))
end

function update_state!(objfun, constraints, state, parents::AbstractVector{IT}, method::GA, options, itr) where {IT}
    populationSize = method.populationSize
    evaltype = options.parallelization
    rng = options.rng
    offspring = similar(parents)

    # select offspring
    selected = method.selection(state.fitpop, populationSize, rng=rng)

    # perform mating
    offspringSize = populationSize - state.eliteSize
    recombine!(offspring, parents, selected, method, offspringSize, rng=rng)

    # Elitism (copy population individuals before they pass to the offspring & get mutated)
    fitidxs = sortperm(state.fitpop)
    for i in 1:state.eliteSize
        subs = offspringSize+i
        offspring[subs] = copy(parents[fitidxs[i]])
    end

    # perform mutation
    mutate!(offspring, method, constraints, offspringSize, rng=rng)

    # calculate fitness of the population
    evaluate!(objfun, state.fitpop, offspring, constraints)

    # select the best individual
    minfit, fitidx = findmin(state.fitpop)
    state.fittest = offspring[fitidx]
    state.fitness = state.fitpop[fitidx]

    # replace population
    parents .= offspring

    return false
end

function recombine!(offspring, parents, selected, method, n=length(selected);
                    rng::AbstractRNG=default_rng())
    mates = ((i,i == n ? i-1 : i+1) for i in 1:2:n)
    for (i,j) in mates
        p1, p2 = parents[selected[i]], parents[selected[j]]
        if rand(rng) < method.crossoverRate
            offspring[i], offspring[j] = method.crossover(p1, p2, rng=rng)
        else
            offspring[i], offspring[j] = p1, p2
        end
    end

end

function mutate!(population, method, constraints, n=length(population);
                 rng::AbstractRNG=default_rng())
    for i in 1:n
        if rand(rng) < method.mutationRate
            method.mutation(population[i], rng=rng)
        end
        apply!(constraints, population[i])
    end
end

function evaluate!(objfun, fitness, population, constraints)
    # calculate fitness of the population
    value!(objfun, fitness, population)
    # apply penalty to fitness
    penalty!(fitness, constraints, population)
end


text:
@info "Train baseline Single Headed Attention Recurrent language model using enwik8 dataset..."

using Knet

include("../src/data.jl")
include("../src/model.jl")
include("../src/train.jl")

function evaluate()
    println()
    @info "Finished training, Starting evaluation ..."
    devloss = loss(model, ddev);
    println("Development set scores:    ", report_lm(devloss))
    testloss = loss(model, dtst);
    println("Test set scores:           ", report_lm(testloss))

    model_name = "small_mel.jld2"
    @info "Saving the model as $(model_name)"
    Knet.save(model_name, "model", model);
end

BATCHSIZE = 3 ; @show BATCHSIZE
BPTT = 256 ; @show BPTT
MEMSIZE = 512 ; @show MEMSIZE
EMSIZE = 512 ; @show EMSIZE

datadir = "../data/enwik8"
jld2dir = "../jld2/enwik8.jld2"
if !isfile(jld2dir)
    println("Reading data from directory: $datadir")
    println("Setting batch size to $BATCHSIZE")
    vocab = Vocab("$datadir/train.txt")
    trainfile = TextReader("$datadir/train.txt", vocab)
    validfile = TextReader("$datadir/valid.txt", vocab)
    testfile = TextReader("$datadir/test.txt", vocab)
    dtrn = TextData(trainfile, batchsize=BATCHSIZE, bptt=BPTT)
    ddev = TextData(validfile, batchsize=BATCHSIZE, bptt=BPTT, randomize = false)
    dtst = TextData(testfile, batchsize=BATCHSIZE, bptt=BPTT, randomize = false)
    println("Saving data to $jld2dir")
    Knet.save(jld2dir, "dtrn", dtrn, "dtst", dtst, "ddev", ddev)
else
    println("Loading data from $jld2dir")
    (dtrn, dtst, ddev) = Knet.load(jld2dir, "dtrn", "dtst", "ddev")
    vocab = dtrn.src.vocab
    if dtrn.batchsize != BATCHSIZE
        changebatchsize!(dtrn, BATCHSIZE)
        changebatchsize!(ddev, BATCHSIZE)
        changebatchsize!(dtst, BATCHSIZE)
    end;
    dtrn.bptt = BPTT
    dtst.bptt = BPTT
    ddev.bptt = BPTT
end;

println()
@info "Initializing the model and collecting training data..."
epochs, em_size, hidden_size, layers = 2, EMSIZE, (EMSIZE*4), 1
println("embedding size: ", em_size)
println("hidden size: ", hidden_size)
println("layers: ", layers)
println("Collecting training data...")
println("epochs: ", epochs)

ctrn = collect(dtrn)
trn = collect(flatten(collect(dtrn) for i in 1:epochs))
dev = collect(ddev)
mintrn = ctrn[1:20];

model = SHARNN(em_size, hidden_size, vocab, layers; num_max_positions=MEMSIZE, tie_weights=true);

# model = Knet.load("small_model.jld2", "model")

@show length(ctrn)

# testloss = loss(model, dtst);
# println("Test set scores:           ", report_lm(testloss))

println()
@info "Starting training, total iteration no: $(length(trn))"
initlamb!(model, length(trn); lr=0.002, warmup=(1200)/length(trn))

#  TODO: stop training at anytime using CTRL-C
model = train!(model, trn, dev, mintrn; report_iter=length(ctrn), update_per_n_batch=1)

atexit(evaluate)


text:
export prove_inclusion, set_includes

# TODO: check es for LogicalSymbols
evaluate(es::ExtensionalSet, ::Dict{LogicalSymbol, Bool}) = es
evaluate(is::IntensionalSet, ::Dict{LogicalSymbol, Bool}) = is

set_includes(es::ExtensionalSet{T}, sub::T) where {T} = sub ∈ elements(es)
function set_includes(sub::T, is::IntensionalSet) where {T <: AbstractExpression}
    pattern_matches = Dict{LogicalSymbol, AbstractExpression}()
    if !find_matches!(sub, transform(is), pattern_matches)
        return false
    end

    @info pattern_matches
    # evaluate(rule(is), pattern_matches)
end

function prove_inclusion(expr::LogicalExpression, substitution::Pair{LogicalSymbol, T}) where {T}
    target_set = right(expr)

    if operation(expr) != set_in || !(target_set isa MathematicalSet) || !(left(expr) isa LogicalSymbol)
        throw(ErrorException("Expression is not a set inclusion statement"))
    end

    if left(expr) != substitution.first
        throw(ErrorException("Substituted symbol and inclusion symbol are different ($(left(expr)) ≠ $(substitution.first))"))
    end

    set_includes(target_set, substitution.second)
end


text:
module NonDimensional

using LinearAlgebra: norm
using PrettyTables

export dynamic_pressure, force_coefficient, moment_coefficient, moment_coefficients, rate_coefficient, pressure_coefficient, aerodynamic_coefficients, print_coefficients, reynolds_number, print_derivatives

"""
    dynamic_pressure(ρ, V)

Compute the dynamic pressure given density ``\\rho`` and speed ``V``.
"""
dynamic_pressure(ρ, V) = 1/2 * ρ * V^2

"""
    force_coefficient(force, q, S)

Compute the non-dimensional force coefficient given force, dynamic pressure ``q``, and area ``S``.
"""
force_coefficient(force, q, S) = force / (q * S)

"""
    moment_coefficient(moment, q, S, L)

Compute the non-dimensional moment coefficient given force, dynamic pressure ``q``, area ``S``, and reference length ``L``.
"""
moment_coefficient(moment, q, S, L) = force_coefficient(moment, q, S) / L

"""
    rate_coefficient(Ω, V, L)

Compute the non-dimensional angular velocity coefficient given angular speed ``\\Omega``, reference speed ``V``, and reference length ``L``.
"""
rate_coefficient(Ω, V, L) = Ω * L / 2V

"""
    pressure_coefficient(force, ρ, V, S)

Compute the pressure coefficient given force, reference density ``\\rho``, speed ``V`` and area ``S``.
"""
pressure_coefficient(force, ρ, V, S) = force_coefficient(force, dynamic_pressure(ρ, V), S)

pressure_coefficient(mag, vels) = 1 - (norm(vels) / mag)^2

moment_coefficient(moment, q, S, b, c) = moment_coefficient.(moment, q, S, [b, c, b])

rate_coefficient(Ω, V, b, c) = rate_coefficient.(Ω, V, [b, c, b])

"""
    aerodynamic_coefficients(force, moment, Ω, V, S, b, c, ρ)

Compute the relevant aerodynamic coefficients given a net force, moment and angular rates ``\\Omega`` with reference speed ``V``, area ``S``, span ``b``, chord ``c``, and density ``\\rho``.
"""
function aerodynamic_coefficients(force, moment, Ω, V, S, b, c, ρ)
    q = dynamic_pressure(ρ, V)

    [ force_coefficient(force, q, S) 		 ;
      moment_coefficient(moment, q, S, b, c) ;
      rate_coefficient(Ω, V, b, c) 			 ]
end

function print_coefficients(nf_coeffs, ff_coeffs, name = ""; browser = false)
    coeffs = [ ifelse(length(nf_coeffs) == 11, ["CD", "CDp"], []); [ "CDi", "CY", "CL", "Cl", "Cm", "Cn", "p̄", "q̄", "r̄" ] ]
    data = [ coeffs nf_coeffs [ ff_coeffs; fill("—", 6) ] ]
    head = [ name, "Nearfield", "Farfield" ]
    h1 = Highlighter( (data,i,j) -> (j == 1), foreground = :blue, bold = true)
    if browser
        pretty_table(String, data, head, alignment = [:c, :c, :c], tf = tf_html_minimalist, backend = :html, highlighters = HTMLHighlighter( (data,i,j) -> (j == 1), HTMLDecoration(font_weight = "bold")), formatters = ft_round(8))
    else
        pretty_table(data, head, alignment = [:c, :c, :c], tf = tf_compact, highlighters = h1, vlines = :none, formatters = ft_round(8))
    end
end

function print_derivatives(derivs, name = ""; browser = false)
    coeffs = ["∂CD", "∂CY", "∂CL", "∂Cl", "∂Cm", "∂Cn"]
    nf_vars = [ "$name" "" "Nearfield" "Stability" "Derivatives" "" ; "" "∂α, 1/rad" "∂β, 1/rad" "∂p̄" "∂q̄" "∂r̄" ]
    nf_rows = [ coeffs derivs ]

    if browser
        pretty_table(String, nf_rows, nf_vars, alignment = :c, tf = tf_html_minimalist, backend = :html, highlighters = HTMLHighlighter( (data,i,j) -> (j == 1), HTMLDecoration(color = "blue", font_weight = "bold")), formatters = ft_round(8))
    else
        pretty_table(nf_rows, nf_vars, alignment = :c, tf = tf_compact, header_crayon = Crayon(bold = true), subheader_crayon = Crayon(foreground = :yellow, bold = true), highlighters = Highlighter( (data,i,j) -> (j == 1), foreground = :blue, bold = true), vlines = :none, formatters = ft_round(8))
    end
end

reynolds_number(ρ, V, c, μ) = ρ * V * c / μ
reynolds_number(V, c, ν) = V * c / ν

force(CF, q, S) = CF * q * S
moment(CM, q, S, c) = CM * q * S * c
moment(CM, q, S, b, c) = moment.(CM, q, S, [b, c, b])

end

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
# ---
# title: 92. Reverse Linked List II
# id: problem92
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Linked List
# link: <https://leetcode.com/problems/reverse-linked-list-ii/description/>
# hidden: true
# ---
#
# Reverse a linked list from position _m_ to _n_. Do it in one-pass.
#
# **Note:  **1 ≤ _m_ ≤ _n_ ≤ length of list.
#
# **Example:**
#
#
#
#     Input: 1->2->3->4->5->NULL, _m_ = 2, _n_ = 4
#     Output: 1->4->3->2->5->NULL
#
#
#
## @lc code=start
using LeetCode

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


text:
using Test
using REPL
using Cthulhu
using Revise

if !isdefined(@__MODULE__, :fake_terminal)
    @eval begin
        includet(joinpath(@__DIR__, "FakeTerminals.jl"))   # FIXME change to include
        using .FakeTerminals
    end
end

cread1(io) = readuntil(io, '↩'; keep=true)
cread(io) = cread1(io) * cread1(io)

# For multi-call sites
module MultiCall
fmulti(::Int32) = 1
fmulti(::Float32) = rand(Float32)
fmulti(::Char) = 3
callfmulti(c) = fmulti(c[])
end

const keydict = Dict(:up => "\e[A",
                     :down => "\e[B",
                     :enter => '\r')

@testset "Terminal" begin
    if isdefined(Base, :active_repl)
        @test Cthulhu.default_terminal() isa REPL.Terminals.TTYTerminal
    end
    colorize(use_color::Bool, c::Char) = Cthulhu.stringify() do io
        use_color ? printstyled(io, c; color=:cyan) : print(io, c)
    end
    # Write a file that we track with Revise. Creating it programmatically allows us to rewrite it with
    # different content
    fn = tempname()
    open(fn, "w") do io
        println(io,
        """
        function simplef(a, b)
            z = a*a
            return z + b
        end
        """)
    end
    includet(fn)

    # Copy the user's current settings and set up the defaults
    CONFIG = deepcopy(Cthulhu.CONFIG)
    config = Cthulhu.CthulhuConfig()
    for fn in fieldnames(Cthulhu.CthulhuConfig)
        setfield!(Cthulhu.CONFIG, fn, getfield(config, fn))
    end

    try
        fake_terminal() do term, in, out
            @async begin
                descend(simplef, Tuple{Float32, Int32}; interruptexc=false, terminal=term)
            end
            lines = cread(out)
            @test occursin("invoke simplef(::Float32,::Int32)::Float32", lines)
            @test occursin(r"Base\.mul_float\(.*, .*\)\u001B\[36m::Float32\u001B\[39m", lines)
            @test occursin('[' * colorize(true, 'o') * "]ptimize", lines)
            @test occursin('[' * colorize(true, 'T') * "]yped", lines)

            @test occursin("\nSelect a call to descend into", lines)   # beginning of the line
            @test occursin('•', lines)
            write(in, 'o')            # switch to unoptimized
            lines = cread(out)
            @test occursin("invoke simplef(::Float32,::Int32)::Float32", lines)
            @test occursin(r"\(z = a \* a\)\u001B\[\d\dm::Float32\u001B\[39m", lines)
            @test occursin('[' * colorize(false, 'o') * "]ptimize", lines)
            @test occursin("\nSelect a call to descend into", lines)   # beginning of the line
            @test occursin("• %1  = *(::Float32,::Float32)::Float32", lines)
            # Call selection
            write(in, keydict[:down])
            write(in, keydict[:enter])
            lines = cread1(out)
            lines = cread(out)
            @test occursin("• %1  = promote(::Float32,::Int32)::Tuple{Float32, Float32}", lines)
            write(in, keydict[:up])
            write(in, keydict[:enter])
            lines = cread1(out)
            lines = cread(out)
            @test occursin("• %1  = *(::Float32,::Float32)::Float32", lines)  # back to where we started
            write(in, 'o')            # back to optimized
            lines = cread(out)
            @test !occursin("Variables", lines)
            @test occursin(r"Base\.mul_float\(.*, .*\)\u001B\[36m::Float32\u001B\[39m", lines)
            write(in, 'i')            # inline cost
            lines = cread(out)
            @test occursin(r"\u001B\[32m \d\u001B\[39m %\d += Base\.mul_float", lines)
            @test occursin('[' * colorize(true, 'i') * "]nlining costs", lines)
            write(in, 'i')
            lines = cread(out)
            write(in, 'o')
            lines = cread(out)
            @test occursin("Variables", lines)
            write(in, 'w')            # unoptimized + warntype
            lines = cread(out)
            @test occursin("Variables", lines)
            @test occursin(r"z.*::Float32", lines)
            @test occursin(r"\nBody.*Float32", lines)
            @test occursin('[' * colorize(true, 'w') * "]arn", lines)
            @test occursin("\nSelect a call to descend into", lines)   # beginning of the line
            @test occursin("• %1  = *(::Float32,::Float32)::Float32", lines)
            # Source view
            write(in, 'S')
            lines = cread(out)
            @test occursin("""
            \n\nfunction simplef(a, b)
                z = a*a
                return z + b
            end
            """, lines)
            @test occursin('[' * colorize(true, 'S') * "]ource", lines)
            # turn on syntax highlighting
            write(in, 's'); cread(out)
            write(in, 'S')
            lines = cread(out)
            @test occursin("simplef", lines)
            @test occursin("\u001B", first(split(lines, '\n')))
            @test occursin('[' * colorize(true, 's') * "]yntax", lines)
            write(in, 's'); cread(out)  # off again
            # Toggling 'o' goes back to typed code, make sure it also updates the selector status
            write(in, 'o')
            lines = cread(out)
            @test occursin('[' * colorize(true, 'T') * "]yped", lines)
            write(in, 'o'); cread(out)   # toggle it back for later tests
            # AST view
            write(in, 'A')
            lines = cread(out)
            @test occursin("Symbol simplef", lines)
            @test occursin("Symbol call", lines)
            @test occursin('[' * colorize(true, 'A') * "]ST", lines)
            # LLVM view
            write(in, 'L')
            lines = cread(out)
            @test occursin(r"sitofp i(64|32)", lines)
            @test occursin("fadd float", lines)
            @test occursin("┌ @ promotion.jl", lines)  # debug info on by default
            @test occursin('[' * colorize(true, 'L') * "]LVM", lines)
            # turn off debug info
            write(in, 'd'); cread(out)
            write(in, 'd'); cread(out)
            write(in, 'L')
            lines = cread(out)
            @test occursin('[' * colorize(false, 'd') * "]ebuginfo", lines)
            @test !occursin("┌ @ promotion.jl", lines)
            # Native-code view
            write(in, 'N')
            lines = cread(out)
            @test occursin(".text", lines) || occursin("__text", lines)
            @test occursin("retq", lines)
            @test occursin('[' * colorize(true, 'N') * "]ative", lines)
            # Typed-view (by selector)
            write(in, 'T')
            lines = cread(out)
            @test occursin(r"\(z \+ b\)\u001B\[\d\dm::Float32\u001B\[39m", lines)
            @test occursin('[' * colorize(true, 'T') * "]yped", lines)
            # Revise
            open(fn, "w") do io
                println(io,
                """
                function simplef(a, b)
                    z = a*b
                    return z + b
                end
                """)
            end
            sleep(0.1)
            write(in, 'R')
            lines = cread(out)
            write(in, 'T')
            lines = cread(out)
            @test occursin("z = a * b", lines)
            write(in, 'q')
        end
        # Multicall & iswarn=true
        fake_terminal() do term, in, out
            @async begin
                descend_code_warntype(MultiCall.callfmulti, Tuple{Any}; interruptexc=false, optimize=false, terminal=term)
            end
            lines = cread(out)
            @test occursin("\nBody\e[91m\e[1m::Union{Float32, $Int}\e[22m\e[39m", lines)
            @test occursin("Base.getindex(c)\e[91m\e[1m::Any\e[22m\e[39m", lines)
            @test occursin("\e[31m%\e[39m2  = call → fmulti(::Any)::Union{Float32, Int64}", lines)
            write(in, keydict[:enter])
            lines = cread1(out)
            @test occursin("%2  = fmulti(::Int32)::Union{Float32, $Int}", lines)
            @test occursin("%2  = fmulti(::Float32)::Union{Float32, $Int}", lines)
            @test occursin("%2  = fmulti(::Char)::Union{Float32, $Int}", lines)
            write(in, 'q')
        end
        # Tasks (see the special handling in `_descend`)
        ftask() = @sync @async show(io, "Hello")
        fake_terminal() do term, in, out
            @async begin
                @descend terminal=term ftask()
            end
            lines = cread(out)
            @test occursin(r"• %\d\d  = task", lines)
            write(in, keydict[:enter])
            lines = cread(out)
            @test occursin("call show(::IO,::String)", lines)
            write(in, 'q')
        end
        # descend with MethodInstances
        mi = Cthulhu.get_specialization(MultiCall.callfmulti, Tuple{typeof(Ref{Any}(1))})
        fake_terminal() do term, in, out
            @async begin
                descend(mi; interruptexc=false, optimize=false, terminal=term)
            end
            lines = cread(out)
            @test occursin("fmulti(::Any)", lines)
            write(in, 'q')
        end
        fake_terminal() do term, in, out
            @async begin
                descend_code_warntype(mi; interruptexc=false, optimize=false, terminal=term)
            end
            lines = cread(out)
            @test occursin("Base.getindex(c)\e[91m\e[1m::Any\e[22m\e[39m", lines)
            write(in, 'q')
        end

        # ascend
        @noinline inner3(x) = 3x
        @inline   inner2(x) = 2*inner3(x)
                  inner1(x) = -1*inner2(x)
        inner1(0x0123)
        mi = Cthulhu.get_specialization(inner3, Tuple{UInt16})
        fake_terminal() do term, in, out
            @async begin
                ascend(term, mi)
            end
            write(in, keydict[:down])
            write(in, keydict[:enter])
            write(in, keydict[:down])
            write(in, keydict[:enter])
            lines = split(cread(out), '\n')
            idx = first(findfirst("inner3", lines[2]))
            @test first(findfirst("inner2", lines[3])) == idx + 2
            ln = occursin(r"caller.*inner3", lines[6]) ? 6 :
                 occursin(r"caller.*inner3", lines[8]) ? 8 : error("not found")
            @test occursin("inner2", lines[ln+1])
            @test any(str -> occursin("Variables", str), lines[ln+2:end])
            write(in, 'q')
            write(in, 'q')
        end
    finally
        # Restore the previous settings
        for fn in fieldnames(Cthulhu.CthulhuConfig)
            setfield!(Cthulhu.CONFIG, fn, getfield(CONFIG, fn))
        end
        rm(fn)
    end
end


text:
include("Packages.jl") # Import the packages
include("Args.jl") # Parsing the command-line

##################################################
# General parameters
##################################################
"""
    PARALLEL

Determining if the code is run in parallel.
Is 'yes' by default.
"""
const PARALLEL = parsed_args["parallel"]
if ((PARALLEL != "yes") && (PARALLEL != "no"))
    error("ERROR: UNKNOWN PARALLEL") # Unknown parallel procedure
end

"""
    eps

Temperature of the disk, defined as the ratio of the internal
energy over the absolut total energy.

A low temperature means that pressure predominates, while a high
temperature means that self-gravity predominates.

Is '0.1' by default.
"""
const eps = parsed_args["eps"]

"""
    N

Number of basis elements in the linear expansion.

Is '100' by default.
"""
const N = parsed_args["N"]

"""
    nbK

Sampling number for numerical integration.

Is '5000' by default.
"""
const nbK = parsed_args["nbK"]

"""
    m

Azimutal wavenumber in the Fourier expansion in angles.

Is '2' by default (mode with fastest growth rates).
"""
const m = parsed_args["m"]

##################################################

include("Constants.jl") # Physical constants
include("Mean.jl") # Useful functions of the system
include("TableValues.jl") # Compute the matrix elements
include("ResponseMatrix.jl") # Compute the linear response matrix
include("Eigenvalues.jl") # Computes the system eigenvalues
include("GrowthRate.jl") # Computes the growth rate of the system

max_stars_repo_name: KerwannTEP/LiRA

max_issues_repo_name: KerwannTEP/LiRA

max_forks_repo_name: KerwannTEP/LiRA


text:
# This file holds functions that compute the upwinding scheme
# for the value function, given inputs for the savings over
# the state space, and solves the HJB equation,
# given the standard form: ρV = u + AV.

# Upwind step for value function
# Exact upwind, assumes If, Ib, and I0 do not conflict
function upwind_value_function(dVf::Matrix{T}, dVb::Matrix{T}, dV0::Matrix{T},
                               sf::Matrix{T}, sb::Matrix{T}) where {T<:Number}
    If = sf .> 0 # positive drift -> forward diff
    Ib = sb .< 0 # negative drift -> backward dift
    I0 = (1 .- If .- Ib)
    V  =  dVf .* If + dVb .* Ib + dV0 .* I0 # Compute upwinded value function difference matrix
    return V, If, Ib, I0
end

# Inexact upwind for vhen we know exact form of value function (e.g. CRRA w/ no labor disutility)
# and for when we are unsure if If, Ib, and I0 will not uniquely identify the upwind scheme.
# This function has not been verified to work yet since we have not used a model with it yet.
# This was based off the following upwind_value_function from OneAssetHANK
function upwind_value_function(dVf::Matrix{T}, dVb::Matrix{T}, dV0::Matrix{T},
                               cf::Matrix{T}, cb::Matrix{T}, c0::Matrix{T},
                               sf::Matrix{T}, sb::Matrix{T}; reverse_sign::Float64 = -1e12) where {T<:Number}
    Vf = (cf .> 0) .* (sf .* dVf) + (cf .<= 0) * (reverse_sign)
    Vb = (cb .> 0) .* (sb .* dVb) + (cb .<= 0) * (reverse_sign)
    V0 = (c0 .<= 0) * (reverse_sign)

    I_neither = (1 - (sf .> 0)) .* (1 - (sb .< 0)) # exactly zero drift
    I_unique = (sb .< 0) .* (1 - (sf .> 0)) + (1 - (sb .< 0)) .* (sf .> 0) # just one direction
    I_both = (sb .< 0) .* (sf .> 0) # both positive and negative

    Ib = I_unique .* (sb .< 0) .* (Vb .> V0) + I_both .* (Vb .== max.(max.(Vb, Vf), V0))
    If = I_unique .* (sf .> 0) .* (Vf .> V0) + I_both .* (Vf .== max.(max.(Vb, Vf), V0))
    I0 = I_neither + (1 - I_neither) .* (V0 .== max.(max.(Vb, Vf), V0))
    I0 = 1 - Ib - If
    V = dVf .* If + dVb .* Ib + dV0 .* I0

    return V, If, Ib, I0
end

# Inexact upwind, assumes positive consumption drift corresponds to positive value function drift,
# written based on model using CRRA utility with labor disutility.
# The upwind is inexact b/c it accomodates for when the savings difference indicates
# we should use both a forward and a backward difference at the same location in the upwind scheme.
# Vaf - value function without labor disutility applied to forward consumption diff
# Vf - value function w/labor disutility, etc., applied to forward consumption diff, forward hours diff, etc.
function upwind_value_function(Vaf::Matrix{T}, Vab::Matrix{T}, Va0::Matrix{T},
                               Vf_orig::Matrix{T}, Vb_orig::Matrix{T}, V0_orig::Matrix{T},
                               cf::Matrix{T}, cb::Matrix{T}, c0::Matrix{T},
                               sf::Matrix{T}, sb::Matrix{T}; reverse_sign::Float64 = -1e12) where {T<:Number}
    # Approximate forward difference in value function? Not entirely sure.
    # The Vf_orig + sf .* Vaf seems to be an approximation of a differencing of the
    # flow utility + ∂_a V * da_t
    # term in the HJB equation. See the determining of indices for another indication
    # why this seems to be some kind of approximation method.
    Vf = (cf .> 0) .* (Vf_orig + sf .* Vaf) + (cf .<= 0) * (reverse_sign)
    Vb = (cb .> 0) .* (Vb_orig + sb .* Vab) + (cb .<= 0) * (reverse_sign)
    V0 = (c0 .> 0) .* V0_orig + (c0 .<= 0) * (reverse_sign)

    I_neither = (1 - (sf .> 0)) .* (1 - (sb .< 0)) # exactly zero drift at this point
    I_unique = (sb .< 0) .* (1 - (sf .> 0)) + (1 - (sb .< 0)) .* (sf .> 0) # just one direction
    I_both = (sb .< 0) .* (sf .> 0) # both positive and negative

    # Determine indices with forward, backward, and no drift
    # Observe that Vb .> V0 is how we are determining when to consider I_unique or not
    Ib = I_unique .* (sb .< 0) .* (Vb .> V0) + I_both .* (Vb .== max.(max.(Vb, Vf), V0))
    If = I_unique .* (sf .> 0) .* (Vf .> V0) + I_both .* (Vf .== max.(max.(Vb, Vf), V0))
    I0 = I_neither + (1 - I_neither) .* (V0 .== max.(max.(Vb, Vf), V0))
    I0 = 1 - Ib - If

    # Compute upwinded value function difference (in wealth direction) matrix
    V = Vaf .* If + Vab .* Ib + Va0 .* I0

    return V, If, Ib, I0
end

# Constructs A matrix that applies to the HJB
# A_switch - summarizes all information due to non-wealth state variables, e.g. income shocks
function upwind_matrix(A_switch::SparseMatrixCSC{R, Int64},
                       sf::Matrix{S}, sb::Matrix{S},
                       f_diff::Matrix{T}, b_diff::Matrix{T},
                       wealth_dim::Int, other_dims::Int;
                       exact::Bool = true,
                       If::Matrix{Int64} = Matrix{Int64}(0, 0),
                       Ib::Matrix{Int64} = Matrix{Int64}(0, 0),
                       I0::Matrix{Int64} = Matrix{Int64}(0, 0)) where {R<:Number, S<:Number, T<:Number}
    if size(sf, 1) != wealth_dim
        error("Dimension of wealth grid, wealth_dim, is incorrect.")
    elseif size(f_diff) != size(sf) || size(b_diff) != size(sf)
        error("Dimension of diff, differenced state space grid, along the wealth dimension is wrong.")
    end

    # Compute diagonals of first-order differenced V matrix
    if exact
        # use min/max b/c assumes we set I0 entries to 1 only if sf/sb exactly zero
        X = -min.(sb, 0) ./ b_diff # If savings drift < 0, re-weight by grid partition
                                   # (to approx a derivative), and set to negative to get the
                                   # right direction b/c lower diagonal (doing a backward
                                   # difference here so set positive)
        Y = -max.(sf, 0) ./ f_diff + min.(sb, 0) ./ b_diff # Same idea here
        Z = max.(sf, 0) ./ f_diff # Same idea here
    else
        X = -Ib .* sb ./ b_diff
        Y = -If .* sf ./ f_diff + Ib .* sb ./ b_diff
        Z =  If .* sf ./ f_diff
    end

    total_dims = wealth_dim * other_dims
    centdiag = vec(Y) # Y is our main diagonal
    updiag = zeros(R, total_dims - 1)
    lowdiag = zeros(R, total_dims - 1)
    for j = 1:other_dims
        for i = 1:wealth_dim
            # When i == wealth_dim, leave entry as zero b/c matrix should be block
            # [nonzero zeros; zeros nonzero] in the case of a wealth_dim x 2 state space
            # Also, this way, we don't save Z[end, end]
            # or X[1, 1] since X and Z are lower/upper diagonals
            if i < wealth_dim
                updiag[wealth_dim * (j-1) + i] = Z[i, j]
                lowdiag[wealth_dim * (j-1) + i] = X[i + 1, j]
            end
        end
    end

    AA = spdiagm((lowdiag, centdiag, updiag), (-1, 0, 1)) # AA * v ≈ ∂_a v * (time derivative of wealth)
    A = AA + A_switch # main A matrix in standard representation of HACT steady state
    return A
end

# Same function as above but with a constant difference in the state space grid
function upwind_matrix(A_switch::SparseMatrixCSC{R, Int64},
                       sf::Matrix{S}, sb::Matrix{S},
                       diff::T, wealth_dim::Int, other_dims::Int;
                       exact::Bool = true,
                       If::Matrix{Int64} = Matrix{Int64}(undef, 0, 0),
                       Ib::Matrix{Int64} = Matrix{Int64}(undef, 0, 0),
                       I0::Matrix{Int64} = Matrix{Int64}(undef, 0, 0)) where {R<:Number, S<:Number, T<:Number}
    if size(sf, 1) != wealth_dim
        error("Dimension of wealth grid is incorrect.")
    end

    # Compute diagonals of first-order differenced V matrix
    if exact
        # use min/max b/c assumes we set I0 entries to 1 only if sf/sb exactly zero
        X = -min.(sb, 0) / diff
        Y = -max.(sf, 0) / diff + min.(sb, 0) / diff
        Z =  max.(sf, 0) / diff
    else
        # Need to multiply with Ib, If to be consistent with how we set the entries of I0
        X = -Ib .* sb ./ diff
        Y = -If .* sf ./ diff + Ib .* sb ./ diff
        Z =  If .* sf ./ diff
    end
    total_dims = wealth_dim * other_dims
    centdiag = vec(Y)
    updiag = zeros(R, total_dims - 1)
    lowdiag = zeros(R, total_dims - 1)
    for j = 1:other_dims
        for i = 1:wealth_dim
            if i < wealth_dim
                updiag[wealth_dim * (j-1) + i] = Z[i, j]
                lowdiag[wealth_dim * (j-1) + i] = X[i + 1, j]
            end
        end
    end

    AA = SparseArrays.spdiagm(0 => centdiag, 1 => updiag, -1 => lowdiag)
    A = AA + A_switch
    return A
end

# Solves hjb via left divide
# Δ_HJB is used to make the left-divide more numerically stable. Derivation is as follows:
# (ρ - A) * V = u
# (ρ - A) * V + V/Δ_HJB = u + V/Δ_HJB
# Since Δ_HJB is large, V/Δ_HJB ≈ 0 but also adds a small perturbation to reflect
# that we are approximating V, so our vector u may not be exactly correct. Re-arranging yields
# the form ((1/Δ_HJB + ρ) - A) V = u + V/Δ_HJB
function solve_hjb(A::SparseMatrixCSC{R, Int64}, ρ::S,
                   Δ_HJB::S, u::Matrix{T}, V::Matrix{T}) where {R<:Number, S<:Number, T<:Number}
    # Stack u and V
    return solve_hjb(A, ρ, Δ_HJB, vec(u), vec(V))
end

function solve_hjb(A::SparseMatrixCSC{R, Int64}, ρ::S,
                   Δ_HJB::S, u::Vector{T}, V::Vector{T}) where {R<:Number, S<:Number, T<:Number}
    B = (1/Δ_HJB + ρ) * SparseMatrixCSC(I, size(A, 1), size(A, 1)) - A
    b = u + V / Δ_HJB
    return vec(B \ b)
end

# Wrapper function combining upwind_value_function and upwind_matrix
# since these do not require intermediate steps. However, b/c
# solve_hjb needs flow utility, which is model-specific,
# we do not wrap all of these steps into one function
function upwind(A_switch::SparseMatrixCSC{R, Int64},
                Vaf::Matrix{S}, Vab::Matrix{S}, Va0::Matrix{S},
                Vf::Matrix{S}, Vb::Matrix{S}, V0::Matrix{S},
                cf::Matrix{S}, cb::Matrix{S}, c0::Matrix{S},
                sf::Matrix{S}, sb::Matrix{S},
                f_diff::Matrix{T}, b_diff::Matrix{T},
                wealth_dim::Int, other_dims::Int;
                reverse_sign = -1e12, exact::Bool = true) where {R<:Number, S<:Number, T<:Number}
    if exact
        dV_upwind, If, Ib, I0 = upwind_value_function(Vaf, Vab, Va0, sf, sb)
        A = upwind_matrix(A_switch, sf, sb, f_diff, b_diff, wealth_dim, other_dims;
                      exact = exact)
    else
        dV_upwind, If, Ib, I0 = upwind_value_function(Vaf, Vab, Va0, Vf, Vb, V0, cf, cb, c0,
                                                  sf, sb; reverse_sign = reverse_sign)
        A = upwind_matrix(A_switch, sf, sb, f_diff, b_diff, wealth_dim, other_dims;
                      exact = exact, If = If, Ib = Ib, I0 = I0)
    end
    return dV_upwind, If, Ib, I0, A
end

function upwind(A_switch::SparseMatrixCSC{R, Int64},
                Vaf::Matrix{S}, Vab::Matrix{S}, Va0::Matrix{S},
                sf::Matrix{S}, sb::Matrix{S},
                f_diff::Matrix{T}, b_diff::Matrix{T},
                wealth_dim::Int, other_dims::Int,
                cf::Matrix{S} = Matrix{S}(0, 0),
                cb::Matrix{S} = Matrix{S}(0, 0),
                c0::Matrix{S} = Matrix{S}(0, 0);
                exact::Bool = true,
                reverse_sign::Float64 = -1e12) where {R<:Number, S<:Number, T<:Number}

    if exact
        dV_upwind, If, Ib, I0 = upwind_value_function(Vaf, Vab, Va0, sf, sb)
        A = upwind_matrix(A_switch, sf, sb, f_diff, b_diff, wealth_dim, other_dims;
                      exact = exact)
    else
        dV_upwind, If, Ib, I0 = upwind_value_function(Vaf, Vab, Va0, cf, cb, c0, sf, sb;
                                                      reverse_sign = reverse_sign)
        A = upwind_matrix(A_switch, sf, sb, f_diff, b_diff, wealth_dim, other_dims;
                      exact = exact, If = If, Ib = Ib, I0 = I0)
    end

    return dV_upwind, If, Ib, I0, A
end

function upwind(A_switch::SparseMatrixCSC{R, Int64},
                Vaf::Matrix{S}, Vab::Matrix{S}, Va0::Matrix{S},
                Vf::Matrix{S}, Vb::Matrix{S}, V0::Matrix{S},
                cf::Matrix{S}, cb::Matrix{S}, c0::Matrix{S},
                sf::Matrix{S}, sb::Matrix{S},
                diff::T, wealth_dim::Int, other_dims::Int;
                reverse_sign = -1e12, exact::Bool = trie) where {R<:Number, S<:Number, T<:Number}
    if exact
        dV_upwind, If, Ib, I0 = upwind_value_function(Vaf, Vab, Va0, sf, sb)
        A = upwind_matrix(A_switch, sf, sb, diff, wealth_dim, other_dims;
                      exact = exact)
    else
        dV_upwind, If, Ib, I0 = upwind_value_function(Vaf, Vab, Va0, Vf, Vb, V0, cf, cb, c0,
                                                  sf, sb; reverse_sign = reverse_sign)
        A = upwind_matrix(A_switch, sf, sb, diff, wealth_dim, other_dims;
                          exact = exact, If = If, Ib = Ib, I0 = I0)
    end
    return dV_upwind, If, Ib, I0, A
end

function upwind(A_switch::SparseMatrixCSC{R, Int64},
                Vaf::Matrix{S}, Vab::Matrix{S}, Va0::Matrix{S},
                sf::Matrix{S}, sb::Matrix{S},
                diff::T, wealth_dim::Int, other_dims::Int;
                exact::Bool = true,
                reverse_sign::Float64 = -1e12) where {R<:Number, S<:Number, T<:Number}
    if exact
        dV_upwind, If, Ib, I0 = upwind_value_function(Vaf, Vab, Va0, sf, sb)
        A = upwind_matrix(A_switch, sf, sb, diff, wealth_dim, other_dims;
                      exact = exact)
    else
        dV_upwind, If, Ib, I0 = upwind_value_function(Vaf, Vab, Va0, cf, cb, c0, sf, sb;
                                                      reverse_sign = reverse_sign)
        A = upwind_matrix(A_switch, sf, sb, diff, wealth_dim, other_dims;
                      exact = exact, If = If, Ib = Ib, I0 = I0)
    end

    return dV_upwind, If, Ib, I0, A
end


text:
RpcApi.SimpleOrganism[RpcApi.SimpleOrganism(0x00000000000021c6,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_1::Int8 = 39
        var_2 = cot(if var_1 === (Int8)(0)
                    (Int8)(1)
                else
                    var_1
                end)
        Creature.stepDown(o)
        var_1 = if 127 >= abs(7059498975635099421)
                Int8(7059498975635099421)
            else
                127
            end
        var_2 = cot(if var_2 === (Float64)(0)
                    (Float64)(1)
                else
                    var_2
                end)
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        return true
    end),[97,53,74,3,69,88,38,55,7],43.3261,886,15.8708,125426,0x0000006a,Dict{Int16,Int16}(),[589,701],28907,4,7),RpcApi.SimpleOrganism(0x0000000000002448,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_1::Int8 = 39
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.stepLeft(o)
        Creature.stepLeft(o)
        Creature.eatUp(c,o,Int(var_1))
        return true
    end),[26,100,9,32,67,86,24,1,24],57.9721,412,58.2612,69743,0x00000066,Dict{Int16,Int16}(),[135,316],22481,64,6),RpcApi.SimpleOrganism(0x0000000000002657,:(function (c::Config.ConfigData,o::Creature.Organism)
        function func_5(var_2::Bool=false,var_3::Int16=-15437,var_4::String="wasUVB6e")
            var_3 = var_3 + var_3
            return var_2
        end
        local var_1::Int8 = 39
        Creature.eatLeft(c,o,Int(var_1))
        var_1 = if 127 >= abs(var_1)
                Int8(var_1)
            else
                127
            end
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        var_1 = var_1 - var_1
        Creature.stepLeft(o)
        return true
    end),[26,100,1,32,7,11,24,1,24],85.7587,863,1.55414,59454,0x00000067,Dict{Int16,Int16}(),[666,354],17203,63,12),RpcApi.SimpleOrganism(0x0000000000002933,:(function (c::Config.ConfigData,o::Creature.Organism)
        function func_4(var_2::String="VLXNQsjK",var_3::Float64=0.393312)
            local var_5::Float64 = 0.862193
            var_3 = sin(-2808405706071492296)
            var_3 = if isempty(var_2)
                    0.0
                else
                    1.0
                end
            var_3 = var_3 + var_3
            var_3 = tan(17974)
            var_3 = tan(4908374754400172539)
            return var_2
        end
        local var_1::Int8 = 77
        Creature.eatRight(c,o,Int(var_1))
        var_1 = (Int8)(var_1) << (Int8)(var_1)
        var_1 = if 127 >= abs(-1189142824310634600)
                Int8(-1189142824310634600)
            else
                127
            end
        for i::Int8 = 1:div(var_1,16)
        end
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.stepDown(o)
        for i::Int8 = 1:div(var_1,16)
        end
        Creature.eatLeft(c,o,Int(var_1))
        Creature.stepUp(o)
        Creature.stepRight(o)
        Creature.stepRight(o)
        Creature.eatLeft(c,o,Int(var_1))
        var_1 = var_1 % if var_1 === (Int8)(0)
                    (Int8)(1)
                else
                    var_1
                end
        var_1 = if 127 >= abs(-13122)
                Int8(-13122)
            else
                127
            end
        return true
    end),[31,97,1,0,96,47,16,1,1],66.2066,876,16.4457,58599,0x0000006c,Dict{Int16,Int16}(),[473,314],9875,50,24),RpcApi.SimpleOrganism(0x00000000000021a2,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_1::Int8 = 77
        function func_5(var_2::Bool=false,var_3::Bool=true,var_4::Int64=5360845696666748569)
            var_4 = Int(Creature.idLeft(o))
            var_3 = Bool(abs(var_3 - var_3))
            return var_2
        end
        Creature.stepRight(o)
        Creature.eatRight(c,o,Int(var_1))
        var_1 = var_1 - var_1
        return true
    end),[50,97,1,69,1,37,47,1,81],2.0,590,17.188,54194,0x00000069,Dict{Int16,Int16}(),[1604,258],29267,100,7),RpcApi.SimpleOrganism(0x00000000000022ad,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_2::Bool = false
        local var_1::Int8 = 39
        var_1 = if 127 >= abs(-9842)
                Int8(-9842)
            else
                127
            end
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.stepUp(o)
        var_2 = Creature.stepDown(o)
        for i::Int8 = 1:div(var_1,16)
            Creature.stepLeft(o)
        end
        Creature.stepLeft(o)
        Creature.stepLeft(o)
        Creature.eatUp(c,o,Int(var_1))
        var_1 = var_1 - var_1
        var_1 = var_1 - var_1
        return true
    end),[80,91,1,32,92,11,24,1,24],83.2223,625,12.8025,51787,0x00000067,Dict{Int16,Int16}(),[871,277],26599,52,14),RpcApi.SimpleOrganism(0x0000000000002801,:(function (c::Config.ConfigData,o::Creature.Organism)
        function func_4(var_2::String="MY3VTiTE",var_3::Int16=6893)
            Creature.stepDown(o)
            Creature.stepLeft(o)
            var_3 = var_3 * var_3
            return var_2
        end
        local var_1::Int8 = 77
        var_1 = var_1 % if var_1 === (Int8)(0)
                    (Int8)(1)
                else
                    var_1
                end
        Creature.eatDown(c,o,Int(var_1))
        Creature.stepRight(o)
        Creature.eatLeft(c,o,Int(var_1))
        var_1 = if 127 >= abs(3392)
                Int8(3392)
            else
                127
            end
        Creature.stepRight(o)
        Creature.eatRight(c,o,Int(var_1))
        return true
    end),[59,35,41,73,91,89,47,2,1],50.7291,887,14.9913,48154,0x00000067,Dict{Int16,Int16}(),[1700,545],12935,50,12),RpcApi.SimpleOrganism(0x0000000000002852,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_1::Int8 = 77
        Creature.eatDown(c,o,Int(var_1))
        var_1 = var_1 % if var_1 === (Int8)(0)
                    (Int8)(1)
                else
                    var_1
                end
        Creature.eatRight(c,o,Int(var_1))
        Creature.stepRight(o)
        Creature.eatUp(c,o,Int(var_1))
        Creature.stepLeft(o)
        return true
    end),[24,0,42,68,1,70,77,2,48],87.4785,780,64.1164,47922,0x00000071,Dict{Int16,Int16}(),[847,577],12124,100,7),RpcApi.SimpleOrganism(0x0000000000002186,:(function (c::Config.ConfigData,o::Creature.Organism)
        function func_4(var_3::String="pz5KgfwO")
            local var_5::Int16 = -20115
            return var_3
        end
        local var_1::Int8 = 77
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        for i::Int8 = 1:div(var_1,16)
            local var_2::String = "uVf1MD5l"
        end
        for i::Int8 = 1:div(var_1,16)
            Creature.stepRight(o)
        end
        Creature.eatDown(c,o,Int(var_1))
        for i::Int8 = 1:div(var_1,16)
        end
        return true
    end),[50,97,1,99,1,47,47,1,79],96.1186,872,96.745,47674,0x0000006a,Dict{Int16,Int16}(),[1666,195],29548,89,13),RpcApi.SimpleOrganism(0x000000000000201e,:(function (c::Config.ConfigData,o::Creature.Organism)
        function func_5(var_2::Float64=0.490693,var_3::Bool=true,var_4::String="GqrOVNUy")
            var_4 = string(var_4)
            var_2 = csc(if var_2 === (Float64)(0)
                        (Float64)(1)
                    else
                        var_2
                    end)
            var_3 = (Bool)(var_3) $ (Bool)(var_3)
            var_3 = Creature.stepDown(o)
            var_2 = pi
            var_4 = string(-60)
            var_2 = pi
            return var_2
        end
        local var_1::Int8 = 39
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.stepDown(o)
        var_1 = (Int8)(var_1) | (Int8)(var_1)
        Creature.stepLeft(o)
        return true
    end),[26,64,1,91,92,11,24,1,0],92.0603,821,11.1109,45014,0x0000006a,Dict{Int16,Int16}(),[265,121],33154,1,17),RpcApi.SimpleOrganism(0x0000000000002882,:(function (c::Config.ConfigData,o::Creature.Organism)
        function func_4(var_2::Int64=-7850605145822392225,var_3::String="9SwUF9QL")
            Creature.stepLeft(o)
            var_2 = Creature.energyDown(o)
            var_2 = Int(Creature.idUp(o))
            var_2 = Creature.energyDown(o)
            if var_2 === var_2
                Creature.stepUp(o)
            end
            return var_2
        end
        local var_1::Int8 = 39
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.stepDown(o)
        Creature.eatRight(c,o,Int(var_1))
        var_1 = var_1 + var_1
        Creature.stepRight(o)
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        for i::Int8 = 1:div(var_1,16)
        end
        Creature.eatDown(c,o,Int(var_1))
        return true
    end),[26,100,1,32,10,11,24,1,24],10.5548,973,13.6002,44008,0x00000066,Dict{Int16,Int16}(),[45,17],11644,10,18),RpcApi.SimpleOrganism(0x0000000000002cc4,:(function (c::Config.ConfigData,o::Creature.Organism)
        function func_4(var_2::String="MY3VTiTE",var_3::Int16=6893)
            return var_2
        end
        local var_1::Int8 = 77
        Creature.eatDown(c,o,Int(var_1))
        if var_1 == var_1
        end
        Creature.stepDown(o)
        if var_1 === var_1
        end
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        return true
    end),[24,0,91,68,61,70,77,2,48],40.1163,942,94.9506,43023,0x00000070,Dict{Int16,Int16}(),[1056,198],740,50,8),RpcApi.SimpleOrganism(0x0000000000002af0,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_2::Float64 = 0.94503
        local var_1::Int8 = 39
        Creature.eatLeft(c,o,Int(var_1))
        var_2 = cot(if -30854 === (Int16)(0)
                    (Int16)(1)
                else
                    -30854
                end)
        Creature.stepLeft(o)
        Creature.eatUp(c,o,Int(var_1))
        return true
    end),[26,100,9,32,67,86,24,1,24],81.1283,374,10.3787,41322,0x00000066,Dict{Int16,Int16}(),[77,315],5419,84,6),RpcApi.SimpleOrganism(0x0000000000001e26,:(function (c::Config.ConfigData,o::Creature.Organism)
        function func_4(var_2::Int64=-8351618089127045604,var_3::Float64=0.81045)
            var_3 = pi
            Creature.stepDown(o)
            return var_2
        end
        local var_1::Int8 = 39
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.stepUp(o)
        return true
    end),[26,35,1,83,76,11,24,1,24],60.9549,887,17.4913,41228,0x00000066,Dict{Int16,Int16}(),[380,938],38192,5,8),RpcApi.SimpleOrganism(0x0000000000002901,:(function (c::Config.ConfigData,o::Creature.Organism)
        function func_4(var_2::String="MY3VTiTE",var_3::Int16=6893)
            if var_2 === var_2
            end
            var_3 = (Int16)(var_3) << (Int16)(var_3)
            var_3 = if 32767 >= abs(-2449718953369934346)
                    Int16(-2449718953369934346)
                else
                    32767
                end
            var_2 = string(4758249667505520108)
            Creature.stepUp(o)
            Creature.stepLeft(o)
            return var_2
        end
        local var_1::Int8 = 77
        Creature.eatDown(c,o,Int(var_1))
        var_1 = if 127 >= abs(3016)
                Int8(3016)
            else
                127
            end
        Creature.stepDown(o)
        for i::Int8 = 1:div(var_1,16)
        end
        Creature.stepDown(o)
        var_1 = (Int8)(var_1) << (Int8)(var_1)
        Creature.eatRight(c,o,Int(var_1))
        return true
    end),[78,35,1,97,11,81,47,1,99],14.5079,768,14.5033,41180,0x0000006b,Dict{Int16,Int16}(),[490,739],10374,65,15),RpcApi.SimpleOrganism(0x0000000000002a75,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_3::Int8 = 2
        local var_1::Int8 = 77
        local var_2::Int16 = -22952
        var_2 = Int16(Creature.eatLeft(c,o,Int(var_3)))
        for i::Int8 = 1:div(var_3,16)
        end
        var_2 = if true
                1
            else
                0
            end
        for i::Int8 = 1:div(var_1,16)
        end
        var_2 = Int16(o.pos.x)
        var_2 = Int16(Creature.eatDown(c,o,Int(var_1)))
        var_2 = Int16(o.pos.y)
        Creature.eatUp(c,o,Int(var_1))
        var_1 = if 127 >= abs(var_2)
                Int8(var_2)
            else
                127
            end
        var_2 = Int16(Creature.eatDown(c,o,Int(var_1)))
        Creature.stepDown(o)
        return true
    end),[50,97,1,89,56,47,47,1,1],75.4864,951,22.1102,40367,0x00000068,Dict(-22952=>-22952),[509,822],6652,50,14),RpcApi.SimpleOrganism(0x0000000000002754,:(function (c::Config.ConfigData,o::Creature.Organism)
        function func_6(var_5::String="FsGapKff")
            return var_5
        end
        local var_1::Int8 = 39
        Creature.stepLeft(o)
        Creature.stepUp(o)
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        var_1 = (Int8)(var_1) | (Int8)(var_1)
        var_1 = if isempty("Mnyn9tG3")
                0
            else
                1
            end
        Creature.stepUp(o)
        Creature.stepRight(o)
        return true
    end),[29,11,82,85,76,98,47,45,5],10.8698,868,6.65285,40254,0x0000006a,Dict{Int16,Int16}(),[625,637],14672,46,10),RpcApi.SimpleOrganism(0x00000000000027d8,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_1::Int8 = 77
        for i::Int8 = 1:div(var_1,16)
        end
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        var_1 = (Int8)(var_1) >> (Int8)(var_1)
        Creature.stepRight(o)
        Creature.stepUp(o)
        return true
    end),[72,48,11,69,81,22,38,10,1],25.3125,711,7.83146,40231,0x00000065,Dict{Int16,Int16}(),[536,46],13345,50,10),RpcApi.SimpleOrganism(0x0000000000002173,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_1::Int8 = 77
        function func_5(var_2::String="hnEOMRsP",var_3::Bool=true,var_4::Bool=false)
            return var_2
        end
        for i::Int8 = 1:div(var_1,16)
        end
        for i::Int8 = 1:div(var_1,16)
        end
        Creature.stepDown(o)
        Creature.eatDown(c,o,Int(var_1))
        return true
    end),[50,97,16,0,1,47,47,1,1],53.9117,930,46.6299,40036,0x00000065,Dict{Int16,Int16}(),[472,482],29744,2,6),RpcApi.SimpleOrganism(0x00000000000029c0,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_1::Int8 = 77
        for i::Int8 = 1:div(var_1,16)
            local var_2::Int64 = -8738709802300115691
        end
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        var_1 = (Int8)(var_1) $ (Int8)(var_1)
        if var_1 == var_1
        end
        var_1 = if 127 >= abs(-31941)
                Int8(-31941)
            else
                127
            end
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.stepUp(o)
        return true
    end),[8,97,1,11,49,85,47,1,81],45.8762,937,10.0402,39910,0x00000069,Dict{Int16,Int16}(),[360,903],8460,92,11)]

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
@testset "dummy conshdlr (always satisfied, no constraint)" begin
    # create an empty problem
    mscip = SCIP.ManagedSCIP()
    SCIP.set_parameter(mscip, "display/verblevel", 0)

    # add the constraint handler
    ch = Dummy.DummyConsHdlr()
    SCIP.include_conshdlr(mscip, ch; needs_constraints=false)

    # solve the problem
    SCIP.@SC SCIP.SCIPsolve(mscip.scip[])

    @test ch.check_called >= 1
    @test ch.enfo_called == 0
    @test ch.lock_called == 1

    # free the problem
    finalize(mscip)
end

@testset "dummy conshdlr (always satisfied, with constraint)" begin
    # create an empty problem
    mscip = SCIP.ManagedSCIP()
    SCIP.set_parameter(mscip, "display/verblevel", 0)

    # add the constraint handler
    ch = Dummy.DummyConsHdlr()
    SCIP.include_conshdlr(mscip, ch; needs_constraints=true)

    # add dummy constraint
    cr = SCIP.add_constraint(mscip, ch, Dummy.DummyCons())

    # solve the problem
    SCIP.@SC SCIP.SCIPsolve(mscip.scip[])

    @test ch.check_called >= 1
    @test ch.enfo_called == 0
    @test ch.lock_called == 1

    # free the problem
    finalize(mscip)
end

@testset "dummy conshdlr (always satisfied, no constraint, but needs it)" begin
    # create an empty problem
    mscip = SCIP.ManagedSCIP()
    SCIP.set_parameter(mscip, "display/verblevel", 0)

    # add the constraint handler
    ch = Dummy.DummyConsHdlr()
    SCIP.include_conshdlr(mscip, ch; needs_constraints=true)

    # solve the problem
    SCIP.@SC SCIP.SCIPsolve(mscip.scip[])

    @test ch.check_called == 0
    @test ch.enfo_called == 0
    @test ch.lock_called == 0

    # free the problem
    finalize(mscip)
end

@testset "never satisfied conshdlr (does not need constraint)" begin
    # create an empty problem
    mscip = SCIP.ManagedSCIP()
    SCIP.set_parameter(mscip, "display/verblevel", 0)

    # add the constraint handler
    ch = NeverSatisfied.NSCH()
    SCIP.include_conshdlr(mscip, ch; needs_constraints=false)

    # solve the problem
    SCIP.@SC SCIP.SCIPsolve(mscip.scip[])

    @test ch.check_called >= 1
    @test ch.enfo_called == 1
    @test ch.lock_called == 1

    # free the problem
    finalize(mscip)
end

@testset "never satisfied conshdlr (needs constraint but does not have it)" begin
    # create an empty problem
    mscip = SCIP.ManagedSCIP()
    SCIP.set_parameter(mscip, "display/verblevel", 0)

    # add the constraint handler
    ch = NeverSatisfied.NSCH()
    SCIP.include_conshdlr(mscip, ch; needs_constraints=true)

    # solve the problem
    SCIP.@SC SCIP.SCIPsolve(mscip.scip[])

    @test ch.check_called == 0
    @test ch.enfo_called == 0
    @test ch.lock_called == 0

    # free the problem
    finalize(mscip)
end

@testset "never satisfied conshdlr (needs constraint and has one)" begin
    # create an empty problem
    mscip = SCIP.ManagedSCIP()
    SCIP.set_parameter(mscip, "display/verblevel", 0)

    # add the constraint handler
    ch = NeverSatisfied.NSCH()
    SCIP.include_conshdlr(mscip, ch; needs_constraints=true)

    # add one constraint
    cr = SCIP.add_constraint(mscip, ch, NeverSatisfied.Cons())

    # solve the problem
    SCIP.@SC SCIP.SCIPsolve(mscip.scip[])

    @test ch.check_called >= 1
    @test ch.enfo_called == 1
    @test ch.lock_called == 1

    # free the problem
    finalize(mscip)
end