text:
for p in ("Knet","ArgParse","Images")
    Pkg.installed(p) == nothing && Pkg.add(p)
end
include(Pkg.dir("Knet","data","mnist.jl"))
include(Pkg.dir("Knet","data","imagenet.jl"))

module CGAN
using Knet
using Images
using ArgParse
using JLD2, FileIO

function main(args)
    o = parse_options(args)
    o[:seed] > 0 && Knet.setseed(o[:seed])

    # load models, data, optimizers
    wd, wg, md, mg = load_weights(o[:atype], o[:zdim], o[:loadfile])
    xtrn,ytrn,xtst,ytst = Main.mnist()
    dtrn = minibatch(xtrn, ytrn, o[:batchsize]; shuffle=true, xtype=o[:atype])
    optd = map(wi->eval(parse(o[:optim])), wd)
    optg = map(wi->eval(parse(o[:optim])), wg)
    z = sample_noise(o[:atype],o[:zdim],prod(o[:gridsize]))
    l0 = map(i->reshape(collect(1:9), 1, 9), 1:9)
    l1 = vec(vcat(l0...))

    if o[:outdir] != nothing && !isdir(o[:outdir])
        mkpath(o[:outdir])
        mkpath(joinpath(o[:outdir],"models"))
        mkpath(joinpath(o[:outdir],"generations"))
    end

    # training
    println("training started..."); flush(STDOUT)
    for epoch = 1:o[:epochs]
        dlossval = glossval = 0
        @time for (x,y) in dtrn
            x0 = 2x-1; x0 = reshape(x0, 784, div(length(x0),784))
            noise = sample_noise(o[:atype],o[:zdim],length(y))
            dlossval += train_discriminator!(wd,wg,md,mg,x0,y,noise,optd,o)
            noise = sample_noise(o[:atype],o[:zdim],length(y))
            glossval += train_generator!(wg,wd,mg,md,noise,y,optg,o)
        end
        dlossval /= length(dtrn); glossval /= length(dtrn)
        println((:epoch,epoch,:dloss,dlossval,:gloss,glossval))
        flush(STDOUT)

        # save models and generations
        if o[:outdir] != nothing
            filename = @sprintf("%04d.png",epoch)
            filepath = joinpath(o[:outdir],"generations",filename)
            plot_generations(
                wg, mg, l1; z=z, savefile=filepath,
                scale=o[:gridscale], gridsize=o[:gridsize])

            filename = @sprintf("%04d.jld2",epoch)
            filepath = joinpath(o[:outdir],"models",filename)
            save_weights(filepath,wd,wg,md,mg)
        end
    end

    return wd,wg,md,mg
end

function parse_options(args)
    s = ArgParseSettings()
    s.description = "Conditional GANs on MNIST."

    @add_arg_table s begin
        ("--atype"; default=(gpu()>=0?"KnetArray{Float32}":"Array{Float32}");
         help="array and float type to use")
        ("--batchsize"; arg_type=Int; default=100; help="batch size")
        ("--zdim"; arg_type=Int; default=100; help="noise dimension")
        ("--epochs"; arg_type=Int; default=20; help="# of training epochs")
        ("--seed"; arg_type=Int; default=-1; help="random seed")
        ("--gridsize"; arg_type=Int; nargs=2; default=[9,9])
        ("--gridscale"; arg_type=Float64; default=2.0)
        ("--optim"; default="Adam(;lr=0.0002, beta1=0.5)")
        ("--loadfile"; default=nothing; help="file to load trained models")
        ("--outdir"; default=nothing; help="output dir for models/generations")
    end

    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true)
    o[:atype] = eval(parse(o[:atype]))
    if o[:outdir] != nothing
        o[:outdir] = abspath(o[:outdir])
    end
    return o
end

function load_weights(atype,zdim,loadfile=nothing)
    if loadfile == nothing
        wd, md = initwd(atype)
        wg, mg = initwg(atype)
    else
        @load loadfile wd wg md mg
        wd = convert_weights(wd, atype)
        wg = convert_weights(wg, atype)
        md = convert_moments(md, atype)
        mg = convert_moments(mg, atype)
    end
    return wd, wg, md, mg
end

function save_weights(savefile,wd,wg,md,mg)
    save(savefile,
         "wd", convert_weights(wd),
         "wg", convert_weights(wg),
         "md", convert_moments(md),
         "mg", convert_moments(mg))
end

function convert_weights(w, atype=Array{Float32})
    w0 = map(wi->convert(atype, wi), w)
    w1 = convert(Array{Any}, w0)
end

function convert_moments(moments,atype=Array{Float32})
    clone = map(mi->bnmoments(), moments)
    for k = 1:length(clone)
        if moments[k].mean != nothing
            clone[k].mean = convert(atype, moments[k].mean)
        end

        if moments[k].var != nothing
            clone[k].var = convert(atype, moments[k].var)
        end
    end
    return convert(Array{Any,1}, clone)
end

function leaky_relu(x, alpha=0.2)
    pos = max(0,x)
    neg = min(0,x) * alpha
    return pos + neg
end

function sample_noise(atype,zdim,nsamples,mu=0.5,sigma=0.5)
    noise = convert(atype, randn(zdim,nsamples))
    normalized = (noise-mu)/sigma
end

function initwd(atype, winit=0.01)
    w = Any[]
    m = Any[]

    # fc layers for inputs (x,y)
    push!(w, winit*randn(1024,784))
    push!(w, zeros(1024,1))
    push!(w, winit*randn(1024,10))
    push!(w, zeros(1024,1))

    # fc layers for concatenated
    push!(w, winit*randn(512,2048))
    push!(w, bnparams(512))
    push!(m, bnmoments())
    push!(w, winit*randn(256,512))
    push!(w, bnparams(256))
    push!(m, bnmoments())

    push!(w, winit*randn(2,256))
    push!(w, zeros(2,1))
    return convert_weights(w,atype), m
end

function dnet(w,x,y,m; training=true, pdrop=0.5, alpha=0.2)
    x1 = leaky_relu.(w[1] * x .+ w[2], alpha)
    x2 = leaky_relu.(w[3][:,y] .+ w[4])
    x3 = vcat(x1,x2)
    x4 = leaky_relu.(batchnorm(w[5]*x3, m[1], w[6]), alpha)
    x5 = leaky_relu.(batchnorm(w[7]*x4, m[2], w[8]), alpha)
    x6 = w[9] * x5 .+ w[10]
end

function dloss(w,m,real_images,real_labels,fake_images,fake_labels,ygold)
    yreal = dnet(w,real_images,ygold,m)
    real_loss = nll(yreal, real_labels)
    yfake = dnet(w,fake_images,ygold,m)
    fake_loss = nll(yfake, fake_labels)
    return real_loss + fake_loss
end

dlossgradient = gradloss(dloss)

function train_discriminator!(wd,wg,md,mg,real_images,ygold,noise,optd,o)
    fake_images = gnet(wg,noise,ygold,mg; training=true)
    nsamples = div(length(real_images),784)
    real_labels = ones(Int64, 1, nsamples)
    fake_labels = 2ones(Int64, 1, nsamples)
    gradients, lossval = dlossgradient(
        wd,md,real_images,real_labels,fake_images,fake_labels,ygold)
    update!(wd, gradients, optd)
    return lossval
end

function initwg(atype, winit=0.01)
    w = Any[]
    m = Any[]

    push!(w, winit*randn(256,100))
    push!(w, bnparams(256))
    push!(m, bnmoments())

    push!(w, winit*randn(256,10))
    push!(w, bnparams(256))
    push!(m, bnmoments())

    push!(w, winit*randn(512,512))
    push!(w, bnparams(512))
    push!(m, bnmoments())

    push!(w, winit*randn(1024,512))
    push!(w, bnparams(1024))
    push!(m, bnmoments())

    push!(w, winit*randn(784,1024))
    push!(w, zeros(784,1))

    return convert_weights(w,atype), m
end

function gnet(w,z,y,m; training=true, pdrop=0.5)
    x0 = relu.(batchnorm(w[1]*z, m[1], w[2]))
    x1 = relu.(batchnorm(w[3][:,y], m[2], w[4]))
    x2 = vcat(x0,x1)
    x3 = relu.(batchnorm(w[5]*x2, m[3], w[6]))
    x4 = relu.(batchnorm(w[7]*x3, m[4], w[8]))
    x5 = tanh.(w[9] * x4 .+ w[10])
end

function gloss(wg,wd,mg,md,noise,ygold,labels)
    fake_images = gnet(wg,noise,labels,mg)
    ypred = dnet(wd,fake_images,labels,md)
    return nll(ypred, ygold)
end

glossgradient = gradloss(gloss)

function train_generator!(wg,wd,mg,md,noise,labels,optg,o)
    ygold = ones(Int64, 1, length(labels))
    gradients, lossval = glossgradient(wg,wd,mg,md,noise,ygold,labels)
    update!(wg,gradients,optg)
    return lossval
end

function plot_generations(
    wg, mg, labels; z=nothing, gridsize=(8,8), scale=1.0, savefile=nothing)
    if z == nothing
        nimg = prod(gridsize)
        zdim = size(wg[1],2)
        atype = wg[1] isa KnetArray ? KnetArray{Float32} : Array{Float32}
        z = sample_noise(atype,zdim,nimg)
    end
    output = Array(0.5*(1+gnet(wg,z,labels,mg; training=false)))
    images = map(i->reshape(output[:,i], 28, 28, 1), 1:size(output,2))
    grid = Main.make_image_grid(images; gridsize=gridsize, scale=scale)
    if savefile == nothing
        display(colorview(Gray, grid))
    else
        save(savefile, grid)
    end
end

splitdir(PROGRAM_FILE)[end] == "conditional-gan.jl" && main(ARGS)

end # module

max_stars_repo_name: ilkerkesen/adversarial

max_issues_repo_name: muhammetfaik/GAN

max_forks_repo_name: muhammetfaik/GAN


text:
module GoogleCodeSearch

using Sockets
using JSON
using HTTP
using GoogleCodeSearch_jll

import Base: show, read, write
export Ctx, index, search, show, indices, clear_indices, paths_indexed, run_http

include("codesearch.jl")
include("server.jl")
include("index.jl")

end # module

ext: jl


text:
export GraphNodeFeaturization, featurize!
export encodable_elements, decode, chunk_vec

"""
    GraphNodeFeaturization(encoded_atom_features::Vector{AtomFeatureDescriptor})
    GraphNodeFeaturization(feature_names, lookup_table, nbins, logspaced, categorical)

A featurization for AtomGraph objects that encodes features associated with each node. Contains a collection of `AtomFeatureDescriptor` objects, and can be initialized by passing those, or by passing parameters for constructing them.

## Notes

## Required Arguments
- `feature_names::Vector{String}`: Names of each feature
- `lookup_table::DataFrame` (optional): source of data for features, if not included in built-in `atom_data_df`

## Keyword Arguments
- `nbins::Union{Vector{Integer},Integer}`: Number of bins to use for one-hot encoding of continuous-valued features. Will be ignored for categorical features.
- `logspaced::Union{Vector{Bool},Bool}`: Whether to logarithmically space the bins
- `categorical::Union{Vector{Bool},Bool}`: Whether each feature is categorical or continous-valued.
"""
struct GraphNodeFeaturization <: AbstractFeaturization
    features::Vector{<:AbstractAtomFeatureDescriptor} # TODO: this might need another name
