text:
import LinearAlgebra.BLAS: @blasfunc
import LinearAlgebra: BlasInt
import LinearAlgebra.LAPACK: stegr!
const liblapack = Base.liblapack_name

mutable struct StegrWork{T<:Real}
    jobz::Char
    range::Char
    dv::Vector{T}
    ev::Vector{T}
    vl::Real
    vu::Real
    il::BlasInt
    iu::BlasInt
    abstol::Vector{T}
    m::Vector{BlasInt}
    w::Vector{T}
    Z::Matrix{T}
    isuppz::Vector{BlasInt}
    work::Vector{T}
    lwork::BlasInt
    iwork::Vector{BlasInt}
    liwork::BlasInt
    info::Vector{BlasInt}
end

for (stegr,elty) in ((:dstegr_,:Float64),
                     (:sstegr_,:Float32))
    @eval begin
        function stegr!(n::BlasInt, sw::StegrWork{$elty})
            ldz = stride(sw.Z, 2)
            ccall((@blasfunc($stegr), liblapack), Cvoid,
                  (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty},
                   Ptr{$elty}, Ref{$elty}, Ref{$elty}, Ref{BlasInt},
                   Ref{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty},
                   Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
                   Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt}, Ptr{BlasInt}),
                  sw.jobz, sw.range, n,
                  sw.dv, sw.ev,
                  sw.vl, sw.vu, sw.il, sw.iu,
                  sw.abstol, sw.m,
                  sw.w, sw.Z, ldz,
                  sw.isuppz, sw.work, sw.lwork, sw.iwork, sw.liwork,
                  sw.info)
        end
    end
end

"""
    stegr!(α, β, sw)

Diagonalize the real-symmetric tridiagonal matrix with `α` on the
diagonal and `β` on the super-/subdiagonal, using the workspaces
allocated in `sw`.
"""
function stegr!(α::AbstractVector{T}, β::AbstractVector{T}, sw::StegrWork{T}) where T
    # @assert length(sw.dv) >= length(α)
    # @assert length(sw.ev) >= length(β)
    copyto!(sw.dv, α)
    copyto!(sw.ev, β)
    stegr!(BlasInt(length(α)), sw)
end

"""
    StegrWork(T, n)

Allocate work arrays for diagonalization of real-symmetric tridiagonal
matrices of sizes up to `n`×`n`.
"""
function StegrWork(::Type{T}, n::BlasInt,
                   jobz::Char = 'V', range::Char = 'A') where T
    dv = Array{T}(undef, n)
    ev = Array{T}(undef, n)
    abstol = Array{T}(undef, 1)
    m = Vector{BlasInt}(undef, 1)
    w = Array{T}(undef, n)
    ldz = jobz == 'N' ? 1 : n
    Z = Array{T}(undef, ldz, n)
    isuppz = Array{BlasInt}(undef, 2n)
    work = Array{T}(undef, 1)
    lwork = -one(BlasInt)
    iwork = Array{BlasInt}(undef, 1)
    liwork = -one(BlasInt)
    info = Array{BlasInt}(undef, 1)
    sw = StegrWork(jobz, range,
                   dv, ev,
                   0.0, 0.0,
                   BlasInt(0), BlasInt(0),
                   abstol, m,
                   w, Z,
                   isuppz,
                   work, lwork,
                   iwork, liwork,
                   info)
    # Call stegr! once to query for necessary workspace sizes.
    stegr!(n, sw)
    sw.lwork = BlasInt(sw.work[1])
    sw.work = Array{T}(undef, sw.lwork)
    sw.liwork = sw.iwork[1]
    sw.iwork = Array{BlasInt}(undef, sw.liwork)
    sw
end

mutable struct StegrCache{T,R<:Real} <: HermitianSubspaceCache{T}
    v::Vector{T} # Subspace-propagated vector
    w::Vector{T}
    sw::StegrWork{R}
    StegrCache(::Type{T}, n::Integer) where T = new{T,real(T)}(
        Vector{T}(undef, n), Vector{T}(undef, n),
        StegrWork(real(T), BlasInt(n)))
end

StegrCache(Ks::KrylovSubspace{B,T,U}) where {B,T,U} =
    StegrCache(T, Ks.maxiter)

"""
    expT!(α, β, t, cache)

Calculate the subspace exponential `exp(t*T)` for a tridiagonal
subspace matrix `T` with `α` on the diagonal and `β` on the
super-/subdiagonal, diagonalizing via `stegr!`.
"""
function expT!(α::AbstractVector{R}, β::AbstractVector{R}, t::Number,
               cache::StegrCache{T,R}) where {T,R<:Real}
    stegr!(α, β, cache.sw)
    sel = 1:length(α)
    @inbounds for i = sel
        cache.w[i] = exp(t*cache.sw.w[i])*cache.sw.Z[1,i]
    end
    mul!(@view(cache.v[sel]), @view(cache.sw.Z[sel,sel]), @view(cache.w[sel]))
end


text:
@testset "3.longest-substring-without-repeating-characters.jl" begin
    @test length_of_longest_substring("abcabcbb") == 3
    @test length_of_longest_substring("b") == 1
    @test length_of_longest_substring("bbbb") == 1
    @test length_of_longest_substring("pwwkew") == 3
    @test length_of_longest_substring("") == 0
    @test length_of_longest_substring("abac") == 3
end

ext: jl


text:
module ITensorGaussianMPS

using Compat
using ITensors
using ITensors.NDTensors
using LinearAlgebra

import LinearAlgebra: Givens

export slater_determinant_to_mps,
  slater_determinant_to_gmps,
  hopping_hamiltonian,
  slater_determinant_matrix,
  slater_determinant_to_gmera

include("gmps.jl")
include("gmera.jl")

end


text:
using NLsolve, LinearAlgebra, RootedTrees
"""
`Base.length(tab::ODERKTableau)`

Defines the length of a Runge-Kutta method to be the number of stages.
"""
Base.length(tab::ODERKTableau) = tab.stages

"""
`stability_region(z,tab::ODERKTableau)`

Calculates the stability function from the tableau at `z`. Stable if <1.

```math
r(z) = 1 + z bᵀ(I - zA)⁻¹ e
```
where e denotes a vector of ones.
"""
function stability_region(z,tab::ODERKTableau; embedded=false)
  A, c = tab.A, tab.c
  b = embedded ? tab.αEEst : tab.α
  e = ones(eltype(A), length(b))
  stages = (I - z*A) \ e
  1 + z * (transpose(b) * stages)
end

"""
`stability_region(tab::ODERKTableau; initial_guess=-3.0)`

Calculates the length of the stability region in the real axis.
"""
function stability_region(tab::ODERKTableau; initial_guess=-3.0)
  residual! = function (resid, x)
    resid[1] = abs(stability_region(x[1], tab)) - 1
  end
  sol = nlsolve(residual!, [initial_guess])
  sol.zero[1]
end

function RootedTrees.residual_order_condition(tab::ODERKTableau, order::Int, reducer=nothing, mapper=x->x^2; embedded=false)
  A, c = tab.A, tab.c
  b = embedded ? tab.αEEst : tab.α
  if reducer === nothing
    resid = map(RootedTreeIterator(order)) do t
      residual_order_condition(t, A, b, c)
    end
  else
    resid = mapreduce(reducer, RootedTreeIterator(order)) do t
      mapper(residual_order_condition(t, A, b, c))
    end
  end
  return resid
end

isfsal(::ExplicitRKTableau{M,V,fsal}) where {M,V,fsal} = fsal
isfsal(::ImplicitRKTableau) = nothing
function check_tableau(tab; tol=10eps(1.0))
  order = all(i -> residual_order_condition(tab, i, +, abs) < tol, 1:tab.order)
  if !order
    error("Tableau's order is not correct.")
  end
  if tab.adaptiveorder != 0
    embedded_order = all(i -> residual_order_condition(tab, i, +, abs; embedded=true) < tol, tab.adaptiveorder)
    if !embedded_order
      error("Tableau's embedded order is not correct.")
    end
  end
  return true
end


text:
#   This file is part of VerTeX.jl. It is licensed under the MIT license
#   Copyright (C) 2018 Michael Reed

export cdpkg, @cdpkg

cdpkg(pkg) = cd(dirname(Base.find_package(string(pkg))))

macro cdpkg(pkg)
    cdpkg(pkg)
    return nothing
end

repos = Dict("julia"=>"~/.julia/vtx/")

manifest = Dict("julia"=>Dict())
dictionary = Dict()

function lookup(ref)
    s = split(ref,':';limit=2)
    if haskey(dictionary,s[1])
        key = length(s) > 1 ? s[2] : ""
        return haskey(dictionary[s[1]],key) ? dictionary[s[1]][key] : [nothing]
    else
        return [nothing]
    end
end

function getdepot()
    repodat = Dict()
    try
        open(joinpath(homedir(),".julia/vtx/Depot.toml"), "r") do f
            repodat = TOML.parse(read(f,String))
        end
    catch
    end
    for key in keys(repos)
        push!(repodat,key=>repos[key])
    end
    return repodat
end

function regdepot(depot,location)
    dep = getdepot()
    if haskey(dep,depot)
        dep[depot] = location
    else
        push!(dep,depot=>location)
    end
    try
        open(joinpath(homedir(),".julia/vtx/Depot.toml"), "w") do f
            write(f, dict2toml(dep))
        end
    catch
    end
end

regpkg(depot,location) = regdepot(string(depot),joinpath(dirname(location),"vtx"))

regcmd = :(println(joinpath(dirname(@__DIR__),"vtx")))

function updateref!(data)
    depot = haskey(data,"depot") ? data["depot"] : "julia"
    !haskey(manifest[depot],data["uuid"]) && push!(manifest[depot],data["uuid"]=>Dict())
    setval!(manifest[depot][data["uuid"]],"dir",data["dir"])
    # update manifest and dictionary
    for cat in ["label","cite"]
        if haskey(data,cat)
            setval!(manifest[depot][data["uuid"]],cat,data[cat])
            cat ≠ "cite" && for ref in data[cat]
                if cat == "label"
                    give = [data["uuid"],haskey(data,"depot") ? data["depot"] : "julia",data["dir"]]
                    s = split(ref,':';limit=2)
                    key = length(s) > 1 ? join(s[2]) : ""
                    haskey(dictionary,s[1]) ? push!(dictionary[s[1]],key=>give) : push!(dictionary,join(s[1])=>Dict(key=>give))
                end
            end
        end
    end
    for item in ["ref","used","deps","show"]
        if haskey(data,item)
            setval!(manifest[depot][data["uuid"]],item,data[item])
        else
            haskey(manifest[depot][data["uuid"]],item) && pop!(manifest[depot][data["uuid"]],item)
        end
    end
    return nothing
end

function updaterefby!(depot,key)
    for cat ∈ ["ref","used","deps","show"]
        # identify cross-references
        if haskey(manifest[depot][key],cat)
            for ref in manifest[depot][key][cat]
                if cat ≠ "show"
                    s = split(ref,':';limit=2)
                    s2 = length(s) > 2 ? s[2] : ""
                    if haskey(dictionary,s[1]) && haskey(dictionary[s[1]],s2)
                        addval!(manifest[dictionary[s[1]][s2][2]][dictionary[s[1]][s2][1]],cat*"by",[key,depot,manifest[depot][key]["dir"]])
                    end
                else
                    if haskey(manifest,ref[2]) && haskey(manifest[ref[2]],ref[1])
                        addval!(manifest[ref[2]][ref[1]],cat*"by",[key,depot,manifest[depot][key]["dir"]])
                    end
                end
            end
        end
        # remove surplus edges
        if haskey(manifest[depot][key],cat*"by")
            amt = length(manifest[depot][key][cat*"by"])
            k = 1
            while k ≤ amt
                ref = manifest[depot][key][cat*"by"][k]
                if haskey(manifest[ref[2]][ref[1]],cat) &&
                        ref ∉ manifest[ref[2]][ref[1]][cat]
                    deleteat!(manifest[depot][key][cat*"by"],k)
                    amt -= 1
                else
                    k += 1
                end
            end
            isempty(manifest[depot][key][cat*"by"]) && pop!(manifest[depot][key],cat*"by")
        end
    end
    return nothing
end

function scan(depot)
    depos = getdepot()
    !haskey(manifest,depot) && push!(manifest,depot=>Dict())
    for (root, dirs, files) in walkdir(checkhome(depos[depot]))
        for dir in dirs
            for file in readdir(joinpath(root,dir))
                data = nothing
                if endswith(file, ".vtx")
                    data = TOML.parsefile(joinpath(root,dir,file))
                    updateref!(data)
                end
            end
        end
    end
end

function scan()
    for depot ∈ keys(manifest)
        for key ∈ keys(manifest[depot])
            updaterefby!(depot,key)
        end
    end
end

