text:
"""
    rodrigues(nv, Θ)

Create a rotation matrix from a normalized axis and an angle (in radian).
"""
function rodrigues(nv, Θ)
    et = eltype(nv)
    #R = nv*nv' + cos(Θ).*(Matrix{et}(I, 3,3) - nv*nv') + sin(Θ).*crossprodtensor(nv)
    R = nv*nv' + cos(Θ).*(Matrix(I, 3,3) - nv*nv')
    pluscrossprod!(R, sin(Θ), nv)
    return R
end

"""
    rodrigues(nv, Θ)

Create a rotation matrix from a normalized axis and an angle (in radian).
"""
function rodrigues(nv::SA, Θ) where {SA<:StaticArray}
    #R = nv*nv' + cos(Θ).*(Matrix{et}(I, 3,3) - nv*nv') + sin(Θ).*crossprodtensor(nv)
    R = MMatrix{3,3}(nv*nv' + cos(Θ).*(Matrix(I, 3,3) - nv*nv'))
    pluscrossprod!(R, sin(Θ), nv)
    return SMatrix{3,3}(R)
end

"""
    pluscrossprod!(A, value, v)

Add a crossproduct tensor to `A`.
Equals to: `A+value .* crossprodtensor(v)`
"""
function pluscrossprod!(A, value, v)
    A[1,2] -= value*v[3]
    A[1,3] += value*v[2]

    A[2,1] += value*v[3]
    A[2,3] -= value*v[1]

    A[3,1] -= value*v[2]
    A[3,2] += value*v[1]

    return A
end

"""
    crossprodtensor(v)

Create the cross product tensor of 3 dimensional vector.
"""
function crossprodtensor(v)
    [0 -v[3] v[2];
    v[3] 0 -v[1];
    -v[2] v[1] 0]
end

"""
    rodriguesrad(nv, ϑ)

Create a rotation matrix from an axis (`nv`) and an angle (`ϑ`) in radians.
"""
function rodriguesrad(nv, ϑ)
    nvn = normalize(nv)
    return rodrigues(nvn, ϑ)
end

"""
    rodriguesdeg(nv, ϑ)

Create a rotation matrix from an axis (`nv`) and an angle (`ϑ`) in degrees.
"""
function rodriguesdeg(nv, ϑ)
    nvn = normalize(nv)
    return rodrigues(nvn, deg2rad(ϑ))
end

## Collection of useful functions.

"""
    arbitrary_orthogonal(vec)

Create an arbitrary orthogonal vector to `vec`.
"""
function arbitrary_orthogonal(vec)
    @assert size(vec,1) == 3 "Implemented only for 3 dimensional."
    v = normalize(vec)
    b0 = (v[1] <  v[2]) && (v[1] <  v[3]);
    b1 = (v[2] <= v[1]) && (v[2] <  v[3]);
    b2 = (v[3] <= v[1]) && (v[3] <= v[2]);
    rv = normalize!([convert(Int, b0), convert(Int, b1), convert(Int, b2)])
    return convert(typeof(vec), cross(v, rv))
end

"""
    arbitrary_orthogonal2(vec)

Create an arbitrary orthogonal vector to `vec`.
"""
function arbitrary_orthogonal2(vec)
    @assert size(vec,1) == 3 "Implemented only for 3 dimensional."
    v = normalize(vec)
    rv = normalize(rand(3))
    return convert(typeof(vec), cross(v, rv))
end

"""
    isparallel(v1, v2, alpharad)

Check if `v1` and `v2` are parallel in an `alpharad` angle
    and point towards the same direction.

The vectors considered to be normalized.
"""
function isparallel(v1, v2, alpharad)
    dot(v1, v2) > cos(alpharad)
end

"""
    findAABB(points)

Find the axis aligned minimum bounding box of the given points.
"""
function findAABB(points)
    minV = convert(Array, deepcopy(points[1]))
    maxV = convert(Array, deepcopy(points[1]))
    for i in 1:length(points)
        a = points[i]
        for j in 1:length(minV)
            minV[j] = minV[j] > a[j] ? a[j] : minV[j]
            maxV[j] = maxV[j] < a[j] ? a[j] : maxV[j]
        end
    end
    return convert(eltype(points), minV), convert(eltype(points), maxV)
end

"""
    chopzaxis(points)

Chop the 3rd element of every point.
"""
function chopzaxis(points)
    @assert length(points[1]) == 3 "Implemented only for 3 long vectors."
    [SVector(a[1], a[2]) for a in points]
end

"""
    unitdisk2square(p)

Transform point on unit disk to unit sqaure.

Maps a unit circle to a unit sqaure [0,1]^2.
Source: shirley1997 - A low distortion map between disk and square
"""
function unitdisk2square(p)
    r = norm(p)
    phi = atan(p[2], p[1])
    if phi < -π/4
        phi += 2*π # in range [-pi/4,7pi/4]
    end
    if phi < π/4 # region 1
        a = r
        b = phi*a/(π/4)
    elseif phi < 3*π/4 # region 2
        b = r
        a = -(phi-π/2)*b/(π/4)
    elseif phi < 5*π/4 # region 3
        a = -r
        b = (phi - π)*a/(π/4)
    else # region 4
        b = -r
        a = -(phi-3*π/2)*b/(π/4)
    end
    # TODO:
    # decision: squre: [-1,1]^2 or [0,1]^2
    convert(typeof(p), [(a+1)/2, (b+1)/2])
end

"""
    smallestdistance(points)

Find the smallest distance between the points.
"""
function smallestdistance(points)
    @assert length(points) > 1 "At least two point is needed for that."
    ld = norm(points[2]-points[1])
    for i in eachindex(points)
        for j in eachindex(points)
            if i!=j
                d = norm(points[i]-points[j])
                ld = d < ld ? d : ld
            end
        end
    end
    ld
end

"""
    minmaxdistance(points)

Find the smallest, largest and average distance between the points.
"""
function minmaxdistance(points)
    @assert length(points) > 1 "At least two point is needed for that."
    ld = norm(points[2]-points[1])
    maxid = norm(points[2]-points[1])
    sumd = 0.
    sumi = 0
    for i in eachindex(points)
        for j in eachindex(points)
            if i!=j
                d = norm(points[i]-points[j])
                sumd += d
                sumi += 1
                ld = d < ld ? d : ld
                maxid = d > maxid ? d : maxid
            end
        end
    end
    avg_d = sumd/sumi
    return (mind=ld, maxd=maxid, avgd=avg_d)
end

"""
    not0minmaxdistance(points)

Find the smallest and largest distance between the points.
Zero distance is not taken into account.
"""
function not0minmaxdistance(points)
    @assert length(points) > 1 "At least two point is needed for that."
    ld = norm(points[2]-points[1])
    maxid = norm(points[2]-points[1])
    for i in eachindex(points)
        for j in eachindex(points)
            if i!=j
                d = norm(points[i]-points[j])
                d == 0 && continue
                ld = d < ld ? d : ld
                maxid = d > maxid ? d : maxid
            end
        end
    end
    return (mind=ld, maxd=maxid)
end

"""
    prob(n, s, N, k)

The probability of successful detection of a shape sized `n`,
from a point cloud size of `N`, with `k` size of minimal sets and `s` draws.

# Arguments:
- `n`: size of the shape.
- `s`: number of candidates that have been drawn.
- `N`: size of the point cloud.
- `k`: size of the minimal set required to define a shape candidate.
"""
prob(n, s, N, k) = 1-(1-(n/N)^k)^s

"""
    havesameelement(A, B)

Return true if `A` and `B` have at least one common element.
"""
function havesameelement(A, B)
    for a in A
        for b in B
            a == b && return true
        end
    end
    false
end

"""
    allisdifferent(a::AbstractArray{T, 1}) where {T}

Check if all the elements of an array are different.
Supports only one dimensional arrays.
"""
function allisdifferent(a::AbstractArray{T, 1}) where {T}
    #TODO: make it more abstract with:
    #https://docs.julialang.org/en/v1/base/iterators/
    size(a, 1) < 2 && return true
    for i in 2:size(a, 1)
        for j in 1:i-1
            a[i] == a[j] && return false
        end
    end
    return true
end

function chooseS(A, k)
    symtoint = Dict(:lengthC=>1, :allcand=>2, :nofminset=>3)
    return A[symtoint[k]]
end

"""
    defaultiterationparameters(shape_types)

Construct a named tuple with the default iteration parameters.
`shape_types` is an array of `FittedShape`s,
that controls which primitives you want to fit to the point cloud.

# Examples

```julia-repl
julia> RANSAC.defaultiterationparameters([FittedPlane])
(iteration = (drawN = 3, minsubsetN = 15, prob_det = 0.9,
shape_types = UnionAll[FittedPlane], τ = 900, itermax = 1000,
extract_s = :nofminset, terminate_s = :nofminset),)

julia> RANSAC.defaultiterationparameters([FittedPlane, FittedSphere, FittedCone])
(iteration = (drawN = 3, minsubsetN = 15, prob_det = 0.9,
shape_types = UnionAll[FittedPlane, FittedSphere, FittedCone], τ = 900, itermax = 1000,
extract_s = :nofminset, terminate_s = :nofminset),)
```

# Implementation
- `drawN`: number of points to be sampled (length of a minimal subset).
- `minsubsetN`: number of minimal sets sampled in one iteration.
- `prob_det`: probability of detection.
- `τ`: minimal shape size.
- `itermax`: maximum number of iteration.
- `shape_types`: shapes that are fitted to the point cloud (array of types).
- `extract_s`, `terminate_s`: they are for easier testing, do not delete or modify it.
"""
function defaultiterationparameters(shape_types)
    # `drawN`: number of points to be sampled (length of a minimal subset)
    # `minsubsetN`: number of minimal sets sampled in one iteration
    # `prob_det`: probability of detection
    # `τ`: minimal shape size
    # `itermax`: maximum number of iteration
    # `shape_types`: shapes that are fitted to the point cloud (array of types)
    # `extract_s`, `terminate_s`:
        # track the number of candidates for the probabilities
        # which to use to terminate the algorithm
    # 1. :lengthC - number of candidates in a given iteration
    # 2. :allcand - number of candidates that have been ever scored
    # 3. :nofminset - number of minimal sets that have been drawn so far
    ip = (drawN=3, minsubsetN=15, prob_det=0.9, shape_types=shape_types, τ=900, itermax=1000, extract_s=:nofminset, terminate_s=:nofminset)
    return (iteration=ip,)
end

"""
    defaultcommonparameters()

Construct a `NamedTuple` with the default common parameters.

# Examples
```julia-repl
julia> defaultcommonparameters()
(common = (collin_threshold = 0.2, parallelthrdeg = 1.0),)
```

# Implementation
This section describes the role of the common parameters.
- `collin_threshold`: 3 points can be nearly collinear,
    in some cases they must be filtered. See the code of:
    ` fit(::Type{FittedPlane}, p, n, params)`.
- `parallelthrdeg`: threshold for two vectos being parallel, in degrees.
    If `abs(dot(a,b))>cosd(parallelthrdeg)`, `a` and `b` are considered to be parallel.
"""
function defaultcommonparameters()
    #`collin_threshold`: threshold of points being collinear
    # `parallelthrdeg`: threshold of two vectors being parallel (in degrees)
    cp = (collin_threshold = 0.2, parallelthrdeg = 1.)
    return (common=cp,)
end

"""
    defaultparameters(shape_types::Vector{T}) where {T}

Construct a `NamedTuple` with the given shape types and the default parameters.

# Examples
```julia-repl
julia> defaultparameters([FittedSphere, FittedPlane])
(iteration = (drawN = 3, minsubsetN = 15, prob_det = 0.9,
shape_types = UnionAll[FittedSphere, FittedPlane], τ = 900, itermax = 1000,
extract_s = :nofminset, terminate_s = :nofminset),
common = (collin_threshold = 0.2, parallelthrdeg = 1.0),
sphere = (ϵ = 0.3, α = 0.08726646259971647, sphere_par = 0.02),
plane = (ϵ = 0.3, α = 0.08726646259971647))
```
"""
function defaultparameters(shape_types::Vector{T}) where {T}
    def_iter_pars = defaultiterationparameters(shape_types)
    def_iter_pars = merge(def_iter_pars, defaultcommonparameters())
    def_pars = defaultshapeparameters.(shape_types)
    for a in def_pars
        def_iter_pars = merge(def_iter_pars, a)
    end
    def_iter_pars
end

"""
    ransacparameters(p::T=DEFAULT_PARAMETERS; kwargs...) where {T<:NamedTuple}

Construct a `NamedTuple` based on a previous one,
defaulting to `DEFAULT_PARAMETERS` and override it with the kwargs.
Check the docs and examples for more.

# Examples
```julia-repl
julia> p1 = ransacparameters()
(iteration = (drawN = 3, minsubsetN = 15, prob_det = 0.9),
plane = (ϵ = 0.3, α = 0.08726646259971647),
cone = (ϵ = 0.3, α = 0.08726646259971647, minconeopang = 0.03490658503988659),
cylinder = (ϵ = 0.3, α = 0.08726646259971647),
sphere = (ϵ = 0.3, α = 0.08726646259971647, sphere_par = 0.1))

julia> p2 = ransacparameters(p1; sphere=(ϵ=0.9, α=deg2rad(1),), plane=(ϵ=1.0,))
(iteration = (drawN = 3, minsubsetN = 15, prob_det = 0.9),
plane = (ϵ = 1.0, α = 0.08726646259971647),
cone = (ϵ = 0.3, α = 0.08726646259971647, minconeopang = 0.03490658503988659),
cylinder = (ϵ = 0.3, α = 0.08726646259971647),
sphere = (ϵ = 0.9, α = 0.017453292519943295, sphere_par = 0.1))
```
"""
function ransacparameters(p::T=DEFAULT_PARAMETERS; kwargs...) where {T<:NamedTuple}
    newp = p
    for a in keys(kwargs)
        oldpar = get(p, a, kwargs[a])
        newpar = merge(oldpar, kwargs[a])
        newp = merge(newp, (a=>newpar,))
    end
    return newp
end

"""
    ransacparameters(p::AbstractArray; kwargs...)

Construct a `NamedTuple` for a given types of shapes
using [`defaultparameters`](@ref) and override it with the kwargs.
Check the docs and examples for more.

# Examples
```julia-repl
julia> p1 = ransacparameters([FittedSphere, FittedCylinder])
(iteration = (drawN = 3, minsubsetN = 15, prob_det = 0.9,
shape_types = UnionAll[FittedSphere, FittedCylinder], τ = 900, itermax = 1000,
extract_s = :nofminset, terminate_s = :nofminset),
common = (collin_threshold = 0.2, parallelthrdeg = 1.0),
sphere = (ϵ = 0.3, α = 0.08726646259971647, sphere_par = 0.02),
cylinder = (ϵ = 0.3, α = 0.08726646259971647))

julia> p2 = ransacparameters([FittedSphere, FittedCylinder], sphere=(ϵ=0.01,), cylinder=(α=0.02,))
(iteration = (drawN = 3, minsubsetN = 15, prob_det = 0.9,
shape_types = UnionAll[FittedSphere, FittedCylinder], τ = 900, itermax = 1000,
extract_s = :nofminset, terminate_s = :nofminset),
common = (collin_threshold = 0.2, parallelthrdeg = 1.0),
sphere = (ϵ = 0.01, α = 0.08726646259971647, sphere_par = 0.02),
cylinder = (ϵ = 0.3, α = 0.02))
```
"""
function ransacparameters(p::AbstractArray; kwargs...)
    dpars = defaultparameters(p)
    return ransacparameters(dpars; kwargs...)
