text:
using DoubleFloats
using LinearAlgebra
#abstract type AbstractExtDataAdv end

using SemiLagrangian:
    Advection,
    sizeall,
    AdvectionData,
    getdata,
    advection!,
    UniformMesh,
    AbstractInterpolation,
    Lagrange,
    Hermite,
    B_SplineLU,
    B_SplineFFT,
    interpolate!,
    gettabmod,
    CachePrecal,
    AbstractExtDataAdv,
    initcoef!,
    getalpha,
    OpTuple,
    TimeAlgorithm,
    NoTimeAlg,
    ABTimeAlg_ip,
    ABTimeAlg_init

import SemiLagrangian: initcoef!, getalpha
# """

#    exact(tf, mesh1, mesh2)

#    Julia function to compute exact solution "

# ```math
# \frac{d f}{dt} +  (y \frac{df}{delta1} - x \frac{df}{delta2}) = 0
# ```

# """
# function exact!(f, mesh1::UniformMesh{T}, mesh2::UniformMesh{T}, tf::T) where {T}
#     for (i, x) in enumerate(mesh1.points), (j, y) in enumerate(mesh2.points)
#         s, c = sincos(tf)
#         xn, yn = c * x - s * y, + s * x + c * y
#         f[i,j] = exp(-5*((xn)^2+(yn+T(6//5))^2))
#     end
# end

# function test_rotation(
#     sz::NTuple{2,Int},
#     interp_sp::AbstractInterpolation{T},
#     interp_v::AbstractInterpolation{T},
#     nbdt::Int
# ) where {T}
#     spmin, spmax, nsp =  T(-5), T(5),  sz[1]
#     vmin, vmax, nv = -T(12//2), T(9//2), sz[2]

#     mesh_sp = UniformMesh( spmin, spmax, nsp)
#     mesh_v = UniformMesh( vmin, vmax, nv)

#     dt = T(2big(pi)/nbdt)
#     adv = Advection1d((mesh_sp,), (mesh_v,), (interp_sp,), (interp_v,), dt; tab_fct=[tan, sin, tan])
# #    adv = Advection1d((mesh_sp,), (mesh_v,), (interp_sp,), (interp_v,), dt; tab_coef=[1], tab_fct=[identity])
#     sz = sizeall(adv)
#     tabref = zeros(T,sz)
#     exact!(tabref, mesh_sp, mesh_v, T(0))

#     pvar = getrotationvar(adv)

#     advdata = Advection1dData(adv, tabref, pvar)

#     diffmax = 0
#     data = getdata(advdata)
#     for ind=1:nbdt
#         while advection!(advdata) end
#         exact!(tabref, mesh_sp, mesh_v, dt*ind)
#         diff = norm(data .- tabref, Inf)
#         diffmax = max(diffmax, diff)
#         @show ind, diff
#     end
#     println("test_rotation sz=$sz interp=$interp_sp, $interp_v nbdt=$nbdt diffmax=$diffmax")
#     return diffmax

# end

function test_swirling(
    sz::NTuple{2,Int},
    ::T,
    interp::AbstractVector{I},
    nbdt::Int,
) where {T, I <: AbstractInterpolation{T}}
    spmin, spmax, nsp = T(0), T(1), sz[ 1]
    vmin, vmax, nv = T(0), T(1), sz[2]

    mesh_sp = UniformMesh(spmin, spmax, nsp)
    mesh_v = UniformMesh(vmin, vmax, nv)

    t_max = T(1.5)
    dt = t_max / nbdt

    # dec1 = -2tan(dt/2) / step(mesh_sp) * mesh_v.points
    # dec2 = 2tan(dt/2) / step(mesh_v) * mesh_sp.points
    dec1 = zeros(T, sz)
    dec2 = zeros(T, sz)
    tgdt = 2tan(dt / 2)
    coef = 1 / (1 + tgdt^2 / 4)
    tabref = zeros(T, sz)
    for i = 1:sz[1], j = 1:sz[2]
        x = mesh_sp.points[i]
        y = mesh_v.points[j]
        v = [
            -cos(T(pi) * x)^2 * sin(2T(pi) * y) / step(mesh_sp),
            cos(T(pi) * y)^2 * sin(2T(pi) * x) / step(mesh_v),
        ]
        dec1[i, j], dec2[i, j] = v
        tabref[i, j] = ((1 - x)^2 + (1 - y)^2 < 0.8) ? 1 : 0
    end

    @show minimum(dec1), maximum(dec1)
    @show minimum(dec2), maximum(dec2)

    #    @show dec1, dec2

    data = zeros(T, sz)

    copyto!(data, tabref)

    diff = norm(data .- tabref)

    @show 0, diff

    tabmod=gettabmod.(sz)
    cache=CachePrecal(interp, one(T))

    buf = zeros(T, sz)
    for ind = 1:nbdt
        coef = dt * cos(T(pi) * (ind-1) / nbdt)
        interpolate!(buf, data, x -> (coef * dec1[x], coef * dec2[x]), interp; tabmod=tabmod, cache=cache)
        copyto!(data, buf)
        diff = norm(data .- tabref)
        @show ind, diff
    end
    println("test_swirling sz=$sz interp=$interp nbdt=$nbdt diff=$diff")
    return diff
end
mutable struct Swirling{T} <: AbstractExtDataAdv
    ref::Array{OpTuple{2,T},2}
    function Swirling(dec::Array{OpTuple{2,T},2}) where {T}
        return new{T}(copy(dec))
    end

end
@inbounds function initcoef!(ext::Swirling{T}, advd::AdvectionData{T}) where{T}
    coef = advd.adv.dt_base * cos(T(pi) * advd.time_cur/T(1.5))
    sz = sizeall(advd.adv)
    if ismissing(advd.bufcur)
        advd.bufcur = zeros(OpTuple{2,T}, sz)
    end
    advd.bufcur .= coef * ext.ref
    missing
end

function test_swirling_adv(
    sz::NTuple{2,Int},
    ::T,
    interp::AbstractVector{I},
    nbdt::Int;
    timealg::TimeAlgorithm=NoTimeAlg,
    ordalg=0
) where {T, I <: AbstractInterpolation{T}}
    spmin, spmax, nsp = T(0), T(1), sz[ 1]
    vmin, vmax, nv = T(0), T(1), sz[2]

    mesh_sp = UniformMesh(spmin, spmax, nsp)
    mesh_v = UniformMesh(vmin, vmax, nv)

    t_max = T(1.5)
    dt = t_max / nbdt

    # dec1 = -2tan(dt/2) / step(mesh_sp) * mesh_v.points
    # dec2 = 2tan(dt/2) / step(mesh_v) * mesh_sp.points
    dec = zeros(OpTuple{2,T}, sz)
    # tgdt = 2tan(dt / 2)
    # coef = 1 / (1 + tgdt^2 / 4)
    tabref = zeros(T, sz)
    for i = 1:sz[1], j = 1:sz[2]
        x = mesh_sp.points[i]
        y = mesh_v.points[j]
        v = (
            -cos(T(pi) * x)^2 * sin(2T(pi) * y) / step(mesh_sp),
            cos(T(pi) * y)^2 * sin(2T(pi) * x) / step(mesh_v),
        )
        dec[i, j] = OpTuple(v)
        tabref[i, j] = exp(-10*((T(1)/4 - x)^2 + (T(1)/2 - y)^2))
    end

    adv = Advection((mesh_sp,mesh_v),interp,dt,[([1,2], 2, 1, false),],tab_coef=[dt,], timealg=timealg, ordalg=ordalg)

    initdatas = missing
    if timealg == ABTimeAlg_init
        initdatas = [copy(tabref) for i=1:(3ordalg-1)]
    end

    advd = AdvectionData(adv, tabref, Swirling(dec),initdatas=initdatas)

    if timealg == ABTimeAlg_init
        advd.time_cur -=  length(initdatas)*dt
    end

    f(v,x) = v.v[x]
    @show extrema(f.(dec,(1,)))
    @show extrema(f.(dec,(2,)))

    diff = norm(getdata(advd) .- tabref)

    @show 0, diff

    for ind = 1:nbdt
        while advection!(advd)
        end
        diff = norm(getdata(advd) .- tabref)
        @show ind, diff
    end
    println("test_swirling_adv sz=$sz interp=$interp nbdt=$nbdt diff=$diff")
    return diff
end

@testset "test swirling" begin
    T = Float64
    @time @test test_swirling((100, 100), one(T), [Lagrange(9,T),Lagrange(9,T)] , 50) < 15
end
@testset "test swirling_adv" begin
    T = Float64
    @time @test test_swirling_adv((100, 100), one(T), [Lagrange(9,T),Lagrange(9,T)] , 50) < 5
    @time @test test_swirling_adv((100, 100), one(T), [B_SplineLU(9,100,T),B_SplineLU(9,100,T)] , 50) < 5
    @time @test test_swirling_adv((100, 100), one(T), [Hermite(9,T),Hermite(9,T)] , 50) < 5

    @time @test test_swirling_adv((100, 100), one(T), [Lagrange(9,T),Lagrange(9,T)] , 50, timealg=ABTimeAlg_init, ordalg=4) < 2
    @time @test test_swirling_adv((100, 100), one(T), [Hermite(9,T),Hermite(9,T)] , 50, timealg=ABTimeAlg_init, ordalg=4) < 2
end


text:
models = (
    GaussianProcessRegressor,
    GaussianProcessClassifier,
    )

fparams = (
    GaussianProcessRegressor=(:X_train, :y_train, :kernel, :L, :alpha, :log_marginal_likelihood_value),
    GaussianProcessClassifier=(:kernel, :log_marginal_likelihood_value, :classes, :n_classes),
    )

@testset "Fit/Predict" begin
    X, y = simple_binaryclf()
    m = GaussianProcessClassifier()
    acc, f = test_clf(m, X, y)
    @test acc >= 0.75
    @test keys(MB.fitted_params(m, f)) == fparams.GaussianProcessClassifier

    X, y, ls = simple_regression()
    m = GaussianProcessRegressor()
    r, f = test_regression(m, X, y, ls)
    @test r < 0.1 # super over-fitting
end

@testset "Docstrings" begin
    for mod in models
        m = mod()
        @test !isempty(MB.docstring(m))
    end
end

ext: jl


text:
@testset "Nodes" begin
    @test typeof(Node(:a)) <: Node
end

@testset "Edges" begin
    @test Edge(Node(:a), Node(:b)) == Edge(Node(:a), Node(:b)) == DirectedEdge(SimpleNode(:a), SimpleNode(:b))
end

# define modifier
struct Weight <: EdgeModifier end

@testset "ModifedNode" begin
    @test ModifyingNode(Node(:a), Weight()) == ModifyingNode(Node(:a), Weight())
end

@testset "ModifedEdges" begin
    # ModifyingNode leads to ModifiedEdge with simple Node
    @test Edge(Node(:a), ModifyingNode(Node(:b), Weight())) == ModifiedEdge(Edge(Node(:a), Node(:b)), Weight())
    # check if the same works with vector
    @test [ModifiedEdge(Edge(Node(:a), Node(:b)), Weight()), Edge(Node(:a), Node(:c))] == Edge(Node(:a), [Node(:b) * Weight() Node(:c)])
end


text:
using GAIO
using StaticArrays
using Test

@testset "exported functionality" begin
    partition = TreePartition(Box(SVector(0.0, 0.0, 0.0, 0.0), SVector(1.0, 1.0, 1.0, 1.0)))
    @testset "basics" begin
        @test depth(partition) == 0
        @test dimension(partition) == 4
    end
    @testset "subdivision" begin
        keys_to_subdivide = [(0, 1), (1, 1), (1, 2), (2, 2), (2, 4)]
        for key in keys_to_subdivide
            subdivide!(partition, key)
        end
        @test depth(partition) == 3
        @test dimension(partition) == 4
    end
    @testset "domain with zero radius" begin
        center = SVector(0.0, 0.0)
        radius = SVector(1.0, 0.0)
        @test_throws DomainError Box(center, radius)
    end
end


text:
function PolyShape(body::Body, verts::Array{Vect}, radius::Real)
    shape_ptr = ccall(dlsym(libchipmunk, :cpPolyShapeAlloc), Ptr{Void}, ())
    PolyShape(ccall(dlsym(libchipmunk, :cpPolyShapeInitRaw), Ptr{Void}, (Ptr{Void}, Ptr{Void}, Int32, Ptr{Vect}, Cdouble,),
                                                                        shape_ptr, body.ptr, length(verts), pointer(verts), radius))
end

function BoxShape(body::Body, width::Real, height::Real, radius::Real)
    PolyShape(ccall(dlsym(libchipmunk, :cpBoxShapeNew), Ptr{Void}, (Ptr{Void}, Cdouble, Cdouble, Cdouble,), body.ptr, width, height, radius))
end

function get_count(shape::PolyShape)
    ccall(dlsym(libchipmunk, :cpPolyShapeGetCount), Int32, (Ptr{Void},), shape.ptr)
end

function get_vert(shape::PolyShape, index::Integer)
    ccall(dlsym(libchipmunk, :cpPolyShapeGetVert), Vect, (Ptr{Void}, Int32,), shape.ptr, index)
end

function get_radius(shape::PolyShape)
    ccall(dlsym(libchipmunk, :cpPolyShapeGetRadius), Cdouble, (Ptr{Void},), shape.ptr)
end

export PolyShape, get_count, get_vert, get_radius


text:
using Test
using LogicCircuits

@testset "Unbatch test" begin # wrapping all the statements into a function also trigers signal (4) due to `XDataset`
    for name in ["msnbc", "jester", "plants", "ad"]
        one_batch_set = dataset(twenty_datasets(name); do_shuffle=false, batch_size=-1)
        batches_set = dataset(twenty_datasets(name); do_shuffle=false, batch_size=100)
        unbatched_set = XDataset(unbatch(batches_set)...)
        for f in [train, valid, test]
            # unweighted data
            @test feature_matrix(f(unbatched_set)) == feature_matrix(f(one_batch_set))

            # weighted data
            w = rand(num_examples(f(one_batch_set)))
            w_one_batch = WXData(f(one_batch_set), w)

            w_batches = batch(w_one_batch, 100)
            w_unbatched = unbatch(w_batches)

            @test w_one_batch.xd.x == w_unbatched.xd.x
            @test w_one_batch.w == w_unbatched.w
        end
    end