resolve(depot) = haskey(getdepot(),depot) && (scan(depot); scan())

function resolve()
    for depot in keys(getdepot())
        scan(depot)
    end
    scan()
end

function save(dat::Dict,path::String;warn=true)
    out = deepcopy(dat)
    !haskey(out,"depot") && push!(out,"depot"=>"julia")
    repo = out["depot"]
    depos = getdepot()
    !haskey(depos,repo) && (@warn "did not save, $repo depot not found"; return dat)
    data = nothing
    try
        data = load(path,repo)
        if (data ≠ nothing)
            cmv = checkmerge(dat["revised"],data,dat["title"],dat["author"],dat["date"],dat["tex"],"Save/Overwrite?")
            if cmv == 0
                throw(error("VerTeX unable to proceed due to merge failure"))
            elseif cmv < 2
                @warn "skipped saving $path"
                return dat
            end
        end
    catch
    end
    way = joinpath(checkhome(depos[repo]),path)
    !isdir(dirname(way)) && mkpath(dirname(way))
    if haskey(dat,"dir") && (dat["dir"] ≠ path)
        #rm(joinpath(checkhome(depos[repo]),dat["dir"]))
        out["dir"] = path
    else
        push!(out,"dir"=>path)
    end
    infotxt = "saving VerTeX: $(out["title"])\n"
    old = data ≠ nothing ? data : dat
    # go through save queue from show list
    if haskey(out,"save")
        for it ∈ out["save"]
            save(it,out["ids"][it["uuid"]][3];warn=warn)
        end
        pop!(out,"save")
    end
    haskey(out,"compact") && pop!(out,"compact")
    updateref!(out)
    updaterefby!(repo,out["uuid"])
    if haskey(out,"edit")
        setval!(out,"revised",out["edit"])
        pop!(out,"edit")
    end
    if data ≠ nothing
        haskey(data,"revised") && pop!(data,"revised")
        compare = deepcopy(out)
        haskey(compare,"revised") && pop!(compare,"revised")
        data == compare && (return out)
    end
    open(way, "w") do f
        write(f, dict2toml(out))
    end
    warn && (@info infotxt*"saved at $path in $(out["depot"])")
    return out
end

function save(dat::Dict;warn=true)
    save(dat, haskey(dat,"dir") ? dat["dir"] : dat["uuid"];warn=warn)
end

function save(dat::Dict,path::String,repo::String;warn=true)
    out = deepcopy(dat)
    if haskey(dat["depot"])
        #rm(joinpath(checkhome(getdepot()[dat["depot"]]),dat["dir"]))
        out["depot"] = repo
        out["dir"] = path
    else
        push!(out,"depot"=>repo,"dir"=>path)
    end
    save(out,path;warn=warn)
end

function load(path::String,repo="julia")
    depos = getdepot()
    !haskey(depos,repo) && (@warn "did not load, $repo depot not found"; return path)
    dat = ""
    open(joinpath(checkhome(depos[repo]),path), "r") do f
        dat = read(f, String)
    end
    return TOML.parse(dat)
end

function loadpath(data::Dict,file::String="/tmp/doc.tex")
    load = ""
    g = getdepot()
    if haskey(data,"dir") && (data["depot"] ∈ keys(g))
        load = joinpath(checkhome(g[data["depot"]]),data["dir"])
        load = replace(load,r".vtx$"=>".tex")
        !occursin(r".tex$",load) && (load = load*".tex")
    else
        load = file
    end
    return load
end

function writetex(data::Dict,file::String="/tmp/doc.tex")
    load = loadpath(data,file)
    # check if file actually exists yet, if not create it.
    open(load, "w") do f
        # check if tex file actually needs to be updated?
        write(f, VerTeX.dict2tex(data))
    end
    return load
end

function readtex(load::String)
    out = ""
    open(load, "r") do f
        out = read(f,String)
    end
    return out
end

function update(data::Dict)
    save(tex2dict(readtex(loadpath(data)),data))
end

function writemanifest(depot)
    depos = getdepot()
    if haskey(manifest,depot) && haskey(depos,depot)
        open(joinpath(checkhome(depos[depot]),"Manifest.toml"), "w") do f
            write(f, dict2toml(manifest[depot]))
        end
    else
        @warn "no $depot manifest found in memory"
    end
    return nothing
end

function writemanifest()
    for key ∈ keys(manifest)
        writemanifest(key)
    end
    return nothing
end

function writedictionary()
    depos = getdepot()
    open(joinpath(checkhome(depos["julia"]),"Dictionary.toml"), "w") do f
        write(f, dict2toml(dictionary))
    end
    return nothing
end

function readmanifest(depot)
    depos = getdepot()
    if haskey(depos,depot)
        dat = ""
        try
            open(joinpath(checkhome(depos[depot]),"Manifest.toml"), "r") do f
                dat = read(f, String)
            end
            setval!(manifest,depot,TOML.parse(dat))
        catch
        end
    else
        @warn "did not load, $depot depot not found"
    end
    return nothing
end

function readmanifest()
    depos = getdepot()
    for depot ∈ keys(depos)
        readmanifest(depot)
    end
    return nothing
end

function readdictionary()
    depos = getdepot()
    dat = ""
    try
        open(joinpath(checkhome(depos["julia"]),"Dictionary.toml"), "r") do f
            dat = read(f, String)
        end
        global dictionary
        dictionary = TOML.parse(dat)
    catch
    end
    return nothing
end


text:
"""
    BackPropagation(NeuralNetwork)

Backpropagation process for the neural network.

    Usage:
```julia-repl
BackPropagation(NeuralNetwork)
```
Arguments:
- `NeuralNetwork` : The NeuralNetwork structure.
"""
function BackPropogation(NeuralNetwork, ActivationFunction::String = "Sigmoid")

    # Check activation function argument

    (ActivationFunction == "Sigmoid" || ActivationFunction == "Tanh" || ActivationFunction == "ReLU" || ActivationFunction == "Linear") || error("`ActivationFunction` should be a String of either 'Sigmoid', 'Tanh', 'ReLU' or 'Linear'.")

    # Use the chain rule to find derivative of loss function w.r.t. w₁ and w₂ (weights)

    if ActivationFunction == "Sigmoid"
        DerivOutput = SigmoidDerivative(NeuralNetwork.output)
        DerivLayer1 = SigmoidDerivative(NeuralNetwork.output)
    elseif ActivationFunction == "Tanh"
        DerivOutput = TanhDerivative(NeuralNetwork.output)
        DerivLayer1 = TanhDerivative(NeuralNetwork.output)
    elseif ActivationFunction == "ReLU"
        DerivOutput = ReLUDerivative(NeuralNetwork.output)
        DerivLayer1 = ReLUDerivative(NeuralNetwork.output)
    elseif ActivationFunction == "Linear"
        DerivOutput = LinearDerivative(NeuralNetwork.output)
        DerivLayer1 = LinearDerivative(NeuralNetwork.output)
    end

    ∂w₂ = transpose(NeuralNetwork.layer1) * ((2 * (NeuralNetwork.y - NeuralNetwork.output)) * DerivOutput)

    ∂w₁ = ((2 * (NeuralNetwork.y - NeuralNetwork.output)) * DerivOutput) * transpose(NeuralNetwork.w₂)

    ∂w₁ = ∂w₁ * DerivLayer1

    ∂w₁ = transpose(NeuralNetwork.input) * ∂w₁

    # Update weights using gradient of the loss function

    NeuralNetwork$w₁ <- NeuralNetwork$w₁ + ∂w₁
    NeuralNetwork$w₂ <- NeuralNetwork$w₂ + ∂w₂

    return NeuralNetwork
end


text:
using StochasticDiffEq
using DiffEqSensitivity
using DiffEqNoiseProcess
using LinearAlgebra, Statistics, Random
using Zygote, ReverseDiff, ForwardDiff
using Test
#using DifferentialEquations
seed = 12345
Random.seed!(seed)

function g(u,p,t)
  sum(u.^2.0/2.0)
end

function dg!(out,u,p,t,i)
  (out.=-u)
end

@testset "noise iip tests" begin
  function f(du,u,p,t,W)
    du[1] = p[1]*u[1]*sin(W[1] - W[2])
    du[2] = p[2]*u[2]*cos(W[1] + W[2])
    return nothing
  end
  dt = 1e-4
  u0 = [1.00;1.00]
  tspan = (0.0,5.0)
  t = tspan[1]:0.1:tspan[2]
  p = [2.0,-2.0]
  prob = RODEProblem(f,u0,tspan,p)

  sol = solve(prob,RandomEM(),dt=dt, save_noise=true)
  # check reversion with usage of Noise Grid
  _sol = deepcopy(sol)
  noise_reverse = NoiseGrid(reverse(_sol.t),reverse(_sol.W.W))
  prob_reverse = RODEProblem(f,_sol[end],reverse(tspan),p,noise=noise_reverse)
  sol_reverse = solve(prob_reverse,RandomEM(),dt=dt)
  @test sol.u ≈ reverse(sol_reverse.u) rtol=1e-3
  @show minimum(sol.u)

  # Test if Forward and ReverseMode AD agree.
  Random.seed!(seed)
  du0ReverseDiff,dpReverseDiff = Zygote.gradient((u0,p)->sum(
    Array(solve(prob,RandomEM(),dt=dt,u0=u0,p=p,saveat=t,sensealg=ReverseDiffAdjoint())).^2/2)
    ,u0,p)
  Random.seed!(seed)
  dForward = ForwardDiff.gradient((θ)->sum(
    Array(solve(prob,RandomEM(),dt=dt,u0=θ[1:2],p=θ[3:4],saveat=t)).^2/2)
    ,[u0;p])

  @info dForward

  @test du0ReverseDiff ≈ dForward[1:2]
  @test dpReverseDiff ≈ dForward[3:4]

  # test gradients
  Random.seed!(seed)
  sol = solve(prob,RandomEM(),dt=dt, save_noise=true, saveat=t)

  ###
  ## BacksolveAdjoint
  ###

  # ReverseDiff
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=BacksolveAdjoint())

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  @info du0, dp'

  # ReverseDiff with compiled tape
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=BacksolveAdjoint(autojacvec=DiffEqSensitivity.ReverseDiffVJP(true)))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=BacksolveAdjoint(autojacvec=DiffEqSensitivity.ZygoteVJP()))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # Tracker
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=BacksolveAdjoint(autojacvec=DiffEqSensitivity.TrackerVJP()))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # isautojacvec = true
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=BacksolveAdjoint(autojacvec=true))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # isautojacvec = false
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=BacksolveAdjoint(autojacvec=false))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # isautojacvec = false and with jac and paramjac
  function jac(J,u,p,t,W)
    J[1,1] = p[1]*sin(W[1] - W[2])
    J[2,1] = zero(u[1])
    J[1,2] = zero(u[1])
    J[2,2] = p[2]*cos(W[1] + W[2])
  end

  function paramjac(J,u,p,t,W)
    J[1,1] = u[1]*sin(W[1] - W[2])
    J[2,1] = zero(u[1])
    J[1,2] = zero(u[1])
    J[2,2] = u[2]*cos(W[1] + W[2])
  end
  Random.seed!(seed)
  faug = RODEFunction(f,jac=jac,paramjac=paramjac)
  prob_aug = RODEProblem{true}(faug,u0,tspan,p)
  sol = solve(prob_aug,RandomEM(),dt=dt, save_noise=true, saveat=t)
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=BacksolveAdjoint(autojacvec=false))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  ###
  ## InterpolatingAdjoint
  ###

  # test gradients with dense solution and no checkpointing
  Random.seed!(seed)
  sol = solve(prob,RandomEM(),dt=dt, save_noise=true, dense=true)

  # ReverseDiff
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(autojacvec=DiffEqSensitivity.ReverseDiffVJP()))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # ReverseDiff with compiled tape
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(autojacvec=DiffEqSensitivity.ReverseDiffVJP(true)))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # Zygote
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(autojacvec=DiffEqSensitivity.ZygoteVJP()))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # Tracker
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(autojacvec=DiffEqSensitivity.TrackerVJP()))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # isautojacvec = true
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(autojacvec=true))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # isautojacvec = false
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(autojacvec=false))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # isautojacvec = false and with jac and paramjac
  Random.seed!(seed)
  faug = RODEFunction(f,jac=jac,paramjac=paramjac)
  prob_aug = RODEProblem{true}(faug,u0,tspan,p)
  sol = solve(prob_aug,RandomEM(),dt=dt, save_noise=true, dense=true)
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(autojacvec=false))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # test gradients with saveat solution and checkpointing
  # need to simulate for dt beyond last tspan to avoid errors in NoiseGrid
  Random.seed!(seed)
  sol = solve(prob,RandomEM(),dt=dt, save_noise=true, dense=true)
  Random.seed!(seed)
  sol_long = solve(remake(prob, tspan=(tspan[1],tspan[2]+10dt)),RandomEM(),dt=dt, save_noise=true, dense=true)

  @test sol_long(t) ≈ sol(t) rtol=1e-12
  @test sol_long.W.W[1:end-10] ≈ sol.W.W[1:end] rtol=1e-12

  # test gradients with saveat solution and checkpointing
  noise = NoiseGrid(sol_long.W.t,sol_long.W.W)
  sol2 = solve(remake(prob,noise=noise,tspan=(tspan[1],tspan[2])),RandomEM(),dt=dt, saveat=t)

  @test sol_long(t) ≈ sol2(t) rtol=1e-12
  @test sol_long.W.W ≈ sol2.W.W rtol=1e-12

  # ReverseDiff
  du0, dp = adjoint_sensitivities(sol2,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(checkpointing=true,autojacvec=DiffEqSensitivity.ReverseDiffVJP()))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # ReverseDiff with compiled tape
  du0, dp = adjoint_sensitivities(sol2,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(checkpointing=true,autojacvec=DiffEqSensitivity.ReverseDiffVJP(true)))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # Zygote
  du0, dp = adjoint_sensitivities(sol2,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(checkpointing=true,autojacvec=DiffEqSensitivity.ZygoteVJP()))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # Tracker
  du0, dp = adjoint_sensitivities(sol2,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(checkpointing=true,autojacvec=DiffEqSensitivity.TrackerVJP()))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # isautojacvec = true
  du0, dp = adjoint_sensitivities(sol2,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(checkpointing=true,autojacvec=true))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # isautojacvec = false
  du0, dp = adjoint_sensitivities(sol2,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(checkpointing=true,autojacvec=false))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # isautojacvec = false and with jac and paramjac
  Random.seed!(seed)
  faug = RODEFunction(f,jac=jac,paramjac=paramjac)
  prob_aug = RODEProblem{true}(faug,u0,tspan,p, noise=noise)
  sol = solve(prob_aug,RandomEM(),dt=dt, save_noise=false, dense=false, saveat=t)
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(checkpointing=true,autojacvec=false))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2
end