end

"""
    push2candidatesandlevels!(candidates, candidate, levels, current_level)

Push `candidate` to `candidates` and `current_level` to `levels`.
Handle the case, when `candidate` is an array of candidates.
Used in [`forcefitshapes!`](@ref).
"""
function push2candidatesandlevels!(candidates, candidate, levels, current_level)
    push!(candidates, candidate)
    push!(levels, current_level)
end

function push2candidatesandlevels!(candidates, candidate::AbstractArray, levels, current_level)
    append!(candidates, candidate)
    append!(levels, fill(current_level, size(candidate)))
end

"""
    setfloattype(nt, T)

Convert every `Real` but not `Integer` value to type `T` in a `NamedTuple` recursively.
"""
function setfloattype(nt, T)
    Keys = []
    Values = []
    for (k,v) in pairs(nt)
        push!(Keys, k)
        v = nt[k]
        if v isa NamedTuple
            push!(Values, setfloattype(v,T))
        elseif (v isa Number) && ! (v isa Integer)
            # v is a Real or a float
            push!(Values, convert(T, v))
        else
            push!(Values, v)
        end
    end
    return (; zip(Keys,Values)...)
end


text:
using DocThemeIndigo
using Test

@testset "install" begin
    indigo = DocThemeIndigo.install(DocThemeIndigo)
    path = joinpath(dirname(dirname(pathof(DocThemeIndigo))), "docs", "src", indigo)
    @test ispath(path)
end

ext: jl


text:
using Pkg; Pkg.activate("GraphFusedElasticNet.jl")
using LinearAlgebra
using Distributions, StatsFuns
using Plots, Images
using Random
using YAML, DotMaps
using GraphFusedElasticNet
using NPZ

## read config
cfg = DotMap(YAML.load_file("simulations/geometries/conf.yaml"))
Random.seed!(cfg.seed)

## cell types

nc, nr = cfg.grid.size
signal = zeros(nr, nc)
attempts = zeros(nr, nc)
successes = zeros(nr, nc)
miss = zeros(Bool, nr, nc)

## background
p = rand(Uniform(cfg.bg.prob...))
N = cfg.bg.obs
signal .= p
attempts .= N
for j in 1:nc, i in 1:nr
    successes[i, j] = rand(Binomial(N, p))
end

## triangles
N = cfg.signal.obs
for k in 1:cfg.signal.triangle.n
    a, b = cfg.signal.triangle.width
    width = rand(DiscreteUniform(a, b))
    a, b = cfg.signal.triangle.height
    height = rand(DiscreteUniform(a, b))
    # rel_height = rand(Uniform(cfg.signal.triangle.rel_height...))
    # height = ceil(Int, width * rel_height)
    c0 = rand(DiscreteUniform(1, nc - width))
    r0 = rand(DiscreteUniform(1, nr - height))
    p = rand(Uniform(cfg.signal.prob...))
    for j in c0:(c0 + width), i in r0:(r0 + height)
        m = 0.5 * height / width
        b = r0 + 0.5 * height
        if abs(i - b) ≤ m * (j - c0)
            signal[i, j] = p
            attempts[i, j] = 1
            successes[i, j] = rand(Binomial(N, p))
        end
    end
end

##
plot(Gray.(1.0 .- signal))
xlabel!("time")
ylabel!("space")

##
savefig("simulations/geometries/geometries_truth.pdf")

## ellipses
for k in 1:cfg.miss.ellipse.n
    width = rand(DiscreteUniform(cfg.miss.ellipse.width...))
    height = rand(DiscreteUniform(cfg.miss.ellipse.height...))
    r0 = rand(DiscreteUniform(1, nr - width))
    c0 = rand(DiscreteUniform(1, nc - height))
    miss_prob = rand(Uniform(cfg.miss.prob...))
    θ = 2π * rand(Uniform(cfg.miss.ellipse.angle...))
    R = [cos(θ) sin(θ); sin(θ) -cos(θ)]
    for j in 1:nc, i in 1:nr
        point = R * [i - r0, j - c0]
        if (point[1] / height)^2 + (point[2] / width)^2 ≤ 1.0
            j
            if rand() < miss_prob
                attempts[i, j] = 0
                successes[i, j] = 0
                miss[i, j] = true
            end
        end
    end
end

## visualize
u = 1.0 .- successes ./ (attempts .+ 1e-6);
z = [
    miss[i, j] ? colorant"#FF9AA2" : RGB(u[i,j], u[i, j], u[i,j])
    for i in 1:nr, j in 1:nc
];

plt = plot(z)
xlabel!("time")
ylabel!("space")

##
savefig("simulations/geometries/geometries_data.pdf")

## spatial neighbor correlation
x = Float64[]
tnb = Float64[]
snb = Float64[]
dnb = Float64[]
for j in 2:nc, i in 2:nr
    push!(x, signal[i, j])
    push!(tnb, signal[i, j - 1])
    push!(snb, signal[i - 1, j])
    push!(dnb, signal[i - 1, j - 1])
end

##
X = [x tnb snb dnb]
cor(X)

## now the smoothing
ptr, brks = grid_trails(nr, nc)
λ1 = 0.75
λ2 = 1.5
lambdasl1 = fill(λ1, size(ptr))
lambdasl2 = fill(λ2, size(ptr))

## MAP estimate
modelopts = Dict{Symbol, Any}(
    :admm_balance_every => 50,
    :admm_init_penalty => 10.0,
    :admm_residual_balancing => true,
    :admm_adaptive_inflation => false,
    :reltol => 1e-3,
    :admm_min_penalty => 0.1,
    :admm_max_penalty => 100.0,
    :abstol => 1e-8,
    :save_norms => true,
    :save_loss => true
)

fitopts = Dict{Symbol, Any}(
    :walltime => 300.0,
    :parallel => true
)

model = BinomialGFEN(ptr, brks, lambdasl1, lambdasl2; modelopts...)

ϵ = 1e-8
succ = vec(successes) .+ ϵ
att = vec(attempts) .+ 2ϵ
@time fit!(model, succ, att; fitopts...)

##
plot(log10.(model.prim_norms))

##
plot(log10.(model.dual_norms))

##
ω = 1. .- logistic.(model.beta)
ω = reshape(ω, nr, nc)
# npzwrite("simulations/geometries/map.npy", Float32.(ω))

# ω = Float64.(npzread("simulations/geometries/map.npy"))
plot(Gray.(ω))
xlabel!("time")
ylabel!("space")

##

# savefig("simulations/geometries/solution_lasso_1-0.pdf")
# savefig("simulations/geometries/solution_ridge_2-0.pdf")
savefig("simulations/geometries/solution_enet_0-75_1-5.pdf")

## uncertainty quantification
edges = Tuple{Int, Int}[]
λ1 = 0.75
λ2 = 1.5
idx = LinearIndices((nr, nc))
for j in 1:nc, i in 1:nr
    (i > 1) && push!(edges, (idx[i, j], idx[i - 1, j]))
    (j > 1) && push!(edges, (idx[i, j], idx[i, j - 1]))
    (i < nr) && push!(edges, (idx[i, j], idx[i + 1, j]))
    (j < nc) && push!(edges, (idx[i, j], idx[i, j + 1]))
end
bmod = BayesianBinomialGFEN(edges, tv1=λ1, tv2=λ2)

##
n = 3000
thinning = 5
burnin = 2000

##
init = model.beta  # zeros(size(mcmc_init))
@time chain = sample_chain(bmod, succ, att, n, burnin=burnin, thinning=thinning, init=init, verbose=true, async=true)

##
z = logistic.(chain)
q75 = mapslices(x -> quantile(x, 0.75), z, dims=2)
q50 = mapslices(x -> quantile(x, 0.50), z, dims=2)
q25 = mapslices(x -> quantile(x, 0.25), z, dims=2)
iqr = q75 .- q25

unc = reshape(iqr, (nr, nc))
unc = (unc .- minimum(unc)) ./ (maximum(unc) .- minimum(unc))
npzwrite("simulations/geometries/bayes_iqr.npy", Float32.(unc))
unc = Float64.(npzread("simulations/geometries/bayes_iqr.npy"))
##

plot(Gray.(1 .- unc))
xlabel!("time")
ylabel!("space")

## posterior covariance  How to order points?

## posterior inverse covariance

##
savefig("simulations/geometries/geometries_iqr.pdf")
##

median_signal = reshape(q50, (nr, nc))
npzwrite("simulations/geometries/bayes_median.npy", Float32.(median_signal))
median_signal = Float64.(npzread("simulations/geometries/bayes_median.npy"))
plot(Gray.(median_signal))
xlabel!("time")
ylabel!("space")

## now with shrinking

##
ϵ = 0.1
succ = vec(successes) .+ ϵ
att = vec(attempts) .+ 2ϵ
rmodel = BinomialGFEN(ptr, brks, lambdasl1, zeros(size(lambdasl2)); modelopts...)
fit!(rmodel, succ, att; fitopts...)
ω = logistic.(rmodel.beta)
ω = reshape(ω, nr, nc)

plot(Gray.(1 .- ω))
xlabel!("time")
ylabel!("space")

##

max_stars_repo_name: mauriciogtec/gfen-reproduce

max_issues_repo_name: mauriciogtec/gfen-reproduce

max_forks_repo_name: mauriciogtec/gfen-reproduce


text:
function helloworld()
    println("Hello world from hello.jl")
end

max_stars_repo_name: PacktPublishing/Learning-Julia

max_issues_repo_name: PacktPublishing/Learning-Julia

max_forks_repo_name: PacktPublishing/Learning-Julia


text:
"""
    test_bqg_rossbywave(stepper, dt, nsteps, dev; kwargs...)

Evolves a Rossby wave and compares with the analytic solution.
"""
function test_bqg_rossbywave(stepper, dt, nsteps, dev::Device=CPU())
    nx = 64
    Lx = 2π
     β = 2.0
     μ = 0.0
     ν = 0.0
     T = Float64

  # the following if statement is called so that all the cases of
  # Problem() fuction are tested
  if stepper=="ForwardEuler"
    eta = zeros(dev, T, (nx, nx))
  else
    eta(x, y) = 0*x
  end

  prob = BarotropicQG.Problem(dev; nx=nx, Lx=Lx, eta=eta, β=β, μ=μ, ν=ν, stepper=stepper, dt=dt)
  sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid

  x, y = gridpoints(g)

  # the Rossby wave initial condition
   ampl = 1e-2
  kwave = 3 * 2π/g.Lx
  lwave = 2 * 2π/g.Ly
      ω = -p.β*kwave/(kwave^2 + lwave^2)
     ζ0 = @. ampl*cos(kwave*x)*cos(lwave*y)
    ζ0h = rfft(ζ0)

  BarotropicQG.set_zeta!(prob, ζ0)

  stepforward!(prob, nsteps)
  dealias!(sol, g)
  BarotropicQG.updatevars!(prob)

  ζ_theory = @. ampl*cos(kwave*(x - ω/kwave*cl.t)) * cos(lwave*y)

  return isapprox(ζ_theory, v.zeta, rtol=g.nx*g.ny*nsteps*1e-12)
end

"""
    test_bqg_stochasticforcingbudgets(dev; kwargs...)

Tests if the energy budget is closed for BarotropicQG problem with stochastic forcing.
"""
function test_bqg_stochasticforcingbudgets(dev::Device=CPU(); n=256, dt=0.01, L=2π, ν=1e-7, nν=2, μ=1e-1, T=Float64)
  n, L  = 256, 2π
  ν, nν = 1e-7, 2
  μ = 1e-1
  dt, tf = 0.005, 0.1/μ
  nt = round(Int, tf/dt)

    gr = TwoDGrid(dev, n, L)
  x, y = gridpoints(gr)

  # Forcing
  kf, dkf = 12.0, 2.0
  ε = 0.1

  CUDA.@allowscalar Kr = ArrayType(dev)([gr.kr[i] for i=1:gr.nkr, j=1:gr.nl])

  forcingcovariancespectrum = zeros(dev, T, (gr.nkr, gr.nl))
  @. forcingcovariancespectrum = exp( -(sqrt(gr.Krsq) - kf)^2 / (2 * dkf^2) )
  CUDA.@allowscalar @. forcingcovariancespectrum[gr.Krsq .< 2^2] = 0
  CUDA.@allowscalar @. forcingcovariancespectrum[gr.Krsq .> 20^2] = 0
  CUDA.@allowscalar @. forcingcovariancespectrum[Kr .< 2π/L] = 0
  ε0 = parsevalsum(forcingcovariancespectrum .* gr.invKrsq / 2, gr) / (gr.Lx * gr.Ly)
  forcingcovariancespectrum .= ε / ε0 * forcingcovariancespectrum

  Random.seed!(1234)

  function calcFq!(Fqh, sol, t, clock, vars, params, grid)
    eta = ArrayType(dev)(exp.(2π * im * rand(T, size(sol))) / sqrt(clock.dt))
    CUDA.@allowscalar eta[1, 1] = 0
    @. Fqh = eta * sqrt(forcingcovariancespectrum)
    return nothing
  end

  prob = BarotropicQG.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, dt=dt,
   stepper="RK4", calcFq=calcFq!, stochastic=true)

  BarotropicQG.set_zeta!(prob, 0*x)

  E = Diagnostic(BarotropicQG.energy,             prob, nsteps=nt)
  D = Diagnostic(BarotropicQG.energy_dissipation, prob, nsteps=nt)
  R = Diagnostic(BarotropicQG.energy_drag,        prob, nsteps=nt)
  W = Diagnostic(BarotropicQG.energy_work,        prob, nsteps=nt)
  diags = [E, D, W, R]

  stepforward!(prob, diags, nt)

  BarotropicQG.updatevars!(prob)

  dEdt_numerical = (E[2:E.i] - E[1:E.i-1]) / prob.clock.dt

  # If the Ito interpretation was used for the work
  # then we need to add the drift term
  # dEdt_computed = W[2:E.i] + ε - D[1:E.i-1] - R[1:E.i-1]      # Ito
  dEdt_computed = W[2:E.i] - D[1:E.i-1] - R[1:E.i-1]        # Stratonovich

  return isapprox(dEdt_numerical, dEdt_computed, rtol=1e-3)