end


text:
function set_velocity_solution!(body::Body)
    state = body.state
    state.vsol[1] = state.v15
    state.vsol[2] = state.v15
    state.ϕsol[1] = state.ϕ15
    state.ϕsol[2] = state.ϕ15
end

function set_previous_configuration!(body::Body{T}, timestep::T) where {T}
    x2, v15, q2, ϕ15 = initial_configuration_velocity(body.state)
    body.state.x1 = next_position(x2, -v15, timestep)
    body.state.q1 = next_orientation(q2, -ϕ15, timestep)
end

function initialize_state!(body::Body{T}, timestep) where T
    set_previous_configuration!(body, timestep)
    state = body.state
    state.F2 = szeros(T,3)
    state.τ2 = szeros(T,3)
end

function update_state!(body::Body{T}, timestep) where T
    state = body.state

    state.x1 = state.x2
    state.q1 = state.q2

    state.v15 = state.vsol[2]
    state.ϕ15 = state.ϕsol[2]

    state.x2 = next_position(state.x2, state.vsol[2], timestep)
    state.q2 = next_orientation(state.q2, state.ϕsol[2], timestep)

    state.F2 = szeros(T,3)
    state.τ2 = szeros(T,3)
end

# maximal
function set_maximal_configurations!(body::Body;
    x::AbstractVector=SA[0.0; 0.0; 0.0],
    q::Quaternion=one(Quaternion))

    body.state.x2 = x
    body.state.q2 = q

    return body.state.x2, body.state.q2
end

function set_maximal_velocities!(body::Body;
    v::AbstractVector=SA[0.0; 0.0; 0.0],
    ω::AbstractVector=SA[0.0; 0.0; 0.0])

    body.state.v15 = v
    body.state.ϕ15 = ω

    return body.state.v15, body.state.ϕ15
end

function set_maximal_configurations!(pbody::Node, cbody::Body;
    parent_vertex::AbstractVector=SA[0.0; 0.0; 0.0],
    child_vertex::AbstractVector=SA[0.0; 0.0; 0.0],
    Δx::AbstractVector=SA[0.0; 0.0; 0.0],
    Δq::Quaternion=one(Quaternion))

    q1 = pbody.state.q2
    q2 = pbody.state.q2 * Δq
    x2 = pbody.state.x2 + vector_rotate(parent_vertex + Δx, q1) - vector_rotate(child_vertex, q2)

    return set_maximal_configurations!(cbody; x=x2, q=q2)
end

function set_maximal_velocities!(pbody::Node, cbody::Body;
    parent_vertex::AbstractVector=SA[0.0; 0.0; 0.0],
    child_vertex::AbstractVector=SA[0.0; 0.0; 0.0],
    Δv::AbstractVector=SA[0.0; 0.0; 0.0],
    Δω::AbstractVector=SA[0.0; 0.0; 0.0])

    x1 = pbody.state.x2
    v1 = pbody.state.v15
    q1 = pbody.state.q2
    ω1 = pbody.state.ϕ15

    x2 = cbody.state.x2
    q2 = cbody.state.q2

    ω2 = vector_rotate(Δω + ω1, inv(q2) * q1)
    ω1w = vector_rotate(ω1, q1)
    ω2w = vector_rotate(ω2, q2)
    Δvw = vector_rotate(Δv, q1)
    cApB_w = (x2 + vector_rotate(child_vertex, q2)) - x1
    pBcB_w = - vector_rotate(child_vertex, q2)
    v2 = copy(v1)
    v2 += skew(ω1w) * cApB_w
    v2 += skew(ω2w) * pBcB_w
    v2 += Δvw

    return set_maximal_velocities!(cbody; v = v2, ω = ω2)
end

# inputs
function set_input!(body::Body;
    F::AbstractVector=SA[0.0; 0.0; 0.0],
    τ::AbstractVector=SA[0.0; 0.0; 0.0],
    p::AbstractVector=SA[0.0; 0.0; 0.0])

    τ += cross(p, F) # in local coordinates
    set_input!(body.state, vector_rotate(F, body.state.q2), τ)
end


text:
module QRupdate

using LinearAlgebra

export qraddcol, qraddrow, qrdelcol, csne

"""
Add a column to a QR factorization without using Q.

`R = qraddcol(A,R,v)`  adds the m-vector `v` to the QR factorization of the
m-by-n matrix `A` without using `Q`. On entry, `R` is a dense n-by-n upper
triangular matrix such that `R'*R = A'*A`.

`R = qraddcol(A,R,v,β)` is similar to the above, except that the
routine updates the QR factorization of
```
[A; β I],   and   R'*R = (A'*A + β^2*I).
```

`A` should have fewer columns than rows.  `A` and `v` may be sparse or
dense.  On exit, `R` is the (n+1)-by-(n+1) factor corresponding to

```
  Anew = [A        V ],    Rnew = [R   u  ],   Rnew'Rnew = Anew'Anew.
         [beta*I     ]            [  gamma]
         [       beta]
```

The new column `v` is assumed to be nonzero.
If `A` has no columns yet, input `A = []`, `R = []`.
"""
function qraddcol(A::AbstractMatrix{T}, Rin::AbstractMatrix{T}, a::Vector{T}, β::T = zero(T)) where {T}

    m, n = size(A)
    anorm  = norm(a)
    anorm2 = anorm^2
    β2  = β^2
    if β != 0
        anorm2 = anorm2 + β2
        anorm  = sqrt(anorm2)
    end

    if n == 0
        return reshape([anorm], 1, 1)
    end

    R = UpperTriangular(Rin)

    c      = A'a           # = [A' β*I 0]*[a; 0; β]
    u      = R'\c
    unorm2 = norm(u)^2
    d2     = anorm2 - unorm2

    if d2 > anorm2 #DISABLE 0.01*anorm2     # Cheap case: γ is not too small
        γ = sqrt(d2)
    else
        z = R\u          # First approximate solution to min ||Az - a||
        r = a - A*z
        c = A'r
        if β != 0
            c = c - β2*z
        end
        du = R'\c
        dz = R\du
        z  += dz          # Refine z
      # u  = R*z          # Original:     Bjork's version.
        u  += du          # Modification: Refine u
        r  = a - A*z
        γ = norm(r)       # Safe computation (we know gamma >= 0).
        if !iszero(β)
            γ = sqrt(γ^2 + β2*norm(z)^2 + β2)
        end
    end

    # This seems to be faster than concatenation, ie:
    # [ Rin         u
    #   zeros(1,n)  γ ]
    Rout = zeros(T, n+1, n+1)
    Rout[1:n,1:n] .= R
    Rout[1:n,n+1] .= u
    Rout[n+1,n+1] = γ

    return Rout
end

"""
Add a row and update a Q-less QR factorization.

`qraddrow!(R, a)` returns the triangular part of a QR factorization of `[A; a]`, where `A = QR` for some `Q`.  The argument `a` should be a row
vector.
"""
function qraddrow(R::AbstractMatrix{T}, a::AbstractMatrix{T}) where {T}

    n = size(R,1)
    @inbounds @simd for k in 1:n
        G, r = givens( R[k,k], a[k], 1, 2 )
        B = G * [ reshape(R[k,k:n], 1, n-k+1)
                  reshape(a[:,k:n], 1, n-k+1) ]
        R[k,k:n] = B[1,:]
        a[  k:n] = B[2,:]
    end
    return R
end

"""
Delete the k-th column and update a Q-less QR factorization.

`R = qrdelcol(R,k)` deletes the k-th column of the upper-triangular
`R` and restores it to upper-triangular form.  On input, `R` is an n x
n upper-triangular matrix.  On output, `R` is an n-1 x n-1 upper
triangle.

    18 Jun 2007: First version of QRdelcol.m.
                 Michael Friedlander (mpf@cs.ubc.ca) and
                 Michael Saunders (saunders@stanford.edu)
                 To allow for R being sparse,
                 we eliminate the k-th row of the usual
                 Hessenberg matrix rather than its subdiagonals,
                 as in Reid's Bartel-Golub LU update and also
                 the Forrest-Tomlin update.
                 (But Matlab will still be pretty inefficient.)
    18 Jun 2007: R is now the exact size on entry and exit.
    30 Dec 2015: Translate to Julia.
"""
function qrdelcol(R::AbstractMatrix{T}, k::Int) where {T}

    m = size(R,1)
    R = R[:,1:m .!= k]          # Delete the k-th column
    n = size(R,2)               # Should have m=n+1

    for j in k:n                # Forward sweep to reduce k-th row to zeros
        G, y = givens(R[j+1,j], R[k,j], 1, 2)
        R[j+1,j] = y
        if j<n && G.s != 0
            @inbounds @simd for i in j+1:n
                tmp = G.c*R[j+1,i] + G.s*R[k,i]
                R[k,i] = G.c*R[k,i] - conj(G.s)*R[j+1,i]
                R[j+1,i] = tmp
            end
        end
    end
    R = R[1:m .!= k, :]         # Delete the k-th row
    return R
end

"""
    x, r = csne(R, A, b)

solves the least-squares problem

    minimize  ||r||_2,  r := b - A*x

using the corrected semi-normal equation approach described by Bjork (1987). Assumes that `R` is upper triangular.
"""
function csne(Rin::AbstractMatrix{T}, A::AbstractMatrix{T}, b::Vector{T}) where {T}

    R = UpperTriangular(Rin)
    q = A'b
    x = R' \ q

    bnorm2 = sum(abs2, b)
    xnorm2 = sum(abs2, x)
    d2 = bnorm2 - xnorm2

    x = R \ x

    # Apply one step of iterative refinement.
    r = b - A*x
    q = A'r
    dx = R' \ q
    dx = R  \ dx
    x += dx
    r = b - A*x
    return (x, r)
end

end # module


text:
module help

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
function my_throw(e::Exception, location_msg::String)
    println("ERROR at ", location_msg)
    throw(e)
end

#=
type exception_recursion <: Exception
    prev_e::Exception
    location_msg::String
    function exception_recursion(prev_e::Exception, location_msg::String)
        return new(prev_e, location_msg)
    end
end

function print(e::exception_recursion)
    println("ERROR at ", e.location_msg)
    print(e.prev_e)
end

function print(e::Exception)
    print(e.msg)
end
=#


text:
immutable Epanechnikov <: ContinuousUnivariateDistribution
    μ::Float64
    σ::Float64

    Epanechnikov(μ::Real, σ::Real) = (@check_args(Epanechnikov, σ > zero(σ)); new(μ, σ))
    Epanechnikov(μ::Real) = new(μ, 1.0)
    Epanechnikov() = new(0.0, 1.0)
end

@distr_support Epanechnikov d.μ - d.σ d.μ + d.σ

## Parameters

location(d::Epanechnikov) = d.μ
scale(d::Epanechnikov) = d.σ
params(d::Epanechnikov) = (d.μ, d.σ)

## Properties
mean(d::Epanechnikov) = d.μ
median(d::Epanechnikov) = d.μ
mode(d::Epanechnikov) = d.μ

var(d::Epanechnikov) = d.σ^2 / 5
skewness(d::Epanechnikov) = 0.0
kurtosis(d::Epanechnikov) = -2.914285714285714  # 3/35-3

## Functions
function pdf(d::Epanechnikov, x::Float64)
    u = abs(x - d.μ) / d.σ
    u >= 1 ? 0.0 : 0.75 * (1 - u^2) / d.σ
end

function cdf(d::Epanechnikov, x::Float64)
    u = (x - d.μ) / d.σ
    u <= -1 ? 0.0 :
    u >= 1 ? 1.0 :
    0.5 + u * (0.75 - 0.25 * u^2)
end

function ccdf(d::Epanechnikov, x::Float64)
    u = (d.μ - x) / d.σ
    u <= -1 ? 1.0 :
    u >= 1 ? 0.0 :
    0.5 + u * (0.75 - 0.25 * u^2)
end

@quantile_newton Epanechnikov

function mgf(d::Epanechnikov, t::Float64)
    a = d.σ * t
    a == 0 ? 1.0 :
    3.0 * exp(d.μ * t) * (cosh(a) - sinh(a) / a) / a^2
end

function cf(d::Epanechnikov, t::Float64)
    a = d.σ * t
    a == 0 ? 1.0+0.0im :
    -3.0 * exp(im * d.μ * t) * (cos(a) - sin(a) / a) / a^2
end


text:
include("euler/Julia/grid.jl")
include("euler/Julia/math.jl")
using .Grid: adjacent_product
using .Math: fast_max

function compute(n::Int, grid::String)::Int
    grid = [parse.(Int, row) for row in split.(split(grid, "\n"), " ")]
    width, height = length(grid[1]), length(grid)
    result = 0
    for y ∈ 1:height
        for x ∈ 1:width
            if x + n - 1 ≤ width result = fast_max(result, adjacent_product(grid, n, (x, y), (1, 0))) end
            if y + n - 1 ≤ height result = fast_max(result, adjacent_product(grid, n, (x, y), (0, 1)))
                if x + n - 1 ≤ width result = fast_max(result, adjacent_product(grid, n, (x, y), (1, 1))) end
                if x - n ≥ 0 result = fast_max(result, adjacent_product(grid, n, (x, y), (-1, 1))) end
            end
        end
    end
    return result
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
abstract LocalNumSys

"""
# Stores local DOF numbering systems

Each local degree of freedom (DOF) has an associated index.
Often, there is an advantage to separating the degrees of
freedom into boundary dof (nonzero support on the boundaries
of the element) and interior degrees of freedom (zero on
boundaries). Depending on how the basis DOF are initially
defined, two local numbering systems are necessary:

 * Sequential numbering, most natural when computing
 * Boundary/Interior numbering where boundary modes are numbered first and then interior modes.

"""
immutable LocalNumSys1d <: LocalNumSys
    "Number of boundary modes"
    nb::Int
    "Number of interior modes"
    ni::Int
    "Indicesof boundary modes"
    bndry::Vector{Int}
    "Indices of interior modes"
    intrr::Vector{Int}
end
LocalNumSys1d(b::Vector{Int}, i::Vector{Int}) = LocalNumSys1d(length(b), length(i), b, i)