@testset "noise oop tests" begin
  function f(u,p,t,W)
    dx = p[1]*u[1]*sin(W[1] - W[2])
    dy = p[2]*u[2]*cos(W[1] + W[2])
    return [dx,dy]
  end
  dt = 1e-4
  u0 = [1.00;1.00]
  tspan = (0.0,5.0)
  t = tspan[1]:0.1:tspan[2]
  p = [2.0,-2.0]
  prob = RODEProblem{false}(f,u0,tspan,p)

  sol = solve(prob,RandomEM(),dt=dt, save_noise=true)
  # check reversion with usage of Noise Grid
  _sol = deepcopy(sol)
  noise_reverse = NoiseGrid(reverse(_sol.t),reverse(_sol.W.W))
  prob_reverse = RODEProblem(f,_sol[end],reverse(tspan),p,noise=noise_reverse)
  sol_reverse = solve(prob_reverse,RandomEM(),dt=dt)
  @test sol.u ≈ reverse(sol_reverse.u) rtol=1e-3
  @show minimum(sol.u)

  # Test if Forward and ReverseMode AD agree.
  Random.seed!(seed)
  du0ReverseDiff,dpReverseDiff = Zygote.gradient((u0,p)->sum(
    Array(solve(prob,RandomEM(),dt=dt,u0=u0,p=p,saveat=t,sensealg=ReverseDiffAdjoint())).^2/2)
    ,u0,p)
  Random.seed!(seed)
  dForward = ForwardDiff.gradient((θ)->sum(
    Array(solve(prob,RandomEM(),dt=dt,u0=θ[1:2],p=θ[3:4],saveat=t)).^2/2)
    ,[u0;p])

  @info dForward

  @test du0ReverseDiff ≈ dForward[1:2]
  @test dpReverseDiff ≈ dForward[3:4]

  # test gradients
  Random.seed!(seed)
  sol = solve(prob,RandomEM(),dt=dt, save_noise=true, saveat=t)

  ###
  ## BacksolveAdjoint
  ###

  # ReverseDiff
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=BacksolveAdjoint(autojacvec=DiffEqSensitivity.ReverseDiffVJP()))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  @info du0, dp'

  # ReverseDiff with compiled tape
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=BacksolveAdjoint(autojacvec=DiffEqSensitivity.ReverseDiffVJP(true)))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # Zygote
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=BacksolveAdjoint(autojacvec=DiffEqSensitivity.ZygoteVJP()))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # Tracker
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=BacksolveAdjoint(autojacvec=DiffEqSensitivity.TrackerVJP()))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # isautojacvec = true
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=BacksolveAdjoint(autojacvec=true))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # isautojacvec = false
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=BacksolveAdjoint(autojacvec=false))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # isautojacvec = false and with jac and paramjac
  function jac(J,u,p,t,W)
    J[1,1] = p[1]*sin(W[1] - W[2])
    J[2,1] = zero(u[1])
    J[1,2] = zero(u[1])
    J[2,2] = p[2]*cos(W[1] + W[2])
  end

  function paramjac(J,u,p,t,W)
    J[1,1] = u[1]*sin(W[1] - W[2])
    J[2,1] = zero(u[1])
    J[1,2] = zero(u[1])
    J[2,2] = u[2]*cos(W[1] + W[2])
  end
  Random.seed!(seed)
  faug = RODEFunction(f,jac=jac,paramjac=paramjac)
  prob_aug = RODEProblem{false}(faug,u0,tspan,p)
  sol = solve(prob_aug,RandomEM(),dt=dt, save_noise=true, saveat=t)
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=BacksolveAdjoint(autojacvec=false))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  ###
  ## InterpolatingAdjoint
  ###

  # test gradients with dense solution and no checkpointing
  Random.seed!(seed)
  sol = solve(prob,RandomEM(),dt=dt, save_noise=true, dense=true)

  # ReverseDiff
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(autojacvec=DiffEqSensitivity.ReverseDiffVJP()))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # ReverseDiff with compiled tape
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(autojacvec=DiffEqSensitivity.ReverseDiffVJP(true)))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # Zygote
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(autojacvec=DiffEqSensitivity.ZygoteVJP()))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # Tracker
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(autojacvec=DiffEqSensitivity.TrackerVJP()))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # isautojacvec = true
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(autojacvec=true))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # isautojacvec = false
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(autojacvec=false))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # isautojacvec = false and with jac and paramjac
  Random.seed!(seed)
  faug = RODEFunction(f,jac=jac,paramjac=paramjac)
  prob_aug = RODEProblem{false}(faug,u0,tspan,p)
  sol = solve(prob_aug,RandomEM(),dt=dt, save_noise=true, dense=true)
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(autojacvec=false))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # test gradients with saveat solution and checkpointing
  # need to simulate for dt beyond last tspan to avoid errors in NoiseGrid
  Random.seed!(seed)
  sol = solve(prob,RandomEM(),dt=dt, save_noise=true, dense=true)
  Random.seed!(seed)
  sol_long = solve(remake(prob, tspan=(tspan[1],tspan[2]+10dt)),RandomEM(),dt=dt, save_noise=true, dense=true)

  @test sol_long(t) ≈ sol(t) rtol=1e-12
  @test sol_long.W.W[1:end-10] ≈ sol.W.W[1:end] rtol=1e-12

  # test gradients with saveat solution and checkpointing
  noise = NoiseGrid(sol_long.W.t,sol_long.W.W)
  sol2 = solve(remake(prob,noise=noise,tspan=(tspan[1],tspan[2])),RandomEM(),dt=dt, saveat=t)

  @test sol_long(t) ≈ sol2(t) rtol=1e-12
  @test sol_long.W.W ≈ sol2.W.W rtol=1e-12

  # ReverseDiff
  du0, dp = adjoint_sensitivities(sol2,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(checkpointing=true,autojacvec=DiffEqSensitivity.ReverseDiffVJP()))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # ReverseDiff with compiled tape
  du0, dp = adjoint_sensitivities(sol2,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(checkpointing=true,autojacvec=DiffEqSensitivity.ReverseDiffVJP(true)))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # Zygote
  du0, dp = adjoint_sensitivities(sol2,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(checkpointing=true,autojacvec=DiffEqSensitivity.ZygoteVJP()))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # Tracker
  du0, dp = adjoint_sensitivities(sol2,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(checkpointing=true,autojacvec=DiffEqSensitivity.TrackerVJP()))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # isautojacvec = true
  du0, dp = adjoint_sensitivities(sol2,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(checkpointing=true,autojacvec=true))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # isautojacvec = false
  du0, dp = adjoint_sensitivities(sol2,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(checkpointing=true,autojacvec=false))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2

  # isautojacvec = false and with jac and paramjac
  Random.seed!(seed)
  faug = RODEFunction(f,jac=jac,paramjac=paramjac)
  prob_aug = RODEProblem{false}(faug,u0,tspan,p,noise=noise)
  sol = solve(prob_aug,RandomEM(),dt=dt, save_noise=false, saveat=t, dense=false)
  du0, dp = adjoint_sensitivities(sol,RandomEM(),dg!,Array(t)
    ,dt=dt,adaptive=false,sensealg=InterpolatingAdjoint(checkpointing=true,autojacvec=false))

  @test du0ReverseDiff ≈ du0 rtol=1e-2
  @test dpReverseDiff ≈ dp' rtol=1e-2
end


text:
# # 2D decaying turbulence
#
# A simulation of decaying two-dimensional turbulence.

using Pkg; Pkg.instantiate()

using FourierFlows, Printf, Random, JLD2

using Random: seed!
using FFTW: rfft, irfft
using LinearAlgebra: ldiv!
using FourierFlows: parsevalsum

import GeophysicalFlows.TwoDNavierStokes
import GeophysicalFlows.TwoDNavierStokes: energy, enstrophy
import GeophysicalFlows: peakedisotropicspectrum

# ## Choosing a device: CPU or GPU

dev = CPU()     # Device (CPU/GPU)

# ## Numerical, domain, and simulation parameters
#
# First, we pick some numerical and physical parameters for our model.

α = 0.3                       # aspect ratio parameter α = Ly / Lx
γ = 1.0                       # forcing amplitude

critical_ν(α, γ) = sqrt((γ * sqrt(1 - α^2)) / (sqrt(2) * (1 + α^2)))

ny, Ly  = 64, 2π             # grid resolution and domain length
nx, Lx  = 128, Ly / α

## Then we pick the time-stepper parameters
    dt = 5e-3  # timestep
 nsubs = 50    # number of steps between each plot

  ν = critical_ν(α, γ) / 10

 k₀ = 1
 tfinal = 10.0 / (ν * k₀^2)
 nsteps = Int(round(tfinal / dt))

 grid = TwoDGrid(dev, nx, Lx, ny, Ly)
 x, y = gridpoints(grid)

 forcing = @. γ * cos(y)
 forcingh = rfft(forcing)

 # Function that constructs the vorticity forcing in Fourier space
 function calcF!(Fh, sol, t, clock, vars, params, grid)
   @. Fh = forcingh
   return nothing
 end

# Initial condition (as function)
initial_vorticity(x, y) = sin(α * x) * sin(y) + 0.05 * (cos(α * 8x) + cos(8y))
initial_vorticity(x, y) = γ / ν * cos(y)

# Setting initial conditions
seed!(1234)

ζ_perturbation = 1 * randn(nx, ny)

ζ₀ = @. initial_vorticity(x, y) + ζ_perturbation

include("simulation.jl")

max_stars_repo_name: navidcy/KolmogorovFlow

max_issues_repo_name: navidcy/KolmogorovFlow

max_forks_repo_name: navidcy/KolmogorovFlow


text:
using Test
using BartIO

include("testIO.jl");


text:
# L!Type is the type of Llama types.
abstract L!Type

# L! is the type of Llama values.
abstract L!
type L!int <: L! val::Int64 end
type L!str <: L! val::String end
type L!Object <: L!
  _type::L!Type
  constr::String
  values::Array{L!, 1}
  attributes::Dict{String, L!}
  extends::Union(L!, Nothing)
end

L!Object(_type::L!Type, constr::String) =
  L!Object(_type, constr, [], Dict(), Nothing())

# A constraint on a type. Asserts `name` exists of type `_type`.
type L!Constraint
  name::String
  _type::L!Type
end

# A trait is a set of one or more constraints.
type L!Trait
  constraints::Set{L!Constraint}
end

# A type variable, for polymorphic types.
type L!TVar <: L!Type name::String end

# A constant type. Defines what constructors and attributes are legal,
# and might have a parent.
type L!TConst <: L!Type
  name::String
  constructors::Array{String, 1}
  parent::Union(L!TConst, Nothing)
  attributes::Dict{String, L!Type}
end