end

function GraphNodeFeaturization(
    feature_names::Vector{String},
    lookup_table::Union{DataFrame,Nothing} = nothing;
    nbins::Union{Vector{Integer},Integer,Nothing} = nothing,
    logspaced::Union{Vector{Bool},Bool,Nothing} = nothing,
    categorical::Union{Vector{Bool},Bool,Nothing} = nothing,
)
    num_features = length(feature_names)
    local lookup_table_here , logspaced_here , categorical_here , nbins_here
    if isnothing(lookup_table)
        lookup_table_here = atom_data_df
    else
        # need to merge them in case some data is in one place and some the other
        lookup_table_here = outerjoin(atom_data_df, lookup_table, on = :Symbol)
    end

    if isnothing(logspaced)
        logspaced_here = map(fn -> default_log(fn, lookup_table_here), feature_names)
    else
        logspaced_here = get_param_vec(logspaced, num_features)
    end

    if isnothing(categorical)
        categorical_here =
            map(fn -> default_categorical(fn, lookup_table_here), feature_names)
    else
        categorical_here = get_param_vec(categorical, num_features)
    end

    if isnothing(nbins)
        nbins_here = [default_nbins for i = 1:num_features]
    else
        nbins_here = get_param_vec(nbins, num_features, pad_val = default_nbins)
    end

    afs = map(zip(feature_names, nbins_here, logspaced_here, categorical_here)) do args
        ElementFeatureDescriptor(
            args[1],
            lookup_table_here,
            nbins = args[2],
            logspaced = args[3],
            categorical = args[4],
        )
    end
    GraphNodeFeaturization(afs)
end

# TODO: function to compute total vector length?

encodable_elements(fzn::GraphNodeFeaturization) =
    intersect([encodable_elements(f) for f in fzn.features]...)

"""
    chunk_vec(vec, nbins)

Helper function that divides up an already-constructed feature vector into "chunks" (one for each feature) of lengths specified by the vector nbins.

Sum of nbins should be equal to the length of vec.

# Examples
```jldoctest
julia> chunk_vec([1,0,0,1,0], [3,2])
2-element Array{Array{Bool,1},1}:
 [1, 0, 0]
 [1, 0]
 ```
"""
function chunk_vec(vec::Vector{<:Real}, nbins::Vector{<:Integer})
    chunks = fill(Bool[], size(nbins, 1))
    @assert length(vec) == sum(nbins) "Total number of bins doesn't match length of feature vector."
    for i = 1:size(nbins, 1)
        if i == 1
            start_ind = 1
        else
            start_ind = sum(nbins[1:i-1]) + 1
        end
        chunks[i] = vec[start_ind:start_ind+nbins[i]-1]
    end
    return chunks
end

function featurize!(ag::AtomGraph, fzn::GraphNodeFeaturization)
    encoded = reduce(vcat, map((x) -> x(ag), fzn.features))
    ag.encoded_features = encoded
    ag.featurization = fzn
end

function decode(fzn::GraphNodeFeaturization, encoded::Matrix{<:Real})
    num_atoms = size(encoded, 2)
    nbins = [f.length for f in fzn.features]
    local decoded = Dict{Integer,Dict{String,Any}}()
    for i = 1:num_atoms
        #println("atom $i")
        chunks = chunk_vec(encoded[:, i], nbins)
        decoded[i] = Dict{String,Any}()
        for (chunk, f) in zip(chunks, fzn.features)
            #println("    $(f.name): $(decode(f, chunk))")
            decoded[i][f.name] = decode(f, chunk)
        end
    end
    return decoded
end

function decode(ag::AtomGraph)
    @assert !(any(isnothing.([ag.featurization, ag.encoded_features])))
    decoded = decode(ag.featurization, ag.encoded_features)
    for (k,v) in decoded
        v["Symbol"] = ag.elements[k]
    end
    return decoded
end


text:
using RandomMatrixDistributions
using Test

@testset "MarchenkoPastur" begin

    @test_throws ErrorException MarchenkoPastur(5)

    dist = MarchenkoPastur(1)

    @test minimum(dist) == 0
    @test maximum(dist) == 4

    @test pdf(dist, -1) == pdf(dist, 5) == 0
    @test pdf(dist, 2) ≈ 1/(2*pi)
end

@testset "Wachter" begin

    @test_throws ErrorException Wachter(10, 10)

    dist = Wachter(1, 0.5)

    @test minimum(dist) == 0
    @test maximum(dist) == 16

    @test pdf(dist, -1) == pdf(dist, 17) == 0
    @test pdf(Wachter(0.5, 0), 1) ≈ pdf(MarchenkoPastur(0.5), 1)
end


text:
module DLPackInterface

import Requires: @require
import ArrayInterface

export DLDeviceType, DLDevice, DLDataTypeCode, DLDataType, DLTensor
export DL_CPU, DL_CUDA, DL_CUDAHost, DL_OpenCL, DL_Vulkan, DL_Metal, DL_VPI, DL_ROCM, DL_ROCMHost, DL_ExtDev, DL_CUDAManaged, DL_OneAPI, DL_WebGPU, DL_Hexagon
export DL_Int, DL_UInt, DL_Float, DL_OpaqueHandle, DL_Bfloat, DL_Complex
export isdleltype, dltensor, dldata, dldataoffset, dldevice, dldatatype, dlndim, dlshape, dlstride, dlstrides, dlsizeof, dlinfo, dllength

### Types
#
# Note: These are similar to but NOT the same as the C structs from dlpack.h.

@enum DLDeviceType begin
    DL_CPU = 1
    DL_CUDA = 2
    DL_CUDAHost = 3
    DL_OpenCL = 4
    DL_Vulkan = 7
    DL_Metal = 8
    DL_VPI = 9
    DL_ROCM = 10
    DL_ROCMHost = 11
    DL_ExtDev = 12
    DL_CUDAManaged = 13
    DL_OneAPI = 14
    DL_WebGPU = 15
    DL_Hexagon = 16
end

"""
    DLDevice(type::DLDeviceType, id=0)

A device of the given type and id.
"""
struct DLDevice
    type::DLDeviceType
    id::Int
end
DLDevice(t) = DLDevice(t, 0)

@enum DLDataTypeCode begin
    DL_Int = 0
    DL_UInt = 1
    DL_Float = 2
    DL_OpaqueHandle = 3
    DL_Bfloat = 4
    DL_Complex = 5
end

"""
    DLDataType(code::DLDataTypeCode, bits, lanes=1)
    DLDataType(T::Type, lanes=1)

A datatype with the specified number of bits.
"""
struct DLDataType
    code::DLDataTypeCode
    bits::Int
    lanes::Int
end
DLDataType(code, bits) = DLDataType(code, bits, 1)

const TYPE_TO_DATATYPE = Dict(
    Int8 => DLDataType(DL_Int, 8),
    Int16 => DLDataType(DL_Int, 16),
    Int32 => DLDataType(DL_Int, 32),
    Int64 => DLDataType(DL_Int, 64),
    UInt8 => DLDataType(DL_UInt, 8),
    UInt16 => DLDataType(DL_UInt, 16),
    UInt32 => DLDataType(DL_UInt, 32),
    UInt64 => DLDataType(DL_UInt, 64),
    Float16 => DLDataType(DL_Float, 16),
    Float32 => DLDataType(DL_Float, 32),
    Float64 => DLDataType(DL_Float, 64),
    ComplexF16 => DLDataType(DL_Complex, 32),
    ComplexF32 => DLDataType(DL_Complex, 64),
    ComplexF64 => DLDataType(DL_Complex, 128),
)

DLDataType(::Type, lanes=1) = error("not a valid DLPack data type")
for (T, t) in TYPE_TO_DATATYPE
    @eval DLDataType(::Type{$T}, lanes=1) = DLDataType($(t.code), $(t.bits), lanes)
end

### Interface functions

"""
    isdleltype(T::Type)

True if `T` is a valid eltype for a DLTensor.
"""
function isdleltype end
@eval isdleltype(::Type{T}) where {T} = T in $(Tuple(keys(TYPE_TO_DATATYPE)))

"""
    dltensor(x)

An object `t` which can be passed to `dldata`, `dldevice`, etc. to determine information
about the tensor `x`.

As long as `t` is not garbarge collected, the pointer `dldata(t)` remains valid.
"""
dltensor(x) = nothing
dltensor(x::T) where {T<:AbstractArray} = isdleltype(eltype(T)) && ArrayInterface.defines_strides(T) && ArrayInterface.device(T) === ArrayInterface.CPUPointer() ? x : nothing
dltensor(x::T) where {T<:StridedArray} = isdleltype(eltype(T)) ? x : nothing

"""
    dldata(t) :: Ptr{Cvoid}

A pointer to the data in tensor `t`.
"""
dldata(t) = convert(Ptr{Cvoid}, pointer(t))

"""
    dldataoffset(t) :: Int

The data in tensor `t` starts at `dldata(t) + dldataoffset(t)`.
"""
dldataoffset(t) = 0

"""
    dldevice(t) :: DLDevice

The device hosting the data in tensor `t`.
"""
dldevice(t) = DLDevice(DL_CPU)

"""
    dldatatype(t) :: DLDataType

The type of data in tensor `t`.
"""
dldatatype(t) = DLDataType(eltype(t))

"""
    dlndim(t) :: Int

The number of dimensions of tensor `t`.
"""
dlndim(t) = Int(ndims(t))

"""
    dlshape(t, i) :: Int

The size of tensor `t` along the `i`th dimension.
"""
dlshape(t, i) = Int(size(t, i))

"""
    dlstride(t, i) :: Int

The distance in memory (in elements, not bytes) between elements of `t` along the `i`th
dimension.
"""
dlstride(t, i) = Int(stride(t, i))

### Utility functions

"""
    dlshape(t) :: NTuple{N,Int}

The size of tensor `t` in each dimension.
"""
dlshape(t) = ntuple(i->dlshape(t, i)::Int, dlndim(t)::Int)

"""
    dlstrides(t) :: NTuple{N,Int}

The distance in memory (in elements, not bytes) between elements of `t` along each
dimension.
"""
dlstrides(t) = ntuple(i->dlstride(t, i)::Int, dlndim(t)::Int)

"""
    dllength(t) :: Int

The number of elements in the tensor `t`.
"""
dllength(t) = prod(dlshape(t))::Int

"""
    dlsizeof(t::DLDataType) :: Int

The size in bytes of an element of this type.
"""
dlsizeof(t::DLDataType) = cld(Int(t.bits) * Int(t.lanes), 8)::Int

"""
    dlsizeof(t) :: Int

The size in bytes of the tensor `t`, if it were stored contiguously.
"""
dlsizeof(t) = dllength(t) * dlsizeof(dldatatype(t)::DLDataType)

"""
    dlinfo(x) :: Union{Nothing,NamedTuple}

If `x` satisfies the DLPack interface, return a NamedTuple of information about it.
Otherwise, return `nothing`.

It has fields `tensor`, `pointer`, `pointeroffset`, `device`, `datatype`, `ndim`, `shape`,
`strides` containing the outputs of the corresponding `dl*` functions.
"""
function dlinfo(x)
    if (t = dltensor(x)) !== nothing
        (
            tensor = t,
            pointer = dldata(t)::Ptr{Cvoid},
            pointeroffset = dldataoffset(t)::Int,
            device = dldevice(t)::DLDevice,
            datatype = dldatatype(t)::DLDataType,
            ndim = dlndim(t)::Int,
            shape = dlshape(t)::Tuple{Vararg{Int}},
            strides = dlstrides(t)::Tuple{Vararg{Int}},
        )
    end