isbndryint(n::LocalNumSys1d) = n.ni > 0

"Total number of modes"
nmodes(n::LocalNumSys1d) = n.nb + n.ni

"Number of boundary modes"
nbndry(n::LocalNumSys1d) = n.nb

"Number of interior modes"
ninterior(n::LocalNumSys1d) = n.ni

"Indices of boundary modes"
bndidx(n::LocalNumSys1d) = n.bndry

"Indices of interior modes"
intidx(n::LocalNumSys1d) = n.intrr

"Convert vector from sequential numbering to boundary/interior"
function seq2bi!{T}(lnum::LocalNumSys1d, x::AbstractVector{T}, y::AbstractVector{T})

    for i = 1:lnum.nb
        y[i] = x[lnum.bndry[i]]
    end
    for i = 1:lnum.ni
        y[i+lnum.nb] = x[lnum.intrr[i]]
    end
    y
end

"Convert vector from sequential numbering to boundary/interior"
seq2bi{T}(lnum::LocalNumSys1d, x::AbstractVector{T}) = seq2bi!(lnum, x, zeros(T, nmodes(lnum)))

"Convert vector from boundary/interior numbering to sequential numbering"
function bi2seq!{T}(lnum::LocalNumSys1d, x::AbstractVector{T}, y::AbstractVector{T})
    for i = 1:lnum.nb
        y[lnum.bndry[i]] = x[i]
    end
    for i = 1:lnum.ni
        y[lnum.intrr[i]] = x[i+lnum.nb]
    end
    y
end

"Convert vector from boundary/interior numbering to sequential numbering"
bi2seq{T}(lnum::LocalNumSys1d, x::AbstractVector{T}) = bi2seq!(lnum, x, zeros(T, nmodes(lnum)))

"Convert vector from boundary/interior numbering to sequential numbering"
function bi2seq!{T}(lnum::LocalNumSys1d, xb::AbstractVector{T}, xi::AbstractVector{T},
                    y::AbstractVector{T})
    for i = 1:lnum.nb
        y[lnum.bndry[i]] = xb[i]
    end
    for i = 1:lnum.ni
        y[lnum.intrr[i]] = xi[i]
    end
    y
end

"Convert vector from boundary/interior numbering to sequential numbering"
bi2seq{T}(lnum::LocalNumSys1d, xb::AbstractVector{T}, xi::AbstractVector{T}) =
    bi2seq!(lnum, xb, xi, zeros(T, nmodes(lnum)))

"Extract boundary modes from modes numbered sequentially"
function seq2b!{T}(lnum::LocalNumSys1d, x::AbstractVector{T}, y::AbstractVector{T})
    for i = 1:lnum.nb
        y[i] = x[lnum.bndry[i]]
    end
    y
end

"Extract boundary modes from modes numbered sequentially"
seq2b{T}(lnum::LocalNumSys1d, x::AbstractVector{T}) = seq2b!(lnum, x, zeros(T, nbndry(lnum)))

"Extract interior modes from modes numbered sequentially"
function seq2i!{T}(lnum::LocalNumSys1d, x::AbstractVector{T}, y::AbstractVector{T})
    for i = 1:lnum.ni
        y[i] = x[lnum.intrr[i]]
    end
    y
end

"Extract interior modes from modes numbered sequentially"
seq2i{T}(lnum::LocalNumSys1d, x::AbstractVector{T}) = seq2i!(lnum, x, zeros(T, ninterior(lnum)))


text:
using Revise
using ZonalFlow
using BenchmarkTools

lx = Float64(2π);
ly = Float64(2π);
nx = 6;
ny = 6;

dt = 0.005;
t_end = 10.0;
savefreq = 20;

β = 10.0 # β = 10.0
μ = 0.01 # μ = 0.01
ν = 0.0
ν₄ = 0.0 # ν₄ = 1.0

kf = 3
dk = 1
ε = 0.01

Λ = 0

@benchmark acoeffs($ny)
@benchmark bcoeffs($lx,$ly,$nx,$ny,$β,$μ,$ν,$ν₄)
@benchmark ccoeffs($lx,$ly,$nx,$ny,$0)
@benchmark fcoeffs($nx,$ny,$0,$kf,$dk,$ε)
@benchmark ic_rand($nx,$ny,1e-3)

A = acoeffs(ny)
B = bcoeffs(lx,ly,nx,ny,β,μ,ν,ν₄)
Cp,Cm = ccoeffs(lx,ly,nx,ny,0)
F = fcoeffs(nx,ny,0,kf,dk,ε).x[2]

u0 = ic_rand(nx,ny,1e-3)
u0 = ic_cumulants(nx,ny,u0)
du = similar(u0)
u = fill!(similar(du),0)

dx = fill!(similar(u0.x[1]),0)
dy = fill!(similar(u0.x[2]),0)
temp = fill!(similar(u0.x[2]),0)

p = [nx,ny,A,B,Cp,Cm,F,dx,dy,temp]
t = (0.0,1.0)

@time ce2_eqs!(du,u,p,t)
@code_warntype ce2_eqs!(du,u,p,t)
@benchmark ce2_eqs!($du,$u,$p,$t)
@code_lowered ce2_eqs!(du,u,p,t)

@code_warntype unit_eqs!(du,u,p,t)
@benchmark unit_eqs!($du,$u,$p,$t)

@benchmark ODEProblem($ce2_eqs!,$u0,$t,$p)
@code_warntype ODEProblem(ce2_eqs!,u0,t,p)
prob = ODEProblem(ce2_eqs!,u0,t,p)
@benchmark solve($prob,$Heun(),save_start=false,save_everystep=false)

@code_warntype solve(prob,Heun(),save_everystep=false)

@code_warntype fcoeffs(nx,ny,0,kf,dk,ε,isotropic=true)
@code_llvm fcoeffs(nx,ny,0,kf,dk,ε,isotropic=true)

# Cp,Cm = ccoeffs(nx,ny)
F = fcoeffs(nx,ny,0,kf,dk,ε,isotropic=isotropic).x[2]


text:
# Kochenderfer, Mykel J.. Algorithms for Optimization (The MIT Press)

function branin(x; a=1, b=5.1/(4π^2), c=5/π, r=6, s=10, t=1/(8π))
    return a*(x[2]-b*x[1]^2+c*x[1]-r)^2 + s*(1-t)*cos(x[1]) + s
end

function ackley(x, a=20, b=0.2, c=2π)
    d = length(x)
    return -a*exp(-b*sqrt(sum(x.^2)/d)) - exp(sum(cos.(c*xi) for xi in x)/d) + a + exp(1)
end

paraboloid(x; a=1, b=1) = x[1]^2/a^2 + x[2]^2/b^2


text:
using Test

m = test_page_model()
include("../src/components/CO2cycle.jl")

add_comp!(m, co2cycle)

set_param!(m, :co2cycle, :e_globalCO2emissions, readpagedata(m, "test/validationdata/e_globalCO2emissions.csv"))
set_param!(m, :co2cycle, :y_year,[2009.,2010.,2020.,2030.,2040.,2050.,2075.,2100.,2150.,2200.])#real values
set_param!(m, :co2cycle, :y_year_0,2008.)#real value
set_param!(m, :co2cycle, :rt_g_globaltemperature, readpagedata(m, "test/validationdata/rt_g_globaltemperature.csv"))

p=load_parameters(m)
set_leftover_params!(m,p)

run(m)

pop=m[:co2cycle,  :c_CO2concentration]

pop_compare=readpagedata(m, "test/validationdata/c_co2concentration.csv")

@test pop ≈ pop_compare rtol=1e-4


text:
mutable struct OrgCache
    doc
    headings::Union{Nothing, Vector{Heading}}
    elements::Union{Nothing, Vector{Element}}
    components::Union{Nothing, Vector{Component}}
    keywords::Union{Nothing, Dict{SubString{String}, Vector{SubString{String}}}}
    macros::Union{Nothing, Dict{SubString{String}, Function}}
    footnotes::Union{Nothing, Dict{Union{SubString{String}, FootnoteReference},
                                   Tuple{Int, Union{FootnoteDefinition, FootnoteReference}}}}
end

OrgCache(doc) = OrgCache(doc, nothing, nothing, nothing,
                         nothing, nothing, nothing)
OrgCache() = OrgCache(nothing)


text:
using DiffEqSensitivity, Flux, OrdinaryDiffEq, Test # , Plots

p = [1.5 1.0;3.0 1.0]
function lotka_volterra(du,u,p,t)
  du[1] = p[1,1]*u[1] - p[1,2]*u[1]*u[2]
  du[2] = -p[2,1]*u[2] + p[2,2]*u[1]*u[2]
end

u0 = [1.0,1.0]
tspan = (0.0,10.0)

prob = ODEProblem(lotka_volterra,u0,tspan,p)
sol = solve(prob,Tsit5())

# plot(sol)

p = [2.2 1.0;2.0 0.4] # Tweaked Initial Parameter Array
ps = Flux.params(p)

function predict_adjoint() # Our 1-layer neural network
  Array(solve(prob,Tsit5(),p=p,saveat=0.0:0.1:10.0))
end

loss_adjoint() = sum(abs2,x-1 for x in predict_adjoint())

data = Iterators.repeated((), 100)
opt = ADAM(0.1)
cb = function () #callback function to observe training
  display(loss_adjoint())
end

predict_adjoint()

# Display the ODE with the initial parameter values.
cb()
Flux.train!(loss_adjoint, ps, data, opt, cb = cb)

@test loss_adjoint() < 1


text:
#export simulate_readcoverage,Para,increment!,hist_add!,vec_hist_add!,simulate_and_get_coveage_hist,plot_readcoverage_simulation

using OffsetArrays
using ProgressMeter
using PyPlot
using LsqFit

Base.@kwdef struct Para
    genome_length::Int = 10_000_000
    read_length::Int = 100
    read_number::Int
end

function increment!(h::Vector{T}, i::Int) where T <: Number
    if i > length(h)
        lh = length(h)
        resize!(h, i)
        for k in (lh + 1) : i
            @inbounds h[k] = zero(T)
        end
    end
    @inbounds h[i] += one(T)
    h
end

function hist_add!(h1::Vector{T}, h2::Vector{T}) where T <: Number
    if length(h1) <= length(h2)
        h = h2
        for i in 1:length(h1)
            @inbounds h[i] += h1[i]
        end
    else
        h = h1
        for i in 1:length(h2)
            @inbounds h[i] += h2[i]
        end
    end
    return h
end

function vec_hist_add!(hs1, hs2)
    @assert length(hs1) == length(hs2)
    map(hist_add!, hs1, hs2)
end

function genome_simulation(p::Para)
    # simulation
    g = fill(zero(Int32), p.genome_length)
    for i in 1:p.read_number
        pos = rand(1:(p.genome_length - p.read_length))
        g[pos] += 1
        g[pos + p.read_length] -= 1
    end
    return g
end

function get_coverage_hist(g,genome_length)
    max_coverage_for_hist = 2
    hists = OffsetArray(map(x->Int[], 1:max_coverage_for_hist+1), 0:max_coverage_for_hist)
    # compute coverages
    g0 = 0 + g[1]
    p0 = 1
    for pos in 2:genome_length
        if g[pos] != 0
            g0 <= max_coverage_for_hist && increment!(hists[g0], pos - p0)
            g0 += g[pos]
            p0 = pos
        end
    end
    hists
end

function simulate_readcoverage(n,genome_length,read_length,no_reads)
    paras = fill(Para(genome_length = genome_length,read_length=read_length,read_number=no_reads), n);
    res1= @showprogress map(genome_simulation,paras)
    res = @showprogress map(get_coverage_hist, res1,genome_length);
    h = reduce(vec_hist_add!, res) ./ n
    return h
end


text:
@testset "CFR" begin
    include("best_response.jl")
    include("nash_conv.jl")
    include("mccfr.jl")
end


text:
using VoronoiCells
using Base.Test

# Corners of the unit square
x = [1.0 ; 0.0 ; 0.0 ; 1.0]
y = [1.0 ; 1.0 ; 0.0 ; 0.0]

@test_approx_eq density(x, y) sqrt(2)/2


text:
using BinaryBuilder, Pkg

# Beefalo!
name = "beefalo"
version = v"1.35.0"
beefalo_ver = "1.35.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/macd/beefalo.c.git", "d48ba216e50d7b0854bf7b6c829f04d919567411")
]

