text:
######################################################################
# Scene Graphs resemble an abstract representation of 2D & 3D scenes.
# SceneGraphs library is designed to plug into existing systems as a
# component.
# -----
# Licensed under Apache License 2.0
module SceneGraphs
using ExtraFun
using GenerateProperties
using StaticArrays
import LinearAlgebra: I
import Base.Threads: @threads

# ===== Typedefs =====

const Optional{T} = Union{T, Nothing}

export Vector2, Vector3, Vector4
const Vector2{T} = SVector{2, T}
const Vector3{T} = SVector{3, T}
const Vector4{T} = SVector{4, T}

export Matrix2, Matrix3, Matrix4
const Matrix2{T} = SMatrix{2, 2, T}
const Matrix3{T} = SMatrix{3, 3, T}
const Matrix4{T} = SMatrix{4, 4, T}

export idmat
idmat(T::Type{<:SMatrix}) = T(I)

# ===== (Abstract) Types =====

export AbstractScene, Scene3D
"""`Scene{E, T}`
 Represents the entire scene. Privided primarily for convenience. It is intended to support cross-dimensional entities
 as well, e.g. 2D entities in a 3D scene. The implementation specifics of the scene are left to the implementer."""
abstract type AbstractScene{E, T<:Real} end

struct Scene3D{E, T<:Real} <: AbstractScene{E, T}
    roots::Set{E}
end
Scene3D{E, T}() where {E, T<:Real} = Scene3D{E, T}(Set{E}())

export Transform
"""`Transform{D, E, T}`
 Represents the composable transformation matrix of a scene node / entity.

 `D` is the transform's dimensionality.

 `E` is the type of the actual scene node. This type must specialize `SceneGraphs.transformof(node::E)`.

 `T` is the underlying numeric type of the transform."""
abstract type Transform{D, E, T<:Real} end

export SpriteTransform
"""`SpriteTransform` is a transform of 2D entities in a 3D scene."""
mutable struct SpriteTransform{E, T<:Real} <: Transform{3, E, T}
    scene::Optional{AbstractScene{E, T}}
    parent::Optional{E}
    children::Vector{E}
    location::Vector3{T}
    rotation::T
    scale::Vector2{T}
    dirty::Bool
    obj2world::Matrix4{T}
    world2obj::Matrix4{T}
end
function SpriteTransform{E, T}(location::Vector3 = Vector3{T}(0, 0, 0), rotation::Number = 0, scale::Vector2 = Vector2{T}(1, 1)) where {E, T}
    SpriteTransform{E, T}(nothing, nothing, Vector(), Vector3{T}(location...), T(rotation), Vector2{T}(scale...), true, idmat(Matrix4{T}), idmat(Matrix4{T}))
end

@generate_properties SpriteTransform begin
    @set location = setlocation!(self, value)
    @set rotation = setrotation!(self, value)
    @set scale    = setscale!(self, value)
end

# ===== Scene Alteration =====

Base.:(∈)(entity::E, scene::AbstractScene{E}) where E = entity ∈ scene.roots
Base.haskey(scene::AbstractScene{E}, entity::E) where E = entity ∈ scene
function Base.push!(scene::AbstractScene{E}, entity::E) where E
    tf = transformof(entity)
    if parentof(tf) !== nothing
        deparent!(tf)
    end
    update_scene!(tf, scene)
    push!(scene.roots, entity)
    scene
end
function Base.delete!(scene::AbstractScene{E}, entity::E) where E
    tf = transformof(entity)
    update_scene!(tf, nothing)
    delete!(scene.roots, entity)
    scene
end

# ===== Matrix Component Setters =====

export setlocation!
setlocation!(tf::Transform, it) = setlocation!(tf, it...)
function setlocation!(tf::Transform{D}, comps::Real...) where D
    setfield!(tf, :dirty, true)
    setfield!(tf, :location, compose_vector(D, tf.location, comps))
    return tf
end

export setrotation!
setrotation!(tf::Transform, it) = setrotation!(tf, it...)
function setrotation!(tf::Transform{D}, comps::Real...) where D
    setfield!(tf, :dirty, true)
    setfield!(tf, :rotation, compose_vector(D, tf.rotation, comps))
    return tf
end
function setrotation!(tf::SpriteTransform, rot::Real)
    setfield!(tf, :dirty, true)
    setfield!(tf, :rotation, rot)
    return tf
end

export setscale!
setscale!(tf::Transform, it) = setscale!(tf, it...)
setscale!(tf::Transform{D}, s::Real) where D = setscale!(tf, (s for _ ∈ 1:D)...)
function setscale!(tf::Transform{D}, comps::Real...) where D
    setfield!(tf, :dirty, true)
    setfield!(tf, :scale, compose_vector(D, tf.scale, comps))
    return tf
end
function setscale!(tf::SpriteTransform{E, T}, scalex::Real, scaley::Real) where {E, T<:Real}
    setfield!(tf, :dirty, true)
    setfield!(tf, :scale, Vector2{T}(scalex, scaley))
    return tf
end

# ===== Matrix Component Getters (object-local & scene-global) =====

export locationof
locationof(entity) = locationof(transformof(entity))
locationof(tf::Transform) = tf.location

export rotationof
rotationof(entity) = rotationof(transformof(entity))
rotationof(tf::Transform) = tf.rotation

export scaleof
scaleof(entity) = scaleof(transformof(entity))
scaleof(tf::Transform) = tf.scale

export worldlocationof
worldlocationof(entity) = worldlocationof(transformof(entity))
worldlocationof(tf::Transform{D}) where D = (mat = obj2world(tf); mat[1:D, lastindex(mat, 2)])

export worldrotationof
worldrotationof(entity) = worldrotationof(transformof(entity))
function worldrotationof(tf::Transform)
    final_rotation = rotationof(tf)
    curr = parentof(tf)
    while curr !== nothing
        final_rotation = final_rotation .+ rotationof(curr)
        curr = parentof(curr)
    end
    return final_rotation
end

export worldscaleof
worldscaleof(entity) = worldscaleof(transformof(entity))
function worldscaleof(tf::Transform)
    final_scale = scaleof(tf)
    curr = parentof(tf)
    while curr !== nothing
        final_scale = final_scale .* scaleof(curr)
        curr = parentof(curr)
    end
    return final_scale
end

# ===== Matrix Component Modifiers =====

export transform!
function transform!(entity; translate = nothing, rotate = nothing, scale = nothing)
    tf = transformof(entity)
    if translate !== nothing
        translate!(tf, translate...)
    end
    if rotate !== nothing
        rotate!(tf, rotate...)
    end
    if scale !== nothing
        scale!(tf, scale...)
    end
    entity
end

export translate!
translate!(tf::Transform, offset) = translate!(tf, offset...)
translate!(tf::Transform{D}, coords::Real...) where D = translate!(tf, compose_vector(D, zeros(D), coords)...)
translate!(tf::Transform{D}, args::Vararg{Real, D}) where D = (tf.location = tf.location .+ args; tf)
translate!(entity, args...) = (translate!(transformof(entity), args...); entity)

export scale!
scale!(tf::Transform, scale) = scale!(tf, scale...)
scale!(tf::Transform{D}, scale::Real...) where D = scale!(tf, compose_vector(D, ones(D), scale)...)
scale!(tf::Transform{D}, scale::Vararg{Real, D}) where D = tf.scale = tf.scale .* scale
scale!(entity, args...) = (scale!(transformof(entity), args...); entity)

scale!(tf::SpriteTransform, scalex::Real, scaley::Real) = tf.scale = tf.scale .* (scalex, scaley)

export rotate!
rotate!(tf::Transform, offset) = rotate!(tf, offset...)
rotate!(tf::Transform{D}, angles::Real...) where D = rotate!(tf, compose_vector(D, zeros(D), angles)...)
rotate!(tf::Transform{D}, angles::Vararg{Real, D}) where D = tf.rotation = tf.rotation .+ angles
rotate!(entity, args...) = (rotate!(transformof(entity), args...); entity)

# Specializations
rotate!(tf::SpriteTransform, rotation::Real) = tf.rotation += rotation

# ===== Scene Graph Alteration =====

export parent!
"""`parent!(child, parent)`
 Assign `child` to its new `parent`. This automatically marks the child as *dirty*. Both `child` and `parent` are
 assumed to have transforms."""
parent!(::Nothing, _) = throw(ArgumentError("child is nothing"))
parent!(_, ::Nothing) = throw(ArgumentError("parent is nothing"))
parent!(child, parent) = (parent!(transformof(child), transformof(parent), child, parent); child)

function parent!(tf_child::Transform, tf_parent::Transform, child, parent)
    if child === parent || tf_child === tf_parent
        throw(ArgumentError("cannot parent an entity to itself"))
    end

    if parentof(tf_child) !== parent
        flagdirty!(child)

        deparent!(tf_child, child)
        push!(tf_parent.children, child)
        tf_child.parent = parent
        update_scene!(tf_child, tf_parent.scene)
    end
    return
end

export deparent!
deparent!(child) = deparent!(transformof(child), child)
deparent!(_, ::Transform) = throw(ArgumentError("missing child entity, only got transform"))
function deparent!(tf::Transform, child)
    if parentof(child) !== nothing
        tf_parent = child |> parentof |> transformof
        if (i = indexof(tf_parent.children, child)) !== nothing
            deleteat!(tf_parent.children, i)
        end

        tf.parent = nothing
        tf.dirty  = true
        update_scene!(tf, nothing)
    end

    return child
end

function update_scene!(tf::Transform, scene::Union{Nothing, AbstractScene})
    tf.dirty = true
    tf.scene = scene
    for child ∈ childrenof(tf)
        update_scene!(transformof(child), scene)
    end
    return tf
end

# ===== Trait Interface =====

export sceneof, parentof, childrenof
sceneof(x)    = transformof(x).scene
parentof(x)   = transformof(x).parent
childrenof(x) = transformof(x).children
childrenof(T::Type, x) = filter(c->isa(c, T), childrenof(x))

function ExtraFun.update!(tf::Transform{3, E, T}, parentmat::Matrix4{T} = idmat(Matrix4{T}), forceupdate::Bool = false) where {E, T}
    if tf.dirty || forceupdate
        tf.obj2world = parentmat * transformmatrix4(tf)
        tf.world2obj = inv(tf.obj2world)

        forceupdate = true
        cleardirty!(tf)
    end

    @threads for child ∈ childrenof(tf)
        update!(transformof(child), tf.obj2world, forceupdate)
    end

    return tf
end

export flagdirty!, cleardirty!
flagdirty!(entity) = (flagdirty!(transformof(entity)); entity)
flagdirty!(tf::Transform) = (tf.dirty = true; tf)
cleardirty!(entity) = (cleardirty!(transformof(entity)); entity)
cleardirty!(tf::Transform) = (tf.dirty = false; tf)

export obj2world, world2obj
obj2world(x) = obj2world(transformof(x))
world2obj(x) = world2obj(transformof(x))
obj2world(transform::Transform) = transform.obj2world
world2obj(transform::Transform) = transform.world2obj

export hastransform
@generated function hastransform(x)
    if hassignature(transformof, x)
        :(true)
    else
        :(false)
    end
end

export transformof, transform_family, transform_dimensionality, transform_entity_type, transform_number_type
transformof(x)::Transform = x.transform
transformof(x::Transform) = x
transform_family(T::Type{<:Transform}) = throw(MethodError(transform_family, (T,)))
transform_family(::Type{<:SpriteTransform}) = SpriteTransform
transform_dimensionality(::Type{<:Transform{D, E, T}}) where {D, E, T} = D
transform_entity_type(::Type{<:Transform{D, E, T}}) where {D, E, T} = E
transform_number_type(::Type{<:Transform{D, E, T}}) where {D, E, T} = T

# ===== Transformation (Matrix) Methods =====

export translationmatrix4, rotationmatrix4, scalematrix4
translationmatrix4(T::Type{<:Real}, x::Real, y::Real, z::Real) = Matrix4{T}(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, x, y, z, 1)
rotationmatrix4(T::Type{<:Real}, p::Real, y::Real, r::Real) = yawmatrix(T, y) * pitchmatrix(T, p) * rollmatrix(T, r)
scalematrix4(T::Type{<:Real}, x::Real, y::Real, z::Real = 1) = Matrix4{T}(x, 0, 0, 0, 0, y, 0, 0, 0, 0, z, 0, 0, 0, 0, 1)
scalematrix4(T::Type{<:Real}, s::Real) = scalematrix4(T, s, s, s)

export pitchmatrix, yawmatrix, rollmatrix
pitchmatrix(T::Type{<:Real}, p::Real) = (sinr = sin(p); cosr = cos(p); Matrix4{T}(1, 0, 0, 0, 0, cosr, sinr, 0, 0, -sinr, cosr, 0, 0, 0, 0, 1))
yawmatrix(  T::Type{<:Real}, y::Real) = (sinr = sin(y); cosr = cos(y); Matrix4{T}(cosr, 0, -sinr, 0, 0, 1, 0, 0, sinr, 0, cosr, 0, 0, 0, 0, 1))
rollmatrix( T::Type{<:Real}, r::Real) = (sinr = sin(r); cosr = cos(r); Matrix4{T}(cosr, sinr, 0, 0, -sinr, cosr, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1))

export transformmatrix4
transformmatrix4(tf::Transform{D, E, T}) where {D, E, T} = translationmatrix4(T, tf.location...) * rotationmatrix4(T, tf.rotation...) * scalematrix4(T, tf.scale...)
transformmatrix4(tf::SpriteTransform{E, T}) where {E, T} = translationmatrix4(T, tf.location...) * rollmatrix(T, tf.rotation) * scalematrix4(T, tf.scale...)

export transform
transform(tf, ::Nothing) = tf
transform(vert::SVector{D}, mat::SMatrix{D, D}) where D = mat * vert
function transform(verts, mat::SMatrix{N, N}) where N
    transformed = [vformat(vert, Vector4(0, 0, 0, 1)) for vert ∈ verts]
    @threads for i ∈ 1:length(transformed)
        transformed[i] = mat * vformat(transformed[i], Vector4(0, 0, 0, 1))
    end
    transformed
end