end

### Integrations

function __init__()
    @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" @eval begin
        dltensor(x::CUDA.CuArray) = x
        dldata(x::CUDA.CuArray) = reinterpret(Ptr{Cvoid}, pointer(x))
        dldevice(x::CUDA.CuArray) = DLDevice(CUDA.device(x))
        DLDevice(d::CUDA.CuDevice) = DLDevice(DL_CUDA, d.handle)
    end
end

end # module


text:
using Reachability: Options, SafeStatesProperty
using MathematicalSystems, LazySets
using DynamicPolynomials, SemialgebraicSets

"""
    vanderpol_MP(; [T], [X0], [variables])

Construct the Van der Pol model using MultivariatePolynomials.

### Input

- `T`         -- (optional, default: `7.0`) the time horizon for the initial
                  value problem
- `X0`        -- (optional, default: `[1.25, 1.55] × [2.35, 2.45]`) set of initial states
- `variables` -- (optional, default: `PolyVar`) the set of polynomial variables that
                  are used in the equations
### Output

The tuple `(𝑃, 𝑂)` where `𝑃` is an initial-value problem and `𝑂` are the options.
"""
function vanderpol_MP(; T=7.0,
                        X0=Hyperrectangle(low=[1.25, 2.35], high=[1.55, 2.45]),
                        variables=@polyvar x₁ x₂)

    𝑂 = Options()
    x₁, x₂ = variables
    𝑂[:variables] = variables

    # instantiate the polynomial system
    f = [x₂, x₂ - x₁ - x₁^2 * x₂]
    𝐹 = PolynomialContinuousSystem(f)

    # instantiate the IVP
    𝑃 = InitialValueProblem(𝐹, X0)

    # time horizon
    𝑂[:T] = T

    # variables to plot
    𝑂[:plot_vars] = [1, 2]

    # safety property
    𝑂[:property] = SafeStatesProperty(HalfSpace([0., 1.], 2.75))

    return (𝑃, 𝑂)
end

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
struct SimpleCurve
    start::Tuple{Number, Number}
    stop::Tuple{Number, Number}
    control::Tuple{Number, Number}

    SimpleCurve(start::Tuple{Number, Number}, stop::Tuple{Number, Number}, control::Tuple{Number, Number}=(start .+ stop) ./ 2) = new(start, stop, control)
end

function getPoint(curve::SimpleCurve, percent::Number)
    x = (1 - percent)^2 * curve.start[1] + 2 * (1 - percent) * percent * curve.control[1] + percent^2 * curve.stop[1]
    y = (1 - percent)^2 * curve.start[2] + 2 * (1 - percent) * percent * curve.control[2] + percent^2 * curve.stop[2]

    return x, y
end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
"""
    PMM_get_problem(fnum; uldim = 5, lldim = 5)

Return upper and lower level objective functions, and the problem configuration in a
dictionary.

### Example

```
julia> using HardTestProblems

julia> F, f, conf = PMM_get_problem(2,uldim=2,lldim=3);

julia> conf
Dict{Symbol, Any} with 15 entries:
  :follower_optimum      => 0.0
  :n_inequality_follower => 0
  :xbest                 => [0.0, 0.0]
  :problem               => "PMM2"
  :n_equality_follower   => 0
  :lldim                 => 3
  :uldim                 => 2
  :n_equality_leader     => 0
  :n_inequality_leader   => 0
  :xmin                  => [-10.0, -10.0]
  :leader_optimum        => 0.0
  :xmax                  => [10.0, 10.0]
  :ymax                  => [10.0, 10.0, 10.0]
  :ymin                  => [-10.0, -10.0, -10.0]
  :ybest                 => [0.0, 1.41421, 1.73205]

```
"""
function PMM_get_problem(fnum; uldim = 5, lldim = 5)
    F = eval(Symbol("PMM"*string(fnum) * "_leader"))
    f = eval(Symbol("PMM"*string(fnum) * "_follower"))

    optimum = eval(Symbol("PMM"*string(fnum) * "_optimum"))

    x, y = optimum(uldim, lldim)

    Fmin,_,_ = F(x,y)
    fmin,_,_ = f(x,y)

    xmin = fill(-10.0, uldim)
    xmax = fill( 10.0, uldim)

    ymin = fill(-10.0, lldim)
    ymax = fill( 10.0, lldim)

    # upperMember, lowerMember
    conf = Dict(:xmin  => xmin,
                :xmax  => xmax,
                :ymin  => ymin,
                :ymax  => ymax,
                :xbest => x,
                :ybest => y,
                :leader_optimum   => Fmin,
                :follower_optimum => fmin,
                :lldim => lldim,
                :uldim => uldim,
                :n_inequality_leader   => 0,
                :n_equality_leader     => 0,
                :n_inequality_follower => 0,
                :n_equality_follower   => 0,
                :problem => "PMM" * string(fnum)
               )

    return F, f, conf
end

export PMM_get_problem


text:
################################################################################
# Roessler SRK for second order weak approx

struct DRI1ConstantCache{T,T2} <: StochasticDiffEqConstantCache
  # hard-coded version
  a021::T
  a031::T
  a032::T

  a121::T
  a131::T
  #a132::T

  #a221::T
  #a231::T
  #a232::T

  b021::T
  b031::T
  #b032::T

  b121::T
  b131::T
  #b132::T

  b221::T
  b222::T
  b223::T
  b231::T
  b232::T
  b233::T

  α1::T
  α2::T
  α3::T

  c02::T2
  c03::T2

  #c11::T2
  c12::T2
  c13::T2

  #c21::T2
  #c22::T2
  #c23::T2

  beta11::T
  beta12::T
  beta13::T

  #beta21::T
  beta22::T
  beta23::T

  beta31::T
  beta32::T
  beta33::T

  #beta41::T
  beta42::T
  beta43::T

  #quantile(Normal(),1/6)
  NORMAL_ONESIX_QUANTILE::T
end