type L!TApply <: L!Type
  left::L!Type
  right::L!Type
end

L!TConst(name::String) = L!TConst(name, [], Nothing(), Dict())

L!TApply(t1::L!Type, t2::L!Type, t3::L!Type) = L!TApply(L!TApply(t1, t2), t3)
function L!TApply(types)
  if length(types) < 2
    error("Not enough types to apply, need at least 2")
  end
  t = L!TApply(types[1], types[2])
  for t_ in types[3:end]
    t = L!TApply(t, t_)
  end
  t
end

L!absfunc!t = L!TConst("->")
L!func!t(from::L!Type, to::L!Type) = L!TApply(L!absfunc!t, from, to)
L!int!t = L!TConst("Int")
L!str!t = L!TConst("Str")

function L!tuple(types::Array{L!Type, 1})
  tup = L!TConst("Tuple($(length(types)))")
  if length(types) == 0
    return tup
  elseif length(types) == 1
    return types[1]
  end
  L!TApply(vcat([tup], types))
end

f = L!func!t(L!int!t, L!str!t)

println(f)
println()
println(L!tuple([L!int!t, L!str!t, f]))

type L!ValueSet
  values::Dict{L!Type, L!}
end

function single(v::L!)
  L!ValueSet({v._type => v})
end

function addval(vset::L!ValueSet, v::L!)
  insert!(vset.values, v._type, v)
end

i = L!int(1)
Maybe = L!TConst("Maybe")

println("Ok!")


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using FMIImport.FMICore: fmi2VariableNamingConventionStructured

myFMU = fmi2Load("SpringFrictionPendulum1D", ENV["EXPORTINGTOOL"], ENV["EXPORTINGVERSION"])
myFMU.executionConfig.assertOnWarning = true

@test fmi2GetVersion(myFMU) == "2.0"
@test fmi2GetTypesPlatform(myFMU) == "default"

@test fmi2GetModelName(myFMU) == "SpringFrictionPendulum1D"
@test fmi2GetVariableNamingConvention(myFMU) == fmi2VariableNamingConventionStructured
@test fmi2IsCoSimulation(myFMU) == true
@test fmi2IsModelExchange(myFMU) == true

@test fmi2GetGUID(myFMU) == "{2e178ad3-5e9b-48ec-a7b2-baa5669efc0c}"
@test fmi2GetGenerationTool(myFMU) == "Dymola Version 2022x (64-bit), 2021-10-08"
@test fmi2GetGenerationDateAndTime(myFMU) == "2022-05-19T06:54:12Z"
@test fmi2GetNumberOfEventIndicators(myFMU) == 24
@test fmi2CanGetSetState(myFMU) == true
@test fmi2CanSerializeFMUstate(myFMU) == true
@test fmi2ProvidesDirectionalDerivative(myFMU) == true
@test fmi2DependenciesSupported(myFMU.modelDescription) == true
@test fmi2DerivativeDependenciesSupported(myFMU.modelDescription) == true

@test fmi2GetDefaultStartTime(myFMU.modelDescription) ≈ 0.0
@test fmi2GetDefaultStopTime(myFMU.modelDescription) ≈ 1.0
@test fmi2GetDefaultTolerance(myFMU.modelDescription) ≈ 1e-4
@test fmi2GetDefaultStepSize(myFMU.modelDescription) === nothing

# comfort getters (dictionaries)

@test length(fmi2GetValueReferencesAndNames(myFMU.modelDescription)) == 42
@test length(fmi2GetValueReferencesAndNames(myFMU)) == 42

@test length(fmi2GetInputNames(myFMU.modelDescription)) == 0
@test length(fmi2GetInputNames(myFMU)) == 0

@test length(fmi2GetOutputNames(myFMU.modelDescription)) == 0
@test length(fmi2GetOutputNames(myFMU)) == 0

@test length(fmi2GetParameterNames(myFMU.modelDescription)) == 0
@test length(fmi2GetParameterNames(myFMU)) == 0

@test length(fmi2GetStateNames(myFMU.modelDescription)) == 2
@test length(fmi2GetStateNames(myFMU)) == 2
@test fmi2GetStateNames(myFMU; mode=:first) == ["mass.s", "mass.v"]
@test fmi2GetStateNames(myFMU; mode=:flat) == ["mass.s", "mass.v", "mass.v_relfric"]
@test fmi2GetStateNames(myFMU; mode=:group) == [["mass.s"], ["mass.v", "mass.v_relfric"]]

@test length(fmi2GetDerivativeNames(myFMU.modelDescription)) == 2
@test length(fmi2GetDerivativeNames(myFMU)) == 2
@test fmi2GetDerivativeNames(myFMU; mode=:first) == ["der(mass.s)", "mass.a_relfric"]
@test fmi2GetDerivativeNames(myFMU; mode=:flat) == ["der(mass.s)", "mass.a_relfric", "mass.a", "der(mass.v)"]
@test fmi2GetDerivativeNames(myFMU; mode=:group) == [["der(mass.s)"], ["mass.a_relfric", "mass.a", "der(mass.v)"]]

@test length(fmi2GetNamesAndDescriptions(myFMU.modelDescription)) == 50
@test length(fmi2GetNamesAndDescriptions(myFMU)) == 50

@test length(fmi2GetNamesAndUnits(myFMU.modelDescription)) == 50
dict = fmi2GetNamesAndUnits(myFMU)
@test length(dict) == 50
@test dict["der(mass.s)"] == "m/s"
@test dict["mass.F_prop"] == "N.s/m"
@test dict["mass.fexp"] == "s/m"
@test dict["der(mass.v)"] == "m/s2"

@test length(fmi2GetNamesAndInitials(myFMU.modelDescription)) == 50
dict = fmi2GetNamesAndInitials(myFMU)
@test length(dict) == 50
@test dict["mass.startForward"] == 0
@test dict["mass.startBackward"] == 0
@test dict["mass.locked"] == 1

fmi2Unload(myFMU)


text:
using ReinforcementLearningBase
using Test

using CommonRLInterface
const CRL = CommonRLInterface

using POMDPs
using POMDPModelTools: Deterministic

@testset "ReinforcementLearningBase" begin
include("CommonRLInterface.jl")
end


text:
using Symata, NumericalMethodsforEngineers
using Test