end

"""
    test_bqg_deterministicforcingbudgets(dev; kwargs...)

Tests if the energy budget is closed for BarotropicQG problem with deterministic forcing.
"""
function test_bqg_deterministicforcingbudgets(dev::Device=CPU(); n=256, dt=0.01, L=2π, ν=1e-7, nν=2, μ=1e-1)
  n, L  = 256, 2π
  ν, nν = 1e-7, 2
  μ = 1e-1
  dt, tf = 0.005, 0.1/μ
  nt = round(Int, tf/dt)

  gr = TwoDGrid(dev, n, L)
  x, y = gridpoints(gr)
  k₀, l₀ = 2π/gr.Lx, 2π/gr.Ly

  # Forcing = 0.01cos(4x)cos(5y)cos(2t)
  f = @. 0.01 * cos(4k₀*x) * cos(5l₀*y)
  fh = rfft(f)

  function calcFq!(Fqh, sol, t, clock, vars, params, grid)
    @. Fqh = fh * cos(2t)
    return nothing
  end

  prob = BarotropicQG.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, dt=dt,
   stepper="RK4", calcFq=calcFq!, stochastic=false)

  BarotropicQG.set_zeta!(prob, 0*x)

  E = Diagnostic(BarotropicQG.energy,             prob, nsteps=nt)
  D = Diagnostic(BarotropicQG.energy_dissipation, prob, nsteps=nt)
  R = Diagnostic(BarotropicQG.energy_drag,        prob, nsteps=nt)
  W = Diagnostic(BarotropicQG.energy_work,        prob, nsteps=nt)
  diags = [E, D, W, R]

  stepforward!(prob, diags, round(Int, nt))

  BarotropicQG.updatevars!(prob)

  dEdt_numerical = (E[3:E.i] - E[1:E.i-2]) / (2 * prob.clock.dt)
  dEdt_computed  = W[2:E.i-1] - D[2:E.i-1] - R[2:E.i-1]

  residual = dEdt_numerical - dEdt_computed

  return isapprox(dEdt_numerical, dEdt_computed, atol=1e-10)
end

"""
    test_bqg_nonlinearadvection(dt, stepper, dev; kwargs...)

Tests the advection term in the BarotropicQG module by timestepping a
test problem with timestep dt and timestepper identified by the string stepper.
The test problem is derived by picking a solution ζf (with associated
streamfunction ψf) for which the advection term J(ψf, ζf) is non-zero. Next, a
forcing Ff is derived according to Ff = ∂ζf/∂t + J(ψf, ζf) - ν∇²ζf. One solution
to the vorticity equation forced by this Ff is then ζf. (This solution may not
be realized, at least at long times, if it is unstable.)
"""
function test_bqg_advection(dt, stepper, dev::Device=CPU(); n=128, L=2π, ν=1e-2, nν=1, μ=0.0)
  n, L  = 128, 2π
  ν, nν = 1e-2, 1
   μ = 0.0
  tf = 1.0
  nt = round(Int, tf/dt)

  gr = TwoDGrid(dev, n, L)
  x, y = gridpoints(gr)

  ψf = @. sin(2x) * cos(2y) + 2sin(x) * cos(3y)
  qf = @. -8sin(2x) * cos(2y) - 20sin(x) * cos(3y)

  Ff = @. -(
    ν*( 64sin(2x) * cos(2y) + 200sin(x) * cos(3y) )
    + 8 * ( cos(x) * cos(3y) * sin(2x) * sin(2y) - 3cos(2x) * cos(2y) * sin(x) * sin(3y) )
    )

  Ffh = rfft(Ff)

  # Forcing
  function calcFq!(Fqh, sol, t, cl, v, p, g)
    Fqh .= Ffh
    return nothing
  end

  prob = BarotropicQG.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, dt=dt, stepper=stepper, calcFq=calcFq!)

  BarotropicQG.set_zeta!(prob, qf)

  stepforward!(prob, round(Int, nt))

  BarotropicQG.updatevars!(prob)

  return isapprox(prob.vars.q, qf, rtol=rtol_barotropicQG)
end

"""
    test_bqg_formstress(dt, stepper, dev; kwargs...)

Tests the form stress term that forces the domain-averaged zonal flow U(t).
"""
function test_bqg_formstress(dt, stepper, dev::Device=CPU(); n=128, L=2π, ν=0.0, nν=1, μ=0.0)
  n, L  = 128, 2π
  ν, nν = 1e-2, 1
  μ = 0.0
  tf = 1
  nt = 1

  gr  = TwoDGrid(dev, n, L)
  x, y = gridpoints(gr)

  ζ₀ = @. -20 * sin(10x) * cos(10y)
  topoPV(x, y) = @. cos(10x) * cos(10y)
  F(t) = 0 #no forcing

  answer = 0.25 # this is what <v*eta> should be

  prob = BarotropicQG.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, dt=dt, stepper=stepper, eta=topoPV, calcFU=F)
  BarotropicQG.set_zeta!(prob, ζ₀)
  BarotropicQG.updatevars!(prob)

  # Step forward
  stepforward!(prob, nt)
  return CUDA.@allowscalar isapprox(prob.timestepper.N[1, 1], answer, rtol=rtol_barotropicQG)
end

"""
    test_bqg_energyenstrophy(dev)

Tests the energy and enstrophy function for a BarotropicQG problem.
"""
function test_bqg_energyenstrophy(dev::Device=CPU())
  nx, Lx  = 64, 2π
  ny, Ly  = 64, 3π
  g  = TwoDGrid(dev, nx, Lx, ny, Ly)
  k₀, l₀ = 2π/Lx, 2π/Ly # fundamental wavenumbers
  x, y = gridpoints(g)

  energy_calc = 29/9
  enstrophy_calc = 2701/162

    eta = @. cos(10k₀*x) * cos(10l₀*y)
   psi0 = @. sin(2k₀*x) * cos(2l₀*y) + 2sin(k₀*x) * cos(3l₀*y)
  zeta0 = @. - ((2k₀)^2+(2l₀)^2) * sin(2k₀*x) * cos(2l₀*y) - (k₀^2+(3l₀)^2) * 2sin(k₀*x) * cos(3l₀*y)

  prob = BarotropicQG.Problem(dev; nx=nx, Lx=Lx, ny=ny, Ly=Ly, eta = eta, stepper="ForwardEuler")
  sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
  BarotropicQG.set_zeta!(prob, zeta0)
  BarotropicQG.updatevars!(prob)

  energyzeta0 = BarotropicQG.energy(prob)
  enstrophyzeta0 = BarotropicQG.enstrophy(prob)

  return isapprox(energyzeta0, energy_calc, rtol=rtol_barotropicQG) && isapprox(enstrophyzeta0, enstrophy_calc, rtol=rtol_barotropicQG) &&
  BarotropicQG.addforcing!(prob.timestepper.N, sol, cl.t, cl, v, p, g)==nothing
end

"""
    test_bqg_energyenstrophy(dev)

Tests the mean flow U(t) energy and enstrophy function for a BarotropicQG problem.
"""
function test_bqg_meanenergyenstrophy(dev::Device=CPU())
  nx, Lx  = 64, 2π
  ny, Ly  = 96, 3π
  g = TwoDGrid(dev, nx, Lx, ny, Ly)
  k₀, l₀ = 2π/Lx, 2π/Ly # fundamental wavenumbers
  x, y = gridpoints(g)

  calcFU(t) = 0.0
  eta(x, y) = @. cos(10x) * cos(10y)
  psi0 = @. sin(2k₀*x) * cos(2l₀*y) + 2sin(k₀*x) * cos(3l₀*y)
 zeta0 = @. - ((2k₀)^2+(2l₀)^2) * sin(2k₀*x) * cos(2l₀*y) - (k₀^2+(3l₀)^2) * 2sin(k₀*x) * cos(3l₀*y)
  β = 10.0
  U = 1.2

  energy_calc = 29/9
  enstrophy_calc = 2701/162

  prob = BarotropicQG.Problem(dev; nx=nx, Lx=Lx, ny=ny, Ly=Ly, β=β, eta=eta, calcFU = calcFU,
                                    stepper="ForwardEuler")

  BarotropicQG.set_zeta!(prob, zeta0)
  BarotropicQG.set_U!(prob, U)
  BarotropicQG.updatevars!(prob)

  energyU = BarotropicQG.meanenergy(prob)
  enstrophyU = BarotropicQG.meanenstrophy(prob)

  energyzeta0 = BarotropicQG.energy(prob)
  enstrophyzeta0 = BarotropicQG.enstrophy(prob)

  return (isapprox(energyU, 0.5*U^2, rtol=rtol_barotropicQG) &&
    isapprox(enstrophyU, β*U, rtol=rtol_barotropicQG) &&
    isapprox(energyzeta0, energy_calc, rtol=rtol_barotropicQG) &&
    isapprox(enstrophyzeta0, enstrophy_calc, rtol=rtol_barotropicQG)
  )
end

"""
    test_bqg_problemtype(dev, T)

Tests the BarotropicQG problem constructor for different DataType `T`.
"""
function test_bqg_problemtype(dev, T)
  prob = BarotropicQG.Problem(dev; T=T)

  A = ArrayType(dev)

  return (typeof(prob.sol)<:A{Complex{T}, 2} && typeof(prob.grid.Lx)==T && eltype(prob.grid.x)==T && typeof(prob.vars.u)<:A{T, 2})
end


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 2a600f08-6660-11eb-3d8a-cdb566c812d5
using XKCD, PlutoUI

# ╔═╡ 38f28954-6660-11eb-298e-cd02923b4fb0
current = Comic();

# ╔═╡ 422e9fe6-6660-11eb-0687-f9fbb67772f5
@bind sel NumberField(1:current.json.num, default=current.json.num)

# ╔═╡ 5579fffa-6660-11eb-2104-a5583bf197b2
Comic(sel)

# ╔═╡ Cell order:
# ╠═2a600f08-6660-11eb-3d8a-cdb566c812d5
# ╠═38f28954-6660-11eb-298e-cd02923b4fb0
# ╠═422e9fe6-6660-11eb-0687-f9fbb67772f5
# ╠═5579fffa-6660-11eb-2104-a5583bf197b2


text:
import Base: depwarn

export Info
# Compatibility with `info`

const IntNum = Union{Integer, Tuple{Integer,Integer}, Tuple{Integer,Integer,Integer}}

mutable struct Info
    n::IntNum
    nsrc::Integer
    nt::Array{Integer,1}
end

function Info(n::IntNum, nsrc::Integer, nt::Integer)
    depwarn("Info is deprecated and will be removed in future versions", :Info; force=true)
    Info(n, nsrc, [nt for i=1:nsrc])
end

for f in [:judiModeling, :judiProjection, :judiWavefield]
    @eval function $f(info::Info, ar...;kw...)
        depwarn("$($f)(info::Info, ar...; kw...) is deprecated, use $($f)(ar...; kw...)", Symbol($f); force=true)
        $f(ar...; kw...)
    end
end

for f in [:judiLRWF, :judiRHS]
    @eval $f(info::Info, ar...;kw...) = throw(ArgumentError("$($f)(info::Info, ...) is deprecated and requires a time sampling rate `dt`"))
end


text:
using AutoHashEquals

"""
    Model{S,D}

Structure for representing the petri net model

represented by states and transition functions
"""
@auto_hash_equals struct Model{S,D}
  S::S  # states
  Δ::D  # transition function
end

Model(s::S, Δ) where S<:UnitRange = Model(collect(s), Δ)
Model(s::S, Δ) where S<:Int = Model(1:s, Δ)

"""
    NullPetri

create a Petri net of no states and no transitions
"""
const NullPetri = Model(Int[], Vector{Tuple{Dict{Int, Number},Dict{Int, Number}}}())

"""
    EmptyPetri(n::Int)

create a Petri net of ``n`` states with no transitions
"""
EmptyPetri(n::Int) = Model(collect(1:n), Vector{Tuple{Dict{Int, Number},Dict{Int, Number}}}())

"""
    EmptyPetri(s::Vector{T})

create a Petri net with states ``s`` with no transitions
"""
EmptyPetri(s::Vector{T}) where T = Model(s, Vector{Tuple{Dict{T, Number},Dict{T, Number}}}())


text:
using OrdinaryDiffEq, DiffEqDevTools, Sundials, ModelingToolkit, ODEInterfaceDiffEq,
      Plots, DASSL, DASKR
using LinearAlgebra

@variables t y₁(t)=1.0 y₂(t)=0.0 y₃(t)=0.0
@parameters k₁=0.04 k₂=3e7 k₃=1e4
D = Differential(t)
eqs = [
  D(y₁) ~ -k₁*y₁ + k₃*y₂*y₃
  D(y₂) ~  k₁*y₁ - k₃*y₂*y₃ - k₂*y₂^2
  0 ~  y₁ + y₂ + y₃ - 1
]
@named sys = ODESystem(eqs)
simpsys = structural_simplify(sys)
mmprob = ODEProblem(sys,[],(0.0,1e5))
daeprob = DAEProblem(sys,[D(y₁)=>-0.04,
                              D(y₂)=>0.04,
                              D(y₃)=>0.0],[],(0.0,1e5))
odaeprob = ODAEProblem(simpsys,[],(0.0,1e5))

ref_sol = solve(daeprob,IDA(),abstol=1/10^14,reltol=1/10^14);
ode_ref_sol = solve(odaeprob,CVODE_BDF(),abstol=1/10^14,reltol=1/10^14);

probs = [mmprob,daeprob,odaeprob]
refs = [ref_sol,ref_sol,ode_ref_sol];

plot(ode_ref_sol, vars = [y₁,y₂,y₃])

abstols = 1.0 ./ 10.0 .^ (5:8)
reltols = 1.0 ./ 10.0 .^ (1:4);
setups = [Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 1, :alg=>Rodas4()),
          #Dict(:prob_choice => 1, :alg=>FBDF()),
          Dict(:prob_choice => 1, :alg=>QNDF()),
          Dict(:prob_choice => 1, :alg=>rodas()),
          Dict(:prob_choice => 1, :alg=>radau()),
          Dict(:prob_choice => 1, :alg=>RadauIIA5()),
          Dict(:prob_choice => 2, :alg=>DFBDF()),
          Dict(:prob_choice => 2, :alg=>IDA()),
]

wp = WorkPrecisionSet(probs,abstols,reltols,setups;
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (6:8)
reltols = 1.0 ./ 10.0 .^ (2:4);
setups = [Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 1, :alg=>Rodas4()),
          Dict(:prob_choice => 2, :alg=>IDA()),
          Dict(:prob_choice => 3, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 3, :alg=>Rodas4()),
          Dict(:prob_choice => 3, :alg=>CVODE_BDF()),
          Dict(:prob_choice => 3, :alg=>TRBDF2()),
          Dict(:prob_choice => 3, :alg=>KenCarp4()),
          ]
