text:
using Test
using Mimi

function page_model()
    Mimi.reset_compdefs()

    include("../src/utils/load_parameters.jl")
    include("../src/utils/mctools.jl")

    m = Model()

    set_dimension!(m, :time, [2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200, 2250, 2300])
    set_dimension!(m, :region, ["EU", "USA", "OECD", "USSR", "China", "SEAsia", "Africa", "LatAmerica"])

    return m
end

function get_scenario(ii)
    Mimi.reset_compdefs()
    if ii == 1
        return "ndcs", "NDCs", true, false
    end
    if ii == 2
        return "2c-saf", "2 degC Target", true, true
    end
end

@testset "mimi-page.jl" begin

##
# Tests not suitable yet for annual version with its own parameters
##
# include("test_climatemodel.jl")
include("test_AbatementCosts.jl")
include("test_AdaptationCosts.jl")
include("test_CH4cycle.jl")
include("test_CH4emissions.jl")
include("test_CH4forcing.jl")
# include("test_ClimateTemperature.jl")
include("test_CO2cycle.jl")
include("test_CO2emissions.jl")
include("test_CO2forcing.jl")
# include("test_Discontinuity.jl")
# include("test_EquityWeighting.jl")
include("test_GDP.jl")
include("test_LGcycle.jl")
include("test_LGemissions.jl")
include("test_LGforcing.jl")
include("test_loadparameters.jl")
# include("test_mainmodel.jl")
# include("test_mainmodel_noperm.jl")
# ##include("test_MarketDamages.jl") # missing data
# include("test_MarketDamagesBurke.jl")
include("test_N2Ocycle.jl")
include("test_N2Oemissions.jl")
include("test_N2Oforcing.jl")
# include("test_NonMarketDamages.jl")
include("test_Population.jl")
include("test_SeaLevelRise.jl")
include("test_SLRDamages.jl")
include("test_SulphateForcing.jl")
include("test_TotalAbatementCosts.jl")
include("test_TotalAdaptationCosts.jl")
include("test_TotalForcing.jl")
include("test_Permafrost.jl")
# #include("test_mcs.jl") # missing data
# include("contrib/test_taxeffect.jl")

end

max_stars_repo_name: jkikstra/PAGE-VAR

max_issues_repo_name: jkikstra/PAGE-VAR

max_forks_repo_name: jkikstra/PAGE-VAR


text:
"""
    body_derivatives(system, surfaces, reference, freestream; kwargs...)

Returns the derivatives of the body forces and moments with respect to the
freestream velocity components `(u, v, w)` and the angular velocity components
`(p, q, r)` in the body frame.

The derivatives are returned as two named tuples: `dCF, dCM`

Note that the derivatives with respect to the freestream variables of the panel
forces must have been previously computed and stored in `system`

# Arguments:
 - `system`: Object of type `System` which holds system properties
"""
function body_derivatives(system)

    CF, CM, dCF, dCM = body_forces_derivatives(system)

    fs = system.freestream[]

    # unpack derivatives
    (CF_a, CF_b, CF_p, CF_q, CF_r) = dCF
    (CM_a, CM_b, CM_p, CM_q, CM_r) = dCM

    # convert derivatives to be with respect to (u, v, w) instead of (alpha, beta)
    u, v, w = freestream_velocity(fs)
    u2, v2, w2 = u^2, v^2, w^2
    # a = atan(w, u)
    a_u = -w/(u2+w2)
    a_w = u/(u2+w2)
    # b = atan(v, u)
    b_u =-v/(u2+v2)
    b_v = u/(u2+v2)

    CF_u = CF_a*a_u + CF_b*b_u
    CF_v = CF_b*b_v
    CF_w = CF_a*a_w

    CM_u = CM_a*a_u + CM_b*b_u
    CM_v = CM_b*b_v
    CM_w = CM_a*a_w

    # pack up derivatives as named tuples
    dCF = (u=CF_u, v=CF_v, w=CF_w, p=CF_p, q=CF_q, r=CF_r)
    dCM = (u=CM_u, v=CM_v, w=CM_w, p=CM_p, q=CM_q, r=CM_r)

    return dCF, dCM
end

"""
    stability_derivatives(system)

Returns the derivatives of the body forces and moments in the stability frame
with respect to the freestream velocity components `(alpha, beta)` and the angular
velocity components `(p, q, r)` in the stability frame.

The derivatives are returned as two named tuples: `dCF, dCM`

Note that the derivatives with respect to the freestream variables of the panel
forces must have been previously computed and stored in `system`

# Arguments:
 - `system`: Object of type `System` which holds system properties
"""
function stability_derivatives(system)

    CFb, CMb, dCFb, dCMb = body_forces_derivatives(system)

    fs = system.freestream[]
    ref = system.reference[]

    # unpack derivatives
    (CFb_a, CFb_b, CFb_pb, CFb_qb, CFb_rb) = dCFb
    (CMb_a, CMb_b, CMb_pb, CMb_qb, CMb_rb) = dCMb

    # rotate forces/moments into the stability frame
    R, R_a = body_to_stability_alpha(fs)

    CFs = R * CFb
    CFs_a = R * CFb_a + R_a * CFb
    CFs_b = R * CFb_b
    CFs_pb = R * CFb_pb
    CFs_qb = R * CFb_qb
    CFs_rb = R * CFb_rb

    CMs = R * CMb
    CMs_a = R * CMb_a + R_a * CMb
    CMs_b = R * CMb_b
    CMs_pb = R * CMb_pb
    CMs_qb = R * CMb_qb
    CMs_rb = R * CMb_rb

    # switch (p,q,r) derivatives to be wrt (p,q,r) of stability frame
    Ωb = fs.Omega # body p, q, r
    Ωs = R * Ωb # stability p, q, r

    # Ωb = R' * Ωs
    pb_ps, qb_ps, rb_ps = R[1,:]
    pb_qs, qb_qs, rb_qs = R[2,:]
    pb_rs, qb_rs, rb_rs = R[3,:]

    CFs_ps = CFs_pb*pb_ps + CFs_qb*qb_ps + CFs_rb*rb_ps
    CFs_qs = CFs_pb*pb_qs + CFs_qb*qb_qs + CFs_rb*rb_qs
    CFs_rs = CFs_pb*pb_rs + CFs_qb*qb_rs + CFs_rb*rb_rs

    CMs_ps = CMs_pb*pb_ps + CMs_qb*qb_ps + CMs_rb*rb_ps
    CMs_qs = CMs_pb*pb_qs + CMs_qb*qb_qs + CMs_rb*rb_qs
    CMs_rs = CMs_pb*pb_rs + CMs_qb*qb_rs + CMs_rb*rb_rs

    # assign outputs, and apply stability derivative normalizations
    CF_a = CFs_a
    CF_b = CFs_b
    CF_p = CFs_ps*2*ref.V/ref.b
    CF_q = CFs_qs*2*ref.V/ref.c
    CF_r = CFs_rs*2*ref.V/ref.b

    CM_a = CMs_a
    CM_b = CMs_b
    CM_p = CMs_ps*2*ref.V/ref.b
    CM_q = CMs_qs*2*ref.V/ref.c
    CM_r = CMs_rs*2*ref.V/ref.b

    # pack up derivatives as named tuples
    dCF = (alpha=CF_a, beta=CF_b, p=CF_p, q=CF_q, r=CF_r)
    dCM = (alpha=CM_a, beta=CM_b, p=CM_p, q=CM_q, r=CM_r)

    return dCF, dCM
end

ext: jl