@sym begin
  ClearAll(xi, yi, Y, Ydot, Ydotdot, ytilde1, ytilde2, R)
  xi = [0, 1//3, 2//3, 1]
  yi = [0, a, b, 1]
  Y(x_) := LagrangePolynomial(xi, yi)
  #
  # Can be re-formulated as Y(x) = F(x) + C1(a) * Ψ1(x) + C2(a) * Ψ2(x)
  #
  # F(x_) := (1/2)*(9x^3 - 9x^2 + 2x)
  # C1(a_) := (27/2)*(a-b)
  Ψ1(x_) := x^3 - x^2
  # C2(a_) := -(9/2)*(2a-b)
  Ψ2(x_) := x^2 - x
  #
  # Y(x_) := F(x) + C1(a)*Ψ1(x) + C2(a)*Ψ2(x)
  #
  R(x_) = Simplify(D(Y(x), x, 2) - 3*x - 4*Y(x))
  SetJ(r, ToString(Simplify(R(x))))
end

println("\n\nExample 7.14: y'' = 3x + 4y, y(0)=0, y(1)=1")
println("Residual for Weighted Residual Method using 2 point Lagragian Polynomial")
@sym Println("\nY(x): ", Expand(Y(x)), "\n")
println("( Example 7.14 gives: ytilde = 1/2*x^3*(27a-27b+9)-x^2*(45a-36b+9)+x*(18a-9b+2) )")
@sym Println("\nR(x) = ", R(x), "\n")

@test r == "-9 - 45a + 36b + 20x - 63b*x + 45a*x + 18*x^2 - 72b*x^2 + 90a*x^2 - 18x^3 - 54a*x^3 + 54b*x^3"


text:
#===============================================================================
# Methods are how functions are implemented over certain types. These are some
# exercises to understand their behavior
#
# Author: Carl Cortright
# Date: 12/30/2016
#
===============================================================================#

#
# Defines the function f for the specific method with argument types of Float64
#
function f(x::Float64, y::Float64)
  return 2x + y
end

#
# Now we define it over more general arguments
#
function f(x::Number, y::Number)
  return 2x - y
end

# Note: the first declaration creates the function declaration, adding creates methods
methods(f)
methods(+)
methods(println)

# it is important to define methods that clear up any ambiguity
g(x::Float64, y) = x + y
g(x, y::Float64) = x - y
g(x::Float64, y::Float64) = x*y

# Parametric functions allow us to define over many types without the catchall
sametype{T}(x::T, y::T) = true
sametype(x,y) = false

sametype("hello", 2)

# Vararg functions can also be parameterized
function bar(a, b, x::Vararg{Any, 2})
  return [a,b,x[1],x[2]]
end

bar(1,2,3,4)

# Functors are functions that are specified by type and not name
type polynomial{R}
  coeffs::Vector{R}
end
function (p::polynomial)(x)
  v = p.coeffs[end]
  for i in (length(p.coeffs)-1):-1:1
    v = v*x + p.coeffs[i]
  end
  return v
end

p = polynomial([1,2,3,4,5,100,10000])
p(10)
p(1)
p(101)

# empty functions can be used to separate interfaces and definitions
function emptyfunction
end

max_issues_repo_name: carlcortright/Playground


text:
import Random
export duplicaterng

# Relevant threads
# https://discourse.julialang.org/t/parallel-mersenne-twister/27567
# https://discourse.julialang.org/t/future-randjump-peculiar-speed/34862
# https://discourse.julialang.org/t/reproducible-multithreaded-monte-carlo-task-local-random/35269

"""
Create's `n` rngs from rng using rand jump

Related threads:

```jldoctest
import Random
seed = 3
rng = Random.MersenneTwister(seed)
rngs = duplicaterng(rng, 4)
julia> rand.(rngs)
4-element Array{Float64,1}:
 0.8116984049958615
 0.282977627325657
 0.13715166493608555
 0.1710542488647211
"""
function duplicaterng(rng::T, n) where T
  jmpamt = big(10)^20
  rngs = Vector{T}(undef, n)
  rngs[1] = rng
  for i in 2:n
    @inbounds rngs[i] = Future.randjump(rngs[i - 1], jmpamt)
  end
  rngs
end

duplicaterng(::Random._GLOBAL_RNG, n) =
  duplicaterng(Random.default_rng(), n)

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
"""
    FOLyapunov()

Computing fractional order Lyapunov exponent of a fractionl order system.

```tex
@article{Danca2018MatlabCF,
  title={Matlab Code for Lyapunov Exponents of Fractional-Order Systems},
  author={Marius-F. Danca and Nikolay V. Kuznetsov},
  journal={Int. J. Bifurc. Chaos},
  year={2018},
  volume={28},
  pages={1850067:1-1850067:14}
}
```
"""
function FOLyapunov(fun, order, t_start, h_norm, t_end, u0, h, out)# TODO: Generate the Lyapunov exponent plot
    ne::Int = length(u0) # System dimension

    Jfdefun(t, u) = jacobian_of_fdefun(fun, t, u)

    tspan = Float64[]
    LE = Float64[]

    # Generate extend system with jacobian
    function extend_fun(t, temp)
        temp=reshape(temp, ne, ne+1)
        result = zeros(ne)
        fun(result, temp[:, 1], nothing, t)
        for i=2:ne+1
            result = [result; Jfdefun(t, temp[:, 1])*temp[:, i]]
        end
        return result
    end
    x = zeros(Float64, ne*(ne+1))
    x0 = zeros(Float64, ne*(ne+1))
    c = zeros(Float64, ne)
    gsc = zeros(Float64, ne)
    zn = zeros(Float64, ne)
    n_it = round(Int, (t_end-t_start)/h_norm)
    x[1:ne] = u0
    q = order*ones(ne*(ne+1))# fractional order of the extend system
    for i=1:ne
        x[(ne+1)*i]=1.0
    end
    t = t_start
    LExp = zeros(ne)
    for it=1:n_it
        (_, Y) = pc(q, extend_fun, t, t+h_norm, x, h) # Solve the extend system
        t = t+h_norm
        Y = Y'

        x = Y[size(Y, 1), :]
        for i=1:ne
            for j=1:ne
                x0[ne*i+j]=x[ne*j+i]
            end
        end
        zn[1] = 0.0
        for j=1:ne
            zn[1] = zn[1]+x0[ne*j+1]^2
        end
        zn[1] = sqrt(zn[1])
        for j=1:ne
            x0[ne*j+1] = x0[ne*j+1]/zn[1]
        end
        for j=2:ne
            for k=1:(j-1)
                gsc[k] = 0.0
                for l=1:ne
                    gsc[k] = gsc[k]+x0[ne*l+j]*x0[ne*l+k]
                end
            end
            for k=1:ne
                for l=1:j-1
                    x0[ne*k+j]=x0[ne*k+j]-gsc[l]*x0[ne*k+l]
                end
            end
            zn[j]=0.0
            for k=1:ne
                zn[j]=zn[j]+x0[ne*k+j]^2
            end
            zn[j]=sqrt(zn[j])
            for k=1:ne
                x0[ne*k+j]=x0[ne*k+j]/zn[j]
            end
        end
        for k=1:ne
            c[k]=c[k]+log(zn[k])
            LExp[k]=c[k]/(t-t_start)
        end

        mod(it, out)==0 ? println(LExp) : nothing

        for i=1:ne
            for j=1:ne
                x[ne*j+i]=x0[ne*i+j]
            end
        end
        LE = [LE; LExp]
        tspan = [tspan; t]
    end
    LE = reshape(LE, ne, :)
    return LE, tspan
end

FOLyapunov(sys::FODESystem, h_norm, h, out) = FOLyapunov(sys.f, sys.α, tspan[1], h_norm, tspan[2], sys.x0, h, out)

mutable struct M
    an
    bn
    a0
    halpha1
    halpha2
    mu
    mu_tol
    r
    index_fft
    an_fft
    bn_fft
end
#TODO: Decouple ABM methods for FODESystems from FractionalSystems.jl to FractionalDiffEq.jl
function pc(alpha, f_fun, t0, T, y0, h)
    METH = M(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
    mu_tol = 1.0e-6
    mu = 1
    alpha = alpha[:]

    # Check compatibility size of the problem with number of fractional orders
    alpha_length = length(alpha)
    problem_size = size(y0, 1)

    m_alpha = ceil.(Int, alpha)
    m_alpha_factorial = zeros(alpha_length, maximum(m_alpha))
    for i = 1 : alpha_length
        for j = 0 : m_alpha[i]-1
            m_alpha_factorial[i, j+1] = factorial(j)
        end
    end

    f_temp = f_vectorfield(t0, y0[:, 1], f_fun)

    r = 16
    N = ceil(Int64, (T-t0)/h)
    Nr = ceil(Int64, (N+1)/r)*r
    Qr::Int = ceil(Int64, log2(Nr/r)) - 1
    NNr = Int64(2^(Qr+1)*r)

    # Preallocation of some variables
    y = zeros(problem_size, N+1)
    fy = zeros(problem_size, N+1)
    zn_pred = zeros(problem_size, NNr+1)
    if mu > 0
        zn_corr = zeros(problem_size, NNr+1)
    else
        zn_corr = 0
    end

    # Evaluation of coefficients of the PECE method
    nvett = 0 : NNr+1
    bn = zeros(alpha_length, NNr+1) ; an = copy(bn); a0 = copy(bn)
    for i_alpha = 1:alpha_length
        #find_alpha = find(alpha[i_alpha]==alpha[1:i_alpha-1])
        find_alpha=[]
        if alpha[i_alpha] == alpha[1:i_alpha-1]
            push!(find_alpha, i_alpha)
        end

        if isempty(find_alpha) == false
            bn[i_alpha, :] = bn[find_alpha[1], :]
            an[i_alpha, :] = an[find_alpha[1], :]
            a0[i_alpha, :] = a0[find_alpha[1], :]
        else
            nalpha = nvett.^alpha[i_alpha]
            nalpha1 = nalpha.*nvett
            bn[i_alpha, :] = nalpha[2:end] - nalpha[1:end-1]
            an[i_alpha, :] = [1; (nalpha1[1:end-2] - 2*nalpha1[2:end-1] + nalpha1[3:end]) ]
            a0[i_alpha, :] = [0; nalpha1[1:end-2]-nalpha[2:end-1].*(nvett[2:end-1].-alpha[i_alpha].-1)]
        end
    end
    METH.bn = bn ; METH.an = an ; METH.a0 = a0
    METH.halpha1 = h.^alpha./gamma.(alpha.+1)
    METH.halpha2 = h.^alpha./gamma.(alpha.+2)
    METH.mu = mu ; METH.mu_tol = mu_tol

    # Evaluation of FFT of coefficients of the PECE method
    METH.r = r
    if Qr >= 0
        index_fft = zeros(Int, 2, Qr+1) ;
        for l = 1 : Qr+1
            if l == 1
                index_fft[1, l] = 1 ; index_fft[2, l] = r*2
            else
                index_fft[1, l] = index_fft[2, l-1]+1 ; index_fft[2, l] = index_fft[2, l-1]+2^l*r
            end
        end

        bn_fft = zeros(Complex, alpha_length, index_fft[2, Qr+1]) ; an_fft = copy(bn_fft)
        for l = 1 : Qr+1
            coef_end = 2^l*r
            for i_alpha = 1 : alpha_length
                #find_alpha = find(alpha(i_alpha)==alpha(1:i_alpha-1)) ;
                find_alpha = []
                if alpha[i_alpha] == alpha[1:i_alpha-1]
                    push!(find_alpha, i_alpha)
                end
                if isempty(find_alpha) == false
                    bn_fft[i_alpha, index_fft[1, l]:index_fft[2, l]] = bn_fft(find_alpha(1),index_fft(1,l):index_fft(2,l)) ;
                    an_fft[i_alpha, index_fft[1, l]:index_fft[2, l]] = an_fft(find_alpha(1),index_fft(1,l):index_fft(2,l)) ;
                else
                    bn_fft[i_alpha, index_fft[1, l]:index_fft[2, l]] = ourfft(METH.bn[i_alpha, 1:coef_end], coef_end)
                    an_fft[i_alpha, index_fft[1, l]:index_fft[2, l]] = ourfft(METH.an[i_alpha, 1:coef_end], coef_end)
                end
            end
        end
        METH.bn_fft = bn_fft ; METH.an_fft = an_fft ; METH.index_fft = index_fft ;
    end

    # Initializing solution and proces of computation
    t = t0 .+ collect(0:N)*h
    y[:, 1] = y0[:, 1]
    fy[:, 1] = f_temp
    (y, fy) = Triangolo(1, r-1, t, y, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, y0, t0, f_fun) ;

    # Main process of computation by means of the FFT algorithm
    ff = zeros(1, 2^(Qr+2)) ; ff[1:2] = [0; 2] ; card_ff = 2
    nx0 = 0 ; ny0 = 0 ;
    for qr = 0 : Qr
        L = 2^qr
        (y, fy) = DisegnaBlocchi(L, ff, r, Nr, nx0+L*r, ny0, t, y, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, y0, t0, f_fun) ;
        ff[1:2*card_ff] = [ff[1:card_ff]; ff[1:card_ff]]
        card_ff = 2*card_ff
        ff[card_ff] = 4*L
    end

    # Evaluation solution in T when T is not in the mesh
    if T < t[N+1]
        c = (T - t[N])/h
        t[N+1] = T
        y[:, N+1] = (1-c)*y[:, N] + c*y[:, N+1]
    end

    t = t[1:N+1]; y = y[:, 1:N+1]
    return t, y

end

function DisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn_pred, zn_corr, N , METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, y0, t0, f_fun)

    nxi = nx0 ; nxf = nx0 + L*r - 1
    nyi = ny0 ; nyf = ny0 + L*r - 1
    is = 1
    s_nxi = zeros(N)
    s_nxf = zeros(N)
    s_nyi = zeros(N)
    s_nyf = zeros(N)
    s_nxi[is] = nxi ; s_nxf[is] = nxf ; s_nyi[is] = nyi ; s_nyf[is] = nyf ;

    i_triangolo = 0 ; stop = false
    while stop == false

        stop = (nxi+r-1 == nx0+L*r-1) || (nxi+r-1>=Nr-1) # Ci si ferma quando il triangolo attuale finisce alla fine del quadrato

        (zn_pred, zn_corr) = Quadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length) ;

        (y, fy) = Triangolo(nxi, nxi+r-1, t, y, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, y0, t0, f_fun) ;
        i_triangolo = i_triangolo + 1 ;

        if stop == false
            if nxi+r-1 == nxf   # Il triangolo finisce dove finisce il quadrato -> si scende di livello
                i_Delta = ff[i_triangolo]
                Delta = i_Delta*r ;
                nxi = s_nxf[is]+1 ; nxf = s_nxf[is]  + Delta ;
                nyi = s_nxf[is] - Delta +1; nyf = s_nxf[is]  ;
                s_nxi[is] = nxi ; s_nxf[is] = nxf ; s_nyi[is] = nyi ; s_nyf[is] = nyf ;
            else # Il triangolo finisce prima del quadrato -> si fa un quadrato accanto
                nxi = nxi + r ; nxf = nxi + r - 1 ; nyi = nyf + 1 ; nyf = nyf + r  ;
                is = is + 1 ;
                s_nxi[is] = nxi ; s_nxf[is] = nxf ; s_nyi[is] = nyi ; s_nyf[is] = nyf ;
            end
        end

    end
    return y, fy
end

function Quadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length)

    #coef_beg = nxi-nyf ;
    coef_end = Int64(nxf-nyi+1)
    i_fft = log2(coef_end/METH.r)
    funz_beg = Int64(nyi+1) ; funz_end = Int64(nyf+1)
    Nnxf = min(N, nxf)

    # Evaluation convolution segment for the predictor
    vett_funz = fy[:, funz_beg:funz_end]
    vett_funz_fft = rowfft(vett_funz, coef_end)
    zzn_pred = zeros(problem_size,coef_end) ;
    for i = 1 : problem_size
        i_alpha = min(alpha_length,i) ;
        Z = METH.bn_fft[i_alpha, METH.index_fft[1, Int64(i_fft)]:METH.index_fft[2, Int64(i_fft)]].*vett_funz_fft[i, :]
        zzn_pred[i, :] = real.(ourifft(Z, coef_end))
    end
    zzn_pred = zzn_pred[:, Int64(nxf-nyf+1-1):end-1]
    zn_pred[:, Int64(nxi+1):Int64(Nnxf+1)] = zn_pred[:, Int64(nxi+1):Int64(Nnxf+1)] + zzn_pred[:, 1:Int64(Nnxf-nxi+1)]

    # Evaluation convolution segment for the corrector
    if METH.mu > 0
        if nyi == 0 # Evaluation of the lowest square
            vett_funz = [zeros(problem_size, 1) fy[:, funz_beg+1:funz_end]]
            vett_funz_fft = rowfft(vett_funz, coef_end)
        end
        zzn_corr = zeros(problem_size, coef_end) ;
        for i = 1 : problem_size
            i_alpha = min(alpha_length,i) ;
            Z = METH.an_fft[i_alpha,METH.index_fft[1, Int64(i_fft)]:METH.index_fft[2, Int64(i_fft)]].*vett_funz_fft[i, :]
            zzn_corr[i, :] = real.(ourifft(Z, coef_end))
        end
        zzn_corr = zzn_corr[:, Int64(nxf-nyf+1):end]
        zn_corr[:, Int64(nxi+1):Int64(Nnxf+1)] = zn_corr[:, Int64(nxi+1):Int64(Nnxf+1)] + zzn_corr[:, 1:Int64(Nnxf-nxi+1)]
    else
        zn_corr = 0
    end
    return zn_pred, zn_corr
end

function Triangolo(nxi, nxf, t, y, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, y0, t0, f_fun)

    for n = Int64(nxi) : Int64(min(N, nxf))

        # Evaluation of the predictor
        Phi = zeros(problem_size, 1)
        if nxi == 1 # Case of the first triangle
            j_beg = 0
        else # Case of any triangle but not the first
            j_beg = nxi
        end
        for j = Int64(j_beg) : Int64(n-1)
            Phi = Phi + METH.bn[1:alpha_length,n-j].*fy[:, j+1]
        end
        St = StartingTerm(t[n+1], y0, m_alpha, t0, m_alpha_factorial)
        y_pred = St + METH.halpha1.*(zn_pred[:, n+1] + Phi)
        f_pred = f_vectorfield(t[n+1], y_pred, f_fun)

        # Evaluation of the corrector
        if METH.mu == 0
            y[:, n+1] = y_pred
            fy[:, n+1] = f_pred
        else
            j_beg = Int64(nxi)
            Phi = zeros(problem_size, 1)
            for j = j_beg : n-1
                Phi = Phi + METH.an[1:alpha_length, n-j+1].*fy[:, j+1]
            end
            Phi_n = St + METH.halpha2.*(METH.a0[1:alpha_length, n+1].*fy[:, 1] + zn_corr[:, n+1] + Phi)
            yn0 = y_pred ; fn0 = f_pred ;
            stop = 0 ; mu_it = 0 ;
            while stop == false
                global yn1 = Phi_n + METH.halpha2.*fn0 ;
                mu_it = mu_it + 1 ;
                if METH.mu == Inf
                    stop = norm(yn1-yn0,inf) < METH.mu_tol ;
                    if mu_it > 100 && ~stop
                        stop = 1 ;
                    end
                else
                    stop = mu_it == METH.mu ;
                end
                global fn1 = f_vectorfield(t[n+1], yn1, f_fun)
                yn0 = yn1 ; fn0 = fn1 ;
            end
            y[:, n+1] = yn1
            fy[:, n+1] = fn1
        end
    end
    return y, fy