# === Base Overrides ===
Base.copy(tf::SpriteTransform{E, T}) where {E, T} = SpriteTransform{E, T}(nothing, nothing, [], tf.location, tf.rotation, tf.scale, true, idmat(Matrix4), idmat(Matrix4))

function Base.show(io::IO, transform::Transform)
    write(io, "$(typeof(transform))(loc: $(transform.location), rot: $(transform.rotation), scale: $(transform.scale)")
    if transform.dirty
        write(io, ", dirty")
    end
    write(io, ")")
end

# === Utilities ===
function compose_vector(D::Integer, missing_defaults, comps)
    if length(comps) < D
        SVector{D}(comps..., (missing_defaults[i] for i ∈ length(comps)+1:D)...)
    else
        SVector{D}((comps[i] for i ∈ 1:D)...)
    end
end

end # module SceneGraphs


text:
"""
struct MCMCResults

Struct for holding the MCMC chain (a vector of NamedTuples which represent the accepted parameter sets), and the MCMC tree statistics
"""
struct MCMCResults
    chain
    logL::Vector{Float64}
    treestatistics::Vector{DynamicHMC.TreeStatisticsNUTS}
end

"""
    mutable struct CoVAreaModel

Struct for holding the fixed data of the area and MCMC results.
"""
Base.@kwdef mutable struct CoVAreaModel
    areaname::String #Name of area
    PCR_cases::Matrix{Int64} = zeros(Int64,365,2) #Daily PCR cases: dim 1: day, dim 2: num pos and total
    sero_cases::Array{Int64,3} = zeros(Int64,365,2,3) #Daily serology tests: dim 1: day, dim 2: num pos and neg, dim 3: 1 = Overall, 2 = lower SES sero tests, 3 = higher SES sero tests
    deaths::Vector{Int64} = zeros(Int64,365) #Daily deaths
    dates::Vector{Date} = [Date(2020,2,20) + Day(k) for k = 1:365]#Dates of each day
    N::Float64 = 4.3e6 #Area's population size
    α::Float64 = 1/3.1 # 1/mean latent period
    γ::Float64 = 1/2.4#Recovery rate
    relative_testing_rate::Vector{Float64} = relative_testing_rate
    prob::ODEProblem = build_prob()
    PCR_array::Vector{Float64} = PCR_array #Relative Sensitivity of PCR test on days post-infection
    PCR_sensitivity::Float64 = 1. #Base sensitivity for RT-PCR test
    PCR_specificity::Float64 = 0.995 #Base specificity for RT_PCR test
    baseline_sero_array::Vector{Float64} = rel_sero_array_26days #Relative sensitivity of serology test on days post-infection
    sero_sensitivity::Float64 = 0.825 #Base sensitivity for serology test
    sero_specificity::Float64 = 0.992 #Base specificity for serology test
    M_BB::Float64 = 40.74 #Shrinkage factor for BetaBinomial distribution that incorporates the uncertainty in the sensitivity (fitted to serology group)
    log_likelihood::Function = θ -> 0.
    log_priors::Function = θ -> 0.
    MCMC_results::Union{MCMCResults,Nothing} = nothing #This field gets added after the MCMC run
end

"""
    struct CoVCountryModel

Struct for holding a collection of `CovAreaModel` objects, fitting to them simultaneously, and a sero-waning rate that covers all of them.
"""
Base.@kwdef struct CoVCountryModel
    name::String = "Kenya"
    areas::Vector{CoVAreaModel}
    p_serowaning::Float64 = 1/180.
    log_priors::Function = θ -> 0.
end

"""
    CoVAreaModel(name::String,loglk_func,log_prior_func;
                        case_data,sero_data,death_data,pop_data::NamedArray)

Constructor for a `CoVAreaModel` struct using country data objects.
"""
function CoVAreaModel(name::String,loglk_func,log_prior_func;
                        case_data,sero_data,death_data,pop_data::NamedArray)
        uprname = uppercase(name)
        PCR_cases = sum(case_data.cases[:,case_data.areas .== uprname,:,:],dims = [2,3])[:,1,1,:]
        sero_cases = sum(sero_data.serodata[:,sero_data.areas .== uprname,:,:],dims = [2,3])[:,1,1,:]
        sero_cases = cat(sero_cases,zeros(Int64,size(sero_cases)),zeros(Int64,size(sero_cases)),dims = 3)
        deaths = death_data.deaths[:,death_data.areas .== uprname][:]
        dates = [Date(2020,2,20) + Day(k) for k = 1:size(PCR_cases,1)]
        N = sum(pop_data[:,name])

        return KenyaCoVSD.CoVAreaModel(areaname = name,
                                        PCR_cases=PCR_cases,
                                        sero_cases=sero_cases,
                                        deaths=deaths,
                                        dates = dates,
                                        N=N,
                                        log_likelihood=loglk_func,
                                        log_priors = log_prior_func)

end

function (model::KenyaCoVSD.CoVAreaModel)(θ)
    model.log_likelihood(θ,model,1/180) + model.log_priors(θ)
end

"""
    function modeldic(areamodel::CoVAreaModel)

Calculate the deviance information criterion (DIC) for the posterior parameter draws in `areamodel`.
"""
function modeldic(areamodel::CoVAreaModel)
        D = -2*areamodel.MCMC_results.logL
        D̄ = mean(D)
        p_D = 0.5*var(D)
        return D̄ + p_D
end

ext: jl


text:
mutable struct TranslationVar{T,N} <: AbstractExtDataAdv
    values::NTuple{N,T}
    valok
end
function gettranslationvar(v::NTuple{N,T}) where {T,N}
    return TranslationVar{T,N}(v, ntuple(x->zero(T),N))
end

function initcoef!(pv::TranslationVar{T,N}, self::AdvectionData{T,N},) where {T,N}
    st = getst(self)
    pv.valok = ntuple( i -> pv.values[st.perm[i]] * getcur_t(self), st.ndims)
end

"""
    getalpha(pv::TranslationVar, self::AdvectionData, i, ind)

Implementation of the interface function that is called before each interpolation in advection

"""
function getalpha(
    pv::TranslationVar{T,N},
    self::AdvectionData{T,N},
    ind,
) where {T,N}
    return pv.valok
end


text:
using FromFile
using Distributed
using LossFunctions
#TODO - eventually move some of these
# into the SR call itself, rather than
# passing huge options at once.
@from "Operators.jl" import plus, pow, mult, sub, div, log_abs, log10_abs, log2_abs, log1p_abs, sqrt_abs, acosh_abs, atanh_clip
@from "Equation.jl" import Node

"""
         build_constraints(una_constraints, bin_constraints,
                           unary_operators, binary_operators)

Build constraints on operator-level complexity from a user-passed dict.
"""
function build_constraints(una_constraints, bin_constraints,
                           unary_operators, binary_operators,
                           nuna, nbin)::Tuple{Array{Int, 1}, Array{Tuple{Int,Int}, 1}}
    # Expect format ((*)=>(-1, 3)), etc.
    # TODO: Need to disable simplification if (*, -, +, /) are constrained?
    #  Or, just quit simplification is constraints violated.

    is_bin_constraints_already_done = typeof(bin_constraints) <: Array{Tuple{Int,Int},1}
    is_una_constraints_already_done = typeof(una_constraints) <: Array{Int,1}

    if typeof(bin_constraints) <: Array && !is_bin_constraints_already_done
        bin_constraints = Dict(bin_constraints)
    end
    if typeof(una_constraints) <: Array && !is_una_constraints_already_done
        una_constraints = Dict(una_constraints)
    end

    if una_constraints == nothing
        una_constraints = [-1 for i=1:nuna]
    elseif !is_una_constraints_already_done
        una_constraints::Dict
        _una_constraints = Int[]
        for (i, op) in enumerate(unary_operators)
            did_user_declare_constraints = haskey(una_constraints, op)
            if did_user_declare_constraints
                constraint::Int = una_constraints[op]
                push!(_una_constraints, constraint)
            else
                push!(_una_constraints, -1)
            end
        end
        una_constraints = _una_constraints
    end
    if bin_constraints == nothing
        bin_constraints = [(-1, -1) for i=1:nbin]
    elseif !is_bin_constraints_already_done
        bin_constraints::Dict
        _bin_constraints = Tuple{Int,Int}[]
        for (i, op) in enumerate(binary_operators)
            did_user_declare_constraints = haskey(bin_constraints, op)
            if did_user_declare_constraints
                constraint::Tuple{Int,Int} = bin_constraints[op]
                push!(_bin_constraints, constraint)
            else
                push!(_bin_constraints, (-1, -1))
            end
        end
        bin_constraints = _bin_constraints
    end

    return una_constraints, bin_constraints
end

function binopmap(op)
    if op == plus
        return +
    elseif op == mult
        return *
    elseif op == sub
        return -
    elseif op == div
        return /
    elseif op == ^
        return pow
    end
    return op
end

function unaopmap(op)
    if op == log
        return log_abs
    elseif op == log10
        return log10_abs
    elseif op == log2
        return log2_abs
    elseif op == log1p
        return log1p_abs
    elseif op == sqrt
        return sqrt_abs
    elseif op == acosh
        return acosh_abs
    elseif op == atanh
        return atanh_clip
    end
    return op
end

struct Options{A,B,C<:Union{SupervisedLoss,Function}}

    binops::A
    unaops::B
    bin_constraints::Array{Tuple{Int,Int}, 1}
    una_constraints::Array{Int, 1}
    ns::Int
    parsimony::Float32
    alpha::Float32
    maxsize::Int
    maxdepth::Int
    fast_cycle::Bool
    migration::Bool
    hofMigration::Bool
    fractionReplacedHof::Float32
    shouldOptimizeConstants::Bool
    hofFile::String
    npopulations::Int
    perturbationFactor::Float32
    annealing::Bool
    batching::Bool
    batchSize::Int
    mutationWeights::Array{Float64, 1}
    warmupMaxsizeBy::Float32
    useFrequency::Bool
    npop::Int
    ncyclesperiteration::Int
    fractionReplaced::Float32
    topn::Int
    verbosity::Int
    probNegate::Float32
    nuna::Int
    nbin::Int
    seed::Union{Int, Nothing}
    loss::C
    progress::Bool
    terminal_width::Union{Int, Nothing}
    optimizer_algorithm::String
    optimize_probability::Float32
    optimizer_nrestarts::Int
    optimizer_iterations::Int

end

