text:
# Run tests on Newton methods & rfi methods

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra

@testset "Newton iterations" begin
    nep=nep_gallery("dep0")

    @bench @testset "Newton and AugNewton" begin
        @info "Newton and AugNewton test"

        # newton and augnewton are equivalent, therefore I expect them
        # to generate identical results
        λ1,x1 = newton(nep,logger=displaylevel,v=ones(size(nep,1)),λ=0,tol=eps()*10)
        λ2,x2 = augnewton(nep,logger=displaylevel,v=ones(size(nep,1)),λ=0,tol=eps()*10)
        @test λ1 ≈ λ2
        @test x1 ≈ x2
        @test compute_resnorm(nep,λ1,x1) < eps()*100
        @test compute_resnorm(nep,λ2,x2) < eps()*100

    end

    @bench @testset "QuasiNewton" begin
    @info "QuasiNewton  test"

        λ,x = quasinewton(nep,logger=displaylevel,v=ones(size(nep,1)),λ=0,tol=1e-12)
        @test compute_resnorm(nep,λ,x) < 1e-11*100

    end

    @bench @testset "Newton QR" begin

        @info "Newton QR test"
        #Run derivative test for left and right eigenvectors returned by newtonqr
        λ3,x3,y3 =  newtonqr(nep, λ=0, v=ones(size(nep,1)), logger=displaylevel,tol=eps()*10)
        @test compute_resnorm(nep,λ3,x3) < eps()*100

        @info "  Testing formula for derivative"
        tau = 1
        Mλ = -I - tau * nep.A[2] * exp(-tau*λ3)
        Mtau = -λ3*nep.A[2]*exp(-tau*λ3)

        # Formula for derivative
        λp=-(y3'*(Mtau)*x3) / (y3'* Mλ*x3)

        δ=0.0001
        nepp = deepcopy(nep)
        nepp.tauv[2] = nep.tauv[2]+δ
        λ3δ,x3,y3 = newtonqr(nepp, λ=0, v=ones(size(nep,1)), logger=displaylevel,tol=eps()*10)
        λp_approx = (λ3δ-λ3)/δ
        @test abs(λp-λp_approx)< (δ*10)
    end

    @bench @testset "Resinv" begin
        @info "resinv + periodicdde"
        nep1=nep_gallery("periodicdde", name="mathieu")
        # basic functionality of resinv. Start close to solution to speed up unit test
        λ4,x4 =  resinv(nep1, λ=-0.2447, v=[0.970208+0.0im, -0.242272+0.0im],
                        logger=displaylevel,tol=eps()*10)
        r4=compute_resnorm(nep1,λ4,x4)

        @test r4 < eps()*100
    end

    @bench @testset "Rayleigh Function Iteration" begin
        @info "rfi"
        nept = DEP([copy(nep.A[1]'), copy(nep.A[2]')],copy(nep.tauv))

        n=size(nep,1)
        u0=ones(n)
        v0=ones(n)

        λ,x,y =rfi(nep, nept, logger=displaylevel, v=ones(n), u=ones(n), tol=1e-15)
        @test compute_resnorm(nep,λ,x) < eps()*100
        @test compute_resnorm(nept,λ,y) < eps()*100

        # Test RFIb

        @info "rfi_b"
        λb,xb,yb =rfi_b(nep, nept, logger=displaylevel, v=v0, u=u0, λ=λ+0.01, tol=1e-15)
        @test λ ≈ λb

        @info "  Testing formula for derivative (with left and right eigvecs)"
        tau = 1
        Mλ = -I - tau*nep.A[2] * exp(-tau*λ)
        Mtau = -λ*nep.A[2]*exp(-tau*λ)

        # Formula for derivative
        λp=-(y'*(Mtau)*x) / (y'* Mλ*x)

        δ=0.0001
        nepp = deepcopy(nep)
        nepp.tauv[2] = nep.tauv[2]+δ
        neptp = deepcopy(nept)
        neptp.tauv[2] = nept.tauv[2]+δ
        λδ,x,y = rfi(nepp,neptp,logger=displaylevel, v=ones(n), u=ones(n))
        λp_approx=(λδ-λ)/δ
        @test abs(λp-λp_approx)< (δ*10)
    end

    @bench @testset "implicitdet" begin
        @info "Implicitdet test"
        nepd=nep_gallery("periodicdde", name="mathieu")
        λ,v=implicitdet(nepd, v=ones(size(nepd,1)), logger=displaylevel)
        @test norm(compute_Mder(nepd,λ)*v) < eps()*100
    end

    @bench @testset "type stability" begin
        dep = nep_gallery("dep0")
        @inferred augnewton(Float64, dep, v = ones(size(dep,1)))
        @inferred augnewton(ComplexF64, dep, v = ones(size(dep,1)))
    end
end


text:
abstract type AbstractSimulator end

# ============================================================================================

"""
    Simulator
# Fields:
    P::StateSpace
    f = (x,p,t) -> x
    y = (x,t)   -> y
"""
struct Simulator{S<:AbstractStateSpace,F<:Function,Y<:Function} <: AbstractSimulator
    P::S
    f::F
    y::Y
end

"""
    Simulator(P::StateSpace, u = (x,t) -> 0)

Used to simulate continuous-time systems. See function `?solve` for additional info.
# Usage:
```
using OrdinaryDiffEq
dt             = 0.1
tfinal         = 20
t              = 0:dt:tfinal
P              = ss(tf(1,[2,1])^2)
K              = 5
reference(x,t) = [1.]
s              = Simulator(P, reference)
x0             = [0.,0]
tspan          = (0.0,tfinal)
sol            = solve(s, x0, tspan, Tsit5())
plot(t, s.y(sol, t)[:], lab="Open loop step response")
```
"""
function Simulator(P::AbstractStateSpace, u::F = (x,t) -> 0) where F
    iscontinuous(P) || throw(ArgumentError("Simulator only supports continuous-time system. See function `lsim` for simulation of discrete-time systems."))
    f = (dx,x,p,t) -> dx .= P.A*x .+ P.B*u(x,t)
    y(x,t) = P.C*x .+ P.D*u(x,t)
    y(sol::ODESolution,t) = P.C*sol(t) .+ P.D*u(sol(t),t)
    Simulator(P, f, y)
end

"""
    sol = solve(s::AbstractSimulator, x0, tspan,  args...; kwargs...)
Simulate the system represented by `s` from initial state `x0` over time span `tspan = (t0,tf)`.
`args` and `kwargs` are sent to the `solve` function from `OrdinaryDiffEq`, e.g., `solve(s, x0, tspan,  Tsit5(), reltol=1e-5)` solves the problem with solver [`Tsit5()`](http://docs.juliadiffeq.org/stable/solvers/ode_solve.html) and relative tolerance 1e-5.

See also `Simulator` `lsim`
"""
DiffEqBase.solve(s::AbstractSimulator, x0, tspan, solver=Tsit5(), args...; kwargs...) = solve(ODEProblem(s.f,x0,tspan), solver, args...; kwargs...)

ext: jl


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

@doc raw"""ManagedFieldsEntry is a workflow-id, a FieldSet and the group version of the resource that the fieldset applies to.

    IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry(;
        apiVersion=nothing,
        fieldsType=nothing,
        fieldsV1=nothing,
        manager=nothing,
        operation=nothing,
        time=nothing,
    )

    - apiVersion::String : APIVersion defines the version of this resource that this field set applies to. The format is \"group/version\" just like the top-level APIVersion field. It is necessary to track the version of a field set because it cannot be automatically converted.
    - fieldsType::String : FieldsType is the discriminator for the different fields format and version. There is currently only one possible value: \"FieldsV1\"
    - fieldsV1::IoK8sApimachineryPkgApisMetaV1FieldsV1 : FieldsV1 holds the first JSON version format as described in the \"FieldsV1\" type.
    - manager::String : Manager is an identifier of the workflow managing these fields.
    - operation::String : Operation is the type of operation which lead to this ManagedFieldsEntry being created. The only valid values for this field are 'Apply' and 'Update'.
    - time::IoK8sApimachineryPkgApisMetaV1Time : Time is timestamp of when these fields were set. It should always be empty if Operation is 'Apply'
"""
mutable struct IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    fieldsType::Any # spec type: Union{ Nothing, String } # spec name: fieldsType
    fieldsV1::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1FieldsV1 } # spec name: fieldsV1
    manager::Any # spec type: Union{ Nothing, String } # spec name: manager
    operation::Any # spec type: Union{ Nothing, String } # spec name: operation
    time::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: time

    function IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry(;apiVersion=nothing, fieldsType=nothing, fieldsV1=nothing, manager=nothing, operation=nothing, time=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry, Symbol("fieldsType"), fieldsType)
        setfield!(o, Symbol("fieldsType"), fieldsType)
        validate_property(IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry, Symbol("fieldsV1"), fieldsV1)
        setfield!(o, Symbol("fieldsV1"), fieldsV1)
        validate_property(IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry, Symbol("manager"), manager)
        setfield!(o, Symbol("manager"), manager)
        validate_property(IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry, Symbol("operation"), operation)
        setfield!(o, Symbol("operation"), operation)
        validate_property(IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry, Symbol("time"), time)
        setfield!(o, Symbol("time"), time)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry

const _property_map_IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("fieldsType")=>Symbol("fieldsType"), Symbol("fieldsV1")=>Symbol("fieldsV1"), Symbol("manager")=>Symbol("manager"), Symbol("operation")=>Symbol("operation"), Symbol("time")=>Symbol("time"))
const _property_types_IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("fieldsType")=>"String", Symbol("fieldsV1")=>"IoK8sApimachineryPkgApisMetaV1FieldsV1", Symbol("manager")=>"String", Symbol("operation")=>"String", Symbol("time")=>"IoK8sApimachineryPkgApisMetaV1Time")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry)
    true
end

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


text:
"""Module containing functions for retrieving substances' data from a database."""
module FunctionsDB

using DataFrames
using DBInterface
using SQLite

export substance_extract, texture_extract, air_extract, cn_extract, cn_list_extract

const dbi = DBInterface
const sql = SQLite

"""
    substance_extract( numCAS::String, fields::Vector{String}, db_path::String = "..\\..\\..\\resources\\Analysis data\\substance.db" )

Extract values that describe a pollutant from the columns indicated by `fields` of the database at `db_path` where the CAS number matches `numCAS`.
"""
function substance_extract( numCAS::String, fields::Vector{String}, db_path::String="" )
    if  isempty(db_path)
        db_path = occursin("src", @__DIR__) ? split(@__DIR__, "src")[1] : *(@__DIR__, "\\..\\")
        db_path *= "resources\\Analysis data\\substance.db"
    end
    # estrazione valori sostanze
    db = sql.DB(db_path)
    sql_fields = join(fields, ", ")
    query_substance = sql.Stmt(db, "SELECT "*sql_fields*" FROM substance WHERE n_CAS LIKE ?")
    results = dbi.execute(query_substance, [numCAS])
    resdf = DataFrame(results)
    return resdf
end

"""
    texture_extract( texture_name::String, fields::Vector{String}, db_path::String = "..\\..\\..\\resources\\Analysis data\\substance.db" )

Extract values that describe a tipe of texture from the columns indicated by `fields` of the database at `db_path` where the column `name` matches the value of `textre_name`.
"""
function texture_extract( texture_name::String, fields::Vector{String}, db_path::String="" )
    if  isempty(db_path)
        db_path = occursin("src", @__DIR__) ? split(@__DIR__, "src")[1] : *(@__DIR__, "\\..\\")
        db_path *= "resources\\Analysis data\\substance.db"
    end
    # estrazione valori sostanze
    db = sql.DB(db_path)
    sql_fields = join(fields, ", ")
    query_texture = sql.Stmt(db, "SELECT "*sql_fields*" FROM texture WHERE nome LIKE ?" )
    results = dbi.execute(query_texture, [texture_name])
    resdf = DataFrame(results)
    return resdf
end

"""
    air_extract( stability_class::String, outdoor_class::String, fields::Vector{String}, db_path::String="..\\..\\..\\resources\\Analysis data\\substance.db" )

Extract values from the columns indicated by `fields` of the database at `db_path` where the outdoor and stability class match `outdoor_class` and `stability_class`.
"""
function air_extract( stability_class::String, outdoor_class::String, fields::Vector{String}, db_path::String="" )
    if  isempty(db_path) || isnothing(db_path)
        db_path = occursin("src", @__DIR__) ? split(@__DIR__, "src")[1] : *(@__DIR__, "\\..\\")
        db_path *= "resources\\Analysis data\\substance.db"
    end
    db = sql.DB(db_path)
    sql_fields = join(fields, ", ")
    query_texture = sql.Stmt(db, "SELECT "*sql_fields*" FROM air_stability WHERE class LIKE ? AND outdoor LIKE ?")
    results = dbi.execute(query_texture, [stability_class, outdoor_class])
    resdf = DataFrame(results)
    return resdf
end