end

function f_vectorfield(t, y, f_fun!)
    f = f_fun!(t, y)
    return f
end

function  StartingTerm(t, y0, m_alpha, t0, m_alpha_factorial)
    ys = zeros(size(y0, 1),1)
    for k = 1 : maximum(m_alpha)
        if length(m_alpha) == 1
            ys = ys + (t-t0)^(k-1)/m_alpha_factorial[k]*y0[:, k]
        else
            i_alpha = findall(x -> x>=k, m_alpha)
            ys[i_alpha, 1] = ys[i_alpha, 1] + (t-t0)^(k-1)*y0[i_alpha, k]./m_alpha_factorial[i_alpha, k]
        end
    end
    return ys
end

function ourfft(x::Vector, n)
    s=length(x)
    x=x[:]
    if s > n
        return fft(x[1:n])
    elseif s < n
        return fft([x; zeros(n-s)])
    else
        return fft(x)
    end
end

function ourifft(x, n)
    s=length(x)
    x=x[:]
    if s > n
        return ifft(x[1:n])
    elseif s < n
        return ifft([x; zeros(n-s)])
    else
        return ifft(x)
    end
end

function rowfft(x::AbstractMatrix, n)
    result = zeros(Complex, size(x)[1], n)
    for i=1:size(x)[1]
        result[i, :] = ourfft(x[i, :], n)
    end
    return result
end

function jacobian_of_fdefun(f, t, y)
    ForwardDiff.jacobian(y) do y
    du = similar(y)
    f(du, y, nothing, t)
    du
    end
end


text:
using Documenter, DescriptorSystems
DocMeta.setdocmeta!(DescriptorSystems, :DocTestSetup, :(using DescriptorSystems); recursive=true)

makedocs(
  modules  = [DescriptorSystems],
  sitename = "DescriptorSystems.jl",
  authors  = "Andreas Varga",
  format   = Documenter.HTML(prettyurls = false),
  pages    = [
     "Home"   => "index.md",
     "Library" => [
        "dss.md",
        "rtf.md",
        "operations.md",
        "operations_rtf.md",
        "conversions.md",
        "order_reduction.md",
        "analysis.md",
        "factorizations.md",
        "advanced_operations.md",
        "model_matching.md"
        ],
     "Utilities" => [
        "dstools.md"
     ],
     "Index" => "makeindex.md"
  ]
)

deploydocs(
  repo = "github.com/andreasvarga/DescriptorSystems.jl.git",
  target = "build",
  devbranch = "main"
)


text:
using DiffEqFlux, CUDA, MLDatasets, OrdinaryDiffEq, Printf, Test
using Flux.Losses: logitcrossentropy
using Flux.Data: DataLoader
using MLDataUtils: LabelEnc, convertlabel, stratifiedobs

CUDA.allowscalar(false)

function loadmnist(batchsize = bs, train_split = 0.9)
    # Use MLDataUtils LabelEnc for natural onehot conversion
    onehot(labels_raw) = convertlabel(LabelEnc.OneOfK, labels_raw,
                                      LabelEnc.NativeLabels(collect(0:9)))
    # Load MNIST
    imgs, labels_raw = MNIST.traindata()
    # Process images into (H,W,C,BS) batches
    x_data = Float32.(reshape(imgs, size(imgs,1), size(imgs,2), 1, size(imgs,3)))
    y_data = onehot(labels_raw)
    (x_train, y_train), (x_test, y_test) = stratifiedobs((x_data, y_data),
                                                         p = train_split)
    return (
        # Use Flux's DataLoader to automatically minibatch and shuffle the data
        DataLoader(gpu.(collect.((x_train, y_train))); batchsize = batchsize,
                   shuffle = true),
        # Don't shuffle the test data
        DataLoader(gpu.(collect.((x_test, y_test))); batchsize = batchsize,
                   shuffle = false)
    )
end

# Main
const bs = 128
const train_split = 0.9
train_dataloader, test_dataloader = loadmnist(bs, train_split)

down = Chain(Conv((3, 3), 1=>64, relu, stride = 1), GroupNorm(64, 64),
             Conv((4, 4), 64=>64, relu, stride = 2, pad=1), GroupNorm(64, 64),
             Conv((4, 4), 64=>64, stride = 2, pad = 1)) |>gpu

dudt = Chain(Conv((3, 3), 64=>64, tanh, stride=1, pad=1),
             Conv((3, 3), 64=>64, tanh, stride=1, pad=1)) |>gpu

fc = Chain(GroupNorm(64, 64), x -> relu.(x), MeanPool((6, 6)),
           x -> reshape(x, (64, :)), Dense(64,10)) |> gpu

nn_ode = NeuralODE(dudt, (0.f0, 1.f0), Tsit5(),
                   save_everystep = false,
                   reltol = 1e-3, abstol = 1e-3,
                   save_start = false) |> gpu

function DiffEqArray_to_Array(x)
    xarr = gpu(x)
    return xarr[:,:,:,:,1]
end

# Build our over-all model topology
model = Chain(down,                 # (28, 28, 1, BS) -> (6, 6, 64, BS)
              nn_ode,               # (6, 6, 64, BS) -> (6, 6, 64, BS, 1)
              DiffEqArray_to_Array, # (6, 6, 64, BS, 1) -> (6, 6, 64, BS)
              fc)                   # (6, 6, 64, BS) -> (10, BS)

# To understand the intermediate NN-ODE layer, we can examine it's dimensionality
img, lab = train_dataloader.data[1][:, :, :, 1:1], train_dataloader.data[2][:, 1:1]

x_d = down(img)

# We can see that we can compute the forward pass through the NN topology
# featuring an NNODE layer.
x_m = model(img)

classify(x) = argmax.(eachcol(x))

function accuracy(model, data; n_batches = 10)
    total_correct = 0
    total = 0
    for (i, (x, y)) in enumerate(data)
        # Only evaluate accuracy for n_batches
        i > n_batches && break
        target_class = classify(cpu(y))
        predicted_class = classify(cpu(model(x)))
        total_correct += sum(target_class .== predicted_class)
        total += length(target_class)
    end
    return total_correct / total
end

# burn in accuracy
accuracy(model, train_dataloader)

loss(x, y) = logitcrossentropy(model(x), y)

# burn in loss
loss(img, lab)

opt = ADAM(0.05)
iter = 0

cb() = begin
    global iter += 1
    # Monitor that the weights do infact update
    # Every 10 training iterations show accuracy
    if iter % 10 == 1
        train_accuracy = accuracy(model, train_dataloader) * 100
        test_accuracy = accuracy(model, test_dataloader;
                                 n_batches = length(test_dataloader)) * 100
        @printf("Iter: %3d || Train Accuracy: %2.3f || Test Accuracy: %2.3f\n",
                iter, train_accuracy, test_accuracy)
    end
end

Flux.train!(loss, Flux.params(down, nn_ode.p, fc), train_dataloader, opt, cb = cb)
@test accuracy(model, test_dataloader; n_batches = length(test_dataloader)) > 0.8