# Bash recipe for building across all platforms
script = raw"""
if [[ "${target}" == *-apple-* ]]; then
    EXFLG="-D MACOS";
elif [[ "${target}" == *-w64-* ]]; then
    EXFLG="-D OSWIN";
else
    EXFLG="";
fi

cd ${WORKSPACE}/srcdir/beefalo.c
./configure  --prefix=$prefix --host=${target} --build=${MACHTYPE} CFLAGS="${CFLAGS} ${EXFLG}"
make
mkdir -p "${bindir}"
cp src/beefalo${exeext} ${bindir}
chmod +x ${bindir}/*
install_license COPYING
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> !(Sys.islinux(p) && arch(p) == "i686" &&
                           libc(p) == "musl"), supported_platforms())

# The products that we will ensure are always built
products = [
    ExecutableProduct("beefalo", :beefalo),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using LinearAlgebra

Δx = 0.04
Lx = 4.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx)

RadC = 1.0
Δs = 1.4*cellsize(g)
body = Circle(RadC,Δs)
X = VectorData(collect(body,endpoints=true))

w = Nodes(Dual,size(g))
ϕ = Nodes(Primal,size(g))
q = Edges(Primal,size(g))
f = ScalarData(X)

angs(n) = range(0,2π,length=n+1)[1:n]

_size(::CartesianGrids.Laplacian{NX,NY}) where {NX,NY} = NX, NY

@testset "Forming a cache" begin
  scache1 = SurfaceScalarCache(body,g,scaling=GridScaling)
  scache2 = SurfaceScalarCache(X,g,scaling=GridScaling)
  @test normals(scache1) ≈ normals(scache2)
  @test areas(scache1) ≈ areas(scache2)
  @test points(scache1) ≈ points(scache2)

  @test size(g) == _size(scache1.L)
end

@testset "Scalar surface ops" begin
  scache = SurfaceScalarCache(body,g,scaling=GridScaling)

  gd = similar_grid(scache)
  @test typeof(gd) == typeof(scache.gdata_cache)

  sd = similar_surface(scache)
  @test typeof(sd) == typeof(scache.sdata_cache)

  f .= sin.(angs(length(body)) .- π/4)
  regularize_normal!(q,f,scache)
  normal_interpolate!(f,q,scache)
  @test maximum(f) ≈ 11 atol = 1e-0
  @test minimum(f) ≈ -11 atol = 1e-0

  A = create_CLinvCT(scache,scale=cellsize(g))

  @test maximum(abs.(eigvals(A))) ≈ 0.2 atol = 1e-1

  A = create_GLinvD(scache,scale=cellsize(g))
  @test maximum(abs.(eigvals(A))) ≈ 0.45 atol = 1e-1

  A = create_nRTRn(scache)
  @test maximum(svdvals(A)) ≈ 11 atol = 1e-0

  # Create additional regularization and interpolation ops
  w = ones_gridcurl(scache)
  Rn = RegularizationMatrix(scache,scache.sdata_cache,w)
  En = InterpolationMatrix(scache,w,scache.sdata_cache)

  f .= En*w

  vd = zeros_gridgradcurl(scache)
  @test typeof(vd) <: Edges{Dual}

  f2 = zeros_griddiv(scache)
  @test typeof(f2) <: Nothing

  f2 = similar_griddiv(scache)
  @test typeof(f2) <: Nothing

  f2 = ones_griddiv(scache)
  @test typeof(f2) <: Nothing

  f2 = x_griddiv(scache)
  @test typeof(f2) <: Nothing

  f2 = y_griddiv(scache)
  @test typeof(f2) <: Nothing

  @test isnothing(scache.Rcurl)
  @test isnothing(scache.Rdiv)

end

@testset "Complex cache" begin
  scache = SurfaceScalarCache(body,g,scaling=GridScaling,dtype=ComplexF64)

  @test eltype(similar_grid(scache)) == ComplexF64
  @test eltype(similar_gridcurl(scache)) == ComplexF64
  @test eltype(similar_gridgrad(scache)) == ComplexF64
  @test eltype(similar_gridgradcurl(scache)) == ComplexF64
  @test eltype(similar_surface(scache)) == ComplexF64

end

dq = EdgeGradient(Primal,size(g))
vs = VectorData(X)

@testset "Vector surface ops" begin
  vcache = SurfaceVectorCache(body,g,scaling=GridScaling)

  vs.u .= sin.(angs(length(body)) .- π/4)

  regularize_normal!(dq,vs,vcache)
  normal_interpolate!(vs,dq,vcache)

  fill!(dq,0.0)
  vs.u .= sin.(angs(length(body)) .- π/4)
  regularize_normal_symm!(dq,vs,vcache)
  normal_interpolate_symm!(vs,dq,vcache)
  @test minimum(vs.u) ≈ -22.5 atol = 1e-0
  @test maximum(vs.u) ≈ 22.5 atol = 1e-0

  surface_divergence!(q,vs,vcache)
  surface_grad!(vs,q,vcache)

  surface_divergence_symm!(q,vs,vcache)
  surface_grad_symm!(vs,q,vcache)

  A = create_GLinvD(vcache,scale=cellsize(g))
  @test maximum(abs.(eigvals(A))) ≈ 0.45 atol = 1e-1

  f2 = zeros_griddiv(vcache)
  @test typeof(f2) <: Nodes{Primal}

  f2 = similar_griddiv(vcache)
  @test typeof(f2) <: Nodes{Primal}

  f2 = ones_griddiv(vcache)
  @test typeof(f2) <: Nodes{Primal}

  f2 = x_griddiv(vcache)
  @test typeof(f2) <: Nodes{Primal}

  f2 = y_griddiv(vcache)
  @test typeof(f2) <: Nodes{Primal}

  @test typeof(vcache.Rcurl) <: RegularizationMatrix
  @test typeof(vcache.Rdiv) <: RegularizationMatrix

  v2 = zeros_grid(vcache)
  f2 = ones_surfacescalar(vcache)

  regularize_normal_cross!(v2,f2,vcache)

  vs2 = ones_surface(vcache,2)
  w2 = zeros_gridcurl(vcache)

  regularize_normal_cross!(w2,vs2,vcache)

end

@testset "Masks" begin

  scache = SurfaceScalarCache(body,g,scaling=GridScaling)

  inner = mask(scache)

  oc = similar(ϕ)
  oc .= 1
  @test abs(dot(oc,inner,g) - π*RadC^2) < 2e-3

  outer = complementary_mask(scache)

end

@testset "Grid operations" begin
  vcache = SurfaceVectorCache(body,g,scaling=GridScaling)

  q .= randn(size(q))
  cv_cache = ConvectiveDerivativeCache(EdgeGradient(Primal,q))
  cv_cache2 = ConvectiveDerivativeCache(vcache)
  qdq = similar(q)
  qdq .= 0.0
  convective_derivative!(qdq,q,vcache,cv_cache2)

  qdq2 = convective_derivative(q,vcache)

  @test qdq == qdq2

  w = zeros_gridcurl(vcache)
  v = zeros_grid(vcache)
  v .= randn(size(v))
  curl!(w,v)
  vdw = zeros_gridcurl(vcache)

  cd_cache = ConvectiveDerivativeCache(zeros_gridgradcurl(vcache))
  convective_derivative!(vdw,v,w,vcache,cd_cache)

end

@testset "Convenience fields" begin
  scache = SurfaceScalarCache(body,g,scaling=GridScaling)

  x = x_grid(scache)
  y = y_grid(scache)
  @test typeof(x) == typeof(y) == typeof(scache.gdata_cache)

  x = x_gridcurl(scache)
  y = y_gridcurl(scache)
  @test typeof(x) == typeof(y) == typeof(scache.gcurl_cache)

  x = x_gridgrad(scache)
  y = y_gridgrad(scache)
  @test typeof(x) == typeof(y) == typeof(scache.gsnorm_cache)

  vcache = SurfaceVectorCache(body,g,scaling=GridScaling)

  x = x_grid(vcache)
  y = y_grid(vcache)
  @test typeof(x) == typeof(y) == typeof(vcache.gdata_cache)

  x = x_gridcurl(vcache)
  y = y_gridcurl(vcache)
  @test typeof(x) == typeof(y) == typeof(vcache.gcurl_cache)

  x = x_gridgrad(vcache)
  y = y_gridgrad(vcache)
  @test typeof(x) == typeof(y) == typeof(vcache.gsnorm_cache)

end

@testset "Helmholtz decomposition" begin

  vcache = SurfaceVectorCache(body,g,scaling=GridScaling)
  veccache = VectorFieldCache(vcache)

  w = zeros_gridcurl(vcache)
  ψ = zeros_gridcurl(vcache)
  d = zeros_griddiv(vcache)
  ϕ = zeros_griddiv(vcache)
  v = zeros_grid(vcache)
  dv = zeros_surface(vcache)

  w .= randn(size(w))
  d .= randn(size(d))

  dv .= randn(size(dv))

  vecfield_helmholtz!(v,w,d,dv,(0.0,0.0),vcache,veccache)

  # Test that curl of the resulting vector field is zero
  w2 = zeros_gridcurl(vcache)
  curl!(w2,v,vcache)
  masked_w = zeros_gridcurl(vcache)
  masked_curlv_from_curlv_masked!(masked_w,w2,dv,vcache,veccache.wcache)
  @test maximum(abs.(masked_w[2:end-1,2:end-1] - w[2:end-1,2:end-1])) < 1e-8

  # Test that divergence of the resulting vector field is zero
  d2 = zeros_griddiv(vcache)
  divergence!(d2,v,vcache)
  masked_d = zeros_griddiv(vcache)
  masked_divv_from_divv_masked!(masked_d,d2,dv,vcache,veccache.dcache)
  @test maximum(abs.(masked_d[2:end-1,2:end-1]-d[2:end-1,2:end-1])) < 1e-8

  vectorpotential_from_masked_curlv!(ψ,w,dv,vcache,veccache.wcache)
  scalarpotential_from_masked_divv!(ϕ,d,dv,vcache,veccache.dcache)

  # Tests that no immersed points lead to no problems
  vcache = SurfaceVectorCache(g,scaling=GridScaling)
  veccache = VectorFieldCache(vcache)

  w = zeros_gridcurl(vcache)
  ψ = zeros_gridcurl(vcache)
  d = zeros_griddiv(vcache)
  ϕ = zeros_griddiv(vcache)
  v = zeros_grid(vcache)
  dv = zeros_surface(vcache)

  vectorpotential_from_masked_curlv!(ψ,w,dv,vcache,veccache.wcache)
  scalarpotential_from_masked_divv!(ϕ,d,dv,vcache,veccache.dcache)
  vecfield_helmholtz!(v,w,d,dv,(0.0,0.0),vcache,veccache)
  @test maximum(abs.(ψ)) == 0.0
  @test maximum(abs.(ϕ)) == 0.0
  @test maximum(abs.(v)) == 0.0

end

@testset "Problem specification" begin

  prob = BasicScalarILMProblem(g,scaling=GridScaling)
  @test typeof(prob.bodies) <: BodyList

  prob = BasicScalarILMProblem(g,body,scaling=GridScaling)
  @test typeof(prob.bodies) <: BodyList

  sys = construct_system(prob)

  scache = SurfaceScalarCache(body,g,scaling=GridScaling)

  @test typeof(sys.base_cache) == typeof(scache)

  vprob = BasicVectorILMProblem(g,body,scaling=GridScaling)
  vsys = construct_system(vprob)

  vcache = SurfaceVectorCache(body,g,scaling=GridScaling)

  @test typeof(vsys.base_cache) == typeof(vcache)

  p = rand()
  prob = BasicScalarILMProblem(g,body,scaling=GridScaling,phys_params=p)
  sys = construct_system(prob)
  @test sys.phys_params == p

  Re = rand()
  p = Dict("Re"=>Re)
  prob = BasicScalarILMProblem(g,body,scaling=GridScaling,phys_params=p)
  sys = construct_system(prob)
  @test sys.phys_params["Re"] == Re

  my_bc_func(x) = x
  prob = BasicScalarILMProblem(g,body,scaling=GridScaling,bc=my_bc_func)
  sys = construct_system(prob)
  @test sys.bc == my_bc_func

  my_f_func(x) = x
  prob = BasicScalarILMProblem(g,body,scaling=GridScaling,forcing=my_f_func)
  sys = construct_system(prob)
  @test sys.forcing == my_f_func

  prob = BasicScalarILMProblem(g,body,scaling=GridScaling,dtype=ComplexF64)
  sys = construct_system(prob)
  @test eltype(sys.base_cache.gdata_cache) == ComplexF64

  prob = BasicVectorILMProblem(g,body,scaling=GridScaling,dtype=ComplexF64)
  sys = construct_system(prob)
  @test eltype(sys.base_cache.gdata_cache) == ComplexF64

end

@testset "Problem regeneration" begin

  my_f_func(x) = x
  prob = BasicScalarILMProblem(g,body,scaling=GridScaling,forcing=my_f_func)
  sys = construct_system(prob)

  prob2 = ImmersedLayers.regenerate_problem(sys,sys.base_cache.bl)
  @test prob.g == prob2.g
  @test prob.bodies == prob2.bodies
  @test prob.phys_params == prob2.phys_params
  @test prob.bc == prob2.bc
  @test prob.forcing == prob2.forcing
  @test prob.timestep_func == prob2.timestep_func

end

@testset "Forcing regions" begin
  Δs = 1.4*cellsize(g)

  scache = SurfaceScalarCache(body,g,scaling=GridScaling)

  phys_params = Dict("areaheater1_flux" => 3.0,
                      "lineheater_flux" => -2.0,
                      "areaheater2_temp" => 1.5,
                       "areaheater2_coeff" => 2.0)

  fregion1 = Circle(0.2,Δs)
  T = RigidTransform((0.5,0.0),0.0)
  T(fregion1)
  function model1!(σ,T,t,fr::AreaRegionCache,phys_params)
      σ .= phys_params["areaheater1_flux"]
   end
   afm = AreaForcingModel(fregion1,model1!)

   fregion2 = Square(0.5,Δs)
   T = RigidTransform((0.0,1.0),0.0)
   T(fregion2)
   function model2!(σ,T,t,fr::LineRegionCache,phys_params)
     σ .= phys_params["lineheater_flux"]
    end
   lfm = LineForcingModel(fregion2,model2!)

   fregion3 = Polygon([0.0,1.0,0.5],[0.0,0.0,0.5],Δs)
   T = RigidTransform((-1.0,-1.0),π/4)
   T(fregion3)
   function model3!(σ,T,t,fr::AreaRegionCache,phys_params)
     σ .= phys_params["areaheater2_coeff"]*(phys_params["areaheater2_temp"] - T)
   end
   afm2 = AreaForcingModel(fregion3,model3!)

   pts = VectorData([-1.2,0.5],[0.5,0.5])
   str = [1.0,-1.0]

   function model4!(σ,T,t,fr::PointRegionCache,phys_params)
     σ .= str
     return σ
   end
   pfm = PointForcingModel(pts,model4!;ddftype=CartesianGrids.M4prime);

   fcache = ForcingModelAndRegion([afm,lfm,afm2,pfm],scache)

   @test typeof(fcache[1].region_cache.mask) <: ScalarGridData

   @test typeof(fcache[4].region_cache.cache.regop.ddf) == DDF{CartesianGrids.M4prime, 1.0}

   T = zeros_grid(scache)
   dT = similar_grid(scache)
   t = 0.0

   apply_forcing!(dT,T,t,fcache,phys_params)

   function point_function(T,t,fr::PointRegionCache,phys_params)
     return VectorData([-1.2,0.5],[0.5,0.5])
   end
   pfm2 = PointForcingModel(point_function,model4!;ddftype=CartesianGrids.M4prime);

   fcache = ForcingModelAndRegion(pfm,scache)
   fcache2 = ForcingModelAndRegion(pfm2,scache)

   T2 = zeros_grid(scache)
   dT2 = similar_grid(scache)
   t = 0.0

   apply_forcing!(dT,T,t,fcache,phys_params)
   apply_forcing!(dT2,T2,t,fcache2,phys_params)

   @test dT == dT2

   str = [1.0,-1.0,2.0]

   @test_throws DimensionMismatch apply_forcing!(dT2,T2,t,fcache2,phys_params)

end


text:
tomlpath = joinpath(@__DIR__, "sbm_config.toml")
config = Wflow.Config(tomlpath)

model = Wflow.initialize_sbm_model(config)

@unpack network = model

min_sto = get(config.model, "min_streamorder", 4)
index_pit = [network.land.order[end]]
@test min_sto == 3

subbas_order, indices_subbas, topo_subbas =
    Wflow.kinwave_set_subdomains(config, network.land.graph, network.land.order, index_pit)

Wflow.close_files(model, delete_output = false)

if nthreads() == 1
    @testset "Nonparallel subdomains kinematic wave (nthreads = 1)" begin
        @test subbas_order == [[1]]
        @test indices_subbas == [[1:length(network.land.order);]]
        @test topo_subbas == [network.land.order]
    end
else
    @testset "Parallel subdomains kinematic wave (nthreads > 1)" begin
        @test length(subbas_order) == 6
        @test length.(subbas_order) == [507, 172, 142, 86, 24, 1]
        @test maximum(subbas_order) == [932]
        @test subbas_order[1][1:4] == [107, 115, 83, 84]
        @test subbas_order[3][1:4] == [97, 109, 56, 68]
        @test length(topo_subbas) == 932
        @test topo_subbas[1][1:4] == [50062, 50063, 50056, 50055]
        @test topo_subbas[end][1:4] == [49884, 49842, 49791, 49779]
        @test length(indices_subbas) == 932
        @test indices_subbas[1][1:4] == [1, 2, 8, 9]
        @test indices_subbas[end][1:4] == [166, 201, 236, 247]
    end
end

@testset "Streamorder and subbasins" begin
    # directed acyclic graph of basin
    g = DiGraph(16)
    add_edge!(g, 1, 3)
    add_edge!(g, 2, 3)
    add_edge!(g, 3, 4)
    add_edge!(g, 4, 5)
    add_edge!(g, 5, 6)
    add_edge!(g, 7, 9)
    add_edge!(g, 8, 9)
    add_edge!(g, 9, 4)
    add_edge!(g, 10, 12)
    add_edge!(g, 11, 12)
    add_edge!(g, 12, 16)
    add_edge!(g, 13, 15)
    add_edge!(g, 14, 15)
    add_edge!(g, 15, 16)
    add_edge!(g, 16, 5)

    # generate streamorder and subbasins (ids, graph and order)
    toposort = topological_sort_by_dfs(g)
    strord = Wflow.stream_order(g, toposort)
    subbas = Wflow.subbasins(g, strord, toposort, 2)
    subbas_fill = Wflow.fillnodata_upstream(g, toposort, subbas, 0)
    graph_subbas = Wflow.graph_from_nodes(g, subbas, subbas_fill)
    toposort_subbas = topological_sort_by_dfs(graph_subbas)
    dist =
        Graphs.Experimental.Traversals.distances(Graph(graph_subbas), toposort_subbas[end])
    max_dist = maximum([dist; 1])
    subbas_order = Wflow.subbasins_order(graph_subbas, toposort_subbas[end], max_dist)

    @test strord == [1, 1, 2, 3, 4, 4, 1, 1, 2, 1, 1, 2, 1, 1, 2, 3]
    @test strord[toposort[end]] == 4
    @test subbas == [0, 0, 5, 6, 0, 7, 0, 0, 4, 0, 0, 2, 0, 0, 1, 3]
    @test subbas_fill == [5, 5, 5, 6, 7, 7, 4, 4, 4, 2, 2, 2, 1, 1, 1, 3]
    @test toposort_subbas == [5, 4, 6, 2, 1, 3, 7]
    @test max_dist == 2
    @test subbas_order == [[1, 2, 4, 5], [3, 6], [7]]
end


text:
function build_array(::Type{MathArrays}, arg0::Field, arg1::jint)
    return jcall(MathArrays, "buildArray", Vector{Object}, (Field, jint), arg0, arg1)
end

function build_array(::Type{MathArrays}, arg0::Field, arg1::jint, arg2::jint)
    return jcall(MathArrays, "buildArray", Vector{Vector{Object}}, (Field, jint, jint), arg0, arg1, arg2)
end

function build_array(::Type{MathArrays}, arg0::Field, arg1::jint, arg2::jint, arg3::jint)
    return jcall(MathArrays, "buildArray", Vector{Vector{Object}}, (Field, jint, jint, jint), arg0, arg1, arg2, arg3)
end

function check_equal_length(::Type{MathArrays}, arg0::Vector{RealFieldElement}, arg1::Vector{RealFieldElement})
    return jcall(MathArrays, "checkEqualLength", void, (Vector{RealFieldElement}, Vector{RealFieldElement}), arg0, arg1)
end

function check_equal_length(::Type{MathArrays}, arg0::Vector{RealFieldElement}, arg1::Vector{RealFieldElement}, arg2::jboolean)
    return jcall(MathArrays, "checkEqualLength", jboolean, (Vector{RealFieldElement}, Vector{RealFieldElement}, jboolean), arg0, arg1, arg2)
end

function check_equal_length(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(MathArrays, "checkEqualLength", void, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function check_equal_length(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jboolean)
    return jcall(MathArrays, "checkEqualLength", jboolean, (Vector{jdouble}, Vector{jdouble}, jboolean), arg0, arg1, arg2)
end

function check_equal_length(::Type{MathArrays}, arg0::Vector{jint}, arg1::Vector{jint})
    return jcall(MathArrays, "checkEqualLength", void, (Vector{jint}, Vector{jint}), arg0, arg1)
end

function check_equal_length(::Type{MathArrays}, arg0::Vector{jint}, arg1::Vector{jint}, arg2::jboolean)
    return jcall(MathArrays, "checkEqualLength", jboolean, (Vector{jint}, Vector{jint}, jboolean), arg0, arg1, arg2)
end

function check_non_negative(::Type{MathArrays}, arg0::Vector{Vector{jlong}})
    return jcall(MathArrays, "checkNonNegative", void, (Vector{Vector{jlong}},), arg0)
end

function check_non_negative(::Type{MathArrays}, arg0::Vector{jlong})
    return jcall(MathArrays, "checkNonNegative", void, (Vector{jlong},), arg0)
end

function check_not_na_n(::Type{MathArrays}, arg0::Vector{jdouble})
    return jcall(MathArrays, "checkNotNaN", void, (Vector{jdouble},), arg0)
end

function check_order(::Type{MathArrays}, arg0::Vector{RealFieldElement})
    return jcall(MathArrays, "checkOrder", void, (Vector{RealFieldElement},), arg0)
end

function check_order(::Type{MathArrays}, arg0::Vector{RealFieldElement}, arg1::MathArrays_OrderDirection, arg2::jboolean)
    return jcall(MathArrays, "checkOrder", void, (Vector{RealFieldElement}, MathArrays_OrderDirection, jboolean), arg0, arg1, arg2)
end

function check_order(::Type{MathArrays}, arg0::Vector{RealFieldElement}, arg1::MathArrays_OrderDirection, arg2::jboolean, arg3::jboolean)
    return jcall(MathArrays, "checkOrder", jboolean, (Vector{RealFieldElement}, MathArrays_OrderDirection, jboolean, jboolean), arg0, arg1, arg2, arg3)
end

function check_order(::Type{MathArrays}, arg0::Vector{jdouble})
    return jcall(MathArrays, "checkOrder", void, (Vector{jdouble},), arg0)
end

function check_order(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::MathArrays_OrderDirection, arg2::jboolean)
    return jcall(MathArrays, "checkOrder", void, (Vector{jdouble}, MathArrays_OrderDirection, jboolean), arg0, arg1, arg2)
end

function check_order(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::MathArrays_OrderDirection, arg2::jboolean, arg3::jboolean)
    return jcall(MathArrays, "checkOrder", jboolean, (Vector{jdouble}, MathArrays_OrderDirection, jboolean, jboolean), arg0, arg1, arg2, arg3)
end

function check_positive(::Type{MathArrays}, arg0::Vector{jdouble})
    return jcall(MathArrays, "checkPositive", void, (Vector{jdouble},), arg0)
end

function check_rectangular(::Type{MathArrays}, arg0::Vector{Vector{jlong}})
    return jcall(MathArrays, "checkRectangular", void, (Vector{Vector{jlong}},), arg0)
end

function concatenate(::Type{MathArrays}, arg0::Vector{Vector{jdouble}})
    return jcall(MathArrays, "concatenate", Vector{jdouble}, (Vector{Vector{jdouble}},), arg0)
end

function convolve(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(MathArrays, "convolve", Vector{jdouble}, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function cos_angle(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(MathArrays, "cosAngle", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function distance(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(MathArrays, "distance", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function distance(::Type{MathArrays}, arg0::Vector{jint}, arg1::Vector{jint})
    return jcall(MathArrays, "distance", jdouble, (Vector{jint}, Vector{jint}), arg0, arg1)
end

function distance1(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(MathArrays, "distance1", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function distance1(::Type{MathArrays}, arg0::Vector{jint}, arg1::Vector{jint})
    return jcall(MathArrays, "distance1", jint, (Vector{jint}, Vector{jint}), arg0, arg1)
end

function distance_inf(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(MathArrays, "distanceInf", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function distance_inf(::Type{MathArrays}, arg0::Vector{jint}, arg1::Vector{jint})
    return jcall(MathArrays, "distanceInf", jint, (Vector{jint}, Vector{jint}), arg0, arg1)
end

function ebe_add(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(MathArrays, "ebeAdd", Vector{jdouble}, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function ebe_divide(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(MathArrays, "ebeDivide", Vector{jdouble}, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function ebe_multiply(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(MathArrays, "ebeMultiply", Vector{jdouble}, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function ebe_subtract(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(MathArrays, "ebeSubtract", Vector{jdouble}, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function equals(::Type{MathArrays}, arg0::Vector{jbyte}, arg1::Vector{jbyte})
    return jcall(MathArrays, "equals", jboolean, (Vector{jbyte}, Vector{jbyte}), arg0, arg1)
end

function equals(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(MathArrays, "equals", jboolean, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function equals(::Type{MathArrays}, arg0::Vector{jfloat}, arg1::Vector{jfloat})
    return jcall(MathArrays, "equals", jboolean, (Vector{jfloat}, Vector{jfloat}), arg0, arg1)
end

function equals(::Type{MathArrays}, arg0::Vector{jint}, arg1::Vector{jint})
    return jcall(MathArrays, "equals", jboolean, (Vector{jint}, Vector{jint}), arg0, arg1)
end

function equals(::Type{MathArrays}, arg0::Vector{jlong}, arg1::Vector{jlong})
    return jcall(MathArrays, "equals", jboolean, (Vector{jlong}, Vector{jlong}), arg0, arg1)
end

function equals(::Type{MathArrays}, arg0::Vector{jshort}, arg1::Vector{jshort})
    return jcall(MathArrays, "equals", jboolean, (Vector{jshort}, Vector{jshort}), arg0, arg1)
end

function equals_including_na_n(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(MathArrays, "equalsIncludingNaN", jboolean, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function equals_including_na_n(::Type{MathArrays}, arg0::Vector{jfloat}, arg1::Vector{jfloat})
    return jcall(MathArrays, "equalsIncludingNaN", jboolean, (Vector{jfloat}, Vector{jfloat}), arg0, arg1)
end

function is_monotonic(::Type{MathArrays}, arg0::Vector{Comparable}, arg1::MathArrays_OrderDirection, arg2::jboolean)
    return jcall(MathArrays, "isMonotonic", jboolean, (Vector{Comparable}, MathArrays_OrderDirection, jboolean), arg0, arg1, arg2)
end

function is_monotonic(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::MathArrays_OrderDirection, arg2::jboolean)
    return jcall(MathArrays, "isMonotonic", jboolean, (Vector{jdouble}, MathArrays_OrderDirection, jboolean), arg0, arg1, arg2)
end

function linear_combination(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(MathArrays, "linearCombination", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function linear_combination(::Type{MathArrays}, arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jdouble)
    return jcall(MathArrays, "linearCombination", jdouble, (jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function linear_combination(::Type{MathArrays}, arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble)
    return jcall(MathArrays, "linearCombination", jdouble, (jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5)
end

function linear_combination(::Type{MathArrays}, arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble)
    return jcall(MathArrays, "linearCombination", jdouble, (jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function natural(::Type{MathArrays}, arg0::jint)
    return jcall(MathArrays, "natural", Vector{jint}, (jint,), arg0)
end

function normalize_array(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::jdouble)
    return jcall(MathArrays, "normalizeArray", Vector{jdouble}, (Vector{jdouble}, jdouble), arg0, arg1)
end

function safe_norm(::Type{MathArrays}, arg0::Vector{jdouble})
    return jcall(MathArrays, "safeNorm", jdouble, (Vector{jdouble},), arg0)
end

function scale(::Type{MathArrays}, arg0::jdouble, arg1::Vector{jdouble})
    return jcall(MathArrays, "scale", Vector{jdouble}, (jdouble, Vector{jdouble}), arg0, arg1)
end

function scale_in_place(::Type{MathArrays}, arg0::jdouble, arg1::Vector{jdouble})
    return jcall(MathArrays, "scaleInPlace", void, (jdouble, Vector{jdouble}), arg0, arg1)
end

function sequence(::Type{MathArrays}, arg0::jint, arg1::jint, arg2::jint)
    return jcall(MathArrays, "sequence", Vector{jint}, (jint, jint, jint), arg0, arg1, arg2)
end

function shuffle(::Type{MathArrays}, arg0::Vector{jint})
    return jcall(MathArrays, "shuffle", void, (Vector{jint},), arg0)
end

function shuffle(::Type{MathArrays}, arg0::Vector{jint}, arg1::RandomGenerator)
    return jcall(MathArrays, "shuffle", void, (Vector{jint}, RandomGenerator), arg0, arg1)
end

function shuffle(::Type{MathArrays}, arg0::Vector{jint}, arg1::jint, arg2::MathArrays_Position)
    return jcall(MathArrays, "shuffle", void, (Vector{jint}, jint, MathArrays_Position), arg0, arg1, arg2)
end

function shuffle(::Type{MathArrays}, arg0::Vector{jint}, arg1::jint, arg2::MathArrays_Position, arg3::RandomGenerator)
    return jcall(MathArrays, "shuffle", void, (Vector{jint}, jint, MathArrays_Position, RandomGenerator), arg0, arg1, arg2, arg3)
end

function sort_in_place(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::MathArrays_OrderDirection, arg2::Vector{Vector{jdouble}})
    return jcall(MathArrays, "sortInPlace", void, (Vector{jdouble}, MathArrays_OrderDirection, Vector{Vector{jdouble}}), arg0, arg1, arg2)
end

function sort_in_place(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::Vector{Vector{jdouble}})
    return jcall(MathArrays, "sortInPlace", void, (Vector{jdouble}, Vector{Vector{jdouble}}), arg0, arg1)
end

function unique(::Type{MathArrays}, arg0::Vector{jdouble})
    return jcall(MathArrays, "unique", Vector{jdouble}, (Vector{jdouble},), arg0)
end

function verify_values(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jint, arg3::jint)
    return jcall(MathArrays, "verifyValues", jboolean, (Vector{jdouble}, Vector{jdouble}, jint, jint), arg0, arg1, arg2, arg3)
end

function verify_values(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jint, arg3::jint, arg4::jboolean)
    return jcall(MathArrays, "verifyValues", jboolean, (Vector{jdouble}, Vector{jdouble}, jint, jint, jboolean), arg0, arg1, arg2, arg3, arg4)
end

function verify_values(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(MathArrays, "verifyValues", jboolean, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function verify_values(::Type{MathArrays}, arg0::Vector{jdouble}, arg1::jint, arg2::jint, arg3::jboolean)
    return jcall(MathArrays, "verifyValues", jboolean, (Vector{jdouble}, jint, jint, jboolean), arg0, arg1, arg2, arg3)
end


text:
Random.seed!(0)

# Dynamics
dt = 0.1
m = 50
A = [1 dt 0 0 ; 0 1 0 0 ; 0 0 1 dt; 0 0 0 1]
B = [0 0; dt/m 0; 0 0; 0 dt/m]
W = 0.1*Matrix{Float64}(I,4,4)

dmodel = LinearDynamicsModel(A,B,W);

# Measurement
C = [0 1.0 0 0; 0 0 0 1.0]
V = 0.5*Matrix{Float64}(I,2,2)

omodel = LinearObservationModel(C,V)

# Build KF
kf = KalmanFilter(dmodel,omodel)

# Simulation
b0 = GaussianBelief([0.0,0.0,0.0,0.0], 2.0*Matrix{Float64}(I,4,4))

times = 0:dt:10
Fmag = 1000
action_sequence = [[Fmag*cos(t), Fmag*sin(t)] for t in times]

sim_states, sim_measurements = simulation(kf, b0,action_sequence);

# Filter
filtered_beliefs = run_filter(kf, b0, action_sequence, sim_measurements)
μ, Σ = unpack(filtered_beliefs);

# Tests
@test kf.d isa LinearDynamicsModel
@test kf.o isa LinearObservationModel
@test length(sim_measurements) == length(action_sequence)

first_update = update(kf, b0, action_sequence[1], sim_measurements[1])
first_predict_measure = measure(kf, predict(kf, b0, action_sequence[1]),
    sim_measurements[1]; u = action_sequence[1])
array_isapprox(first_update.μ, first_predict_measure.μ)
array_isapprox(first_update.Σ, first_predict_measure.Σ)

@test length(filtered_beliefs[end].μ) == length(b0.μ)
@test size(filtered_beliefs[end].Σ) == size(b0.Σ)
@test size(μ) == (length(filtered_beliefs), length(filtered_beliefs[1].μ))
@test size(Σ) == (length(filtered_beliefs), size(filtered_beliefs[1].Σ)...)

# change with algorithm / rng changes
array_isapprox(sim_states[end], [59.9223, -7.17527, 178.149, 30.4653];atol=0.001)
array_isapprox(sim_measurements[end], [-7.01552, 30.4259];atol=0.001)
array_isapprox(filtered_beliefs[end].μ, [60.1627, -7.13376, 183.227, 30.3895];
                atol=0.001)
array_isapprox(filtered_beliefs[end].Σ,   [12.6069 0.0320871 0.0 0.0;
                                            0.0320871 0.179129 0.0 0.0;
                                            0.0 0.0 12.6069 0.0320871;
                                            0.0 0.0 0.0320871 0.179129]; atol=0.001)


text:
##################################################################
#    IMPORTANT: ANALYTICAL QUERIES THAT ARE CALLED THROUGHOUT    #
##################################################################

#delimiter to separate ACS fields that need to be aggregated
delim_agg_fields = ";;"
#overwrite data if they exist?
overwrite_data_q = true
#overwrite fp_csv_ca_pl_unused if it exists?
overwrite_unused_fields_q = true
#use existing dataset in "build restore_existing_database_q" notebook?
restore_existing_database_q = false

#########################
#    SET DIRECTORIES    #
#########################

#set the acs year
acs_yr = 2010
census_yr = 2010

#high level
dir_cur = dirname(pwd())
dir_data = joinpath(dirname(dir_cur), "data")
dir_bin = joinpath(dir_cur, "bin")
dir_acs = joinpath(dir_data, "ACS_$(acs_yr)")
dir_ref = joinpath(dir_cur, "ref")
dir_ref_acs = joinpath(dir_ref, "ACS_$(acs_yr)")
dir_out = joinpath(dir_cur, "out")
dir_out_img = joinpath(dir_out, "img")

#acs subdirectories
dir_acs_templates = joinpath(dir_acs, "$(acs_yr)_5yr_Summary_FileTemplates")
dir_acs_data = joinpath(dir_acs, "California_Tracts_Block_Groups_Only")
dir_data_estimates = joinpath(dir_data, "preliminary_downscaled_estimates", "ACS-$(acs_yr)_TO_CENSUS-$(census_yr)")
#california
dir_cvap = joinpath(dir_data, "CA_block_level")
#set directory for state blocks
dir_sblks = joinpath(dir_data, "state_blocks_$(census_yr)_shp")
#pl 94 for california
dir_capl = joinpath(dir_data, "ca$(census_yr).pl")
dir_casf = joinpath(dir_data, "ca$(census_yr).sf1")

# check that directories exist/make them if need be
for d_check in [dir_out, dir_out_img, dir_data_estimates]
    if !ispath(d_check)
        mkdir(d_check)
    end
end

#######################
#    SET FILEPATHS    #
#######################

##  CSVs
fp_csv_acs_field_sequence_index = joinpath(dir_acs, "acs_field_sequence_index_$(acs_yr).csv")
fp_csv_acs_geodat = joinpath(dir_acs_data, "g$(acs_yr)5ca.csv")
fp_csv_adjacecy_index = joinpath(dir_data, "adjacency_index.csv")

fp_csv_aers_counts_stage_1 = joinpath(dir_out, "aers_counts_acs-$(acs_yr)_to_census-$(census_yr)_stage-1.csv")
fp_csv_aers_counts_stage_2 = joinpath(dir_out, "aers_counts_acs-$(acs_yr)_to_census-$(census_yr)_stage-2.csv")
fp_csv_aers_termination_status_stage_1 = joinpath(dir_out, "aers_termination_status_stage-1.csv")
fp_csv_aers_termination_status_stage_2 = joinpath(dir_out, "aers_termination_status_stage-2.csv")

fp_csv_attribute_census_sf1_table = joinpath(dir_ref, "attribute_census_sf1_table_id.csv")
fp_csv_attribute_demographic = joinpath(dir_ref, "attribute_demographic_id.csv")
fp_csv_attribute_household = joinpath(dir_ref_acs, "attribute_household_id.csv")
fp_csv_attribute_geo_unit = joinpath(dir_ref, "attribute_geo_unit.csv")
fp_csv_attribute_variable = joinpath(dir_ref_acs, "attribute_variable_id.csv")

fp_csv_ca_pl94_template = joinpath(dir_capl, "ca0000##SEGMENT##$(census_yr).pl")
fp_csv_ca_pl_unused = joinpath(dir_data_estimates, "california_pl_94_data_codes-unused.csv")
fp_csv_composite_data_block_group = joinpath(dir_data_estimates, "composite_data-block_group.csv")
fp_csv_composite_data_tract = joinpath(dir_data_estimates, "composite_data-tract.csv")
fp_csv_codemap_sf1_aggregation = joinpath(dir_ref, "codemap_sf1_aggregation.csv")
fp_csv_codemap_sf1_to_acs = joinpath(dir_ref_acs, "codemap_sf1_to_acs.csv")
fp_csv_extraction_fields_acs = joinpath(dir_ref_acs, "extraction_fields-acs_$(acs_yr).csv")
fp_csv_extraction_fields_file_reduce = joinpath(dir_ref, "extraction_fields_2020_pl94_file_reduce.csv")
fp_csv_geo_ca = joinpath(dir_acs, "geo_ca.csv")
fp_csv_out_constraint_matrix = joinpath(dir_data_estimates, "base_dasymmetric_constraints.csv")
fp_csv_preliminary_block_household_variable_estimates = joinpath(dir_data_estimates, "california_preliminary_block_household_variable_estimates.csv")
fp_csv_preliminary_block_population_variable_estimates = joinpath(dir_data_estimates, "california_preliminary_block_population_variable_estimates.csv")
fp_csv_ref_constraint_matrix_manual = joinpath(dir_ref_acs, "manual_dasymmetric_constraints.csv")
fp_csv_summing_matrix = joinpath(dir_ref, "summing_matrix-demographic_id.csv")

##  SHAPEFILEs
fp_shp_state_blocks = joinpath(dir_sblks, "state_blks$(census_yr - 2000).shp")

##  OTHER TEXT FILES
fp_fwtxt_ca_pl94_geoheader = joinpath(dir_capl, "cageo$(census_yr).pl")
fp_fwtxt_ca_sf1_geoheader = joinpath(dir_casf, "cageo$(census_yr).sf1")

##  EXCEL FILES
fp_xlsx_acs_geodat_header = joinpath(dir_acs_templates, "$(acs_yr)_SFGeoFileTemplate.xlsx")


text:
module CellMapApply

using Gridap
using Gridap.Helpers
using Gridap.CachedArrays
using Gridap.CachedValues
using Gridap.CellNumberApply: _getvalues
using Gridap.MapApply: _compute_S, _compute_M, setinputs!
using Gridap.MapApply: MapFromKernel
using Gridap.Kernels: _compute_N, _compute_T

import Gridap: evaluate
import Gridap.MapApply: _stype, _m
import Gridap: apply
import Base: iterate
import Base: length
import Base: size
import Base: getindex

function apply(k::ArrayKernel,m::CellMap,v::Vararg{<:CellValue})
  CellMapFromKernel(k,m,v...)
end

function apply(k::ArrayKernel,m::CellValue,v::CellMap)
  CellMapFromKernel(k,m,v)
end

function apply(k::ArrayKernel,m::CellMap,v::CellMap)
  CellMapFromKernel(k,m,v)
end

function apply(k::ArrayKernel,m::IndexCellMap,v::Vararg{<:IndexCellValue})
  IndexCellMapFromKernel(k,m,v...)
end

function apply(k::ArrayKernel,m::IndexCellValue,v::IndexCellMap)
  IndexCellMapFromKernel(k,m,v)
end

function apply(k::ArrayKernel,m::IndexCellMap,v::IndexCellMap)
  IndexCellMapFromKernel(k,m,v)
end

struct CellMapFromKernel{S,M,T,N,R,K,V} <: IterCellMap{S,M,T,N,R}
  kernel::K
  cellvalues::V
end

function CellMapFromKernel(k::ArrayKernel,v::Vararg{<:CellValue})
  S = _compute_S(v)
  M = _compute_M(v)
  T = _compute_T(k,v)
  N = _compute_N(k,v)
  R = _compute_R(k,v)
  K = typeof(k)
  V = typeof(v)
  CellMapFromKernel{S,M,T,N,R,K,V}(k,v)
end

function length(self::CellMapFromKernel)
  vi, = self.cellvalues
  length(vi)
end

@inline function iterate(self::CellMapFromKernel{T,N}) where {T,N}
  zipped = zip(self.cellvalues...)
  znext = iterate(zipped)
  if znext === nothing; return nothing end
  a, zstate = znext
  v = MapFromKernel(self.kernel,a...)
  _iterate(self,znext,zipped,v)
end

@inline function iterate(self::CellMapFromKernel,state)
  zipped, zstate, v = state
  znext = iterate(zipped,zstate)
  _iterate(self,znext,zipped,v)
end

@inline function _iterate(self::CellMapFromKernel,znext,zipped,v)
  if znext === nothing; return nothing end
  a, zstate = znext
  setinputs!(v,a...)
  state = (zipped,zstate,v)
  (v, state)
end

function evaluate(
  m::CellMapFromKernel{S,M,T,N},a::CellArray{<:S,M}) where {S,M,T,N}
  v = [ _eval(mi,a) for mi in m.cellvalues ]
  apply(m.kernel,v...)
end

_eval(m::CellMap,a) = evaluate(m,a)

_eval(m::CellArray,a) = m

_eval(m::CellNumber,a) = m

function _compute_R(k,v)
  m = tuple([eltype(vi) for vi in v]...)
  S = _compute_S(m)
  M = _compute_M(m)
  T = _compute_T(k,m)
  N = _compute_N(k,m)
  K = typeof(k)
  V = Tuple{m...}
  C = Tuple{[ _cache_type(mi) for mi in m]...}
  MapFromKernel{S,M,T,N,K,V,C}
end

_cache_type(::Type{<:Map{S,M,T,N}}) where {S,M,T,N} = CachedArray{T,N,Array{T,N}}

_cache_type(::Type{<:AbstractArray{T,N}}) where {T,N} = CachedArray{T,N,Array{T,N}}

_cache_type(::Type{T}) where T<:Number = Vector{T}

_stype(v::CellMap{S}) where S = S

_stype(v::CellArray) = nothing

_stype(v::CellNumber) = nothing

_m(v::CellMap{S,M}) where {S,M} = M

_m(v::CellArray) = nothing

_m(v::CellNumber) = nothing

_stype(v::Type{<:Map{S}}) where S = S

_stype(v::Type{<:AbstractArray}) = nothing

_stype(v::Type{<:Number}) = nothing

_m(v::Type{<:Map{S,M}}) where {S,M} = M

_m(v::Type{<:AbstractArray}) = nothing

_m(v::Type{<:Number}) = nothing

struct IndexCellMapFromKernel{S,M,T,N,R<:Map{S,M,T,N},K,V,F} <: IndexCellValue{R,1}
  kernel::K
  cellvalues::V
  cache::CachedValue{F}
end

function IndexCellMapFromKernel(k::ArrayKernel,v::Vararg{<:CellValue})
  S = _compute_S(v)
  M = _compute_M(v)
  T = _compute_T(k,v)
  N = _compute_N(k,v)
  R = _compute_R(k,v)
  K = typeof(k)
  V = typeof(v)
  F = Union{Nothing,R}
  cache = CachedValue{F}(nothing)
  IndexCellMapFromKernel{S,M,T,N,R,K,V,F}(k,v,cache)
end

function length(self::IndexCellMapFromKernel)
  vi, = self.cellvalues
  length(vi)
end

size(self::IndexCellMapFromKernel) = (length(self),)

function getindex(self::IndexCellMapFromKernel,i::Integer)
  vals = _getvalues(i,self.cellvalues...)
  if self.cache.value === nothing
    self.cache.value = MapFromKernel(self.kernel,vals...)
  else
  setinputs!(self.cache.value,vals...)
  end
  self.cache.value
end

function evaluate(
  m::IndexCellMapFromKernel{S,M,T,N},a::IndexCellArray{<:S,M}) where {S,M,T,N}
  v = [ _eval(mi,a) for mi in m.cellvalues ]
  apply(m.kernel,v...)
end

function evaluate(
  m::IndexCellMapFromKernel{S,M,T,N},a::CellArray{<:S,M}) where {S,M,T,N}
  v = [ _eval(mi,a) for mi in m.cellvalues ]
  apply(m.kernel,v...)
end

end # module


text:
@testset "hold" begin
    @testset "basic" begin
        @system SHold(Controller) begin
            a ~ hold
        end
        @test_throws ErrorException instance(SHold)
    end
end


text:
using DelayDiffEq, DiffEqProblemLibrary.DDEProblemLibrary
using Test

DDEProblemLibrary.importddeproblems()

# Check that numerical solutions approximate analytical solutions,
# independent of problem structure

@testset "standard history" begin
  # standard algorithm
  alg = MethodOfSteps(BS3(); constrained = false)

  alg1 = MethodOfSteps(Tsit5(); constrained = false,
                       fpsolve = NLFunctional(; max_iter = 100))
  alg2 = MethodOfSteps(DP8(); constrained = false,
                       fpsolve = NLFunctional(; max_iter = 10))
  alg3 = MethodOfSteps(Tsit5(); constrained = true)
  alg4 = MethodOfSteps(DP5(); constrained = false,
                       fpsolve = NLFunctional(; max_iter = 100))

  ## Single constant delay
  @testset "single constant delay" begin
    @testset "short time span" begin
      ### Scalar function
      sol_scalar = solve(DDEProblemLibrary.prob_dde_constant_1delay_scalar, alg)

      @test sol_scalar.errors[:l∞] < 3.0e-5
      @test sol_scalar.errors[:final] < 2.1e-5
      @test sol_scalar.errors[:l2] < 1.2e-5

      ### Out-of-place function
      sol_oop = solve(DDEProblemLibrary.prob_dde_constant_1delay_oop, alg)

      @test sol_scalar.t ≈ sol_oop.t && sol_scalar.u ≈ sol_oop[1, :]

      ### In-place function
      sol_ip = solve(DDEProblemLibrary.prob_dde_constant_1delay_ip, alg)

      @test sol_scalar.t ≈ sol_ip.t && sol_scalar.u ≈ sol_ip[1, :]
    end

    @testset "long time span" begin
      prob = DDEProblemLibrary.prob_dde_constant_1delay_long_scalar

      sol1 = solve(prob, alg1; abstol = 1e-12, reltol = 1e-12)
      sol2 = solve(prob, alg2; abstol = 1e-8, reltol = 1e-10)
      sol3 = solve(prob, alg3; abstol = 1e-8, reltol = 1e-10)
      sol4 = solve(prob, alg4; abstol = 1e-12, reltol = 1e-12)

      # relaxed tests to prevent floating point issues
      @test abs(sol1[end] - sol2[end]) < 2.4e-8
      @test abs(sol1[end] - sol3[end]) < 3.7e-8
      @test abs(sol1[end] - sol4[end]) < 9.0e-11 # 9.0e-13
    end
  end

  ## Two constant delays
  @testset "two constant delays" begin
    @testset "short time span" begin
      ### Scalar function
      sol_scalar = solve(DDEProblemLibrary.prob_dde_constant_2delays_scalar, alg)

      @test sol_scalar.errors[:l∞] < 2.4e-6
      @test sol_scalar.errors[:final] < 2.1e-6
      @test sol_scalar.errors[:l2] < 1.2e-6

      ### Out-of-place function
      sol_oop = solve(DDEProblemLibrary.prob_dde_constant_2delays_oop, alg)

      @test sol_scalar.t ≈ sol_oop.t && sol_scalar.u ≈ sol_oop[1, :]

      ### In-place function
      sol_ip = solve(DDEProblemLibrary.prob_dde_constant_2delays_ip, alg)

      @test sol_scalar.t ≈ sol_ip.t && sol_scalar.u ≈ sol_ip[1, :]
    end

    @testset "long time span" begin
      prob = DDEProblemLibrary.prob_dde_constant_2delays_long_scalar

      sol1 = solve(prob, alg1; abstol = 1e-12, reltol = 1e-12)
      sol2 = solve(prob, alg2; abstol = 1e-8, reltol = 1e-10)
      sol3 = solve(prob, alg3; abstol = 1e-8, reltol = 1e-10)
      sol4 = solve(prob, alg4; abstol = 1e-12, reltol = 1e-12)

      # relaxed tests to prevent floating point issues
      @test abs(sol1[end] - sol3[end]) < 1.2e-13 # 1.2e-15
      @test abs(sol1[end] - sol4[end]) < 3.1e-13 # 3.1e-15
    end
  end
end

## Non-standard history functions
@testset "non-standard history" begin
  alg = MethodOfSteps(Tsit5(); constrained = false,
                      fpsolve = NLFunctional(; max_iter = 100))

  @testset "idxs" begin
    function f(du,u,h,p,t)
      du[1] = -h(p, t-0.2;idxs=1) + u[1]
    end
    h(p, t; idxs = nothing) = idxs isa Number ? 0.0 : [0.0]

    prob = DDEProblem(f, [1.0], h, (0.0, 100.0); constant_lags = [0.2])
    solve(prob, alg; abstol = 1e-12, reltol = 1e-12)
  end

  @testset "in-place" begin
    function f(du,u,h,p,t)
      h(du, p, t-0.2)
      du[1] = -du[1] + u[1]
    end
    h(val, p, t) = (val .= 0.0)

    prob = DDEProblem(f, [1.0], h, (0.0, 100.0); constant_lags = [0.2])
    solve(prob, alg; abstol = 1e-12, reltol = 1e-12)
  end
end


text:
function fib05()
#
#  test polynomial suggested by Goedecker
#
    n = 5;
    p = fib(n)
    z = [-.6783507129699967-.4585361872731445*im   1;
         -.6783507129699967+.4585361872731445*im   1;
         .19537659464725405-.8488536405462456*im   1;
         .19537659464725405+.8488536405462456*im   1;
         1.9659482366454853                       1];

    p, PolyZeros(z)
end


text:
module ReduxImGui

Base.Experimental.@optlevel 1

using Redux
using CImGui

include("forward_decls.jl")

## Buttons
include("buttons/Buttons.jl")
using .Buttons
import .Buttons: Button
export Buttons

include("buttons/OnOffButtons.jl")
using .OnOffButtons
import .OnOffButtons: OnOffButton
export OnOffButtons

include("buttons/ColorButtons.jl")
using .ColorButtons
import .ColorButtons: ColorButton
export ColorButtons

include("buttons/RadioButtons.jl")
using .RadioButtons
import .RadioButtons: RadioButton
export RadioButtons

include("buttons/RadioButtonGroups.jl")
using .RadioButtonGroups
import .RadioButtonGroups: RadioButtonGroup
export RadioButtonGroups

include("buttons/ArrowButtons.jl")
using .ArrowButtons
import .ArrowButtons: ArrowButton
export ArrowButtons

## Checkboxes
include("checkboxes/Checkboxes.jl")
using .Checkboxes
import .Checkboxes: Checkbox
export Checkboxes

## Combos
include("combos/Combos.jl")
using .Combos
import .Combos: Combo
export Combos

## Inputs
include("inputs/InputTexts.jl")
using .InputTexts
import .InputTexts: InputText
export InputTexts

include("inputs/InputTextWithHints.jl")
using .InputTextWithHints
import .InputTextWithHints: InputTextWithHint
export InputTextWithHints

include("inputs/InputInts.jl")
using .InputInts
import .InputInts: InputInt
export InputInts

include("inputs/InputFloats.jl")
using .InputFloats
import .InputFloats: InputFloat
export InputFloats

include("inputs/InputDoubles.jl")
using .InputDoubles
import .InputDoubles: InputDouble
export InputDoubles

## Drags
include("drags/DragInts.jl")
using .DragInts
import .DragInts: DragInt
export DragInts

include("drags/DragFloats.jl")
using .DragFloats
import .DragFloats: DragFloat
export DragFloats

## Sliders
include("sliders/SliderInts.jl")
using .SliderInts
import .SliderInts: SliderInt
export SliderInts

include("sliders/SliderFloats.jl")
using .SliderFloats
import .SliderFloats: SliderFloat
export SliderFloats

include("sliders/SliderAngles.jl")
using .SliderAngles
import .SliderAngles: SliderAngle
export SliderAngles

include("sliders/SliderStrings.jl")
using .SliderStrings
import .SliderStrings: SliderString
export SliderStrings

## ColorEdits
include("color_edits/ColorEdits.jl")
using .ColorEdits
import .ColorEdits: ColorEdit
export ColorEdits

## ListBoxes
include("listboxes/ListBoxes.jl")
using .ListBoxes
import .ListBoxes: ListBox, get_items, get_item
export ListBoxes

include("menus/includes.jl")

## Others
include("trees/TreeNodes.jl")
using .TreeNodes
import .TreeNodes: TreeNode, TreeNodeEx
export TreeNodes

include("tooltips/Tooltips.jl")
using .Tooltips
import .Tooltips: Tooltip
export Tooltips

end # module


text:
# This file was generated, do not modify it. # hide
y1_oos = vcat(y11, y12, y13);


text:
"""
    getplotvectors(pp::PoissonProcess, intervallength::Real) -> (Tuple{Vector{Float64},
                                                                 Vector{Int}})