text:
function interaction_effect_probabilities(ch, data)
    γpost_x = collect(hcat(data.Xunique...)')
    γpost_mins = [minimum(γpost_x[:, i]) for i in 1:size(γpost_x, 2)]
    γpost_x = γpost_x - ones(size(γpost_x, 1)) * (γpost_mins .- 1)'

    ch_gamma_mat = collect(hcat(ch.gamma...)')
    nvar = size(γpost_x, 2)
    p_effects = zeros(nvar * (nvar - 1) ÷ 2)
    var1 = zeros(Int, nvar * (nvar - 1) ÷ 2)
    var2 = zeros(Int, nvar * (nvar - 1) ÷ 2)
    row = 1
    for i in 1:(nvar - 1), j = (i + 1):nvar
        subset = deepcopy(γpost_x)
        subset[:, i] = 2 .- (subset[:, i] .!= 1)
        subset[:, j] = 2 .- (subset[:, j] .!= 1)
        subset = maximum(subset, dims = 2)[:, 1] .== 1
        if sum(subset) > 0
            p_effects[row] = mean(maximum(ch_gamma_mat[:, subset], dims = 2))
            var1[row] = i
            var2[row] = j
            row += 1
        end
    end
    return DataFrame(
        var1 = "x" .* string.(var1),
        var2 = "x" .* string.(var2),
        prob = p_effects
    )[1:(row - 1), :]
end


text: add(a, b) = ([a b] * [1,1])[1]

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

const OptionalLADJ = Union{Real,Missing}

"""
    abstract type AbstractVariateTransform <: Function

Abstract type for change-of-variables transformations.

Subtypes (e.g. `SomeTrafo <: AbstractVariateTransform`) must support (with
`trafo::SomeTrafo`):

```julia
    (trafo)(v_prev::SomeVariate) == v_new
    (trafo)(v_prev::SomeVariate, ladj_prev::Union{Real,Missing})) == (v = v_new, ladj = ladj_new)
    (trafo)(s_prev::DensitySample)::DensitySample
    ((trafo2 ∘ trafo1)(v)::AbstractVariateTransform)(v) == trafo2(trafo1(v))
    inv(trafo)(trafo(v)) == v
    inv(inv(trafo)) == trafo

    ValueShapes.varshape(trafo)::ValueShapes.AbstractValueShape
    ValueShapes.valshape(trafo)::ValueShapes.AbstractValueShape
```

with `valshape(v_prev) == varshape(trafo)` and
`valshape(trafo(v_prev)) == valshape(trafo)`

`ladj` must be `logabsdet(jacobian(trafo, v))`.
"""
abstract type AbstractVariateTransform <: Function end
export AbstractVariateTransform

ValueShapes.unshaped(trafo::AbstractVariateTransform) =
    _generic_unshaped_impl(trafo, varshape(trafo), valshape(trafo))

_generic_unshaped_impl(trafo::AbstractVariateTransform, ::ArrayShape{<:Real,1}, ::ArrayShape{<:Real,1}) =
    trafo

# ToDo: Add `UnshapedVariateTransform` that shapes at input and unshaped at
# output, to return from
# `_generic_unshaped_impl(trafo::AbstractVariateTransform, ::AbstractValueShape, ::AbstractValueShape)`.

"""
    ladjof(r::NamedTuple{(...,:ladj,...)})::Union{Real,Missing}

Extract the `log(abs(det(jacobian)))` value that is part of a result `r`.

Examples:

```julia
ladjof((..., ladj = some_ladj, ...)) == some_ladj
ladjof(trafo)(v) = trafo(v, )
```
"""
function ladjof end
export ladjof

ladjof(x::NamedTuple) = x.ladj

struct LADJOfVarTrafo{T<:AbstractVariateTransform} <: Function
    trafo::T
end

(ladjof_trafo::LADJOfVarTrafo)(v::Any) = ladjof(trafo(v, 0))
(ladjof_trafo::LADJOfVarTrafo)(v::Any, prev_ladj::OptionalLADJ) = ladjof(trafo(v, prev_ladj))

"""
    ladjof(trafo::AbstractVariateTransform)::Function

Returns a function that computes the `log(abs(det(jacobian)))` of `trafo` for
a given variate `v`:

```julia
    ladjof(trafo)(v) == ladjof(trafo(v, 0))
    ladjof(trafo)(v, prev_ladj) == ladjof(trafo(v, prev_ladj))
```
"""
ladjof(trafo::AbstractVariateTransform) = LADJOfVarTrafo(trafo)

function _transform_density_sample(trafo::AbstractVariateTransform, s::DensitySample)
    r = trafo(s.v, zero(Float32))
    v = r.v
    logd = s.logd - r.ladj
    DensitySample(v, logd, s.weight, s.info, s.aux)
end

(trafo::AbstractVariateTransform)(s::DensitySample) = _transform_density_sample(trafo, s)

# Custom broadcast(::AbstractVariateTransform, DensitySampleVector), multithreaded:
function Base.copy(
    instance::Base.Broadcast.Broadcasted{
        <:Base.Broadcast.AbstractArrayStyle{1},
        <:Any,
        <:AbstractVariateTransform,
        <:Tuple{<:Union{ArrayOfSimilarVectors{<:Real},ShapedAsNTArray}}
    }
)
    trafo = instance.f
    v_src = instance.args[1]
    vs_trg = valshape(trafo)
    R = eltype(unshaped(trafo(first(v_src)), vs_trg))
    v_src_us = unshaped.(v_src)
    trafo_us = unshaped(trafo)

    n = length(eachindex(v_src_us))
    v_trg_unshaped = nestedview(similar(flatview(v_src_us), R, totalndof(vs_trg), n))
    @assert axes(v_trg_unshaped) == axes(v_src)
    @assert v_trg_unshaped isa ArrayOfSimilarArrays
    @threads for i in eachindex(v_trg_unshaped, v_src)
        v_trg_unshaped[i] = trafo_us(v_src_us[i])
    end
    vs_trg.(v_trg_unshaped)
end

function Base.copy(
    instance::Base.Broadcast.Broadcasted{
        <:Base.Broadcast.AbstractArrayStyle,
        <:Any,
        <:AbstractVariateTransform,
        <:Tuple{DensitySampleVector}
    }
)
    trafo = instance.f
    s_src = instance.args[1]
    vs_trg = valshape(trafo)
    R = eltype(unshaped(trafo(first(s_src.v)), vs_trg))
    s_src_us = unshaped.(s_src)
    trafo_us = unshaped(trafo)

    n = length(eachindex(s_src_us))
    s_trg_unshaped = DensitySampleVector((
        nestedview(similar(flatview(s_src_us.v), R, totalndof(vs_trg), n)),
        zero(s_src_us.logd),
        deepcopy(s_src_us.weight),
        deepcopy(s_src_us.info),
        deepcopy(s_src_us.aux),
    ))
    @assert axes(s_trg_unshaped) == axes(s_src)
    @assert s_trg_unshaped.v isa ArrayOfSimilarArrays
    @threads for i in eachindex(s_trg_unshaped, s_src)
        r = trafo_us(s_src_us.v[i], zero(Float32))
        s_trg_unshaped.v[i] .= r.v
        s_trg_unshaped.logd[i] = s_src_us.logd[i] - r.ladj
    end
    vs_trg.(s_trg_unshaped)
end

function _combined_trafo_ladj(trafo_ladj::OptionalLADJ, prev_ladj::OptionalLADJ, trg_v_isinf::Bool)
    if ismissing(trafo_ladj) || ismissing(prev_ladj)
        missing
    else
        ladj_sum = trafo_ladj + prev_ladj
        R = typeof(ladj_sum)
        if !isnan(ladj_sum)
            ladj_sum
        else
            # Should be safe to assume that target dist goes to zero at infinity, should win out over infinite prev_ladj:
            ladjs_should_cancel = (trafo_ladj == R(-Inf) && prev_ladj == R(+Inf) && trg_v_isinf)
            ladjs_should_cancel ? zero(R) : ladj_sum
        end
    end
end

function var_trafo_result(trg_v::Any, src_v::Any, trafo_ladj::OptionalLADJ, prev_ladj::OptionalLADJ)
    # TODO: Preserve precision of src_v while keeping dual-number properties of trg_v
    trg_ladj = _combined_trafo_ladj(trafo_ladj, prev_ladj, any_isinf(trg_v))
    (v = trg_v, ladj = trg_ladj)
end

function ChainRulesCore.rrule(::typeof(var_trafo_result), trg_v::Any, src_v::Any, trafo_ladj::OptionalLADJ, prev_ladj::OptionalLADJ)
    result = var_trafo_result(trg_v, src_v, trafo_ladj, prev_ladj)
    function var_trafo_result_pullback(thunked_ΔΩ)
        ΔΩ = ChainRulesCore.unthunk(thunked_ΔΩ)
        (NoTangent(), ΔΩ.v, ZeroTangent(),
            ismissing(trafo_ladj) ? missing : ΔΩ.ladj, ismissing(prev_ladj) ? missing : ΔΩ.ladj)
    end
    return result, var_trafo_result_pullback
end

# ToDo: Remove intermediate type `VariateTransform`?

"""
    abstract type VariateTransform{VT<:AbstractValueShape,VF<:AbstractValueShape}

*BAT-internal, not part of stable public API.*

Abstract parameterized type for change-of-variables transformations.

Subtypes (e.g. `SomeTrafo <: VariateTransform`) must implement:

* `BAT.apply_vartrafo_impl(trafo::SomeTrafo, v)`
* `BAT.apply_vartrafo_impl(inv_trafo::InverseVT{SomeTrafo}, v)`
* `ValueShapes.varshape(trafo::SomeTrafo)`

for real values and/or real-valued vectors `v`.
"""
abstract type VariateTransform{
    VT<:AbstractValueShape,VF<:AbstractValueShape
} <: AbstractVariateTransform end

function apply_vartrafo end

function apply_vartrafo_impl end

apply_vartrafo(trafo::VariateTransform{<:Any,<:ScalarShape{T}}, v::T, prev_ladj::OptionalLADJ) where {T<:Real} =
    apply_vartrafo_impl(trafo, v, prev_ladj)

function apply_vartrafo(trafo::VariateTransform{<:Any,<:ScalarShape{T}}, v::AbstractArray{<:T,0}, prev_ladj::OptionalLADJ) where {T<:Real}
    r = apply_vartrafo_impl(trafo, v[], prev_ladj)
    (v = fill(r.v), ladj = r.ladj)
end

apply_vartrafo(trafo::VariateTransform{<:Any,<:ArrayShape{T,N}}, v::AbstractArray{<:T,N}, prev_ladj::OptionalLADJ) where {T<:Real,N} =
    apply_vartrafo_impl(trafo, v, prev_ladj)

apply_vartrafo(trafo::VariateTransform{<:Any,<:ValueShapes.NamedTupleShape{names}}, v::NamedTuple{names}, prev_ladj::OptionalLADJ) where names =
    apply_vartrafo_impl(trafo, v, prev_ladj)

apply_vartrafo(trafo::VariateTransform{<:Any,<:ValueShapes.NamedTupleShape{names}}, v::ShapedAsNT{<:NamedTuple{names}}, prev_ladj::OptionalLADJ) where names =
    apply_vartrafo_impl(trafo, v, prev_ladj)

(trafo::VariateTransform)(v::Any) = apply_vartrafo(trafo, v, missing).v
(trafo::VariateTransform)(v::Any, prev_ladj::OptionalLADJ) = apply_vartrafo(trafo, v, prev_ladj)
(trafo::VariateTransform)(s::DensitySample) = _transform_density_sample(trafo, s)

struct IdentityVT{
    VTF <: AbstractValueShape
} <: VariateTransform{VTF,VTF}
    varshape::VTF
end

Base.inv(trafo::IdentityVT) = trafo

ValueShapes.varshape(trafo::IdentityVT) = trafo.varshape
ValueShapes.valshape(trafo::IdentityVT) = trafo.varshape

ValueShapes.unshaped(trafo::IdentityVT{<:ArrayShape{<:Any,1}}) = trafo

import Base.∘
@inline ∘(a::AbstractVariateTransform, b::IdentityVT) = a
@inline ∘(a::IdentityVT, b::IdentityVT) = a
@inline ∘(a::IdentityVT, b::AbstractVariateTransform) = b

@inline apply_vartrafo_impl(trafo::IdentityVT, v::Any, prev_ladj::OptionalLADJ) = (v = v, ladj = prev_ladj)

(trafo::IdentityVT)(s::DensitySample) = s

# Custom broadcast(::IdentityVT, DensitySampleVector), multithreaded:

function Base.copy(
    instance::Base.Broadcast.Broadcasted{
        <:Base.Broadcast.AbstractArrayStyle,
        <:Any,
        <:IdentityVT,
        <:Tuple{<:Union{ArrayOfSimilarVectors{<:Real},ShapedAsNTArray,DensitySampleVector}}
    }
)
    deepcopy(instance.args[1])
end


text:
include("gradient_hessian/gradient_hessian.jl")
include("log_likelihood.jl")
include("simulation.jl")
include("transition_probabilities.jl")
include("mle.jl")


text:
struct RestrictedSequencingNode{T<:Real}<:Node
    path                    ::Vector{Int}
    visited                 ::BitVector
    value                   ::T
    domain                  ::BitVector
end

Base.show(io::IO,node::RestrictedSequencingNode) = Base.print(
    io,
    "\n   path: ", node.path, "\n",
    "   value: ", node.value, "\n",
    "   visited: ", node.visited, "\n",
    "   domain: ", node.domain, "\n"
)

"""
    deepCopy(node::RestrictedSequencingNode)

Make a deep copy of a RestrictedSequencingNode

# Arguments
- `node::RestrictedSequencingNode`: The node to copy
"""
function deepCopy(node::RestrictedSequencingNode)
    return RestrictedSequencingNode(
        copy(node.path),
        copy(node.visited),
        node.value,
        copy(node.domain)
    )
end

"""
    getState(node::RestrictedSequencingNode)

Get the last location visited in the path leading to the given node

# Arguments
- `node::RestrictedSequencingNode`: The node
"""
function getState(node::RestrictedSequencingNode)
    return last(node.path)
end

"""
    inDomain(node::RestrictedSequencingNode, i::Int)

Check if a value is in the domain of a node (check if a location has been visited or not)

# Arguments
- `node::RestrictedSequencingNode`: The node whose domain is being checked
- `i::Int`: The location/value to check for in the domain
"""
function inDomain(node::RestrictedSequencingNode, i::Int)
    return node.domain[i]
end

"""
    remove!(node::RestrictedSequencingNode, i::Int)

Remove a value from the domain of a node

# Arguments
- `node::RestrictedSequencingNode`: the node
- `i::Int`: the value to remove
"""
function remove!(node::RestrictedSequencingNode, i::Int)
    node.domain[i] = false
end

"""
    add!(node::RestrictedSequencingNode, i::Int)

Add a a value to the domain of a node

# Arguments
- `node::RestrictedSequencingNode`: the domain
- `i::Int`: the value to add
"""
function add!(node::RestrictedSequencingNode, i::Int)
    node.domain[i] = true
end

"""
    initializeDomain(model::SequencingModel)

Initalize the domain for a given sequencing problem

# Arguments
- `model::SequencingModel`: the problem instance
"""
function initializeDomain(model::SequencingModel)
    d = trues(model.length)

    #handle precedence constraints
    for follower in keys(model.precedenceF2P)
        d[follower]=false
    end

    return d
end

"""
    createRestrictedRootNode(model::SequencingModel)

Create the root node for a problem using the given model

# Arguments
- `model::SequencingModel`: The problem instance
"""
function createRestrictedRootNode(model::SequencingModel)
    if hasStart(model)
        v = falses(model.length)
        v[model.seqStart] = true
        d = initializeDomain(model)
        d[model.seqStart] = false
        return RestrictedSequencingNode(Int[model.seqStart], v,0, d)
    else
        return RestrictedSequencingNode(Int[], falses(model.length),0, initializeDomain(model))
    end
end


text:
# Julia wrapper for header: /usr/include/scip/scip_disp.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeDisp(scip, name, desc, header, dispstatus, dispcopy, dispfree, dispinit, dispexit, dispinitsol, dispexitsol, dispoutput, dispdata, width, priority, position, stripline)
    ccall((:SCIPincludeDisp, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring, Cstring, SCIP_DISPSTATUS, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{SCIP_DISPDATA}, Cint, Cint, Cint, UInt32), scip, name, desc, header, dispstatus, dispcopy, dispfree, dispinit, dispexit, dispinitsol, dispexitsol, dispoutput, dispdata, width, priority, position, stripline)
end

function SCIPfindDisp(scip, name)
    ccall((:SCIPfindDisp, libscip), Ptr{SCIP_DISP}, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPgetDisps(scip)
    ccall((:SCIPgetDisps, libscip), Ptr{Ptr{SCIP_DISP}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNDisps(scip)
    ccall((:SCIPgetNDisps, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPautoselectDisps(scip)
    ccall((:SCIPautoselectDisps, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPchgDispMode(disp, mode)
    ccall((:SCIPchgDispMode, libscip), Cvoid, (Ptr{SCIP_DISP}, SCIP_DISPMODE), disp, mode)
end


text:
function get_volume(dimensions)
    reduce(*, dimensions)
end

function read_lines_from_cwd(filename)
    input_path = joinpath(dirname(Base.source_path()), filename)
    f = open(input_path)
    lines = readlines(f)
    close(f)
    lines
end

function main()
    lines = read_lines_from_cwd("input.txt")
    paper = 0
    ribbon = 0
    for line in lines
        dimensions = split(line, "x")
        dimensions = map(x->parse(Int32, x), dimensions)
        sort!(dimensions)  # sort in-place
        x, y, z = dimensions
        a = 2 * x * y
        b = 2 * x * z
        c = 2 * y * z
        surface = a + b + c
        smallest = Int32(minimum([a, b, c]) / 2)
        paper += surface + smallest

        volume = get_volume(dimensions)
        perimeter = 2 * (x + y)
        ribbon += volume + perimeter
    end
    println("Paper: ", paper)
    println("Ribbon: ", ribbon)
end

main()


text:
## Accurate Summation
# See what happens in two floatint Point
# precisions
using Printf
""" Run an example where we compute
a sum in a low-precision and compare with
a high-precision sum. You can use BigFloat
in case you need to check Float64. """
function sum_example_1(ns, Tlow, Thigh=Float64)
  for n in ns
    x = randn(Tlow, n)
    y = Thigh.(x) # convert to high-precision
    s1 = sum(x) # sum in low-precision
    s2 = sum(y) # sum in high-precision
    @printf("%10i - %20s - %20s\n", n, "$s1", "$s2")
  end
end
sum_example_1([10^3 10^4 10^5 10^6], Float16)
## Let's do the same example but where we control the sum
# and don't allocate vectors.
using Printf
function sum_example(ns, Tlow, Thigh=Float64)
  for n in ns
    s1 = zero(Tlow)
    s2 = zero(Thigh)
    for i =1:n
      x = randn(Tlow)
      s1 += x
      s2 += Thigh.(x)
      x = rand(Tlow)
      x *= 2*rand(Bool)-1 # random sign
      s1 += x
      s2 += Thigh.(x)
    end
    @printf("%10i - %20s - %20s\n", n, "$s1", "$s2")
  end
end
sum_example([10^3 10^4 10^5 10^6], Float16)

## Kahan Summation

"""
kahan_sum implements one step of the Kahan summation procedure.
function KahanSum(input)
  variables sum,c,y,t,i          // Local to the routine.
    sum = 0.0                    // Prepare the accumulator.
    c = 0.0                      // A running compensation for lost low-order bits.
    for i = 1 to input.length do // The array input has elements indexed input[1] to input[input.length].
        y = input[i] - c         // c is zero the first time around.
        t = sum + y              // Alas, sum is big, y small, so low-order digits of y are lost.
        c = (t - sum) - y        // (t - sum) cancels the high-order part of y; subtracting y recovers negative (low part of y)
        sum = t                  // Algebraically, c should always be zero. Beware overly-aggressive optimizing compilers!
    next i                       // Next time around, the lost low part will be added to y in a fresh attempt.
    return sum
"""
function kahan_sum(x,sum,c)
  y = x-c
  t = sum+y
  c = (t-sum)-y
  sum = t
  return sum,c
end

function ksum(xs::AbstractArray)
  sum = zero(eltype(xs))
  c = zero(eltype(xs))
  for x in xs
    sum,c = kahan_sum(x, sum, c)
  end
  return sum
end
ksum(randn(Float16,1000))

##
using Printf
function sum_example_kahan(ns, Tlow, Thigh=Float64)
  for n in ns
    s0 = zero(Tlow)
    s1 = zero(Tlow)
    c1 = zero(Tlow)
    s2 = zero(Thigh)
    for i =1:n
      x = randn(Tlow)
      s0 += x
      s1,c1 = kahan_sum(x, s1, c1)
      s2 += Thigh.(x)
      x = rand(Tlow)
      x *= 2*rand(Bool)-1 # random sign
      s0 += x
      s1,c1 = kahan_sum(x, s1, c1)
      s2 += Thigh.(x)
    end
    @printf("%10i - %8s - %8s - %20s\n", n, "$s0","$s1", "$s2")
  end
end
using Random
Random.seed!(0)
sum_example_kahan([10^3 10^4 10^5 10^6], Float16)


text:
using Random

"""

Arguments:
    - `f`: Function to be optimized
    - `g`: Gradient function for `f`
    - `c`: Constraints function. Evaluates all the constraint equations and return a Vector of values
    - `x0`: (Vector) Initial position to start from
    - `n`: (Int) Number of evaluations allowed. Remember `g` costs twice of `f`
    - `prob`: (String) Name of the problem. So you can use a different strategy for each problem
"""
function optimize(f, g, c, x0, n, prob)
    x_best = x0
    y_best = f(x0)
    cs = c(x0)

    for i in 1 : floor(n/3)-1
        x_next = x_best + randn(length(x_best)) .* g(x_best)
        y_next = f(x_next)
        if y_next < y_best && all(c(x_next) .<= 0)
            x_best, y_best = x_next, y_next
        end
    end
    return x_best
end


text:
#------------------------------
# Non-linear ODE system with two unkown parameters
#------------------------------

#------------------------------
# Setup
#------------------------------

using Pkg;

# all paths are relative to the repository main folder

Pkg.activate(".")
Pkg.instantiate()
Pkg.status()

using Distributions
using Random
using Flux
using DiffEqFlux
using OrdinaryDiffEq
using Plots
using LaTeXStrings
gr()

include("../src/simulation.jl")
include("../src/model.jl")
include("../src/plotting.jl")

#------------------------------
# Define and obtain ground-truth developments
#------------------------------

true_u0 = Float32[2, 2]
tspan = (0.0f0, 10.0f0)
true_odeparams_group1 = Float32[1.0, 1.0, 1.0, 0.5]
true_odeparams_group2 = Float32[0.5, 1.0, 1.0, 2.0]

lvprob1 = ODEProblem(lotka_volterra, true_u0, tspan, true_odeparams_group1)
lvprob2 = ODEProblem(lotka_volterra, true_u0, tspan, true_odeparams_group2)

# true trajectory of each group
dt=0.1
sol_group1 = solve(lvprob1, Tsit5(), saveat = dt);
sol_group2 = solve(lvprob2, Tsit5(), saveat = dt);
#plot(sol_group1)
#plot(sol_group2)

#------------------------------
# Simulate data based on true trajectories
#------------------------------

n = 200
p = 10
q, q_info = 50, 20
baseline_simulation = "groupsonly"

Random.seed!(1234)
xs, tvals, group1, group2 = generate_xs(n, p, true_u0, sol_group1, sol_group2, t_start=0.0, maxntps=1, σ_var=0.1, σ_ind=0.1); # default vals: t_start=1.5, t_end=10, maxntps=10, dt=0.1, σ_var=0.1, σ_ind=0.5
if baseline_simulation == "trueparams"
    x_baseline = generate_baseline(n, q, q_info, group1, true_odeparams_group1, true_odeparams_group2); # defaul vals: σ_info=0.1, σ_noise=0.1
elseif baseline_simulation == "groupsonly"
    x_baseline = generate_baseline(n, q, q_info, group1,  σ_info=0.5, σ_noise=0.5); # default vals: σ_info=1, σ_noise=1
else
    error("Please select either 'trueparams' or 'groupsonly' as mode for simulation of baseline variables")
end

# look at data:
data = simdata(xs, x_baseline, tvals, group1, group2);

plot(plot_truesolution(2, data, sol_group1, sol_group2, showdata=false),
    plot_truesolution(1,  data, sol_group1, sol_group2, showdata=false),
    plot_truesolution(2, data, sol_group1, sol_group2, showdata=true),
    plot_truesolution(1,  data, sol_group1, sol_group2, showdata=true),
    layout = (2,2),
    legend = false
)

#------------------------------
# Define and train model
#------------------------------

zdim = nODEparams = 2
m = init_vae(p, q, zdim, nODEparams, lvprob2, seed=84)
L = loss_wrapper(m)
ps = getparams(m)
opt = ADAM(0.0001)
trainingdata = zip(xs, x_baseline, tvals);
evalcb() = @show(mean(L.(xs, x_baseline, tvals)))
evalcb_zs() = eval_z_trajectories(xs, x_baseline, tvals, group1, sol_group1, sol_group2, m, dt)

for epoch in 1:40
    Flux.train!(L, ps, trainingdata, opt)
    evalcb()
    evalcb_zs()
end

#------------------------------
# Look at results
#------------------------------

plot(allindsplot(2, data, m, sol_group1, sol_group2),
    allindsplot(1, data, m, sol_group1, sol_group2, showlegend=false),
    layout = (1,2)
)

# individual plots
p1 = plot_individual_solutions(46, xs, x_baseline, tvals, group1, sol_group1, sol_group2, m, dt)
p2 = plot_individual_solutions(182, xs, x_baseline, tvals, group1, sol_group1, sol_group2, m, dt, showlegend=false)

# solutions across the entire dataset
p3 = allindsplot(2, data, m, sol_group1, sol_group2)
p4 = allindsplot(1, data, m, sol_group1, sol_group2, showlegend=false)
plot(p3,p4, layout = (1,2))


text:
module Jl2Py

export jl2py

using PyCall

macro __unaryop(jl_expr, op)
    quote
        operand = __jl2py($(esc(jl_expr)).args[2])
        if isa(operand, Vector)
            operand = operand[1]
        end
        unaryop = AST[].UnaryOp($op(), operand)
        return [unaryop]
    end
end

macro __binop(jl_expr, op)
    quote
        jl_expr = $(esc(jl_expr))
        left = __jl2py(jl_expr.args[2])
        if isa(left, Vector)
            left = left[1]
        end
        right = __jl2py(jl_expr.args[3])
        if isa(right, Vector)
            right = right[1]
        end
        binop = AST[].BinOp(left, $op(), right)
        return [binop]
    end
end

macro __compareop_from_comparison(jl_expr)
    quote
        jl_expr = $(esc(jl_expr))
        elts = __jl2py(jl_expr.args[1:2:end])
        ops = map(x -> OP_DICT[x](), jl_expr.args[2:2:end])
        compareop = AST[].Compare(elts[1], ops, elts[2:end])
        return [compareop]
    end
end

macro __compareop_from_call(jl_expr, op)
    quote
        jl_expr = $(esc(jl_expr))
        elts = __jl2py(jl_expr.args[2:end])
        compareop = AST[].Compare(elts[1], [$op()], elts[2:end])
        return [compareop]
    end
end

macro __multiop(jl_expr, op)
    quote
        jl_expr = $(esc(jl_expr))
        left = __jl2py(jl_expr.args[2])
        if isa(left, Vector)
            left = left[1]
        end
        right = __jl2py(jl_expr.args[3])
        if isa(right, Vector)
            right = right[1]
        end
        binop = AST[].BinOp(left, $op(), right)
        argcnt = length(jl_expr.args)
        for i in 4:argcnt
            right = __jl2py(jl_expr.args[i])
            if isa(right, Vector)
                right = right[1]
            end
            binop = AST[].BinOp(binop, $op(), right)
        end
        return [binop]
    end
end

const AST = Ref{PyObject}()
const OP_DICT = Dict{Symbol,PyObject}()

function __jl2py(args::AbstractVector)
    return map(x -> isa(x, Vector) ? x[1] : x, map(__jl2py, args))
end

function __jl2py(jl_constant::Union{Number,String})
    return AST[].Constant(jl_constant)
end

function __jl2py(jl_symbol::Symbol)
    return AST[].Name(jl_symbol)
end

function __jl2py(jl_expr::Expr)
    if jl_expr.head == :block || jl_expr.head == :toplevel
        py_exprs = PyObject[]
        for expr in jl_expr.args
            if isa(expr, LineNumberNode)
                continue
            end

            py_expr = __jl2py(expr)
            if isa(py_expr, Vector)
                append!(py_exprs, py_expr)
            else
                push!(py_exprs, AST[].Expr(py_expr))
            end
        end

        return py_exprs
    elseif jl_expr.head == :call
        if jl_expr.args[1] == :+
            if length(jl_expr.args) == 2
                @__unaryop(jl_expr, AST[].UAdd)
            else
                @__multiop(jl_expr, OP_DICT[jl_expr.args[1]])
            end
        elseif jl_expr.args[1] == :-
            if length(jl_expr.args) == 2
                @__unaryop(jl_expr, AST[].USub)
            else
                @__binop(jl_expr, OP_DICT[jl_expr.args[1]])
            end
        elseif jl_expr.args[1] == :*
            @__multiop(jl_expr, OP_DICT[jl_expr.args[1]])
        elseif jl_expr.args[1] ∈ [:/, :÷, :div, :%, :mod, :^]
            @__binop(jl_expr, OP_DICT[jl_expr.args[1]])
        elseif jl_expr.args[1] ∈ [:(==), :(===), :≠, :(!=), :(!==), :<, :<=, :>, :>=]
            @__compareop_from_call(jl_expr, OP_DICT[jl_expr.args[1]])
        else
            # FIXME: handle keyword arguments
            func = __jl2py(jl_expr.args[1])
            parameters = __jl2py(jl_expr.args[2:end])
            call = AST[].Call(func, parameters, [])
            return [AST[].Expr(call)]
        end
    elseif jl_expr.head == :comparison
        @__compareop_from_comparison(jl_expr)
    elseif jl_expr.head == :(=)
        targets = PyObject[]
        curr = jl_expr.args
        while isa(curr[2], Expr) && curr[2].head == :(=)
            push!(targets, __jl2py(curr[1]))
            curr = curr[2].args
        end
        last_target, value = __jl2py(curr)
        push!(targets, last_target)
        return [AST[].fix_missing_locations(AST[].Assign(targets, value))]
    elseif jl_expr.head == :vect
        listop = AST[].List(__jl2py(jl_expr.args))
        return [listop]
    else
        return PyObject[]
    end
end

function jl2py(jl_str::String)
    jl_ast = Meta.parse(jl_str)
    py_body = PyObject[]

    if isa(jl_ast, Union{Number,String})
        push!(py_body, __jl2py(jl_ast))
    elseif isa(jl_ast, Symbol)
        push!(py_body, __jl2py(jl_ast))
    elseif isa(jl_ast, Expr)
        append!(py_body, __jl2py(jl_ast))
    end

    _module = AST[].Module(py_body, [])
    py_str = AST[].unparse(_module)
    return py_str
end

function __init__()
    AST[] = pyimport("ast")
    OP_DICT[:+] = AST[].Add
    OP_DICT[:-] = AST[].Sub
    OP_DICT[:*] = AST[].Mult
    OP_DICT[:/] = AST[].Div
    OP_DICT[:÷] = AST[].FloorDiv
    OP_DICT[:div] = AST[].FloorDiv
    OP_DICT[:%] = AST[].Mod
    OP_DICT[:mod] = AST[].Mod
    OP_DICT[:^] = AST[].Pow
    OP_DICT[:(==)] = AST[].Eq
    OP_DICT[:(===)] = AST[].Eq
    OP_DICT[:≠] = AST[].NotEq
    OP_DICT[:(!=)] = AST[].NotEq
    OP_DICT[:(!==)] = AST[].NotEq
    OP_DICT[:<] = AST[].Lt
    OP_DICT[:<=] = AST[].LtE
    OP_DICT[:>=] = AST[].GtE
    OP_DICT[:>] = AST[].Gt
    OP_DICT[:in] = AST[].In
    OP_DICT[:∈] = AST[].In
    OP_DICT[:∉] = AST[].NotIn
end

end


text:
# functions to set and get associated user data in variable `instance` of `DataType`.

"""
    settypevar!(t::Type, value)

Associate type with a value. The second call for the same type is gnored.
The value can be retrieved by calling `gettypevar(t)`.
"""
function settypevar!(t::Type, value)
    get!(TypeVariablesTemp, t, value)
end
"""
    settypevar!(f::Function, t::Type)

Associate type with value `f()`. `f` is called only once per type.
The value can be retrieved by calling `gettypevar(t)`.
"""
function settypevar!(f::Function, t::Type)
    get!(f, TypeVariablesTemp, t)
end
const TypeVariablesTemp = IdDict{DataType, Any}()

"""
    gettypevar(t::Type)

Return value, which has previously been associated with this type
"""
@generated function gettypevar(::Type{R}) where {T,R<:Ring{T}}
    tv = TypeVariablesTemp[R]
    :( $tv )
end

"""
    new_class(t::Type{<:Ring{T}}, args...) where T<:RingClass

Store type variable `T(args...)` for `t` and return `t`.
The return value can be used to create ring elements of the new class.
Example:
```
    ZZp = new_class(ZZmod{:p,BigInt}, 1000000000000000000000000000057)
    element = ZZp(123456)^10000000
```
"""
function new_class(t::Type{<:Ring{T}}, args...) where T
    typevar() = T(args...)
    settypevar!(typevar, t)
    t
end
function new_class(f::Function, t::Type{<:Ring{T}}, args...) where T
    settypevar!(t) do
        T(f(args...)...)
    end
    t
end

const IdentSymbols = Union{Symbol,Base.BitInteger,BigInt}
"""
    sintern(a)

Return a symbol, which uniquly identifies the argument.
"""
sintern(m::IdentSymbols) = m
sintern(m::BigInt) = Symbol(m)
sintern(a::IdentSymbols...) = Symbol(tuple(a...))
sintern(a::AbstractVector{<:IdentSymbols}) = Symbol(a)
sintern(a::Tuple{Vararg{<:IdentSymbols}}) = Symbol(a)
sintern(a) = Symbol(hash(a))


text:
module Combinatorics

using Polynomials

include("numbers.jl")
include("factorials.jl")
include("combinations.jl")
include("permutations.jl")
include("partitions.jl")
include("multinomials.jl")
include("youngdiagrams.jl")

end #module


text:
using Colors, GtkUtilities, Gtk.ShortNames, Graphics

include("julia_set.jl")

# Create the window and canvas
win = @Window("julia", 400, 400)
c = @Canvas()
push!(win, c)

# Set up computation buffers
s, z_r, z_i, col = alloc(Float32, width(win))

# This is called for every drawing event (resizing the window,
# changing the zoom/pan)
@guarded Gtk.ShortNames.draw(c) do widget
    # Retrieve the display region
    xview, yview = guidata[widget, :xview], guidata[widget, :yview]
    set_coords(getgc(widget), xview, yview)
    # Render the image for this region
    iterate!(s, z_r, z_i, (xview.min,xview.max), (yview.min,yview.max))
    colorize!(col, s, z_r, z_i)
    # Blit it to the screen
    copy!(widget, col)
end

# Initialize panning & zooming
lim = (-2.0,2.0)
panzoom(c, lim, lim)
panzoom_mouse(c, factor=1.1)
panzoom_key(c)

# Activate the application
showall(win)


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -2.6308206135564163 <= q <= 3.997051239074799)

                     add_NL_constraint(m, :((0.8157606255059098 + 0.2914814987498171*$(x[1]) + -0.962674206030298*$(x[2]))/(1 + exp(-(0.8157606255059098 + 0.2914814987498171*$(x[1]) + -0.962674206030298*$(x[2])))) + (0.7907664778929373 + 0.4020245577506052*$(x[1]) + -0.1936195018990734*$(x[2]))/(1 + exp(-(0.7907664778929373 + 0.4020245577506052*$(x[1]) + -0.1936195018990734*$(x[2])))) + (-0.5840361484149175 + 0.5201746807648506*$(x[1]) + 0.3251106634438119*$(x[2]))/(1 + exp(-(-0.5840361484149175 + 0.5201746807648506*$(x[1]) + 0.3251106634438119*$(x[2])))) + (-0.33937564222473826 + -0.3486039500062561*$(x[1]) + 0.27024686767089534*$(x[2]))/(1 + exp(-(-0.33937564222473826 + -0.3486039500062561*$(x[1]) + 0.27024686767089534*$(x[2])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
as_struct(dt::Type{T} where T,x::BSON_VALUE_PRIMITIVE)=x
as_struct(dt::Type{T} where T<:Enum,x::BSON_VALUE_PRIMITIVE)=dt(x)
as_struct(dt::Type{T} where T<:Symbol,x::BSON_VALUE_PRIMITIVE)=dt(x)

as_struct(dt::Type{T} where T, arr::Array)=map(x->as_struct(Any,x),arr)

function as_struct(dt::Type{T} where T<:AbstractDict,x)
    _type=haskey(x,"_type") ? str_to_type(x["_type"]) : Dict{Any,Any}
    ret=_type()
    for (k,v) in x
        k=="_type" ? continue : nothing
        if v isa AbstractDict && haskey(v,"_k") && haskey(v,"_v")
            if v["_k"] isa AbstractDict && haskey(v["_k"],"_type")
               kc=as_struct(str_to_type(v["_k"]["_type"]),v["_k"])
            else
               kc=as_struct(Any,v["_k"])
            end

            if v["_v"] isa AbstractDict && haskey(v["_v"],"_type")
               vc=as_struct(str_to_type(v["_v"]["_type"]),v["_v"])
            else
                vc=as_struct(Any,v["_v"])
            end
            ret[kc]=vc
        elseif v isa AbstractDict && haskey(v,"_type")
            ret[k]=as_struct(str_to_type(v["_type"]),v)
        else
            ret[k]=as_struct(Any,v)
        end
    end

    return ret
end

function iter_data_types!(ret::Vector{DataType},arr::Vector)
    for r in arr
        if isconcretetype(r)
           push!(ret,r)
        else
            n_ret=get_concrete_types(r)
            length(n_ret)==0 ? nothing : append!(ret,n_ret)
        end
    end
end

function get_concrete_types(dt::Type)
    ret=Vector{DataType}()
    if typeof(dt)==Union
         iter_data_types!(ret,Base.uniontypes(dt))
    else
        iter_data_types!(ret,InteractiveUtils.subtypes(dt))
    end
    return ret
end

function str_to_type(str::AbstractString)
    ex=Meta.parse(str)

    if ex isa Symbol
        curr_module = isdefined(parentmodule(@__MODULE__), ex) ? parentmodule(@__MODULE__) : Main
        return getfield(curr_module,ex)
    elseif ex isa Expr
        if ex.head==:.
            namespace_arr=split(string(ex.args[1]),".")
            curr_module = isdefined(parentmodule(@__MODULE__), Symbol(first(namespace_arr))) ? parentmodule(@__MODULE__) : Main
            for n in namespace_arr
                curr_module=getfield(curr_module,Symbol(n))
            end
            return getfield(curr_module,ex.args[2].value)
        # Parametric
        elseif ex.head==:curly
            dt=str_to_type(string(ex.args[1]))
            parameters=Vector{Type}()
            for e in ex.args[2:end]
               push!(parameters,str_to_type(string(e)))
            end
            return dt{parameters...}
        end
    end
end

function as_struct(dt::Type,x)

    _type=haskey(x,"_type") ? x["_type"] : nothing

    if dt<:Type
        return str_to_type(get(x,"_value",nothing))
    elseif dt<:AbstractDict
        return as_struct(dt,x)
    elseif dt<:MongocUtils.BSON_PRIMITIVE
        return x
    elseif _type=="Symbol"
        return Symbol(x["_value"])
    elseif isconcretetype(dt)
        return dt([as_struct(fieldtype(dt,k),x[string(k)]) for k in fieldnames(dt)]...)
    elseif dt==Any
        if _type==nothing
            if x isa AbstractDict
                return as_struct(Dict,x)
            else
                return x
            end
        else
            return as_struct(str_to_type(_type),x)
        end
    ## dt is abstract
    else
        if _type==nothing
            types_arr=get_concrete_types(dt)
            try
               for t in types_arr
                  return  as_struct(t,x)
               end
            catch;
               return x
            end
        elseif _type=="Type"
            return str_to_type(get(x,"_value",nothing))
        else
            stored_type=str_to_type(_type)
            @assert stored_type<:dt "Stored $(string(stored_type)) is not a subtype of abstract type $(string(dt))"
            return as_struct(stored_type,x)
        end
    end
end


text:
export Data

"""
Data structure to store the solution of the problem along time
"""
mutable struct Data

    U
    datasize :: Int
    datalength :: Int

    function Data( v )

        (datalength , datasize ) = size(v)

        U = [copy(v)]

        new(U, datasize, datalength)

    end

end

import Base.length

length( data :: Data ) = length( data.U )


text:
module Runner

using Compat; import Compat.String

import ..Utilities: Utilities, Head, @H_str, issymbol, isexpr

export State, withscope, withref, findmethods, findtuples, findvcats

include("state.jl")
include("lookup.jl")

end


text:
import ThorinDistributions
import Polynomials
import DynamicPolynomials
import Combinatorics
using Test

@testset "kStats.jl tests" begin
    @testset "better_stirling correctness" begin
        @test ThorinDistributions._better_stirling(10,5)[10,5]==Combinatorics.stirlings2(10,5)
    end

    @testset "_fd correctness" begin
        @test ThorinDistributions._fd(3,0) == Polynomials.Polynomial([0,-6,11,-6,1])
        @test ThorinDistributions._fd(3,1) == Polynomials.Polynomial([-6,11,-6,1])
        @test ThorinDistributions._fd(3,2) == Polynomials.Polynomial([6,-5,1])
    end

    @testset "_compute_k correctness" begin
        DynamicPolynomials.@polyvar N
        DynamicPolynomials.@polyvar S[1:9]
        result = ThorinDistributions._compute_k(N,S);
        Test.@test result[1] == S[1]/N
        Test.@test result[2] == (N * S[2] - S[1]^2)/(N*(N-1))
        Test.@test result[3] == (2S[1]^3 - 3N*S[1]*S[2] +N^2 * S[3]) / (N*(N-1)*(N-2))
        Test.@test result[4] == (-6*S[1]^4 + 12 * N * S[1]^2 * S[2] - 3 * N * (N-1) * S[2]^2 - 4 * N * (N+1) * S[1] * S[3] + N^2*(N+1)*S[4])/(N*(N-1)*(N-2)*(N-3))
    end
end


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Main parser file

let parser_dir = "parser"
    include(joinpath(parser_dir, "scene.jl"))
    include(joinpath(parser_dir, "utils.jl"))
    include(joinpath(parser_dir, "expectations.jl"))
    include(joinpath(parser_dir, "parse_constructors.jl"))
    include(joinpath(parser_dir, "parse_commands.jl"))
    include(joinpath(parser_dir, "parse_scene.jl"))
    include(joinpath(parser_dir, "cli_parsing.jl"))
end


text:
"""
    OptOrtho(loss, grad, s)

This function returns the opimizer of loss, whose feasible
set is O(s), using the GradientDescent method from
Manifold optimization package Optim

# Arguments
-   `loss`: loss function whose input is a orthogonal matrix
-   `grad`: Euclidean gradient function of loss function w.r.t
    to the orthogonal matrix
-   `s`: The feasible set is O(s)
-   `opt_config`: configuration passed to Optim

# Outputs:
-   `Xopt`: orthogonal matrix optmize loss function
"""
function OptOrtho(loss, grad, s; opt_config=nothing)
    # default optim config
    if opt_config === nothing
        opt_config = Optim.Options(f_tol = 0.0,
                                    g_tol = 1e-8,
                                    iterations = 1_000)
    end

    # use Stiefel manifold
    manif = Optim.Stiefel()

    # initialize a orthogonal matrix of size s-by-s
    X0 = complex(Matrix(qr(randn(s,s)).Q))

    # find the minimum of loss function using the
    # gradient descent method, the Euclidean gradient
    # is given by grad
    optimres = Optim.optimize(loss, grad, X0,
        Optim.GradientDescent(manifold=manif), opt_config; inplace = false)

    # extract and return the minimizer
    Xopt = Optim.minimizer(optimres)
    return real(Xopt)
end

"""
    OptOrtho(loss, grad, s)

This function returns the opimizer of loss, whose feasible
set is S(s - 1), using the GradientDescent method from
Manifold optimization package Optim

# Arguments
-   `loss`: loss function whose input is a orthogonal matrix
-   `grad`: Euclidean gradient function of loss function w.r.t
    to the orthogonal matrix
-   `s`: The feasible set is S(s-1)
-   `opt_config`: configuration passed to Optim

# Outputs:
-   `Xopt`: orthogonal matrix optmize loss function
"""
function OptSphere(loss, grad, s; opt_config=nothing)
    # default optim config
    if opt_config === nothing
        opt_config = Optim.Options(f_tol = 0.0,
                                    g_tol = 1e-8,
                                    iterations = 1_000)
    end

    # use Sphere manifold
    manif = Optim.Sphere()

    # initialize a unit vector of size s
    x0 = randn(s)
    x0 = complex(x0/sqrt(x0'*x0))

    # find the minimum of loss function using the
    # gradient descent method, the Euclidean gradient
    # is given by grad
    optimres = Optim.optimize(loss, grad, x0,
        Optim.GradientDescent(manifold=manif), opt_config; inplace = false)

    # extract and return the minimizer
    xopt = Optim.minimizer(optimres)
    return real(xopt)
end


text:
"""
    load_rda(rda_file)

Load header dictionary and spectroscopy data from a Siemens .rda file.  Returns
`header,data`, where `header` is a simple dictionary containing the ASCII
header, and `data` is the time domain induced magnetization signal.

Single voxel spectroscopy data will be loaded as a vector containing the single
magnentization signal.

For spectroscopy with a spatial component (that is, CSI data with at least one
header field `CSIMatrixSize[i]` not equal to 1), an array of size N×I×K×J will
be returned where I,J,K are the CSI indices 0,1,2 and the FID length is N.
"""
function load_rda(io::IO)
    firstline = readline(io)
    if firstline != ">>> Begin of header <<<"
        error("Expected RDA header, got: $firstline")
    end
    # Read header
    header = OrderedDict{String,String}()
    while true
        if eof(io)
            throw(EOFError("RDA file ended while reading header"))
        end
        # TODO: Is the encoding correct here?
        line = readline(io)
        if line == ">>> End of header <<<"
            break
        end
        # To split array indices...
        # r"^(\w+)(?:\[(\d+)\])?:\s*(.*)$",
        m = match(r"^(\S+)\s*:\s*(.*)$", line)
        if m == nothing
            @error "BUG: could not match RDA header line \"$line\""
        end
        header[m[1]] = m[2]
    end
    vecsize = parse(Int, header["VectorSize"])
    csi_dims = [parse(Int, get(header, "CSIMatrixSize[$i]", "1")) for i=0:2]
    # RDA should always be channel combined (& averaged?), so expect a single
    # FID per voxel.
    data_dims = prod(csi_dims) == 1 ? (vecsize,) : (vecsize, csi_dims...)
    raw_data = read(io)
    if sizeof(raw_data) != prod(data_dims)*sizeof(ComplexF64)
        error("""
              Wrong amount of raw data in .rda file.
              Expected $data_dims ComplexF64's. Got $(length(raw_data)) bytes instead.""")
    end
    data = Array{ComplexF64}(undef, data_dims)
    copy!(data, reshape(reinterpret(ComplexF64, raw_data), data_dims))
    header, data
end

load_rda(fname::AbstractString) = open(load_rda, fname)

"""
    save_rda(file, header, data)

Save a single FID as Siemens format .rda.  `header` should be a dictionary of
key values for the rda ASCII header.
"""
function save_rda(io::IO, header::AbstractDict, data)
    if parse(Int, header["VectorSize"]) != size(data,1)
        error("rda header VectorSize mismatches with length(data)")
    end
    for dim in 2:4
        csi_key = "CSIMatrixSize[$(dim-2)]"
        h_size = parse(Int, get(header, csi_key, "1"))
        d_size = size(data, dim)
        if h_size != d_size
            error("rda header entry $csi_key=$h_size is inconsistent with data dimension $dim of length $d_size")
        end
    end
    write(io, ">>> Begin of header <<<\r\n")
    for (key,val) in header
        write(io, "$key: $val\r\n")
    end
    write(io, ">>> End of header <<<\r\n")
    write(io, ComplexF64.(data))
    nothing
end

save_rda(fname::String, header, data) = open(io->save_rda(io, header, data), fname, "w")

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
struct BabyPOMG
    babyPOMDP::CryingBaby
end

function MultiCaregiverCryingBaby()
    BabyPOMDP = CryingBaby()
    return BabyPOMG(BabyPOMDP)
end

n_agents(pomg::BabyPOMG) = 2

ordered_states(pomg::BabyPOMG) = [SATED, HUNGRY]
ordered_actions(pomg::BabyPOMG, i::Int) = [FEED, IGNORE, SING]
ordered_joint_actions(pomg::BabyPOMG) = vec(collect(Iterators.product([ordered_actions(pomg, i) for i in 1:n_agents(pomg)]...)))

n_actions(pomg::BabyPOMG, i::Int) = length(ordered_actions(pomg, i))
n_joint_actions(pomg::BabyPOMG) = length(ordered_joint_actions(pomg))

ordered_observations(pomg::BabyPOMG, i::Int) = [CRYING, QUIET]
ordered_joint_observations(pomg::BabyPOMG) = vec(collect(Iterators.product([ordered_observations(pomg, i) for i in 1:n_agents(pomg)]...)))

n_observations(pomg::BabyPOMG, i::Int) = length(ordered_observations(pomg, i))
n_joint_observations(pomg::BabyPOMG) = length(ordered_joint_observations(pomg))

function transition(pomg::BabyPOMG, s, a, s′)
    # Regardless, feeding makes the baby sated.
    if a[1] == FEED || a[2] == FEED
        if s′ == SATED
            return 1.0
        else
            return 0.0
        end
    else
        # If neither caretaker fed, then one of two things happens.
        # First, a baby that is hungry remains hungry.
        if s == HUNGRY
            if s′ == HUNGRY
                return 1.0
            else
                return 0.0
            end
        # Otherwise, it becomes hungry with a fixed probability.
        else
            probBecomeHungry = 0.5 #pomg.babyPOMDP.p_become_hungry
            if s′ == SATED
                return 1.0 - probBecomeHungry
            else
                return probBecomeHungry
            end
        end
    end
end

function joint_observation(pomg::BabyPOMG, a, s′, o)
    # If at least one caregiver sings, then both observe the result.
    if a[1] == SING || a[2] == SING
        # If the baby is hungry, then the caregivers both observe crying/silent together.
        if s′ == HUNGRY
            if o[1] == CRYING && o[2] == CRYING
                return pomg.babyPOMDP.p_cry_when_hungry_in_sing
            elseif o[1] == QUIET && o[2] == QUIET
                return 1.0 - pomg.babyPOMDP.p_cry_when_hungry_in_sing
            else
                return 0.0
            end
        # Otherwise the baby is sated, and the baby is silent.
        else
            if o[1] == QUIET && o[2] == QUIET
                return 1.0
            else
                return 0.0
            end
        end
    # Otherwise, the caregivers fed and/or ignored the baby.
    else
        # If the baby is hungry, then there's a probability it cries.
        if s′ == HUNGRY
            if o[1] == CRYING && o[2] == CRYING
                return pomg.babyPOMDP.p_cry_when_hungry
            elseif o[1] == QUIET && o[2] == QUIET
                return 1.0 - pomg.babyPOMDP.p_cry_when_hungry
            else
                return 0.0
            end
        # Similarly when it is sated.
        else
            if o[1] == CRYING && o[2] == CRYING
                return pomg.babyPOMDP.p_cry_when_not_hungry
            elseif o[1] == QUIET && o[2] == QUIET
                return 1.0 - pomg.babyPOMDP.p_cry_when_not_hungry
            else
                return 0.0
            end
        end
    end
end

function joint_reward(pomg::BabyPOMG, s, a)
    r = [0.0, 0.0]

    # Both caregivers do not want the child to be hungry.
    if s == HUNGRY
        r += [pomg.babyPOMDP.r_hungry, pomg.babyPOMDP.r_hungry]
    end

    # One caregiver prefers to feed.
    if a[1] == FEED
        r[1] += pomg.babyPOMDP.r_feed / 2.0
    elseif a[1] == SING
        r[1] += pomg.babyPOMDP.r_sing
    end

    # One caregiver prefers to sing.
    if a[2] == FEED
        r[2] += pomg.babyPOMDP.r_feed
    elseif a[2] == SING
        r[2] += pomg.babyPOMDP.r_sing / 2.0
    end

    # Note that caregivers only experience a cost if they do something.

    return r
end

joint_reward(pomg::BabyPOMG, b::Vector{Float64}, a) = sum(joint_reward(pomg, s, a) * b[s] for s in ordered_states(pomg))

function POMG(pomg::BabyPOMG)
    return POMG(
        pomg.babyPOMDP.γ,
        vec(collect(1:n_agents(pomg))),
        ordered_states(pomg),
        [ordered_actions(pomg, i) for i in 1:n_agents(pomg)],
        [ordered_observations(pomg, i) for i in 1:n_agents(pomg)],
        (s, a, s′) -> transition(pomg, s, a, s′),
        (a, s′, o) -> joint_observation(pomg, a, s′, o),
        (s, a) -> joint_reward(pomg, s, a)
    )
end


text:
# This file was generated, do not modify it. # hide
y_wind = data.Wind_gen
X = data[:, [:Windspeed, :Temperature, :Radiation_dir, :Radiation_dif]];


text:
using GardenOptim
using Test

@testset "classification" begin
    clf = GardenOptim.Classification(
        Dict(
            "type"=>"Life",
            "name"=>"God",
            "children"=>[
                Dict(
                    "type"=>"Family",
                    "name"=>"Homo",
                    "children"=>[
                        Dict(
                            "type"=>"Species",
                            "name"=>"Human",
                            "bio"=>"Homo sapiens",
                            "children"=>[]
                        )
                    ]
                )
            ]
        )
    )
    @testset "constructor" begin
        @test_throws UndefRefError clf.parent
        @test clf.children[1].parent === clf
        @test length(clf.children) == 1
        @test clf.children[1].type == :family
        @test clf.children[1].children[1].bio == "Homo sapiens"
        @test clf.children[1].children[1].name == :human
        @test clf.children[1].children[1].children == []
    end
    @testset "getfirstparent" begin
        @test GardenOptim.getfirstparent(:human, clf) === clf.children[1]
        @test GardenOptim.getfirstparent(:homo, clf) === clf.children[1]
        @test_throws UndefRefError GardenOptim.getfirstparent(:god, clf)
    end
end

max_stars_repo_name: dlozeve/GardenOptim

max_issues_repo_name: dlozeve/GardenOptim

max_forks_repo_name: dlozeve/GardenOptim


text:
# Copyright 2020 Google LLC

# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

#     https://www.apache.org/licenses/LICENSE-2.0

# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
@enum State floor empty occupied

function Base.show(io::Core.IO, state::State)
    if state == empty
        print(io, 'L')
    elseif state == occupied
        print(io, '#')
    else
        print(io, '.')
    end
end

function Base.show(io::Core.IO, state::AbstractArray{State, 2})
    height, width = size(state)
    for i in 1:height
        for j in 1:width
            show(io, state[i, j])
        end
        print(io, '\n')
    end
end

function parsechar(c)::State
    if c == '.'
        return floor
    elseif c == 'L'
        return empty
    elseif c == '#'
        return occupied
    end
    throw(ErrorException("Unrecognized character"))
end

function adjacents(state::AbstractArray{State, 2}, i, j)
    height, width = size(state)
    [(-1, -1), (-1, 0), (-1, +1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)] |>
        dirs->map(c->(x=i+c[1], y=j+c[2]), dirs) |>
        coords->filter(c -> !(c.x < 1 || c.y < 1 || c.x > height || c.y > width), coords) |>
        coords->filter(c -> state[c.x, c.y] != floor, coords) |>
        coords->map(c->(c.x, c.y), coords)
end

function nextstate(state::State, neighbors::AbstractArray{State}, threshold_to_leave::Int)
    if state == empty
        if count(s->s==occupied, neighbors) == 0
            return occupied
        end
    elseif state == occupied
        if count(s->s==occupied, neighbors) >= threshold_to_leave
            return empty
        end
    end
    return state
end

function coords_visible(state::AbstractArray{State, 2}, i, j)
    height, width = size(state)
    neighbor_dirs = [(-1, -1), (-1, 0), (-1, +1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)]
    coords = Tuple{Int, Int}[]
    for dir ∈ neighbor_dirs
        x = i
        y = j
        while true
            x = x + dir[1]
            y = y + dir[2]
            if x < 1 || x > height || y < 1 || y > width
                break
            end
            if state[x, y] != floor
                push!(coords, (x, y))
                break
            end
        end
    end
    coords
end

function find_relevant_coords(state::AbstractArray{State, 2}, relevancy_function)
    height, width = size(state)
    relevant_coords = Dict{Tuple{Int, Int}, AbstractArray{Tuple{Int, Int}}}()
    for i ∈ 1:height
        for j ∈ 1:width
            relevant_coords[(i, j)] = relevancy_function(state, i, j)
        end
    end
    relevant_coords
end

function parseinput(lines)
    height = length(lines)
    width = length(lines[1])
    state = Array{State}(undef, (height, width))
    for i ∈ 1:height
        for j ∈ 1:width
            state[i, j] = parsechar(lines[i][j])
        end
    end
    state
end

function evolve(state::AbstractArray{State, 2}, relevant_coords, threshold_to_leave::Int)
    next = copy(state)
    for i ∈ 1:size(state)[1]
        for j ∈ 1:size(state)[2]
            neighbors = map(ind->state[ind[1], ind[2]], relevant_coords[(i, j)])
            next[i, j] = nextstate(state[i, j], neighbors, threshold_to_leave)
        end
    end
    next
end

function part1(filename="day11input.txt")
    state = parseinput(readlines(filename))
    relevant_coords = find_relevant_coords(state, adjacents)

    prev = nothing
    while prev != state
        prev = state
        state = evolve(state, relevant_coords, 4)
    end
    return count(s->s==occupied, state)
end

function part2(filename="day11input.txt")
    state = parseinput(readlines(filename))
    relevant_coords = find_relevant_coords(state, coords_visible)
    prev = nothing
    while prev != state
        prev = state
        state = evolve(state, relevant_coords, 5)
    end
    return count(s->s==occupied, state)
end


text:
#=
Use Clang.jl to wrap the GDAL C API.

Several custom transformations are applied that should make using this package more convenient.

- docstrings are added, created from GDAL Doxygen XML output
- function names are lowercased following Julia conventions
- functions that return a Cstring are wrapped in unsafe_string to return a String
- functions that return a Ptr{Cstring} are wrapped in unsafe_loadstringlist to return a Vector{String}
- functions that return a Ptr{Cvoid} are wrapped in failsafe to catch errors
=#

using Clang
using Clang.Generators
using MacroTools
using Logging
using EzXML
using GDAL_jll
using JuliaFormatter

const output_dir = joinpath(@__DIR__, "..", "src")
const xmlpath = joinpath(@__DIR__, "doxygen.xml")

# set up a global logger to log.txt to store the large amount of logging
loghandle = open(joinpath(@__DIR__, "log.txt"), "w")
logger = SimpleLogger(loghandle)
global_logger(logger)

# several functions for building docstrings
include(joinpath(@__DIR__, "doc.jl"))

includedir = joinpath(GDAL_jll.artifact_dir, "include")
headerfiles =
    joinpath.(
        includedir,
        [
            "cpl_conv.h",
            "cpl_error.h",
            "cpl_minixml.h",
            "cpl_port.h",
            "cpl_progress.h",
            "cpl_virtualmem.h",
            "cpl_vsi.h",
            "gdal.h",
            "gdal_alg.h",
            "gdalwarper.h",
            "gdal_vrt.h",
            "gdal_utils.h",
            "ogr_api.h",
            "ogr_core.h",
            "ogr_srs_api.h",
        ],
    )

for headerfile in headerfiles
    if !isfile(headerfile)
        error("Header file missing `($headerfile)")
    end
end

# functions that return a Cstring, of which we know we must free them
# because for instance the documentation says: "Free with CPLFree() or VSIFree()"
# use the renamed (lowercased) function names here
# if it is not in this list, it will not be freed, add function names to make it more complete
const string_free_me = Symbol[:gdalinfo]

"""
Custom rewriter for Clang.jl's C wrapper

Gets called with all expressions in a header file, or all expressiong in a common file.
If available, it adds docstrings before every expression, such that Clang.jl prints them
on top of the expression. The expressions themselves get sent to `rewriter(::Expr)`` for
further treatment.
"""
function rewriter(xs::Vector)
    rewritten = Any[]
    for x in xs
        # Clang.jl inserts strings like "# Skipping MacroDefinition: X"
        # keep these to get a sense of what we are missing
        if x isa String
            push!(rewritten, x)
            continue
        end
        @assert x isa Expr

        name = cname(x)
        node = findnode(name, doc)
        docstr = node === nothing ? "" : build_docstring(node)
        isempty(docstr) || push!(rewritten, addquotes(docstr))
        x2 = rewriter(x)
        push!(rewritten, x2)
    end
    rewritten
end

"Rewrite expressions in the ways listed at the top of this file."
function rewriter(x::Expr)
    if @capture(x, function f_(fargs__)
        ccall(fname_, rettype_, argtypes_, argvalues__)
    end)
        # it is a function wrapper around a ccall

        # lowercase the function name
        f2 = Symbol(lowercase(String(f)))

        # bind the ccall such that we can easily wrap it
        cc = :(ccall($fname, $rettype, $argtypes, $(argvalues...)))

        cc2 = if occursin(r"^cpl.*error", String(f2))
            # We will take care of GDAL error states and memory freeing by
            # calling `aftercare`, except for `cpl*error*` functions, as these
            # will be called by `aftercare` itself, so we prevent a loop.
            if rettype == :Cstring
                :(unsafe_string($cc))
            else
                cc
            end
        elseif rettype == :Cstring
            # boolean to let aftercare know if it should call vsifree
            free = f2 in string_free_me
            :(aftercare($cc, $free))
        else
            :(aftercare($cc))
        end

        # stitch the modified function expression back together
        :(function $f2($(fargs...))
            $cc2
        end) |> prettify
    else
        # do not modify expressions that are no ccall function wrappers
        x
    end
end

# parse GDAL's Doxygen XML file
const doc = readxml(xmlpath)

# custom function that prevents overwriting GDAL.jl by mapping gdal.h to gdal_h.jl
function header_outputfile(h)
    stem = splitext(basename(h))[1]
    if stem == "gdal"
        joinpath(output_dir, "gdal_h.jl")
    else
        joinpath(output_dir, stem * ".jl")
    end
end

function rewriter(dag::Clang.ExprDAG)
    for node in get_nodes(dag)
        # Macros are not generated, so no need to skip
        map!(rewriter, node.exprs, node.exprs)
    end
end

# do not wrap, handled in prologue.jl
@add_def stat
@add_def _stat64
@add_def time_t

options = load_options(joinpath(@__DIR__, "generator.toml"))

# add compiler flags, e.g. "-DXXXXXXXXX"
args = get_default_args()
push!(args, "-I$includedir")

# create context
ctx = create_context(headerfiles, args, options)

# run generator
build!(ctx, BUILDSTAGE_NO_PRINTING)
rewriter(ctx.dag)
build!(ctx, BUILDSTAGE_PRINTING_ONLY)

add_doc(joinpath(@__DIR__, "..", "src", "libgdal.jl"))

format(joinpath(@__DIR__, ".."))

close(loghandle)


text:
# Define the mechanism by which code snippet is evaluated
#
export initialize_boxing_module

#+Evaluate, Obsolete
# -> evaluate_new
# Parse and execute a Julia code snippet
#
function evaluate(code::Union{SubString{String},String},
                  m::Module)::String
    result=m.eval(Meta.parse(code,raise=true))
    io=IOBuffer()
    show(io,"text/plain",result)
    return String(take!(io))
end

# +Evaluate, API   L:initialize_boxing_module
#
# Initialize a boxing module. This module is used to run Julia comment
# code snippet (tagged by "#!" or by "# !")
#
# *Example:*
# #+BEGIN_SRC julia :eval never :exports code
# initialize_boxing_module(boxingModule="MyBoxing",
#                          usedModules=["RequiredPackage_1",
#                                       "RequiredPackage_2",...])
# #+END_SRC
#
# creates
#
# #+BEGIN_SRC julia :eval never :exports code
# module MyBoxing
# using RequiredPackage_1,RequiredPackage_2,...
# end
# #+END_SRC
#
# and future "# !" statements are executed after using MyBoxing:
# #+BEGIN_SRC julia :eval never :exports code
# using MyBoxing
# # !statements
# #+END_SRC
function initialize_boxing_module(;
                                  boxingModule::String="BoxingModule",
                                  usedModules::Vector{String}=String[],
                                  force::Bool=false)::Nothing
    # Common errors
    @assert boxingModule!=""
    @assert usedModules!=String[""]

    # initialize a module (for boxing)
    module_exists = isdefined(J4Org,Symbol(boxingModule))

    if force||(!module_exists)
        if length(usedModules)>0
            usedModules_asString = "using $(join(usedModules,", ")) "
        else
            usedModules_asString = ""
        end
        eval(Meta.parse("module $(boxingModule) $(usedModules_asString) end"))
    else
        # force=false && module_exists = true
        # -> nothing is done... interpreted as an error if usedModules if different
        # --> however for the moment we do not know how to get module "used" module_name
        #     thus we trigger an error if usedModules is not empty
        @assert isempty(usedModules) "Tried to define an already existing module, maybe use the force=true kwarg"
    end
    nothing
end

# +Evaluate
#
# Parse and execute a Julia code snippet
#
# The format is raw comment
# #+BEGIN_EXAMPLE
# # some text
# #!5+6
# #!rand(5)
# # some text
# #+END_EXAMPLE
#
# Example:
# !5+6
# !rand(5)
#
function with_hash_evaluate(comment::String,
                            boxingModule::String)::String
    # Common error
    @assert boxingModule!=""
    @assert isdefined(J4Org,Symbol(boxingModule)) "Module $(boxingModule) does not exist, create a new one with initialize_boxing_module()"

    boxingModule_asType=eval(Meta.parse("J4Org.$(boxingModule)"))

    # process comment, line by line
    output=Array{String,1}()
    comment_line_by_line=split(comment,"\n")
    n_comment_line_by_line=length(comment_line_by_line)
    code_to_execute=r"^#[ ]?!(.*)$"
    i=1
    while i<=n_comment_line_by_line

        # only process line beginning with "#!" or by "# !"
        # (other lines are forwarded without modification)
        if !occursin(code_to_execute,comment_line_by_line[i])
            push!(output,comment_line_by_line[i])
            i=i+1
            continue
        end

        # Here process code
        @assert occursin(code_to_execute,comment_line_by_line[i])

        output_code=Array{String,1}()
        output_result=Array{String,1}()

        push!(output_code,"# #+BEGIN_SRC julia")
        push!(output_result,"# #+BEGIN_SRC julia")

        while occursin(code_to_execute,comment_line_by_line[i])
            code_local = match(code_to_execute,comment_line_by_line[i])[1]
            push!(output_code,"# "*code_local)
            show_code_local_result = !occursin(r";\s*$",code_local) # Caveat: do not take into account final comment

            try
                result_line_by_line=evaluate(code_local,boxingModule_asType)
            catch e
                io=IOBuffer()
                show(io,"text/plain",e)
                result_line_by_line="***ERROR***  $(String(take!(io)))"
                show_code_local_result = true; # if an error occurred, show it!
                warning_message("An error occurred evaluating $(code_local)")
            end
            result_line_by_line=split(result_line_by_line,"\n")

            if show_code_local_result
                result_line_by_line=map(x->"# "*x,result_line_by_line)
                output_result=vcat(output_result,result_line_by_line)
            end
            i=i+1
        end

        push!(output_code,"# #+END_SRC")
        push!(output_result,"# #+END_SRC")

        output=vcat(output,output_code,output_result)
    end

#    return foldr((x,y)->x*"\n"*y,output,"")
    return join(output,"\n")
end


text:
module BudAutomata

using StaticArrays

include("bud_levenshtein_automata.jl")

end # module


text:
struct Model
    fw
    fc
    fs
end

function Model{T}(wordembeds::Matrix{T}, charembeds::Matrix{T}, ntags::Int)
    fw = @graph x begin
        Lookup(wordembeds)(x)
    end

    d = size(charembeds, 1)
    fc = @graph x begin
        x = Lookup(charembeds)(x)
        x = Conv1D(T,5,d,5d,2,1)(x)
        max(x, 2)
    end

    d = size(wordembeds,1) + size(charembeds,1)*5
    dh = 300
    fs = @graph x begin
        x = Conv1D(T,5,d,dh,2,1)(x)
        x = relu(x)

        x1 = dropout(x, 0.2)
        x1 = Conv1D(T,5,dh,dh,2,1)(x1)
        x1 = relu(x1)
        x1 = Standardize(T,(dh,10))(x1)
        x += x1

        Linear(T,dh,ntags)(x)
    end
    Model(fw, fc, fs)
end

function (m::Model)(data::Dataset)
    w = m.fw(data.w)
    c = m.fc(data.c)
    c.batchdims = w.batchdims
    x = cat(1, w, c)
    y = m.fs(x)
    if Merlin.config.train
        softmax_crossentropy(data.t, y)
    else
        vec(argmax(y.data,1))
    end
end

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
using ExecutableSpecifications.Gherkin: parsescenario, issuccessful, Given, When, Then, ByLineParser, ScenarioStep

@testset "Scenario Outline" begin
    @testset "Outline has a Given step; Step is parsed" begin
        text = """
        Scenario Outline: This is one scenario outline
            Given a precondition with field <Foo>

        Examples:
            | Foo |
            | 1   |
            | 2   |
        """
        byline = ByLineParser(text)

        result = parsescenario(byline)

        @test issuccessful(result)
        scenario = result.value
        @test scenario.steps == ScenarioStep[Given("a precondition with field <Foo>")]
    end

    @testset "Scenario Outline has description; Description is parsed" begin
        text = """
        Scenario Outline: This is one scenario outline
            Given a precondition with field <Foo>

        Examples:
            | Foo |
            | 1   |
            | 2   |
        """
        byline = ByLineParser(text)

        result = parsescenario(byline)

        @test issuccessful(result)
        scenario = result.value
        @test scenario.description == "This is one scenario outline"
    end

    @testset "Scenario Outline has tags; Tags are parsed" begin
        text = """
        @tag1 @tag2
        Scenario Outline: This is one scenario outline
            Given a precondition with field <Foo>

        Examples:
            | Foo |
            | 1   |
            | 2   |
        """
        byline = ByLineParser(text)

        result = parsescenario(byline)

        @test issuccessful(result)
        scenario = result.value
        @test scenario.tags == ["@tag1", "@tag2"]
    end

    @testset "Scenario Outline Examples" begin
        @testset "Outline has three placeholders; The placeholders are parsed" begin
            text = """
            Scenario Outline: This is one scenario outline
                Given a precondition with placeholders <Foo>, <Bar>, <Baz>

            Examples:
                | Foo | Bar | Baz |
                | 1   | 2   | 3   |
                | 1   | 2   | 3   |
            """
            byline = ByLineParser(text)

            result = parsescenario(byline)

            @test issuccessful(result)
            scenario = result.value
            @test scenario.placeholders == ["Foo", "Bar", "Baz"]
        end

        @testset "Two examples with three placeholders are provided; Examples array is 3x2" begin
            text = """
            Scenario Outline: This is one scenario outline
                Given a precondition with placeholders <Foo>, <Bar>, <Baz>

            Examples:
                | Foo | Bar | Baz |
                | 1   | 2   | 3   |
                | 1   | 2   | 3   |
            """
            byline = ByLineParser(text)

            result = parsescenario(byline)

            @test issuccessful(result)
            scenario = result.value
            @test size(scenario.examples) == (3,2)
        end

        @testset "Three examples with four placeholders are provided; Examples array is 4x3" begin
            text = """
            Scenario Outline: This is one scenario outline
                Given a precondition with placeholders <Foo>, <Bar>, <Baz>, <Quux>

            Examples:
                | Foo | Bar | Baz | Quux |
                | 1   | 2   | 3   | 4    |
                | 1   | 2   | 3   | 4    |
                | 1   | 2   | 3   | 4    |
            """
            byline = ByLineParser(text)

            result = parsescenario(byline)

            @test issuccessful(result)
            scenario = result.value
            @test size(scenario.examples) == (4,3)
        end

        @testset "Two examples with three placeholders are provided; Examples array has all values" begin
            text = """
            Scenario Outline: This is one scenario outline
                Given a precondition with placeholders <Foo>, <Bar>, <Baz>

            Examples:
                | Foo | Bar | Baz |
                | 1   | 2   | 3   |
                | 4   | 5   | 6   |
            """
            byline = ByLineParser(text)

            result = parsescenario(byline)

            @test issuccessful(result)
            scenario = result.value
            @test scenario.examples[:,1] == ["1", "2", "3"]
            @test scenario.examples[:,2] == ["4", "5", "6"]
        end

        @testset "Examples with spaces; Examples are split on | not on spaces" begin
            text = """
            Scenario Outline: This is one scenario outline
                Given a precondition with placeholder <Foo>

            Examples:
                | Foo       |
                | word      |
                | two words |
            """
            byline = ByLineParser(text)

            result = parsescenario(byline)

            @test issuccessful(result)
            scenario = result.value
            @test scenario.examples[:,1] == ["word"]
            @test scenario.examples[:,2] == ["two words"]
        end
    end

    @testset "Malformed Scenario Outlines" begin
        @testset "Invalid step definition; Parse is unsuccessful" begin
            text = """
            Scenario Outline: This is one scenario outline
                NotAStepDefinition with placeholders <Foo>, <Bar>, <Baz>

            Examples:
                | Foo | Bar | Baz |
                | 1   | 2   | 3   |
                | 4   | 5   | 6   |
            """
            byline = ByLineParser(text)

            result = parsescenario(byline)

            @test !issuccessful(result)
        end
    end
end


text:
function mvuniform(offset::Integer, lb::Real, ub::Real, nrows::Integer, ncols::Integer)
    iid(Float64, i->uniform(i, lb,ub), nrows, ncols, offset = offset)
end

"Multivariate Uniform Vector"
mvuniform(lb::Real, ub::Real, n::Integer) = RandArray(RandVar{Float64}[uniform(lb,ub) for i = 1:n])
"Multivariate Uniform Matrix"
mvuniform(lb::Real, ub::Real, nrows::Integer, ncols::Integer) =
  RandArray(RandVar{Float64}[uniform(lb,ub) for i = 1:nrows,j=1:ncols])

"Multivariate Logistic Vector"
mvlogistic(μ, s, n::Integer) = RandArray(RandVar{Float64}[logistic(μ, s) for i = 1:n])
"Multivariate Logistic Matrix"
mvlogistic(μ, s, nrows::Integer, ncols::Integer) =
  RandArray(RandVar{Float64}[logistic(μ, s) for i = 1:nrows,j=1:ncols])

"Multivariate Exponential Vector"
mvexponential(λ, n::Integer) = RandArray(RandVar{Float64}[exponential(λ) for i = 1:n])

"Multivariate Exponential Matrix"
mvexponential(λ, nrows::Integer, ncols::Integer) =
  RandArray(RandVar{Float64}[exponential(λ) for i = 1:nrows,j=1:ncols])

"Multivariate Normal Vector"
mvnormal(μ, σ, nrows::Integer) =
  RandArray(RandVar{Float64}[normal(μ, σ) for i = 1:nrows])

"Multivariate Normal Matrix"
mvnormal(μ, σ, nrows::Integer, ncols::Integer) =
  RandArray(RandVar{Float64}[normal(μ, σ) for i = 1:nrows, j=1:ncols])

## Independent RandVars
## ====================
function iid(T::DataType, c::Function,
             nrows::Int64, ncols::Int64; offset::Int64 = 0)
  a = RandVar{T}[c(i+(j-1)*(nrows) + offset) for i = 1:nrows, j = 1:ncols]
  RandArray{T,2}(a)
end

function iid(T::DataType, c::Function, nrows::Int64; offset = 0)
  a = RandVar{T}[c(i + offset) for i = 1:nrows]
  RandArray{T,1}(a)
end


text:
import Base: run

"""
    Action(func, args...; kwargs...)

An `Action` is a structure (a functor in fact) which stores function, arguments and keyword arguments.

An `Action` can be run (in fact it's run internally by a scheduler when a `Job` is triggered.)
"""
struct Action
    func::Function
    args
    kwargs

    Action(func, args...; kwargs...) = new(func, args, kwargs)
end

"""
    run(action::Action)

Run `action`.

This function shouldn't be called directly. It's called by scheduler when a job is triggered.
"""
run(action::Action) = action.func(action.args...; action.kwargs...)


text:
#############################################################################
#############################################################################
#
# This file implements polynomial multiplication
#
#############################################################################
#############################################################################
using Primes

"""
Multiply two polynomials.
"""

#Multiplication in Z_p[x]

function *(p1::PolynomialModP, p2::PolynomialModP) #Fastest Method

    @assert p1.mod == p2.mod
    prime = p1.mod
    iszero(p1) && return zero(PolynomialModP,prime)
    iszero(p2) && return zero(PolynomialModP,prime)

    degree(p1) == 0 && return p2*leading(p1).coeff
    degree(p2) == 0 && return p1*leading(p2).coeff

    length(p1) == 0 && return p1.polynomial.terms[1]*p2
    length(p2) == 0 && return p2.polynomial.terms[1]*p1

    #Converts Polynomials into vector of coefficients
    p1 = coeffvector(p1.polynomial)
    p2 = coeffvector(p2.polynomial)

    n1 = length(p1)
    n2 = length(p2)
    output = zeros(n1+n2)
    outpoly = Term[]

    #Expands terms like with old multiplication method
    for i in 1:n1
        if p1[i] == 0
            i = i
        else
            for j in 1:n2
                if p2[j] == 0
                    j = j
                else
                    output[i+j] = mod(output[i+j] + mod(p1[i]*p2[j],prime),prime)
                end
            end
        end
    end
    output = Int.(output)

    #converting back to Polynomial
    for i in 1:n1+n2
        push!(outpoly, Term(output[i],(n1+n2)-i))
    end

    return PolynomialModP(outpoly,prime)

end

function mult(p1::PolynomialModP, p2::PolynomialModP)::PolynomialModP #Old Method
    @assert p1.mod == p2.mod
    p = p1.mod
    p_out = PolynomialModP(p)
    for t in p1
        p_out = p_out + (t * p2)
    end
    return p_out
end

function mult2(p1::PolynomialModP, p2::PolynomialModP) #Splitting method V2 (Demonstated in lectures)
    @assert p1.mod == p2.mod
    prime = p1.mod

    (p1 == 0 || p2 == 0) && return zero(prime)

    n1 = length(p1)
    n2 = length(p2)
    p1pairs = [(0,0) for _ in 1:n1]
    p2pairs = [(0,0) for _ in 1:n2]
    i = 1
    j = 1
    for t in p1.polynomial.terms
        p1pairs[i] = (t.coeff,t.degree)
        i = i+1
    end
    for t in p2.polynomial.terms
        p2pairs[j] = (t.coeff,t.degree)
        j = j+1
    end

    #Converts the polynomials to vetcors of tuples

    return split(p1pairs,p2pairs,prime)

end

function mult3(p1::PolynomialModP, p2::PolynomialModP) #splitting method V1, same as mult 3 except with polynomials instead of vectors
    @assert p1.mod == p2.mod
    prime = p1.mod

    (p1 == 0 || p2 == 0) && return zero(prime)
    n1 = length(p1.polynomial.terms)
    n2 = length(p2.polynomial.terms)

    (n1 == 1 && n2 == 1) && return PolynomialModP(leading(p1)*leading(p2),prime)

    n1 == 1 && return p2*leading(p1)
    n2 == 1 && return p1*leading(p2)

    p1a = PolynomialModP(Polynomial(p1.polynomial.terms[1:Int(floor(n1/2))],true),prime,true)
    p1b = PolynomialModP(Polynomial(p1.polynomial.terms[Int(floor(n1/2+1)):n1],true),prime,true)
    p2a = PolynomialModP(Polynomial(p2.polynomial.terms[1:Int(floor(n2/2))],true),prime,true)
    p2b = PolynomialModP(Polynomial(p2.polynomial.terms[Int(floor(n2/2+1)):n2],true),prime,true)

    return (p1a*p2a+p1b*p2b+p1b*p2a+p2b*p1a)

end

# Multiplication in Z[x]
height(p::Polynomial)::Int = maximum(abs.(coeff.(p.terms)))

nonzero(p::Polynomial)::Int = length(p.terms) #number of non-zero terms (the # operator)

required_M(p1::Polynomial, p2::Polynomial)::Int = 2 * height(p1) * height(p2) * minimum([nonzero(p1),nonzero(p2)]) + 1
#An upper bound on the product of primes used for CRT, M

function fullcrt(u::Vector{Int},m::Vector{Int}) #An implimentation of the full CRT algorthm
    n = length(u)
    v = Array{Int}(undef, n)
    v[1] = u[1]
    out = v[1]
    if n > 1
        v[2] = ((u[2]- v[1])*int_inverse_mod(m[1],m[2]))%m[2]
        modprod = m[1]
        a = v[1]
        out  = out+v[2]*modprod
        if n > 2
            for i in 2:n-1
                a = a + v[i]*modprod
                modprod = modprod*m[i]
                v[i+1] = ((u[i+1] - a)*int_inverse_mod(modprod,m[i+1]))%m[i+1]
                out  = out+v[i+1]*modprod
            end
        end
    end
    return symmod(out,modprod*last(m))
end

function crt(u::Vector{Int},m::Vector{Int}) #An implimentation of CRT for just two cases
    v = ((u[2]- u[1])*int_inverse_mod(m[1],m[2]))%m[2]
    out  = u[1]+v*m[1]
    return symmod(out,m[1]*m[2])
end

function split(a::Vector{Tuple{Int64, Int64}}, b::Vector{Tuple{Int64, Int64}}, p::Int)
    n1 = length(a)
    n2 = length(b)
    #used for mult2  and 3
    #Impliments the splitting multiplication demonstated in lecturs

    if (n1 == 1 && n2 == 1) #these are mostly edge cases

        return PolynomialModP([Term(mod(a[1][1]*b[1][1],p),a[1][2]+b[1][2])],p)

    elseif (n1 == 0 || n2 == 0)

         return zero(PolynomialModP,p)

    elseif (n1 == 1)

        coef = a[1][1]
        deg = a[1][2]
        i = 1
        for t in a
            t = (mod(t[1]*coef,p),t[2]+deg)
        end
        output = Term.(b)
        return PolynomialModP(output,p)

    elseif (n2 == 1)

        coef = b[1][1]
        deg = b[1][2]
        i = 1
        for t in b
            t = ( mod(t[1]*coef,p), t[2]+deg )
        end
        output = Term.(b)
        return PolynomialModP(output,p)

    else #splits vector in half and expands the terms uisng a recusive call
        a1 = a[1:Int(floor(n1/2))]
        a2 = a[Int(floor(n1/2+1)):n1]
        b1 = b[1:Int(floor(n2/2))]
        b2 = b[Int(floor(n2/2+1)):n2]

        return split(a1,b1,p) + split(a2,b2,p) + split(a1,b2,p) + split(a2,b1,p)

    end
end
function *(p1::Polynomial, p2::Polynomial) #CRT multiplication
    iszero(p1)||iszero(p2) && return zero(Polynomial)
    degree(p1) == 0 && return p2*leading(p1).coeff
    degree(p2) == 0 && return p1*leading(p2).coeff
    (degree(p1) + degree(p2) < 6) && (return p1*̄p2) #The old multimplication is more efficien for small degrees
    M = required_M(p1,p2)
    n = degree(p1)+degree(p2)
    prime = 1
    pprod = 3
    nextpolyterms = coeffvector(PolynomialModP(p1,3)*PolynomialModP(p2,3), n+1)

    while pprod < M

        prime = nextprime(pprod+1)
        nextpolyterms = mod.(nextpolyterms,pprod)

        p1mod = PolynomialModP(p1,prime)
        p2mod = PolynomialModP(p2,prime)
        p = p1mod*p2mod

        a = coeffvector(p,n+1)

        for j in 1:n+1
            t = crt([a[j],nextpolyterms[j]],[prime,pprod])
            nextpolyterms[j] = t
        end
        pprod = prime*pprod
    end

    output = Array{Term}(undef,n+1)
    i = n+1

    for t in nextpolyterms
        output[i] = Term(t,i-1)
        i = i-1
    end

    return Polynomial(output)
end

function *̄(p1::Polynomial, p2::Polynomial)::Polynomial #Old Multiplication
    p_out = Polynomial()
    for t in p1
        p_out = p_out + (t * p2)
    end
    return p_out
end

"""
Power of a polynomial.
"""

function pow(p::PolynomialModP, n::Int) #old ^ function,
    n < 0 && error("No negative power")
    out = one(p)
    for _ in 1:n
        out *= p
    end
    return out
end

function ^(p::Polynomial,n::Int) #The new optimised Power Function, which uses
    n == 0 && return one(Polynomial) #the binary expansion of the power to reduce the number of operations needed
    n == 1 && return p
    iszero(p) && return zero(Polynomial)

    a = bitstring(n)
    twopowers = zeros(64)
    maxdegree = 1
    result = one(Polynomial)

    for i in 1:64
        (a[65-i] == '1') && (maxdegree = i)
        twopowers[i] = Int(a[i])-48
    end
    #Makes a  vector that represents the binary expansion of the power

    (twopowers[64] == 1) && (result = result*p)

    for i in 2:maxdegree # loops through and multiplies together each of powers of two corresponding
                        # to the binary expansion
        p = p*p
        (twopowers[65-i] == 1) && (result = result*p)
    end

    return result
end

function ^(p::PolynomialModP,n::Int) #Same as above implimentation accept with mod p
    n == 0 && return one(PolynomialModP,p.mod)
    n == 1 && return p
    iszero(p) && return zero(PolynomialModP,p.mod)
    prime = p.mod
    a = bitstring(n)
    twopowers = zeros(64)
    maxdegree = 1
    for i in 1:64

        if (a[65-i] == '1')
            maxdegree = i

        end
        twopowers[i] = Int(a[i])-48
    end

    result = one(PolynomialModP,prime)

    for i in 1:maxdegree
        (twopowers[65-i] == 1) && (result = result*p)
        p = p*p #everything gets reduced mod prime at this step

    end
    return result
end


text:
using OffsetArrays
using ImageFiltering
using Distances
using PaddedViews

"""Apply the opening operation to `A` with window size `ω`."""
function opening(A::Array{T,2}, ω::Integer) where {T<:Real}
    A = mapwindow(minimum, A, (ω, ω))  # erosion
    A = mapwindow(maximum, A, (ω, ω))  # dilation
    A
end

"""Apply the opening operation to `A` with window size `ω`."""
function opening!(A::Array{T,2}, ω::Integer, out::Array{T,2}) where {T<:Real}
    mapwindow_sep!(minimum, A, ω, out, Inf)  # erosion
    mapwindow_sep!(maximum, out, ω, A, -Inf)  # dilation
    A
end

"""Apply the opening operation to `A` with window size `ω`."""
function opening_circ!(A::Array{T,2}, ω::Integer, out::Array{T,2}) where {T<:Real}
    mapwindowcirc!(minimum_mask, A, ω, out, Inf)  # erosion
    mapwindowcirc!(maximum_mask, out, ω, A, -Inf)  # dilation
    A
end

function circmask(n::Integer)
    # TODO This could be precomputed for first N integers
    kern = falses(-n:n, -n:n)
    for I in CartesianIndices(kern)
        i, j = Tuple(I)
        kern[I] = i^2 + j^2 <= n^2
    end
    return kern.parent
end

function opening_circ(A::Array{T,2}, ω::Integer) where {T<:Real}
    m = circmask(ω >> 1)
    A = mapwindow(x -> minimum(x[m]), A, (ω, ω))  # erosion
    A = mapwindow(x -> maximum(x[m]), A, (ω, ω))  # dilation
    A
end

# First discussed here https://github.com/JuliaImages/ImageFiltering.jl/issues/179
function mapwindow!(f, img, window, out)
    R = CartesianIndices(img)
    I_first, I_last = first(R), last(R)
    Δ = CartesianIndex(ntuple(x -> window ÷ 2, ndims(img)))
    @inbounds for I in R
        patch = max(I_first, I - Δ):min(I_last, I + Δ)
        out[I] = f(view(img, patch))
    end
    out
end

function mapwindow_stack!(f, img, window, out)
    R = CartesianIndices(img)
    I_first, I_last = first(R), last(R)
    Δ = CartesianIndex(1, 1)
    out2 = copy(img)
    iterations = window:-2:3
    @inbounds for _ in iterations  # repeat 3x3 window
        for I in R
            patch = max(I_first, I - Δ):min(I_last, I + Δ)
            out[I] = f(view(out2, patch))
        end
        out2 .= out
    end
    out
end

function mapwindow_sep!(f, img, window, out, fill = Inf)
    Δ = window ÷ 2

    w, h = size(img)
    A = PaddedView(fill, img, (-Δ+1:w+Δ, -Δ+1:h+Δ))
    out2 = copy(out)

    # Maximum/minimum is seperable into 1d
    @inbounds for i in 1:h, j in 1:w
        out2[j, i] = f(@view A[j-Δ:j+Δ, i])
    end
    A = PaddedView(fill, out2, (-Δ+1:w+Δ, -Δ+1:h+Δ))
    @inbounds for j in 1:w, i in 1:h
        out[j, i] = f(@view A[j, i-Δ:i+Δ])
    end
    out
end

function mapwindowcirc!(f, img, window, out, fill = Inf)
    R = CartesianIndices(img)
    Δ = CartesianIndex(ntuple(x -> window ÷ 2, ndims(img)))

    w, h = size(img)
    A = PaddedView(fill, img, (-Δ[1]+1:w+Δ[1], -Δ[2]+1:h+Δ[2]))

    m = euclidean.(Tuple.(-Δ:Δ), Ref((0, 0))) .<= Δ[1]
    @inbounds for I in R
        patch = (I-Δ):(I+Δ)
        out[I] = f(view(A, patch), m)
    end
    out
end

function mapwindowcirc_approx!(f, img, window, out, fill = Inf)
    R = CartesianIndices(img)
    Δ = CartesianIndex(1, 1)

    w, h = size(img)

    iterations = window:-2:3
    A = PaddedView(fill, img, (-Δ[1]+1:w+Δ[1], -Δ[2]+1:h+Δ[2]))

    m = euclidean.(Tuple.(-Δ:Δ), Ref((0, 0))) .<= Δ[1]
    @inbounds for _ in iterations  # repeat 3x3 window
        for I in R
            patch = (I-Δ):(I+Δ)
            out[I] = f(view(A, patch), m)
        end
        img .= out
    end
    out
end

@inline function maximum_mask(x, m)
    o = -Inf
    @inbounds for I in eachindex(x)
        o = ifelse(m[I], max(o, x[I]), o)
    end
    o
end

@inline function minimum_mask(x, m)
    o = Inf
    @inbounds for I in eachindex(x)
        o = ifelse(m[I], min(o, x[I]), o)
    end
    o
end