"""
    Options(;kws...)

Construct options for `EquationSearch` and other functions.

# Arguments
- `binary_operators=(div, plus, mult)`: Tuple of binary
    operators to use. Each operator should be defined for two input scalars,
    and one output scalar. All operators need to be defined over the entire
    real line (excluding infinity - these are stopped before they are input).
    Thus, `log` should be replaced with `log_abs`, etc.
    For speed, define it so it takes two reals
    of the same type as input, and outputs the same type. For the SymbolicUtils
    simplification backend, you will need to define a generic method of the
    operator so it takes arbitrary types.
- `unary_operators=(exp, cos)`: Same, but for
    unary operators (one input scalar, gives an output scalar).
- `constraints=nothing`: Array of pairs specifying size constraints
    for each operator. The constraints for a binary operator should be a 2-tuple
    (e.g., `(-1, -1)`) and the constraints for a unary operator should be an `Int`.
    A size constraint is a limit to the size of the subtree
    in each argument of an operator. e.g., `[(^)=>(-1, 3)]` means that the
    `^` operator can have arbitrary size (`-1`) in its left argument,
    but a maximum size of `3` in its right argument. Default is
    no constraints.
- `batching=false`: Whether to evolve based on small mini-batches of data,
    rather than the entire dataset.
- `batchSize=50`: What batch size to use if using batching.
- `loss=L2DistLoss()`: What loss function to use. Can be one of
    the following losses, or any other loss of type
    `SupervisedLoss`. You can also pass a function that takes
    a scalar target (left argument), and scalar predicted (right
    argument), and returns a scalar. This will be averaged
    over the predicted data. If weights are supplied, your
    function should take a third argument for the weight scalar.
    Included losses:
        Regression:
            - `LPDistLoss{P}()`,
            - `L1DistLoss()`,
            - `L2DistLoss()` (mean square),
            - `LogitDistLoss()`,
            - `HuberLoss(d)`,
            - `L1EpsilonInsLoss(ϵ)`,
            - `L2EpsilonInsLoss(ϵ)`,
            - `PeriodicLoss(c)`,
            - `QuantileLoss(τ)`,
        Classification:
            - `ZeroOneLoss()`,
            - `PerceptronLoss()`,
            - `L1HingeLoss()`,
            - `SmoothedL1HingeLoss(γ)`,
            - `ModifiedHuberLoss()`,
            - `L2MarginLoss()`,
            - `ExpLoss()`,
            - `SigmoidLoss()`,
            - `DWDMarginLoss(q)`.
- `npopulations=nothing`: How many populations of equations to use. By default
    this is set equal to the number of cores
- `npop=1000`: How many equations in each population.
- `ncyclesperiteration=300`: How many generations to consider per iteration.
- `ns=10`: Number of equations in each subsample during regularized evolution.
- `topn=10`: Number of equations to return to the host process, and to
    consider for the hall of fame.
- `alpha=0.100000f0`: The probability of accepting an equation mutation
    during regularized evolution is given by exp(-delta_loss/(alpha * T)),
    where T goes from 1 to 0. Thus, alpha=infinite is the same as no annealing.
- `maxsize=20`: Maximum size of equations during the search.
- `maxdepth=nothing`: Maximum depth of equations during the search, by default
    this is set equal to the maxsize.
- `parsimony=0.000100f0`: A multiplicative factor for how much complexity is
    punished.
- `useFrequency=false`: Whether to use a parsimony that adapts to the
    relative proportion of equations at each complexity; this will
    ensure that there are a balanced number of equations considered
    for every complexity.
- `fast_cycle=false`: Whether to thread over subsamples of equations during
    regularized evolution. Slightly improves performance, but is a different
    algorithm.
- `migration=true`: Whether to migrate equations between processes.
- `hofMigration=true`: Whether to migrate equations from the hall of fame
    to processes.
- `fractionReplaced=0.1f0`: What fraction of each population to replace with
    migrated equations at the end of each cycle.
- `fractionReplacedHof=0.1f0`: What fraction to replace with hall of fame
    equations at the end of each cycle.
- `shouldOptimizeConstants=true`: Whether to use NelderMead optimization
    to periodically optimize constants in equations.
- `optimizer_nrestarts=3`: How many different random starting positions to consider
    when using NelderMead optimization.
- `hofFile=nothing`: What file to store equations to, as a backup.
- `perturbationFactor=1.000000f0`: When mutating a constant, either
    multiply or divide by (1+perturbationFactor)^(rand()+1).
- `probNegate=0.01f0`: Probability of negating a constant in the equation
    when mutating it.
- `mutationWeights=[10.000000, 1.000000, 1.000000, 3.000000, 3.000000, 0.010000, 1.000000, 1.000000]`:
- `annealing=true`: Whether to use simulated annealing.
- `warmupMaxsize=0`: Whether to slowly increase the max size from 5 up to
    `maxsize`. If nonzero, specifies how many cycles (populations*iterations)
    before increasing by 1.
- `verbosity=convert(Int, 1e9)`: Whether to print debugging statements or
    not.
- `bin_constraints=nothing`:
- `una_constraints=nothing`:
- `seed=nothing`: What random seed to use. `nothing` uses no seed.
- `progress=false`: Whether to use a progress bar output (`verbosity` will
    have no effect).
"""
function Options(;
    binary_operators::NTuple{nbin, Any}=(div, plus, mult),
    unary_operators::NTuple{nuna, Any}=(exp, cos),
    constraints=nothing,
    loss=L2DistLoss(),
    ns=10, #1 sampled from every ns per mutation
    topn=10, #samples to return per population
    parsimony=0.000100f0,
    alpha=0.100000f0,
    maxsize=20,
    maxdepth=nothing,
    fast_cycle=false,
    migration=true,
    hofMigration=true,
    fractionReplacedHof=0.1f0,
    shouldOptimizeConstants=true,
    hofFile=nothing,
    npopulations=nothing,
    perturbationFactor=1.000000f0,
    annealing=true,
    batching=false,
    batchSize=50,
    mutationWeights=[10.000000, 1.000000, 1.000000, 3.000000, 3.000000, 0.010000, 1.000000, 1.000000],
    warmupMaxsizeBy=0f0,
    useFrequency=false,
    npop=1000,
    ncyclesperiteration=300,
    fractionReplaced=0.1f0,
    verbosity=convert(Int, 1e9),
    probNegate=0.01f0,
    seed=nothing,
    bin_constraints=nothing,
    una_constraints=nothing,
    progress=false,
    terminal_width=nothing,
    warmupMaxsize=nothing,
    optimizer_algorithm="NelderMead",
    optimizer_nrestarts=3,
    optimize_probability=0.1f0,
    optimizer_iterations=100,
    nrestarts=nothing,
   ) where {nuna,nbin}

    if nrestarts != nothing
        optimizer_nrestarts = nrestarts
    end
    if warmupMaxsize != nothing
        error("warmupMaxsize is deprecated. Please use warmupMaxsizeBy, and give the time at which the warmup will end as a fraction of the total search cycles.")
    end

    if hofFile == nothing
        hofFile = "hall_of_fame.csv" #TODO - put in date/time string here
    end

    @assert maxsize > 3
    @assert warmupMaxsizeBy >= 0f0

    constraints::Union{Tuple,Array{Pair{Any,Any}, 1},Nothing}

    if typeof(constraints) <: Tuple
        constraints = collect(constraints)
    end
    if constraints !== nothing
        @assert bin_constraints == nothing
        @assert una_constraints == nothing
        # TODO: This is redundant with the checks in EquationSearch
        for op in binary_operators
            @assert !(op in unary_operators)
        end
        for op in unary_operators
            @assert !(op in binary_operators)
        end
        bin_constraints = constraints
        una_constraints = constraints
    end

    una_constraints, bin_constraints = build_constraints(una_constraints, bin_constraints,
                                                         unary_operators, binary_operators,
                                                         nuna, nbin)

    if maxdepth == nothing
        maxdepth = maxsize
    end

    if npopulations == nothing
        npopulations = nworkers()
    end

    binary_operators = map(binopmap, binary_operators)
    unary_operators = map(unaopmap, unary_operators)

    mutationWeights = map((x,)->convert(Float64, x), mutationWeights)
    if length(mutationWeights) != 8
        error("Not the right number of mutation probabilities given")
    end

    for (op, f) in enumerate(map(Symbol, binary_operators))
        _f = if f == Symbol(pow)
            Symbol(^)
        else
            f
        end
        if !isdefined(Base, _f)
            continue
        end
        @eval begin
            Base.$_f(l::Node, r::Node)::Node = (l.constant && r.constant) ? Node($f(l.val, r.val)::AbstractFloat) : Node($op, l, r)
            Base.$_f(l::Node, r::AbstractFloat)::Node =        l.constant ? Node($f(l.val, r)::AbstractFloat)     : Node($op, l, r)
            Base.$_f(l::AbstractFloat, r::Node)::Node =        r.constant ? Node($f(l, r.val)::AbstractFloat)     : Node($op, l, r)
        end
    end

    for (op, f) in enumerate(map(Symbol, unary_operators))
        if !isdefined(Base, f)
            continue
        end
        @eval begin
            Base.$f(l::Node)::Node = l.constant ? Node($f(l.val)::AbstractFloat) : Node($op, l)
        end
    end

    if progress
        verbosity = 0
    end

    Options{typeof(binary_operators),typeof(unary_operators), typeof(loss)}(binary_operators, unary_operators, bin_constraints, una_constraints, ns, parsimony, alpha, maxsize, maxdepth, fast_cycle, migration, hofMigration, fractionReplacedHof, shouldOptimizeConstants, hofFile, npopulations, perturbationFactor, annealing, batching, batchSize, mutationWeights, warmupMaxsizeBy, useFrequency, npop, ncyclesperiteration, fractionReplaced, topn, verbosity, probNegate, nuna, nbin, seed, loss, progress, terminal_width, optimizer_algorithm, optimize_probability, optimizer_nrestarts, optimizer_iterations)
end


text:
# Copyright (c) 2022 Code Komali
#
# This software is released under the MIT License.
# https://opensource.org/licenses/MIT

using IndianTrafficSim

fig = plot_environment!()

fig


text:
# This file was generated, do not modify it. # hide
using MLJ, DataFrames, Random
MLJ.color_off() # hide
@load RidgeRegressor pkg=MultivariateStats

Random.seed!(5) # for reproducibility
x1 = rand(300)
x2 = rand(300)
x3 = rand(300)
y = exp.(x1 - x2 -2x3 + 0.1*rand(300))
X = DataFrame(x1=x1, x2=x2, x3=x3)

test, train = partition(eachindex(y), 0.8);

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
export FacElem, FacElemMon, simplify, factor_coprime

export transform

################################################################################
#
#  Multiplicative representation
#
################################################################################

# abstract nonsense

const FacElemMonDict = IdDict()

function (x::FacElemMon{S})() where S
  z = FacElem{elem_type(S), S}()
  z.fac = Dict{elem_type(S), fmpz}()
  z.parent = x
  return z
end
@doc Markdown.doc"""
    FacElem{B}(base::Array{B, 1}, exp::Array{fmpz, 1}) -> FacElem{B}
Returns the element $\prod b_i^{e_i}$, un-expanded.
"""
function FacElem(base::Array{B, 1}, exp::Array{fmpz, 1}) where B

  length(base) == 0 && error("Array must not be empty")

  length(base) != length(exp) && error("not the same length")

  z = FacElem{B, typeof(parent(base[1]))}()

  for i in 1:length(base)
    if exp[i] == 0
      continue
    end
    if haskey(z.fac, base[i])
      z.fac[base[i]] += exp[i]
    else
      z.fac[base[i]] = exp[i]
    end
  end

  z.parent = FacElemMon(parent(base[1]))
  return z
end

@doc Markdown.doc"""
    FacElem{B}(d::Dict{B, fmpz}) -> FacElem{B}
    FacElem{B}(d::Dict{B, Integer}) -> FacElem{B}
Returns the element $\prod b^{d[p]}$, un-expanded.
"""
function FacElem(d::Dict{B, fmpz}) where B

  length(d) == 0 && error("Dictionary must not be empty")

  z = FacElem{B, typeof(parent(first(keys(d))))}()
  z.fac = d

  z.parent = FacElemMon(parent(first(keys(d))))
  return z
end

function FacElem(R::Ring)
  z = FacElem{elem_type(R), typeof(R)}()
  z.fac = Dict{elem_type(R), fmpz}()
  z.parent = FacElemMon(R)
  return z
end

function FacElem(d::Dict{B, T}) where {B, T <: Integer}

  length(d) == 0 && error("Dictionary must not be empty")

  z = FacElem{B, typeof(parent(first(keys(d))))}()
  z.fac = Dict{B, fmpz}((k,fmpz(v)) for (k,v) = d)

  z.parent = FacElemMon(parent(first(keys(d))))
  return z
end

parent(x::FacElem) = x.parent

base_ring(x::FacElemMon) = x.base_ring

base_ring(x::FacElem) = base_ring(parent(x))

base(x::FacElem) = keys(x.fac)

function Base.deepcopy_internal(x::FacElem{B, S}, dict::IdDict) where {B, S}
  z = FacElem{B, S}()
  z.fac = Base.deepcopy_internal(x.fac, dict)
  if isdefined(x, :parent)
    z.parent = x.parent
  end
#  z.hash = x.hash # this needs to be deleted in ! operations
  return z
end

function Base.copy(x::FacElem{B, S}) where {B, S}
  z = FacElem{B, S}()
  z.fac = Base.copy(x.fac)
  if isdefined(x, :parent)
    z.parent = x.parent
  end
#  z.hash = x.hash
  return z
end

Base.iterate(a::FacElem) = Base.iterate(a.fac)

Base.iterate(a::FacElem, state) = Base.iterate(a.fac, state)

Base.length(a::FacElem) = Base.length(a.fac)

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, x::FacElemMon)
  print(io, "Factored elements over $(x.base_ring)")
end

function show(io::IO, x::FacElem)
  print(io, "Factored element with data\n$(x.fac)")
end

################################################################################
#
#  Inverse
#
################################################################################
function inv!(x::FacElem)
  for a in base(x)
    x.fac[a] = -x.fac[a]
  end
end

function inv(x::FacElem)
  y = copy(x)
  inv!(y)
  return y
end

################################################################################
#
#  Exponentiation
#
################################################################################

function pow!(z::FacElem, x::FacElem, y::T) where T <: Union{fmpz, Integer}
  z.fac = copy(x.fac)
  for a in base(x)
    # this should be inplace ... not sure anymore: using copy, inplace is bad
    z.fac[a] = y*x.fac[a]
  end
end

function pow!(x::FacElem, y::T) where T <: Union{fmpz, Integer}
  for a in base(x)
    # this should be inplace ... not sure anymore: using copy, inplace is bad
    x.fac[a] = y*x.fac[a]
  end
end

# ^(x::FacElem, y::Union{fmpz, Integer}) is ambiguous
for T in [:Integer, fmpz]
  @eval begin
    function ^(x::FacElem, y::($T))
      z = parent(x)()
      if y == 0
        return z
      end
      if y == 1
        return copy(x)
      else
        z.fac = copy(x.fac)
        for a in base(x)
          # this should be inplac
          z.fac[a] = y*x.fac[a]
        end
        return z
      end
    end
  end
end

################################################################################
#
#  Multiplication
#
################################################################################

function mul!(z::FacElem{B, S}, x::FacElem{B, S}, y::FacElem{B, S}) where {B, S}
  z.fac = copy(x.fac)
  for a in base(y)
    if haskey(x.fac, a)
      z.fac[a] = z.fac[a] + y.fac[a]
    else
      z.fac[a] = y.fac[a]
    end
  end
end

function *(x::FacElem{B, S}, y::FacElem{B, S}) where {B, S}
  if length(x.fac) == 0
    return copy(y)
  end

  if length(y.fac) == 0
    return copy(x)
  end

  z = copy(x)
  for a in base(y)
    if haskey(x.fac, a)
      z.fac[a] = z.fac[a] + y.fac[a]
    else
      z.fac[a] = y.fac[a]
    end
  end
  return z
end

function *(x::FacElem{B}, y::B) where B
  z = copy(x)
  if haskey(x.fac, y)
    z.fac[y] = z.fac[y] + 1
  else
    z.fac[y] = 1
  end
  return z
end

function *(y::B, x::FacElem{B}) where B
  z = copy(x)
  if haskey(x.fac, y)
    z.fac[y] = z.fac[y] + 1
  else
    z.fac[y] = 1
  end
  return z
end

function div(x::FacElem{B}, y::FacElem{B}) where B
  z = copy(x)
  for a in base(y)
    if haskey(x.fac, a)
      z.fac[a] = z.fac[a] - y.fac[a]
    else
      z.fac[a] = -y.fac[a]
    end
  end
  return z
end

################################################################################
#
#  Apply transformation
#
################################################################################

# return (x1,...,xr)*y
function _transform(x::Array{FacElem{T, S}, 1}, y::fmpz_mat) where {T, S}
  length(x) != nrows(y) &&
              error("Length of array must be number of rows of matrix")

  z = Array{FacElem{T, S}}(undef, ncols(y))

  t = parent(x[1])()

  for i in 1:ncols(y)
    z[i] = x[1]^y[1,i]
    for j in 2:nrows(y)
      if y[j, i] == 0
        continue
      end
      pow!(t, x[j], y[j, i])
      mul!(z[i], z[i], t)
    end
  end
  return z
end

function transform(x::Array{FacElem{S, T}, 1}, y::fmpz_mat) where {S, T}
  return _transform(x, y)
end

function transform_left!(x::Array{FacElem{S, T}, 1}, y::TrafoSwap{fmpz}) where {S, T}
  x[y.i], x[y.j] = x[y.j], x[y.i]
  nothing
end

function transform_left!(x::Array{FacElem{S, T}, 1}, y::TrafoAddScaled{fmpz}) where {S, T}
  x[y.j] = x[y.j] * x[y.i]^y.s
  nothing
end

function transform_left!(x::Array{FacElem{S, T}, 1}, y::TrafoPartialDense{R}) where {S, T, R}
  z = view(deepcopy(x), y.rows)
  xx = view(x, y.rows)
  for i in 1:nrows(y.U)
    xx[i] = z[1]^Int(y.U[i, 1])
    for j in 2:ncols(y.U)
      xx[i] *= z[j]^Int(y.U[i, j])
    end
  end