Simulate `pp` over the interval `[0, intervallength]`, and return the vectors needed to
plot.

We duplicate vector entries in the returned tuple in order to help us plot an orthogonal
polygonal chain.
"""
function getplotvectors(pp::PoissonProcess, intervallength::Real)
    # simulate
    arrivaltimes = rand(pp, intervallength)
    if length(arrivaltimes) == 0
        return Float64[], Int[]
    end
    # allocate memory
    paddedarrivals = Vector{Float64}(undef, 2*length(arrivaltimes))
    paddednarrivals = Vector{Int}(undef, 2*length(arrivaltimes))
    # duplicate vector entries for plotting purposes (we want an orthogonal polygonal chain)
    for (i, value) in enumerate(arrivaltimes)
        paddedarrivals[2*i-1], paddedarrivals[2*i] = value, value
    end
    # pattern: 0, 1, 1, 2, 2, ..., length(arrivaltimes)
    paddednarrivals[1] = 0
    for i in 1:(length(arrivaltimes)-1)
        paddednarrivals[2*i], paddednarrivals[2*i+1] = i, i
    end
    paddednarrivals[length(paddednarrivals)] = length(arrivaltimes)
    return paddedarrivals, paddednarrivals
end

"""
    plot(pp::PoissonProcess, intervallength::Real = 30.0, nrealizations::Integer = 10) -> (
    nothing)