function DRI1ConstantCache(::Type{T}, ::Type{T2}) where {T,T2}

  a021 = convert(T, 1//2)
  a031 = convert(T, -1)
  a032 = convert(T, 2)

  a121 = convert(T, 342//491)
  a131 = convert(T, 342//491)

  b021 = convert(T, (6-sqrt(6))/10)
  b031 = convert(T, (3+2*sqrt(6))/5)

  b121 = convert(T, 3*sqrt(38//491))
  b131 = convert(T, -3*sqrt(38/491))

  b221 = convert(T, -214//513*sqrt(1105//991))
  b222 = convert(T, -491//513*sqrt(221//4955))
  b223 = convert(T, -491//513*sqrt(221//4955))
  b231 = convert(T, 214//513*sqrt(1105//991))
  b232 = convert(T, 491//513*sqrt(221//4955))
  b233 = convert(T, 491//513*sqrt(221//4955))

  α1 = convert(T, 1//6)
  α2 = convert(T, 2//3)
  α3 = convert(T, 1//6)

  c02 = convert(T2, 1//2)
  c03 = convert(T2, 1)

  c12 = convert(T2, 342//491)
  c13 = convert(T2, 342//491)

  beta11 = convert(T, 193//684)
  beta12 = convert(T, 491//1368)
  beta13 = convert(T, 491//1368)

  beta22 = convert(T, 1//6*sqrt(491//38))
  beta23 = convert(T, -1//6*sqrt(491//38))

  beta31 = convert(T, -4955//7072)
  beta32 = convert(T, 4955//14144)
  beta33 = convert(T, 4955//14144)

  beta42 = convert(T, -1//8*sqrt(4955//221))
  beta43 = convert(T, 1//8*sqrt(4955//221))

  NORMAL_ONESIX_QUANTILE = convert(T, -0.9674215661017014)

  DRI1ConstantCache(a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::DRI1,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  DRI1ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function RI1ConstantCache(T::Type, T2::Type)
  a021 = convert(T, 2//3) # convert(T, 2//3)
  a031 = convert(T, -1//3)
  a032 = convert(T, 1)

  a121 = convert(T, 1)
  a131 = convert(T, 1)

  b021 = convert(T, 1)
  b031 = convert(T, 0)

  b121 = convert(T, 1)
  b131 = convert(T, -1)

  b221 = convert(T, 1)
  b222 = convert(T, 0)
  b223 = convert(T, 0)
  b231 = convert(T, -1)
  b232 = convert(T, 0)
  b233 = convert(T, 0)

  α1 = convert(T, 1//4)
  α2 = convert(T, 1//2)
  α3 = convert(T, 1//4)

  c02 = convert(T2, 2//3)
  c03 = convert(T2, 2//3)

  c12 = convert(T2, 1)
  c13 = convert(T2, 1)

  beta11 = convert(T, 1//2)
  beta12 = convert(T, 1//4)
  beta13 = convert(T, 1//4)

  beta22 = convert(T, 1//2)
  beta23 = convert(T, -1//2)

  beta31 = convert(T, -1//2)
  beta32 = convert(T, 1//4)
  beta33 = convert(T, 1//4)

  beta42 = convert(T, 1//2)
  beta43 = convert(T, -1//2)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  DRI1ConstantCache(a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::RI1,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  RI1ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function RI3ConstantCache(T::Type, T2::Type)
  a021 = convert(T, 1)
  a031 = convert(T, 1//4)
  a032 = convert(T, 1//4)

  a121 = convert(T, 1)
  a131 = convert(T, 1)

  b021 = convert(T, (3-2*sqrt(6))/5)
  b031 = convert(T, (6+sqrt(6))/10)

  b121 = convert(T, 1)
  b131 = convert(T, -1)

  b221 = convert(T, 1)
  b222 = convert(T, 0)
  b223 = convert(T, 0)
  b231 = convert(T, -1)
  b232 = convert(T, 0)
  b233 = convert(T, 0)

  α1 = convert(T, 1//6)
  α2 = convert(T, 1//6)
  α3 = convert(T, 2//3)

  c02 = convert(T2, 1)
  c03 = convert(T2, 1//2)

  c12 = convert(T2, 1)
  c13 = convert(T2, 1)

  beta11 = convert(T, 1//2)
  beta12 = convert(T, 1//4)
  beta13 = convert(T, 1//4)

  beta22 = convert(T, 1//2)
  beta23 = convert(T, -1//2)

  beta31 = convert(T, -1//2)
  beta32 = convert(T, 1//4)
  beta33 = convert(T, 1//4)

  beta42 = convert(T, 1//2)
  beta43 = convert(T, -1//2)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  DRI1ConstantCache(a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::RI3,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  RI3ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function RI5ConstantCache(T::Type, T2::Type)
  a021 = convert(T, 1)
  a031 = convert(T, 25//144)
  a032 = convert(T, 35//144)

  a121 = convert(T, 1//4)
  a131 = convert(T, 1//4)

  b021 = convert(T, 1//3)
  b031 = convert(T, -5//6)

  b121 = convert(T, 1//2)
  b131 = convert(T, -1//2)

  b221 = convert(T, 1)
  b222 = convert(T, 0)
  b223 = convert(T, 0)
  b231 = convert(T, -1)
  b232 = convert(T, 0)
  b233 = convert(T, 0)

  α1 = convert(T, 1//10)
  α2 = convert(T, 3//14)
  α3 = convert(T, 24//35)

  c02 = convert(T2, 1)
  c03 = convert(T2, 5//12)

  c12 = convert(T2, 1//4)
  c13 = convert(T2, 1//4)

  beta11 = convert(T, 1)
  beta12 = convert(T, -1)
  beta13 = convert(T, -1)

  beta22 = convert(T, 1)
  beta23 = convert(T, -1)

  beta31 = convert(T, 1//2)
  beta32 = convert(T, -1//4)
  beta33 = convert(T, -1//4)

  beta42 = convert(T, 1//2)
  beta43 = convert(T, -1//2)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  DRI1ConstantCache(a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::RI5,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  RI5ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function RI6ConstantCache(T::Type, T2::Type)
  a021 = convert(T, 1)
  a031 = convert(T, 0)
  a032 = convert(T, 0)

  a121 = convert(T, 1)
  a131 = convert(T, 1)

  b021 = convert(T, 1)
  b031 = convert(T, 0)

  b121 = convert(T, 1)
  b131 = convert(T, -1)

  b221 = convert(T, 1)
  b222 = convert(T, 0)
  b223 = convert(T, 0)
  b231 = convert(T, -1)
  b232 = convert(T, 0)
  b233 = convert(T, 0)

  α1 = convert(T, 1//2)
  α2 = convert(T, 1//2)
  α3 = convert(T, 0)

  c02 = convert(T2, 1)
  c03 = convert(T2, 0)

  c12 = convert(T2, 1)
  c13 = convert(T2, 1)

  beta11 = convert(T, 1//2)
  beta12 = convert(T, 1//4)
  beta13 = convert(T, 1//4)

  beta22 = convert(T, 1//2)
  beta23 = convert(T, -1//2)

  beta31 = convert(T, -1//2)
  beta32 = convert(T, 1//4)
  beta33 = convert(T, 1//4)

  beta42 = convert(T, 1//2)
  beta43 = convert(T, -1//2)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  DRI1ConstantCache(a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::RI6,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  RI6ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function RDI2WMConstantCache(T::Type, T2::Type)
  a021 = convert(T, 1)
  a031 = convert(T, 0)
  a032 = convert(T, 0)

  a121 = convert(T, 2//3)
  a131 = convert(T, 2//3)

  b021 = convert(T, 1)
  b031 = convert(T, 0)

  b121 = convert(T, sqrt(2//3))
  b131 = convert(T, -sqrt(2//3))

  b221 = convert(T, sqrt(2))
  b222 = convert(T, 0)
  b223 = convert(T, 0)
  b231 = convert(T, -sqrt(2))
  b232 = convert(T, 0)
  b233 = convert(T, 0)

  α1 = convert(T, 1//2)
  α2 = convert(T, 1//2)
  α3 = convert(T, 0)

  c02 = convert(T2, 1)
  c03 = convert(T2, 0)

  c12 = convert(T2, 2//3)
  c13 = convert(T2, 2//3)

  beta11 = convert(T, 1//4)
  beta12 = convert(T, 3//8)
  beta13 = convert(T, 3//8)

  beta22 = convert(T, sqrt(6)/4)
  beta23 = convert(T, -sqrt(6)/4)

  beta31 = convert(T, -1//4)
  beta32 = convert(T, 1//8)
  beta33 = convert(T, 1//8)

  beta42 = convert(T, sqrt(2)/4)
  beta43 = convert(T, -sqrt(2)/4)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  DRI1ConstantCache(a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::RDI2WM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  RDI2WMConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function RDI3WMConstantCache(T::Type, T2::Type)
  a021 = convert(T, 1//2)
  a031 = convert(T, 0)
  a032 = convert(T, 3//4)

  a121 = convert(T, 2//3)
  a131 = convert(T, 2//3)

  b021 = convert(T, (9-2*sqrt(15))/14)
  b031 = convert(T, (18+3*sqrt(15))/28)

  b121 = convert(T, sqrt(2//3))
  b131 = convert(T, -sqrt(2//3))

  b221 = convert(T, sqrt(2))
  b222 = convert(T, 0)
  b223 = convert(T, 0)
  b231 = convert(T, -sqrt(2))
  b232 = convert(T, 0)
  b233 = convert(T, 0)

  α1 = convert(T, 2//9)
  α2 = convert(T, 1//3)
  α3 = convert(T, 4//9)

  c02 = convert(T2, 1//2)
  c03 = convert(T2, 3//4)

  c12 = convert(T2, 2//3)
  c13 = convert(T2, 2//3)

  beta11 = convert(T, 1//4)
  beta12 = convert(T, 3//8)
  beta13 = convert(T, 3//8)

  beta22 = convert(T, sqrt(6)/4)
  beta23 = convert(T, -sqrt(6)/4)

  beta31 = convert(T, -1//4)
  beta32 = convert(T, 1//8)
  beta33 = convert(T, 1//8)

  beta42 = convert(T, sqrt(2)/4)
  beta43 = convert(T, -sqrt(2)/4)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  DRI1ConstantCache(a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::RDI3WM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  RDI3WMConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function RDI4WMConstantCache(T::Type, T2::Type)
  a021 = convert(T, 1//2)
  a031 = convert(T, -1)
  a032 = convert(T, 2)

  a121 = convert(T, 2//3)
  a131 = convert(T, 2//3)

  b021 = convert(T, (6-sqrt(6))/10)
  b031 = convert(T, (3+2*sqrt(6))/5)

  b121 = convert(T, sqrt(2//3))
  b131 = convert(T, -sqrt(2//3))

  b221 = convert(T, sqrt(2))
  b222 = convert(T, 0)
  b223 = convert(T, 0)
  b231 = convert(T, -sqrt(2))
  b232 = convert(T, 0)
  b233 = convert(T, 0)

  α1 = convert(T, 1//6)
  α2 = convert(T, 2//3)
  α3 = convert(T, 1//6)

  c02 = convert(T2, 1//2)
  c03 = convert(T2, 1)

  c12 = convert(T2, 2//3)
  c13 = convert(T2, 2//3)

  beta11 = convert(T, 1//4)
  beta12 = convert(T, 3//8)
  beta13 = convert(T, 3//8)

  beta22 = convert(T, sqrt(6)/4)
  beta23 = convert(T, -sqrt(6)/4)

  beta31 = convert(T, -1//4)
  beta32 = convert(T, 1//8)
  beta33 = convert(T, 1//8)

  beta42 = convert(T, sqrt(2)/4)
  beta43 = convert(T, -sqrt(2)/4)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  DRI1ConstantCache(a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::RDI4WM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  RDI4WMConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

@cache struct DRI1Cache{uType,randType,MType1,tabType,rateNoiseType,rateType,possibleRateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType

  _dW::randType
  _dZ::randType
  chi1::randType
  Ihat2::MType1

  tab::tabType

  g1::rateNoiseType
  g2::Vector{rateNoiseType}
  g3::Vector{rateNoiseType}

  k1::rateType
  k2::rateType
  k3::rateType

  H02::uType
  H03::uType
  H12::Vector{uType}
  H13::Vector{uType}
  H22::Vector{uType}
  H23::Vector{uType}

  tmp1::possibleRateType
  tmpg::rateNoiseType

end

function alg_cache(alg::DRI1,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = DRI1ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = [zero(noise_rate_prototype) for k=1:m]
  g3 = [zero(noise_rate_prototype) for k=1:m]
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = Vector{typeof(u)}()
  H13 = Vector{typeof(u)}()
  H22 = Vector{typeof(u)}()
  H23 = Vector{typeof(u)}()

  for k=1:m
    push!(H12,zero(u))
    push!(H13,zero(u))
    push!(H22,zero(u))
    push!(H23,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  DRI1Cache(u,uprev,_dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,k1,k2,k3,H02,H03,H12,H13,H22,H23,tmp1,tmpg)
end

function alg_cache(alg::RI1,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = RI1ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = [zero(noise_rate_prototype) for k=1:m]
  g3 = [zero(noise_rate_prototype) for k=1:m]
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = Vector{typeof(u)}()
  H13 = Vector{typeof(u)}()
  H22 = Vector{typeof(u)}()
  H23 = Vector{typeof(u)}()

  for k=1:m
    push!(H12,zero(u))
    push!(H13,zero(u))
    push!(H22,zero(u))
    push!(H23,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  DRI1Cache(u,uprev,_dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,k1,k2,k3,H02,H03,H12,H13,H22,H23,tmp1,tmpg)
end

function alg_cache(alg::RI3,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = RI3ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = [zero(noise_rate_prototype) for k=1:m]
  g3 = [zero(noise_rate_prototype) for k=1:m]
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = Vector{typeof(u)}()
  H13 = Vector{typeof(u)}()
  H22 = Vector{typeof(u)}()
  H23 = Vector{typeof(u)}()

  for k=1:m
    push!(H12,zero(u))
    push!(H13,zero(u))
    push!(H22,zero(u))
    push!(H23,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  DRI1Cache(u,uprev,_dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,k1,k2,k3,H02,H03,H12,H13,H22,H23,tmp1,tmpg)
end

function alg_cache(alg::RI5,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = RI5ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = [zero(noise_rate_prototype) for k=1:m]
  g3 = [zero(noise_rate_prototype) for k=1:m]
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = Vector{typeof(u)}()
  H13 = Vector{typeof(u)}()
  H22 = Vector{typeof(u)}()
  H23 = Vector{typeof(u)}()

  for k=1:m
    push!(H12,zero(u))
    push!(H13,zero(u))
    push!(H22,zero(u))
    push!(H23,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  DRI1Cache(u,uprev,_dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,k1,k2,k3,H02,H03,H12,H13,H22,H23,tmp1,tmpg)
end

function alg_cache(alg::RI6,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = RI6ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = [zero(noise_rate_prototype) for k=1:m]
  g3 = [zero(noise_rate_prototype) for k=1:m]
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = Vector{typeof(u)}()
  H13 = Vector{typeof(u)}()
  H22 = Vector{typeof(u)}()
  H23 = Vector{typeof(u)}()

  for k=1:m
    push!(H12,zero(u))
    push!(H13,zero(u))
    push!(H22,zero(u))
    push!(H23,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  DRI1Cache(u,uprev,_dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,k1,k2,k3,H02,H03,H12,H13,H22,H23,tmp1,tmpg)
end

function alg_cache(alg::RDI2WM,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = RDI2WMConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = [zero(noise_rate_prototype) for k=1:m]
  g3 = [zero(noise_rate_prototype) for k=1:m]
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = Vector{typeof(u)}()
  H13 = Vector{typeof(u)}()
  H22 = Vector{typeof(u)}()
  H23 = Vector{typeof(u)}()

  for k=1:m
    push!(H12,zero(u))
    push!(H13,zero(u))
    push!(H22,zero(u))
    push!(H23,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  DRI1Cache(u,uprev,_dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,k1,k2,k3,H02,H03,H12,H13,H22,H23,tmp1,tmpg)
end

function alg_cache(alg::RDI3WM,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = RDI3WMConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = [zero(noise_rate_prototype) for k=1:m]
  g3 = [zero(noise_rate_prototype) for k=1:m]
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = Vector{typeof(u)}()
  H13 = Vector{typeof(u)}()
  H22 = Vector{typeof(u)}()
  H23 = Vector{typeof(u)}()

  for k=1:m
    push!(H12,zero(u))
    push!(H13,zero(u))
    push!(H22,zero(u))
    push!(H23,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  DRI1Cache(u,uprev,_dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,k1,k2,k3,H02,H03,H12,H13,H22,H23,tmp1,tmpg)
end

function alg_cache(alg::RDI4WM,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = RDI4WMConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = [zero(noise_rate_prototype) for k=1:m]
  g3 = [zero(noise_rate_prototype) for k=1:m]
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = Vector{typeof(u)}()
  H13 = Vector{typeof(u)}()
  H22 = Vector{typeof(u)}()
  H23 = Vector{typeof(u)}()

  for k=1:m
    push!(H12,zero(u))
    push!(H13,zero(u))
    push!(H22,zero(u))
    push!(H23,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  DRI1Cache(u,uprev,_dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,k1,k2,k3,H02,H03,H12,H13,H22,H23,tmp1,tmpg)
end

# Roessler SRK for first order weak approx
struct RDI1WMConstantCache{T,T2} <: StochasticDiffEqConstantCache
  # hard-coded version
  a021::T

  #a121::T

  #a221::T

  b021::T

  #b121::T

  #b221::T

  α1::T
  α2::T

  c02::T2

  #c11::T2
  #c12::T2

  #c21::T2
  #c22::T2

  beta11::T
  #beta12::T

  #beta21::T
  #beta22::T

  #beta31::T
  #beta32::T

  #beta41::T
  #beta42::T

  #quantile(Normal(),1/6)
  NORMAL_ONESIX_QUANTILE::T
end

function RDI1WMConstantCache(::Type{T}, ::Type{T2}) where {T,T2}

  a021 = convert(T, 2//3)

  b021 = convert(T, 2//3)

  α1 = convert(T, 1//4)
  α2 = convert(T, 3//4)

  c02 = convert(T2, 2//3)

  beta11 = convert(T, 1)

  NORMAL_ONESIX_QUANTILE = convert(T, -0.9674215661017014)

  RDI1WMConstantCache(a021,b021,α1,α2,c02,beta11,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::RDI1WM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  RDI1WMConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

@cache struct RDI1WMCache{uType,randType,MType1,tabType,rateNoiseType,rateType,possibleRateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType

  _dW::randType
  _dZ::randType
  chi1::randType
  Ihat2::MType1

  tab::tabType

  g1::rateNoiseType

  k1::rateType
  k2::rateType

  H02::uType

  tmp1::possibleRateType
end

function alg_cache(alg::RDI1WM,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = RDI1WMConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  k1 = zero(rate_prototype); k2 = zero(rate_prototype)

  H02 = zero(u)

  tmp1 = zero(rate_prototype)

  RDI1WMCache(u,uprev,_dW,_dZ,chi1,Ihat2,tab,g1,k1,k2,H02,tmp1)
end

# Stratonovich sense

struct RSConstantCache{T,T2} <: StochasticDiffEqConstantCache
  # hard-coded version
  a021::T
  a031::T
  a032::T

  a131::T
  a141::T

  b031::T
  b032::T

  b121::T
  b131::T
  b132::T
  b141::T
  b142::T
  b143::T

  b221::T
  b231::T

  b331::T
  b332::T
  b341::T
  b342::T

  α1::T
  α2::T
  α3::T
  α4::T

  c02::T2
  c03::T2

  c13::T2
  c14::T2

  beta11::T
  beta12::T
  beta13::T
  beta14::T

  beta22::T
  beta23::T

  #quantile(Normal(),1/6)
  NORMAL_ONESIX_QUANTILE::T
end

function RS1ConstantCache(T::Type, T2::Type)
  a021 = convert(T, 0)
  a031 = convert(T, 1)
  a032 = convert(T, 0)

  a131 = convert(T, 1)
  a141 = convert(T, 1)

  b031 = convert(T, 1//4)
  b032 = convert(T, 3//4)

  b121 = convert(T, 2//3)
  b131 = convert(T, 1//12)
  b132 = convert(T, 1//4)
  b141 = convert(T, -5//4)
  b142 = convert(T, 1//4)
  b143 = convert(T, 2)

  b221 = convert(T, 1)
  b231 = convert(T, -1)

  b331 =  convert(T, 1//4)
  b332 =  convert(T, 3//4)
  b341 =  convert(T, 1//4)
  b342 =  convert(T, 3//4)

  α1 = convert(T, 0)
  α2 = convert(T, 0)
  α3 = convert(T, 1//2)
  α4 = convert(T, 1//2)

  c02 = convert(T2, 0)
  c03 = convert(T2, 1)

  c13 = convert(T2, 1)
  c14 = convert(T2, 1)

  beta11 = convert(T, 1//8)
  beta12 = convert(T, 3//8)
  beta13 = convert(T, 3//8)
  beta14 = convert(T, 1//8)

  beta22 = convert(T, -1//4)
  beta23 = convert(T, 1//4)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  RSConstantCache(a021,a031,a032,a131,a141,b031,b032,b121,b131,b132,b141,b142,b143,b221,b231,b331,b332,b341,b342,α1,α2,α3,α4,c02,c03,c13,c14,beta11,beta12,beta13,beta14,beta22,beta23,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::RS1,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  RS1ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function RS2ConstantCache(T::Type, T2::Type)
  a021 = convert(T, 2//3)
  a031 = convert(T, 1//6)
  a032 = convert(T, 1//2)

  a131 = convert(T, 1)
  a141 = convert(T, 1)

  b031 = convert(T, 1//4)
  b032 = convert(T, 3//4)

  b121 = convert(T, 2//3)
  b131 = convert(T, 1//12)
  b132 = convert(T, 1//4)
  b141 = convert(T, -5//4)
  b142 = convert(T, 1//4)
  b143 = convert(T, 2)

  b221 = convert(T, 1)
  b231 = convert(T, -1)

  b331 =  convert(T, 1//4)
  b332 =  convert(T, 3//4)
  b341 =  convert(T, 1//4)
  b342 =  convert(T, 3//4)

  α1 = convert(T, 1//4)
  α2 = convert(T, 1//4)
  α3 = convert(T, 1//2)
  α4 = convert(T, 0)

  c02 = convert(T2, 2//3)
  c03 = convert(T2, 2//3)

  c13 = convert(T2, 1)
  c14 = convert(T2, 1)

  beta11 = convert(T, 1//8)
  beta12 = convert(T, 3//8)
  beta13 = convert(T, 3//8)
  beta14 = convert(T, 1//8)

  beta22 = convert(T, -1//4)
  beta23 = convert(T, 1//4)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  RSConstantCache(a021,a031,a032,a131,a141,b031,b032,b121,b131,b132,b141,b142,b143,b221,b231,b331,b332,b341,b342,α1,α2,α3,α4,c02,c03,c13,c14,beta11,beta12,beta13,beta14,beta22,beta23,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::RS2,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  RS2ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

@cache struct RSCache{uType,randType,MType1,tabType,rateNoiseType,rateType,possibleRateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType

  _dW::randType
  _dZ::randType
  chi1::randType
  Ihat2::MType1

  tab::tabType

  g1::rateNoiseType
  g2::Vector{rateNoiseType}
  g3::Vector{rateNoiseType}
  g4::Vector{rateNoiseType}

  k1::rateType
  k2::rateType
  k3::rateType

  H02::uType
  H03::uType
  H12::Vector{uType}
  H13::Vector{uType}
  H14::Vector{uType}
  H22::Vector{uType}
  H23::Vector{uType}

  tmp1::possibleRateType
  tmpg::rateNoiseType

end

function alg_cache(alg::RS1,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = RS1ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = [zero(noise_rate_prototype) for k=1:m]
  g3 = [zero(noise_rate_prototype) for k=1:m]
  g4 = [zero(noise_rate_prototype) for k=1:m]
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = Vector{typeof(u)}()
  H13 = Vector{typeof(u)}()
  H14 = Vector{typeof(u)}()
  H22 = Vector{typeof(u)}()
  H23 = Vector{typeof(u)}()

  for k=1:m
    push!(H12,zero(u))
    push!(H13,zero(u))
    push!(H14,zero(u))
    push!(H22,zero(u))
    push!(H23,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  RSCache(u,uprev,_dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,g4,k1,k2,k3,H02,H03,H12,H13,H14,H22,H23,tmp1,tmpg)
end

function alg_cache(alg::RS2,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = RS2ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = [zero(noise_rate_prototype) for k=1:m]
  g3 = [zero(noise_rate_prototype) for k=1:m]
  g4 = [zero(noise_rate_prototype) for k=1:m]
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = Vector{typeof(u)}()
  H13 = Vector{typeof(u)}()
  H14 = Vector{typeof(u)}()
  H22 = Vector{typeof(u)}()
  H23 = Vector{typeof(u)}()

  for k=1:m
    push!(H12,zero(u))
    push!(H13,zero(u))
    push!(H14,zero(u))
    push!(H22,zero(u))
    push!(H23,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  RSCache(u,uprev,_dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,g4,k1,k2,k3,H02,H03,H12,H13,H14,H22,H23,tmp1,tmpg)
end


text:
struct CheckResults <: Results
    variables::Dict{Symbol, DataFrames.DataFrame}
    total_cost::Dict
    optimizer_log::Dict
    time_stamp::DataFrames.DataFrame
    check_sum::Array
end
"""This function creates the correct results struct for the context"""
function _make_results(variables::Dict,
                      total_cost::Dict,
                      optimizer_log::Dict,
                      time_stamp::DataFrames.DataFrame,
                      check_sum::Array)
    return CheckResults(variables, total_cost, optimizer_log, time_stamp, check_sum)
end
# sums the power generation across all variables into one number to verify that results
# have not been tampered with in writing and reading back
function _sum_variable_results(results::Results)
    variable_sums = []
    for (k,v) in results.variables
        variable_sums = vcat(variable_sums,sum([sum(v[!,i]) for i in 1 : size(v,2)]))
    end
    return total_sum = [sum(variable_sums)]
end

function _sum_variable_results(variables::Dict)
    variable_sums = []
    for (k,v) in variables
        variable_sums = vcat(variable_sums,sum([sum(v[!,i]) for i in 1 : size(v,2)]))
    end
    return total_sum = [sum(variable_sums)]
end


text:
function maxtrisum(values::Array{T})::T where T <: Number
    sum(unique!(sort(values))[end-2:end])
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
# Method when using Boykov-Kolmogorov as a subroutine
# Kishimoto algorithm

@traitfn function kishimoto(
    flow_graph::::IsDirected,                       # the input graph
    source::Integer,                           # the source vertex
    target::Integer,                           # the target vertex
    capacity_matrix::AbstractMatrix,      # edge flow capacities
    flow_algorithm::BoykovKolmogorovAlgorithm, # keyword argument for algorithm
    routes::Int                                # keyword argument for routes
    )
    # Initialisation
    flow, F, labels = maximum_flow(flow_graph, source, target,
    capacity_matrix, algorithm = flow_algorithm)
    restriction = flow / routes
    flow, F, labels = maximum_flow(flow_graph, source, target, capacity_matrix,
    algorithm = flow_algorithm, restriction = restriction)

    # Loop condition : approximatively not equal is enforced by floating precision
    i = 1
    while flow < routes * restriction && flow ≉ routes * restriction
        restriction = (flow - i * restriction) / (routes - i)
        i += 1
        flow, F, labels = maximum_flow(flow_graph, source, target, capacity_matrix,
        algorithm = flow_algorithm, restriction = restriction)
    end

    # End
    return flow, F, labels
end

"""
    kishimoto(flow_graph, source, target, capacity_matrix, flow_algorithm, routes)

Compute the maximum multiroute flow (for an integer number of `route`s)
between `source` and `target` in `flow_graph` with capacities in `capacity_matrix`
using the [Kishimoto algorithm](http://dx.doi.org/10.1109/ICCS.1992.255031).
Return the value of the multiroute flow as well as the final flow matrix,
along with a multiroute cut if Boykov-Kolmogorov is used as a subroutine.
"""
function kishimoto end
@traitfn function kishimoto(
    flow_graph::::IsDirected,                   # the input graph
    source::Integer,                       # the source vertex
    target::Integer,                       # the target vertex
    capacity_matrix::AbstractMatrix,  # edge flow capacities
    flow_algorithm::AbstractFlowAlgorithm, # keyword argument for algorithm
    routes::Int                            # keyword argument for routes
    )
    # Initialisation
    flow, F = maximum_flow(flow_graph, source, target,
    capacity_matrix, algorithm = flow_algorithm)
    restriction = flow / routes

    flow, F = maximum_flow(flow_graph, source, target, capacity_matrix,
    algorithm = flow_algorithm, restriction = restriction)

    # Loop condition : approximatively not equal is enforced by floating precision
    i = 1
    while flow < routes * restriction && flow ≉ routes * restriction
        restriction = (flow - i * restriction) / (routes - i)
        i += 1
        flow, F = maximum_flow(flow_graph, source, target, capacity_matrix,
        algorithm = flow_algorithm, restriction = restriction)
    end

    # End
    return flow, F
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
using Scallop.Backtest

@assert 12 == xbacktest(6)


text:
export confusion_params, confusion_matrix, class_confusion

#Confusion matrix related classes

using LinearAlgebra
using ..utils: convert_1d

"""
    confusion_params(matrix::Array{Number,2})

Return the true positives, true negatives, false positives, false negatives arrays
from the given n x n matrix. If the provided matrix is not n x n, an assertion
exception: "Given matrix is not n x n" will be raised. As a visualization for the inner
calculation of the function, [this page](https://devopedia.org/images/article/208/6541.1566280388.jpg) may be visited

"""
function confusion_params(matrix::Array{Number,2})
    @assert size(matrix)[1] == size(matrix)[2] "Given matrix is not n x n "
    tp = []; tn = []; fp = []; fn = []
    matrix_sum = sum(matrix)
     @inbounds for i in 1:size(matrix)[1]
        push!(tp, matrix[i,i])
        push!(fn, sum(matrix[i,:]) - tp[i] )
        push!(fp, sum(matrix[:,i]) -tp[i])
        push!(tn, (matrix_sum - tp[i] - fn[i] - fp[i]))
    end
    return tp, tn, fp, fn
end

"""
A struct for representing confusion matrix and related computations

## Fields
**`true_positives`** : An array that contains the true positive values of each label. For binary case,
`true_positives` is a single value. For multiclass, ith element in the array corresponds to the
`true_positives` of the ith class in the labels list.

**`true_negatives`** : An array that contains the true negative values of each label. For binary case,
`true_negatives` is a single value. For multiclass, ith element in the array corresponds to the
`true_negatives` of the ith class in the labels list.

**`false_positives`** : An array that contains the false positive values of each label. For binary case,
`false_positives` is a single value. For multiclass, ith element in the array corresponds to the
`false_positives` of the ith class in the labels list.

**false_negatives** : An array that contains the false negative values of each label. For binary case,
`false_negatives` is a single value. For multiclass, ith element in the array corresponds to the
`false_negatives` of the ith class in the labels list.

**`matrix`** : an n x n matrix where n is the length of labels. It represents the actual confusion matrix
from which true positives, true negatives, false positives and false negatives are calculated.

**`Labels`** : an array representing the labels which are used for printing and visualization purposes

**`zero division`** :
    \n\t"warn" => all NaN and Inf values are replaced with zeros and user is warned by @warn macro in the
    \tprocess

    \t"0" => all NaN and Inf values are replaced with zeros but the user is not warned by @warn macro in the
    \tprocess

    \t"1" => all NaN and Inf values are replaced with ones but the user is not warned by @warn macro in the
    \tprocess

"""
struct confusion_matrix
    true_positives::Array{Int}
    true_negatives::Array{Int}
    false_positives::Array{Int}
    false_negatives::Array{Int}
    matrix::Array{Number,2}
    Labels::Array{Union{Int,AbstractString}}
    zero_division::String
end

"""
## Constructors

```confusion_matrix(expected, predicted; keywords)```

Return a confusion matrix object constructed by the expected and predicted arrays. Expected and predicted arrays
must be of size (n,1) or or vector type. Lengths of the expected and predicted arrays must be equal; thus,
there should be a prediction and a ground truth for each classification.

## Arguments
- `expected::Vector` : Ground truth values for the classification

- `predicted::Vector` : Predictions of the classifier

## Keywords

- `labels::Vector = nothing` : List of labels to index the matrix. If nothing is given, those that appear at least in expected or predicted are used in sorted order.
- `normalize::Bool = nothing` : Determines whether or not the confusion matrix (matrix field of the produced confusion matrix) will be normalized.
- `sample_weight::Number = nothing` : Sample weights which will be filled in the matrix before confusion params function is called
- `zero_division::String = "warn"` : "warn", "0", "1", default = "warn"

See: `confusion_matrix`

## Example
\n
```julia-repl
julia> y_true = [1,1,1,2,3,3,1,2,1,1,2,1];

julia> y_pred = [1,3,2,1,2,3,1,1,2,3,2,1];

julia> x = confusion_matrix(y_true, y_pred)
\n┌ Warning: No labels provided, constructing a label set by union of the unique elements in Expected
and Predicted arrays\n

                  1      2      3
            _____________________
                  3      2      2   │1
                  2      1      0   │2     Predicted
                  0      1      1   │3

julia> y_true = ["emirhan", "knet", "metrics", "confusion", "knet", "confusion", "emirhan", "metrics", "confusion"];

julia> y_pred = ["knet", "knet", "confusion", "confusion", "knet", "emirhan", "emirhan", "knet", "confusion"];

julia> x = confusion_matrix(y_true, y_pred, labels = ["emirhan", "knet", "confusion", "metrics"])

Expected

emirhan      knet      confusion        metrics
____________________________________________________________
1              1              0              0   │emirhan
0              2              0              0   │knet
1              0              2              0   │confusion       Predicted
0              1              1              0   │metrics

```
## References
   [1] [Wikipedia entry for the Confusion matrix]
           (https://en.wikipedia.org/wiki/Confusion_matrix)
           (Note: Wikipedia and other references may use a different
           convention for axes)

See: `confusion_params`  \n
"""
function confusion_matrix(expected, predicted; labels = nothing, normalize = false, sample_weight = 0, zero_division = "warn", type = Number, suppresswarning = false)
    expected = convert_1d(expected)
    predicted = convert_1d(predicted)
    @assert length(expected) == length(predicted) "Sizes of the expected and predicted values do not match"
    @assert eltype(expected) <: Union{Int, String} &&  eltype(predicted) <: Union{Int, String} "Expected and Predicted arrays must either be integers or strings"
    @assert eltype(expected) == eltype(predicted) "Element types of Expected and Predicted arrays do not match"

    if labels != nothing; @assert length(labels) != 0 "Labels array must contain at least one value"; end;
    @assert zero_division in ["warn", "0", "1"] "Unknown zero division behaviour specification"
    if labels == nothing
        if !suppresswarning; @info "No labels provided, constructing a label set by union of the unique elements in Expected and Predicted arrays"; end
        labels = union(unique(expected),unique(predicted))
        if eltype(labels) == Int
            sort!(labels)
        end
    end
    dictionary = Dict{eltype(labels),Int}()
    for i in 1:length(labels)  #faster than passing a generator or array of pairs into dict?
         dictionary[labels[i]] = i
    end
    matrix = Array{type, 2}(undef, length(labels), length(labels))
    fill!(matrix, type(sample_weight))
    for i in 1:length(expected)
        matrix[dictionary[expected[i]],dictionary[predicted[i]]] += 1
    end
    tp, tn, fp, fn = confusion_params(matrix)
    if !suppresswarning
        if 0 in tp
            @warn "There are elements of value 0 in the true positives array. This may lead to false values for some functions"
        end
        if 0 in tn
            @warn "There are elements of value 0 in the true negatives array. This may lead to false values for some functions"
        end
        if 0 in fp
           @warn "There are elements of value 0 in the false positives array. This may lead to false values for some functions"
        end
        if 0 in fn
            @warn "There are elements of value 0 in the false negatives array. This may lead to false values for some functions"
        end
    end
    if normalize
       matrix = [round(i, digits = 3) for i in LinearAlgebra.normalize(matrix)]
    end
    return confusion_matrix(tp,tn,fp,fn,matrix,labels, zero_division)
end

confusion_params(c::confusion_matrix) = c.true_positives, c.true_negatives, c.false_positives, c.false_negatives

"""
```class_confusion(c::confusion_matrix; class_name = nothing, ith_class = nothing)```

\nReturn a binary confusion matrix for the class denoted by `class_name` or `ith_class` arguments.

## Keywords

- `ith_class::Int = nothing` : Return the binary confusion matrix of the ith class in the Labels array. This will be ignored if class_name is not `nothing`
- `class_name::Int = nothing` : Return the binary confusion matrix of the class of given value if exists in the Labels array.

## Example
\n
```julia-repl
julia> y_true = [1,1,1,2,3,3,1,2,1,1,2,1];

julia> y_pred = [1,3,2,1,2,3,1,1,2,3,2,1];

julia> x = confusion_matrix(y_true, y_pred);

julia> class_confusion(x, ith_class = 2)
2×2 Array{Int64,2}:
 1  3
 2  6

julia> class_confusion(x, class_name = 2)
2×2 Array{Int64,2}:
 1  3
 2  6
```
"""
function class_confusion(c::confusion_matrix; class_name = nothing, ith_class = nothing)
    index = check_index(c.Labels, false ,class_name = class_name, ith_class = ith_class)
    return [c.true_positives[index] c.false_positives[index]; c.false_negatives[index] c.true_negatives[index]]
end

function Base.show(io::IO, ::MIME"text/plain", c::confusion_matrix)
    printer = Int(round(size(c.matrix)[1] / 2)) +1
    label_len = maximum([length(string(i)) for i in c.Labels])[1] + 6
    label_size = length(c.Labels)
    println(io, lpad("Expected\n", printer* label_len ))
    println(io, [lpad(i,label_len) for i in c.Labels]...)
    println(io, repeat("_", length(c.Labels) * label_len))
    for i in 1:size(c.matrix)[1]
        println(io,  [lpad(string(i),label_len) for i in c.matrix[i,:]]..., "   │", c.Labels[i], i == printer ? "\tPredicted" : " ")
    end
end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# =========================================================================================
"""
    struct PipelineTask
        event::Condition
        done::Ref{Bool}
        task::Task
    end

Encapsulates all the data needed to monitor and stop a task reading the data from
Tesseract.

__Values:__

| Name   | Description
| :----- | :----------
| event  | The event to trigger when there is likely more data or the done flag is set.
| done   | The flag to set when the Tesseract process is done.
| task   | The background task that is actually executing.

__Constructors:__

    PipelineTask(event, done, task)

Create a new instance of the object.  Sets a finalizer to ensure that the background task
is terminated when this object is garbage collected.
"""
struct PipelineTask
    event::Condition
    done::Ref{Bool}
    task::Task
end

# =========================================================================================
"""
    pipeline_stop_task(
        task::PipelineTask,
        waitFor::Bool = true
    )::Nothing

Stop the background task and wait for it to terminate before returning.

__Arguments:__

| T | Name    | Default | Description
|---| :------ | :------ | :----------
| R | task    |         | The task to terminate.
| O | waitFor | `true`  | Wait for the tasks to terminate?
"""
function pipeline_stop_task(
            task::PipelineTask,
            waitFor::Bool = true
        )::Nothing
    task.done[] = true
    notify(task.event)
    if waitFor == true
        wait(task.task)
    end
    nothing
end

# =========================================================================================
"""
    pipeline_start_task(
        source::AbstractString,
        target::AbstractString
    )::PipelineTask

Start a background task that waits for the job to complete then renames the output file.

__Arguments:__

| T | Name   | Default | Description
|---| :----- | :------ | :----------
| R | source |         | The file name being generated by Tesseract.
| R | target |         | The target filename to rename the file to.
"""
function pipeline_start_task(
            source::AbstractString,
            target::AbstractString
        )::PipelineTask
    local event = Condition()
    local done  = Ref{Bool}(false)
    local task = @async begin
        while(done[] == false)
            wait(event)
        end
        mv(source, target; force = true)
    end

    return PipelineTask(event, done, task)
end

# =========================================================================================
"""
    pipeline_start_task(
        path::AbstractString,
        filename::AbstractString,
        output::AbstractString
)::PipelineTask

Start a background task that transcodes the file generate by Tesseract to UTF-8.

__Arguments:__

| T | Name     | Default | Description
|---| :------- | :------ | :----------
| R | path     |         | The directory the temporary file is being created in.
| R | filename |         | The name of the file being created by Tesseract.
| R | output   |         | The filename we want to write to.
"""
function pipeline_start_task(
            path::AbstractString,
            filename::AbstractString,
            output::AbstractString
        )::PipelineTask
    local event   = Condition()
    local done    = Ref{Bool}(false)
    local outFile = open(output; create=true, write=true)

    local task
    try
        task = start_job(event, done, path, filename, true) do job
            job_copy(job, outFile)
            close(outFile)
        end
    catch ex
        isopen(outFile) && close(outFile)
        rethrow(ex)
    end

    return PipelineTask(event, done, task)
end

# =========================================================================================
"""
    pipeline_start_task(
        path:AbstractString,
        filename::AbstractString,
        output::TessOutput{UInt8},
        transcode::Bool = false
    )::PipelineTask

Start a background task to read a file created by Tesseract and convert it to a byte array.

__Arguments:__

| T | Name      | Default | Description
|---| :-------- | :------ | :----------
| R | path      |         | The directory to use for the task.
| R | filename  |         | The filename that Tesseract will be writing to.
| R | output    |         | The object to write the result to.
| O | transcode | `false` | Should the file be transcoded from Latin-1 to UTF-8?
"""
function pipeline_start_task(
            path::AbstractString,
            filename::AbstractString,
            output::TessOutput{Vector{UInt8}},
            transcode::Bool = false
        )::PipelineTask
    local event = Condition()
    local done  = Ref{Bool}(false)
    local task  = start_job(event, done, path, filename, transcode) do job
        job_create_byte_array(job, output)
    end

    return PipelineTask(event, done, task)
end

# =========================================================================================
"""
    pipeline_start_task(
        path::AbstractString,
        filename::AbstractString,
        output::TessOutput{String},
        transcode::Bool = false
    )::PipelineTask

Start a background task to read a file created by Tesseract and convert it to a string.

__Arguments:__

| T | Name      | Default | Description
|---| :-------- | :------ | :----------
| R | path      |         | The directory to use for the task.
| R | filename  |         | The filename that Tesseract will be writing to.
| R | output    |         | The object to write the result to.
| O | transcode | `false` | Should the file be transcoded from Latin-1 to UTF-8?
"""
function pipeline_start_task(
            path::AbstractString,
            filename::AbstractString,
            output::TessOutput{String},
            transcode::Bool = false
        )::PipelineTask
    local event = Condition()
    local done  = Ref{Bool}(false)
    local task  = start_job(event, done, path, filename, transcode) do job
        job_create_string(job, output)
    end

    return PipelineTask(event, done, task)
end

# =========================================================================================
"""
    pipeline_start_task(
        path::AbstractString,
        filename::AbstractString,
        dispatch::Function,
        pageSeparator::AbstractString,
        transcode::Bool = false
    )::PipelineTask

Start a background task to read a file created by Tesseract and dispatch all the lines to
the specified function.

__Arguments:__

| T | Name      | Default | Description
|---| :-------- | :------ | :----------
| R | path      |         | The directory to use for the task.
| R | filename  |         | The filename that Tesseract will be writing to.
| R | dispatch  |         | The function to dispatch the lines to.
| R | separator |         | The string that Tesseract is using to separate pages.
| O | transcode | `false` | Should the file be transcoded from Latin-1 to UTF-8?
"""
function pipeline_start_task(
            path::AbstractString,
            filename::AbstractString,
            dispatch::Function,
            separator::AbstractString,
            transcode::Bool = false
        )::PipelineTask
    local event = Condition()
    local done  = Ref{Bool}(false)
    local task  = start_job(event, done, path, filename, transcode) do job
        job_dispatch(job, dispatch, separator)
    end

    return PipelineTask(event, done, task)
end


text:
abstract type AbstractBitstringProblem <: AbstractProblem end

import Base: length

function fitness(x::Array{Bool,1}, problem)
    fitness(BitArray(x), problem)
end

function fitness(
    x::BitArray,
    problem::T,
    fitnessvalue::Number,
    flipsposition,
) where {T<:AbstractProblem}
    y = copy(x)
    for item ∈ flipsposition
        flip!(y, item)
    end
    fitness(y, problem)
end
"""
    struct OneMax <: AbstractBitstringProblem
        n::Integer
    end

Maximization of the number of `true` in a `BitArray`. The mathematical definition of this problem is
maximizing the number of one-bits in a bitstring, which can be seen in the following formula:

```math
OneMax(x_1,\\dots,x_n)\\coloneqq\\sum_{i=1}^{n}x_i
```

## Examples
```jldoctest; setup = :(using Pidoh)
julia> problem = OneMax(4)
OneMax(4)

julia> fitness(BitArray([true,true,false,true]), problem)
3

julia> optimum(problem)
4

```
"""
struct OneMax <: AbstractBitstringProblem
    n::Integer
end
fitness(x::BitArray, problem::OneMax) = sum(x)
optimum(problem::OneMax) = problem.n
length(problem::OneMax) = problem.n

struct ZeroMax <: AbstractBitstringProblem
    n::Integer
end

fitness(x::BitArray, problem::ZeroMax) = problem.n - sum(x)
optimum(problem::ZeroMax) = 0

struct Jump <: AbstractBitstringProblem
    n::Integer
    jumpsize::Integer
end

function fitness(x::BitArray, problem::Jump)
    onemax = sum(x)
    n = problem.n
    k = problem.jumpsize
    jump = -1
    if ((onemax ≤ n - k) || (onemax == n))
        jump = k + onemax
    else
        jump = n - onemax
    end
    jump
end

optimum(problem::Jump) = problem.n + problem.jumpsize
length(problem::Jump) = problem.n

# function fitness(x::BitArray, problem::OneMax, fitnessvalue:: Number, flipsposition)
#     if length(flipsposition) == 0
#         return fitnessvalue
#     end
#     fit = fitnessvalue
#     for item ∈ flipsposition
#         @inbounds if x[item]
#             fit -= 1
#         else
#             fit += 1
#         end
#     end
#     fit
# end

# function fitness(x::BitArray, problem::Jump, fitnessvalue:: Number, flipsposition)
#     if length(flipsposition) == 0
#         return fitnessvalue
#     end
#     n=problem.n
#     k = problem.jumpsize
#
#     onemax = fitnessvalue
#
#     if fitnessvalue < k
#         onemax = n - fitnessvalue
#     else
#         onemax = fitnessvalue - k
#     end
#
#     if ((onemax ≤ n - k) || (onemax == n))
#         jump =k + onemax
#     else
#         jump = n - onemax
#     end
#
#     for item ∈ flipsposition
#         @inbounds if x[item]
#             onemax -= 1
#         else
#             onemax += 1
#         end
#     end
#
#     jump = -1
#     if ((onemax ≤ n - k) || (onemax == n))
#         jump =k + onemax
#     else
#         jump = n - onemax
#     end
#     jump
# end


text:
export Conv2d

"""
    Conv2d(T, filtersize, kwargs...)

* W: (W1,W2,...,I,O)
* X: (X1,X2,...,I,N)
* Y: (Y1,Y2,...,O,N)

where
* I: number of input channels
* O: number of output channels
* N: batch size

```julia
T = Float32
conv = Conv2d(T, (1,1,3,2))
x = CuArray{T}(5,5,3,3)
y = conv(x)
```
"""
mutable struct Conv2d
    w::Var
    pad::NTuple{2,Int}
    stride::NTuple{2,Int}
    dilation::NTuple{2,Int}
end

function Conv2d(::Type{T}, filtersize::Tuple;
    pad=0, stride=1, dilation=1, init_w=Xavier(), init_b=Fill(0)) where T

    N = length(filtersize) - 2
    isa(pad,Int) && (pad = ntuple(_ -> pad, N))
    isa(stride,Int) && (stride = ntuple(_ -> stride, N))
    isa(dilation,Int) && (dilation = ntuple(_ -> dilation, N))

    w = init_w(T, filtersize...)
    Conv(param(w), pad, stride, dilation)
end

function (f::Conv2d)(x::Var)
    configure!(f.w, x)
    y, work = f(x.data)
    Var(y, (f,x,work))
end

function (f::Conv2d)(x::Array{T,4}) where T
    hdims = ntuple(2) do d
        k = (size(w,d)-1) * dilation[d] + 1
        1 + (size(x,d) + 2pad[d] - k) ÷ stride[d]
    end
    h = similar(x, hdims)
    #im2col(x, h, )
    y = linear(h, f.W.data, f.b.data)
    Var(y, (f,x,f.W,f.b,batchdims,h))
end

function (f::Conv2d)(x::CuArray)
    CUDNN.convolution(f.w.data, x, f.pad, f.stride, f.dilation)
end

function addgrad!(y::Var, f::Conv2d, x::Var, work)
    isvoid(x.grad) && return
    ∇conv!(y.grad, f, x.data, x.grad, work)
end

function ∇conv!(gy::CuArray, f::Conv2d, x::CuArray, gx, convdesc)
    w, gw = f.w.data, f.w.grad
    isvoid(gw) || CUDNN.∇convolution_filter!(convdesc, x, gy, gw)
    isvoid(gx) || CUDNN.∇convolution_data!(convdesc, w, gy, gx)
end

function im2col{T}(x::Array{T,4}, kernel::NTuple{4,Int}, pad::NTuple{2,Int}, stride::NTuple{2,Int}, dilation::NTuple{2,Int})
    width, height, channels, num = size(x)
    kernel_w, kernel_h = kernel[1], kernel[2]
    pad_w, pad_h = pad
    stride_w, stride_h = stride

    y_w = (width + 2pad_w - (dilation_w * (kernel_w-1) + 1)) ÷ stride_w + 1
    y_h = (height + 2pad_h - (dilation_h * (kernel_h-1) + 1)) ÷ stride_h + 1
    y = similar(x, y_w, y_h, kernel[4], num)

    height_col = div(height + 2pad_h - kernel_h, stride_h) + 1
    width_col = div(width + 2pad_w - kernel_w, stride_w) + 1
    channels_col = channels * kernel_h * kernel_w

    for c = 0:channels_col-1
        w_offset = c % kernel_w
        h_offset = div(c, kernel_w) % kernel_h
        c_im = div(c, kernel_h * kernel_w) # channel
        for h = 0:height_col-1
            for w = 0:width_col-1
                h_pad = h*stride_h - pad_h + h_offset
                w_pad = w*stride_w - pad_w + w_offset
                if (h_pad >= 0 && h_pad < height && w_pad >= 0 && w_pad < width)
                    @inbounds col[1 + (c*height_col+h) * width_col + w] =
                        img[1 + (c_im * height + h_pad) * width + w_pad]
                else
                    @inbounds col[1 + (c*height_col+h) * width_col + w] = 0
                end
            end
        end
    end
end

function col2im{T}(col::Array{T}, img::Array{T}, width::Int, height::Int, channels::Int,
    kernel::NTuple{2,Int}, pad::NTuple{2,Int}, stride::NTuple{2,Int})

    kernel_w, kernel_h = kernel
    pad_w, pad_h = pad
    stride_w, stride_h = stride

    height_col = div(height + 2pad_h - kernel_h, stride_h) + 1
    width_col = div(width + 2pad_w - kernel_w, stride_w) + 1
    channels_col = channels * kernel_h * kernel_w

    fill!(img, 0)
    for c = 0:channels_col-1
        w_offset = c % kernel_w
        h_offset = div(c, kernel_w) % kernel_h
        c_im = div(c, kernel_w * kernel_h)
        for h = 0:height_col-1
            for w = 0:width_col-1
                h_pad = h * stride_h - pad_h + h_offset
                w_pad = w * stride_w - pad_w + w_offset
                if h_pad >= 0 && h_pad < height && w_pad >= 0 && w_pad < width
                    @inbounds img[1 + (c_im * height + h_pad) * width + w_pad] +=
                        col[1 + (c * height_col + h) * width_col + w]
                end
            end
        end
    end
end


text:
import ..SUN_generator:Generator,lie2matrix!,matrix2lie!

abstract type TA_Gaugefields{NC,Dim}  <: AbstractGaugefields{NC,Dim} #Traceless antihermitian matrix
end

include("./4D/TA_gaugefields_4D.jl")
include("./2D/TA_gaugefields_2D.jl")

function Base.:*(x::Array{<: TA_Gaugefields{NC,Dim},1},y::Array{<: TA_Gaugefields{NC,Dim},1})  where {NC,Dim}
    s = 0.0
    for μ=1:Dim
        s += x[μ]*y[μ]
    end

    return s
end

function initialize_TA_Gaugefields(U::Array{<:AbstractGaugefields{NC,Dim},1}) where {NC,Dim}
    F = Array{TA_Gaugefields{NC,Dim},1}(undef,Dim)
    for μ=1:Dim
        F[μ] = initialize_TA_Gaugefields(U[μ])
    end
    return F
end

function initialize_TA_Gaugefields(u::AbstractGaugefields{NC,Dim}) where {NC,Dim}
    mpi = u.mpi
    if mpi
        if Dim==4
            return TA_Gaugefields_4D_mpi(u)
        elseif Dim==2
            return TA_Gaugefields_2D_mpi(u)
            #error("Dim = $Dim is not supoorted")
        else
            error("Dim = $Dim is not supoorted")
        end

    else
        if Dim==4
            return TA_Gaugefields(NC,u.NX,u.NY,u.NZ,u.NT)
        elseif Dim==2
            return TA_Gaugefields(NC,u.NX,u.NT)
        else
            error("Dim = $Dim is not supoorted")
        end
    end
end

#function gauss_distribution!(p::T) where T<: TA_Gaugefields
#    error("gauss_distribution!(p) is not implemented in type $(typeof(p)) ")
#end

function gauss_distribution!(p::T;σ=1.0) where T<: TA_Gaugefields
    error("gauss_distribution!(p) is not implemented in type $(typeof(p)) ")
end

function gauss_distribution!(p::Array{<:TA_Gaugefields{NC,Dim},1};σ=1.0) where {NC,Dim}
    for μ=1:Dim
        gauss_distribution!(p[μ],σ=σ)
    end
end

function Base.setindex!(x::T,v,i...)  where T<: TA_Gaugefields
    error("setindex! is not implemented in type $(typeof(x)) ")
    x.a[i...] = v
end

function Base.getindex(x::T,i...) where T<: TA_Gaugefields
    error("setindex! is not implemented in type $(typeof(x)) ")
    return x.a[i...]
end

function Base.similar(U::T) where T<: TA_Gaugefields
    error("similar! is not implemented in type $(typeof(U)) ")
end

function Traceless_antihermitian_add!(U::T,factor,temp1) where T <: TA_Gaugefields
    error("Traceless_antihermitian_add! is not implemented in type $(typeof(U)) ")
end

function Traceless_antihermitian!(vout::T,vin::T) where T <: TA_Gaugefields
    error("Traceless_antihermitian! is not implemented in type $(typeof(vout)) ")
end


text:
struct Space{T}
    net::MetNet
    box::BoxGrid
    vec::Vector{Vector{T}}
    iter::Base.Iterators.Zip

    function Space{elT}(net::MetNet, freeids::Vector{Symbol}, δ::Int;
            nthrs = nthreads(), filter::Function = (net, v) -> true
        ) where {elT  <: AbstractFloat}

        box = BoxGrid(net, freeids, δ)
        vec = subspace(net, box, freeids; nthrs, filter)
        any(length.(vec) .!= length(first(vec))) && error("All dimensions must have equal length")
        new(net, box, vec, zip(vec...))
    end

    Space(net::MetNet, freeids::Vector{Symbol}, δ::Int; kwargs...) = Space{Float64}(net, freeids, δ; kwargs...)

end

## ------------------------------------------------------
function Vi(V::Space, free::Symbol)
    ci = free_ci(V, free)
    V.vec[ci]
end

## ------------------------------------------------------
free_ids(V::Space) = V.box.dims_ids
free_ids(V::Space, free) = V.box.dims_ids[free_ci(V, free)]
free_ci(V::Space, freeid::Symbol) = findfirst(isequal(freeid), free_ids(V))
free_ci(V::Space, free::Int) = findfirst(isequal(free), eachindex(free_ids(V)))
free_ci(V::Space, frees) = free_ci.([V], frees)
Base.vec(V::Space) = V.vec

## ------------------------------------------------------
# Iterator Interface
Base.IteratorSize(::Space) = Base.HasLength()
Base.IteratorEltype(::Space{T}) where {T} = Base.HasEltype()
Base.eltype(V::Space{T}) where {T} = Vector{T}
Base.length(V::Space) = length(first(V.vec))
Base.size(V::Space) = tuple(length(V))
Base.size(V::Space, dim) = size(V)[dim]
Base.iterate(V::Space) = iterate(V.iter)
Base.iterate(V::Space, state) = iterate(V.iter, state)

## ------------------------------------------------------
# Array interface
Base.eachindex(V::Space) = eachindex(first(V.vec))
Base.firstindex(V::Space) = firstindex(first(V.vec))
Base.lastindex(V::Space) = lastindex(first(V.vec))

Base.getindex(V::Space, i::Int) = getindex.(V.vec, i)
Base.getindex(V::Space, i) = getindex.([V], i)

## ------------------------------------------------------
function Base.show(io::IO, V::Space{T}) where {T}
    println(io, "Space{", T, "}")
    println(io, "frees: [", join(string.(free_ids(V)), ", "), "]")
    println(io, "len: ", length(V))
    println(io, V.net)
end
Base.show(V::Space) = show(stdout, V)

## ------------------------------------------------------
Base.empty!(V::Space) = (empty!(V.vec); empty!(V.sortperms); V)


text:
using BenchmarkTools
using DataFrames
using DelimitedFiles
using CSV
using XLSX
using Downloads

# downloading files

P = Downloads.download("https://raw.githubusercontent.com/nassarhuda/easy_data/master/programming_languages.csv",
    "data/programming_languages.csv")

matrix, header = readdlm("data/programming_languages.csv",',';header=true)
# writedlm("data/programminglanguages_dlm.txt", matrix, '-')

data = CSV.read("data/programming_languages.csv", DataFrame)

@show data[1:10,:]
@show names(data)
# @show typeof(data)
# @show data.year
@show describe(data)

# write a csv
# @btime CSV.write("data/programming_languages.csv", DataFrame(matrix, :auto))

G = XLSX.readtable("data/zillow_data_download_april2020.xlsx","Sale_counts_city")
data = DataFrame(G...)  #splat unwraps the array G[1], G[2]

@show data[1:10, 1:5]

max_stars_repo_name: lejinvarghese/julia-datascience

max_issues_repo_name: lejinvarghese/julia-datascience

max_forks_repo_name: lejinvarghese/julia-datascience


text:
@init @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin
    using .CUDA
    import .CUDA.CUBLAS

    for (fname, elty) in
        ((:cublasDgemmStridedBatched,:Float64),
         (:cublasSgemmStridedBatched,:Float32),
         (:cublasHgemmStridedBatched,:Float16),
         (:cublasZgemmStridedBatched,:ComplexF64),
         (:cublasCgemmStridedBatched,:ComplexF32))
        @eval begin

            @inline function gemm_strided_batched!(
                transA::Char, transB::Char,
                alpha::($elty), A::CuArray{$elty, 3},
                B::CuArray{$elty, 3}, beta::($elty),
                C::CuArray{$elty, 3})
                return CUBLAS.gemm_strided_batched!(transA, transB, alpha, A, B, beta, C)
            end

            @inline function gemm_strided_batched_impl!(
                transA::Char, transB::Char,
                m::Int, n::Int, k::Int,
                alpha::($elty), A::CuArray{$elty}, lda::Int, strideA::Int,
                B::CuArray{$elty}, ldb::Int, strideB::Int, beta::($elty),
                C::CuArray{$elty}, ldc::Int, strideC::Int, batchCount::Int)

                CUBLAS.$fname(CUBLAS.handle(),
                              transA, transB, m, n, k,
                              alpha, A, lda, strideA,
                              B, ldb, strideB, beta,
                              C, ldc, strideC, batchCount)

                return C
            end

        end
    end

end


text:
@testset "Map" begin
  include("Map/FqPolyRing.jl")
end