text:
@testset "soc problems" begin
    @testset "soc1_test" begin
    #=
   primal
       max 0x + 1y + 1z
   s.t.
       x == 1          :w_4
       x >= ||(y,z)||  :w_1, w_2, w_3

   dual
       min -w_4
   s.t.
       w_2 == -1
       w_3 == -1
       w_1 + w_4 == 0
       w_1 >= ||(w_2, w_3)||
    =#
       primal_model = soc1_test()
       dual_model, primal_dual_map = dual_model_and_map(primal_model)

       @test MOI.get(dual_model, MOI.NumberOfVariables()) == 4
       list_of_cons =  MOI.get(dual_model, MOI.ListOfConstraints())
       @test Set(list_of_cons) == Set([
           (SAF{Float64}, MOI.EqualTo{Float64})
           (VVF, MOI.SecondOrderCone)
       ])
       @test MOI.get(dual_model, MOI.NumberOfConstraints{VVF, MOI.SecondOrderCone}()) == 1
       @test MOI.get(dual_model, MOI.NumberOfConstraints{SAF{Float64}, MOI.EqualTo{Float64}}()) == 3
       obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
       @test obj_type == SAF{Float64}
       obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
       @test MOI.constant(obj) == 0.0
       @test MOI.coefficient.(obj.terms) == [-1.0]

       eq_con1_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(2))
       eq_con1_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(2))
       @test MOI.coefficient.(eq_con1_fun.terms) == [1.0; 1.0]
       @test MOI.constant.(eq_con1_fun) == 0.0
       @test MOI.constant(eq_con1_set) == 0.0
       eq_con2_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(3))
       eq_con2_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(3))
       @test MOI.coefficient.(eq_con2_fun.terms) == [1.0]
       @test MOI.constant.(eq_con2_fun) == 0.0
       @test MOI.constant(eq_con2_set) == -1.0
       eq_con3_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(4))
       eq_con3_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(4))
       @test MOI.coefficient.(eq_con3_fun.terms) == [1.0]
       @test MOI.constant.(eq_con3_fun) == 0.0
       @test MOI.constant(eq_con3_set) == -1.0

       soc_con = MOI.get(dual_model, MOI.ConstraintFunction(), CI{VVF, MOI.SecondOrderCone}(1))
       @test soc_con.variables == VI.(1:3)

       primal_con_dual_var = primal_dual_map.primal_con_dual_var
       @test primal_con_dual_var[CI{VAF{Float64}, MOI.Zeros}(1)] == [VI(4)]
       @test primal_con_dual_var[CI{VVF, MOI.SecondOrderCone}(2)] == [VI(1); VI(2); VI(3)]

       primal_var_dual_con = primal_dual_map.primal_var_dual_con
       @test primal_var_dual_con[VI(1)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(2)
       @test primal_var_dual_con[VI(2)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(3)
       @test primal_var_dual_con[VI(3)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(4)
   end

   @testset "soc2_test" begin
    #=
   primal
       max 0x + 1y + 1z
   s.t.
       x == 1          :w_4
       x >= ||(y,z)||  :w_1, w_2, w_3

   dual
       min -w_4
   s.t.
       w_2 == -1
       w_3 == -1
       w_1 + w_4 == 0
       w_1 >= ||(w_2, w_3)||
   =#
       primal_model = soc2_test()
       dual_model, primal_dual_map = dual_model_and_map(primal_model)

       @test MOI.get(dual_model, MOI.NumberOfVariables()) == 4
       list_of_cons =  MOI.get(dual_model, MOI.ListOfConstraints())
       @test Set(list_of_cons) == Set([
           (SAF{Float64}, MOI.EqualTo{Float64})
           (VVF, MOI.SecondOrderCone)
       ])
       @test MOI.get(dual_model, MOI.NumberOfConstraints{VVF, MOI.SecondOrderCone}()) == 1
       @test MOI.get(dual_model, MOI.NumberOfConstraints{SAF{Float64}, MOI.EqualTo{Float64}}()) == 3
       obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
       @test obj_type == SAF{Float64}
       obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
       @test MOI.constant(obj) == 0.0
       @test MOI.coefficient.(obj.terms) == [-1.0]

       eq_con1_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(2))
       eq_con1_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(2))
       @test MOI.coefficient.(eq_con1_fun.terms) == [1.0; 1.0]
       @test MOI.constant.(eq_con1_fun) == 0.0
       @test MOI.constant(eq_con1_set) == 0.0
       eq_con2_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(3))
       eq_con2_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(3))
       @test MOI.coefficient.(eq_con2_fun.terms) == [1.0]
       @test MOI.constant.(eq_con2_fun) == 0.0
       @test MOI.constant(eq_con2_set) == -1.0
       eq_con3_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(4))
       eq_con3_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(4))
       @test MOI.coefficient.(eq_con3_fun.terms) == [1.0]
       @test MOI.constant.(eq_con3_fun) == 0.0
       @test MOI.constant(eq_con3_set) == -1.0

       soc_con = MOI.get(dual_model, MOI.ConstraintFunction(), CI{VVF, MOI.SecondOrderCone}(1))
       @test soc_con.variables == VI.(2:4)

       primal_con_dual_var = primal_dual_map.primal_con_dual_var
       @test primal_con_dual_var[CI{VAF{Float64}, MOI.Zeros}(1)] == [VI(1)]
       @test primal_con_dual_var[CI{VAF{Float64}, MOI.SecondOrderCone}(2)] == [VI(2); VI(3); VI(4)]

       primal_var_dual_con = primal_dual_map.primal_var_dual_con
       @test primal_var_dual_con[VI(1)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(2)
       @test primal_var_dual_con[VI(2)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(3)
       @test primal_var_dual_con[VI(3)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(4)
   end
end


text:
# File iteration
Tables.rowaccess(::Type{<:File}) = true
Tables.rows(f::File) = f

struct Row
    file::File
    row::Int
end

getfile(r::Row) = getfield(r, :file)
getrow(r::Row) = getfield(r, :row)

Base.propertynames(r::Row) = getnames(getfile(r))

function Base.show(io::IO, r::Row)
    println(io, "CSV.Row($(getrow(r))) of:")
    show(io, getfile(r))
end

Base.eltype(f::File) = Row
Base.length(f::File) = getrows(f)

@inline function Base.iterate(f::File, st=1)
    st > length(f) && return nothing
    return Row(f, st), st + 1
end

@noinline badcolumnerror(name) = throw(ArgumentError("`$name` is not a valid column name"))

@inline function Base.getproperty(row::Row, name::Symbol)
    f = getfile(row)
    i = findfirst(x->x===name, getnames(f))
    i === nothing && badcolumnerror(name)
    return getcell(f, gettypes(f)[i], i, getrow(row))
end

@inline Base.getproperty(row::Row, ::Type{T}, col::Int, name::Symbol) where {T} =
    getcell(getfile(row), T, col, getrow(row))

# internal method for getting a cell value
@inline function getcell(f::File, ::Type{T}, col::Int, row::Int) where {T}
    ind = metaind(row)
    @inbounds offlen = gettape(f, col)[ind]
    missingvalue(offlen) && return missing
    return getvalue(Base.nonmissingtype(T), f, ind, offlen, col)
end

function getvalue(::Type{T}, f, indexoffset, offlen, col) where {T}
    @inbounds x = reinterp_func(T)(gettape(f, col)[indexoffset + 1])
    return x
end

function getvalue(::Type{Float64}, f, indexoffset, offlen, col)
    @inbounds x = gettape(f, col)[indexoffset + 1]
    return ifelse(intvalue(offlen), Float64(int64(x)), float64(x))
end

getvalue(::Type{Missing}, f, indexoffset, offlen, col) = missing

function getvalue(::Type{PooledString}, f, indexoffset, offlen, col)
    @inbounds x = gettape(f, col)[indexoffset + 1]
    @inbounds str = getrefs(f, col)[x]
    return str
end

function getvalue(::Type{String}, f, indexoffset, offlen, col)
    s = PointerString(pointer(getbuf(f), getpos(offlen)), getlen(offlen))
    return escapedvalue(offlen) ? unescape(s, gete(f)) : String(s)
end


text:
using LinearAlgebra: I

const POST_PRED_YVAL = -3.5
const POST_PRED_DIMENSION = 10
const POST_PRED_Y_OBSERVED = POST_PRED_YVAL * ones(POST_PRED_DIMENSION) / sqrt(POST_PRED_DIMENSION)
const POST_PRED_GROUND_TRUTH = pdf(
    MvNormal(POST_PRED_Y_OBSERVED / 2, I), -POST_PRED_Y_OBSERVED
)

@expectation function post_pred(y)
    x ~ MvNormal(zeros(length(y)), I)
    y ~ MvNormal(x, I)
    return pdf(MvNormal(-y, 0.5*I), x)
end


text:
using ..Parameters: transform_to_unconstrained, build_parameters_named_tuple

using Oceananigans.Utils: prettysummary

struct IterationSummary{P, M, C, V, E, O}
    parameters :: P     # constrained
    ensemble_mean :: M  # constrained
    ensemble_cov :: C   # constrained
    ensemble_var :: V
    mean_square_errors :: E
    objective_values :: O
    iteration :: Int
    pseudotime :: Float64
    pseudo_Δt :: Float64
end

"""
    eki_objective(eki, θ::AbstractVector, G::AbstractVector; constrained = false)

Given forward map `G` and parameters `θ`, return a tuple `(Φ₁, Φ₂)`
of terms in the ensemble Kalman inversion (`eki`) regularized objective
function, where,

```code
Φ = Φ₁ + Φ₂
```

`Φ₁` measures output misfit `½ || Γy^(-¹/₂) * (y .- G(θ)) ||²` and
`Φ₂` measures prior misfit `½ || Γθ^(-¹/₂) * (θ .- μθ) ||²`, where `y` is the observation
map, `G(θ)` is the forward map, `Γy` is the observation noise covariance, `Γθ` is
the prior covariance, and `μθ` represents the prior means. Note that `Γ^(-¹/₂) =
inv(sqrt(Γ))`.

When keyword argument `constrained` is provided with `true` then input `θ`
is assumed to represent constrained parameters.
"""
function eki_objective(eki, θ::AbstractVector, G::AbstractVector; constrained = false)
    y = eki.mapped_observations
    Γy = eki.noise_covariance

    fp = eki.inverse_problem.free_parameters
    priors = fp.priors
    unconstrained_priors = [unconstrained_prior(priors[name]) for name in fp.names]
    μθ = getproperty.(unconstrained_priors, :μ)
    Γθ = diagm( getproperty.(unconstrained_priors, :σ).^2 )

    if constrained
        θ = [transform_to_unconstrained(priors[name], θ[i])
                for (i, name) in enumerate(keys(priors))]
    end

    # Φ₁ = (1/2)*|| Γy^(-½) * (y - G) ||²
    Φ₁ = (1/2) * norm(inv(sqrt(Γy)) * (y .- G))^2
    # Φ₂ = (1/2)*|| Γθ^(-½) * (θ - μθ) ||²
    Φ₂ = (1/2) * norm(inv(sqrt(Γθ)) * (θ .- μθ))^2
    return (Φ₁, Φ₂)
end

"""
    IterationSummary(eki, X, forward_map_output=nothing)

Return the summary for ensemble Kalman inversion `eki`
with unconstrained parameters `X` and `forward_map_output`.
"""
function IterationSummary(eki, X, forward_map_output=nothing)
    priors = eki.inverse_problem.free_parameters.priors

    ensemble_mean = mean(X, dims=2)[:]
    constrained_ensemble_mean = transform_to_constrained(priors, ensemble_mean)

    ensemble_covariance = cov(X, dims=2)
    constrained_ensemble_covariance = inverse_covariance_transform(values(priors), X, ensemble_covariance)

    constrained_ensemble_variance = build_parameters_named_tuple(eki.inverse_problem.free_parameters,
                                                                 diag(constrained_ensemble_covariance))

    constrained_parameters = transform_to_constrained(priors, X)

    G = forward_map_output
    if !isnothing(forward_map_output)
        Nobs, Nens= size(forward_map_output)
        y = eki.mapped_observations
        mean_square_errors = [mapreduce((x, y) -> (x - y)^2, +, y, view(G, :, k)) / Nobs for k = 1:Nens]
    else
        mean_square_errors = nothing
    end

    # Vector of (Φ₁, Φ₂) pairs, one for each ensemble member at the current iteration
    objective_values = [eki_objective(eki, X[:, j], G[:, j]) for j in 1:size(G, 2)]

    return IterationSummary(constrained_parameters,
                            constrained_ensemble_mean,
                            constrained_ensemble_covariance,
                            constrained_ensemble_variance,
                            mean_square_errors,
                            objective_values,
                            eki.iteration,
                            eki.pseudotime,
                            eki.pseudo_Δt)
end

function finitefind(a, val, find)
    finite_a = deepcopy(a)
    finite_a[.!isfinite.(a)] .= val
    return find(finite_a)
end

finitefindmin(a) = finitefind(a, Inf, findmin)
finitefindmax(a) = finitefind(a, -Inf, findmax)

function Base.show(io::IO, is::IterationSummary)
    max_error, imax = finitefindmax(is.mean_square_errors)
    min_error, imin = finitefindmin(is.mean_square_errors)

    names = keys(is.ensemble_mean)
    parameter_matrix = [is.parameters[k][name] for name in names, k = 1:length(is.parameters)]
    min_parameters = minimum(parameter_matrix, dims=2)
    max_parameters = maximum(parameter_matrix, dims=2)

    print(io, summary(is), '\n')

    print(io, "                      ", param_str.(keys(is.ensemble_mean))..., '\n',
              "       ensemble_mean: ", param_str.(values(is.ensemble_mean))..., '\n',
              particle_str("best", is.mean_square_errors[imin], is.parameters[imin]), '\n',
              particle_str("worst", is.mean_square_errors[imax], is.parameters[imax]), '\n',
              "             minimum: ", param_str.(min_parameters)..., '\n',
              "             maximum: ", param_str.(max_parameters)..., '\n',
              "   ensemble_variance: ", param_str.(values(is.ensemble_var))...)

    return nothing
end

Base.summary(is::IterationSummary) = string("IterationSummary(",
                                            "iteration=", is.iteration,
                                            ", pseudotime=", prettysummary(is.pseudotime),
                                            ", pseudo_Δt=", prettysummary(is.pseudo_Δt), ") ",
                                            "for ", length(is.parameters),
                                            " particles and ", length(keys(is.ensemble_mean)),
                                            " parameters")

function param_str(p::Symbol)
    p_str = string(p)
    length(p_str) > 9 && (p_str = p_str[1:9])
    return @sprintf("% 10s | ", p_str)
end

param_str(p::Number) = @sprintf("% -1.3e | ", p)

particle_str(particle, error, parameters) =
    @sprintf("% 11s particle: ", particle) *
    string(param_str.(values(parameters))...) *
    @sprintf("error = %.6e", error)


text:
import NNlib.σ

# define some functions
loss(x,y) = x - y

function σ(x)
    x = 1+exp.(-x)
    for i = 1:length(x)
        x[i] = 1/x[i]
    end
    return x
end

function σ_derivative(x)
    res = 1-x
    for i = 1:length(x)
        res[i] *= x[i]
    end
    return res
end

function vecprod(a,b)
    if length(a) == length(b)
        res = reshape(zeros(length(a)),1,length(a))
        for i = 1:length(a)
            res[i] = a[i]*b[i]
        end
        return res
    end
end

function print_status(err, label, pred)
    println("Error: ", err)
    println("True: ", label)
    println("Pred: ", pred)
end

# input variables
bin_length = 8
max_num = 2^bin_length

α = 0.1                     # learning rate
input_dim = 2               # two input numbers
hidden_dim = 9              # hidden neurons
output_dim = 1              # the output is sum of the two numbers
set_size = 15000            # size of the training set

# training data
num_x = reshape([convert(Int16,rand(0:max_num/2)) for i = 1:set_size*2], 2, set_size)
num_y = reshape([num_x[1,i] + num_x[2,i] for i = 1:set_size], 1, set_size)

# convert to bin
x, y = [], []
for i = 1:set_size
    a = bin(num_x[1,i],bin_length)
    b = bin(num_x[2,i],bin_length)
    c = bin(num_y[i],bin_length)

    in_bin = zeros(2,bin_length)
    out_bin = zeros(1,bin_length)
    for j = 1:bin_length
        in_bin[1,j] = convert(Float64, a[j])-48
        in_bin[2,j] = convert(Float64, b[j])-48
        out_bin[j] = convert(Float64, c[j])-48
    end
    push!(x,in_bin)
    push!(y,out_bin)
end

# weight matrices
syn_0 = 2*randn(input_dim,hidden_dim)
syn_1 = 2*randn(hidden_dim,output_dim)
syn_r = 2*randn(hidden_dim,hidden_dim)

syn_0_update = zeros(size(syn_0))
syn_1_update = zeros(size(syn_1))
syn_r_update = zeros(size(syn_r))

# training
for i = 1:set_size

    prediction = reshape(zeros(bin_length),1,bin_length)          # result of the net
    overallError = 0

    layer_2_Δs = []
    layer_1_values = []
    push!(layer_1_values, reshape(zeros(hidden_dim),1,hidden_dim))

    # calculate output and error
    for pos = bin_length:-1:1
        # MODEL -->
        # hidden layer (input ~+ prev_hidden)
        layer_1 = σ(x[i][1:end, pos]'*syn_0 + layer_1_values[end]*syn_r)

        # output layer (new binary representation)
        layer_2 = σ(layer_1*syn_1)

        layer_2_error = loss(y[i][pos],layer_2)
        push!(layer_2_Δs, layer_2_error*σ_derivative(layer_2))     # save derivatives for each timestep
        overallError += abs(layer_2_error[1])

        prediction[pos] = round(layer_2[1,1])

        push!(layer_1_values, deepcopy(layer_1))                   # store hidden layer for next timestep
        # MODEL --<
    end

    future_layer_1_Δ = reshape(zeros(hidden_dim),1,hidden_dim)

    # backpropagation
    for pos = 1:bin_length
        layer_1 = layer_1_values[end-pos+1]        # selecting the current hidden layer
        prev_layer_1 = layer_1_values[end-pos]     # selecting the previous hidden layer

        # error at output layer
        layer_2_Δ = layer_2_Δs[end-pos+1]
        # error at hidden layer (given the error at the future hidden layer and the current output error)
        layer_1_Δ = vecprod((future_layer_1_Δ*syn_r' + layer_2_Δ*syn_1'),σ_derivative(layer_1))

        # let's update all our weights so we can try again
        syn_1_update += layer_1'*layer_2_Δ
        syn_r_update += prev_layer_1'*layer_1_Δ
        syn_0_update += x[i][1:end, pos]*layer_1_Δ

        future_layer_1_Δ = layer_1_Δ
    end

    # backpropagation is calculated, now update the weights
    syn_0 += syn_0_update * α
    syn_1 += syn_1_update * α
    syn_r += syn_r_update * α

    syn_0_update *= 0
    syn_1_update *= 0
    syn_r_update *= 0

    if i % (0.1*set_size) == 0
        print_status(overallError, y[i], prediction)
    end
end


text:
### Model Parameters

mutable struct DiscriminantParameters{T}
    "Model fit indicator - `true` if model has been fit"
    fit::Bool
    "Dimension along which observations are stored (1 for rows, 2 for columns)"
    dims::Int
    "Number of classes"
    m::Int
    "Number of predictors"
    p::Int
    "Matrix of class centroids (one per row if `dims` is 1 or 1 per column if `dims` is 2)"
    M::Matrix{T}
    "Prior-weighted overall centroid"
    μ::Vector{T}
    "Vector of class prior probabilities"
    π::Vector{T}
    "Vector of class counts"
    nₘ::Vector{Int}
    "Shrinkage parameter towards average eigenvalue of overall covariance matrix"
    γ::Union{Nothing,T}
    "Overall covariance matrix without `γ` regularization"
    Σ::Union{Nothing,Matrix{T}}
    function DiscriminantParameters{T}() where T
        new{T}(false, 0, 0, 0, Matrix{T}(undef,0,0), Vector{T}(undef,0), Vector{T}(undef,0),
               Vector{Int}(undef,0), nothing, nothing)
    end
end

function set_dimensions!(Θ::DiscriminantParameters{T},
                         y::Vector{<:Integer},
                         X::Matrix{T},
                         dims::Integer) where {T}
    n, p = check_dims(X, dims=dims)
    m = maximum(y)

    n₂ = length(y)
    n₂ == n || throw(DimensionMismatch("observation count along length of class index " *
                                       "vector y must match dimension $(dims) of data " *
                                       "matrix X (got $(n₂) and $(n))"))

    Θ.dims = dims
    Θ.m = m
    Θ.p = p

    return Θ
end

function set_gamma!(Θ::DiscriminantParameters{T}, γ::Union{Nothing,Real}) where T
    if γ !== nothing
        0 ≤ γ ≤ 1 || throw(DomainError(γ, "γ must be in the interval [0,1]"))
    end
    Θ.γ = γ

    return Θ
end

function set_statistics!(Θ::DiscriminantParameters{T},
                         y::Vector{<:Integer},
                         X::Matrix{T},
                         centroids::Union{Nothing,AbstractMatrix}) where T
    dims = Θ.dims
    m = Θ.m
    p = Θ.p

    # Compute centroids and class counts from data if not specified
    if centroids === nothing
        Θ.M = dims == 1 ? Matrix{T}(undef, m, p) : Matrix{T}(undef, p, m)
        Θ.nₘ = Vector{Int}(undef, m)

        class_statistics!(Θ.M, Θ.nₘ, X, y, dims=dims)
    else
        m₂, p₂ = check_dims(centroids, dims=dims)
        if m₂ != m
            throw(DimensionMismatch("class count along dimension $(dims) of centroid " *
                                    "matrix M must match maximum class index found in " *
                                    "class index vector y (got $(m₂) and $(m))"))
        elseif p₂ != p
            altdims = dims == 1 ? 2 : 1
            throw(DimensionMismatch("predictor count along dimension $(altdims) of " *
                                    "centroid matrix M must match dimension $(dims) of " *
                                    "data matrix X (got $(p₂) and $(p))"))
        end

        Θ.M = copyto!(similar(centroids, T), centroids)
        Θ.nₘ = class_counts!(Vector{Int}(undef, m), y)
    end

    n = 0
    for (k, nₖ) in enumerate(Θ.nₘ)
        n += nₖ
        nₖ > 1 || error("class count at class index $(k) must be greater than 1")
    end

    n₂ = size(X, dims)
    n == n₂ || error("class count total must observation count along dimension $(dims) " *
                     "of X (got $(n) and $(n₂))")

    return Θ
end

function set_priors!(Θ::DiscriminantParameters{T},
                     π::Union{Nothing,AbstractVector}) where T
    m = Θ.m

    # Compute priors from class frequencies in data if not specified
    if π === nothing
        Θ.π = broadcast!(/, Vector{T}(undef, m), Θ.nₘ, sum(Θ.nₘ))
    else
        m₂ = length(π)
        if m₂ != m
            throw(DimensionMismatch("class count along length of class prior probability " *
                                    "vector π must match maximum class index found in " *
                                    "class index vector y (got $(m₂) and $(m))"))
        end
        validate_priors(π)

        Θ.π = copyto!(similar(π, T), π)
    end

    return Θ
end


text:
using Pkg
Pkg.activate(string(@__DIR__, "/../../src/julia/"))

using Plots
using SQLite
using DataFrames

db = SQLite.DB(string(@__DIR__,"/../../data/KuraBenchXeon.db"))
_df = DataFrame(DBInterface.execute(db, "select * from experiments"))

df = _df[_df.solver .!= "radau", :]

# The database contains experiments on 30 systems of coupled ODEs
# Each system is defined on a unique network topology
# Systems 1:10  have 10 nodes
# Systems 11:20 have 100 nodes
# Systems 21:30 have 1000 nodes
# Every system is simulated with 10 different solver tolerances [0.001, 0.0001, ..., 10e-12]
# In total we end up with 300 different experiments
# These are repeated for 5 runtimes [scipy, jl_sparse, jl_nd, jitcode, Fortran]
#
# The reported results are averages of all systems with the same solver tolerances

const ntols = 10 # number of different tolerance settings

"""
      function system_size_mean(df, system_ids, runtime, column)

      Selects all rows in df that match the specified `runtime` and one of the `system_ids` and computes their mean of the given `column` with respect to a fixed solver tolerance
"""
function system_size_mean(df, system_ids, runtime, column)
      nsystems = length(system_ids)
      rows = df[(df.runtime .== runtime) .& map(x -> x ∈ system_ids, df.system), :]
      return sum(reshape(rows[:,column], ntols, nsystems), dims = 2) ./ nsystems
end

# Simple Test. Average of tolerance should give back the same tolerance.
@assert system_size_mean(df, 1:10, "jitcode", "rtol") ≈ df[1:10,:].rtol

# We produce different plots for each system size
ids_per_size = [1:10, 11:20, 21:30]
i = 1
saving = true

## Error vs tolerance
begin
      plot(system_size_mean(df, ids_per_size[i], "jitcode", "rtol"),
           system_size_mean(df, ids_per_size[i], "jitcode", "err_v_ref"),
           xlabel = "Requested tolerance",
           ylabel = "Error",
           scale = :log10,
           label = "jitcode",
           color_palette = :seaborn_colorblind6,
           legend = :bottomright)

      for runtime in ["networkdynamics.jl", "DifferentialEquations.jl", "scipy", "Fortran"]
            Plots.display(
                  plot!(system_size_mean(df, ids_per_size[i], runtime, "rtol"),
                        system_size_mean(df, ids_per_size[i], runtime, "err_v_ref"),
                        label = runtime))
      end
end

## WPD

for i = 1:3
      begin
            plot(system_size_mean(df, ids_per_size[i], "jitcode", "err_v_ref"),
                 system_size_mean(df, ids_per_size[i], "jitcode","tintegration"),
                 xlabel = "Error",
                 ylabel = "Integration time (s)",
                 title = "N = $(10^i)",
                 scale = :log10,
                 marker = :diamond,
                 markersize = 5,
                 label = "JiTCODE",
                 color_palette =:seaborn_colorblind6,
                 legend = :topright)

            runtime = "networkdynamics.jl"
            plot!(system_size_mean(df, ids_per_size[i], runtime, "err_v_ref"),
                  system_size_mean(df, ids_per_size[i], runtime, "tintegration"),
                  marker = :star5,
                  markersize = 5,
                  label = "NetworkDynamics.jl")

            runtime = "DifferentialEquations.jl"
            plot!(system_size_mean(df, ids_per_size[i], runtime, "err_v_ref"),
                  system_size_mean(df, ids_per_size[i], runtime, "tintegration"),
                  marker = :utriangle,
                  markersize = 5,
                  label = "SparseArrays.jl")

            runtime = "Fortran"
            Plots.display(plot!(system_size_mean(df, ids_per_size[i], runtime, "err_v_ref"),
                        system_size_mean(df, ids_per_size[i], runtime, "tintegration"),
                        marker = :x,
                        markersize = 5,
                        label = runtime))

            if i == 1
                  runtime = "scipy"

                  Plots.display(plot!(system_size_mean(df, ids_per_size[i], runtime, "err_v_ref"),
                              system_size_mean(df, ids_per_size[i], runtime, "tintegration"),
                              marker = :o,
                              markersize = 5,
                              label = "SciPy"))
            end
            saving ? savefig(string(@__DIR__, "/WPD$(10^i).pdf")) : nothing
      end
end
## WPD jit

begin

      plot(system_size_mean(df, ids_per_size[i], "jitcode", "err_v_ref"),
            system_size_mean(df, ids_per_size[i], "jitcode","tjit"),
            xlabel = "Error",
            ylabel = "Preparation + integration time (s)",
            scale = :log10,
            label = "JiTCODE",
            color_palette = :seaborn_colorblind6,
            legend = :topright)

      runtime = "networkdynamics.jl"
      plot!(system_size_mean(df, ids_per_size[i], runtime, "err_v_ref"),
            system_size_mean(df, ids_per_size[i], runtime, "tjit"),
            label = "NetworkDynamics.jl")

      runtime = "DifferentialEquations.jl"
      plot!(system_size_mean(df, ids_per_size[i], runtime, "err_v_ref"),
            system_size_mean(df, ids_per_size[i], runtime, "tjit"),
            label = "SparseArrays.jl")
end

# Scaling of compute time vs system size
# For Jitcode the inital improvement is porbably due to SymEngine replacing SymPy (?),
# the later increase due to increased cost of symbolic preprocessing

begin

      plot([10, 100, 1000],
           [system_size_mean(df, sys_ids, "jitcode","ttotal")[4] for sys_ids in ids_per_size],
           xlabel = "Number of nodes",
           ylabel = "Preparation + integration time (s)",
           scale = :log10,
           color_palette = :seaborn_colorblind6,
           label = "JiTCODE",
           marker = :diamond,
           markersize = 5,
           legend = :topleft)

      runtime = "networkdynamics.jl"
      plot!([10, 100, 1000],
            [system_size_mean(df, sids, runtime,"ttotal")[4] for sids in ids_per_size],
            marker = :star5,
            markersize = 5,
            label = "NetworkDynamics.jl")

      runtime = "DifferentialEquations.jl"
      plot!([10, 100, 1000],
            [system_size_mean(df, sids, runtime,"ttotal")[4] for sids in ids_per_size],
            marker = :utriangle,
            markersize = 5,
            label = "SparseArrays.jl")

      Plots.display(
            plot!([10, 100, 1000],
                  [system_size_mean(df, sids, "scipy","ttotal")[4] for sids in ids_per_size],
                  marker = :circle,
                  markersize = 5,
                  label = "SciPy"))

      saving ?  savefig(string(@__DIR__, "/size_vs_jit.pdf")) : nothing
end

max_stars_repo_name: PIK-ICoNe/NetworkDynamicsBenchmarks

max_issues_repo_name: PIK-ICoNe/NetworkDynamicsBenchmarks

max_forks_repo_name: PIK-ICoNe/NetworkDynamicsBenchmarks


text:
export OrderedPlan

"Solution type for fully ordered plans."
struct OrderedPlan <: OrderedSolution
    plan::Vector{Term}
end

get_action(sol::OrderedPlan, t::Int) = sol.plan[t]

Base.iterate(sol::OrderedPlan) = iterate(sol.plan)
Base.iterate(sol::OrderedPlan, istate) = iterate(sol.plan, istate)
Base.getindex(sol::OrderedPlan, i::Int) = getindex(sol.plan, i)
Base.length(sol::OrderedPlan) = length(sol.plan)