Simulate and plot `nrealizations` realizations of a `PoissonProcess`.
"""
function plot(pp::PoissonProcess, intervallength::Real = 30.0, nrealizations::Integer = 10)
    mytitle = string(nrealizations, " realizations of a Poisson process with arrival rate ",
                     pp.rate)

    if nrealizations == 1
        mylw, myalpha, mytitle = 3, 0.9, string("1 realization of a Poisson process with ",
                                                "arrival rate ", pp.rate)
    elseif nrealizations <= 3
        mylw, myalpha = 3, 0.9
    elseif nrealizations <= 10
        mylw, myalpha = 3, 0.6
    elseif nrealizations <= 20
        mylw, myalpha = 2, 0.5
    else
        mylw, myalpha = 1, 0.5
    end

    xs, ys = getplotvectors(pp, intervallength)
    myplot = Plots.plot(xs, ys, title=mytitle, lw=mylw, alpha=myalpha, xlabel="Time",
                        leg=false)

    # display and return if only one realization is to be drawn
    if nrealizations == 1
        display(myplot)
        return nothing
    end
    # mutate myplot with additional realizations
    for step in 2:nrealizations
        xs, ys = getplotvectors(pp, intervallength)
        Plots.plot!(myplot, xs, ys, lw=mylw, alpha=myalpha, leg=false)
    end
    display(myplot)
    return nothing
end

max_stars_repo_name: LoganDGraham/StochasticProcesses

max_issues_repo_name: LoganDGraham/StochasticProcesses

max_forks_repo_name: LoganDGraham/StochasticProcesses


text:
export computeXtX, computeBeta, computeFits, estimateD, callPeaks, generatePeakFile, generateUnfusedPeakFile

#########################################################################################################
# Computes XtX for linear regression
#
# Example usage:
# out = computeXtX(MatrixReader("/scratch/hiranumn/forward.data100", 10000))
# out1 = computeXtX(MatrixReader("/scratch/hiranumn/reverse.data100", 10000))
# JLD.save("../data/xtxs.jld", "XtX1-f", out[1], "XtX2-f", out[2], "XtX1-r", out1[1], "XtX2-r", out1[2])
#########################################################################################################

function computeXtX(mr::MatrixReader; num_chroms=0, verbose=0, binsize=100)

    ##############################
    # Use all chroms for default #
    ##############################
    if num_chroms > length(ReferenceContigs_hg38.sizes) || num_chroms < 1
        num_chroms = length(ReferenceContigs_hg38.sizes)
    end
    training_limit = Int(ceil(sum(ReferenceContigs_hg38.sizes[1:num_chroms])/binsize))

    ################
    # Compute XtXs #
    ################
    XtX1 = zeros(Float64, mr.expsize+1, mr.expsize+1)
    XtX2 = zeros(Float64, mr.expsize+1, mr.expsize+1)

    advance!(mr)
    count = 0
    while !eof(mr) && count*mr.blocksize < training_limit

        # compute
        ctrl = convert(Array{Float64,2}, addConstColumn(mr.data)')
        if count%2==0
            BLAS.syrk!('U', 'N', 1.0, ctrl, 1.0, XtX1)
        else
            BLAS.syrk!('U', 'N', 1.0, ctrl, 1.0, XtX2)
        end

        # report progress
        if verbose>0 && count % (verbose) == 0
            println(count, ":", training_limit)
        end

        # update
        advance!(mr)
        count += 1
    end

    #############################
    # Converting to full matrix #
    #############################
    XtX1 = XtX1+XtX1'
    XtX2 = XtX2+XtX2'
    for i in 1:size(XtX1)[1]
       XtX1[i, i] ./= 2
       XtX2[i, i] ./= 2
    end

    XtX1, XtX2
end

#########################################################################################################
# Computes Beta for specific target
#
# Example usage:
# verbosity = 100
#
# _mr = MatrixReader("/scratch/hiranumn/forward.data100", 10000)
# _br = BinnedReader("/scratch/hiranumn/target_data_nodup/ENCFF000YRS.bam.fbin100")
#
# w_f = computeBeta(_mr, _br, "f", verbose=verbosity, xtxfile="../data/xtxs.jld")
#
# _mr = MatrixReader("/scratch/hiranumn/reverse.data100", 10000)
# _br = BinnedReader("/scratch/hiranumn/target_data_nodup/ENCFF000YRS.bam.rbin100")
#
# w_r = computeBeta(_mr, _br, "r", verbose=verbosity, xtxfile="../data/xtxs.jld")
#
# JLD.save("ENCFF000YRS.jld", "w1-f", w_f[1], "w2-f", w_f[2], "w1-r", w_r[1], "w2-r", w_r[2])
#########################################################################################################

function computeBeta(mr::MatrixReader, br::BinnedReader, direction::String; binsize=100, num_chroms=0, verbose=0, mask=[], xtxfile="../data/xtx.jld")

    ##############################
    # Use all chroms for default #
    ##############################
    if num_chroms > length(ReferenceContigs_hg38.sizes) || num_chroms < 1
        num_chroms = length(ReferenceContigs_hg38.sizes)
    end
    training_limit = Int(ceil(sum(ReferenceContigs_hg38.sizes[1:num_chroms])/binsize))

    ###########################################
    # Prepare matrices for weight calculation #
    ###########################################
    datapath = joinpath(@__DIR__, "..", "data")
    XtX1 = load(joinpath(datapath, xtxfile))["XtX1-$(direction)"]
    XtX2 = load(joinpath(datapath, xtxfile))["XtX2-$(direction)"]
    if length(mask)>0
        # mask input if necessary
        @assert mr.expsize+1==length(mask)
        XtX1 = filter2d(XtX1, mask, mask)
        XtX2 = filter2d(XtX2, mask, mask)
    end
    Xty1 = zeros((size(XtX1)[1],1))
    Xty2 = zeros((size(XtX2)[1],1))

    ###################
    # Compute Weights #
    ###################

    count = 0
    for target in denseblocks([br], mr.blocksize, constantColumn=false, loop=true)

        # update
        count += 1
        if count*mr.blocksize > training_limit break end

        # load control
        advance!(mr)
        ctrl = convert(Array{Float64,2}, addConstColumn(mr.data)')
        if length(mask)>0 ctrl = ctrl[mask, :] end

        # compute
        if count % 2 == 0
            Xty1 .+= ctrl*target
        else
            Xty2 .+= ctrl*target
        end

        # report progress
	if verbose>0 && count % (verbose) == 0
        	progress = Int(floor((count*mr.blocksize/training_limit)*1000))/10
        	printString = "$(progress)% completed ($(count*mr.blocksize)/$(training_limit))"
        	if direction=="f"
            		printString = printString*" on forward signals."
        	else
            		printString = printString*" on reverse signals."
        	end
        	println(printString)
    	end
    end

    m = size(XtX1)[1]
    beta1 = inv(XtX1 + 0.00001*Matrix(1.0I, m, m))*Xty1
    beta2 = inv(XtX2 + 0.00001*Matrix(1.0I, m, m))*Xty2

    beta1, beta2
end

###############################################################################
# Computes fits for specific target
#
# Example usage:
# verbosity = 100
#
# _mr = MatrixReader("/scratch/hiranumn/forward.data100", 10000)
# ff = computeFits(_mr, "ENCFF000YRS.jld", "f", verbose=verbosity)
#
# _mr = MatrixReader("/scratch/hiranumn/reverse.data100", 10000)
# fr = computeFits(_mr, "ENCFF000YRS.jld", "r", verbose=verbosity)
#
# JLD.save("ENCFF000YRS_fit.jld", "fit-f", ff, "fit-r", fr)
###############################################################################

function computeFits(mr::MatrixReader, weightfile::String, direction::String; binsize=100, num_chroms=0, verbose=0, mask=[])

    ##############################
    # Use all chroms for default #
    ##############################
    if num_chroms > length(ReferenceContigs_hg38.sizes) || num_chroms < 1
        num_chroms = length(ReferenceContigs_hg38.sizes)
    end
    training_limit = Int(ceil(sum(ReferenceContigs_hg38.sizes[1:num_chroms])/binsize))

    weight1 = load(weightfile)["w1-$(direction)"]
    weight2 = load(weightfile)["w2-$(direction)"]

    ##########################
    # Compute regression fit #
    ##########################
    regression_fit = zeros(Float16, training_limit)

    advance!(mr)

    count = 0
    binpos = 0
    while !eof(mr) && count*mr.blocksize < training_limit

        # get data
        ctrl = convert(Array{Float64,2}, addConstColumn(mr.data)')
        if length(mask)>0 ctrl = ctrl[mask, :] end

        # compute
        if count % 2 == 0
            pred = ctrl'*weight2
        else
            pred = ctrl'*weight1
        end

        # record
        for j in 1:length(pred)
            binpos +=1
            try
                regression_fit[binpos] = pred[j]
            catch
            end
        end

        # report progress
	if verbose>0 && (count+1) % (verbose) == 0
        	progress = Int(floor(((count+1)*mr.blocksize/training_limit)*1000))/10
        	printString = "$(progress)% completed ($((count+1)*mr.blocksize)/$(training_limit))"
        	if direction=="f"
            		printString = printString*" on forward signals."
        	else
            		printString = printString*" on reverse signals."
        	end
        	println(printString)
    	end

        advance!(mr)
        # update
        count += 1
    end

    regression_fit
end

############################################################################
# Estimates distance between forward and reverse reads for specific target #
############################################################################

function estimateD(forwardtarget, reversetarget; binsize=100)

    #########################
    # Load in target values #
    #########################
    b = BinnedReader(forwardtarget)
    targetf = zeros(Float32, Int(ceil(sum(ReferenceContigs_hg38.sizes)/binsize)))
    while !eof(b) && position(b) < length(targetf)
        targetf[position(b)] = value(b)
        advance!(b)
    end
    close(b)
    b = BinnedReader(reversetarget)
    targetr = zeros(Float32, Int(ceil(sum(ReferenceContigs_hg38.sizes)/binsize)))
    while !eof(b) && position(b) < length(targetr)
        targetr[position(b)] = value(b)
        advance!(b)
    end
    close(b)

    ##################################################
    # Figure out distance that gives minimum overlap #
    ##################################################
    d = []
    for i in 0:4
        f = vcat([0 for j in  1:i], targetf[1:end-i])
        @assert f!=targetr
        @assert length(f)==length(targetr)
        push!(d, sum(abs.(f-targetr)))
    end
    argmin(d)-1
end

# Calls peak for both forward and reverse strands
function callPeaks(br::BinnedReader, fitfile::String, direction::String; num_chroms=0, verbose=0, binsize=100, base=1, smoothing=true)

    if num_chroms > length(ReferenceContigs_hg38.sizes) || num_chroms < 1
        num_chroms = length(ReferenceContigs_hg38.sizes)
    end
    training_limit = Int(ceil(sum(ReferenceContigs_hg38.sizes[1:num_chroms])/binsize))

    # fill in target vector
    target = zeros(Float32, Int(ceil(sum(ReferenceContigs_hg38.sizes[1:num_chroms])/binsize)))
    while !eof(br) && position(br) < length(target)
        target[position(br)] = value(br)
        advance!(br)
    end
    close(br)
    regfit = load(fitfile)["fit-$(direction)"]

    if verbose>0 println("Loaded peak signals.") end

    # Do smoothing if necessary
    if smoothing
        smooth1000 = smooth(regfit, 10)
        smooth5000 = smooth(regfit, 50)
        smooth10000 = smooth(regfit, 100)
    end
    m = mean(target)

    #if verbose>0 println("smoothed.") end

    # Recording vector
    pvals = zeros(Float32, Int(ceil(sum(ReferenceContigs_hg38.sizes[1:num_chroms])/binsize)))
    folds = zeros(Float32, Int(ceil(sum(ReferenceContigs_hg38.sizes[1:num_chroms])/binsize)))
    lambdas = zeros(Float32, Int(ceil(sum(ReferenceContigs_hg38.sizes[1:num_chroms])/binsize)))

    # Compute p-values and folds
    for i in 1:length(regfit)
        if smoothing
            lambda = maximum([regfit[i], smooth1000[i], smooth5000[i], smooth10000[i], base])
        else
            lambda = maximum([regfit[i], base])
        end
        pval = -1*log(10, 1-cdf(Poisson(lambda), target[i]))
        fold = target[i]/lambda
        if pval == Inf
            pval = Float32(typemax(Int64))
        end

        #This version of the code will assign pval to individual bins.
        pvals[i] = pval
        folds[i] = fold
        lambdas[i] = lambda

	# report progress
	if verbose>0 && (i/10000) % (verbose) == 0
        	progress = Int(floor((i/training_limit)*1000))/10
        	printString = "$(progress)% completed ($(i)/$(training_limit))"
        	if direction=="f"
            		printString = printString*" on forward signals."
        	else
            		printString = printString*" on reverse signals."
        	end
        	println(printString)
    	end

    end

    pvals, folds, target, lambdas
end

# combines reverse and forward peaks
function generatePeakFile(pfile::String, name::String; th=1.5, binsize=100)

    # create a vector final folds and pvals
    pvals = zeros(Float32, Int(ceil(sum(ReferenceContigs_hg38.sizes)/binsize)))
    folds = zeros(Float32, Int(ceil(sum(ReferenceContigs_hg38.sizes)/binsize)))

    # get signals for each direction
    saved_data = load(pfile)
    forward_p = saved_data["p-f"]
    reverse_p = saved_data["p-r"]
    folds_f = saved_data["fold-f"]
    folds_r = saved_data["fold-r"]
    offset = saved_data["offset"]

    # Figuring out forward and reverse offset
    forward_offset = Int(ceil(offset/2))
    reverse_offset = Int(floor(offset/2))

    #Write into pvals and folds vectors from forward signals
    for i in 1:length(forward_p)
        try
            pvals[i+forward_offset] = forward_p[i]
            folds[i+forward_offset] = folds_f[i]
        catch
        end
    end

    #For reverse signals
    for i in 1:length(reverse_p)
        try
            if reverse_p[i] < pvals[i-reverse_offset]
                pvals[i-reverse_offset] = reverse_p[i]
            end
            if folds_r[i] < folds[i-reverse_offset]
                folds[i-reverse_offset] = folds_r[i]
            end
        catch
        end
    end

    #Get peaks to write on files
    processed_peaks = sortPeaks(pvals, folds, th)
    pw = PeakWriter(open("$(name).narrowPeak", "w"), ReferenceContigs_hg38)
    written = []
    for p in processed_peaks
        push!(written, writePeak(pw, 100, Int(p[1]), Int(p[2]), p[3], p[4]))
    end
    close(pw)
    written
end

function generateUnfusedPeakFile(pfile::String, name::String; th=1.5, binsize=100)

    # Data loading
    saved_data = load(pfile)
    forward_p = saved_data["p-f"]
    reverse_p = saved_data["p-r"]
    folds_f = saved_data["fold-f"]
    folds_r = saved_data["fold-r"]
    offset = saved_data["offset"]

    pvals = zeros(Float32, Int(ceil(sum(ReferenceContigs_hg38.sizes)/binsize)))
    folds = zeros(Float32, Int(ceil(sum(ReferenceContigs_hg38.sizes)/binsize)))

    forward_offset = Int(ceil(offset/2))
    reverse_offset = Int(floor(offset/2))

    for i in 1:length(forward_p)
        try
            pvals[i+forward_offset] = forward_p[i]
            folds[i+forward_offset] = folds_f[i]
        catch
        end
    end

    for i in 1:length(reverse_p)
        try
            if reverse_p[i] < pvals[i-reverse_offset]
                pvals[i-reverse_offset] = reverse_p[i]
            end
            if folds_r[i] < folds[i-reverse_offset]
                folds[i-reverse_offset] = folds_r[i]
            end
        catch
        end
    end

    fout = open("$(name).narrowPeak","w")
    pw = PeakWriter_unfused(fout, ReferenceContigs_hg38)
    for i in 1:length(pvals)
        #This version of the code will assign pval to individual bins.
        pval = pvals[i]
        if pval >= th
            WritePeak_unfused(pw, 100, i, i, pval, folds[i])
        end
    end
    close(fout)
    pvals
end