"""
    cn_extract( cnl::String, id_soil::Int64, dbloc::String="..\\..\\..\\resources\\Analysis data\\substance.db" )

Extract the value of class `cnl` of a substance from the database at `db_path` where the id matches `soil`.
"""
function cn_extract( cnl::String, id_soil::Int64, dbloc::String="" )
    if  isempty(db_path) || isnothing(db_path)
        db_path = occursin("src", @__DIR__) ? split(@__DIR__, "src")[1] : *(@__DIR__, "\\..\\")
        db_path *= "resources\\Analysis data\\substance.db"
    end
    db = sql.DB(dbloc*"substance.db")
    classecn = "cn_$cnl"
    query_cn = sql.Stmt( db, "SELECT "*classecn*" FROM cn WHERE id = ?" )
    results = dbi.execute(query_cn, [id_soil])
    resdf = DataFrame(results)
    return resdf
end

"""
    cn_list_extract( db_path::String="..\\..\\..\\resources\\Analysis data\\substance.db" )

Extract all values from table `cn` of the database at `db_path`.
"""
function cn_list_extract( db_path::String="" )
    if  isempty(db_path) || isnothing(db_path)
        db_path = occursin("src", @__DIR__) ? split(@__DIR__, "src")[1] : *(@__DIR__, "\\..\\")
        db_path *= "resources\\Analysis data\\substance.db"
    end
	db = sql.DB(db_path)
    query_cn = sql.Stmt( db, "SELECT * FROM cn" )
    results = dbi.execute(query_cn)
    resdf = DataFrame(results)
    return resdf
end

end # module

#=Tables
"family":
    Chemical family of substances, referenced by table "substance"
"volatility":
    Classes of volatility, referenced by table "substance"
"substance":
    Chemical substances, with various information on each
"air_stability":
    Stability classes, outdoor classes, various sigma values
"sqlite_sequence":
    Not relevant
"texture":
    Kind of terrain textures, with informations on them like permeability and such
"sqlite_stat1"
    Not relevant
"cn":
    Related to runoff (cn -> runoff curve number), also has clc (-> CORINE Land Cover ?)
=#
#= CHECK A TABLE TO FIND THE VALUE TO BE USED FOR THE CONDITION
db_path = occursin("src", @__DIR__) ? split(@__DIR__, "src")[1] : *(@__DIR__, "\\..\\")
db_path *= "resources\\Analysis data\\substance.db"
db = sql.DB(db_path)
query = sql.Stmt( db, "SELECT * FROM substance" )
results = dbi.execute(query)
resdf = DataFrame(results)

query2 = sql.Stmt(db, "SELECT rfd_ing, rfd_inal, rfc FROM substance WHERE n_CAS LIKE ?")
result2 = dbi.execute(query2, ["16065-83-1"])
resdf2 = DataFrame(result2)
=#

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
using ModelingToolkit, StaticArrays, LinearAlgebra
using StochasticDiffEq, SparseArrays
using Random, Test

# Define some variables
@parameters t σ ρ β
@variables x(t) y(t) z(t)
D = Differential(t)

eqs = [D(x) ~ σ * (y - x),
    D(y) ~ x * (ρ - z) - y,
    D(z) ~ x * y - β * z]

noiseeqs = [0.1 * x,
    0.1 * y,
    0.1 * z]

# ODESystem -> SDESystem shorthand constructor
@named sys = ODESystem(eqs, t, [x, y, z], [σ, ρ, β])
@test SDESystem(sys, noiseeqs, name = :foo) isa SDESystem

@named de = SDESystem(eqs, noiseeqs, t, [x, y, z], [σ, ρ, β])
f = eval(generate_diffusion_function(de)[1])
@test f(ones(3), rand(3), nothing) == 0.1ones(3)

f = SDEFunction(de)
prob = SDEProblem(SDEFunction(de), f.g, [1.0, 0.0, 0.0], (0.0, 100.0), (10.0, 26.0, 2.33))
sol = solve(prob, SRIW1(), seed = 1)

probexpr = SDEProblem(SDEFunction(de), f.g, [1.0, 0.0, 0.0], (0.0, 100.0),
                      (10.0, 26.0, 2.33))
solexpr = solve(eval(probexpr), SRIW1(), seed = 1)

@test all(x -> x == 0, Array(sol - solexpr))

# Test no error
@test_nowarn SDEProblem(de, nothing, (0, 10.0))

noiseeqs_nd = [0.01*x 0.01*x*y 0.02*x*z
               σ 0.01*y 0.02*x*z
               ρ β 0.01*z]
@named de = SDESystem(eqs, noiseeqs_nd, t, [x, y, z], [σ, ρ, β])
f = eval(generate_diffusion_function(de)[1])
@test f([1, 2, 3.0], [0.1, 0.2, 0.3], nothing) == [0.01*1 0.01*1*2 0.02*1*3
       0.1 0.01*2 0.02*1*3
       0.2 0.3 0.01*3]

f = eval(generate_diffusion_function(de)[2])
du = ones(3, 3)
f(du, [1, 2, 3.0], [0.1, 0.2, 0.3], nothing)
@test du == [0.01*1 0.01*1*2 0.02*1*3
             0.1 0.01*2 0.02*1*3
             0.2 0.3 0.01*3]

f = SDEFunction(de)
prob = SDEProblem(SDEFunction(de), f.g, [1.0, 0.0, 0.0], (0.0, 100.0), (10.0, 26.0, 2.33),
                  noise_rate_prototype = zeros(3, 3))
sol = solve(prob, EM(), dt = 0.001)

u0map = [
    x => 1.0,
    y => 0.0,
    z => 0.0,
]

parammap = [
    σ => 10.0,
    β => 26.0,
    ρ => 2.33,
]

prob = SDEProblem(de, u0map, (0.0, 100.0), parammap)
@test size(prob.noise_rate_prototype) == (3, 3)
@test prob.noise_rate_prototype isa Matrix
sol = solve(prob, EM(), dt = 0.001)

prob = SDEProblem(de, u0map, (0.0, 100.0), parammap, sparsenoise = true)
@test size(prob.noise_rate_prototype) == (3, 3)
@test prob.noise_rate_prototype isa SparseMatrixCSC
sol = solve(prob, EM(), dt = 0.001)

# Test eval_expression=false
function test_SDEFunction_no_eval()
    # Need to test within a function scope to trigger world age issues
    f = SDEFunction(de, eval_expression = false)
    @test f([1.0, 0.0, 0.0], (10.0, 26.0, 2.33), (0.0, 100.0)) ≈ [-10.0, 26.0, 0.0]
end
test_SDEFunction_no_eval()

# modelingtoolkitize and Ito <-> Stratonovich sense
seed = 10
Random.seed!(seed)

# simple 2D diagonal noise
u0 = rand(2)
t = randn()
trange = (0.0, 100.0)
p = [1.01, 0.87]
f1!(du, u, p, t) = (du .= p[1] * u)
σ1!(du, u, p, t) = (du .= p[2] * u)

prob = SDEProblem(f1!, σ1!, u0, trange, p)
# no correction
sys = modelingtoolkitize(prob)
fdrift = eval(generate_function(sys)[1])
fdif = eval(generate_diffusion_function(sys)[1])
@test fdrift(u0, p, t) == p[1] * u0
@test fdif(u0, p, t) == p[2] * u0
fdrift! = eval(generate_function(sys)[2])
fdif! = eval(generate_diffusion_function(sys)[2])
du = similar(u0)
fdrift!(du, u0, p, t)
@test du == p[1] * u0
fdif!(du, u0, p, t)
@test du == p[2] * u0