wp = WorkPrecisionSet(probs,abstols,reltols,setups;
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (6:8)
reltols = 1.0 ./ 10.0 .^ (3:5);
setups = [Dict(:prob_choice => 3, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 3, :alg=>Rodas4()),
          Dict(:prob_choice => 2, :alg=>IDA()),
          Dict(:prob_choice => 2, :alg=>DASSL.dassl()),
          Dict(:prob_choice => 2, :alg=>DASKR.daskr()),
          ]
wp = WorkPrecisionSet(probs,abstols,reltols,setups;
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (5:8)
reltols = 1.0 ./ 10.0 .^ (1:4);
setups = [Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 1, :alg=>Rodas4()),
          #Dict(:prob_choice => 1, :alg=>FBDF()), # too slow
          Dict(:prob_choice => 1, :alg=>QNDF()),
          Dict(:prob_choice => 1, :alg=>rodas()),
          Dict(:prob_choice => 1, :alg=>radau()),
          Dict(:prob_choice => 1, :alg=>RadauIIA5()),
          Dict(:prob_choice => 2, :alg=>DFBDF()),
          Dict(:prob_choice => 2, :alg=>IDA()),
          ]
wp = WorkPrecisionSet(probs,abstols,reltols,setups;error_estimate = :l2,
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (6:8)
reltols = 1.0 ./ 10.0 .^ (2:4);
setups = [Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 1, :alg=>Rodas4()),
          Dict(:prob_choice => 2, :alg=>IDA()),
          Dict(:prob_choice => 3, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 3, :alg=>Rodas4()),
          Dict(:prob_choice => 3, :alg=>CVODE_BDF()),
          Dict(:prob_choice => 3, :alg=>TRBDF2()),
          Dict(:prob_choice => 3, :alg=>KenCarp4()),
          ]
wp = WorkPrecisionSet(probs,abstols,reltols,setups;error_estimate = :l2,
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (7:12)
reltols = 1.0 ./ 10.0 .^ (4:9)

setups = [Dict(:prob_choice => 1, :alg=>Rodas5()),
          Dict(:prob_choice => 3, :alg=>Rodas5()),
          Dict(:prob_choice => 1, :alg=>Rodas4()),
          Dict(:prob_choice => 3, :alg=>Rodas4()),
          #Dict(:prob_choice => 1, :alg=>FBDF()),
          Dict(:prob_choice => 1, :alg=>QNDF()),
          Dict(:prob_choice => 1, :alg=>rodas()),
          Dict(:prob_choice => 1, :alg=>radau()),
          Dict(:prob_choice => 1, :alg=>RadauIIA5()),
          Dict(:prob_choice => 2, :alg=>DFBDF()),
          Dict(:prob_choice => 2, :alg=>IDA()),
          ]

wp = WorkPrecisionSet(probs,abstols,reltols,setups;
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

wp = WorkPrecisionSet(probs,abstols,reltols,setups;error_estimate = :l2,
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (7:9)
reltols = 1.0 ./ 10.0 .^ (4:6)

setups = [Dict(:prob_choice => 1, :alg=>Rodas5()),
          Dict(:prob_choice => 3, :alg=>Rodas5()),
          Dict(:prob_choice => 1, :alg=>Rodas4()),
          Dict(:prob_choice => 3, :alg=>Rodas4()),
          Dict(:prob_choice => 2, :alg=>DFBDF()),
          Dict(:prob_choice => 2, :alg=>IDA()),
          Dict(:prob_choice => 2, :alg=>DASKR.daskr()),
          ]
gr()
wp = WorkPrecisionSet(probs,abstols,reltols,setups;
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

using SciMLBenchmarks
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

max_stars_repo_name: SciML/SciMLBenchmarksOutput

max_issues_repo_name: SciML/SciMLBenchmarksOutput

max_forks_repo_name: SciML/SciMLBenchmarksOutput


text:
module XClipboard

import Base.Multimedia.@textmime
export xclipboard, xclipboard_targets, decodemime

include("XLib.jl")
include("xmime.jl")

# Get a list of clipboard targets, as X11 Atoms
xclipboard_target_atoms(timeout::Int=2) =
    reinterpret(Atom, xclipboard_request(XA_TARGETS, timeout))

# Get a list of targets as MIME types
xclipboard_targets(targets::Array{Atom}) =
    [MIME(XGetAtomName(atom)) for atom in targets]
xclipboard_targets(timeout::Int=2) =
    xclipboard_targets(xclipboard_target_atoms(timeout))

# Get a Dict of clipboard targets as (MIME_type => Atom)[]
function xclipboard_targets_atoms(timeout::Int=2)
    atoms = xclipboard_target_atoms(timeout)
    return Dict(xclipboard_targets(atoms), atoms)
end

# Main xclipboard function
function xclipboard(mimetypes::Array{MIME}, timeout::Int=2)

    atoms = xclipboard_targets_atoms(timeout)
    local request_atom, mimetype
    for mimetype in mimetypes
        if (request_atom = get(atoms, mimetype, NoneAtom)) != NoneAtom
            break
        end
    end

    if request_atom == NoneAtom
        error("Requested clipboard target(s) not found.")
    end

    data = xclipboard_request(request_atom, timeout)

    return decodemime(mimetype, data)
end
xclipboard(mimetype::MIME, timeout::Int=2) = xclipboard([mimetype], timeout)
xclipboard(mimetype::String...) = xclipboard([MIME(m) for m in mimetype])
xclipboard() = xclipboard("UTF8_STRING", "STRING", "TEXT")

end # module


text:
module FieldValuesTests

using Test
using Gridap
using TensorValues
using LinearAlgebra: dot
using LinearAlgebra: norm

a = 1
@test isa(a,FieldValue)

a = 1.0
@test isa(a,FieldValue)

a = 1.0 - 3.0im
@test isa(a,FieldValue)

a = VectorValue(1,2,3)
@test isa(a,FieldValue)

p = Point(1,2,3)

@test isa(p,VectorValue{3,Int})

p = Point{3,Float64}(1,2,3)

@test isa(p,VectorValue{3,Float64})

v = MultiValue{Tuple{1,2}}(10,20)
n = normalvec(v)
@test n == VectorValue(20,-10)
@test meas(v) == sqrt(500)

v = MultiValue{Tuple{2,3}}(1,0,0,1,0,0)
n = normalvec(v)
@test n == VectorValue(0,0,1)
@test meas(v) ≈ 1.0

v = MultiValue{Tuple{2,3}}(1,0,0,1,1,0)
n = normalvec(v)
@test n == VectorValue(-1,0,1)
@test meas(v) ≈ sqrt(2)

t = TensorValue(1,2,3,4)
@test trace(t) == 5
@test tr(t) == 5

t = TensorValue(1,2,3,4,5,6,7,8,9)
@test trace(t) == 15
@test tr(t) == 15

@test symmetic_part(t) == TensorValue(1.0, 3.0, 5.0, 3.0, 5.0, 7.0, 5.0, 7.0, 9.0)

a = TensorValue(1,2,3,4)
b = a'
@test b == TensorValue(1,3,2,4)
@test a*b == TensorValue(10,14,14,20)

u = VectorValue(1.0,2.0)
v = VectorValue(2.0,3.0)
@test dot(u,v) ≈ inner(u,v)
@test norm(u) ≈ sqrt(inner(u,u))

end # module


text:
module HDSolver

# ----------
# Navier–Stokes Solver for 3D Hydrodynamics Problem
# ----------

export
	UᵢUpdate!,
	HDcalcN_advection!,
	HDupdatevars!

using
  CUDA,
  Reexport,
  DocStringExtensions

@reexport using FourierFlows

using LinearAlgebra: mul!, ldiv!
using FourierFlows: parsevalsum

# δ function
δ(a::Int,b::Int) = ( a == b ? 1 : 0 );

function UᵢUpdate!(N, sol, t, clock, vars, params, grid;direction="x")

  if direction == "x"

  	# a = {1,2,3} -> {x,y,z} direction
  	a    = 1;
  	kₐ   = grid.kr;
  	k⁻²  = grid.invKrsq;
  	uᵢh  = vars.uxh;
  	∂uᵢh∂t = @view N[:,:,:,params.ux_ind];

  elseif direction == "y"

  	a    = 2;
  	kₐ   = grid.l;
  	k⁻²  = grid.invKrsq;
  	uᵢh  = vars.uyh;
  	∂uᵢh∂t = @view N[:,:,:,params.uy_ind];

  elseif direction == "z"

  	a    = 3;
  	kₐ   = grid.m;
  	k⁻²  = grid.invKrsq;
  	uᵢh  = vars.uzh;
  	∂uᵢh∂t = @view N[:,:,:,params.uz_ind];

  else

  	@warn "Warning : Unknown direction is declerad"

  end

  @. ∂uᵢh∂t*= 0;

  for (uᵢ,kᵢ) ∈ zip([vars.ux,vars.uy,vars.uz],[grid.kr,grid.l,grid.m])
        for (uⱼ,kⱼ,j) ∈ zip([vars.ux,vars.uy,vars.uz],[grid.kr,grid.l,grid.m],[1, 2, 3])

          # Initialization
          @. vars.nonlin1 *= 0;
          uᵢuⱼ  = vars.nonlin1;
          uᵢuⱼh = vars.nonlinh1;

          # Pre-Calculation in Real Space
          @. uᵢuⱼ = uᵢ*uⱼ;

          # Fourier transform
          mul!(uᵢuⱼh, grid.rfftplan, uᵢuⱼ);

          # Perform the actual calculation
          @. ∂uᵢh∂t += -im*kᵢ*(δ(a,j)-kₐ*kⱼ*k⁻²)*uᵢuⱼh;

        end
    end

    #Compute the diffusion term  - νk^2 u_i
    @. ∂uᵢh∂t += -grid.Krsq*params.ν*uᵢh;
    return nothing

end

function HDcalcN_advection!(N, sol, t, clock, vars, params, grid)

  #Update V  Fourier Conponment
  @. vars.uxh = sol[:, :, :, params.ux_ind];
  @. vars.uyh = sol[:, :, :, params.uy_ind];
  @. vars.uzh = sol[:, :, :, params.uz_ind];

  #Update V + B Real Conponment
  ldiv!(vars.ux, grid.rfftplan, deepcopy(vars.uxh))
  ldiv!(vars.uy, grid.rfftplan, deepcopy(vars.uyh))
  ldiv!(vars.uz, grid.rfftplan, deepcopy(vars.uzh))

  #Update V Advection
  UᵢUpdate!(N, sol, t, clock, vars, params, grid;direction="x")
  UᵢUpdate!(N, sol, t, clock, vars, params, grid;direction="y")
  UᵢUpdate!(N, sol, t, clock, vars, params, grid;direction="z")

  return nothing
end

function HDupdatevars!(prob)
  vars, grid, sol, params = prob.vars, prob.grid, prob.sol, prob.params

  dealias!(sol, grid)

  #Update V Fourier Conponment
  @. vars.uxh = sol[:, :, :, params.ux_ind];
  @. vars.uyh = sol[:, :, :, params.uy_ind];
  @. vars.uzh = sol[:, :, :, params.uz_ind];

  #Update V + B Real Conponment
  ldiv!(vars.ux, grid.rfftplan, deepcopy(vars.uxh)) # deepcopy() since inverse real-fft destroys its input
  ldiv!(vars.uy, grid.rfftplan, deepcopy(vars.uyh)) # deepcopy() since inverse real-fft destroys its input
  ldiv!(vars.uz, grid.rfftplan, deepcopy(vars.uzh)) # deepcopy() since inverse real-fft destroys its input

  return nothing
end

end

max_stars_repo_name: MHDFlows/MHDFlows

max_issues_repo_name: MHDFlows/MHDFlows

max_forks_repo_name: MHDFlows/MHDFlows


text:
using JSON, Plots, StatsPlots

#j_data = open("COVID19_Finland.json", "r")

j_data = read("COVID19-Finland.json", String)

jobj = JSON.Parser.parse(j_data)

println(jobj[1])

cases = map(x -> get(x, "Cases", 0), jobj)
plot(cases, xlabel = "Number of Cases")


text:
const HistoryMaximum = 400
const HistoryMultiply = 32
const HistoryDivide = 512

"""
    updatehistory!(thread::Thread, quiets::MoveStack, ply::Int, depthbonus::Int)

A function to update the history heuristics stored in `thread`, by passing a list of attempted quiet moves, the current `ply`, and a depth bonus.
"""
function updatehistory!(thread::Thread, quietstried::MoveStack, ply::Int, depthbonus::Int)::Nothing
    colour = thread.board.turn.val
    @inbounds best_move = quietstried[quietstried.idx]

    # Extract counter move information.
    # If we haven't made a ply, we can't extract this.
    if ply > 0
        @inbounds counter = thread.movestack[ply]
        @inbounds cm_piece = thread.piecestack[ply]
        cm_to = to(counter)
    else
        counter = MOVE_NONE
        cm_piece = VOID
        cm_to = zero(UInt16)
    end

    # Extract move from 2 moves ago.
    # Must have made more than 1 ply to be able to do this.
    if ply > 1
        @inbounds follow = thread.movestack[ply - 1]
        @inbounds fm_piece = thread.piecestack[ply - 1]
        fm_to = to(follow)
    else
        follow = MOVE_NONE
        fm_piece = VOID
        fm_to = zero(UInt16)
    end

    # Set the killer moves.
    killer1 = thread.killer1s[ply + 1]
    if killer1 !== best_move
        @inbounds thread.killer2s[ply + 1] = killer1
        @inbounds thread.killer1s[ply + 1] = best_move
    end

    # Set the counter move.
    if (counter !== MOVE_NONE) && (counter !== NULL_MOVE)
        thread.cmtable[(!thread.board.turn).val][cm_piece.val][cm_to] = best_move
    end

    if quietstried.idx == 1 && depthbonus <= 1
        return
    end

    bonus = min(depthbonus, HistoryMaximum)

    # Below, the decision to split the cases up into 4 functions is so that the compiler can allow for SIMD optimisations.
    if (counter !== MOVE_NONE) && (counter !== NULL_MOVE) && (follow !== MOVE_NONE) && (follow !== NULL_MOVE)
        updatehistory_internal_countfollow!(thread, quietstried, bonus, best_move, colour, counter, cm_piece, cm_to, follow, fm_piece, fm_to)
    elseif (counter !== MOVE_NONE) && (counter !== NULL_MOVE)
        updatehistory_internal_count!(thread, quietstried, bonus, best_move, colour, counter, cm_piece, cm_to)
    elseif (follow !== MOVE_NONE) && (follow !== NULL_MOVE)
        updatehistory_internal_follow!(thread, quietstried, bonus, best_move, colour, follow, fm_piece, fm_to)
    else
        updatehistory_internal!(thread, quietstried, bonus, best_move, colour)
    end

    return
end

# Update the history heuristics for the case where we a counter move.
function updatehistory_internal_count!(thread::Thread, quietstried::MoveStack, bonus::Int, best_move::Move, colour::UInt8,
    counter::Move, cm_piece::PieceType, cm_to::Integer)::Nothing
    abs_δ = abs(bonus)
    @inbounds for move in quietstried
        δ = (move == best_move) ? bonus : -bonus
        sqr_to = to(move)
        sqr_from = from(move)
        move_piece = type(thread.board[sqr_from])

        mul_common = HistoryMultiply * δ

        # Update the butterfly table. (History heuristics)
        entry = thread.history[colour][sqr_from]
        entry[sqr_to] += mul_common - fld(entry[sqr_to] * abs_δ, HistoryDivide)

        # Update the counter move history table.
        entry_c = thread.counterhistory[cm_piece.val][cm_to][move_piece.val]
        entry_c[sqr_to] += mul_common - fld(entry_c[sqr_to] * abs_δ, HistoryDivide)
    end
    return
end

# Function to update history heuristics when there is a follow up move.
function updatehistory_internal_follow!(thread::Thread, quietstried::MoveStack, bonus::Int, best_move::Move, colour::UInt8,
    follow::Move, fm_piece::PieceType, fm_to::Integer)::Nothing
    abs_δ = abs(bonus)
    @inbounds for move in quietstried
        δ = (move == best_move) ? bonus : -bonus
        sqr_to = to(move)
        sqr_from = from(move)
        move_piece = type(thread.board[sqr_from])

        mul_common = HistoryMultiply * δ

        # Update the butterfly table. (History heuristics)
        entry = thread.history[colour][sqr_from]
        entry[sqr_to] += mul_common - fld(entry[sqr_to] * abs_δ, HistoryDivide)

        # Update the follow up move history table.
        entry_f = thread.followhistory[fm_piece.val][fm_to][move_piece.val]
        entry_f[sqr_to] += mul_common - fld(entry_f[sqr_to] * abs_δ, HistoryDivide)
    end
    return
end

# Internals to update the history heuristics when we have both a counter and follow up move.
function updatehistory_internal_countfollow!(thread::Thread, quietstried::MoveStack, bonus::Int, best_move::Move, colour::UInt8,
    counter::Move, cm_piece::PieceType, cm_to::Integer, follow::Move, fm_piece::PieceType, fm_to::Integer)::Nothing
    abs_δ = abs(bonus)
    @inbounds for move in quietstried
        δ = (move == best_move) ? bonus : -bonus
        sqr_to = to(move)
        sqr_from = from(move)
        move_piece = type(thread.board[sqr_from])

        mul_common = HistoryMultiply * δ

        # Update the butterfly table. (History heuristics)
        entry = thread.history[colour][sqr_from]
        entry[sqr_to] += mul_common - fld(entry[sqr_to] * abs_δ, HistoryDivide)

        # Update the counter move history table.
        entry_c = thread.counterhistory[cm_piece.val][cm_to][move_piece.val]
        entry_c[sqr_to] += mul_common - fld(entry_c[sqr_to] * abs_δ, HistoryDivide)

        # Update the follow up move history table.
        entry_f = thread.followhistory[fm_piece.val][fm_to][move_piece.val]
        entry_f[sqr_to] += mul_common - fld(entry_f[sqr_to] * abs_δ, HistoryDivide)
    end
    return
end

# Update the history heuristics for the case where we have a counter move, and no follow move.
function updatehistory_internal!(thread::Thread, quietstried::MoveStack, bonus::Int, best_move::Move, colour::UInt8)::Nothing
    abs_δ = abs(bonus)
    @inbounds for move in quietstried
        δ = (move == best_move) ? bonus : -bonus
        sqr_to = to(move)
        sqr_from = from(move)

        mul_common = HistoryMultiply * δ

        # Update the butterfly table. (History heuristics)
        entry = thread.history[colour][sqr_from]
        entry[sqr_to] += mul_common - fld(entry[sqr_to] * abs_δ, HistoryDivide)
    end
    return
end

"""
    gethistoryscores!(thread::Thread, moves::MoveStack, scores::Vector{Int}, idx_start::Int, idx_end::Int, ply::Int)

A function to extract the scores of quiet moves using the history heuristics. `idx_start` and `idx_end` determing the starting and ending indexes of the quiet moves in `moves`.
"""
function gethistoryscores!(thread::Thread, moves::MoveStack, scores::Vector{Int32}, idx_start::Int, idx_end::Int, ply::Int)::Nothing
    board = thread.board

    # extract one move ago
    if ply > 0
        @inbounds counter = thread.movestack[ply]
        @inbounds cm_piece = thread.piecestack[ply]
        cm_to = to(counter)
    else
        counter = MOVE_NONE
        cm_piece = VOID
        cm_to = zero(UInt16)
    end

    # extract move from two moves ago
    if ply > 1
        @inbounds follow = thread.movestack[ply - 1]
        @inbounds fm_piece = thread.piecestack[ply - 1]
        fm_to = to(follow)
    else
        follow = MOVE_NONE
        fm_piece = VOID
        fm_to = zero(UInt16)
    end

    # Below, the decision to split the cases up into 4 functions is so that the compiler can allow for SIMD optimisations.
    if (counter !== MOVE_NONE) && (counter !== NULL_MOVE) && (follow !== MOVE_NONE) && (follow !== NULL_MOVE)
        gethistoryscores_internal_histcountfollow!(thread, moves, scores, idx_start, idx_end, counter, cm_piece, cm_to, follow, fm_piece, fm_to)
        return
    elseif (counter !== MOVE_NONE) && (counter !== NULL_MOVE)
        gethistoryscores_internal_histcount!(thread, moves, scores, idx_start, idx_end, counter, cm_piece, cm_to)
        return
    elseif (follow !== MOVE_NONE) && (follow !== NULL_MOVE)
        gethistoryscores_internal_histfollow!(thread, moves, scores, idx_start, idx_end, follow, fm_piece, fm_to)
        return
    else
        gethistoryscores_internal_history!(thread, moves, scores, idx_start, idx_end)
        return
    end
end

# Internals for the case where we have no counter or follow up move.
function gethistoryscores_internal_history!(thread::Thread, moves::MoveStack, scores::Vector{Int32}, idx_start::Int, idx_end::Int)::Nothing
    @inbounds thist = thread.history[thread.board.turn.val]
    @inbounds for i in idx_start:idx_end
        # Extract useful move information.
        sqr_to = to(moves[i])
        sqr_from = from(moves[i])
        # Add the scores
        scores[i] = thist[sqr_from][sqr_to]
    end
    return
end

# Internals for the case where we have a counter move.
function gethistoryscores_internal_histcount!(thread::Thread, moves::MoveStack, scores::Vector{Int32}, idx_start::Int, idx_end::Int,
    counter::Move, cm_piece::PieceType, cm_to::Integer)::Nothing
    @inbounds thist = thread.history[thread.board.turn.val]
    @inbounds tchist = thread.counterhistory[cm_piece.val][cm_to]
    @inbounds for i in idx_start:idx_end
        # Extract useful move information.
        sqr_to = to(moves[i])
        sqr_from = from(moves[i])
        move_piece = type(thread.board[sqr_from])
        # Add the scores
        scores[i] = thist[sqr_from][sqr_to]
        scores[i] += tchist[move_piece.val][sqr_to]
    end
    return
end

# Internals for the case where we have a counter and a follow up move.
function gethistoryscores_internal_histcountfollow!(thread::Thread, moves::MoveStack, scores::Vector{Int32}, idx_start::Int, idx_end::Int,
    counter::Move, cm_piece::PieceType, cm_to::Integer, follow::Move, fm_piece::PieceType, fm_to::Integer)::Nothing
    @inbounds thist = thread.history[thread.board.turn.val]
    @inbounds tchist = thread.counterhistory[cm_piece.val][cm_to]
    @inbounds tfhist = thread.followhistory[fm_piece.val][fm_to]
    @inbounds for i in idx_start:idx_end
        # Extract useful move information.
        sqr_to = to(moves[i])
        sqr_from = from(moves[i])
        move_piece = type(thread.board[sqr_from])
        # Add the scores
        scores[i] = thist[sqr_from][sqr_to]
        scores[i] += tchist[move_piece.val][sqr_to]
        scores[i] += tfhist[move_piece.val][sqr_to]
    end
    return
end

# Internals for the case where we have a follow up move.
function gethistoryscores_internal_histfollow!(thread::Thread, moves::MoveStack, scores::Vector{Int32}, idx_start::Int, idx_end::Int,
    follow::Move, fm_piece::PieceType, fm_to::Integer)::Nothing
    @inbounds thist = thread.history[thread.board.turn.val]
    @inbounds tfhist = thread.followhistory[fm_piece.val][fm_to]
    @inbounds for i in idx_start:idx_end
        # Extract useful move information.
        sqr_to = to(moves[i])
        sqr_from = from(moves[i])
        move_piece = type(thread.board[sqr_from])
        # Add the scores
        scores[i] = thist[sqr_from][sqr_to]
        scores[i] += tfhist[move_piece.val][sqr_to]
    end
    return
end

function gethistory(thread::Thread, move::Move, ply::Int)
    board = thread.board

    sqr_to = to(move)
    sqr_from = from(move)
    move_piece = type(thread.board[sqr_from])

    # extract one move ago
    if ply > 0
        @inbounds counter = thread.movestack[ply]
        @inbounds cm_piece = thread.piecestack[ply]
        cm_to = to(counter)
    else
        counter = MOVE_NONE
        cm_piece = VOID
        cm_to = zero(UInt16)
    end

    # extract move from two moves ago
    if ply > 1
        @inbounds follow = thread.movestack[ply - 1]
        @inbounds fm_piece = thread.piecestack[ply - 1]
        fm_to = to(follow)
    else
        follow = MOVE_NONE
        fm_piece = VOID
        fm_to = zero(UInt16)
    end

    @inbounds hist = thread.history[thread.board.turn.val][sqr_from][sqr_to]

    if (counter == MOVE_NONE || counter == NULL_MOVE)
        cmhist = 0
    else
        @inbounds cmhist = thread.counterhistory[cm_piece.val][cm_to][move_piece.val][sqr_to]
    end

    if (follow == MOVE_NONE || follow == NULL_MOVE)
        fmhist = 0
    else
        @inbounds fmhist = thread.followhistory[fm_piece.val][fm_to][move_piece.val][sqr_to]
    end

    hist, cmhist, fmhist
end

max_stars_repo_name: J-Revell/Astellarn

max_issues_repo_name: J-Revell/Astellarn

max_forks_repo_name: J-Revell/Astellarn


text:
using DandelionWebSockets
using MsgPack
using UUIDs: uuid4

import DandelionWebSockets: on_text, on_binary
import DandelionWebSockets: state_connecting, state_open,
                            state_closing, state_closed
import DandelionWebSockets: getrequestheaders, HTTPHandshakeLogic

# include("CPU.jl")

get_request_id() = string(uuid4())[1:7]
client_id = get_request_id()

Base.@kwdef mutable struct PlutoClient <: WebSocketHandler
        connection::Union{WebSocketConnection, Nothing}
        stop_channel::Channel{Any}
        notebook_id::String
end
PlutoClient(;notebook_id) = PlutoClient(nothing, Channel{Any}(3), notebook_id)

handle_message(handler::PlutoClient, message) = nothing

on_text(handler::PlutoClient, text::String) = println("received: ", text)
function on_binary(handler::PlutoClient, data::Vector{UInt8})
        println("received binary data")

        message = unpack(data)

        handle_message(handler, message)
end
function send_msg(handler::PlutoClient, payload)
        msg = merge(payload, Dict("request_id" => get_request_id(), "client_id" => client_id)) |> pack
        send_binary(handler.connection, msg)
end

function ping(handler::PlutoClient)
        send_msg(handler, Dict("request_id" => get_request_id(), "type" => "ping", "client_id" => client_id, "body" => Dict()))
end
state_connecting(handler::PlutoClient, connection::WebSocketConnection) =
        handler.connection = connection

set_bond_value(handler::PlutoClient, path, value) =
        send_msg(handler, Dict(
                "notebook_id" => handler.notebook_id,
                "type" => "update_notebook",
                "body" => Dict("updates" => [Dict("value" => Dict("value" => value), "op" => "replace", "path" => split(path, "/"))])))

function state_open(handler::PlutoClient)
        @async begin
                send_msg(handler, Dict(
                        "body" => Dict(),
                        "type" => "connect",
                        "notebook_id" => client.notebook_id))
                ping(handler)
                send_msg(handler, Dict(
                        "body" => Dict("updates" => []),
                        "type" => "update_notebook",
                        "notebook_id" => client.notebook_id))

                while true
                        ping(handler)
                        sleep(28)
                end
        end
end
state_closing(handler::PlutoClient) = @info "Closing connection..."
state_closed(handler::PlutoClient) = put!(handler.stop_channel, true)

function connect(handler::PlutoClient, url)
        wsclient = WSClient()
        url = "ws://$(url)/?id=$(handler.notebook_id)"
        wsconnect(wsclient, url, handler)
end


text:
using LinearAlgebra
using SparseArrays
using SuiteSparse

struct IndGraphSparse{T <: RealOrComplex, Ti} <: IndGraph
  m::Int
  n::Int
  A::SparseMatrixCSC{T, Ti}
  F::SuiteSparse.CHOLMOD.Factor{T} #LDL factorization

  tmp::Array{T, 1}
  tmpx::SubArray{T, 1, Array{T, 1}, Tuple{UnitRange{Int}}, true}
  res::Array{T, 1}
end

function IndGraphSparse(A::SparseMatrixCSC{T,Ti}) where {T <: RealOrComplex, Ti}

  m, n = size(A)
  K = [SparseMatrixCSC{T}(I, n, n) A'; A -SparseMatrixCSC{T}(I, m, m)]

  F = ldlt(K)

  tmp = Array{T, 1}(undef, m + n)
  tmpx = view(tmp, 1:n)
  tmpy = view(tmp, (n + 1):(n + m)) #second part is always zeros
  fill!(tmpy, 0)

  res = Array{T,1}(undef, m + n)
  return IndGraphSparse(m, n, A, F, tmp, tmpx, res)
end

function prox!(
    x::AbstractVector{T},
    y::AbstractVector{T},
    f::IndGraphSparse,
    c::AbstractVector{T},
    d::AbstractVector{T},
    gamma=1.0
  ) where {T <: RealOrComplex}
  #instead of res = [c + f.A' * d; zeros(f.m)]
  mul!(f.tmpx, adjoint(f.A), d)
  f.tmpx .+= c
  # A_ldiv_B!(f.res, f.F, f.tmp) #is not working
  f.res .= f.F \ f.tmp
  # f.res .= f.F \ f.tmp #note here f.tmp which is m+n array
  copyto!(x, 1, f.res, 1, f.n)
  copyto!(y, 1, f.res, f.n + 1, f.m)
  return 0.0
end

function (f::IndGraphSparse)(x::AbstractVector{T}, y::AbstractVector{T}) where
    {T <: RealOrComplex}
  # the tolerance in the following line should be customizable
  tmpy = view(f.res, 1:f.m) # the res is rewritten in prox!
  mul!(tmpy, f.A, x)
  tmpy .-= y
  if norm(tmpy, Inf) <= 1e-12
    return 0.0
  end
  return +Inf
end

fun_name(f::IndGraphSparse) = "Indicator of an operator graph defined by sparse matrix"
# fun_dom(f::IndGraph) = "AbstractArray{Real,1}, AbstractArray{Complex,1}"
# fun_expr(f::IndGraph) = "x,y ↦ 0 if Ax = y, +∞ otherwise"
# fun_params(f::IndGraph) =
#   string( "A = ", typeof(f.A), " of size ", size(f.A))

function prox_naive(
    f::IndGraphSparse,
    c::AbstractVector{T},
    d::AbstractVector{T},
    gamma=1.0
  ) where {T <: RealOrComplex}

  tmp = f.A'*d
  tmp .+= c
  res = [tmp; zeros(f.m)]
  xy = f.F \ res
  return xy[1:f.n], xy[f.n + 1:f.n + f.m], 0.0
end

function (f::IndGraphSparse)(xy::AbstractVector{T}) where
  {T <: RealOrComplex}
  x, y = splitinput(f, xy)
  return f(x, y)
end

(f::IndGraphSparse)(xy::Tuple{AbstractVector{T}, AbstractVector{T}}) where
  {T <: RealOrComplex} = f(xy[1], xy[2])


text:
#=
Convert extensive form to matrix form
=#

infodict(sol::AbstractCFRSolver{K,G,I}) where {K,G,I} = getfield(sol,:I)::Dict{K,I}
strategy(I::AbstractInfoState) = getfield(I,:σ)::Vector{Float64}

function pure_strategies(sol::CFRSolver, p::Int)
    game = sol.game
    K = infokeytype(sol.game)
    strategies = Vector{Pair{K,Vector{Float64}}}[]
    key_map = Tuple{K,Int}[]
    for (k,I) in infodict(sol)
        if player(game,k) === p
            push!(key_map, (k, length(I.σ)))
        end
    end
    add_next!(strategies, Pair{K,Vector{Float64}}[], key_map)

    return strategies
end

function add_next!(
        strategies::Vector{Vector{Pair{K,Vector{Float64}}}},
        cur_vec::Vector{Pair{K,Vector{Float64}}},
        key_map::Vector{Tuple{K,Int}} ) where K

    L = length(cur_vec)
    Lf = length(key_map)
    if L === Lf
        push!(strategies, copy(cur_vec))
    else
        key, a_length = key_map[L+1]
        for i in 1:a_length
            v = copy(cur_vec)
            σ = zeros(Float64, a_length)
            σ[i] = 1.0
            push!(v, (key => σ))
            add_next!(strategies, v, key_map)
        end
    end
end

function Base.Matrix(game::Game)
    # find all info states
    solver = CFRSolver(game)
    train!(solver, 1)

    # enumerate all pure strategies
    p1_strats = pure_strategies(solver, 1)
    p2_strats = pure_strategies(solver, 2)

    # evaluate all pure strategies
    mat = Matrix{NTuple{2,Float64}}(undef,length(p1_strats), length(p2_strats))
    for (i,s1) in enumerate(p1_strats)
        for (j,s2) in enumerate(p2_strats)
            for (k1,σ1) in s1
                solver.I[k1].s .= σ1
            end
            for (k2,σ2) in s2
                solver.I[k2].s .= σ2
            end
            eval = evaluate(solver)
            mat[i,j] = eval
        end
    end

    return mat
end

max_stars_repo_name: WhiffleFish/HelloCFR

max_issues_repo_name: WhiffleFish/HelloCFR

max_forks_repo_name: WhiffleFish/HelloCFR


text:
arf = inv(ArbFloat(3))
arb = inv(ArbReal(3))
acb = inv(ArbComplex(3))

@testset "Arithmetic" begin
    @test Float64(arf + 2*arf) == 1.0
    @test Float64(arb + 2*arb) == 1.0
    @test Float64(acb + 2*acb) == 1.0
    @test Float64(arf / arf) == 1.0
    @test Float64(arb / arb) == 1.0
    @test Float64(acb / acb) == 1.0
end


text:
# ---
# title: 1046. Last Stone Weight
# id: problem1046
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Heap, Greedy
# link: <https://leetcode.com/problems/last-stone-weight/description/>
# hidden: true
# ---
#
# We have a collection of stones, each stone has a positive integer weight.
#
# Each turn, we choose the two **heaviest**  stones and smash them together.
# Suppose the stones have weights `x` and `y` with `x <= y`.  The result of this
# smash is:
#
#   * If `x == y`, both stones are totally destroyed;
#   * If `x != y`, the stone of weight `x` is totally destroyed, and the stone of weight `y` has new weight `y-x`.
#
# At the end, there is at most 1 stone left.  Return the weight of this stone
# (or 0 if there are no stones left.)
#
#
#
# **Example 1:**
#
#
#
#     Input: [2,7,4,1,8,1]
#     Output: 1
#     Explanation:
#     We combine 7 and 8 to get 1 so the array converts to [2,4,1,1,1] then,
#     we combine 2 and 4 to get 2 so the array converts to [2,1,1,1] then,
#     we combine 2 and 1 to get 1 so the array converts to [1,1,1] then,
#     we combine 1 and 1 to get 0 so the array converts to [1] then that's the value of last stone.
#
#
#
# **Note:**
#
#   1. `1 <= stones.length <= 30`
#   2. `1 <= stones[i] <= 1000`
#
#
## @lc code=start
using LeetCode

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


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""RunAsGroupStrategyOptions defines the strategy type and any options used to create the strategy.

    IoK8sApiPolicyV1beta1RunAsGroupStrategyOptions(;
        ranges=nothing,
        rule=nothing,
    )

    - ranges::Vector{IoK8sApiPolicyV1beta1IDRange} : ranges are the allowed ranges of gids that may be used. If you would like to force a single gid then supply a single range with the same start and end. Required for MustRunAs.
    - rule::String : rule is the strategy that will dictate the allowable RunAsGroup values that may be set.
"""
mutable struct IoK8sApiPolicyV1beta1RunAsGroupStrategyOptions <: SwaggerModel
    ranges::Any # spec type: Union{ Nothing, Vector{IoK8sApiPolicyV1beta1IDRange} } # spec name: ranges
    rule::Any # spec type: Union{ Nothing, String } # spec name: rule

    function IoK8sApiPolicyV1beta1RunAsGroupStrategyOptions(;ranges=nothing, rule=nothing)
        o = new()
        validate_property(IoK8sApiPolicyV1beta1RunAsGroupStrategyOptions, Symbol("ranges"), ranges)
        setfield!(o, Symbol("ranges"), ranges)
        validate_property(IoK8sApiPolicyV1beta1RunAsGroupStrategyOptions, Symbol("rule"), rule)
        setfield!(o, Symbol("rule"), rule)
        o
    end
end # type IoK8sApiPolicyV1beta1RunAsGroupStrategyOptions

const _property_map_IoK8sApiPolicyV1beta1RunAsGroupStrategyOptions = Dict{Symbol,Symbol}(Symbol("ranges")=>Symbol("ranges"), Symbol("rule")=>Symbol("rule"))
const _property_types_IoK8sApiPolicyV1beta1RunAsGroupStrategyOptions = Dict{Symbol,String}(Symbol("ranges")=>"Vector{IoK8sApiPolicyV1beta1IDRange}", Symbol("rule")=>"String")
Base.propertynames(::Type{ IoK8sApiPolicyV1beta1RunAsGroupStrategyOptions }) = collect(keys(_property_map_IoK8sApiPolicyV1beta1RunAsGroupStrategyOptions))
Swagger.property_type(::Type{ IoK8sApiPolicyV1beta1RunAsGroupStrategyOptions }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiPolicyV1beta1RunAsGroupStrategyOptions[name]))}
Swagger.field_name(::Type{ IoK8sApiPolicyV1beta1RunAsGroupStrategyOptions }, property_name::Symbol) =  _property_map_IoK8sApiPolicyV1beta1RunAsGroupStrategyOptions[property_name]

function check_required(o::IoK8sApiPolicyV1beta1RunAsGroupStrategyOptions)
    (getproperty(o, Symbol("rule")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiPolicyV1beta1RunAsGroupStrategyOptions }, name::Symbol, val)
end


text:
using RandomUseful
using Test

include("test_fourier_diff.jl")
include("test_cheb_diff.jl")
include("test_delaydifferentialequation.jl")


text:
const plotsock = ZMQ.Socket(ctx, ZMQ.SUB)
ZMQ.subscribe(plotsock)
ZMQ.connect(plotsock, "tcp://127.0.0.1:50002")
Plots.plotlyjs()

const context = Dict{Symbol,Any}()

function hear(context)
    while true
        msg = ZMQ.recv(s)
        out = convert(IOStream, msg)
        seekstart(out)
        d = deserialize(out)
        handle(d, context)
    end
end

function handle(x::PlotSetup, context)
    context[:setup] = x
    A = (x.arrtype)(x.size...)
    A[:] = zero(eltype(A))
    context[:array] = A
    context[:plot] = plot(A, show=true)
    yield()
end

function handle(pt::PlotPoint, context)
    A = context[:array]
    A[pt.inds...] = pt.v

    p = context[:plot]
    x = p[1][1][:x]
    p[1] = (x, A)
end


text:
"""
Adapted from sandreza/Learning/sandbox/gaussian_process.jl
https://github.com/sandreza/Learning/blob/master/sandbox/gaussian_process.jl
Changed handling of kernel functions, some variable names, and added log marginal likelihood function.
"""

using LinearAlgebra

"""
    GP{K, 𝒮, 𝒮2, 𝒰, 𝒱, 𝒜}

Data structure for typical Gaussian process regression computations.
"""
struct GP{K, 𝒮, 𝒮2, 𝒰, 𝒱, 𝒜}
     kernel :: K  # a Kernel object
    x_train :: 𝒮  # an array of vectors (n-length array of D-length vectors)
          α :: 𝒮2 # an array
          K :: 𝒰  # matrix or sparse matrix
         CK :: 𝒱  # cholesky factorization of K
      cache :: 𝒜
end

"""
    GPmodel(x_train, y_train, kernel, zavg; sparsity_threshold = 0.0, robust = true, entry_threshold = sqrt(eps(1.0)))

Constructs the posterior distribution for a Gaussian process. In other words this does the "training" automagically.

# Arguments
- `x_train`: (array). training inputs (predictors), must be an array of states.
                      length-n array of D-length vectors, where D is the length of each input n is the number of training points.
- `y_train`: (array). training outputs (prediction), must have the same number as x_train
                      length-n array of D-length vectors.
- `kernel`: (Kernel). Kernel object. See kernels.jl.
                      kernel_function(kernel)(x,x') maps predictor x predictor to real numbers.
# Keyword Arguments
- `z`: (vector). values w.r.t. which to derivate the state (default none).

- `normalize`: (bool). whether to normalize the data during preprocessing and reverse the scaling for postprocessing. Can lead to better performance.

- `hyperparameters`: (array). default = []. hyperparameters that enter into the kernel

- `sparsity_threshold`: (number). default = 0.0. a number between 0 and 1 that determines when to use sparse array format. The default is to never use it

- `robust`: (bool). default = true. This decides whether to uniformly scale the diagonal entries of the Kernel Matrix. This sometimes helps with Cholesky factorizations.

- `entry_threshold`: (number). default = sqrt(eps(1.0)). This decides whether an entry is "significant" or not. For typical machines this number will be about 10^(-8) * largest entry of kernel matrix.
"""
function GPmodel(x_train, y_train, kernel, zavg; sparsity_threshold = 0.0, robust = true, entry_threshold = sqrt(eps(1.0)))

    # get k(x,x') function from kernel object
    kernel = kernel_function(kernel; z=zavg)

    # fill kernel matrix with values
    K = compute_kernel_matrix(kernel, x_train)

    # get the maximum entry for scaling and sparsity checking
    mK = maximum(K)

    # make Cholesky factorization work by adding a small amount to the diagonal
    if robust
        K += mK*sqrt(eps(1.0))*I
    end

    # check sparsity
    bools = K .> entry_threshold * mK
    sparsity = sum(bools) / length(bools)
    if sparsity < sparsity_threshold
        sparse_K = similar(K) .* 0
        sparse_K[bools] = sK[bools]
        K = sparse(Symmetric(sparse_K))
        CK = cholesky(K)
    else
        CK = cholesky(K)
    end

    # get prediction weights FIX THIS SO THAT IT ALWAYS WORKS
    y = hcat(y_train...)'
    α = CK \ y # α = K + σ_noise*I

    # construct struct
    return GP(kernel, x_train, α', K, Array(CK), zeros(length(x_train)))
end

"""
    GPmodel(𝒱; kernel::Kernel = Kernel())

Create an instance of GP using data from ProfileData object 𝒟.

# Arguments
- 𝒱::VData, Data for training the GP
"""
function GPmodel(𝒱; kernel::Kernel = Kernel())
    x_train = [pair[1] for pair in 𝒱.training_data]
    y_train = [pair[2] for pair in 𝒱.training_data]

    # create instance of GP using data from ProfileData object
    return GPmodel(x_train, y_train, kernel, 𝒱.z);
end

"""
prediction(x, 𝒢::GP)

Given state x, GP 𝒢, returns the mean GP prediction

# Arguments
- `x`: single scaled state
- `𝒢`: GP object with which to make the prediction

# Return
- `y`: scaled prediction
"""
function model_output(x, 𝒢::GP)
    for i in 1:length(𝒢.cache)
        𝒢.cache[i] = 𝒢.kernel(x, 𝒢.x_train[i])
    end
    return 𝒢.α * 𝒢.cache
end

"""
    uncertainty(x, 𝒢::GP)

Given state x and GP 𝒢, output the variance at a point

# Arguments
- `x`: state

# Return
- `var`: variance
"""
function uncertainty(x, 𝒢::GP)
    tmpv = zeros(size(𝒢.x_train)[1])
    for i in eachindex(𝒢.x_train)
        tmpv[i] = 𝒢.kernel(x, 𝒢.x_train[i])
    end
    # no ldiv for suitesparse
    tmpv2 = 𝒢.CK \ tmpv
    var = k(x, x) .- tmpv'*tmpv2
    return var
end

"""
    compute_kernel_matrix(kernel, x)

Computes the kernel matrix for GPR.

# Arguments
- `k` : (Kernel) kernel function k(a,b).
- `x` : (array of predictors). x[1] is a vector

# Return
- `sK`: (symmetric matrix). A symmetric matrix with entries sK[i,j] = k(x[i], x[j]). This is only meaningful if k(x,y) = k(y,x) (it should)
"""
function compute_kernel_matrix(k, x)
    K = [k(x[i], x[j]) for i in eachindex(x), j in eachindex(x)]

    if typeof(K[1,1]) <: Number
        sK = Symmetric(K)
    else
        sK = K
    end
    return sK
end

"""
    mean_log_marginal_loss(y_train, 𝒢::GP; add_constant=false)

Computes log marginal loss for each element in the output and averages the results.

Assumes noise-free observations

    log(p(y|X)) = -(1/2) * (y'*α + 2*sum(Diagonal(CK)) + n*log(2*pi))

where n is the number of training points and

# Arguments
- `y_train`: (Array). training outputs (prediction), must have the same number as x_train
- `𝒢`: (GP).

# Keyword Arguments
- `add_constant`: (bool). whether to give the exact value of the loss or leave out an added constant for efficiency.
"""
function mean_log_marginal_loss(y_train, 𝒢::GP; add_constant=false)
    n = length(𝒢.x_train)
    D = length(𝒢.x_train[1])

    ys = hcat(y_train...)' # n x D

    if add_constant
        c = sum([log(𝒢.CK[i,i]) for i in 1:n]) + 0.5*n*log(2*pi)
        total_loss=0.0
        for i in 1:D
            total_loss -= 0.5*ys[:,i]'*𝒢.α[:,i] + c
        end
    else
        total_loss=0.0
        for i in 1:D
            total_loss -= 0.5*ys[:,i]'*𝒢.α[:,i]
        end
    end

    return total_loss / D
end


text:
module ModelBuilds

using TestSpecDef, JuMP, CPLEX, Dictionaries, NamedDims

export arrayModel, arrayModelLongName, arrayModelLongBasename, denseAxisArrayModel, denseAxisArrayModelLongBasename, denseAxisArrayModelLongName, denseAxisArrayModelLongNameShortBasename, denseAxisArrayModelStr20, dictModel, dictModelStr, dictionaryModel, namedDimsModel

#core function that includes model build and optimization for various problem specifications
function arrayModel(coeff::Array{Float64,2}, bound::Array{Float64,1}, n::Int64, directMode::Bool, optimization::Bool, optimizer::Optimizers)
	#DataType A

	#coeff 				matrix: for defintion of the linear constraints of the form coeff*x .<= bound
	#bound 				array: for definition of constraints (see above)
	#n 					integer: showing problem size
	#directMode 		bool: info for model build
	#optimization		bool: shall the JuMP model be optimized

	if directMode == false
		m = Model();
		set_optimizer(m, CPLEX.Optimizer)

	else
		m = direct_model(CPLEX.Optimizer());
	end

	@variable(m,0<= x[1:n]<=1)
	@constraint(m, con[i = 1:n], sum(coeff[i,j]*x[j] for j = 1:n) <= bound[i] )

	@objective(m, Max, sum(x))

	if optimization == true
		JuMP.optimize!(m)
	end

	return(m)

end

function arrayModelLongName(coeff::Array{Float64,2}, bound::Array{Float64,1}, n::Int64, directMode::Bool, optimization::Bool)
	#DataType AL

	#coeff 				matrix: for defintion of the linear constraints of the form coeff*x .<= bound
	#bound 				array: for definition of constraints (see above)
	#n 					integer: showing problem size
	#directMode 		bool: info for model build
	#optimization		bool: shall the JuMP model be optimized

	if directMode == false
		m = Model();
		set_optimizer(m, CPLEX.Optimizer)
	else
		m = direct_model(CPLEX.Optimizer());
	end

	@variable(m, 0<= verylooooooooongname[1:n]<=1,base_name = "x")
	@constraint(m, con[i = 1:n], sum(coeff[i,j]*verylooooooooongname[j] for j in 1:n) <= bound[i] )

	@objective(m, Max, sum(verylooooooooongname))

	if optimization == true
		JuMP.optimize!(m)
	end

	return(m)

end

function arrayModelLongBasename(coeff::Array{Float64,2}, bound::Array{Float64,1}, n::Int64, directMode::Bool, optimization::Bool)
	#DataType ALB

	#coeff 				matrix: for defintion of the linear constraints of the form coeff*x .<= bound
	#bound 				array: for definition of constraints (see above)
	#n 					integer: showing problem size
	#directMode 		bool: info for model build
	#optimization		bool: shall the JuMP model be optimized

	if directMode == false
		m = Model();
		et_optimizer(m, CPLEX.Optimizer)
	else
		m = direct_model(CPLEX.Optimizer());
	end

	@variable(m, 0<= x[1:n]<=1,base_name = "VERYLOOOOOOOOONGNAME")
	@constraint(m, con[i = 1:n], sum(coeff[i,j]*x[j] for j in 1:n) <= bound[i] )

	obj = @objective(m, Max, sum(x))

	if optimization == true
		JuMP.optimize!(m)
	end

	return(m)

end

function denseAxisArrayModel(coeff::JuMP.Containers.DenseAxisArray{Float64,2,Tuple{UnitRange{Int64},UnitRange{Int64}},Tuple{Dict{Int64,Int64},Dict{Int64,Int64}}},
	bound::JuMP.Containers.DenseAxisArray{Float64,1,Tuple{UnitRange{Int64}},Tuple{Dict{Int64,Int64}}}, n::Int64, directMode::Bool, optimization::Bool)
	#DataType DAA

	#coeff 				matrix: for defintion of the linear constraints of the form coeff*x .<= bound
	#bound 				array: for definition of constraints (see above)
	#n 					integer: showing problem size
	#directMode 		bool: info for model build
	#optimization		bool: shall the JuMP model be optimized

	if directMode == false
		m = Model();
		set_optimizer(m, CPLEX.Optimizer)
	else
		m = direct_model(CPLEX.Optimizer());
	end

	@variable(m,0<= x[1:n]<=1, container = DenseAxisArray)
	@constraint(m, con[i = 1:n], sum(coeff[i,j]*x[j] for j in 1:n) <= bound[i], container = DenseAxisArray)

	@objective(m, Max, sum(x[j] for j in 1:n))

	if optimization == true
		JuMP.optimize!(m)
	end

	return(m)
end

function denseAxisArrayModelLongName(coeff::JuMP.Containers.DenseAxisArray{Float64,2,Tuple{UnitRange{Int64},UnitRange{Int64}},Tuple{Dict{Int64,Int64},Dict{Int64,Int64}}},
	bound::JuMP.Containers.DenseAxisArray{Float64,1,Tuple{UnitRange{Int64}},Tuple{Dict{Int64,Int64}}}, n::Int64, directMode::Bool, optimization::Bool)
	#DataType DAAL

	#coeff 				matrix: for defintion of the linear constraints of the form coeff*x .<= bound
	#bound 				array: for definition of constraints (see above)
	#n 					integer: showing problem size
	#directMode 		bool: info for model build
	#optimization		bool: shall the JuMP model be optimized

	if directMode == false
		m = Model();
		set_optimizer(m, CPLEX.Optimizer)
	else
		m = direct_model(CPLEX.Optimizer());
	end

	@variable(m,0<= verylooooooongname_x[1:n]<=1, container = DenseAxisArray)

	@constraint(m, verylooooongname_con[i=1:n], sum(coeff[i,j]*verylooooooongname_x[j] for j in 1:n) <= bound[i], container = DenseAxisArray)

	@objective(m, Max, sum(verylooooooongname_x[j] for j in 1:n))

	if optimization == true
		JuMP.optimize!(m)
	end

	return(m)
end

function denseAxisArrayModelLongNameShortBasename(coeff::JuMP.Containers.DenseAxisArray{Float64,2,Tuple{UnitRange{Int64},UnitRange{Int64}},Tuple{Dict{Int64,Int64},Dict{Int64,Int64}}},
	bound::JuMP.Containers.DenseAxisArray{Float64,1,Tuple{UnitRange{Int64}},Tuple{Dict{Int64,Int64}}}, n::Int64, directMode::Bool, optimization::Bool)
	#DataType DAALS

	#coeff 				matrix: for defintion of the linear constraints of the form coeff*x .<= bound
	#bound 				array: for definition of constraints (see above)
	#n 					integer: showing problem size
	#directMode 		bool: info for model build
	#optimization		bool: shall the JuMP model be optimized

	if directMode == false
		m = Model();
		set_optimizer(m, CPLEX.Optimizer)
	else
		m = direct_model(CPLEX.Optimizer());
	end

	@variable(m,0<= verylooooooongname_x[1:n]<=1, container = DenseAxisArray, base_name = "x")

	@constraint(m, verylooooongname_con[i=1:n], sum(coeff[i,j]*verylooooooongname_x[j] for j in 1:n) <= bound[i], container = DenseAxisArray)

	@objective(m, Max, sum(verylooooooongname_x[j] for j in 1:n))

	if optimization == true
		JuMP.optimize!(m)
	end

	return(m)
end

function denseAxisArrayModelLongBasename(coeff::JuMP.Containers.DenseAxisArray{Float64,2,Tuple{UnitRange{Int64},UnitRange{Int64}},Tuple{Dict{Int64,Int64},Dict{Int64,Int64}}},
	bound::JuMP.Containers.DenseAxisArray{Float64,1,Tuple{UnitRange{Int64}},Tuple{Dict{Int64,Int64}}}, n::Int64, directMode::Bool, optimization::Bool)
	#DAALB

	#coeff 				matrix: for defintion of the linear constraints of the form coeff*x .<= bound
	#bound 				array: for definition of constraints (see above)
	#n 					integer: showing problem size
	#directMode 		bool: info for model build
	#optimization		bool: shall the JuMP model be optimized

	if directMode == false
		m = Model();
		set_optimizer(m, CPLEX.Optimizer)
	else
		m = direct_model(CPLEX.Optimizer());
	end

	@variable(m,0<= x[1:n]<=1, container = DenseAxisArray, base_name = "VERYLOOOOOOOOONGNAME")

	@constraint(m, con[i=1:n], sum(coeff[i,j]*x[j] for j in 1:n) <= bound[i], container = DenseAxisArray)

	@objective(m, Max, sum(x[j] for j in 1:n))

	if optimization == true
		JuMP.optimize!(m)
	end

	return(m)
end

function denseAxisArrayModelStr20(coeff::JuMP.Containers.DenseAxisArray{Float64,2,Tuple{Array{String,1},Array{String,1}},Tuple{Dict{String,Int64},Dict{String,Int64}}},
	bound::JuMP.Containers.DenseAxisArray{Float64,1,Tuple{Array{String,1}},Tuple{Dict{String,Int64}}},
	str_obj::Array{String,1}, str_con::Array{String,1}, directMode::Bool, optimization::Bool)
	#DataType DAAS20

	#coeff 				matrix: for defintion of the linear constraints of the form coeff*x .<= bound
	#bound 				array: for definition of constraints (see above)
	#n 					integer: showing problem size
	#directMode 		bool: info for model build
	#optimization		bool: shall the JuMP model be optimized

	if directMode == false
		m = Model();
		set_optimizer(m, CPLEX.Optimizer)

	else
		m = direct_model(CPLEX.Optimizer());
	end

	@variable(m,0<= x[str_obj]<=1)

	@constraint(m, con[name_con in str_con], sum(coeff[name_con,name_obj]*x[name_obj] for name_obj in str_obj) <= bound[name_con])

	@objective(m, Max, sum(x[name_obj] for name_obj in str_obj))

	if optimization == true
		JuMP.optimize!(m)
	end

	return(m)
end

function dictModel(coeff::Dict{Any,Any}, bound::Dict{Any,Any}, n::Int64, directMode::Bool, optimization::Bool)
	#DataType D

	#coeff 				matrix: for defintion of the linear constraints of the form coeff*x .<= bound
	#bound 				array: for definition of constraints (see above)
	#n 					integer: showing problem size
	#directMode 		bool: info for model build
	#optimization		bool: shall the JuMP model be optimized

	if directMode == false
		m = Model();
		set_optimizer(m, CPLEX.Optimizer)
	else
		m = direct_model(CPLEX.Optimizer());
	end

	x = Dict()
	for j = 1:n
		x[j] = @variable(m, lower_bound = 0, upper_bound =1)
	end

	con = Dict()
	for i = 1:n
		con[i] = @constraint(m, sum(coeff[i,j]*x[j] for j in 1:n) <= bound[i])
	end

	@objective(m, Max, sum(x[j] for j in 1:n))

	if optimization == true
		JuMP.optimize!(m)
	end

	return(m)
end

function dictModelStr(coeff::Dict{Any,Any}, bound::Dict{Any,Any}, str_obj::Array{String,1}, str_con::Array{String,1}, directMode::Bool, optimization::Bool)
	#DataType DS

	#coeff 				matrix: for defintion of the linear constraints of the form coeff*x .<= bound
	#bound 				array: for definition of constraints (see above)
	#n 					integer: showing problem size
	#directMode 		bool: info for model build
	#optimization		bool: shall the JuMP model be optimized

	if directMode == false
		m = Model();
		set_optimizer(m, CPLEX.Optimizer)
	else
		m = direct_model(CPLEX.Optimizer());
	end

	x = Dict()
	for name_obj in str_obj
		x[name_obj] = @variable(m, lower_bound = 0, upper_bound = 1)
	end

	con = Dict()
	for name_con in str_con
		con[name_con] = @constraint(m, sum(coeff[name_con, name_obj]*x[name_obj] for name_obj in str_obj) <= bound[name_con])
	end

	@objective(m, Max, sum(x[name_obj] for name_obj in str_obj))

	if optimization == true
		JuMP.optimize!(m)
	end

	return(m)
end

function dictionaryModel(coeff::Dictionary{CartesianIndex{2},Float64}, bound::Dictionary{Int64,Float64}, n::Int64, directMode::Bool, optimization::Bool)
	#DataType DD

	#coeff 				matrix: for defintion of the linear constraints of the form coeff*x .<= bound
	#bound 				array: for definition of constraints (see above)
	#n 					integer: showing problem size
	#directMode 		bool: info for model build
	#optimization		bool: shall the JuMP model be optimized

	if directMode == false
		m = Model();
		set_optimizer(m, CPLEX.Optimizer)
	else
		m = direct_model(CPLEX.Optimizer());
	end

	x = Dictionary(1:n,@variable(m, 0 <= x[1:n] <= 1))
	con = Dictionary(1:n, @constraint(m, con[i = 1:n], sum(coeff[CartesianIndex(i ,j)]*x[j] for j in 1:n) <= bound[i]))

	@objective(m, Max, sum(x[j] for j in 1:n))

	if optimization == true
		JuMP.optimize!(m)
	end

	return(m)
end

function namedDimsModel(coeff::NamedDimsArray{(:x1, :x2),Float64,2,Array{Float64,2}}, bound::NamedDimsArray{(:x1,),Float64,1,Array{Float64,1}}, n::Int64, directMode::Bool, optimization::Bool, optimizer::Optimizers)
	#DataType ND

	#coeff 				matrix: for defintion of the linear constraints of the form coeff*x .<= bound
	#bound 				array: for definition of constraints (see above)
	#n 					integer: showing problem size
	#directMode 		bool: info for model build
	#optimization		bool: shall the JuMP model be optimized

	if directMode == false
		m = Model();
		set_optimizer(m, CPLEX.Optimizer)
	else
		m = direct_model(CPLEX.Optimizer());
	end

	x = NamedDimsArray{(:x2,)}(@variable(m, 0 <= x[1:n] <= 1))

	con = NamedDimsArray{(:x1,)}(@constraint(m, con[i = 1:n], sum(coeff[x1 = i ,x2 = j]*x[x2 = j] for j in 1:n) <= bound[i] ))

	@objective(m, Max, sum(x[j] for j in 1:n))

	if optimization == true
		JuMP.optimize!(m)
	end

	return(m)
end

end

max_stars_repo_name: SonntagK/JuMPContainerBenchmarks

max_issues_repo_name: SonntagK/JuMPContainerBenchmarks

max_forks_repo_name: SonntagK/JuMPContainerBenchmarks


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()
Pkg.update()
Pkg.add("SearchLight")
Pkg.add("SearchLightPostgreSQL")
Pkg.add("Intervals")
Pkg.add("TimeZones")
Pkg.add("Dates")
Pkg.add("ToStruct")
import BitemporalPostgres


text:
@testset "gbarray.jl" begin
    @testset "construction and getindex" begin
        @testset "dense" begin
            #Construction and indexing correct for dense matrices.
            x = rand(Float64, 100, 100)
            m = GBMatrix(x)
            @test m[50, 64] == x[50, 64]
            @test Matrix(m) == x
            #Construction and indexing correct for dense vectors
            x = rand(Bool, 10000)
            v = GBVector(x)
            m = GBMatrix(x)
            @test m[20] == x[20] == v[20]
            @test v == x

            #Indexing tests
            x = sprand(Int64, 100, 100, 0.05)
            m = GBMatrix(x)
            @test m[1, 2] === nothing
            @test m[:, 2] == GBMatrix(x[:, 2])
            @test m[2, :] == copy(GBMatrix(x[2, :])')
            @test m[:, :] == m
            @test m[1:2:5, 1:2] == GBMatrix(x[1:2:5, 1:2])
            @test m[1:2:5, :] == GBMatrix(x[1:2:5, :])

        end
        @testset "sparse" begin
            #Construction and indexing correct for sparse matrices
            x = sprand(Int32, 1000, 1000, 0.001)
            m = GBMatrix(x)
            @test nnz(x) == nnz(m)
            mnz = findnz(m)
            xnz = findnz(x)
            @test mnz == xnz

            x = sprand(UInt8, 1000, 0.05)
            v = GBVector(x)
            @test nnz(v) == nnz(x)
            vnz = findnz(v)
            xnz = findnz(x)
            @test vnz == xnz
        end
    end

    @testset "setindex and clear!" begin
        x = sprand(UInt16, 10, 10, 0.1)
        m = GBMatrix(x)
        clear!(m)
        @test nnz(m) == 0
        #steprange and scalar
        m[1:2:10, 1] = [1, 2, 3, 4, 5]
        @test m[1:2:10, 1] == GBMatrix([1,2,3,4,5])

        #range and range
        m[8:10, 8:10] = GBMatrix([[1,2,3] [4,5,6] [7,8,9]])
        @test Matrix(m[8:10, 8:10]) == [[1,2,3] [4,5,6] [7,8,9]]

        #range, range, accum, and mask
        mask = GBMatrix([[true, true, false] [false, true, true] [true, false,true]])
        m[8:10, 8:10, mask = mask, accum = BinaryOps.TIMES, desc = Descriptor(replace_output=true)] =
            fill(10, 3, 3)
        @test m[9, 10] === nothing
        @test m[10, 10] == 90

        #vectors
        x = sprand(Float32, 100, 0.35)
        v = GBVector(x)
        clear!(v)
        @test nnz(v) == 0

        #steprange
        v[10:10:100] = collect(1:10)
        @test v[100] == 10

        #steprange, mask, accum
        v[10:10:100, mask = GBVector(rand(Bool, 10)), accum = BinaryOps.ISEQ] =
        collect(1:10)
        @test v[10] == 1 && v[60] == 6 && v[100] == 1
    end
end


text:
using AlgebraicDynamics.UWDDynam
using Catlab.WiringDiagrams,  Catlab.Programs
using DelayDiffEq
using Test

const UWD = UndirectedWiringDiagram

@testset "UWDDynam" begin
  dx(x, p, t) = [x[1]^2, 2*x[1]-x[2]]
  dy(y, p, t) = [1 - y[1]^2]

  r = ContinuousResourceSharer{Float64}(2, dx)
  s = ContinuousResourceSharer{Float64}(1, dy)
  t = ContinuousResourceSharer{Float64}(1, 2, dx, [2])
  @test all(eltype.([r, s, t]) .== Float64)

  x0 = [10.0, 7.5]

  @testset "Identity & merge" begin
    d = @relation (x,y) begin
      f(x,y)
    end
    r2 = oapply(d, [r])

    @test nstates(r) == nstates(r2)
    @test nports(r) == nports(r2)
    @test portmap(r) == portmap(r2)
    @test eval_dynamics(r, x0) == eval_dynamics(r2, x0)
    @test exposed_states(r, x0) == exposed_states(r2, x0)
    @test ports(r) == [1, 2]
    @test ports(r2) == [:x, :y]

    h = 0.1
    drs = oapply(d, [euler_approx(r, h)])
    drs2 = euler_approx(r2, h)
    @test eval_dynamics(drs, x0) == eval_dynamics(drs2, x0)
    drs3 = euler_approx(r2)
    drs4 = oapply(d, [euler_approx(r)])
    @test eval_dynamics(drs, x0) == eval_dynamics(drs3, x0, [h])
    @test eval_dynamics(drs, x0) == eval_dynamics(drs4, x0, [h])

    d = @relation (x,) begin
      f(x,x)
    end
    r2 = oapply(d, [r])

    @test nstates(r2) == 1
    @test nports(r2) == 1
    @test portmap(r2) == [1]
    @test eval_dynamics(r2, [5.0]) == [30.0]
    @test exposed_states(r2, [5.0]) == [5.0]
    @test ports(r2) == [:x]

    let r = ContinuousResourceSharer{Float64}(2, (u,p,t) -> [u[1]*p[1], u[2]*p[2] + t])
    r2 = oapply(d, [r])
    @test eval_dynamics(r2, [5.0], [1.0, 2.0], 10.0) == [25.0]
    end
  end

  @testset "Copy" begin
    d = UWD(2)
    add_box!(d, 1)
    add_junctions!(d, 1)
    set_junction!(d, [1])
    set_junction!(d, [1,1], outer = true)

    let r = t
    r2 = oapply(d, r)
    @test nstates(r2) == 2
    @test nports(r2) == 2
    @test portmap(r2) == [2,2]
    @test eval_dynamics(r2, x0) == eval_dynamics(r, x0)
    @test exposed_states(r2, x0) == [x0[2], x0[2]]
    end
  end

  @testset "Copy & merge" begin
    d = @relation (x,) begin
      f(x,x)
    end

    let r = ContinuousResourceSharer{Float64}(2, 2, dx, [1,1])
    r2 = oapply(d, r)
    @test nstates(r2) == 2
    @test nports(r2) == 1
    @test ports(r2) == [:x]
    @test portmap(r2) == [1]
    @test eval_dynamics(r2, x0) == eval_dynamics(r, x0)
    @test exposed_states(r2, x0) == [x0[1]]
    end
  end

  @testset "Duplicate & add state" begin
    let r = s
    rcopy = ContinuousResourceSharer{Float64}(2, 1, dy, [1,1])
    d = HypergraphDiagram{Nothing, Symbol}(2)
    add_box!(d, 1, name = :r); add_box!(d, 2, name = :copy); add_box!(d, 1, name = :r)
    add_junctions!(d, 2)
    set_junction!(d, [1,1,2,2])
    set_junction!(d, [1,2], outer = true)
    xs = Dict(:r => r, :copy => rcopy)
    r2 = oapply(d, xs)
    @test nstates(r2) == 1
    @test nports(r2) == 2
    @test portmap(r2) == [1,1]
    @test eval_dynamics(r2, [7.0]) == [3 *(1 - 7.0^2)]
    @test exposed_states(r2, [7.0]) == [7.0, 7.0]

    d = UWD(2)
    add_box!(d, 1)
    add_junctions!(d, 2)
    set_junction!(d, [1])
    set_junction!(d, [1,2], outer = true)
    r2 = oapply(d, [r])
    @test nstates(r2) == 2
    @test nports(r2) == 2
    @test portmap(r2) == [1,2]
    @test eval_dynamics(r2, [7.0, 11.0]) == [-48.0, 0.0]
    @test exposed_states(r2, [7.0, 11.0]) == [7.0, 11.0]
    end
  end

  @testset "Boxes" begin
    xs = Dict(:r => r, :s => s)
    d = HypergraphDiagram{Nothing, Symbol}(5)
    add_box!(d, 2, name = :r); add_box!(d, 1, name = :s); add_box!(d, 2, name = :r)
    add_junctions!(d, 4)
    set_junction!(d, [1,1,1,4,2])
    set_junction!(d, [1,1,2,3,3], outer = true)
    r2 = oapply(d, xs)
    @test nstates(r2) == 4
    @test nports(r2) == 5
    @test portmap(r2) == [1,1, 3, 4,4]
    x0 = [2.0, 7.0, 3.0, 5.0]
    @test eval_dynamics(r2, x0) == [3.0, 49.0, 11.0, 0.0]
    @test exposed_states(r2, x0) == [2.0, 2.0, 3.0, 5.0, 5.0]

    h = 0.1
    dr = oapply(d, euler_approx([r,s,r], h))
    dr2 = euler_approx(r2, h)
    dr3 = oapply(d, euler_approx([r,s,r]))
    dr4 = euler_approx(r2)
    @test eval_dynamics(dr, x0) == [2.3, 11.9, 4.1, 5.0]
    @test eval_dynamics(dr, x0) == eval_dynamics(dr2, x0)
    @test eval_dynamics(dr, x0) == eval_dynamics(dr3, x0, [h], 0)
    @test eval_dynamics(dr, x0) == eval_dynamics(dr4, x0, [h], 0)

    h = 0.25
    dr = oapply(d, euler_approx(xs, h))
    dr2 = euler_approx(r2, h)
    dr3 = oapply(d, euler_approx(xs))
    dr4 = euler_approx(r2)
    @test eval_dynamics(dr, x0) == eval_dynamics(dr2, x0)
    @test eval_dynamics(dr, x0) == eval_dynamics(dr3, x0, [h], 0)
    @test eval_dynamics(dr, x0) == eval_dynamics(dr4, x0, [h], 0)
  end

  @testset "Commutativity: substitute & oapply" begin
    d = UWD(4)
    add_box!(d, 1); add_box!(d, 2)
    add_junctions!(d, 3)
    set_junction!(d, [1,2,1])
    set_junction!(d, [1,2,3,3], outer = true)

    din = UWD(2)
    add_box!(din, 1); add_box!(din, 2)
    add_junctions!(din, 2)
    set_junction!(din, [1,1,2])
    set_junction!(din, [1,2], outer = true)

    dtot = ocompose(d, 2, din)

    let s = t
    r1 = oapply(dtot, [s,s,r])
    r2 = oapply(d, [s, oapply(din, [s,r])])

    @test nstates(r1) == nstates(r2)
    @test nports(r1) == nports(r2)
    @test portmap(r1) == portmap(r2)
    x0 = [2.0, 3.0, 5.0, 7.0, 11.0]
    @test eval_dynamics(r1, x0) == eval_dynamics(r2, x0)
    end
  end
end

@testset "DDE Problems" begin
  # delay differential equation with analytic solution via method of steps
  # dx/dt = -x(t-1) with history x(t) = 10 for t <= 0
  u0 = 10.0
  x0 = [u0]
  p = nothing
  hist(p, t) = [u0] # history

  alg = MethodOfSteps(Tsit5())

  # solve over 0<t<1; we have x(t-1)=10 over this interval therefore
  # dx/dt = -10; dx = -10dt; \int{dx} = \int{-10dt}; x + c1 = -10t + c2
  # x(t) = c - 10t; x(t) = 10 - 10t
  #
  # solve over 1<t<2; we have x(t-1)=10-10(t-1) over this interval therefore
  # dx/dt = -(10 - 10(t-1)); \int{dx} = \int{-(10 - 10(t-1))dt};
  # x + c1 = -10(t-1) + 10\int{(t-1) dt}; x + c1 = -10t + 10((t-1)^2/2) + c2;
  # x(t) = -10(t-1) + 5(t-1)^2
  ref_sol(t) = t < 1 ? 10.0 - 10.0*t : -10.0*(t-1.0) + 5.0*(t-1.0)^2
  ts = [0.5, 0.9, 1.5, 1.9]
  consistent(sol) = all(vcat(sol.(ts)...) .≈ ref_sol.(ts))

  @testset "Imperative diagram" begin
    df(u, h, p, t) = -h(p, t - 1.0)
    r = DelayResourceSharer{Float64}(1, df)
    prob = DDEProblem(r, [u0], hist, (0.0, 3.0), p)
    sol = solve(prob, alg,abstol=1e-12,reltol=1e-12)
    @test consistent(sol)

    # test oapply with UWD
    d = UWD(1)
    add_box!(d, 1)
    add_junctions!(d, 1)
    set_junction!(d, [1])
    set_junction!(d, [1], outer=true)

    r2 = oapply(d, [r])

    @test nstates(r) == nstates(r2)
    @test nports(r) == nports(r2)
    @test portmap(r) == portmap(r2)

    @test eval_dynamics(r, x0, hist, p, 0.0) == eval_dynamics(r2, x0, hist, p, 0.0)
    @test exposed_states(r, x0) == exposed_states(r2, x0)

    prob = DDEProblem(r2, [u0], hist, (0.0, 3.0), p)
    sol = solve(prob, alg,abstol=1e-12,reltol=1e-12)
    @test consistent(sol)
  end

  @testset "Relational diagram" begin
    df(u, h, p, t) = -0.5*h(p, t - 1.0)
    r = DelayResourceSharer{Float64}(1, df)
    d = @relation (x,) begin
      f(x)
      g(x)
    end
    r2 = oapply(d, r)

    prob = DDEProblem(r2, [u0], hist, (0.0, 3.0), p)
    sol = solve(prob, alg, abstol = 1e-12, reltol=1e-12)

    @test eval_dynamics(r2, x0, hist, p, 0.0) == [-10.0]
    @test consistent(sol)
  end
end