end

function transform_left!(x::Array{T, 1}, t::TrafoDeleteZero{S}) where {S, T}
  # move ith position to the back
  for j in t.i:length(x)-1
    r = x[j]
    x[j] = x[j + 1]
    x[j + 1] = r
  end
end

function transform_left!(x::Array{FacElem{S, T}, 1}, y::TrafoParaAddScaled{R}) where {S, T, R}
  # [ Ri; Rj] <- [a b; c d]* [ Ri; Rj ]
  ri = deepcopy(x[y.i])
  rj = deepcopy(x[y.j])
  x[y.i] = ri^y.a * rj^y.b
  x[y.j] = ri^y.c * rj^y.d
end

################################################################################
#
#  Evaluate
#
################################################################################

function evaluate(x::FacElem{NfOrdIdl, NfOrdIdlSet}; coprime::Bool = false)
  O = order(base_ring(x))
  if !coprime
    x = simplify(x) # the other method won't work due to one()
  end
  if length(x.fac)==0
    return frac_ideal(O, O(1))
  end
  # still doesn't work
  D = collect(x.fac)
  A = (D[1][1]//1)^Int(D[1][2])
  for i in 2:length(D)
    A = A * (D[i][1]//1)^Int(D[i][2])
  end
  return A
end

function _ev(d::Dict{nf_elem, fmpz}, oe::nf_elem)
  z = deepcopy(oe)
  if length(d)==0
    return z
  elseif length(d)==1
    x = first(d)
    return x[1]^x[2]
  end
  b = empty(d)
  for (k, v) in d
    if iszero(v)
      continue
    end
    if v>-10 && v<10
      mul!(z, z, k^Int(v))
    else
      r = isodd(v) ? 1 : 0
      vv = div(v-r, 2)
      if vv != 0
        b[k] = vv
      end
      if r != 0
        mul!(z, z, k)
      end
    end
  end
  if isempty(b)
    return z
  end
  res = _ev(b, oe)
  mul!(res, res, res)
  mul!(z, z, res)
  return z
end

function _ev(d::Dict{fq_nmod, fmpz}, z::fq_nmod)
  Fq = parent(z)
  if length(d) == 0
    return z
  elseif length(d)==1
    x = first(d)
    return x[1]^Int(x[2])
  end
  b = empty(d)
  for (k, v) in d
    if iszero(v)
      continue
    end
    if abs(v) < 10
      if v >0
        kv = Fq()
        ccall((:fq_nmod_pow, :libflint), Nothing, (Ref{fq_nmod}, Ref{fq_nmod}, Ref{fmpz}, Ref{FqNmodFiniteField}), kv, k, v, Fq)
        mul!(z, z, kv)
      else
        kv = Fq()
        ccall((:fq_nmod_inv, :libflint), Nothing, (Ref{fq_nmod}, Ref{fq_nmod}, Ref{FqNmodFiniteField}), kv, k, Fq)
        ccall((:fq_nmod_pow, :libflint), Nothing, (Ref{fq_nmod}, Ref{fq_nmod}, Ref{fmpz}, Ref{FqNmodFiniteField}), kv, kv, -v, Fq)
        mul!(z, z, kv)
      end
    else
      r = isodd(v) ? 1 : 0
      vv = div(v-r, 2)
      if vv != 0
        b[k] = vv
      end
      if r != 0
        mul!(z, z, k)
      end
    end
  end
  if isempty(b)
    return z
  end
  res = _ev(b, one(Fq))
  mul!(res, res, res)
  mul!(z, z, res)
  return z
end

function _ev(d::Dict{T, fmpz}, oe::T) where T
  z = copy(oe)
  if length(d)==0
    return z
  elseif length(d)==1
    x = first(d)
    return x[1]^x[2]
  end
  b = empty(d)
  for (k,v) in d
    if v>-10 && v<10
      z *= k^Int(v)
    else
      r = isodd(v) ? 1 : 0
      vv = div(v-r, 2)
      if vv!=0
        b[k] = vv
      end
      if r!=0
        z*= k
      end
    end
  end
  return _ev(b, oe)^2*z
end

function one(A::NfOrdFracIdlSet)
  return ideal(order(A), 1)//1
end

function copy(A::NfOrdFracIdl)
  return deepcopy(A)
end

function ^(A::NfOrdFracIdl, d::fmpz)
  return A^Int(d)
end

@doc Markdown.doc"""
  evaluate{T}(x::FacElem{T}) -> T

Expands or evaluates the factored element, i.e. actually computes the
value.
Does "square-and-multiply" on the exponent vectors.
"""
function evaluate(x::FacElem{T}) where T
  return _ev(x.fac, one(base_ring(x)))
end

@doc Markdown.doc"""
    evaluate(x::FacElem{fmpq}) -> fmpq
    evaluate(x::FacElem{fmpz}) -> fmpz

Expands or evaluates the factored element, i.e. actually computes the
the element.
Works by first obtaining a simplified version of the power product
into coprime base elements.
"""
function evaluate(x::FacElem{fmpq})
  return evaluate_naive(simplify(x))
end

function evaluate(x::FacElem{fmpz})
  return evaluate_naive(simplify(x))
end
@doc Markdown.doc"""
    simplify(x::FacElem{fmpq}) -> FacElem{fmpq}
    simplify(x::FacElem{fmpz}) -> FacElem{fmpz}

Simplfies the factored element, i.e. arranges for the base to be coprime.
"""
function simplify(x::FacElem{fmpq})
  y = deepcopy(x)
  simplify!(y)
  return y
end

function simplify(x::FacElem{fmpz})
  y = deepcopy(x)
  simplify!(y)
  return y
end

function simplify!(x::FacElem{fmpq})
  if length(x.fac) <= 1
    return
  end
  cp = coprime_base(vcat([denominator(y) for y = base(x)], [numerator(y) for y=base(x)]))
  ev = Dict{fmpq, fmpz}()
  for p = cp
    if p == 1 || p == -1
      continue
    end
    v = fmpz(0)
    for b = base(x)
      v += valuation(b, abs(p))*x.fac[b]
    end
    if v != 0
      ev[fmpq(abs(p))] = v
    end
  end
  f = b -> b < 0 && isodd(x.fac[b]) ? -1 : 1
  s = prod((f(v) for v in base(x)))
  if s == -1
    ev[fmpq(-1)] = 1
  else
    if length(ev) == 0
      ev[fmpq(1)] = 0
    end
  end
  x.fac = ev
  nothing
end

function simplify!(x::FacElem{fmpz})
  if length(x.fac) == 0
    x.fac[fmpz(1)] = 0
    return
  end
  if length(x.fac) <= 1
    k,v = first(x.fac)
    if isone(k)
      x.fac[k] = 0
    elseif k == -1
      if isodd(v)
        x.fac[k] = 1
      else
        delete!(x.fac, k)
        x.fac[fmpz(1)] = 0
      end
    end
    return
  end
  cp = coprime_base(collect(base(x)))
  ev = Dict{fmpz, fmpz}()
  for p = cp
    if p == 1 || p == -1
      continue
    end
    v = fmpz(0)
    for b = base(x)
      v += valuation(b, abs(p))*x.fac[b]
    end
    if v < 0
      throw(DomainError(v, "Negative valuation in simplify!"))
    end
    if v != 0
      ev[abs(p)] = v
    end
  end
  f = b -> b < 0 && isodd(x.fac[b]) ? -1 : 1
  s = prod(f(v) for v in base(x))
  if s == -1
    ev[-1] = 1
  else
    if length(ev) == 0
      ev[fmpz(1)] = 0
    end
  end
  x.fac = ev
  nothing
end

@doc Markdown.doc"""
    isone(x::FacElem{fmpq}) -> Bool
    isone(x::FacElem{fmpz}) -> Bool
Tests if $x$ represents $1$ without an evaluation.
"""
function isone(x::FacElem{fmpq})
  y = simplify(x)
  return all(iszero, values(y.fac))
end

function isone(x::FacElem{fmpz})
  y = simplify(x)
  return all(iszero, values(y.fac))
end

#TODO: expand the coprime stuff to automatically also get the exponents
@doc Markdown.doc"""
    simplify(x::FacElem{NfOrdIdl, NfOrdIdlSet}) -> FacElem
    simplify(x::FacElem{NfOrdFracIdl, NfOrdFracIdlSet}) -> FacElem

Uses ```coprime_base``` to obtain a simplified version of $x$, ie.
in the simplified version all base ideals will be pariwise coprime
but not neccessarily prime!.
"""
function simplify(x::FacElem{NfOrdIdl, NfOrdIdlSet})
  z = deepcopy(x)
  simplify!(z)
  return z
end

function simplify!(x::FacElem{NfOrdIdl, NfOrdIdlSet})
  if length(x.fac) <= 1
    p = first(keys(x.fac))
    x.fac =  Dict(p => x.fac[p])
    return
  end
  cp = coprime_base(collect(base(x)))
  ev = Dict{NfOrdIdl, fmpz}()
  for p = cp
    if isone(p)
      continue
    end
    assure_2_normal(p)
    v = fmpz(0)
    for b = base(x)
      v += valuation(b, p)*x.fac[b]
    end
    if v != 0
      ev[p] = v
    end
  end
  x.fac = ev
end

function simplify(x::FacElem{NfOrdFracIdl, NfOrdFracIdlSet})
  z = deepcopy(x)
  simplify!(z)
  return z
end

function simplify!(x::FacElem{NfOrdFracIdl, NfOrdFracIdlSet})
  de = factor_coprime(x)
  if length(de)==0
    de = Dict(ideal(order(base_ring(parent(x))), 1) => fmpz(1))
  end
  x.fac = Dict((i//1, k) for (i,k) = de)
end

@doc Markdown.doc"""
    factor_coprime(x::FacElem{NfOrdIdl, NfOrdIdlSet}) -> Dict{NfOrdIdl, Int}
Computed a partial factorisation of $x$, ie. writes $x$ as a product
of pariwise coprime integral ideals.
"""
function factor_coprime(x::FacElem{NfOrdIdl, NfOrdIdlSet})
  x = deepcopy(x)
  simplify!(x)
  return Dict(p=>Int(v) for (p,v) = x.fac)
end

@doc Markdown.doc"""
    factor_coprime(x::FacElem{fmpz}) -> Fac{fmpz}
Computed a partial factorisation of $x$, ie. writes $x$ as a product
of pariwise coprime integers.
"""
function factor_coprime(x::FacElem{fmpz})
  x = deepcopy(x)
  simplify!(x)
  d = Dict(abs(p) => Int(v) for (p,v) = x.fac)
  if haskey(d, fmpz(-1))
    delete!(d, fmpz(-1))
    return Fac(fmpz(-1), d)
  else
    return Fac(fmpz(1), d)
  end
end

@doc Markdown.doc"""
  evaluate_naive{T}(x::FacElem{T}) -> T

Expands or evaluates the factored element, i.e. actually computes the
value. Uses the obvious naive algorithm. Faster for input in finite rings.
"""
function evaluate_naive(x::FacElem{T}) where T
  z = one(base_ring(x))
  d = x.fac
  if length(d)==0
    return z
  end
  for (k,v) in d
    z *= k^v
  end
  return z
end

function ^(a::fmpz, k::fmpz)
  if a == 0
    if k == 0
      return fmpz(1)
    end
    return fmpz(0)
  end

  if a == 1
    return fmpz(1)
  end
  if a == -1
    if isodd(k)
      return fmpz(-1)
    else
      return fmpz(1)
    end
  end
  return a^Int(k)
end

function ^(a::fmpq, k::fmpz)
  if a == 0
    if k == 0
      return fmpq(1)
    end
    return fmpq(0)
  end

  if a == 1
    return fmpq(1)
  end
  if a == -1
    if isodd(k)
      return fmpq(-1)
    else
      return fmpq(1)
    end
  end
  return a^Int(k)
end

#################################################################################
@doc Markdown.doc"""
    max_exp(a::FacElem)
Finds the largest exponent in the factored element $a$
"""
function max_exp(a::FacElem)
  return maximum(values(a.fac))
end

@doc Markdown.doc"""
    min_exp(a::FacElem)
Finds the smallest exponent in the factored element $a$
"""
function min_exp(a::FacElem)
  return minimum(values(a.fac))
end

@doc Markdown.doc"""
    maxabs_exp(a::FacElem)
Finds the largest exponent by absolute value the factored element $a$
"""
function maxabs_exp(a::FacElem)
  return maximum(abs, values(a.fac))
end

function Base.hash(a::FacElem, u::UInt)
  if a.hash == UInt(0)
    h = hash(u, UInt(3127))
    for (k,v) = a.fac
      h = hash(k, hash(v, h))
    end
    a.hash = h
  end
  return a.hash
end

#used (hopefully) only inside the class group
function FacElem(A::Array{nf_elem_or_fac_elem, 1}, v::Array{fmpz, 1})
  if typeof(A[1]) == nf_elem
    B = FacElem(A[1])
  else
    B = A[1]
  end
  B = B^v[1]
  for i=2:length(A)
    if typeof(A[i]) == nf_elem
      C = FacElem(A[i])^v[i]
    else
      C = A[i]^v[i]
    end
    B *= C
  end
  return B
end

#################################################################################
##
##  Auxillary deep copy functions
##
#################################################################################
#
#function _deepcopy{K, V}(x::Dict{K, V})
#  z = typeof(x)()
#  for a in keys(x)
#    z[a] = deepcopy(x[a])
#  end
#  return z
#end

################################################################################
#
#  Parent object overloading
#
################################################################################

function (F::FacElemMon{T})(a::T) where T
  z = F()
  z.fac[a] = fmpz(1)
  return z
end

#function order(A::FacElemMon{IdealSet})
#  return order(A.base_ring)
#end

function order(A::FacElemMon)
  return order(A.base_ring)
end

(::Type{FacElem{T}})(a::FacElem{T}) where {T} = a

(::Type{FacElem})(a::FacElem{T}) where {T} = a


text:
// passing doubles to printInt().

int main() {
	printInt(1.0);
	return 0 ;
}

max_stars_repo_name: jwilk/jtc

max_issues_repo_name: jwilk/jtc

max_forks_repo_name: jwilk/jtc


text:
export COC,DNC,DND,DES1500,CL1500,SDES1500,SCL1500,SDES2500,SCL2500, stateType, actType, acts, discount_f, hMin, hMax,hs, vMin, vMax, vowns, vints, tauMin, tauMax, numTau, taus, resps, interp, accels, velRanges, allowedTrans

# ADVISORY INDICES
COC=0
DNC=1
DND=2
DES1500 = 3
CL1500 = 4
SDES1500=5
SCL1500=6
SDES2500=7
SCL2500=8

# State Type:
stateType = Tuple{Float64,Float64,Float64,Int,Bool}
actType = Int
acts = [0,1,2,3,4,5,6,7,8]

# Default parameters
discount_f = 1.0

### STATE CUTPOINTS ###
# H: Altitude
hMin = -8000.0
hMax = 8000.0
hs   = vcat(LinRange(-8000,-4000,5),LinRange(-3000,-1250,8),LinRange(-1000,-800,3),LinRange(-700,-150,12),LinRange(-100,100,9),LinRange(150,700,12),LinRange(800,1000,3),LinRange(1250,3000,8),LinRange(4000,8000,5))

# Velocities
vMin = -100.
vMax = 100.
vels = vcat(LinRange(-100,-60,5),LinRange(-50,-35,4),LinRange(-30,30,21),LinRange(35,50,4),LinRange(60,100,5))
vowns = vels #LinRange(vMin,vMax,6)
vints = vels #LinRange(vMin,vMax,6)

# Taus
tauMin = 0.0
tauMax = 40.0
numTau = 41
taus  = LinRange(tauMin,tauMax,numTau)

# Responding
resps = [false,true]

interp = LocalGIFunctionApproximator(RectangleGrid(hs,vowns,vints,acts,resps))  # Create the local function approximator using the grid

### Dictionaires to define transitions ###
# Tuple of probabilities and corresponding acceleration in ft/s^2
accels = Dict(COC=>([0.5,0.25,0.25],[0.0,3.0,-3.0]),
              DNC=>([0.5,0.25,0.25],[-8.33,-9.33,-7.33]),
              DND=>([0.5,0.25,0.25],[8.33,9.33,7.33]),
              DES1500=>([0.5,0.25,0.25],[-8.33,-9.33,-7.33]),
              CL1500=>([0.5,0.25,0.25],[8.33,9.33,7.33]),
              SDES1500=>([0.5,0.25,0.25],[-10.7,-11.7,-9.7]),
              SCL1500=>([0.5,0.25,0.25],[10.7,11.7,9.7]),
              SDES2500=>([0.5,0.25,0.25],[-10.7,-11.7,-9.7]),
              SCL2500=>([0.5,0.25,0.25],[10.7,11.7,9.7]))

# Velocity range where aircraft is NON-compliant with advisory
velRanges = Dict(COC=>(-100.0,100.0),
                DNC=>(0.0,100.0),
                DND=>(-100.0,0.0),
                DES1500=>(-25.0,100.0),
                CL1500=>(-100.0,25.0),
                SDES1500=>(-25.0,100.0),
                SCL1500=>(-100.0,25.0),
                SDES2500=>(-41.67,100.0),
                SCL2500=>(-100.0,41.67))

# Allowed transitions between advisories
allowedTrans = Dict(COC=>[1,1,1,1,1,0,0,0,0],
                   DNC=>[1,1,1,1,1,0,0,0,0],
                   DND=>[1,1,1,1,1,0,0,0,0],
                   DES1500=>[1,1,1,1,1,1,1,0,0],
                   CL1500=>[1,1,1,1,1,1,1,0,0],
                   SDES1500=>[1,1,1,1,1,1,1,1,1],
                   SCL1500=>[1,1,1,1,1,1,1,1,1],
                   SDES2500=>[1,1,1,1,1,1,1,1,1],
                   SCL2500=>[1,1,1,1,1,1,1,1,1])

max_stars_repo_name: alisamiii/VerticalCAS

max_issues_repo_name: alisamiii/VerticalCAS

max_forks_repo_name: bethelhall/CoCoCAS


text:
###############################
# CUDA Array and memory types #
###############################

# Abstract CUDA array types
abstract AbstractCudaArray{T,N}
typealias AbstractCudaVector{T} AbstractCudaArray{T,1}
typealias AbstractCudaMatrix{T} AbstractCudaArray{T,2}

typealias HostOrDevArray{T,N} Union{AbstractArray{T,N}, AbstractCudaArray{T,N}}

# copy method for AbstractCudaArray
copy(a::AbstractCudaArray; stream=null_stream) = copy!(similar(a),a;stream=stream)

# Set the following to true to store a backtrace
# at allocation time, to identify where each array comes from
const debugMemory = false

# Contiguous arrays on the GPU
if !debugMemory
    type CudaArray{T,N} <: AbstractCudaArray{T,N}
        ptr::CudaPtr{T}
        dims::NTuple{N,Int}
        dev::Int
    end
else
    type CudaArray{T,N} <: AbstractCudaArray{T,N}
        ptr::CudaPtr{T}
        dims::NTuple{N,Int}
        dev::Int
        bt

        function CudaArray(ptr::CudaPtr{T}, dims::NTuple{N,Int}, dev::Integer)
            new(ptr, dims, dev, backtrace())
        end
    end
end

# Vector and matrix aliases
typealias CudaVector{T} CudaArray{T,1}
typealias CudaMatrix{T} CudaArray{T,2}
typealias CudaVecOrMat{T} Union{CudaVector{T}, CudaMatrix{T}}

# Layout-optimized 1-, 2-, and 3-dimensional arrays
if !debugMemory
    type CudaPitchedArray{T,N} <: AbstractCudaArray{T,N}
        ptr::rt.cudaPitchedPtr
        dims::NTuple{N,Int}
        dev::Int
    end
else
    type CudaPitchedArray{T,N} <: AbstractCudaArray{T,N}
        ptr::rt.cudaPitchedPtr
        dims::NTuple{N,Int}
        dev::Int
        bt

        function CudaPitchedArray(ptr::rt.cudaPitchedPtr, dims::NTuple{N,Int}, dev::Integer)
            new(ptr, dims, dev, backtrace())
        end
    end
end

# Pinned memory on the host
type HostArray{T,N} <: AbstractArray{T,N}
    ptr::Ptr{Void}
    data::Array{T,N}
end

typealias CdArray{T,N} Union{DenseArray{T,N},HostArray{T,N},AbstractCudaArray{T,N}}
typealias ContiguousArray{T,N} Union{Array{T,N},HostArray{T,N},CudaArray{T,N}}

###################
# Implementations #
###################

length(g::AbstractCudaArray) = prod(g.dims)
size(g::AbstractCudaArray) = g.dims
function size{T,N}(g::AbstractCudaArray{T,N}, d::Integer)
    d >= 1 ? (d <= N ? g.dims[d] : 1) : error("Invalid index of dimension.")
end
ndims{T,N}(g::AbstractCudaArray{T,N}) = N
eltype{T,N}(g::AbstractCudaArray{T,N}) = T
device(A::AbstractCudaArray) = A.dev
device(A::AbstractArray) = -1  # for host

pointer(g::AbstractCudaArray) = g.ptr

to_host{T}(g::AbstractCudaArray{T}) = copy!(Array(T, size(g)), g)

summary(g::AbstractCudaArray) = string(g)

######################################
# CudaArray: contiguous array on GPU #
######################################

CudaArray(T::Type, dims::Integer...) = CudaArray(T, dims)

function CudaArray(T::Type, dims::Dims)
    n = prod(dims)
    p = malloc(T, n)
    CudaArray{T,length(dims)}(p, dims, device())
end

CudaArray{T,N}(a::Array{T,N}; stream=null_stream) = copy!(CudaArray(T, size(a)), a; stream=stream)
CudaArray{T,N}(a::HostArray{T,N}; stream=null_stream) = copy!(CudaArray(T, size(a)), a; stream=stream)
CudaArray{T,N}(a::AbstractArray{T,N}) = CudaArray(convert(Array{T,N}, a))

stride(g::CudaArray, dim::Integer) = prod(size(g)[1:dim-1])

similar(g::CudaArray, T, dims::Int...) = similar(g,T,dims)
similar(g::CudaArray, T, dims::Dims) = CudaArray(T, dims)
similar(g::CudaArray) = CudaArray(eltype(g), size(g))
similar(g::CudaArray, T) = CudaArray(T, size(g))
similar(g::CudaArray, dims::Dims) = CudaArray(eltype(g), dims)

unsafe_convert{T}(::Type{Ptr{T}}, g::CudaArray) = unsafe_convert(Ptr{T}, pointer(g))

# convert{T,N}(::Type{CudaArray{T,N}}, g::CudaArray{T,N}) = g
# convert{R,S,N}(::Type{CudaArray{R,N}}, g::CudaArray{S,N}) = CudaArray{R,N}(g.ptr, g.dims, g.dev)
reinterpret{T}(::Type{T}, g::CudaArray{T}) = g
function reinterpret{R,S}(::Type{R}, g::CudaArray{S})
    if sizeof(R) != sizeof(S)
        error("result shape not specified")
    end
    CudaArray{R,ndims(g)}(g.ptr, g.dims, g.dev)
end
function reinterpret{R,S,N}(::Type{R}, g::CudaArray{S}, dims::NTuple{N,Int})
    lenR = div(sizeof(S)*length(g),sizeof(R))
    if prod(dims) != lenR
        throw(DimensionMismatch("New dimensions $dims must be consistent with array of length $lenR"))
    end
    CudaArray{R,N}(g.ptr, dims, g.dev)
end

free(g::CudaArray) = free(pointer(g))

if debugMemory
    function show(io::IO, g::CudaArray)
        println(io, typeof(g))
        println(io, "  ptr: ", g.ptr)
        println(io, "  dims: ", g.dims)
        println(io, "  dev: ", g.dev)
        print(io, "Allocated from:")
        Base.show_backtrace(io, g.bt)
    end
end

# To limit the likelihood of ambiguity warnings with other packages,
# these should be the only two-argument definitions of copy!
copy!(dst::AbstractArray, src::HostOrDevArray; stream=null_stream) = _copy!(dst, src, stream)
copy!(dst::HostOrDevArray, src::HostOrDevArray; stream=null_stream) = _copy!(dst, src, stream)

function _copy!{T}(dst::ContiguousArray{T}, src::ContiguousArray{T}, stream)
    if length(dst) != length(src)
        throw(ArgumentError("Inconsistent array length."))
    end
    nbytes = length(src) * sizeof(T)
    rt.cudaMemcpyAsync(dst, src, nbytes, cudamemcpykind(dst, src), stream)
    return dst
end
_copy!{T}(dst::ContiguousArray{T}, src::ContiguousArray, stream) = _copy!(dst, to_eltype(T, src), stream)
_copy!{T}(dst::AbstractCudaArray{T}, src, stream) = _copy!(dst, copy!(Array(T, size(src)), src), stream)

function fill!{T}(X::CudaArray{T}, val; stream=null_stream)
    valT = convert(T, val)
    func = ptxdict[device()].fns["fill_contiguous", T]
    nsm = attribute(device(), rt.cudaDevAttrMultiProcessorCount)
    mul = min(32, ceil(Int, length(X)/(256*nsm)))
    cudacall(func, mul*nsm, 256, (Ptr{T}, Csize_t, T), X, length(X), valT;
             stream=convert(CuStream, stream))
    X
end

#############################################################################
# CudaPitchedArray: layout-optimized 1-, 2- and 3-dimensional arrays on GPU #
#############################################################################

CudaPitchedArray(T::Type, dims::Integer...) = CudaPitchedArray(T, dims)

function CudaPitchedArray(T::Type, dims::Dims)
    nd = length(dims)
    1 <= nd <= 3 || error("Supports only dimensions 1, 2, or 3")
    p = Array(rt.cudaPitchedPtr, 1)
    ext = CudaExtent(T, dims)
    rt.cudaMalloc3D(p, ext)
    pp = p[1]
    pp = rt.cudaPitchedPtr(pp, dims[1]) # correct the xsize
    cptr = pp.ptr
    dev = device()
    cuda_ptrs[WeakRef(cptr)] = dev
    g = CudaPitchedArray{T,length(dims)}(pp, dims, dev)
    finalizer(g, free)
    g
end

CudaPitchedArray{T,N}(a::Array{T,N}) = copy!(CudaPitchedArray(T, size(a)), a)
CudaPitchedArray{T,N}(a::SharedArray{T,N}) = CudaPitchedArray(sdata(a))
CudaPitchedArray{T,N}(a::AbstractArray{T,N}) = CudaPitchedArray(convert(Array{T,N}, a))

stride(g::CudaPitchedArray, dim::Integer) = dim == 1 ? 1 :
                                            dim == 2 ? pitchel(g) :
                                            dim == 3 ? pitchel(g)*size(g,2) : pitchel(g)*size(g,2)*size(g,3)

similar(g::CudaPitchedArray, T, dims::Dims) = CudaPitchedArray(T, dims)
similar(g::CudaPitchedArray) = CudaPitchedArray(eltype(g), size(g))
similar(g::CudaPitchedArray, T) = CudaPitchedArray(T, size(g))
similar(g::CudaPitchedArray, dims::Dims) = CudaPitchedArray(eltype(g), dims)

unsafe_convert(::Type{rt.cudaPitchedPtr}, g::CudaPitchedArray) = pointer(g)
unsafe_convert{T}(::Type{Ptr{T}}, g::CudaPitchedArray{T}) = unsafe_convert(Ptr{T}, rawpointer(g))

# convert{T,N}(::Type{CudaPitchedArray{T,N}}, g::CudaPitchedArray{T,N}) = g
# convert{R,S,N}(::Type{CudaPitchedArray{R,N}}, g::CudaPitchedArray{S,N}) = CudaPitchedArray{R,N}(g.ptr, g.dims, g.dev)

reinterpret{T}(::Type{T}, g::CudaPitchedArray{T}) = g
function reinterpret{R,S}(::Type{R}, g::CudaPitchedArray{S})
    if sizeof(R) != sizeof(S)
        error("result shape not specified")
    end
    CudaPitchedArray{R,ndims(g)}(g.ptr, g.dims, g.dev)
end
function reinterpret{R,S}(::Type{R}, g::CudaPitchedArray{S}, dims::Dims)
    if length(g) == 0 && prod(dims) == 0
        return CudaPitchedArray{R,length(dims)}(g.ptr, dims, g.dev)
    end
    # Because of the pitch, we don't need the first dimension to agree exactly
    if dims[1]*sizeof(R) > size(g,1)*sizeof(S)
        throw(DimensionMismatch("First dimension $(dims[1]) is too large for $(size(g,1)*sizeof(S)) bytes"))
    end
    # Require agreement for dimensions 2 and 3
    nd = length(dims)
    lg, ld = 1, 1
    for d = 2:nd
        lg *= size(g,d)
        ld *= dims[d]
    end
    for d = nd+1:ndims(g)
        lg *= size(g,d)
    end
    if ld != lg
        throw(DimensionMismatch("New dimensions $dims must be consistent---in all dimensions except the first---with array of dimensions with product $lg"))
    end
    pptr = g.ptr
    newptr = rt.cudaPitchedPtr(pptr.ptr, pptr.pitch, dims[1], nd>1?dims[2]:1)
    CudaPitchedArray{R,length(dims)}(newptr, dims, g.dev)
end

function free(g::CudaPitchedArray)
    p = g.ptr.ptr
    if p != C_NULL && haskey(cuda_ptrs, p)
        delete!(cuda_ptrs, p)
        rt.cudaFree(p)
        g.ptr = rt.cudaPitchedPtr()
    end
end

if debugMemory
    function show(io::IO, g::CudaPitchedArray)
        println(io, typeof(g))
        println(io, "  ptr: ", g.ptr)
        println(io, "  dims: ", g.dims)
        println(io, "  dev: ", g.dev)
        print(io, "Allocated from:")
        Base.show_backtrace(io, g.bt)
    end
end

pitchedptr(g::CudaPitchedArray) = pointer(g)
function pitchedptr(a::ContiguousArray)
    1 <= ndims(a) <= 3 || error("Supports only dimensions 1, 2, or 3")
    rt.cudaPitchedPtr(pointer(a), size(a,1)*sizeof(eltype(a)), size(a,1), size(a,2))
end
function pitchedptr{T}(p::Ptr{T}, sz)
    1 <= length(sz) <= 3 || error("Supports only dimensions 1, 2, or 3")
    rt.cudaPitchedPtr(p, sz[1]*sizeof(T), sz[1], length(sz) > 1 ? sz[2] : 1)
end

pitchbytes(g::CudaArray) = size(g,1)*sizeof(eltype(g))
pitchbytes(g::CudaPitchedArray) = g.ptr.pitch

pitchel(g::AbstractCudaArray) = div(pitchbytes(g),sizeof(eltype(g)))

rawpointer(g::CudaArray) = g.ptr
rawpointer(g::CudaPitchedArray) = g.ptr.ptr

CudaExtent{T}(a::AbstractCudaArray{T}) = CudaExtent(pitchbytes(a), size(a,2), size(a,3))

_copy!{T}(dst::AbstractCudaArray{T}, src::AbstractCudaArray{T}, stream) = cudacopy!(dst, src; stream=stream)
_copy!{T}(dst::CdArray{T}, src::AbstractCudaArray{T}, stream) = cudacopy!(dst, src; stream=stream)
_copy!{T}(dst::AbstractCudaArray{T}, src::CdArray{T}, stream) = cudacopy!(dst, src, stream=stream)
_copy!{T}(dst::AbstractCudaArray{T}, src::CdArray, stream) = cudacopy!(dst, to_eltype(T, src), stream=stream)

function cudacopy!{T}(dst::CdArray{T}, src::CdArray{T}; kwargs...)
    if size(dst) != size(src)
        throw(DimensionMismatch("Size $(size(dst)) of dst is not equal to $(size(src)) of src"))
    end
    copy!(dst, map(d->1:d, size(dst)), src, map(d->1:d, size(src)); kwargs...)
end

copy!{T}(dst::AbstractCudaArray{T}, dstI::Tuple{Vararg{Union{Int,UnitRange{Int}}}}, src::AbstractCudaArray{T}, srcI::Tuple{Vararg{Union{Int,UnitRange{Int}}}}; stream=null_stream) =
    cudacopy!(dst, dstI, src, srcI, stream=stream)
copy!{T}(dst::CdArray{T}, dstI::Tuple{Vararg{Union{Int,UnitRange{Int}}}}, src::AbstractCudaArray{T}, srcI::Tuple{Vararg{Union{Int,UnitRange{Int}}}}; stream=null_stream) =
    cudacopy!(dst, dstI, src, srcI, stream=stream)
copy!{T}(dst::AbstractCudaArray{T}, dstI::Tuple{Vararg{Union{Int,UnitRange{Int}}}}, src::CdArray{T}, srcI::Tuple{Vararg{Union{Int,UnitRange{Int}}}}; stream=null_stream) =
    cudacopy!(dst, dstI, src, srcI, stream=stream)
function cudacopy!{T}(dst::CdArray{T}, dstI::Tuple{Vararg{Union{Int,UnitRange{Int}}}}, src::CdArray{T}, srcI::Tuple{Vararg{Union{Int,UnitRange{Int}}}}; stream=null_stream)
    nd = length(srcI)
    if length(dstI) != nd
        throw(DimensionMismatch("Destination of $(length(dstI)) dimensions differs from dimensionality $nd of src"))
    end
    for i = 1:nd
        length(dstI[i]) == length(srcI[i]) || throw(DimensionMismatch("In dimension $i, destination range of length $(length(dstI[i])) but source range has length $(length(srcI[i]))"))
        first(dstI[i]) >= 1 && last(dstI[i]) <= size(dst, i) || throw(DimensionMismatch("In dimension $i, destination range of $(dstI[i]) is not within array dimension of $(size(dst,i))"))
        first(srcI[i]) >= 1 && last(srcI[i]) <= size(src, i) || throw(DimensionMismatch("In dimension $i, source range of $(srcI[i]) is not within array dimension of $(size(src,i))"))
    end
    checkstrides_pitched(dst)
    checkstrides_pitched(src)
    ext = CudaExtent(eltype(src),map(length, srcI))
    srcpos = CudaPos(eltype(src),map(first, srcI))
    dstpos = CudaPos(eltype(dst),map(first, dstI))
    params = [rt.cudaMemcpy3DParms(C_NULL, srcpos, pitchedptr(src), C_NULL, dstpos, pitchedptr(dst), ext, cudamemcpykind(dst, src))]
    rt.cudaMemcpy3DAsync(params, stream)
    dst
end

function copy!{T}(dst::CdArray{T}, src::AbstractCudaArray{T}, srcI::Tuple{Vararg{Union{Int,UnitRange{Int}}}}; stream=null_stream)
    nd = length(srcI)
    for i = 1:nd
        first(srcI[i]) >= 1 && last(srcI[i]) <= size(src, i) || throw(DimensionMismatch("In dimension $i, source range of $(srcI[i]) is not within array dimension of $(size(src,i))"))
        length(srcI[i]) == size(dst, i) || throw(DimensionMismatch("In dimension $i, the size $(size(dst, i)) of the destination is inconsistent with the size $(length(srcI[i])) being copied"))
    end
    checkstrides_pitched(dst)
    sz = map(length, srcI)
    ext = CudaExtent(eltype(src),sz)
    srcpos = CudaPos(eltype(src),map(first, srcI))
    dstpos = CudaPos(0, 0, 0)
    pdst = pitchedptr(dst)
    params = [rt.cudaMemcpy3DParms(C_NULL, srcpos, pitchedptr(src), C_NULL, dstpos, pdst, ext, rt.cudaMemcpyDeviceToHost)]
    rt.cudaMemcpy3DAsync(params, stream)
    dst
end

function cudaget!{T}(dst::CdArray{T}, src::CdArray{T}, srcI::Tuple{Vararg{Union{Int,UnitRange{Int}}}}, default; kwargs...)
    nd = length(srcI)
    if ndims(dst) != nd
        throw(DimensionMismatch("Destination of $(ndims(dst)) dimensions differs from dimensionality $nd of src"))
    end
    for i = 1:nd
        size(dst,i) == length(srcI[i]) || throw(DimensionMismatch("In dimension $i, destination has dimension $(size(dst,i)) but source range has length $(length(srcI[i]))"))
    end
    srcII = ntuple(i->isa(srcI[i], Int) ? (srcI[i]:srcI[i]) : srcI[i], nd)
    dstIc, srcIc = Base.indcopy(size(src), srcII)
    if srcIc != srcII
        fill!(dst, default; kwargs...)
        copy!(dst, dstIc, src, srcIc; kwargs...)
    else
        copy!(dst, dstIc, src, srcIc; kwargs...)
    end
    dst
end
get!{T}(dst::AbstractCudaArray{T}, src::AbstractCudaArray{T}, srcI::Tuple{Vararg{Union{Int,UnitRange{Int}}}}, default; kwargs...) = cudaget!(dst, src, srcI, default; kwargs...)
get!{T}(dst::AbstractCudaArray{T}, src::CdArray{T}, srcI::Tuple{Vararg{Union{Int,UnitRange{Int}}}}, default; kwargs...) = cudaget!(dst, src, srcI, default; kwargs...)
get!{T}(dst::CdArray{T}, src::AbstractCudaArray{T}, srcI::Tuple{Vararg{Union{Int,UnitRange{Int}}}}, default; kwargs...) = cudaget!(dst, src, srcI, default; kwargs...)

function fill!{T}(X::CudaPitchedArray{T}, val; stream=null_stream)
    valT = convert(T, val)
    func = ptxdict[device()].fns["fill_pitched", T]
    nsm = attribute(device(), rt.cudaDevAttrMultiProcessorCount)
    mul = min(32, ceil(Int, length(X)/(256*nsm)))
    blockspergrid, threadsperblock = ndims(X) == 1 ? (mul*nsm, 256) : (mul*nsm, (16,16))
    cudacall(func, blockspergrid, threadsperblock, (Ptr{T}, Csize_t, Csize_t, Csize_t, Csize_t, T),
             X, size(X,1), size(X,2), size(X,3), pitchel(X), valT;
             stream=convert(CuStream, stream))
    X
end

#####################################
# HostArray (pinned memory on host) #
#####################################

HostArray{T}(::Type{T}, sz...; flags::Integer=rt.cudaHostAllocDefault) = HostArray(T, sz, flags=flags)
function HostArray{T}(::Type{T}, sz::Dims; flags::Integer=rt.cudaHostAllocDefault)
    p = Ptr{Void}[C_NULL]
    rt.cudaHostAlloc(p, prod(sz)*sizeof(T), flags)
    ptr = p[1]
    if VERSION < v"0.5.0-dev+4597"
        data = pointer_to_array(unsafe_convert(Ptr{T}, ptr), sz, false)
    else
        data = unsafe_wrap(Array, unsafe_convert(Ptr{T}, ptr), sz, false)
    end
    ha = HostArray(ptr, data)
    finalizer(ha, free)
    cuda_ptrs[WeakRef(ptr)] = device()
    ha
end

function free(ha::HostArray)
    if ha.ptr != C_NULL && haskey(cuda_ptrs, ha.ptr)
        rt.cudaFreeHost(ha.ptr)
        ha.ptr = C_NULL
        ha.data = Array(eltype(ha), ntuple(d->0, ndims(ha)))
    end
end

size(ha::HostArray) = size(ha.data)
size(ha::HostArray, d) = size(ha.data, d)
ndims{T,N}(ha::HostArray{T,N}) = N
eltype{T,N}(ha::HostArray{T,N}) = T
getindex(ha::HostArray, i::Real) = getindex(ha.data, i)
getindex(ha::HostArray, i::Union{Real,AbstractVector}) = getindex(ha.data, i)
getindex(ha::HostArray, i::Union{Real,AbstractVector}, j::Union{Real,AbstractVector}) = getindex(ha.data, i, j)
getindex(ha::HostArray, i::Union{Real,AbstractVector}...) = getindex(ha.data, i...)
setindex!(ha::HostArray, val, i::Real) = setindex!(ha.data, val, i)
setindex!(ha::HostArray, val, i::Union{Real,AbstractVector}) = setindex!(ha.data, val, i)
setindex!(ha::HostArray, val, i::Union{Real,AbstractVector}, j::Union{Real,AbstractVector}) = setindex!(ha.data, val, i, j)
setindex!(ha::HostArray, val, i...) = setindex!(ha.data, val, i...)
pointer(ha::HostArray)    = pointer(ha.data)
rawpointer(ha::HostArray) = pointer(ha)
pitch(ha::HostArray) = size(ha,1)*sizeof(eltype(ha))
unsafe_convert{T}(::Type{Ptr{Void}}, ha::HostArray{T}) = ha.ptr
fill!(ha::HostArray, val) = fill!(ha.data, val)

function checkstrides_pitched(A)
    stride(A, 1) == 1 || error("A must have a stride of 1 along the first dimension")
    stride(A, 3) == stride(A, 2) * size(A, 2) || error("A must be contiguous for dimension 2")
end

vec{T, N}(v::CudaArray{T,N}) = CudaArray{T, 1}(v.ptr, (length(v),), v.dev)
vec{T}(v::CudaArray{T,1}) = v

to_eltype{T}(::Type{T}, A::AbstractArray{T}) = A
to_eltype{T,S,N}(::Type{T}, A::AbstractArray{S,N}) = convert(Array{T,N}, A)


text:
module day14
using Test

function readinput(filename)
    input = readlines(filename)
    inputs = []
    mask_re = r"^mask = (.{31})$"
    mem_re = r"^mem\[([0-9]*)] = ([0-9]*)$"
    for i in input
        if occursin("mask", i)
            mask = [(idx, parse(Int, c)) for (idx, c) in enumerate(i[8:end]) if c ∈ ['0', '1']]
            push!(inputs, mask)
        end
        m = match(mem_re, i)
        if !(m isa Nothing)
            push!(inputs, [parse(Int, m[1]), parse(Int, m[2])])
        end
    end
    inputs
end

genval(x) = digits(x, base=2, pad=36) |> reverse
bitarray2int(x) = sum(2^(i-1)*σ for (i, σ) in enumerate(reverse(x)))
function maskit(x, mask)
    b = genval(x)
    for (idx, val) in mask
        b[idx] = val
    end
    bitarray2int(b)
end

function solve(filename)
    input = readinput(filename)
    mem = Dict{Int64,Int64}()
    current_mask = []
    for i in input
        if isa(i[1], Int)
            mem[i[1]] = maskit(i[2], current_mask)
        else
            current_mask = i
            continue
        end
    end
    sum(values(mem))
end

@testset "solve" begin
    arg = "input/day14.txt"
    @test solve(arg) == 11179633149677
end

@testset "genval" begin
    @test genval(11)[end-4:end] == [0, 1, 0, 1, 1]
end

@testset "readinput" begin
    arg = "input/day14-test.txt"
    res = [
        [(30, 1), (35, 0)],
        [8, 11],
        [7, 101],
        [8, 0]
    ]
    @test readinput(arg) == res
end

end


text:
using DelayDiffEq, DiffEqProblemLibrary.DDEProblemLibrary
using Test

DDEProblemLibrary.importddeproblems()

const prob = DDEProblemLibrary.prob_dde_constant_1delay_scalar

@testset "Basic iterator" begin
  # compute the solution of the DDE
  sol = solve(prob, MethodOfSteps(BS3());
              dt = 1//2^(4), tstops = [1.5], saveat = 0.01, save_everystep = true)

  # initialize integrator
  integrator = init(prob, MethodOfSteps(BS3());
                    dt = 1//2^(4), tstops = [1.5], saveat = 0.01, save_everystep = true)

  # perform one step
  step!(integrator)
  @test integrator.iter == 1

  # move to next grid point
  integrator.opts.advance_to_tstop = true
  step!(integrator)
  @test integrator.t == 1.5

  # solve the DDE
  solve!(integrator)
  @test integrator.t == 10
  @test integrator.sol(9) ≈ sol(9)

  # move to next grid point
  push!(integrator.opts.tstops, 15)
  step!(integrator)
  @test integrator.t == 15

  # move just one step
  integrator.opts.advance_to_tstop = false
  step!(integrator)
  @test integrator.t > 15
end

@testset "Advanced iterators" begin
  # move to grid point in one step
  integrator1 = init(prob, MethodOfSteps(Tsit5());
                    dt = 1//2^(4), tstops = [0.5], advance_to_tstop = true)
  for (u, t) in tuples(integrator1)
    @test 0.5 ≤ t ≤ 10
  end

  # move to grid point in one step and stop there
  integrator2 = init(prob, MethodOfSteps(Tsit5());
                     dt = 1//2^(4), tstops = [0.5], advance_to_tstop = true,
                     stop_at_next_tstop = true)
  for (u, t) in tuples(integrator2)
    @test t == 0.5
  end
  integrator2([10; 20])

  # step and show intervals
  integrator3 = init(prob, MethodOfSteps(Tsit5()); dt = 1//2^(4), tstops = [0.5])
  for (uprev, tprev, u, t) in intervals(integrator3)
    @show tprev, t
  end
  integrator3([10; 20])

  # iterator for chosen time points
  integrator4 = init(prob, MethodOfSteps(Tsit5()); dt = 1//2^(4))
  ts = 1:10
  us = Float64[]
  for (u, t) in TimeChoiceIterator(integrator4, ts)
    push!(us, u)
  end
  @test us ≈ integrator4.sol(ts)
end

@testset "iter" begin
  integrator = init(prob, MethodOfSteps(RK4()); dt = 1//2^(9), adaptive = false)
  for i in Base.Iterators.take(integrator, 12) end
  @test integrator.iter == 12
  for i in Base.Iterators.take(integrator, 12) end
  @test integrator.iter == 24
end

max_stars_repo_name: JuliaDiffEq/DelayDiffEq

max_issues_repo_name: JuliaDiffEq/DelayDiffEq

max_forks_repo_name: JuliaDiffEq/DelayDiffEq


text:
using Base: Float64
using Dates: string
println("---> START <---")
# Set up environment
using Pkg
Pkg.activate(".")
using JSON
using Glob
using CSV
using Random
using DataFrames

using Statistics
using StatsBase

using Distributions
using DataStructures
using Parameters

using Plots
using Colors
using ColorSchemes
using Dates
using Plotly

include("post_simulation_plots_functions.jl")

# SETUP
outputs_path = "/Users/kingasiuta/Desktop/results/results_5_07"
vars = ["p0"] # what has to be constant
vars_heatmap = ["N_agents", "p0"] # what has to be constant
points_in_time = [[1*24*3600, 2*24*3600, 3*24*3600]]
points_in_time_sing = [1*24*3600, 2*24*3600, 3*24*3600]

# reading possible parameter values

params_to_plot = CSV.read("$outputs_path/sim_params_comb_df.csv", DataFrame)

file_list = glob("simul_res_parset_id_*.json", outputs_path)
    for id in params_to_plot[:,:parset_id]
        if string(outputs_path,"/simul_res_parset_id_",id,".json") ∉ file_list
            filter!(r -> !in(r.parset_id, [id]), params_to_plot)
        end
    end
CSV.write("$outputs_path/sim_params_comb_df.csv", params_to_plot)

values_to_loop = Dict()
for var in vars
    values_to_loop[var] = unique(params_to_plot[:,var])
end

# actual loop
for n_ag in values_to_loop[vars[1]], #  THIS HAS TO BE CHANGED IF VAR NUMBER IS NOT 3
    speed in values_to_loop[vars[2]],
    p0 in values_to_loop[vars[3]]

    full_output, params_sets = read_results(outputs_path,
        control_vars = vars,
        control_val_equal = [n_ag, speed, p0],
        vars_to_check = ["TTC_car_freq"])

    develplot_by_var("total_infected", "TTC_car_freq", full_output,
    with_std = false,
    total_indic = "N_agents",
    div = 60,
    var_unit = "min",
    max_it = nothing,
    x_lab = "iteration",
    y_lab = "% infected",
    title = "% infected vs freq when ag = $n_ag, speed = $speed, p0 = $p0", #COSMETIC
    out_name = "total_infected_S/ag_$n_ag&_sp_$speed&_p0_$p0", #COSMETIC
    with_export = true)

    develplot_by_var("TTC_infected", "TTC_car_freq", full_output,
    with_std = false,
    div = 60,
    var_unit = "min",
    max_it = nothing,
    x_lab = "iteration",
    y_lab = "people infected",
    title = "TTC infected vs freq when ag = $n_ag, speed = $speed, p0 = $p0", #COSMETIC
    out_name = "TTC_infected_S/ag_$n_ag&_sp_$speed&_p0_$p0", #COSMETIC
    with_export = true)

    develplot_by_var("TTC_infected", "TTC_car_freq", full_output,
    with_std = false,
    div = 60,
    var_unit = "min",
    max_it = nothing,
    x_lab = "iteration",
    y_lab = "people infected",
    total_indic = "prct_of_agents_used_TTC",
    title = "% TTC infected vs freq when ag = $n_ag, speed = $speed, p0 = $p0", #COSMETIC
    out_name = "perc_used/ag_$n_ag&_sp_$speed&_p0_$p0", #COSMETIC
    with_export = true)

    for (index, value) in enumerate(points_in_time)

        point_in_time_by_var(["total_infected"], "TTC_car_freq", full_output,
            points_in_time = value,
            total_indic = "N_agents",
            div = 60,
            time_div = 3600,
            time_unit = "h",
            x_lab = "Public Transport Frequency",
            y_lab = "% infected",
            title = "Optimal Public Transport Frequency = 9 min", #COSMETIC
            out_name = "do_networks", #COSMETIC
            with_export = true)
        point_in_time_by_var(["total_infected"], "TTC_car_freq", full_output,
            points_in_time = value,
            total_indic = "N_agents",
            div = 60,
            time_div = 3600,
            time_unit = "h",
            x_lab = nothing,
            y_lab = "% infected",
            title = "%infected when ag = $n_ag, speed = $speed, p0 = $p0", #COSMETIC
            out_name = "points_in_time/ag_$n_ag&_sp_$speed&_p0_$p0&_$index", #COSMETIC
            with_export = true)

        point_in_time_by_var(["TTC_infected"], "TTC_car_freq", full_output,
            points_in_time = value,
            total_indic = "N_agents",
            div = 60,
            time_div = 3600,
            time_unit = "h",
            x_lab = nothing,
            y_lab = "% infected",
            title = "%infected when ag = $n_ag, speed = $speed, p0 = $p0", #COSMETIC
            out_name = "TTC_points_in_time/ag_$n_ag&_sp_$speed&_p0_$p0&_$index", #COSMETIC
            with_export = true)

        point_in_time_by_var(["TTC_infected"], "TTC_car_freq", full_output,
            points_in_time = value,
            total_indic = "prct_of_agents_used_TTC",
            div = 60,
            time_div = 3600,
            time_unit = "h",
            x_lab = nothing,
            y_lab = "% infected",
            title = "%infected when ag = $n_ag, speed = $speed, p0 = $p0", #COSMETIC
            out_name = "perc_used/ag_$n_ag&_sp_$speed&_p0_$p0&_$index", #COSMETIC
            with_export = true)
    end

    for (index, value) in enumerate(points_in_time_sing)
        point_in_time_by_var(["total_infected", "TTC_infected", "street_infected"],
        "TTC_car_freq", full_output,
        points_in_time = [value],
        total_indic = "N_agents",
        div = 60,
        time_div = 3600,
        time_unit = "h",
        x_lab = "TTC frequency",
        y_lab = "infected",
        title = "Odsetek zarażonych vs. częstotliwość kursowania TTC, EOD $index", #COSMETIC
        out_name = "multiple_variables/ag_$n_ag&_sp_$speed&_p0_$p0&_EOD$index", #COSMETIC
        with_export = true)

        plot_XY("prct_of_agents_used_TTC", "total_infected", full_output,
            point_in_time = value,
            total_indic = "N_agents",
            x_lab = "TTC users",
            y_lab = "infected",
            title = "Odsetek zarażonych vs. odsetek osób korzystających z TTC, EOD $index", #COSMETIC
            out_name = "XY_plots/ag_$n_ag&_sp_$speed&_p0_$p0&_EOD$index", #COSMETIC
            with_export = true)

        plot_XYY("prct_of_agents_used_TTC", "total_infected", "TTC_car_freq", full_output,
                point_in_time = value,
                divx = 60,
                total_indic = "N_agents",
                x_lab = "TTC frequency",
                y1_lab = "TTC_users",
                y2_lab = "infected",
                title = "Infected population & TTC usage by TTC frequency, EOD $index",
                out_name = "XYY_plots/ag_$n_ag&_sp_$speed&_p0_$p0&_EOD$index",
                with_export = true)
    end

    # full_output, params_sets = read_results(outputs_path,
    #     control_vars = vars_heatmap,
    #     control_val_equal = [n_ag, p0],
    #     vars_to_check = ["TTC_car_freq", "agents_speed_in_sec_per_m"])

    # make_heatmap("prct_of_agents_used_TTC", "TTC_car_freq",
    #     "agents_speed_in_sec_per_m", full_output, params_sets,
    #     div1 = 60,
    #     div2 = 1,
    #     out_name = "heatmaps/ag_$n_ag&_p0_$p0",
    #     with_export = true)

    println("parameters $n_ag agents, $speed speed, $p0 probability - finished" ) #COSMETIC
end

println("plotting finished")

#### plotting scenarios for p0

total_indic = "N_agents"
response_vars= ["total_infected"]
variable = "TTC_car_freq"
points_in_time = [1*24*3600, 2*24*3600, 3*24*3600, 4*24*3600, 5*24*3600]
div = 60

for p0 in values_to_loop["p0"]
    full_output, params_sets = read_results(outputs_path,
        control_vars = vars,
        control_val_equal = [p0],
        vars_to_check = ["TTC_car_freq"])

    point_in_time_by_var(["total_infected"], "TTC_car_freq", full_output,
    points_in_time = points_in_time,
    total_indic = "N_agents",
    div = div,
    time_div = 3600,
    time_unit = "h",
    x_lab = nothing,
    y_lab = "% infected",
    title = "%infected", #COSMETIC
    out_name = "more_iters_$p0", #COSMETIC
    with_export = true)
end

max_stars_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_issues_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_forks_repo_name: NykPol/EpidemicInUrbanNetworkToronto


text:
using LinearAlgebra

@testset "Quick bounds update smp" begin
    Sigma = [Matrix(1.0I, 4, 4) .+ 1e-4]
    s = discreteswitchedsystem(Sigma)

    qlb, qub = quickb(s, 1)

    @test qlb > 1.

    @test qub > 1.

    @test hassmp(s)

    @test getlb(s) > 1.

    @test getub(s) > 1.
end


text:
module XMem

export load, save, xmem, dfmerge!, getflags

using CSV, DataFrames

include("flags.jl")

# Include the X-Mem executable
const XMEM = joinpath(@__DIR__, "..", "deps", "X-Mem", "bin", "xmem-linux-x64_avx")

load(path) = CSV.File(path) |> DataFrame
save(path, df) = CSV.write(path, df)

# To handle cases where columns were inferred differently.
function dfmerge!(A, B)
    columns = names(A)

    # Collect the promoted types for the columns of A and B
    promote_types = [promote_type(eltype(A[col]), eltype(B[col])) for col in columns]

    for (col, typ) in zip(columns, promote_types)
        A[col] = convert(Vector{typ}, A[col])
        B[col] = convert(Vector{typ}, B[col])
    end
    return append!(A, B)
end

# Single child process versions
xmem(flags::Vector{String}; wait = true, stdout = devnull) = run(pipeline(`$XMEM $flags`; stdout = stdout); wait = wait)

function xmem(flagses::Vector{Vector{T}}) where {T}

    first = true
    local df
    for (index, flags) in enumerate(flagses)
        outfile = getoutfile(flags)
        # Print progress
        if mod(index, 5) == 0
            println("On index $index of $(length(flagses))")
        end

        xmem(flags)
        if first
            df = load(outfile)
            first = false
        else
            dfmerge!(df, load(outfile))
        end
    end
    return df
end

# Dual child process
function xmem(latency_flags::Vector{Vector{T}}, load_flags::Vector{Vector{S}}, outfile) where {T,S}
    if length(latency_flags) != length(load_flags)
        throw(error("atency flags and load flags must be the same length!"))
    end
    first = true
    local df
    for (index, (latency, load)) in enumerate(zip(latency_flags, load_flags))
        if mod(index, 5) == 0
            println("On index $index of $(length(latency_flags))")
        end

        if first
            df = dual_xmem(latency, load, outfile)
            first = false
        else
            dfmerge!(df, dual_xmem(latency, load, outfile))
        end
    end
    save(outfile, df)
    return df
end

# Don't give outfile flag
function dual_xmem(latency_flags, load_flags, outfile)
    # Throw error if outfile flag is given in either load flags
    f(x) = startswith(x, flag_outfile(""))
    if any(f, latency_flags)
        throw(error("Do not include `outfile` in latency flags"))
    end
    if any(f, load_flags)
        throw(error("Do not include `outfile` in load flags"))
    end

    latency_outfile = ".a$outfile"
    load_outfile = ".b$outfile"
    _latency_flags = vcat(latency_flags, [flag_outfile(latency_outfile)])
    _load_flags = vcat(load_flags, [flag_outfile(load_outfile)])
    # Run the two workloads
    latency_process = runxmem(_latency_flags; wait = false)
    while process_running(latency_process)
        runxmem(_load_flags; wait = true)
    end

    # Read and combine both files
    dfa = load(latency_outfile)
    dfb = load(load_outfile)
    dfmerge!!(dfa, dfb)

    # Remove the two temporary files
    rm(latency_outfile)
    rm(load_outfile)
    save(outfile, dfa)

    return dfa
end

end #module


text:
# using Sys built-in
# using command line
# julia> is the prompt

# on a 64-bit OS
julia> Sys.WORD_SIZE
64

# version info
julia> versioninfo()
Julia Version 1.0.0
Commit 5d4eaca0c9 (2018-08-08 20:58 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i5-6300U CPU @ 2.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.0 (ORCJIT, skylake)


text:
#=

Maybe do generated functions to reduce dispatch times?

f(x) = x
g(x,i) = f(x[i])
g{i}(x,::Type{Val{i}}) = f(x[i])
@generated function gg(tup::Tuple, num)
                  N = length(tup.parameters)
                  :(@nif $(N+1) i->(i == num) i->(f(tup[i])) i->error("unreachable"))
              end
h(i) = g((1,1.0,"foo"), i)
h2{i}(::Type{Val{i}}) = g((1,1.0,"foo"), Val{i})
h3(i) = gg((1,1.0,"foo"), i)
@benchmark h(1)
mean time:        31.822 ns (0.00% GC)
@benchmark h2(Val{1})
mean time:        1.585 ns (0.00% GC)
@benchmark h3(1)
mean time:        6.423 ns (0.00% GC)

@generated function foo(tup::Tuple, num)
  N = length(tup.parameters)
  :(@nif $(N+1) i->(i == num) i->(tup[i]) i->error("unreachable"))
end

@code_typed foo((1,1.0), 1)

@generated function perform_step!(integrator, cache::CompositeCache, repeat_step=false)
  N = length(cache.parameters)
  :(@nif $(N+1) i->(i == num) i->(tup[i]) i->error("unreachable"))
end

=#

function initialize!(integrator, cache::CompositeCache)
  cache.current = cache.choice_function(integrator)
  if cache.current == 1
    initialize!(integrator, @inbounds(cache.caches[1]))
  elseif cache.current == 2
    initialize!(integrator, @inbounds(cache.caches[2]))
  else
    initialize!(integrator, @inbounds(cache.caches[cache.current]))
  end
  resize!(integrator.k, integrator.kshortsize)
end

function perform_step!(integrator, cache::CompositeCache, repeat_step=false)
  if cache.current == 1
    perform_step!(integrator, @inbounds(cache.caches[1]), repeat_step)
  elseif cache.current == 2
    perform_step!(integrator, @inbounds(cache.caches[2]), repeat_step)
  else
    perform_step!(integrator, @inbounds(cache.caches[cache.current]), repeat_step)
  end
end

choose_algorithm!(integrator,cache::OrdinaryDiffEqCache) = nothing
function choose_algorithm!(integrator,cache::CompositeCache)
  new_current = cache.choice_function(integrator)
  @inbounds if new_current != cache.current
    if new_current == 1
      initialize!(integrator, @inbounds(cache.caches[1]))
    elseif new_current == 2
      initialize!(integrator, @inbounds(cache.caches[2]))
    else
      initialize!(integrator, @inbounds(cache.caches[new_current]))
    end
    if cache.current == 1 && new_current == 2
      reset_alg_dependent_opts!(integrator,integrator.alg.algs[1],integrator.alg.algs[2])
      transfer_cache!(integrator,integrator.cache.caches[1],integrator.cache.caches[2])
    elseif cache.current == 2 && new_current == 1
      reset_alg_dependent_opts!(integrator,integrator.alg.algs[2],integrator.alg.algs[1])
      transfer_cache!(integrator,integrator.cache.caches[2],integrator.cache.caches[1])
    else
      reset_alg_dependent_opts!(integrator,integrator.alg.algs[cache.current],integrator.alg.algs[new_current])
      transfer_cache!(integrator,integrator.cache.caches[cache.current],integrator.cache.caches[new_current])
    end
    cache.current = new_current
  end
end

"""
If no user default, then this will change the default to the defaults
for the second algorithm.
"""
function reset_alg_dependent_opts!(integrator,alg1,alg2)
  integrator.dtchangeable = isdtchangeable(alg2)
  if integrator.opts.adaptive == isadaptive(alg1)
    integrator.opts.adaptive = isadaptive(alg2)
  end
  if integrator.opts.qmin == qmin_default(alg1)
    integrator.opts.qmin = qmin_default(alg2)
  end
  if integrator.opts.qmax == qmax_default(alg1)
    integrator.opts.qmax == qmax_default(alg2)
  end
  reset_alg_dependent_opts!(integrator.opts.controller, alg1, alg2)
end

# Write how to transfer the cache variables from one cache to the other
# Example: send the history variables from one multistep method to another

transfer_cache!(integrator,alg1,alg2) = nothing


text:
# Jonas August, 1/2021

"Denoises image v by minimizing DTV(u) + λ|u-v|^2 using the 'digital total variation' proposed by Chan, Osher, and Shen (COS) in 'The Digital TV Filter and Nonlinear Denoising', Trans. IEEE Image Process, Feb. 2001.  Returns sequence of iterates."
function chan_osher_shen_tv(v::Array{T,2};λ::S=S(8),itermax::Int=100,devmin²::S=S((2^-10)^2)) where {S<:AbstractFloat,T<:Union{S,Gray{S}}}

    # Pixels should be a floating point type to avoid artifacts, e.g., subtraction a-b where a<b, a>0, b>0.
    # devmin² could be even smaller (e.g., eps(T)^2), but that would require subnormal numbers which may be slow.
    # Also, they might be flushed to zero, depending on details.  I'm mostly only assuming 0 < devmin² << 1, and
    # especially 0 < sqrt(devmin²).

    m,n = size(v)
    u = zeros(eltype(v),0:m+1,0:n+1)  # pad by 1
    u[1:m,1:n] .= v  # Initial condition

    # Enforce 0 outward derivative for u BC.
    u[0,:] .= u[1,:];  u[:,0] .= u[:,1];  u[end,:] .= u[end-1,:];  u[:,end] .= u[:,end-1]

    # devinv = 1 / local deviation
    devinv = similar(u)

    # Dummy BCs for devinv (ultimately ignored).
    devinv[0,:] .= 1;  devinv[:,0] .= 1;  devinv[end,:] .= 1;  devinv[:,end] .= 1

    us = zeros(eltype(v),m,n,itermax)  # TODO: grow history instead of preallocate?
    for iter in 1:itermax
        # 1/(local deviation) image is precomputed for speedup, so not
        # exactly the same behavior as in the COS paper.
        # TODO: don't precompute, so that we have exact COS behavior.  Timing?
        @inbounds for i=1:m, j=1:n
            devinv[i,j] = 1/sqrt(devmin² + sum(y->y^2,(u[i+1,j],u[i,j+1],u[i-1,j],u[i,j-1]) .- u[i,j]))
        end

        @inbounds for i=1:m, j=1:n
            we = devinv[i,j+1]+devinv[i,j]
            wn = devinv[i-1,j]+devinv[i,j]
            ww = devinv[i,j-1]+devinv[i,j]
            ws = devinv[i+1,j]+devinv[i,j]

            # Kill weights along boundary
            # TODO: not much speedup if we pull this BC handling outside i,j loop? looks less than 20% time vs. no BC
            wn = (i == 1) ? 0 : wn
            ws = (i == m) ? 0 : ws
            ww = (j == 1) ? 0 : ww
            we = (j == n) ? 0 : we

            u[i,j] = (λ*v[i,j] + we*u[i,j+1] + wn*u[i-1,j] + ww*u[i,j-1] + ws*u[i+1,j])/(λ+we+wn+ww+ws)
        end

        # Enforce 0 outward derivative for u BC.
        u[0,:] .= u[1,:]; u[:,0] .= u[:,1]; u[end,:] .= u[end-1,:]; u[:,end] .= u[:,end-1]

        us[:,:,iter] = u[1:m,1:n]
    end
    us  # sequence of iterations is returned
end

max_stars_repo_name: jonasaugust/total-variation-methods

max_issues_repo_name: jonasaugust/total-variation-methods

max_forks_repo_name: jonasaugust/total-variation-methods


text:
module SpineRuntimes

include("../gen/Lib.jl")

end # module SpineRuntimes


text:
######### Estimators for fBm and related processes #########

# abstract type AbstractEstimator end
# abstract type AbstractRollingEstimator <: AbstractEstimator end
# abstract type AbstractfBmEstimator <: AbstractEstimator end

include("Estimator_Variogram.jl")
include("Estimator_Scalogram.jl")
include("Estimator_MLE.jl")

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
@testset "complete substitution" begin
    @polyvar x y z

    p = x^2 * y + y^2 + x
    @test (@inferred p(x=>1, y=>2)) == 7
    @test (@inferred p(x=>5)(y=>3)) == p(x=>5, y=>3)
    @test @wrappedallocs(p(x=>1, y=>2)) == 0 # Works on Julia v0.6 but fails on Julia v0.7
    @test (@inferred x(x=>5)) == 5
    @test (@inferred x(y=>2)) == x
    @test (@inferred x(x=>y)) == y

end

@testset "partial substitution" begin
    @polyvar x y z

    p = x^2 * y + y^2 + x
    p2 = @inferred p(x=>4)
    @test length(terms(p2)) == 3
    @test terms(p2)[1].coefficient == 1
    @test degree(terms(p2)[1].monomial) == 2
    @test exponents(monomial(terms(p2)[1])) == (2,)
    @test terms(p2)[2].coefficient == 16
    @test degree(terms(p2)[2].monomial) == 1
    @test exponents(monomial(terms(p2)[2])) == (1,)
    @test terms(p2)[3].coefficient == 4
    @test degree(terms(p2)[3].monomial) == 0
end

@testset "more inference" begin
    @polyvar x y z

    @test @inferred(subs(x, x=>1)) == 1
    @test @inferred(subs(x, x=>1, y=>2)) == 1
    @test @inferred(subs(x^2, x=>3.0)) == 9.0
    @test @inferred(subs(x^2 * y, x=>1)) == y
    @test @inferred(subs(x^2 * y, y=>5)) == 5x^2
    @test @inferred(subs(x^2 * y, y=>1.0, x=>2)) == 4
    @test @inferred(subs(1x, x=>2)) == 2
    @test @inferred(subs(x + 1, x=>1)) == 2
    @test @inferred(subs(x + 2y + z, y=>2.0, z=>π)) ≈ π + 4.0 + x
end

@testset "more call overloads" begin
    @polyvar x y z
    t = 5x^2 * z * 3 * y * z
    @test @inferred(t(x=>2, y=>2)) == @inferred(15 * 2^2 * z^2 * 2)
    m = z * z * z
    @test @inferred(m(x=>5)) == z^3
    @test @inferred(m(z=>2)) == 8
end

@testset "monomial substitution" begin
    @polyvar x y z

    @test @inferred(subs(x^2 * y, y=>z^2)) == x^2 * z^2
end

@testset "tuple substitution" begin
    @polyvar x y z

    @test @inferred(subs(x^2 * y * z^3, (x, y)=>(2, 3))) == (2^2 * 3 * z^3)
    @test @inferred(subs(5x, (x,)=>(6,))) == 30
    @test @inferred(subs(x * y * z^2, (y, x) => (z, z))) == z^4
    @test @inferred(subs(x^2 * y, (x, y) => (z^3, 5))) == 5 * z^6

    @testset "call overloads" begin
        @test @inferred(x((x,)=>(3,))) == 3

        m = z^1 * y^2
        @test @inferred(m((y, z) => (2, 3))) == 2^2 * 3

        t = 3 * m
        @test @inferred(t((y, z) => (2, 3))) == 3 * 2^2 * 3

        p = 3 * m + 1
        @test @inferred(p((y, z) => (2, 3))) == 3 * 2^2 * 3 + 1
    end
end

@testset "swapping variables" begin
    @polyvar x y z

    @test subs(x^2 * y, x=>y^2, y=>x) == y^4 * x
    # From example 1 in MultivariatePolynomials README
    p = 2x + 3.0x*y^2 + y
    @test p((x, y)=>(y, x)) == 2y + 3y*x^2 + x
    @test p((x, y)=>[1, 2]) == 16
end

@testset "non-variable substitution" begin
    @polyvar x y z

    @test @inferred(subs(5, x=>2)) == 5
    @test @inferred(subs(11, (x, y)=>(y, x))) == 11
end


text: Dict{String, Any}("infinity" => Dict{String, Any}("value" => "inf", "type" => "float"), "infinity_neg" => Dict{String, Any}("value" => "-inf", "type" => "float"), "nan_neg" => Dict{String, Any}("value" => "nan", "type" => "float"), "nan" => Dict{String, Any}("value" => "nan", "type" => "float"), "nan_plus" => Dict{String, Any}("value" => "nan", "type" => "float"), "infinity_plus" => Dict{String, Any}("value" => "+inf", "type" => "float"))