# Ito -> Strat
sys2 = stochastic_integral_transform(sys, -1 // 2)
fdrift = eval(generate_function(sys2)[1])
fdif = eval(generate_diffusion_function(sys2)[1])
@test fdrift(u0, p, t) == p[1] * u0 - 1 // 2 * p[2]^2 * u0
@test fdif(u0, p, t) == p[2] * u0
fdrift! = eval(generate_function(sys2)[2])
fdif! = eval(generate_diffusion_function(sys2)[2])
du = similar(u0)
fdrift!(du, u0, p, t)
@test du == p[1] * u0 - 1 // 2 * p[2]^2 * u0
fdif!(du, u0, p, t)
@test du == p[2] * u0

# Strat -> Ito
sys2 = stochastic_integral_transform(sys, 1 // 2)
fdrift = eval(generate_function(sys2)[1])
fdif = eval(generate_diffusion_function(sys2)[1])
@test fdrift(u0, p, t) == p[1] * u0 + 1 // 2 * p[2]^2 * u0
@test fdif(u0, p, t) == p[2] * u0
fdrift! = eval(generate_function(sys2)[2])
fdif! = eval(generate_diffusion_function(sys2)[2])
du = similar(u0)
fdrift!(du, u0, p, t)
@test du == p[1] * u0 + 1 // 2 * p[2]^2 * u0
fdif!(du, u0, p, t)
@test du == p[2] * u0

# somewhat complicated 1D without explicit parameters but with explicit time-dependence
f2!(du, u, p, t) = (du[1] = sin(t) + cos(u[1]))
σ2!(du, u, p, t) = (du[1] = pi + atan(u[1]))

u0 = rand(1)
prob = SDEProblem(f2!, σ2!, u0, trange)
# no correction
sys = modelingtoolkitize(prob)
fdrift = eval(generate_function(sys)[1])
fdif = eval(generate_diffusion_function(sys)[1])
@test fdrift(u0, p, t) == @. sin(t) + cos(u0)
@test fdif(u0, p, t) == pi .+ atan.(u0)
fdrift! = eval(generate_function(sys)[2])
fdif! = eval(generate_diffusion_function(sys)[2])
du = similar(u0)
fdrift!(du, u0, p, t)
@test du == @. sin(t) + cos(u0)
fdif!(du, u0, p, t)
@test du == pi .+ atan.(u0)

# Ito -> Strat
sys2 = stochastic_integral_transform(sys, -1 // 2)
fdrift = eval(generate_function(sys2)[1])
fdif = eval(generate_diffusion_function(sys2)[1])
@test fdrift(u0, p, t) == @. sin(t) + cos(u0) - 1 // 2 * 1 / (1 + u0^2) * (pi + atan(u0))
@test fdif(u0, p, t) == pi .+ atan.(u0)
fdrift! = eval(generate_function(sys2)[2])
fdif! = eval(generate_diffusion_function(sys2)[2])
du = similar(u0)
fdrift!(du, u0, p, t)
@test du == @. sin(t) + cos(u0) - 1 // 2 * 1 / (1 + u0^2) * (pi + atan(u0))
fdif!(du, u0, p, t)
@test du == pi .+ atan.(u0)

# Strat -> Ito
sys2 = stochastic_integral_transform(sys, 1 // 2)
fdrift = eval(generate_function(sys2)[1])
fdif = eval(generate_diffusion_function(sys2)[1])
@test fdrift(u0, p, t) ≈ @. sin(t) + cos(u0) + 1 // 2 * 1 / (1 + u0^2) * (pi + atan(u0))
@test fdif(u0, p, t) == pi .+ atan.(u0)
fdrift! = eval(generate_function(sys2)[2])
fdif! = eval(generate_diffusion_function(sys2)[2])
du = similar(u0)
fdrift!(du, u0, p, t)
@test du ≈ @. sin(t) + cos(u0) + 1 // 2 * 1 / (1 + u0^2) * (pi + atan(u0))
fdif!(du, u0, p, t)
@test du == pi .+ atan.(u0)

# 2D diagonal noise with mixing terms (no parameters, no time-dependence)
u0 = rand(2)
t = randn()
function f3!(du, u, p, t)
    du[1] = u[1] / 2
    du[2] = u[2] / 2
    return nothing
end
function σ3!(du, u, p, t)
    du[1] = u[2]
    du[2] = u[1]
    return nothing
end

prob = SDEProblem(f3!, σ3!, u0, trange, p)
# no correction
sys = modelingtoolkitize(prob)
fdrift = eval(generate_function(sys)[1])
fdif = eval(generate_diffusion_function(sys)[1])
@test fdrift(u0, p, t) == u0 / 2
@test fdif(u0, p, t) == reverse(u0)
fdrift! = eval(generate_function(sys)[2])
fdif! = eval(generate_diffusion_function(sys)[2])
du = similar(u0)
fdrift!(du, u0, p, t)
@test du == u0 / 2
fdif!(du, u0, p, t)
@test du == reverse(u0)

# Ito -> Strat
sys2 = stochastic_integral_transform(sys, -1 // 2)
fdrift = eval(generate_function(sys2)[1])
fdif = eval(generate_diffusion_function(sys2)[1])
@test fdrift(u0, p, t) == u0 * 0
@test fdif(u0, p, t) == reverse(u0)
fdrift! = eval(generate_function(sys2)[2])
fdif! = eval(generate_diffusion_function(sys2)[2])
du = similar(u0)
fdrift!(du, u0, p, t)
@test du == u0 * 0
fdif!(du, u0, p, t)
@test du == reverse(u0)

# Strat -> Ito
sys2 = stochastic_integral_transform(sys, 1 // 2)
fdrift = eval(generate_function(sys2)[1])
fdif = eval(generate_diffusion_function(sys2)[1])
@test fdrift(u0, p, t) == u0
@test fdif(u0, p, t) == reverse(u0)
fdrift! = eval(generate_function(sys2)[2])
fdif! = eval(generate_diffusion_function(sys2)[2])
du = similar(u0)
fdrift!(du, u0, p, t)
@test du == u0
fdif!(du, u0, p, t)
@test du == reverse(u0)

# simple 2D diagonal noise oop
u0 = rand(2)
t = randn()
p = [1.01, 0.87]
f1(u, p, t) = p[1] * u
σ1(u, p, t) = p[2] * u

prob = SDEProblem(f1, σ1, u0, trange, p)
# no correction
sys = modelingtoolkitize(prob)
fdrift = eval(generate_function(sys)[1])
fdif = eval(generate_diffusion_function(sys)[1])
@test fdrift(u0, p, t) == p[1] * u0
@test fdif(u0, p, t) == p[2] * u0
fdrift! = eval(generate_function(sys)[2])
fdif! = eval(generate_diffusion_function(sys)[2])
du = similar(u0)
fdrift!(du, u0, p, t)
@test du == p[1] * u0
fdif!(du, u0, p, t)
@test du == p[2] * u0

# Ito -> Strat
sys2 = stochastic_integral_transform(sys, -1 // 2)
fdrift = eval(generate_function(sys2)[1])
fdif = eval(generate_diffusion_function(sys2)[1])
@test fdrift(u0, p, t) == p[1] * u0 - 1 // 2 * p[2]^2 * u0
@test fdif(u0, p, t) == p[2] * u0
fdrift! = eval(generate_function(sys2)[2])
fdif! = eval(generate_diffusion_function(sys2)[2])
du = similar(u0)
fdrift!(du, u0, p, t)
@test du == p[1] * u0 - 1 // 2 * p[2]^2 * u0
fdif!(du, u0, p, t)
@test du == p[2] * u0

# Strat -> Ito
sys2 = stochastic_integral_transform(sys, 1 // 2)
fdrift = eval(generate_function(sys2)[1])
fdif = eval(generate_diffusion_function(sys2)[1])
@test fdrift(u0, p, t) == p[1] * u0 + 1 // 2 * p[2]^2 * u0
@test fdif(u0, p, t) == p[2] * u0
fdrift! = eval(generate_function(sys2)[2])
fdif! = eval(generate_diffusion_function(sys2)[2])
du = similar(u0)
fdrift!(du, u0, p, t)
@test du == p[1] * u0 + 1 // 2 * p[2]^2 * u0
fdif!(du, u0, p, t)
@test du == p[2] * u0

# non-diagonal noise
u0 = rand(2)
t = randn()
p = [1.01, 0.3, 0.6, 1.2, 0.2]
f4!(du, u, p, t) = du .= p[1] * u
function g4!(du, u, p, t)
    du[1, 1] = p[2] * u[1]
    du[1, 2] = p[3] * u[1]
    du[2, 1] = p[4] * u[1]
    du[2, 2] = p[5] * u[2]
    return nothing
end

prob = SDEProblem(f4!, g4!, u0, trange, noise_rate_prototype = zeros(2, 2), p)
# no correction
sys = modelingtoolkitize(prob)
fdrift = eval(generate_function(sys)[1])
fdif = eval(generate_diffusion_function(sys)[1])
@test fdrift(u0, p, t) == p[1] * u0
@test fdif(u0, p, t) == [p[2]*u0[1] p[3]*u0[1]
                         p[4]*u0[1] p[5]*u0[2]]
fdrift! = eval(generate_function(sys)[2])
fdif! = eval(generate_diffusion_function(sys)[2])
du = similar(u0)
fdrift!(du, u0, p, t)
@test du == p[1] * u0
du = similar(u0, size(prob.noise_rate_prototype))
fdif!(du, u0, p, t)
@test du == [p[2]*u0[1] p[3]*u0[1]
             p[4]*u0[1] p[5]*u0[2]]

# Ito -> Strat
sys2 = stochastic_integral_transform(sys, -1 // 2)
fdrift = eval(generate_function(sys2)[1])
fdif = eval(generate_diffusion_function(sys2)[1])
@test fdrift(u0, p, t) ≈ [
    p[1] * u0[1] - 1 // 2 * (p[2]^2 * u0[1] + p[3]^2 * u0[1]),
    p[1] * u0[2] - 1 // 2 * (p[2] * p[4] * u0[1] + p[5]^2 * u0[2]),
]
@test fdif(u0, p, t) == [p[2]*u0[1] p[3]*u0[1]
                         p[4]*u0[1] p[5]*u0[2]]
fdrift! = eval(generate_function(sys2)[2])
fdif! = eval(generate_diffusion_function(sys2)[2])
du = similar(u0)
fdrift!(du, u0, p, t)
@test du ≈ [
    p[1] * u0[1] - 1 // 2 * (p[2]^2 * u0[1] + p[3]^2 * u0[1]),
    p[1] * u0[2] - 1 // 2 * (p[2] * p[4] * u0[1] + p[5]^2 * u0[2]),
]
du = similar(u0, size(prob.noise_rate_prototype))
fdif!(du, u0, p, t)
@test du == [p[2]*u0[1] p[3]*u0[1]
             p[4]*u0[1] p[5]*u0[2]]

# Strat -> Ito
sys2 = stochastic_integral_transform(sys, 1 // 2)
fdrift = eval(generate_function(sys2)[1])
fdif = eval(generate_diffusion_function(sys2)[1])
@test fdrift(u0, p, t) ≈ [
    p[1] * u0[1] + 1 // 2 * (p[2]^2 * u0[1] + p[3]^2 * u0[1]),
    p[1] * u0[2] + 1 // 2 * (p[2] * p[4] * u0[1] + p[5]^2 * u0[2]),
]
@test fdif(u0, p, t) == [p[2]*u0[1] p[3]*u0[1]
                         p[4]*u0[1] p[5]*u0[2]]
fdrift! = eval(generate_function(sys2)[2])
fdif! = eval(generate_diffusion_function(sys2)[2])
du = similar(u0)
fdrift!(du, u0, p, t)
@test du ≈ [
    p[1] * u0[1] + 1 // 2 * (p[2]^2 * u0[1] + p[3]^2 * u0[1]),
    p[1] * u0[2] + 1 // 2 * (p[2] * p[4] * u0[1] + p[5]^2 * u0[2]),
]
du = similar(u0, size(prob.noise_rate_prototype))
fdif!(du, u0, p, t)
@test du == [p[2]*u0[1] p[3]*u0[1]
             p[4]*u0[1] p[5]*u0[2]]

# non-diagonal noise: Torus -- Strat and Ito are identical
u0 = rand(2)
t = randn()
p = rand(1)
f5!(du, u, p, t) = du .= false
function g5!(du, u, p, t)
    du[1, 1] = cos(p[1]) * sin(u[1])
    du[1, 2] = cos(p[1]) * cos(u[1])
    du[1, 3] = -sin(p[1]) * sin(u[2])
    du[1, 4] = -sin(p[1]) * cos(u[2])
    du[2, 1] = sin(p[1]) * sin(u[1])
    du[2, 2] = sin(p[1]) * cos(u[1])
    du[2, 3] = cos(p[1]) * sin(u[2])
    du[2, 4] = cos(p[1]) * cos(u[2])
    return nothing
end

prob = SDEProblem(f5!, g5!, u0, trange, noise_rate_prototype = zeros(2, 4), p)
# no correction
sys = modelingtoolkitize(prob)
fdrift = eval(generate_function(sys)[1])
fdif = eval(generate_diffusion_function(sys)[1])
@test fdrift(u0, p, t) == 0 * u0
@test fdif(u0, p, t) ==
      [cos(p[1])*sin(u0[1]) cos(p[1])*cos(u0[1]) -sin(p[1])*sin(u0[2]) -sin(p[1])*cos(u0[2])
       sin(p[1])*sin(u0[1]) sin(p[1])*cos(u0[1]) cos(p[1])*sin(u0[2]) cos(p[1])*cos(u0[2])]
fdrift! = eval(generate_function(sys)[2])
fdif! = eval(generate_diffusion_function(sys)[2])
du = similar(u0)
fdrift!(du, u0, p, t)
@test du == 0 * u0
du = similar(u0, size(prob.noise_rate_prototype))
fdif!(du, u0, p, t)
@test du ==
      [cos(p[1])*sin(u0[1]) cos(p[1])*cos(u0[1]) -sin(p[1])*sin(u0[2]) -sin(p[1])*cos(u0[2])
       sin(p[1])*sin(u0[1]) sin(p[1])*cos(u0[1]) cos(p[1])*sin(u0[2]) cos(p[1])*cos(u0[2])]

# Ito -> Strat
sys2 = stochastic_integral_transform(sys, -1 // 2)
fdrift = eval(generate_function(sys2)[1])
fdif = eval(generate_diffusion_function(sys2)[1])
@test fdrift(u0, p, t) == 0 * u0
@test fdif(u0, p, t) ==
      [cos(p[1])*sin(u0[1]) cos(p[1])*cos(u0[1]) -sin(p[1])*sin(u0[2]) -sin(p[1])*cos(u0[2])
       sin(p[1])*sin(u0[1]) sin(p[1])*cos(u0[1]) cos(p[1])*sin(u0[2]) cos(p[1])*cos(u0[2])]
fdrift! = eval(generate_function(sys2)[2])
fdif! = eval(generate_diffusion_function(sys2)[2])
du = similar(u0)
fdrift!(du, u0, p, t)
@test du == 0 * u0
du = similar(u0, size(prob.noise_rate_prototype))
fdif!(du, u0, p, t)
@test du ==
      [cos(p[1])*sin(u0[1]) cos(p[1])*cos(u0[1]) -sin(p[1])*sin(u0[2]) -sin(p[1])*cos(u0[2])
       sin(p[1])*sin(u0[1]) sin(p[1])*cos(u0[1]) cos(p[1])*sin(u0[2]) cos(p[1])*cos(u0[2])]

# Strat -> Ito
sys2 = stochastic_integral_transform(sys, 1 // 2)
fdrift = eval(generate_function(sys2)[1])
fdif = eval(generate_diffusion_function(sys2)[1])
@test fdrift(u0, p, t) == 0 * u0
@test fdif(u0, p, t) ==
      [cos(p[1])*sin(u0[1]) cos(p[1])*cos(u0[1]) -sin(p[1])*sin(u0[2]) -sin(p[1])*cos(u0[2])
       sin(p[1])*sin(u0[1]) sin(p[1])*cos(u0[1]) cos(p[1])*sin(u0[2]) cos(p[1])*cos(u0[2])]
fdrift! = eval(generate_function(sys2)[2])
fdif! = eval(generate_diffusion_function(sys2)[2])
du = similar(u0)
fdrift!(du, u0, p, t)
@test du == 0 * u0
du = similar(u0, size(prob.noise_rate_prototype))
fdif!(du, u0, p, t)
@test du ==
      [cos(p[1])*sin(u0[1]) cos(p[1])*cos(u0[1]) -sin(p[1])*sin(u0[2]) -sin(p[1])*cos(u0[2])
       sin(p[1])*sin(u0[1]) sin(p[1])*cos(u0[1]) cos(p[1])*sin(u0[2]) cos(p[1])*cos(u0[2])]

# issue #819
@testset "Combined system name collisions" begin
    @variables t
    eqs_short = [D(x) ~ σ * (y - x),
        D(y) ~ x * (ρ - z) - y,
    ]
    sys1 = SDESystem(eqs_short, noiseeqs, t, [x, y, z], [σ, ρ, β], name = :sys1)
    sys2 = SDESystem(eqs_short, noiseeqs, t, [x, y, z], [σ, ρ, β], name = :sys1)
    @test_throws ArgumentError SDESystem([sys2.y ~ sys1.z], t, [], [], [],
                                         systems = [sys1, sys2], name = :foo)
end

# observed variable handling
@variables t x(t) RHS(t)
@parameters τ
D = Differential(t)
@named fol = SDESystem([D(x) ~ (1 - x) / τ], [x], t, [x], [τ];
                       observed = [RHS ~ (1 - x) / τ])
@test isequal(RHS, @nonamespace fol.RHS)
RHS2 = RHS
@unpack RHS = fol
@test isequal(RHS, RHS2)


text:
module SimpleNeuralNets

    using Reexport
    @reexport using Distances

    import Base: size, show, split, shuffle! #, slice
    export NN, DataSet,
           size, depth, inputsize, nsamples,
           split, shuffle!,
           calc,
           Sigmoid, Tanh, Relu, SoftRelu,
           sgd

    include("nodefunc.jl")
    include("net.jl")
    include("state.jl")
    include("metrics.jl")
    include("dataset.jl")
    include("solvers.jl")

end # module

max_stars_repo_name: fredo-dedup/SimpleNeuralNets

max_issues_repo_name: fredo-dedup/SimpleNeuralNets

max_forks_repo_name: fredo-dedup/SimpleNeuralNets


text:
using Test, Logging, Random
using DecisionProgramming

@info "Testing ChandeNode"
@test isa(ChanceNode(1, Node[]), ChanceNode)
@test isa(ChanceNode(2, Node[1]), ChanceNode)
@test_throws DomainError ChanceNode(1, Node[1])
@test_throws DomainError ChanceNode(1, Node[2])

@info "Testing DecisionNode"
@test isa(DecisionNode(1, Node[]), DecisionNode)
@test isa(DecisionNode(2, Node[1]), DecisionNode)
@test_throws DomainError DecisionNode(1, Node[1])
@test_throws DomainError DecisionNode(1, Node[2])

@info "Testing ValueNode"
@test isa(ValueNode(1, Node[]), ValueNode)
@test isa(ValueNode(2, Node[1]), ValueNode)
@test_throws DomainError ValueNode(1, Node[1])
@test_throws DomainError ValueNode(1, Node[2])

@info "Testing State"
@test isa(States([1, 2, 3]), States)
@test_throws DomainError States([0, 1])
@test States([(2, [1, 3]), (3, [2, 4, 5])]) == States([2, 3, 2, 3, 3])

@info "Testing validate_influence_diagram"
@test_throws DomainError validate_influence_diagram(
    States([1]),
    [ChanceNode(1, Node[])],
    [DecisionNode(2, Node[])],
    [ValueNode(3, [1, 2])]
)
@test_throws DomainError validate_influence_diagram(
    States([1, 1]),
    [ChanceNode(1, Node[])],
    [DecisionNode(1, Node[])],
    [ValueNode(3, [1, 2])]
)
@test_throws DomainError validate_influence_diagram(
    States([1, 1]),
    [ChanceNode(1, Node[])],
    [DecisionNode(2, Node[])],
    [ValueNode(2, [1])]
)
@test_throws DomainError validate_influence_diagram(
    States([1, 1]),
    [ChanceNode(1, Node[])],
    [DecisionNode(2, Node[])],
    [ValueNode(3, [2]), ValueNode(4, [3])]
)
# Test redundancy
@test validate_influence_diagram(
    States([1, 1]),
    [ChanceNode(1, Node[])],
    [DecisionNode(2, Node[])],
    [ValueNode(3, Node[])]
) === nothing

@info "Testing paths"
@test vec(collect(paths(States([2, 3])))) == [(1, 1), (2, 1), (1, 2), (2, 2), (1, 3), (2, 3)]
@test vec(collect(paths(States([2, 3]), Dict(1=>2)))) == [(2, 1), (2, 2), (2, 3)]

@info "Testing Probabilities"
@test isa(Probabilities(1, [0.4 0.6; 0.3 0.7]), Probabilities)
@test isa(Probabilities(1, [0.0, 0.4, 0.6]), Probabilities)
@test_throws DomainError Probabilities(1, [1.1, 0.1])

@info "Testing DefaultPathProbability"
P = DefaultPathProbability(
    [ChanceNode(1, Node[]), ChanceNode(2, [1])],
    [Probabilities(1, [0.4, 0.6]), Probabilities(2, [0.3 0.7; 0.9 0.1])]
)
@test isa(P, DefaultPathProbability)
@test P((1, 2)) == 0.4 * 0.7

@info "Testing Consequences"
@test isa(Consequences(1, [-1.1, 0.0, 2.7]), Consequences)
@test isa(Consequences(1, [-1.1 0.0; 2.7 7.0]), Consequences)

@info "Testing DefaultPathUtility"
U = DefaultPathUtility(
    [ValueNode(3, [2]), ValueNode(4, [1, 2])],
    [Consequences(3, [1.0, 1.4]), Consequences(4, [1.0 1.5; 0.6 3.4])]
)
@test isa(U, DefaultPathUtility)
@test U((2, 1)) == 1.0 + 0.6

@info "Testing LocalDecisionStrategy"
@test_throws DomainError LocalDecisionStrategy(1, [0, 0, 2])
@test_throws DomainError LocalDecisionStrategy(1, [0, 1, 1])


text:
#--- continuum ---#
prim = [1.0, 0.0, 1.0]
KitBase.prim_conserve(prim, 3.0)
KitBase.prim_conserve(prim[1], prim[2], prim[3], 3.0)

mprim = hcat(prim, prim)
KitBase.mixture_prim_conserve(mprim, 3.0)

KitBase.conserve_prim(prim[1])
KitBase.conserve_prim(prim[1], 1.0)
KitBase.conserve_prim(prim, 3.0)
KitBase.conserve_prim(prim[1], prim[2], prim[3], 3.0)
KitBase.mixture_conserve_prim(mprim, 3.0)

# Rykov
KitBase.prim_conserve([1.0, 0.0, 1.0, 1.0, 1.0], 5 / 3, 2)
KitBase.prim_conserve([1.0, 0.0, 0.0, 1.0, 1.0, 1.0], 5 / 3, 2)
KitBase.conserve_prim([1.0, 0.0, 1.0, 0.1], 5 / 3, 2)
KitBase.conserve_prim([1.0, 0.0, 0.0, 1.0, 0.1], 5 / 3, 2)

prim = [1.0, 0.2, 0.3, -0.1, 1.0]
mprim = hcat(prim, prim)
KitBase.em_coefficients(mprim, randn(3), randn(3), 100, 0.01, 0.01, 0.001)

KitBase.advection_flux(1.0, -0.1)
KitBase.burgers_flux(1.0)
KitBase.euler_flux(prim, 3.0)
KitBase.euler_jacobi(prim, 3.0)

#--- thermo ---#
KitBase.heat_capacity_ratio(2.0, 1)
KitBase.heat_capacity_ratio(2.0, 2)
KitBase.heat_capacity_ratio(2.0, 3)
KitBase.heat_capacity_ratio(2.0, 2, 1)
KitBase.heat_capacity_ratio(2.0, 2, 2)
KitBase.heat_capacity_ratio(2.0, 2, 3)

KitBase.sound_speed(1.0, 5 / 3)
KitBase.sound_speed([1.0, 0.0, 1.0], 5 / 3)
KitBase.sound_speed(rand(3, 2), 5 / 3)

#--- atom ---#
KitBase.pdf_slope(1.0, 0.1)
KitBase.pdf_slope(prim, randn(5), 0.0)
KitBase.mixture_pdf_slope(mprim, randn(5, 2), 0.0)

u = collect(-5:0.2:5)
ω = ones(51) .* 0.2
prim = [1.0, 0.0, 1.0]

M = KitBase.maxwellian(u, prim)
KitBase.maxwellian(randn(16, 16), randn(16, 16), [1.0, 0.0, 0.0, 1.0])
KitBase.maxwellian(
    randn(8, 8, 8),
    randn(8, 8, 8),
    randn(8, 8, 8),
    [1.0, 0.0, 0.0, 0.0, 1.0],
)

KitBase.energy_maxwellian(M, prim, 2)

KitBase.maxwellian!(rand(16), rand(16), prim)
KitBase.maxwellian!(rand(16, 16), randn(16, 16), randn(16, 16), [1.0, 0.0, 0.0, 1.0])
KitBase.maxwellian!(
    rand(8, 8, 8),
    randn(8, 8, 8),
    randn(8, 8, 8),
    randn(8, 8, 8),
    [1.0, 0.0, 0.0, 0.0, 1.0],
)

mprim = hcat(prim, prim)
KitBase.mixture_maxwellian(hcat(u, u), mprim)
KitBase.mixture_maxwellian(randn(8, 8, 2), randn(8, 8, 2), rand(4, 2))
KitBase.mixture_maxwellian(
    randn(8, 8, 8, 2),
    randn(8, 8, 8, 2),
    randn(8, 8, 8, 2),
    rand(5, 2),
)

KitBase.mixture_maxwellian!(randn(8, 2), randn(8, 2), rand(3, 2))
KitBase.mixture_maxwellian!(randn(8, 8, 2), randn(8, 8, 2), randn(8, 8, 2), rand(4, 2))
KitBase.mixture_maxwellian!(
    randn(8, 8, 8, 2),
    randn(8, 8, 8, 2),
    randn(8, 8, 8, 2),
    randn(8, 8, 8, 2),
    rand(5, 2),
)

KitBase.shakhov(u, M, 0.01, prim, 1.0)
KitBase.shakhov(u, M, M, 0.01, prim, 1.0, 2.0)
KitBase.shakhov(randn(16, 16), randn(16, 16), rand(16, 16), rand(2), [1.0, 0.0, 1.0], 1.0)
KitBase.shakhov(
    randn(16, 16),
    randn(16, 16),
    rand(16, 16),
    rand(16, 16),
    rand(2),
    [1.0, 0.0, 0.0, 1.0],
    1.0,
    1.0,
)
KitBase.shakhov(
    randn(8, 8, 8),
    randn(8, 8, 8),
    rand(8, 8, 8),
    rand(8, 8, 8),
    rand(3),
    [1.0, 0.0, 0.0, 0.0, 1.0],
    1.0,
)

KitBase.shakhov!(randn(16), randn(16), rand(16), 0.01, prim, 1.0)
KitBase.shakhov!(randn(16), randn(16), randn(16), rand(16), rand(16), 0.01, prim, 1.0, 2.0)
KitBase.shakhov!(
    randn(16, 16),
    randn(16, 16),
    randn(16, 16),
    rand(16, 16),
    rand(2),
    [1.0, 0.0, 1.0],
    1.0,
)
KitBase.shakhov!(
    randn(16, 16),
    randn(16, 16),
    randn(16, 16),
    randn(16, 16),
    rand(16, 16),
    rand(16, 16),
    rand(2),
    [1.0, 0.0, 0.0, 1.0],
    1.0,
    1.0,
)
KitBase.shakhov!(
    randn(8, 8, 8),
    randn(8, 8, 8),
    randn(8, 8, 8),
    rand(8, 8, 8),
    rand(8, 8, 8),
    rand(3),
    [1.0, 0.0, 0.0, 0.0, 1.0],
    1.0,
)

# Rykov
KitBase.maxwellian!(
    zeros(16),
    zeros(16),
    zeros(16),
    zeros(16),
    zeros(16),
    zeros(16),
    randn(16),
    rand(5),
    4,
    2,
)
KitBase.maxwellian!(
    zeros(8, 8),
    zeros(8, 8),
    zeros(8, 8),
    zeros(8, 8),
    zeros(8, 8),
    zeros(8, 8),
    randn(8, 8),
    randn(8, 8),
    rand(6),
    4,
    2,
)

KitBase.rykov!(
    zeros(16),
    zeros(16),
    zeros(16),
    zeros(16),
    zeros(16),
    zeros(16),
    randn(16),
    rand(16),
    rand(16),
    rand(16),
    rand(16),
    rand(16),
    rand(16),
    rand(2),
    rand(5),
    0.72,
    4,
    1 / 1.55,
    0.2354,
    0.3049,
)
KitBase.rykov!(
    zeros(8, 8),
    zeros(8, 8),
    zeros(8, 8),
    zeros(8, 8),
    zeros(8, 8),
    zeros(8, 8),
    randn(8, 8),
    randn(8, 8),
    rand(8, 8),
    rand(8, 8),
    rand(8, 8),
    rand(8, 8),
    rand(8, 8),
    rand(8, 8),
    rand(4),
    rand(6),
    0.72,
    4,
    1 / 1.55,
    0.2354,
    0.3049,
)

KitBase.reduce_distribution(randn(16, 51), ω, 1)
KitBase.reduce_distribution(randn(16, 24, 24), ones(24, 24), 1)
KitBase.reduce_distribution(
    randn(16, 24, 24),
    randn(16, 24, 24),
    randn(16, 24, 24),
    ones(24, 24),
    1,
)
KitBase.full_distribution(M, M, u, ω, ones(51, 24, 24), ones(51, 24, 24), 1.0, 3.0)
KitBase.full_distribution(M, M, u, ω, ones(51, 24, 24), ones(51, 24, 24), prim, 3.0)

KitBase.ref_vhs_vis(1.0, 1.0, 0.5)
KitBase.vhs_collision_time(prim[1], prim[end], 1e-3, 0.81)
KitBase.vhs_collision_time(prim, 1e-3, 0.81)

KitBase.rykov_zr(100, 91.5, 18.1)

KitBase.hs_boltz_kn(1e-3, 1.0)
vs = VSpace3D(-5, 5, 16, -5, 5, 16, -5, 5, 16)
fsm = KitBase.fsm_kernel(vs, 1e-3)
phi, psi, phipsi =
    KitBase.kernel_mode(5, 5.0, 5.0, 5.0, 0.1, 0.1, 0.1, 16, 16, 16, 1.0, quad_num = 16)
KitBase.kernel_mode(5, 5.0, 5.0, 5.0, 16, 16, 16, 1.0, quad_num = 16)
KitBase.kernel_mode(5, 5.0, 5.0, 0.1, 0.1, 16, 16, quad_num = 16)
KitBase.boltzmann_fft(rand(16, 16, 16), fsm)
KitBase.boltzmann_fft!(rand(16, 16, 16), rand(16, 16, 16), fsm)

KitBase.boltzmann_ode!(zeros(16, 16, 16), rand(16, 16, 16), (1.0, 5, phi, psi, phipsi), 0.0)
KitBase.bgk_ode!(zeros(16, 16, 16), rand(16, 16, 16), (rand(16, 16, 16), 1e-2), 0.0)

vs = KitBase.VSpace3D(-5.0, 5.0, 16, -5.0, 5.0, 16, -5.0, 5.0, 16, "algebra")
u, v, w = vs.u[:, 1, 1], vs.v[1, :, 1], vs.w[1, 1, :]
vnu = hcat(vs.u[:], vs.v[:], vs.w[:])
uuni1d = linspace(vs.u[1, 1, 1], vs.u[end, 1, 1], 16)
vuni1d = linspace(vs.v[1, 1, 1], vs.v[1, end, 1], 16)
wuni1d = linspace(vs.w[1, 1, 1], vs.w[1, 1, end], 16)
u13d = [uuni1d[i] for i = 1:16, j = 1:16, k = 1:16]
v13d = [vuni1d[j] for i = 1:16, j = 1:16, k = 1:16]
w13d = [wuni1d[k] for i = 1:16, j = 1:16, k = 1:16]
vuni = hcat(u13d[:], v13d[:], w13d[:])
KitBase.boltzmann_nuode!(
    zeros(16, 16, 16),
    rand(16, 16, 16),
    (5.0, 5, phi, psi, phipsi, u, v, w, vnu, uuni1d, vuni1d, wuni1d, vuni),
    0.0,
)

τ = KitBase.aap_hs_collision_time(mprim, 1.0, 0.5, 0.5, 0.5, 1.0)
KitBase.aap_hs_prim(mprim, τ, 1.0, 0.5, 0.5, 0.5, 1.0)
KitBase.aap_hs_prim(rand(4, 2), rand(2), 1.0, 0.5, 0.5, 0.5, 1e-2)
KitBase.aap_hs_prim(rand(5, 2), rand(2), 1.0, 0.5, 0.5, 0.5, 1e-2)

KitBase.aap_hs_diffeq!(
    similar(mprim),
    mprim,
    [τ[1], τ[2], 1.0, 0.5, 0.5, 0.5, 1.0, 3.0],
    0.0,
)
KitBase.shift_pdf!(M, 1.0, 1e-4, 1e-4)
KitBase.shift_pdf!(rand(16, 2), randn(2), rand(2), 1e-4)

KitBase.chapman_enskog(rand(16), [1.0, 0.0, 1.0], rand(3), rand(3), 0.1)
KitBase.chapman_enskog(rand(16), [1.0, 0.0, 1.0], zeros(3), 0, 0.1)
KitBase.chapman_enskog(
    rand(16, 16),
    rand(16, 16),
    [1.0, 0.0, 0.0, 1.0],
    rand(4),
    rand(4),
    rand(4),
    0.1,
)
KitBase.chapman_enskog(
    rand(16, 16),
    rand(16, 16),
    [1.0, 0.0, 0.0, 1.0],
    zeros(4),
    zeros(4),
    0.0,
    0.1,
)
KitBase.chapman_enskog(
    rand(16, 16, 16),
    rand(16, 16, 16),
    rand(16, 16, 16),
    [1.0, 0.0, 0.0, 0.0, 1.0],
    rand(5),
    rand(5),
    rand(5),
    rand(5),
    0.1,
)
KitBase.chapman_enskog(
    rand(16, 16, 16),
    rand(16, 16, 16),
    rand(16, 16, 16),
    [1.0, 0.0, 0.0, 0.0, 1.0],
    zeros(5),
    zeros(5),
    zeros(5),
    0.0,
    0.1,
)


text:
using QuasiArrays, DomainSets, StaticArrays, Test

@testset "Kron" begin
    @testset "InclusionKron" begin
        a = Inclusion(0:0.5:2)
        b = Inclusion([1,3,4])
        c = quasikron(a,b)
        @test eltype(c) == SVector{2,Float64}
        @test c[SVector(0.5,3)] == SVector(0.5,3)
        @test_throws BoundsError c[SVector(0.6,3)]
        @test_throws BoundsError c[SVector(0.5,2)]
    end

    @testset "QuasiKron" begin
        a = QuasiVector(randn(5), 0:0.5:2)
        b = QuasiVector([5,6,8], [1,3,4])
        K = quasikron(a, b)
        @test axes(K) == (quasikron(axes(a,1),axes(b,1)),)
        @test K[SVector(0.5,3)] == a[0.5]b[3]
        @test_throws BoundsError K[SVector(0.6,3)]
        @test_throws BoundsError K[SVector(0.5,2)]
    end
end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*--------------------------------------------------------------------
#* Model of cstr reactor
#*--------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id$
#*--------------------------------------------------------------------
type cstr
	cstr()=begin
		new(
			cstr_basic(),
			fill(reaction_mol(Dict{Symbol,Any}(
				:Brief=>"Molar Reaction Rate"
			)),(NReac)),
			fill(heat_reaction(Dict{Symbol,Any}(
				:Brief=>"Heat Reaction"
			)),(NReac)),
			[
				:(diff(_base_1.Outlet.z*M) = (_base_1.Inlet.F*_base_1.Inlet.z - _base_1.Outlet.F*_base_1.Outlet.z) + sumt(_base_1.stoic*r)*_base_1.Vr),
				:(diff(M*_base_1.Outlet.h) = _base_1.Inlet.F*_base_1.Inlet.h - _base_1.Outlet.F*_base_1.Outlet.h +sum(Hr*sum(_base_1.stoic*r))*_base_1.Vr - q),
			],
			[
				"Component Molar Balance","Reactor Energy Balance",
			],
			[:r,:Hr,]
		)
	end
	_base_1::cstr_basic
	r::Array{reaction_mol}
	Hr::Array{heat_reaction}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export cstr
function setEquationFlow(in::cstr)
	addEquation(1)
	addEquation(2)
	#FIXME sum(sum())

end
function atributes(in::cstr,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/cstr"
	fields[:Brief]="Model of a generic CSTR"
	fields[:Info]="
Requires the information of:
* Reaction values
* Heat of reaction
"
	drive!(fields,_)
	return fields
end
cstr(_::Dict{Symbol,Any})=begin
	newModel=cstr()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
#==============================================================================================#
#                            functions.jl -- TypeTree functionality                            #
#==============================================================================================#

using InteractiveUtils: subtypes

"""
`function tt(TY::Type, pref=("",); uni=true, concrete=true, mod=nothing)`\n
Returns a `Vector{String}` of type entries formatted as a tree for human (programmer)
inspection.

- The `TY` argument is the type for which to produce the tree;
- The `pref` optional argument is a text prefix;
- The `uni` optional keyword argument controls whether Unicode characters will be used;
- The `concrete` optional keyword argument controls whether to filter off concrete subtypes.
- The `mod` optional keyword argument removes module prefixes from type names.
"""
function tt(TY::Type, pref=("",); uni=true, concrete=true, mod=nothing)
    ret = Array{String,1}()
    ELL, FRK, BAR = uni ? (" └─ ", " ├─ ", " │  ") : (" \\-- ", " +-- ", " |  ")
    SPC = "    "
    LEN = length(pref) - 1
    PREF = LEN > 0 ? ( i == ELL ? SPC : i == FRK ? BAR : i for i in pref[1:LEN]) : ("",)
    append!(ret, ["$(join((PREF..., pref[end])))$(typename(TY, mod))\n"])
    ST = [i for i in subtypes(TY) if isabstracttype(i) || concrete]
    LE = length(ST)
    for i in 1:LE
        append!(
            ret, tt(ST[i], (pref..., i < LE ? FRK : ELL); uni, concrete, mod)
        )
    end
    ret
end

"""
    typename(SomeModule.SomeType) -> "SomeModule.SomeType"
    typename(SomeModule.SomeType, SomeModule) -> "SomeType"

Returns the name of SomeType in the context of SomeModule.
"""
function typename(TY::Type, mod)
    io = IOBuffer()
    show(mod == nothing ? io : IOContext(io, :module => mod), TY)
    String(take!(io))
end


text:
const IN_PLACE = MPI_IN_PLACE

"""
    Barrier(comm::Comm)

Blocks until `comm` is synchronized.

If `comm` is an intracommunicator, then it blocks until all members of the group have called it.

If `comm` is an intercommunicator, then it blocks until all members of the other group have called it.
"""
function Barrier(comm::Comm)
    # int MPI_Barrier(MPI_Comm comm)
    @mpichk ccall((:MPI_Barrier, libmpi), Cint, (MPI_Comm,), comm)
end

"""
    Bcast!(buf[, count=length(buf)], root, comm::Comm)

Broadcast the first `count` elements of the buffer `buf` from `root` to all processes.
"""
function Bcast!(buffer::MPIBuffertype{T}, count::Integer,
                root::Integer, comm::Comm) where T

    # int MPI_Bcast(void* buffer, int count, MPI_Datatype datatype, int root,
    #               MPI_Comm comm)
    @mpichk ccall((:MPI_Bcast, libmpi), Cint,
                  (Ptr{T}, Cint, MPI_Datatype, Cint, MPI_Comm),
                  buffer, count, mpitype(T), root, comm)
    buffer
end

function Bcast!(buffer::Array{T}, root::Integer, comm::Comm) where T
    Bcast!(buffer, length(buffer), root, comm)
end

function Bcast!(buffer::SubArray{T}, root::Integer, comm::Comm) where T
    @assert Base.iscontiguous(buffer)
    Bcast!(buffer, length(buffer), root, comm)
end

#=
function Bcast{T}(obj::T, root::Integer, comm::Comm)
    buf = [T]
    Bcast!(buf, root, comm)
    buf[1]
end
=#

function bcast(obj, root::Integer, comm::Comm)
    isroot = Comm_rank(comm) == root
    count = Array{Cint}(undef, 1)
    if isroot
        buf = MPI.serialize(obj)
        count[1] = length(buf)
    end
    Bcast!(count, root, comm)
    if !isroot
        buf = Array{UInt8}(undef, count[1])
    end
    Bcast!(buf, root, comm)
    if !isroot
        obj = MPI.deserialize(buf)
    end
    obj
end

"""
    Reduce!(sendbuf, recvbuf[, count=length(sendbuf)], op, root, comm)

Performs `op` reduction on the first `count` elements of the  buffer `sendbuf`
and stores the result in `recvbuf` on the process of rank `root`.

On non-root processes `recvbuf` is ignored.

To perform the reduction in place, see [`Reduce_in_place!`](@ref).

To handle allocation of the output buffer, see [`Reduce`](@ref).
"""
function Reduce!(sendbuf::MPIBuffertype{T}, recvbuf::MPIBuffertype{T},
                 count::Integer, op::Union{Op,MPI_Op}, root::Integer,
                 comm::Comm) where T
    isroot = Comm_rank(comm) == root
    isroot && typeof(recvbuf) <: AbstractArray && @assert length(recvbuf) >= count
    # int MPI_Reduce(const void* sendbuf, void* recvbuf, int count,
    #                MPI_Datatype datatype, MPI_Op op, int root, MPI_Comm comm)
    @mpichk ccall((:MPI_Reduce, libmpi), Cint,
                  (Ptr{T}, Ptr{T}, Cint, MPI_Datatype, MPI_Op, Cint, MPI_Comm),
                  sendbuf, recvbuf, count, mpitype(T), op, root, comm)
    isroot ? recvbuf : nothing
end

# Convert user-provided functions to MPI.Op
Reduce!(sendbuf::MPIBuffertype{T}, recvbuf::MPIBuffertype{T},
        count::Integer, opfunc, root::Integer, comm::Comm) where {T} =
    Reduce!(sendbuf, recvbuf, count, Op(opfunc, T), root, comm)

function Reduce!(sendbuf::MPIBuffertype{T}, recvbuf::MPIBuffertype{T},
                 op, root::Integer, comm::Comm) where T
    Reduce!(sendbuf, recvbuf, length(sendbuf), op, root, comm)
end

"""
    Reduce(sendbuf, count, op, root, comm)

Performs `op` reduction on the buffer `sendbuf` and stores the result in
an output buffer allocated on the process of rank `root`. An empty array
will be returned on all other processes.

To specify the output buffer, see [`Reduce!`](@ref).

To perform the reduction in place, see [`Reduce_in_place!`](@ref).
"""
function Reduce(sendbuf::MPIBuffertype{T}, count::Integer,
                op, root::Integer, comm::Comm) where T
    isroot = Comm_rank(comm) == root
    recvbuf = Array{T}(undef, isroot ? count : 0)
    Reduce!(sendbuf, recvbuf, count, op, root, comm)
end

function Reduce(sendbuf::Array{T,N}, op,
    root::Integer, comm::Comm) where {T,N}
    isroot = Comm_rank(comm) == root
    recvbuf = Array{T,N}(undef, isroot ? size(sendbuf) : Tuple(zeros(Int, ndims(sendbuf))))
    Reduce!(sendbuf, recvbuf, length(sendbuf), op, root, comm)
end

function Reduce(sendbuf::SubArray{T}, op, root::Integer, comm::Comm) where T
    @assert Base.iscontiguous(sendbuf)
    Reduce(sendbuf, length(sendbuf), op, root, comm)
end

function Reduce(object::T, op
                , root::Integer, comm::Comm) where T
    isroot = Comm_rank(comm) == root
    sendbuf = T[object]
    recvbuf = Reduce(sendbuf, op, root, comm)
    isroot ? recvbuf[1] : nothing
end

"""
    Reduce_in_place!(buf, count, op, root, comm)

Performs `op` reduction on the first `count` elements of the buffer `buf` and
stores the result on `buf` of the `root` process in the group.

This is equivalent to calling
```julia
if root == MPI.Comm_rank(comm)
    Reduce!(MPI.IN_PLACE, buf, count, root, comm)
else
    Reduce!(buf, C_NULL, count, root, comm)
end
```

To handle allocation of the output buffer, see [`Reduce`](@ref).

To specify a separate output buffer, see [`Reduce!`](@ref).
"""
function Reduce_in_place!(buf::MPIBuffertype{T}, count::Integer,
                          op::Union{Op,MPI_Op}, root::Integer,
                          comm::Comm) where T
    if Comm_rank(comm) == root
        typeof(buf) <: AbstractArray && @assert length(buf) >= count
        @mpichk ccall((:MPI_Reduce, libmpi), Cint,
                      (Ptr{T}, Ptr{T}, Cint, MPI_Datatype, MPI_Op, Cint, MPI_Comm),
                      MPI_IN_PLACE, buf, count, mpitype(T), op, root, comm)
    else
        @mpichk ccall((:MPI_Reduce, libmpi), Cint,
                      (Ptr{T}, Ptr{T}, Cint, MPI_Datatype, MPI_Op, Cint, MPI_Comm),
                      buf, C_NULL, count, mpitype(T), op, root, comm)
    end
    buf
end

# Convert to MPI.Op
Reduce_in_place!(buf::MPIBuffertype{T}, count::Integer, op,
                 root::Integer, comm::Comm) where T =
                    Reduce_in_place!(buf, count, Op(op,T), root, comm)

"""
    Allreduce!(sendbuf, recvbuf[, count=length(sendbuf)], op, comm)

Performs `op` reduction on the first `count` elements of the buffer
`sendbuf` storing the result in the `recvbuf` of all processes in the
group.

All-reduce is equivalent to a [`Reduce!`](@ref) operation followed by
a [`Bcast!`](@ref), but can lead to better performance.

If `sendbuf==MPI.IN_PLACE` the data is read from `recvbuf` and then overwritten
with the results.

To handle allocation of the output buffer, see [`Allreduce`](@ref).
"""
function Allreduce!(sendbuf::MPIBuffertypeOrConst{T}, recvbuf::MPIBuffertype{T},
                   count::Integer, op::Union{Op,MPI_Op}, comm::Comm) where T
    typeof(recvbuf) <: AbstractArray && @assert length(recvbuf) >= count
    # int MPI_Allreduce(const void* sendbuf, void* recvbuf, int count,
    #                   MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
    @mpichk ccall((:MPI_Allreduce, libmpi), Cint,
                  (Ptr{T}, Ptr{T}, Cint, MPI_Datatype, MPI_Op, MPI_Comm),
                  sendbuf, recvbuf, count, mpitype(T), op, comm)
    recvbuf
end

# Convert user-provided functions to MPI.Op
Allreduce!(sendbuf::MPIBuffertypeOrConst{T}, recvbuf::MPIBuffertype{T},
           count::Integer, opfunc, comm::Comm) where {T} =
    Allreduce!(sendbuf, recvbuf, count, Op(opfunc,T), comm)

function Allreduce!(sendbuf::MPIBuffertypeOrConst{T}, recvbuf::MPIBuffertype{T},
                   op, comm::Comm) where T
    Allreduce!(sendbuf, recvbuf, length(recvbuf), op, comm)
end

"""
    Allreduce!(buf, op, comm)

Performs `op` reduction in place on the buffer `sendbuf`, overwriting it with
the results on all the processes in the group.

Equivalent to calling `Allreduce!(MPI.IN_PLACE, buf, op, comm)`
"""
function Allreduce!(buf::MPIBuffertype{T}, op, comm::Comm) where T
    Allreduce!(MPI.IN_PLACE, buf, length(buf), op, comm)
end

"""
    Allreduce(sendbuf, op, comm)

Performs `op` reduction on the buffer `sendbuf`, allocating and returning the
output buffer in all processes of the group.

To specify the output buffer or perform the operation in pace, see [`Allreduce!`](@ref).
"""
function Allreduce(sendbuf::MPIBuffertype{T}, op, comm::Comm) where T

  recvbuf = similar(sendbuf)
  Allreduce!(sendbuf, recvbuf, length(recvbuf), op, comm)
end

function Allreduce(sendbuf::Array{T, N}, op, comm::Comm) where {T, N}
    recvbuf = Array{T,N}(undef, size(sendbuf))
    Allreduce!(sendbuf, recvbuf, length(sendbuf), op, comm)
end

function Allreduce(obj::T, op, comm::Comm) where T
    objref = Ref(obj)
    outref = Ref{T}()
    Allreduce!(objref, outref, 1, op, comm)

    outref[]
end

# Deprecation warning for lowercase allreduce that was used until v. 0.7.2
# Should be removed at some point in the future
function allreduce(sendbuf::MPIBuffertype{T}, op,
                   comm::Comm) where T
    @warn "`allreduce` is deprecated, use `Allreduce` instead."
    Allreduce(sendbuf, op, comm)
end

"""
    Scatter!(sendbuf, recvbuf, count, root, comm)

Splits the buffer `sendbuf` in the `root` process into `Comm_size(comm)` chunks
and sends the j-th chunk to the process of rank j into the `recvbuf` buffer,
which must be of length at least `count`.

`count` should be the same for all processes. If the number of elements varies between
processes, use [`Scatter!`](@ref) instead.

To perform the reduction in place, see [`Scatter_in_place!`](@ref).

To handle allocation of the output buffer, see [`Scatter`](@ref).
"""
function Scatter!(sendbuf::MPIBuffertype{T}, recvbuf::MPIBuffertype{T},
                  count::Integer, root::Integer,
                  comm::Comm) where T
    typeof(recvbuf) <: AbstractArray && @assert length(recvbuf) >= count
    isroot = Comm_rank(comm) == root
    isroot && typeof(sendbuf) <: AbstractArray && @assert length(sendbuf) >= count*Comm_size(comm)

    # int MPI_Scatter(const void* sendbuf, int sendcount, MPI_Datatype sendtype,
    #                 void* recvbuf, int recvcount, MPI_Datatype recvtype, int root,
    #                 MPI_Comm comm)
    @mpichk ccall((:MPI_Scatter, libmpi), Cint,
                  (Ptr{T}, Cint, MPI_Datatype, Ptr{T}, Cint, MPI_Datatype, Cint, MPI_Comm),
                  sendbuf, count, mpitype(T), recvbuf, count, mpitype(T), root, comm)
    recvbuf
end

"""
    Scatter_in_place!(buf, count, root, comm)

Splits the buffer `buf` in the `root` process into `Comm_size(comm)` chunks
and sends the j-th chunk to the process of rank j. No data is sent to the `root`
process.

This is functionally equivalent to calling
```
if root == MPI.Comm_rank(comm)
    Scatter!(buf, MPI.IN_PLACE, count, root, comm)
else
    Scatter!(C_NULL, buf, count, root, comm)
end
```

To specify a separate output buffer, see [`Scatter!`](@ref).

To handle allocation of the output buffer, see [`Scatter`](@ref).
"""
function Scatter_in_place!(buf::MPIBuffertype{T},
                  count::Integer, root::Integer,
                  comm::Comm) where T
    if Comm_rank(comm) == root
        @mpichk ccall((:MPI_Scatter, libmpi), Cint,
                      (Ptr{T}, Cint, MPI_Datatype, Ptr{T}, Cint, MPI_Datatype, Cint, MPI_Comm),
                      buf, count, mpitype(T), MPI_IN_PLACE, count, mpitype(T), root, comm)
    else
        @mpichk ccall((:MPI_Scatter, libmpi), Cint,
                      (Ptr{T}, Cint, MPI_Datatype, Ptr{T}, Cint, MPI_Datatype, Cint, MPI_Comm),
                      C_NULL, count, mpitype(T), buf, count, mpitype(T), root, comm)
    end
    buf
end

"""
    Scatter(sendbuf, count, root, comm)

Splits the buffer `sendbuf` in the `root` process into `Comm_size(comm)` chunks
and sends the j-th chunk to the process of rank j, allocating the output buffer.
"""
function Scatter(sendbuf::MPIBuffertype{T}, count::Integer, root::Integer,
                 comm::Comm) where T
    recvbuf = Array{T}(undef, count)
    Scatter!(sendbuf, recvbuf, count, root, comm)
end

"""
    Scatterv!(sendbuf, recvbuf, counts, root, comm)

Splits the buffer `sendbuf` in the `root` process into `Comm_size(comm)` chunks
of length `counts[j]` and sends the j-th chunk to the process of rank j into the
`recvbuf` buffer, which must be of length at least `count`.

To perform the reduction in place refer to [`Scatterv_in_place!`](@ref).
"""
function Scatterv!(sendbuf::MPIBuffertype{T}, recvbuf::MPIBuffertype{T},
                  counts::Vector{Cint}, root::Integer, comm::Comm) where T
    recvcnt = counts[Comm_rank(comm) + 1]
    disps = accumulate(+, counts) - counts
    typeof(recvbuf) <: AbstractArray && @assert length(recvbuf) >= recvcnt
    # int MPI_Scatterv(const void* sendbuf, const int sendcounts[],
    #                  const int displs[], MPI_Datatype sendtype, void* recvbuf,
    #                  int recvcount, MPI_Datatype recvtype, int root, MPI_Comm comm)
    @mpichk ccall((:MPI_Scatterv, libmpi), Cint,
                  (Ptr{T}, Ptr{Cint}, Ptr{Cint}, MPI_Datatype, Ptr{T}, Cint, MPI_Datatype, Cint, MPI_Comm),
                  sendbuf, counts, disps, mpitype(T), recvbuf, recvcnt, mpitype(T), root, comm)
    recvbuf
end

"""
    Scatterv(sendbuf, counts, root, comm)

Splits the buffer `sendbuf` in the `root` process into `Comm_size(comm)` chunks
of length `counts[j]` and sends the j-th chunk to the process of rank j, which
allocates the output buffer
"""
function Scatterv(sendbuf::MPIBuffertype{T}, counts::Vector{Cint},
                   root::Integer, comm::Comm) where T
    recvbuf = Array{T}(undef, counts[Comm_rank(comm) + 1])
    Scatterv!(sendbuf, recvbuf, counts, root, comm)
end

"""
    Scatterv_in_place(buf, counts, root, comm)

Splits the buffer `buf` in the `root` process into `Comm_size(comm)` chunks
of length `counts[j]` and sends the j-th chunk to the process of rank j into the
`buf` buffer, which must be of length at least `count`. The `root` process sends
nothing to itself.

This is functionally equivalent to calling
```
if root == MPI.Comm_rank(comm)
    Scatterv!(buf, MPI.IN_PLACE, counts, root, comm)
else
    Scatterv!(C_NULL, buf, counts, root, comm)
end
```
"""
function Scatterv_in_place!(buf::MPIBuffertype{T}, counts::Vector{Cint},
                           root::Integer, comm::Comm) where T
    recvcnt = counts[Comm_rank(comm) + 1]
    disps = accumulate(+, counts) - counts
    typeof(buf) <: AbstractArray && @assert length(buf) >= recvcnt

    if Comm_rank(comm) == root
        @mpichk ccall((:MPI_Scatterv, libmpi), Cint,
                      (Ptr{T}, Ptr{Cint}, Ptr{Cint}, MPI_Datatype, Ptr{T}, Cint, MPI_Datatype, Cint, MPI_Comm),
                      buf, counts, disps, mpitype(T), MPI_IN_PLACE, recvcnt, mpitype(T), root, comm)
    else
        @mpichk ccall((:MPI_Scatterv, libmpi), Cint,
                      (Ptr{T}, Ptr{Cint}, Ptr{Cint}, MPI_Datatype, Ptr{T}, Cint, MPI_Datatype, Cint, MPI_Comm),
                      C_NULL, counts, disps, mpitype(T), buf, recvcnt, mpitype(T), root, comm)
    end
    buf
end

"""
    Gather!(sendbuf, recvbuf, count, root, comm)

Each process sends the first `count` elements of the buffer `sendbuf` to the
`root` process. The `root` process stores elements in rank order in the buffer
buffer `recvbuf`.

`count` should be the same for all processes. If the number of elements varies between
processes, use [`Gatherv!`](@ref) instead.

To perform the reduction in place refer to [`Gather_in_place!`](@ref).
"""
function Gather!(sendbuf::MPIBuffertype{T}, recvbuf::MPIBuffertype{T},
                count::Integer, root::Integer, comm::Comm) where T
    typeof(sendbuf) <: AbstractArray && @assert length(sendbuf) >= count
    isroot = Comm_rank(comm) == root
    isroot && typeof(recvbuf) <: AbstractArray && @assert length(recvbuf) >= count*Comm_size(comm)

    # int MPI_Gather(const void* sendbuf, int sendcount, MPI_Datatype sendtype,
    #                void* recvbuf, int recvcount, MPI_Datatype recvtype, int root,
    #                MPI_Comm comm)
    @mpichk ccall((:MPI_Gather, libmpi), Cint,
                  (Ptr{T}, Cint, MPI_Datatype, Ptr{T}, Cint, MPI_Datatype, Cint, MPI_Comm),
                  sendbuf, count, mpitype(T), recvbuf, count, mpitype(T), root, comm)
    isroot ? recvbuf : nothing
end

"""
    Gather(sendbuf[, count=length(sendbuf)], root, comm)

Each process sends the first `count` elements of the buffer `sendbuf` to the
`root` process. The `root` allocates the output buffer and stores elements in
rank order.
"""
function Gather(sendbuf::MPIBuffertype{T}, count::Integer,
                root::Integer, comm::Comm) where T
    isroot = Comm_rank(comm) == root
    recvbuf = Array{T}(undef, isroot ? Comm_size(comm) * count : 0)
    Gather!(sendbuf, recvbuf, count, root, comm)
end

function Gather(sendbuf::Array{T}, root::Integer, comm::Comm) where T
    Gather(sendbuf, length(sendbuf), root, comm)
end

function Gather(sendbuf::SubArray{T}, root::Integer, comm::Comm) where T
    @assert Base.iscontiguous(sendbuf)
    Gather(sendbuf, length(sendbuf), root, comm)
end

function Gather(object::T, root::Integer, comm::Comm) where T
    isroot = Comm_rank(comm) == root
    sendbuf = T[object]
    recvbuf = Gather(sendbuf, root, comm)
    isroot ? recvbuf : nothing
end

"""
    Gather_in_place!(buf, count, root, comm)

Each process sends the first `count` elements of the buffer `buf` to the
`root` process. The `root` process stores elements in rank order in the buffer
buffer `buf`, sending no data to itself.

This is functionally equivalent to calling
```
if root == MPI.Comm_rank(comm)
    Gather!(MPI.IN_PLACE, buf, count, root, comm)
else
    Gather!(buf, C_NULL, count, root, comm)
end
```
"""
function Gather_in_place!(buf::MPIBuffertype{T}, count::Integer, root::Integer,
                          comm::Comm) where T
    if Comm_rank(comm) == root
        typeof(buf) <: AbstractArray && @assert length(buf) >= count*Comm_size(comm)
        @mpichk ccall((:MPI_Gather, libmpi), Cint,
                      (Ptr{T}, Cint, MPI_Datatype, Ptr{T}, Cint, MPI_Datatype, Cint, MPI_Comm),
                      MPI_IN_PLACE, count, mpitype(T), buf, count, mpitype(T), root, comm)
    else
        @mpichk ccall((:MPI_Gather, libmpi), Cint,
                      (Ptr{T}, Cint, MPI_Datatype, Ptr{T}, Cint, MPI_Datatype, Cint, MPI_Comm),
                      buf, count, mpitype(T), C_NULL, count, mpitype(T), root, comm)
    end
    buf
end

"""
    Allgather!(sendbuf, recvbuf, count, comm)

Each process sends the first `count` elements of `sendbuf` to the
other processes, who store the results in rank order into
`recvbuf`.

If `sendbuf==MPI.IN_PLACE` the input data is assumed to be in the
area of `recvbuf` where the process would receive it's own
contribution.
"""
function Allgather!(sendbuf::MPIBuffertypeOrConst{T}, recvbuf::MPIBuffertype{T},
                    count::Integer, comm::Comm) where T
    typeof(recvbuf) <: AbstractArray && @assert length(recvbuf) >= Comm_size(comm)*count
    # int MPI_Allgather(const void* sendbuf, int sendcount,
    #                   MPI_Datatype sendtype, void* recvbuf, int recvcount,
    #                   MPI_Datatype recvtype, MPI_Comm comm)
    @mpichk ccall((:MPI_Allgather, libmpi), Cint,
                  (Ptr{T}, Cint, MPI_Datatype, Ptr{T}, Cint, MPI_Datatype, MPI_Comm),
                  sendbuf, count, mpitype(T), recvbuf, count, mpitype(T), comm)
    recvbuf
end

"""
    Allgather!(buf, count, comm)

Equivalent to `Allgather!(MPI.IN_PLACE, buf, count, comm)`.
"""
function Allgather!(buf::MPIBuffertype{T}, count::Integer,
                   comm::Comm) where T
    Allgather!(MPI.IN_PLACE, buf, count, comm)
end

"""
    Allgather(sendbuf, count, comm)

Each process sends the first `count` elements of `sendbuf` to the
other processes, who store the results in rank order allocating
the output buffer.
"""
function Allgather(sendbuf::MPIBuffertype{T}, count::Integer,
                   comm::Comm) where T
    recvbuf = Array{T}(undef, Comm_size(comm) * count)
    Allgather!(sendbuf, recvbuf, count, comm)
end

function Allgather(sendbuf::Array{T}, comm::Comm) where T
    Allgather(sendbuf, length(sendbuf), comm)
end

function Allgather(sendbuf::SubArray{T}, comm::Comm) where T
    @assert Base.iscontiguous(sendbuf)
    Allgather(sendbuf, length(sendbuf), comm)
end

function Allgather(object::T, comm::Comm) where T
    sendbuf = T[object]
    recvbuf = Allgather(sendbuf, comm)
    recvbuf
end

"""
    Gatherv!(sendbuf, recvbuf, counts, root, comm)

Each process sends the first `counts[rank]` elements of the buffer `sendbuf` to
the `root` process. The `root` stores elements in rank order in the buffer
`recvbuf`.

To perform the reduction in place refer to [`Gatherv_in_place!`](@ref).
"""
function Gatherv!(sendbuf::MPIBuffertype{T}, recvbuf::MPIBuffertype{T},
                  counts::Vector{Cint}, root::Integer, comm::Comm) where T
    isroot = Comm_rank(comm) == root
    displs = accumulate(+, counts) - counts
    sendcnt = counts[Comm_rank(comm) + 1]
    isroot && typeof(recvbuf) <: AbstractArray && @assert length(recvbuf) >= sum(counts)
    # int MPI_Gatherv(const void* sendbuf, int sendcount, MPI_Datatype sendtype,
    #                 void* recvbuf, const int recvcounts[], const int displs[],
    #                 MPI_Datatype recvtype, int root, MPI_Comm comm)
    @mpichk ccall((:MPI_Gatherv, libmpi), Cint,
                  (Ptr{T}, Cint, MPI_Datatype, Ptr{T}, Ptr{Cint}, Ptr{Cint}, MPI_Datatype, Cint, MPI_Comm),
                  sendbuf, sendcnt, mpitype(T), recvbuf, counts, displs, mpitype(T), root, comm)
    isroot ? recvbuf : nothing
end

"""
    Gatherv(sendbuf, counts, root, comm)

Each process sends the first `counts[rank]` elements of the buffer `sendbuf` to
the `root` process. The `root` allocates the output buffer and stores elements
in rank order.
"""
function Gatherv(sendbuf::MPIBuffertype{T}, counts::Vector{Cint},
                 root::Integer, comm::Comm) where T
    isroot = Comm_rank(comm) == root
    recvbuf = Array{T}(undef, isroot ? sum(counts) : 0)
    Gatherv!(sendbuf, recvbuf, counts, root, comm)
end

"""
    Gatherv_in_place!(buf, counts, root, comm)

Each process sends the first `counts[rank]` elements of the buffer `buf` to
the `root` process. The `root` allocates the output buffer and stores elements
in rank order.

This is functionally equivalent to calling
```
if root == MPI.Comm_rank(comm)
    Gatherv!(MPI.IN_PLACE, buf, counts, root, comm)
else
    Gatherv!(buf, C_NULL, counts, root, comm)
end
```
"""
function Gatherv_in_place!(buf::MPIBuffertype{T}, counts::Vector{Cint},
                           root::Integer, comm::Comm) where T
    isroot = Comm_rank(comm) == root
    displs = accumulate(+, counts) - counts
    sendcnt = counts[Comm_rank(comm) + 1]

    if isroot
        typeof(buf) <: AbstractArray && @assert length(buf) >= sum(counts)
        @mpichk ccall((:MPI_Gatherv, libmpi), Cint,
                      (Ptr{T}, Cint, MPI_Datatype, Ptr{T}, Ptr{Cint}, Ptr{Cint}, MPI_Datatype, Cint, MPI_Comm),
                      MPI_IN_PLACE, sendcnt, mpitype(T), buf, counts, displs, mpitype(T), root, comm)
    else
        @mpichk ccall((:MPI_Gatherv, libmpi), Cint,
                      (Ptr{T}, Cint, MPI_Datatype, Ptr{T}, Ptr{Cint}, Ptr{Cint}, MPI_Datatype, Cint, MPI_Comm),
                      buf, sendcnt, mpitype(T), C_NULL, counts, displs, mpitype(T), root, comm)
    end
    buf
end

"""
    Allgatherv!(sendbuf, recvbuf, counts, comm)

Each process sends the first `counts[rank]` elements of the buffer `sendbuf` to
all other process. Each process stores the received data in rank order in the
buffer `recvbuf`.

if `sendbuf==MPI.IN_PLACE` every process takes the data to be sent is taken from
the interval of `recvbuf` where it would store it's own data.
"""
function Allgatherv!(sendbuf::MPIBuffertypeOrConst{T}, recvbuf::MPIBuffertype{T},
	                     counts::Vector{Cint}, comm::Comm) where T
    typeof(recvbuf) <: AbstractArray && @assert length(recvbuf) >= sum(counts)
    displs = accumulate(+, counts) - counts
    sendcnt = counts[Comm_rank(comm) + 1]
    # int MPI_Allgatherv(const void* sendbuf, int sendcount,
    #                    MPI_Datatype sendtype, void* recvbuf, const int recvcounts[],
    #                    const int displs[], MPI_Datatype recvtype, MPI_Comm comm)
    @mpichk ccall((:MPI_Allgatherv, libmpi), Cint,
              (Ptr{T}, Cint, MPI_Datatype, Ptr{T}, Ptr{Cint}, Ptr{Cint}, MPI_Datatype, MPI_Comm),
              sendbuf, sendcnt, mpitype(T), recvbuf, counts, displs, mpitype(T), comm)
    recvbuf
end

"""
    Allgatherv(sendbuf, counts, comm)

Each process sends the first `counts[rank]` elements of the buffer `sendbuf` to
all other process. Each process allocates an output buffer and stores the
received data in rank order.
"""
function Allgatherv(sendbuf::MPIBuffertype{T}, counts::Vector{Cint},
                    comm::Comm) where T
    recvbuf = Array{T}(undef, sum(counts))
    Allgatherv!(sendbuf, recvbuf, counts, comm)
end

"""
    Alltoall!(sendbuf, recvbuf, count, comm)

Every process divides the buffer `sendbuf` into `Comm_size(comm)` chunks of
length `count`, sending the `j`-th chunk to the `j`-th process.
Every process stores the data received from the `j`-th process in the `j`-th
chunk of the buffer `recvbuf`.

```
rank    send buf                        recv buf
----    --------                        --------
 0      a,b,c,d,e,f       Alltoall      a,b,A,B,α,β
 1      A,B,C,D,E,F  ---------------->  c,d,C,D,γ,ψ
 2      α,β,γ,ψ,η,ν                     e,f,E,F,η,ν
```

If `sendbuf==MPI.IN_PLACE`, data is sent from the `recvbuf` and then
overwritten.
"""
function Alltoall!(sendbuf::MPIBuffertypeOrConst{T}, recvbuf::MPIBuffertype{T},
                   count::Integer, comm::Comm) where T

    # int MPI_Alltoall(const void* sendbuf, int sendcount, MPI_Datatype sendtype,
    #                  void* recvbuf, int recvcount, MPI_Datatype recvtype,
    #                  MPI_Comm comm)
    @mpichk ccall((:MPI_Alltoall, libmpi), Cint,
                  (Ptr{T}, Cint, MPI_Datatype, Ptr{T}, Cint, MPI_Datatype, MPI_Comm),
                  sendbuf, count, mpitype(T), recvbuf, count, mpitype(T), comm)
    recvbuf
end

"""
    Alltoall(sendbuf, count, comm)

Every process divides the buffer `sendbuf` into `Comm_size(comm)` chunks of
length `count`, sending the `j`-th chunk to the `j`-th process.
Every process allocates the output buffer and stores the data received from the
`j`-th process in the `j`-th chunk.

```
rank    send buf                        recv buf
----    --------                        --------
 0      a,b,c,d,e,f       Alltoall      a,b,A,B,α,β
 1      A,B,C,D,E,F  ---------------->  c,d,C,D,γ,ψ
 2      α,β,γ,ψ,η,ν                     e,f,E,F,η,ν
```
"""
function Alltoall(sendbuf::MPIBuffertype{T}, count::Integer,
                  comm::Comm) where T
    recvbuf = Array{T}(undef, Comm_size(comm)*count)
    Alltoall!(sendbuf, recvbuf, count, comm)
end

"""
    Alltoallv!(sendbuf::T, recvbuf::T, scounts, rcounts, comm)

`MPI.IN_PLACE` is not supported for this operation.
"""
function Alltoallv!(sendbuf::MPIBuffertype{T}, recvbuf::MPIBuffertype{T},
                   scounts::Vector{Cint}, rcounts::Vector{Cint},
                   comm::Comm) where T
    typeof(recvbuf) <: AbstractArray && @assert length(recvbuf) == sum(rcounts)

    sdispls = accumulate(+, scounts) - scounts
    rdispls = accumulate(+, rcounts) - rcounts
    # int MPI_Alltoallv(const void* sendbuf, const int sendcounts[],
    #                   const int sdispls[], MPI_Datatype sendtype, void* recvbuf,
    #                   const int recvcounts[], const int rdispls[],
    #                   MPI_Datatype recvtype, MPI_Comm comm)
    @mpichk ccall((:MPI_Alltoallv, libmpi), Cint,
                  (Ptr{T}, Ptr{Cint}, Ptr{Cint}, MPI_Datatype, Ptr{T}, Ptr{Cint}, Ptr{Cint}, MPI_Datatype, MPI_Comm),
                  sendbuf, scounts, sdispls, mpitype(T), recvbuf, rcounts, rdispls, mpitype(T), comm)
    recvbuf
end

function Alltoallv(sendbuf::MPIBuffertype{T}, scounts::Vector{Cint},
                   rcounts::Vector{Cint}, comm::Comm) where T
    recvbuf = Array{T}(undef, sum(rcounts))
    Alltoallv!(sendbuf, recvbuf, scounts, rcounts, comm)
end

function Scan(sendbuf::MPIBuffertype{T}, count::Integer,
              op::Union{Op,MPI_Op}, comm::Comm) where T
    recvbuf = Array{T}(undef, count)
    # int MPI_Scan(const void* sendbuf, void* recvbuf, int count,
    #              MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
    @mpichk ccall((:MPI_Scan, libmpi), Cint,
                  (Ptr{T}, Ptr{T}, Cint, MPI_Datatype, MPI_Op, MPI_Comm),
                  sendbuf, recvbuf, count, mpitype(T), op, comm)
    recvbuf
end

function Scan(object::T, op::Union{Op,MPI_Op}, comm::Comm) where T
    sendbuf = T[object]
    Scan(sendbuf,1,op,comm)
end

function Exscan(sendbuf::MPIBuffertype{T}, count::Integer,
                op::Union{Op,MPI_Op}, comm::Comm) where T
    recvbuf = Array{T}(undef, count)
    # int MPI_Exscan(const void* sendbuf, void* recvbuf, int count,
    #                MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
    @mpichk ccall((:MPI_Exscan, libmpi), Cint,
          (Ptr{T}, Ptr{T}, Cint, MPI_Datatype, MPI_Op, MPI_Comm),
          sendbuf, recvbuf, count, mpitype(T), op, comm)
    recvbuf
end

function Exscan(object::T, op::Union{Op,MPI_Op}, comm::Comm) where T
    sendbuf = T[object]
    Exscan(sendbuf,1,op,comm)
end


text:
module Discretization

using LinearAlgebra
using SparseArrays
using DataStructures

import sbp

struct DiscretizationType

        sizes::OrderedDict
        boundary_points::Array
        boundary_conditions::Array
        normal::Array
        boundaries::Array
        mat_boundaries::Array

        function DiscretizationType(names, sizes, boundary_points,
                                           boundary_conditions, normal)

        num_bnd = length(boundary_points)
        @assert num_bnd == length(boundary_conditions)

        vars = variables(names, sizes)

        for indices in boundary_points
        end

        end
end

struct SettingsType
        nx
        dt
        order
        c
        diss
        nu
        nt

        function SettingsType(nx, nu, order, c, T, diss)
                dt = nu / (nx - 1)
                nt = Int64(floor(T/dt))
                return new(nx, dt, order, c, diss, nu, nt)
        end
end

struct OperatorsType

        Dx::AbstractArray
        x::Array
        h::Float64
        Hx::AbstractArray
        Hxi::AbstractArray
        nx::Int64

        function OperatorsType(nx, order)
                Dx, x, h, Hx, Hxi = sbp_operators(nx, order)
                new(Dx, x, h, Hx, Hxi, nx)
        end

        function OperatorsType(nx, order, grid)
                Dx, x, h, Hx, Hxi, mx = staggered_sbp_operators(nx, order, grid)
                return new(Dx, x, h, Hx, Hxi, mx)
        end

end

function normal(nx)
        nl = spzeros(nx)
        nl[1] = -1
        nl[nx] = 1
        return nl
end

function variables(names, sizes)
        d = OrderedDict()

        for (name, size) in zip(names, sizes)
                d[name] = size
        end

        return d
end

function group(vars)
        d = OrderedDict()

        for var in vars
                for (k, v) in var
                        d[k] = v
                end
        end

        return d
end

function boundary(vars::OrderedDict, boundary_points)
        bnd = OrderedDict()
        for (k, v) in vars
                bnd[k] = boundary_points
        end
        return block_diag(restrict_all(vars, bnd))
end

function boundary_each(vars::OrderedDict, boundary_points)
        bnd = OrderedDict()
        i = 1
        for (k, v) in vars
                bnd[k] = boundary_points[i]
                i += 1
        end
        return block_diag(restrict_all(vars, bnd))
end

function restrict_all(vars::OrderedDict, indices::OrderedDict)
        d = OrderedDict()

        for (k, v) in vars
                d[k] = restrict(v, indices[k])
        end

        return d
end

function restrict(vsize, indices)
        res = spzeros(length(indices), vsize)
        for (j, idx) in enumerate(indices)
                res[j, idx] = 1
        end
        return res
end

function block_diag(vars::OrderedDict)

        rows = zeros(Int64, length(vars))
        cols = zeros(Int64, length(vars))

        i = 1
        for (k, v) in vars
                rows[i] = size(v, 1)
                cols[i] = size(v, 2)
                i += 1
        end

        B = sbp.Sparse.block_matrix(rows, cols)

        i = 1
        for (k, v) in vars
                B = sbp.Sparse.block_matrix_insert(B,rows,cols,i,i,v)
                i += 1
        end

        return B
end

function projection_matrix(A, R1, dR1)
        bnd = sbp.BoundaryOperator
        Xp, Lp = bnd.eigen_positive(A)
        Xm, Lm = bnd.eigen_negative(A)
        return projection_matrix_eig(Xp, Xm, R1, dR1)
end

function projection_matrix_eig(Xp, Xm, R1, dR1)
        bnd = sbp.BoundaryOperator
        L = bnd.boundary_operator(Xp, Xm, R1)
        #dLp = bnd.virtual_complementary_boundary_operator(Xp, Xm, dR1)
        dLp = bnd.boundary_operator(-Xp, Xm, dR1)
        P = bnd.virtual_solution_matrix(dLp, L)
        return P
end

function check_energy_stability(H, M)
        Q = H * M + (H * M)'
        return real(eigen(Q).values)
end

function check_stability(M)
        return real(eigen(M).values)
end

function sbp_operators(nx, order)
        # SBP operators
        st=sbp.Strand1994
        desc_D, desc_H = st.load_description(order=order)
        x, h, Dx, Hx = st.build_operators(desc_D, desc_H, nx)
        Hxi = spdiagm(0 => 1 ./ diag(Hx))
        return Dx, x, h, Hx, Hxi
end

function staggered_sbp_operators(nx, order, grid)
        # staggered SBP operators
        st=sbp.OP2019
        desc = st.load_description()
        xp, xm, h, Dp, Dm, Hp, Hm, Pp, Pm = st.build_operators(desc, nx, true)
        if grid == "p"
                Hxi = spdiagm(0 => 1 ./ diag(Hp))
                return Dp, xp, h, Hp, Hxi, nx
        elseif grid == "m"
                Hxi = spdiagm(0 => 1 ./ diag(Hm))
                return Dm, xm, h, Hm, Hxi, (nx + 1)
        else
                error("Unknown grid type")
        end
end

function sat(Lift, A, B, ops, P)
        sat_term_data = kron(Lift, ops.Hxi) *  B' * (A * P)
        sat_term = sat_term_data * B
        return sat_term, sat_term_data
end

function discretization(Ax, ops::OperatorsType)
        M = kron(Ax, -ops.Dx)
        return M
end

function refine(nx)
        return 2 * (nx - 1) + 1
end

function time_step(sys, q, dt, nt, functor=nothing)
        dq = 0 * q
        t = 0.0

        for step in 1:nt
                q = sbp.LSRK4.lsrk4(sys, q, dq, t, dt)
                t = step * dt

                if functor != nothing
                        functor(t, step, q)
                end
        end

        return q, t
end

function time_step_staggered(f, g, u, v, dt, nt::Int64, functor=nothing)
        t = 0.0

        for step in 1:nt
                u, v = sbp.StaggeredRK4.staggered_rk4(f, g, u, v, t, dt)
                t = step * dt
                if functor != nothing
                        functor(t, step, u, v)
                end
        end

        return u, v, t
end

function time_step_staggered_linear(Ap, Av, p, v, dt, nt::Int64, functor=nothing)
        t = 0.0

        for step in 1:nt
                p, v = sbp.StaggeredRK4.staggered_rk4_linear(Ap, Av, p, v, dt)
                t = step * dt
                if functor != nothing
                        functor(t, step, p, v)
                end
        end

        return p, v, t
end

function functional(u, norm)
        return sum(norm * u)
end

function dissipation(A, P)
        diss = P' * A * P
        diss_eig = maximum(abs.(eigen(diss).values))
        return diss_eig
end

function spectral_radius(M, h)
        rho = h * maximum(abs.(eigen(M).values))
        return rho
end

function convergence_study(nx0, nt0, dt0, nrefine, errors, rates, init_fcn, disc_fcn, err_fcn, info_fcn=nothing)
        err = []
        nx = nx0
        nt = nt0
        dt = dt0
        grid_points = []

        for i in 1:nrefine
                f, ops = disc_fcn(nx)
                u = init_fcn(ops, nx)
                u, t = time_step(f, u, dt, nt)
                err_fcn(errors, u, ops, t)

                if info_fcn != nothing
                        info_fcn(i, nx, nt, dt, ops, errors)
                end
                append!(grid_points, nx)

                nx = refine(nx)
                nt = refine(nt)
                dt = 0.5 * dt
        end

        for (k, v) in errors
            rates[k] = sbp.MMS.log2_convergence_rates(v)
        end

        return grid_points
end

function convergence_study_staggered(nx0, nt0, dt0, nrefine, errors, rates, init_fcn, disc_fcn, err_fcn, info_fcn=nothing)
        err = []
        nx = nx0
        nt = nt0
        dt = dt0
        grid_points = []

        for i in 1:nrefine
                f, g, ops = disc_fcn(nx)
                u, v = init_fcn(ops, nx)
                u, v, t = time_step_staggered(f, g, u, dt, nt)
                err_fcn(errors, u, v, ops, t)

                if info_fcn != nothing
                        info_fcn(i, nx, nt, dt, ops, errors)
                end
                append!(grid_points, nx)

                nx = refine(nx)
                nt = refine(nt)
                dt = 0.5 * dt
        end

        for (k, v) in errors
            rates[k] = sbp.MMS.log2_convergence_rates(v)
        end

        return grid_points
end

end

max_stars_repo_name: ooreilly/sbp

max_issues_repo_name: ooreilly/sbp

max_forks_repo_name: ooreilly/sbp


text:
Random.seed!(2019)

# Dynamics
dt=0.1
function dynamics(x,u)
    v,ϕ = u
    xp = x
    xp += dt*([v*cos(x[3]);v*sin(x[3]); ϕ])
    xp[3] = xp[3]
    return xp
end

W = dt*0.05*Matrix{Float64}(I,3,3)

dmodel = NonlinearDynamicsModel(dynamics,W)

# Measurement
observe(x,u) = [norm(x[1:2])]

V = 0.1*Matrix{Float64}(I,1,1)

omodel = NonlinearObservationModel(observe,V)

# Build UKF
ukf = UnscentedKalmanFilter(dmodel,omodel);

# Simulation
times = 0:dt:5
action_sequence = [[1.0, 2*sin(pi/2*t)] for t in times]

b0 = GaussianBelief([0.0,0.0,0.0], 0.01*Matrix{Float64}(I,3,3))

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

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

# Tests
@test ukf.d isa NonlinearDynamicsModel
@test ukf.o isa NonlinearObservationModel
@test length(sim_measurements) == length(action_sequence)

points, w_μ, w_Σ = @inferred unscented_transform(b0, ukf.λ, ukf.α, ukf.β)
bp = @inferred unscented_transform_inverse(points, w_μ, w_Σ)

first_update = update(ukf, b0, action_sequence[1], sim_measurements[1])
first_predict = predict(ukf, b0, action_sequence[1])
first_predict_measure = measure(ukf, first_predict, 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], [0.16459, 3.34923, 2.65721];atol=0.001)
array_isapprox(sim_measurements[end], [3.17086];atol=0.001)
array_isapprox(filtered_beliefs[end].μ, [1.83418, 2.78963, 1.37579];
                atol=0.001)
array_isapprox(filtered_beliefs[end].Σ, [0.739291 -0.472135 -0.282821;
                                        -0.472135 0.342115 0.197056;
                                        -0.282821 0.197056 0.185713]; atol=0.001)