text:
include("EnsembleMethods.jl")

d_raw = readtable("../data/concrete_data.csv")
dset = as_floats(d_raw)

X = dset[:, 1:8]
y = convert(Array, dset[:, 9])

@time fm1 = build_forest_df(y, X, 3, 5; nthreads = 1, n_surrogates = 0);
Profile.clear_malloc_data()
@time fm1 = build_forest_df(y, X, 3, 50; nthreads = 1, n_surrogates = 0);

n = 1000
p = 20
X = DataFrame(randn(n, p));
y = rand([true, false], n);
wgt = ones(Int, n);

@code_warntype _split_classifcation_error_loss(y, X, collect(1:n), collect(1:p), wgt)

using DecisionTree

X = convert(Array, dset[:, 1:8])
y = convert(Array, dset[:, 9])
@time fm2 = build_forest(y, X, 3, 50);


text:
using ChebyshevQuantum
using Test
using SpecialFunctions

tol_var=1e-8

@testset "bessel eq test" begin

    a=0.0;
    b= 5.520078110588564; #second zero of besselj0

    #solve bessels equation between 0 and 2nd zero of besselj0 with Direchlet b.c.'s.
    t = ChebyshevQuantum.DiffEq.solve(bc1 = [:a,0, 1.0], bc2 = [:b,0,0.0], A0=x->x^2, A1= x->x, A2 = x->(x^2) ,B = 0.0, N=20,  a=a, b= b);

    x = a:.01:b

    @test maximum(abs.(besselj0.(x) - t.(x))) < tol_var
    @test maximum(abs.(x.^2 .* t''.(x) + x .* t'.(x) + (x.^2 ).*t.(x))) < tol_var #test diffeq directly

    a=1.0;
 #   b= 3.8319091796875; #first zero of besselj1
    b=   3.8317059218883514
    #solve bessels equation α = 1.0  between 1 and zero of besselj1 with Direchlet b.c.'s.
    t = ChebyshevQuantum.DiffEq.solve(bc1 = [:a,0, besselj1(a)], bc2 = [:b,0,besselj1(b)], A0=x->(x^2 - 1.0), A1= x->x, A2 = x->x^2  ,B = 0.0, N=20,  a=a, b= b);

    x = a:.01:b

    @test maximum(abs.(besselj1.(x) - t.(x))) < tol_var
    @test maximum(abs.(x.^2 .* t''.(x) + x .* t'.(x) + (x.^2 .- 1 ).*t.(x))) < tol_var #test diffeq directly

end

@testset "heat eq test" begin

    a=0.0;
    b= 1.0;

    #solve heat equation
    t = ChebyshevQuantum.DiffEq.solve(bc1 = [:a,0,1.0], bc2 = [:b,0,0.5], A2 = 0.5 , N=20,  a=a, b= b);

    x = a:.01:b

    function f(x) #solution is a line
        1.0 - 0.5*x
    end

    @test maximum(abs.(f.(x) - t.(x))) < tol_var

end

@testset "SHO test" begin

    a=0.0;
    b=2*pi;

    #solve simple harmonic oscillator equation, neuman b.c.'s
    t = ChebyshevQuantum.DiffEq.solve(bc1 = [:a,0,1.0], bc2 = [:a,1,0.0], A2 = 1.0, A0=1.0 , N=20,  a=a, b= b);

    x = a:.01:b

    function f(x) #solution is a cosine
        cos(x)
    end

    @test maximum(abs.(f.(x) - t.(x))) < tol_var

end

@testset "Decay Exp test" begin

    a=0.0;
    b=2*pi;

    #solve decaying exp equation, neuman b.c.'s
    t = ChebyshevQuantum.DiffEq.solve(bc1 = [:a,0,1.0], bc2 = [:a,1,-1.0], A2 = 1.0, A0=-1.0 , N=20,  a=a, b= b);

    x = a:.01:b

    function f(x) #solution is a decaying exp
        exp(-x)
    end

    @test maximum(abs.(f.(x) - t.(x))) < tol_var

end

@testset "First order Decay Exp" begin

    a=0.0;
    b=10.0;

    #solve decaying exp equation, first order
    t = ChebyshevQuantum.DiffEq.solve(bc1 = [:a,0,1.0],  A1 = 1.0, A0=1.0 , N=20,  a=a, b= b);

    x = a:.01:b

    function f(x) #solution is a decaying exp
        exp(-x)
    end

    @test maximum(abs.(f.(x) - t.(x))) < tol_var
    @test maximum(abs.(  t'.(x) .+ t.(x))) < tol_var #test diffeq directly

end

ext: jl


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>20.0, :alg=>GLGM06(δ=0.01)),
    "X0"   => Hyperrectangle(low=[fill(0.0002, 10); zeros(10903)],
                            high=[fill(0.00025, 10); zeros(10903)]),
    "info" => Dict(
        "name" => "mna5",
        "dim" => 10913,
        "linear" => true,
        "hybrid" => false
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
using Documenter, DistributedArrays

makedocs(
    modules = [DistributedArrays],
    format = :html,
    sitename = "DistributedArrays.jl",
    pages = [
        "Introduction" => "index.md"
        "API" => "api.md"
    ],
    doctest = true
)

deploydocs(
    repo = "github.com/JuliaParallel/DistributedArrays.jl.git",
    julia = "1.0",
    # no need to build anything here, re-use output of `makedocs`
    target = "build",
    deps = nothing,
    make = nothing
)


text:
using Revise

using RCall
using FastGroupBy, BenchmarkTools, DataFrames

const N = Int(1e8); const K = 100;

R"""
memory.limit(2e9)
library(data.table)
N=$N; K=$K
set.seed(1)
DT <- data.table(
  id_bool = sample(K, N, TRUE),                      # bool
  v3 =  sample(round(runif(100,max=100),4), N, TRUE) # numeric e.g. 23.5749
)
cat("GB =", round(sum(gc()[,2])/1024, 3), "\n")
pt = proc.time()
DT[, sum(v3), keyby=id_bool]
y = timetaken(pt)
rm(DT); gc()
"""

@rget y;

srand(1);
id_bool = rand(Bool, N);
val = rand(round.(rand(K)*100,4), N);
y1 = @elapsed fastby!(sum, id_bool, val);

y
y1

r_speed = parse(Float64, y[1:end-3])

fgb_speed = y1/parse(Float64, y[1:end-3])

srand(1);
using StatsBase;
id_bool = rand(Bool, N);
val = rand(round.(rand(K)*100,4), N);
y2 = @elapsed countmap(id_bool, weights(val));

cm_speed = y2/parse(Float64, y[1:end-3])

using Plots
bar([1, fgb_speed, cm_speed], title = "Bool fastby! Speed",
xaxis = (""))
hline!(1)


text:
function next_boolean(obj::IntRandomGenerator)
    return jcall(obj, "nextBoolean", jboolean, ())
end

function next_bytes(obj::IntRandomGenerator, arg0::Vector{jbyte})
    return jcall(obj, "nextBytes", void, (Vector{jbyte},), arg0)
end

function next_bytes(obj::IntRandomGenerator, arg0::Vector{jbyte}, arg1::jint, arg2::jint)
    return jcall(obj, "nextBytes", void, (Vector{jbyte}, jint, jint), arg0, arg1, arg2)
end

function next_double(obj::IntRandomGenerator)
    return jcall(obj, "nextDouble", jdouble, ())
end

function next_float(obj::IntRandomGenerator)
    return jcall(obj, "nextFloat", jfloat, ())
end

function next_long(obj::IntRandomGenerator)
    return jcall(obj, "nextLong", jlong, ())
end


text:
using Yao
using QuDiffEq
using LinearAlgebra
using BitBasis
using OrdinaryDiffEq
using Plots
using Test
# vertex and egde values
vertx = 7
ege = 8

#Constructing incidence matrix B
B = zeros(vertx,ege)
@inbounds for i in 1:vertx
    B[i,i] = -1
    B[i,i+1] = 1
end
B[1,1] = 1

de = 0.5
sn = sin.((0.0:de:(vertx-1)*de)*2*pi/((vertx-1)*de))
u0 = ComplexF64.(sn)

# Intial conditions (stationary to begin with)
u1 = [u0; zero(u0); 0.0; 0.0]
u_ = Float64[u0 zero(u0)]

k = 2   # order in Taylor expansion
t = 1e-2 # time step
B_t = transpose(B)
n = 11  # number of steps
a = 1e-1 # spactial discretization

function make_hamiltonian(B,B_t,a)
    vertx,ege = size(B)
    n = nextpow(2,vertx+ege)
    H = zeros(n,n)
    @inbounds H[1:vertx,vertx+1:vertx+ege] = B
    @inbounds H[vertx+1:vertx+ege,1:vertx] = B_t
    H = -im/a*H
    return H
end

function do_pde(ϕ,B,B_t,k,t,a,n)
    vertx, = size(B)
    H = make_hamiltonian(B,B_t,a)
    res = Array{Array{ComplexF64,1},1}(undef,n)
    res[1] = @view ϕ[1:vertx]
    for i in 2:n
    r, N = taylorsolve(H,ϕ,k,t)
    ϕ = N*vec(state(r))
    res[i] = @view ϕ[1:vertx]
    end
    res_real = real(res)
    return res_real
end

#Dirichlet
res1 = do_pde(u1,B,B_t,k,t,a,n)

const D1 = -1/(a*a)*B*B_t # Laplacian Operator
function f(du,u,p,t)
    buffer, D = p
    u1 = @view(u[:,1])
    u2 = @view(u[:,2])
    mul!(buffer, D, u2)
    Du = buffer

    du[:,1] = Du
    du[:,2] = u1
end

tspan = (0.0,0.1)
prob1 = ODEProblem(f,u_,tspan,(zero(u_[:,1]), D1))
sol1 = solve(prob1,Tsit5(),dt=0.01,adaptive = false)
s1 = Array{Array{Float64,1},1}(undef,n)
for i in 1:n
    s1[i] = @view (sol1.u[i][:,1])
end

@test isapprox.(res1,s1, atol = 1e-2) |> all


text:
############################################################
##Interface to the ABC-PMC algorithms. Mainly documentation.
############################################################

"
Perform a version of ABC-PMC: Algorithm 2 of the paper which is essentially the original Toni et al approach.
Its arguments are:

* `abcinput`
An `ABCInput` variable.

* `N`
Number of accepted particles in each iteration.

* `α`
A tuning parameter between 0 and 1 determining how fast the acceptance threshold is reduced. (In more detail, the acceptance threshold in iteration t is the α quantile of distances from particles which meet the acceptance criteria of the previous iteration.)

* `maxsims`
The algorithm terminates once this many simulations have been performed.

* `nsims_for_init`
How many simulations are stored to initialise the distance function (by default 10,000).

* `store_init` (optional)
Whether to store the parameters and simulations used to initialise the distance function in each iteration (useful to produce many of the paper's figures).
These are stored for every iteration even if `adaptive` is false.
By default false.

* `diag_perturb` (optional)
Whether perturbation should be based on a diagonalised estimate of the covariance.
By default this is false, giving the perturbation described in Section 2.2 of the paper.

* `silent` (optional)
When true no progress bar or messages are shown during execution.
By default false.

* `h1` (optional)
The acceptance threshold for the first iteration, by default `Inf` (accept everything).
If `initialise_dist` is true, then this must be left at its default value.

The output is a `ABCPMCOutput` object.
"
function abcPMC2(abcinput::ABCInput, N::Integer, α::Float64, maxsims::Integer, nsims_for_init=10000; store_init=false, diag_perturb=false, silent=false, h1=Inf)
    abcPMC_23(abcinput, N, α, maxsims, false, nsims_for_init, store_init=store_init, diag_perturb=diag_perturb, silent=silent, h1=h1)
end

"
Perform a version of ABC-PMC: Algorithm 3 of the paper which chooses the distance function in the initial iteration.
Its arguments are:

* `abcinput`
An `ABCInput` variable.

* `N`
Number of accepted particles in each iteration.

* `α`
A tuning parameter between 0 and 1 determining how fast the acceptance threshold is reduced. (In more detail, the acceptance threshold in iteration t is the α quantile of distances from particles which meet the acceptance criteria of the previous iteration.)

* `maxsims`
The algorithm terminates once this many simulations have been performed.

* `nsims_for_init`
How many simulations are stored to initialise the distance function (by default 10,000).

* `store_init` (optional)
Whether to store the parameters and simulations used to initialise the distance function in each iteration (useful to produce many of the paper's figures).
These are stored for every iteration even if `adaptive` is false.
By default false.

* `diag_perturb` (optional)
Whether perturbation should be based on a diagonalised estimate of the covariance.
By default this is false, giving the perturbation described in Section 2.2 of the paper.

* `silent` (optional)
When true no progress bar or messages are shown during execution.
By default false.

* `h1` (optional)
The acceptance threshold for the first iteration, by default `Inf` (accept everything).
If `initialise_dist` is true, then this must be left at its default value.

The output is a `ABCPMCOutput` object.
"
function abcPMC3(abcinput::ABCInput, N::Integer, α::Float64, maxsims::Integer, nsims_for_init=10000; store_init=false, diag_perturb=false, silent=false, h1=Inf)
    abcPMC_23(abcinput, N, α, maxsims, true, nsims_for_init, store_init=store_init, diag_perturb=diag_perturb, silent=silent, h1=h1)
end

"
Perform a version of ABC-PMC.
This is the variant of Algorithm 3 mentioned in Section 4.2 of the paper.
As for Algorithm 3, the distance function is chosen in the initial iteration.
The variation is that the acceptance threshold for iteration t is chosen during that iteration.
Its arguments are as follows:

* `abcinput`
An `ABCInput` variable.

* `N`
Number of accepted particles in each iteration.

* `α`
A tuning parameter between 0 and 1 determining how fast the acceptance threshold is reduced. (In more detail, the acceptance threshold in iteration t is the α quantile of distances from particles which meet the acceptance criteria of the previous iteration.)

* `maxsims`
The algorithm terminates once this many simulations have been performed.

* `nsims_for_init`
How many simulations are stored to initialise the distance function (by default 10,000).

* `store_init` (optional)
Whether to store the parameters and simulations used to initialise the distance function in each iteration (useful to produce many of the paper's figures).
These are stored for every iteration even if `adaptive` is false.
By default false.

* `diag_perturb` (optional)
Whether perturbation should be based on a diagonalised estimate of the covariance.
By default this is false, giving the perturbation described in Section 2.2 of the paper.

* `silent` (optional)
When true no progress bar or messages are shown during execution.
By default false.

The output is a `ABCPMCOutput` object.
"
function abcPMC3V(abcinput::ABCInput, N::Integer, α::Float64, maxsims::Integer, nsims_for_init=10000; store_init=false, diag_perturb=false, silent=false)
    abcPMC_3V5(abcinput, N, α, maxsims, false, nsims_for_init, store_init=store_init, diag_perturb=diag_perturb, silent=silent)
end

"
Perform a version of ABC-PMC proposed in the paper, Algorithm 5.
This chooses the distance function for iteration t during that iteration.
Its arguments are as follows:

* `abcinput`
An `ABCInput` variable.

* `N`
Number of accepted particles in each iteration.

* `α`
A tuning parameter between 0 and 1 determining how fast the acceptance threshold is reduced. (In more detail, the acceptance threshold in iteration t is the α quantile of distances from particles which meet the acceptance criteria of the previous iteration.)

* `maxsims`
The algorithm terminates once this many simulations have been performed.

* `nsims_for_init`
How many simulations are stored to initialise the distance function (by default 10,000).

* `store_init` (optional)
Whether to store the parameters and simulations used to initialise the distance function in each iteration (useful to produce many of the paper's figures).
These are stored for every iteration even if `adaptive` is false.
By default false.

* `diag_perturb` (optional)
Whether perturbation should be based on a diagonalised estimate of the covariance.
By default this is false, giving the perturbation described in Section 2.2 of the paper.

* `silent` (optional)
When true no progress bar or messages are shown during execution.
By default false.

The output is a `ABCPMCOutput` object.
"
function abcPMC5(abcinput::ABCInput, N::Integer, α::Float64, maxsims::Integer, nsims_for_init=10000; store_init=false, diag_perturb=false, silent=false)
    abcPMC_3V5(abcinput, N, α, maxsims, true, nsims_for_init, store_init=store_init, diag_perturb=diag_perturb, silent=silent)
end


text:
include("Continuum.jl")
export SmallStrainContinuum

@doc raw"""
    SmallStrainContinuum

Constructs a small strain element.

- `eledim`: spatial dimension of the element (default = 2).
- `mat`: constitutive law, a length `#elem` vector of materials such as [`PlaneStress`](@ref)
- `elnodes`: the node indices in this finite element, an integer array
- `props`: property dictionary
- `coords`: coordinates of the vertices of the element
- `dhdx`: list of `ngp` shape functions for first order derivatives $\nabla \phi(x)$ (`ndof×2`) on the Gaussian points
- `weights`: weight vector of length `n_gauss_points`, for numerical quadrature
- `hs`: list of `ngp` shape functions for function values $\phi(x)$ (length `ndof` vectors) on the Gaussian points
- `stress`: stress at each quadrature points; this field is reserved for visualization.

# Example
```julia
#   Local degrees of freedom
#   4 ---- 3
#
#   1 ---- 2

nx = 10
ny = 5
h = 0.1
element = SmallStrainContinuum[]
prop = Dict("name"=> "PlaneStrain", "rho"=> 0.0876584, "E"=>0.07180760098, "nu"=>0.4)
for j = 1:ny
    for i = 1:nx
        n = (nx+1)*(j-1) + (i-1)+1
        elnodes = [n, n + 1, n + 1 + (nx + 1), n + (nx + 1)]
        ngp = 3 # 3 x 3 Gauss points per element
        coords = [(i-1)*h (j-1)*h
                    i*h (j-1)*h
                    i*h j*h
                    (i-1)*h j*h]
        push!(element, SmallStrainContinuum(coords,elnodes, prop, ngp))
    end
end
```

"""
mutable struct SmallStrainContinuum <: Continuum
    eledim::Int64
    mat  # constitutive law
    elnodes::Array{Int64}   # the node indices in this finite element
    props::Dict{String, Any}
    coords::Array{Float64}
    dhdx::Array{Array{Float64}}  # 4nPointsx2 matrix
    weights::Array{Float64}
    hs::Array{Array{Float64}}
    stress::Array{Array{Float64}}
end

Base.show(io::IO, z::SmallStrainContinuum) = print(io, "SmallStrainContinuum with $(length(z.weights)) Gauss quadrature points and $(z.props["name"]) material.")

"""
    SmallStrainContinuum(coords::Array{Float64}, elnodes::Array{Int64}, props::Dict{String, Any}, ngp::Int64=2)
"""
function SmallStrainContinuum(coords::Array{Float64}, elnodes::Array{Int64}, props::Dict{String, Any}, ngp::Int64=2)
    eledim = 2
    # @show coords, ngp

    dhdx, weights, hs = get2DElemShapeData( coords, ngp )
    nGauss = length(weights)
    name = props["name"]
    if name=="PlaneStrain"
        mat = [PlaneStrain(props) for i = 1:nGauss]
    elseif name=="Scalar1D"
        mat = [Scalar1D(props) for i = 1:nGauss]
    elseif name=="ViscoelasticityMaxwell"
        mat = [ViscoelasticityMaxwell(props) for i = 1:nGauss]
    elseif name=="PlaneStrainViscoelasticityProny"
        mat = [PlaneStrainViscoelasticityProny(props) for i = 1:nGauss]
    elseif name=="PlaneStressViscoelasticityProny"
        mat = [PlaneStressViscoelasticityProny(props) for i = 1:nGauss]
    elseif name=="PlaneStress"
        mat = [PlaneStress(props) for i = 1:nGauss]
    elseif name=="PlaneStressPlasticityLawBased"
        mat = [PlaneStressPlasticityLawBased(props) for i = 1:nGauss]
    elseif name=="PlaneStressPlasticity"
        mat = [PlaneStressPlasticity(props) for i = 1:nGauss]
    elseif name=="NeuralNetwork2D"
        mat = [NeuralNetwork2D(props) for i = 1:nGauss]
    elseif name=="PlaneStressIncompressibleRivlinSaunders"
        mat = [PlaneStressIncompressibleRivlinSaunders(props) for i = 1:nGauss]
    else
        error("Not implemented yet: $name")
    end
    strain = Array{Array{Float64}}(undef, length(weights))
    SmallStrainContinuum(eledim, mat, elnodes, props, coords, dhdx, weights, hs, strain)
end

"""
    getStiffAndForce(elem::SmallStrainContinuum, state::Array{Float64}, Dstate::Array{Float64}, Δt::Float64)

Returns the internal force term and the stiffness matrix. `state` and `Dstate` are restriction of full state variables to this element.
"""
function getStiffAndForce(elem::SmallStrainContinuum, state::Array{Float64}, Dstate::Array{Float64}, Δt::Float64)
    ndofs = dofCount(elem);
    nnodes = length(elem.elnodes)
    fint = zeros(Float64, ndofs)
    stiff = zeros(Float64, ndofs,ndofs)
    out = Array{Float64}[]
    u = state[1:nnodes]; v = state[nnodes+1:2*nnodes]
    Du = Dstate[1:nnodes]; Dv = Dstate[nnodes+1:2*nnodes]
    # #@show "u ", u, " Du ", Du
    # #@show "v " v, " Dv " Dv

    for k = 1:length(elem.weights)
        # #@show "Gaussian point ", k
        g1 = elem.dhdx[k][:,1]; g2 = elem.dhdx[k][:,2]

        ux = u'*g1; uy = u'*g2; vx = v'*g1; vy = v'*g2
        Dux = Du'*g1; Duy = Du'*g2; Dvx = Dv'*g1; Dvy = Dv'*g2
        #@show "gauss", k , u, Du, ux, Dux

        # compute  ∂E∂u.T, 8 by 3 array
        ∂E∂u = [g1   zeros(nnodes)    g2;
                zeros(nnodes)    g2   g1;]

        E = [ux; vy; uy+vx]
        DE = [Dux; Dvy; Duy+Dvx]

        # #@show E, DE
        S, dS_dE = getStress(elem.mat[k], E, DE, Δt)

        # @info "gauss ", k, " E ", E, " S ", S

        elem.stress[k] = S
        # @show size(S), size(∂E∂u)
        fint += ∂E∂u * S * elem.weights[k] # 1x8

        stiff += (∂E∂u * dS_dE * ∂E∂u')*elem.weights[k] # 8x8

    end

    return fint, stiff
end

"""
    getStiffAndForce(elem::SmallStrainContinuum, state::Array{Float64}, Dstate::Array{Float64}, Δt::Float64)

Returns the internal force term and the stiffness matrix. `state` and `Dstate` are restriction of full state variables to this element.
"""
function getStiffAndForce1(elem::SmallStrainContinuum, state::Array{Float64}, Dstate::Array{Float64}, Δt::Float64)
    ndofs = length(elem.elnodes);
    nnodes = length(elem.elnodes)
    fint = zeros(Float64, ndofs)
    stiff = zeros(Float64, ndofs,ndofs)
    out = Array{Float64}[]
    u = state[1:nnodes]
    Du = Dstate[1:nnodes]
    # #@show "u ", u, " Du ", Du
    # #@show "v " v, " Dv " Dv

    for k = 1:length(elem.weights)
        # #@show "Gaussian point ", k
        g1 = elem.dhdx[k][:,1]; g2 = elem.dhdx[k][:,2]

        ux = u'*g1; uy = u'*g2
        Dux = Du'*g1; Duy = Du'*g2
        #@show "gauss", k , u, Du, ux, Dux

        # compute  ∂E∂u.T, 4 x 2
        ∂E∂u = [g1   g2]

        E = [ux; uy]
        DE = [Dux; Duy]

        # #@show E, DE
        S, dS_dE = getStress(elem.mat[k], E, DE, Δt) # S: 2 , dS_dE: 2 x 2

        # @info "gauss ", k, " E ", E, " S ", S

        elem.stress[k] = S
        # @show size(S), size(∂E∂u)
        fint += ∂E∂u * S * elem.weights[k] # 4

        stiff += (∂E∂u * dS_dE * ∂E∂u')*elem.weights[k] # 4 x 4

    end

    return fint, stiff
end

"""
    getInternalForce(elem::SmallStrainContinuum, state::Array{Float64}, Dstate::Array{Float64}, Δt::Float64)

Returns the internal force term. `state` and `Dstate` are restriction of full state variables to this element.
"""
function getInternalForce(elem::SmallStrainContinuum, state::Array{Float64}, Dstate::Array{Float64}, Δt::Float64)
    n = dofCount(elem)
    nnodes = length(elem.elnodes)
    fint = zeros(Float64,n)
    out = Array{Float64}[]
    u = state[1:nnodes]; v = state[nnodes+1:end]
    Du = Dstate[1:nnodes]; Dv = Dstate[nnodes+1:end]
    for k = 1:length(elem.weights)
        g1 = elem.dhdx[k][:,1]; g2 = elem.dhdx[k][:,2]

        ux = u'*g1; uy = u'*g2; vx = v'*g1; vy = v'*g2
        Dux = Du'*g1; Duy = Du'*g2; Dvx = Dv'*g1; Dvy = Dv'*g2
        # compute  ∂E∂u.T, 8 by 3 array
        ∂E∂u = [g1   zeros(nnodes)    g2;
                zeros(nnodes)    g2   g1;]

        E = [ux; vy; uy+vx]
        DE = [Dux; Dvy; Duy+Dvx]

        S, dS_dE = getStress(elem.mat[k], E, DE, Δt)

        elem.stress[k] = S

        fint += ∂E∂u * S * elem.weights[k] # 1x8
    end
    return fint
end

"""
    getStrain(elem::SmallStrainContinuum, state::Array{Float64})

Returns the strain of this element.  `state` is restricted to this variable.
"""
function getStrain(elem::SmallStrainContinuum, state::Array{Float64})
    n = dofCount(elem)
    nnodes = length(elem.elnodes)
    u = state[1:nnodes]; v = state[nnodes+1:end]
    nGauss = length(elem.weights)
    E = zeros(nGauss, 3)
    w∂E∂u = zeros(nGauss, 2nnodes, 3)
    for k = 1:nGauss
        g1 = elem.dhdx[k][:,1]; g2 = elem.dhdx[k][:,2]
        ux = u'*g1; uy = u'*g2; vx = v'*g1; vy = v'*g2
        # compute  ∂E∂u.T, 8 by 3 array
        w∂E∂u[k,:,:] = [g1   zeros(nnodes)    g2;
                zeros(nnodes)    g2   g1;] * elem.weights[k]
        E[k,:] = [ux; vy; uy+vx]
    end
    return E, w∂E∂u
end

# for the adjoint solver

# Return:
#    strain{ngps_per_elem, nstrain}
#    dstrain_dstate_tran{neqs_per_elem, ngps_per_elem*nstrain}
function getStrainState(elem::SmallStrainContinuum, state::Array{Float64})
    n = dofCount(elem)
    nnodes = length(elem.elnodes)
    u = state[1:nnodes]; v = state[nnodes+1:end]
    nGauss = length(elem.weights)
    nStrain = 3
    E = zeros(nGauss, nStrain)
    ∂E∂u = zeros(nGauss * nStrain, 2nnodes)

    for k = 1:nGauss
        g1 = elem.dhdx[k][:,1]; g2 = elem.dhdx[k][:,2]
        ux = u'*g1; uy = u'*g2; vx = v'*g1; vy = v'*g2
        E[k,:] = [ux; vy; uy+vx]

        # compute  ∂E∂u
        ∂E∂u[(k-1)*nStrain+1:k*nStrain, :] =
                [g1   zeros(nnodes)    g2;
                zeros(nnodes)    g2   g1;]'

    end
    return E, ∂E∂u
end

function getStiffAndForce(elem::SmallStrainContinuum, state::Array{Float64},
                          stress::Array{Float64,2}, dstress_dstrain_T::Array{Float64,3})
    ndofs = dofCount(elem);
    nnodes = length(elem.elnodes)
    fint = zeros(Float64, ndofs)
    stiff = zeros(Float64, ndofs,ndofs)

    u = state[1:nnodes]; v = state[nnodes+1:2*nnodes]

    for k = 1:length(elem.weights)
        g1 = elem.dhdx[k][:,1]; g2 = elem.dhdx[k][:,2]

        # compute  ∂E∂u.T, 8 by 3 array
        ∂E∂u = [g1   zeros(nnodes)    g2;
                zeros(nnodes)    g2   g1;]

        S, dS_dE_T = stress[k, :], dstress_dstrain_T[k,:,:]

        elem.stress[k] = S
        # @show size(S), size(∂E∂u)
        fint += ∂E∂u * S * elem.weights[k] # 1x8

        stiff += (∂E∂u * dS_dE_T' * ∂E∂u')*elem.weights[k] # 8x8
    end
    return fint, stiff
end

function  getStiffAndDforceDstress(elem::SmallStrainContinuum, state::Array{Float64},
    stress::Array{Float64,2}, dstress_dstrain_T::Array{Float64,3})
    ndofs = dofCount(elem);
    nnodes = length(elem.elnodes)
    nStrain = 3
    nGauss = length(elem.weights)
    stiff = zeros(Float64, ndofs,ndofs)
    dfint_dstress = zeros(Float64,  ndofs, nGauss * nStrain)

    u = state[1:nnodes]; v = state[nnodes+1:2*nnodes]

    for k = 1:length(elem.weights)
        g1 = elem.dhdx[k][:,1]; g2 = elem.dhdx[k][:,2]

        # compute  ∂E∂u, 3 by 2nnodes array
        ∂E∂u = [g1   zeros(nnodes)    g2;
                zeros(nnodes)    g2   g1;]

        # #@show E, DE
        S, dS_dE_T = stress[k, :], dstress_dstrain_T[k,:,:]

        elem.stress[k] = S

        #fint += ∂E∂u * S * elem.weights[k] # 1x8

        dfint_dstress[:, (k-1)*nStrain+1:k*nStrain] = ∂E∂u * elem.weights[k]

        stiff += (∂E∂u * dS_dE_T' * ∂E∂u')*elem.weights[k] # 8x8
    end
    return stiff , dfint_dstress
end

function  getForceAndDforceDstress(elem::SmallStrainContinuum, state::Array{Float64},
    stress::Array{Float64,2})
    ndofs = dofCount(elem);
    nnodes = length(elem.elnodes)
    nStrain = 3
    nGauss = length(elem.weights)
    fint = zeros(Float64, ndofs)
    dfint_dstress = zeros(Float64,  ndofs, nGauss * nStrain)

    u = state[1:nnodes]; v = state[nnodes+1:2*nnodes]

    for k = 1:length(elem.weights)
        g1 = elem.dhdx[k][:,1]; g2 = elem.dhdx[k][:,2]

        # compute  ∂E∂u, 3 by 2nnodes array
        ∂E∂u = [g1   zeros(nnodes)    g2;
                zeros(nnodes)    g2   g1;]

        # #@show E, DE
        S = stress[k, :]

        elem.stress[k] = S

        fint += ∂E∂u * S * elem.weights[k] # 1x8

        dfint_dstress[:, (k-1)*nStrain+1:k*nStrain] = ∂E∂u * elem.weights[k]

    end
    return fint , dfint_dstress
end


text:
## Get unary and binary functions

function is_unary(f::Function, t::Type = Number)
    f ∈ [+, -, *, /, ^] && return false
    for m in methods(f, (t, ))
        m.nargs - 1 > 1 && return false
    end
    return true
end

function is_binary(f::Function, t::Type = Number)
    f ∈ [+, -, *, /, ^] && return true
    !is_unary(f, t)
end

function ariety(f::Function, t::Type = Number)
    is_unary(f, t) && return 1
    is_binary(f, t) && return 2
    return 0
end

function sort_ops(f::Vector{Function})
    U = Function[]
    B = Function[]
    for fi in f
        is_unary(fi) ? push!(U, fi) : push!(B, fi)
    end
    return U, B
end

## Create linear independent basis

count_operation(x::Number, op::Function, nested::Bool = true) = 0
count_operation(x::Sym, op::Function, nested::Bool = true) = 0
count_operation(x::Num, op::Function, nested::Bool = true) = count_operation(value(x), op, nested)

function count_operation(x, op::Function, nested::Bool = true)
    if operation(x)== op
        if is_unary(op)
            # Handles sin, cos and stuff
            nested && return 1 + count_operation(arguments(x), op)
            return 1
        else
            # Handles +, *
            nested && length(arguments(x))-1 + count_operation(arguments(x), op)
            return length(arguments(x))-1
        end
    elseif nested
        return count_operation(arguments(x), op, nested)
    end
    return 0
end

function count_operation(x, ops::AbstractArray, nested::Bool = true)
    return sum([count_operation(x, op, nested) for op in ops])
end

function count_operation(x::AbstractArray, op::Function, nested::Bool = true)
    sum([count_operation(xi, op, nested) for xi in x])
end

function count_operation(x::AbstractArray, ops::AbstractArray, nested::Bool = true)
    counter = 0
    @inbounds for xi in x, op in ops
        counter += count_operation(xi, op, nested)
    end
    counter
end

function split_term!(x::AbstractArray, o, ops::AbstractArray = [+])
    if istree(o)
        n_ops = count_operation(o, ops, false)
        c_ops = 0
        @views begin
            if n_ops == 0
                x[begin]= o
            else
                counter_ = 1
                for oi in arguments(o)
                    c_ops = count_operation(oi, ops, false)
                    split_term!(x[counter_:counter_+c_ops], oi, ops)
                    counter_ += c_ops + 1
                end
            end
        end
    else
        x[begin] = o
    end
    return
end

split_term!(x::AbstractArray,o::Num, ops::AbstractArray = [+]) = split_term!(x, value(o), ops)

remove_constant_factor(x::Num) = remove_constant_factor(value(x))
remove_constant_factor(x::Number) = one(x)

function remove_constant_factor(x)
    # Return, if the function is nested
    istree(x) || return x
    # Count the number of operations
    n_ops = count_operation(x, [*], false)+1
    # Create a new array
    ops = Array{Any}(undef, n_ops)
    @views split_term!(ops, x, [*])
    filter!(x->!isa(x, Number), ops)
    return Num(prod(ops))
end

function remove_constant_factor(o::AbstractArray)
    oi = Array{Any}(undef, size(o))
    for i in eachindex(o)
        oi[i] = remove_constant_factor(o[i])
    end
    return Num.(oi)
end

function create_linear_independent_eqs(ops::AbstractVector, simplify_eqs::Bool = false)
    o = simplify.(ops)
    o = remove_constant_factor(o)
    n_ops = [count_operation(oi, +, false) for oi in o]
    n_x = sum(n_ops) + length(o)
    u_o = Array{Any}(undef, n_x)
    ind_lo, ind_up = 0, 0
    for i in eachindex(o)
        ind_lo = i > 1 ? sum(n_ops[1:i-1]) + i : 1
        ind_up = sum(n_ops[1:i]) + i

        @views split_term!(u_o[ind_lo:ind_up], o[i], [+])
    end
    u_o = remove_constant_factor(u_o)
    unique!(u_o)
    return simplify_eqs ? simplify.(Num.(u_o)) : Num.(u_o)
end


text:
using PhysConsts
#using Base.Test

if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@testset "Data File" begin
    datafile = DataFile("test.url", "directory", "filename")

    @test url(datafile) == "test.url"
    @test dirname(datafile) == abspath("directory")
    @test filename(datafile) == "filename"
    @test filepath(datafile) == joinpath(dirname(datafile), "filename")
end

import PhysConsts: e

@testset "Base" begin
    @test e.quantity == "atomic unit of charge"
    @test e.value == 1.602176634e-19
end


text:
function unsafe_string_maybe(ptr::Cstring)
    if ptr == C_NULL
        return ""
    else
        return unsafe_string(ptr)
    end
end

function NsError_maybe(ptr::MtlError)
    if ptr === C_NULL
        return nothing
    else
        return MtlError(ptr)
    end
end

export @enum_without_prefix

## redeclare enum values without a prefix

# this is useful when enum values from an underlying C library, typically prefixed for the
# lack of namespacing in C, are to be used in Julia where we do have module namespacing.
macro enum_without_prefix(enum, prefix)
    if isa(enum, Symbol)
        mod = __module__
    elseif Meta.isexpr(enum, :(.))
        mod = getfield(__module__, enum.args[1])
        enum = enum.args[2].value
    else
        error("Do not know how to refer to $enum")
    end
    enum = getfield(mod, enum)
    prefix = String(prefix)

    ex = quote end
    for instance in instances(enum)
        name = String(Symbol(instance))
        @assert startswith(name, prefix)
        push!(ex.args, :(const $(Symbol(name[length(prefix) + 1:end])) = $(mod).$(Symbol(name))))
    end

    return esc(ex)
end

##
Base.convert(::Type{MtResourceOptions}, val::UInt32) =
    MtResourceOptions(val)

##
"""
    @mtlthrows error_var function(..., error_var)

Marks that this Metal function has an argument error_var which
must be passed by reference in the underlying ccall, and when
the function returns checks that no error has been set.

Expands roughly to
```julia
error_var = Ref{MTLError}()
result = function(..., error_var)
error[] != C_NULL && throw(MtlError(error[]))
```
"""
macro mtlthrows(error, fun)
    expr = quote
        $error = Ref{MTLError}()
        result = $fun
        if $error[] != C_NULL
            throw(MtlError($(error)[]))
        end
        result
    end
    return esc(expr)
end


text:
struct RewardvsNSim
    mdp
    bandit_t
    n_iters::Int
    seed::Int
    max_steps::Int
    metadata::Dict
end
function RewardvsNSim(mdp, bandit_t, n_iters, seed, max_steps)
    RewardvsNSim(mdp, bandit_t, n_iters, seed, max_steps, Dict(:sim_algorithm=>string(bandit_t), :sim_mdp=>string(mdp), :sim_n_ters=>n_iters, :sim_seed=>seed, :sim_max_steps=>max_steps))
end
function POMDPs.simulate(sim::RewardvsNSim)
    rng = MersenneTwister(sim.seed)
    mdp = sim.mdp
    b = sim.bandit_t(mdp)
    solver = ModularSolver(mdp, b, sim.seed; n_iterations=sim.n_iters)
    policy = solve(solver, mdp)
    s0 = initial_state(mdp, rng)
    hr = HistoryRecorder(; max_steps=sim.max_steps, rng=rng)
    h = simulate(hr, mdp, policy, s0)
    r_total = sum(h.reward_hist)
    r_total
end

@with_kw mutable struct RewardvsN
    mdp = Walk1D()
    bandits=[DPWBandit,RandomBandit,CBTSBandit]
    n_seeds=5
    n_iters=[10,25,50]
    max_steps=20
end
struct RewardvsNResult
    title::String
    data::DataFrame
end
function RewardvsN(::Type{Walk1D})
    mdp = Walk1D()
    RewardvsN(mdp, [DPWBandit,RandomBandit,CBTSDPWBandit], 200, [10,20,50,100,200,400,700,1000,1250], mdp.p.t_max)
end
RewardvsN(::Type{VDPTagMDP}) = RewardvsN(VDPTagMDP(), [DPWBandit,RandomBandit,CBTSDPWBandit], 500, [10,50,100,200,500,1000,1500,2000,2500,3000], 10)
RewardvsN(::Type{InvertedPendulum}) = RewardvsN(InvertedPendulum(), [DPWBandit,RandomBandit,CBTSDPWBandit], 20, [10,50,100,200,500,1000], 100)
RewardvsN(::Type{MountainCar}) = RewardvsN(MountainCar(), [DPWBandit,RandomBandit,CBTSDPWBandit], 20, [100,500,1000,2000,5000,10000], 20)
RewardvsN(::Type{LightDark2D}) = RewardvsN(LightDark2D(), [DPWBandit,RandomBandit,CBTSDPWBandit], 200, [50,100,200,500,1000,1500,2000], 10)

function generate_sim_q(study::RewardvsN)
    q = []
    for b in study.bandits, n in study.n_iters, i in 1:study.n_seeds
        push!(q, RewardvsNSim(study.mdp, b, n, i, study.max_steps))
    end
    q
end

function run_study(study::RewardvsN)
    t1 = time()
    df = DataFrame([String, Int, Int, Float64], [:alg, :seed, :n_iters, :r_total], 0)
    q = generate_sim_q(study)
    result = pmap(POMDPs.simulate, q)
    for (r, sim) in zip(result, q)
        push!(df, [string(sim.bandit_t), sim.seed, sim.n_iters, r])
    end
    println("Elapsed time: $(time()-t1)")
    RewardvsNResult(string(study.mdp), df)
end
@recipe function plot(result::RewardvsNResult, more_results::RewardvsNResult...)
    title := result.title
    xlabel := "number of iterations"
    ylabel := "reward"
    data = vcat(result.data, [r.data for r in more_results]...)
    df = aggregate(data, [:alg, :n_iters], [mean, std])[[:alg, :n_iters, :r_total_mean, :r_total_std]]
    for dd in groupby(df, :alg)
        @series begin
            dd = sort(dd[:], :n_iters)
            label := split(dd[1,:alg], ".")[end]
            legend := :bottomright
            err := dd[:r_total_std] ./ sqrt(maximum(result.data[:seed]))
            dd[:n_iters], dd[:r_total_mean]
        end
    end
end


text:
using YaoTensorNetwork, Yao

nx = 4
ny = 4
depth = 8
c = rand_supremacy2d(nx, ny, depth);
reg = rand_state(nx*ny)
pl = reg |> c |> probs

using Plots
N = 2^(nx*ny)
y = exp.(-N.*pl)
x = N.*pl
order = sortperm(x)
x = x[order]
y = y[order]
plot(x, y, yscale=:log10)


text:
# Static Car Escape
T = Float64;

# model
model = Dynamics.DubinsCar()
n,m = size(model)
x0 = @SVector [2.5,2.5,0.]
xf = @SVector [7.5,2.5,0.]
N = 101
tf = 3.0

# cost
Q = (1e-3)*Diagonal(@SVector ones(n))
R = (1e-2)*Diagonal(@SVector ones(m))
Qf = 100.0*Diagonal(@SVector ones(n))
obj = LQRObjective(Q,R,Qf,xf,N)

# constraints
r = 0.5
s1 = 30; s2 = 50; s3 = 15

circles_escape = NTuple{3,Float64}[]

for i in range(0,stop=5,length=s1)
    push!(circles_escape,(0.,i,r))
end
for i in range(0,stop=5,length=s1)
    push!(circles_escape,(5.,i,r))
end
for i in range(0,stop=5,length=s1)
    push!(circles_escape,(10.,i,r))
end
for i in range(0,stop=10,length=s2)
    push!(circles_escape,(i,0.,r))
end
for i in range(0,stop=3,length=s3)
    push!(circles_escape,(i,5.,r))
end
for i in range(5,stop=8,length=s3)
    push!(circles_escape,(i,5.,r))
end

n_circles_escape = 3*s1 + s2 + 2*s3

circles_escape
x,y,r = collect(zip(circles_escape...))
x = SVector{n_circles_escape}(x)
y = SVector{n_circles_escape}(y)
r = SVector{n_circles_escape}(r)

obs = CircleConstraint(n,x,y,r)
con_obs = ConstraintVals(obs, 2:N-1)

bnd = BoundConstraint(n,m,u_min=-5.,u_max=5.)
con_bnd = ConstraintVals(bnd, 1:N-1)

goal = GoalConstraint(xf)
con_xf = ConstraintVals(goal, N:N)

conSet = ConstraintSet(n,m,[con_obs, con_bnd, con_xf], N)

# Build problem
U0 = [@SVector ones(m) for k = 1:N-1]

car_escape_static = Problem(model, obj, xf, tf;
    constraints=conSet, x0=x0)
initial_controls!(car_escape_static, U0);

X_guess = [2.5 2.5 0.;
           4. 5. .785;
           5. 6.25 0.;
           7.5 6.25 -.261;
           9 5. -1.57;
           7.5 2.5 0.]
X0_escape = interp_rows(N,tf,Array(X_guess'))
initial_states!(car_escape_static, X0_escape)


text:
# %% Interactive bifurcations for 1D energy balance
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("colorscheme.jl"))
using DynamicalSystems, OrdinaryDiffEq, Random
using InteractiveDynamics: record_interaction
import GLMakie

D = 25
omega_vec = randn(D)

ds = Systems.kuramoto(D)
diffeq = (alg = Tsit5(), adaptive = false, dt = 0.05)
integ = integrator(ds; diffeq...)

fig = GLMakie.Figure()
display(fig)
cmap = :curl
axku = GLMakie.Axis(fig[1,1])

# Other static elements
GLMakie.hidedecorations!(axku)
axku.aspect = GLMakie.DataAspect()
GLMakie.lines!(axku, cos.(0:0.001:2π+0.002), sin.(0:0.001:2π+0.002);color = :black)

# Plot balls
phases = GLMakie.Observable(copy(integ.u))
balls = GLMakie.lift(u -> [GLMakie.Point2f0(cos(φ), sin(φ)) for φ in u], phases)
GLMakie.scatter!(axku, balls;
    markersize  = 12, color = 1:D, colormap = cmap, strokewidth=2, strokecolor = :black
)

# Plot arrows
arrows_end = GLMakie.lift(b -> 0.8b, balls)
arrows_start = [GLMakie.Point2f0(0,0) for i in 1:D]
GLMakie.arrows!(axku, arrows_start, arrows_end;
    color = 1:D, colormap = (cmap, 0.5), arrowsize = 10, linewidth = 4,
)

# Plot mean field Vector
# add sliders and buttons
fig[2, 1] = controllayout = GLMakie.GridLayout(tellwidth = false)
run = controllayout[1, 1] = GLMakie.Button(fig; label = "run")
Ks = 0.0:0.1:10.0
kslider = GLMakie.labelslider!(fig, "K =", Ks; sliderkw = Dict(:startvalue => ds.p.K))
controllayout[1, 2] = kslider.layout

# run button functionality
isrunning = GLMakie.Observable(false)
GLMakie.on(run.clicks) do clicks; isrunning[] = !isrunning[]; end
GLMakie.on(run.clicks) do clicks
    @async while isrunning[]
        GLMakie.isopen(fig.scene) || break # ensures computations stop if closed window
        step!(integ)
        phases[] = integ.u
        yield()
        sleep(0.001) # or `yield()` instead
    end
end

GLMakie.on(kslider.slider.value) do val
    integ.p.K = val
    u_modified!(integ, true)
end

# %%
record_interaction(fig, string(@__FILE__)[1:end-2]*"mp4"; total_time = 15, sleep_time = 2)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
function standard_normal_gausshermite(n::Int)
    ϵᵢ, wᵢ = gausshermite(n) # approximates exp(-x²)
    ϵᵢ   .*= sqrt(2.)        # Normalize ϵᵢ and wᵢ nodes to approximate standard normal
    wᵢ   ./= sqrt(π)

    return ϵᵢ, wᵢ
end

"""
```
gausshermite_expectation(f, μ, σ, n = 10)
gausshermite_expectation(f, μ, Σ, n = 10)
gausshermite_expectation(f, μ, Σ, ns)
```

calculates the expectation of a function of a Gaussian random variable/vector.
The first method evalulates ``\\mathbb{E}[f(X)]`` where ``X \\sim N(\\mu, \\sigma)``,
while the other two methods evaluate ``\\mathbb{E}[f(X)]`` where
``X \\sim \\mathcal{N}(\\mu, \\Sigma)`` and ``\\Sigma`` is diagonal.
The latter two methods differ in that the first assumes the same number of
quadrature points in every dimension while the second does not.

### Inputs
- `f::Function`: some function of a random variable. If `f(x) = x`, then
    `gausshermite_expectation(f, μ, σ)` calculates the mean of ``N(\\mu, \\sigma)``
    using 10-point Gauss-Hermite quadrature.
- `μ::Number` or `μ::AbstractVector`: mean of the Gaussian random variable/vector.
- `σ::Number`: standard deviation of the Gaussian random variable.
- `Σ::AbstractVector`: diagonal of the variance-covariance matrix of
     the Gaussian random vector.
- `n::Int`: number of quadrature points to use
- `ns::AbstractVector{Int}` or `ns::NTuple{N, Int} where N`: number of quadrature points to use
    in each dimension of the Gaussian random vector.
"""
function gausshermite_expectation(f::Function, μ::Number, σ::Number, n::Int = 10)
    ϵᵢ, wᵢ = gausshermite(n)
    ϵᵢ   .*= sqrt(2.)  # Normalize ϵᵢ and wᵢ nodes to approximate standard normal
    # wᵢ   ./= sqrt(π) # This step done later to reduce number of computations

    if μ ≈ 0.
        return sum([wᵢ[i] * f(ϵᵢ[i] * σ) for i in 1:n]) / sqrt(π)
    else
        return sum([wᵢ[i] * f(ϵᵢ[i] * σ + μ) for i in 1:n]) / sqrt(π)
    end

end

function gausshermite_expectation(f::Function, μ::AbstractVector{S},
                                  Σ::AbstractVector{<: Number}, n::Int = 10) where {S <: Number}

    d = length(μ)
    @assert length(Σ) == d "The length of μ and Σ must be the same."
    ϵ, w = gausshermite(n)
    ϵ   .*= sqrt(2.)  # Normalize ϵ and w nodes to approximate standard normal
    # w   ./= sqrt(π) # This step done later to reduce number of computations

    # Evaluate over the tensor grid
    feval = Array{S}(undef, (n for i in 1:d)...)
    allCI = CartesianIndices(feval)
    if all(μ .≈ 0.)
        @simd for CI in allCI
            feval[CI] = f([ϵ[i] for i in Tuple(CI)] .* Σ)
        end
    else
        @simd for CI in allCI
            feval[CI] = f([ϵ[i] for i in Tuple(CI)] .* Σ + μ)
        end
    end

    for n_dim in 1:(d - 1)
        # Iteratively integrate out each dimension, i.e. law of iterated expectations
        iter = CartesianIndices(tuple(Tuple(1:n for i in 1:(d - n_dim))...,
                                      Tuple(1:1 for i in 1:n_dim)...)) # Create CartesianIndices for all remaining dimensions
        # ((1:n for i in 1:(d - n_dim + 1))..., (1 for i in 1:(n_dim - 1))...) creates a Tuple of 1:n for the dimensions
        # that are not to be integrated out and uses 1s for the remaining dimensions. We want to use each dimension of feval
        # from 1 to (d - n_dim) (inclusive). So on the first iteration, the tuple should be (1:n, 1:n).
        # We then assign it to the dimensions of feval from 1 to (d - n_dim - 1) (inclusive) to avoid allocations
        feval[iter] .= dropdims(sum(mapslices(fᵢ -> fᵢ .* w, (@view feval[((1:n for i in 1:(d - n_dim + 1))...,
                                                                           (1 for i in 1:(n_dim - 1))...)...]),
                                              dims = (d - n_dim) + 1), dims = (d - n_dim) + 1), dims = (d - n_dim) + 1)
    end

    # Handle final integration on its own
    return sum(w .* (@view feval[:, (1 for i in 1:(d - 1))...])) / π^(d / 2)
end

function gausshermite_expectation(f::Function, μ::AbstractVector{S},
                                  Σ::AbstractVector{<: Number}, ns::AbstractVector{Int}) where {S <: Number}

    d = length(μ)
    @assert length(Σ) == d "The length of μ and Σ must be the same."
    ϵ = Dict{Int, Vector{S}}()
    w = Dict{Int, Vector{S}}()
    for i in 1:d
        ϵ[i], w[i] = gausshermite(ns[i])
        ϵ[i]     .*= sqrt(2.)  # Normalize ϵ and w nodes to approximate standard normal
        # w[i]   ./= sqrt(π) # This step done later to reduce number of computations
    end

    # Evaluate over the tensor grid
    feval = Array{S}(undef, (n for n in ns)...)
    allCI = CartesianIndices(feval)
    if all(μ .≈ 0.)
        @simd for CI in allCI
            feval[CI] = f([ϵ[n_dim][gridᵢ] for (n_dim, gridᵢ) in enumerate(Tuple(CI))] .* Σ)
        end
    else
        @simd for CI in allCI
            feval[CI] = f([ϵ[n_dim][gridᵢ] for (n_dim, gridᵢ) in enumerate(Tuple(CI))] .* Σ + μ)
        end
    end

    # Iteratively integrate out each dimension, i.e. law of iterated expectations
    for n_dim in 1:(d - 1)
        iter = CartesianIndices(tuple(Tuple(1:ns[i] for i in 1:(d - n_dim))...,
                                      Tuple(1:1 for i in 1:n_dim)...))
        feval[iter, 1] .= dropdims(sum(mapslices(fᵢ -> fᵢ .* w[d - n_dim + 1], (@view feval[((1:ns[i] for i in 1:(d - n_dim + 1))...,
                                                                                (1 for i in 1:(n_dim - 1))...)...]),
                                                 dims = (d - n_dim) + 1), dims = (d - n_dim) + 1), dims = (d - n_dim) + 1)
    end

    # Handle final integration on its own
    return sum(w[1] .* (@view feval[:, (1 for i in 1:(d - 1))...])) / π^(d / 2)
end

function gausshermite_expectation(f::Function, μ::AbstractVector{S},
                                  Σ::AbstractVector{<: Number}, ns::NTuple{N, Int}) where {S<: Number, N}

    d = length(μ)
    @assert length(Σ) == d "The length of μ and Σ must be the same."
    ϵ = Dict{Int, Vector{S}}()
    w = Dict{Int, Vector{S}}()
    for i in 1:d
        ϵ[i], w[i] = gausshermite(ns[i])
        ϵ[i]     .*= sqrt(2.)  # Normalize ϵ and w nodes to approximate standard normal
        # w[i]   ./= sqrt(π) # This step done later to reduce number of computations
    end

    # Evaluate over the tensor grid
    feval = Array{S}(undef, (n for n in ns)...)
    allCI = CartesianIndices(feval)
    if all(μ .≈ 0.)
        @simd for CI in allCI
            feval[CI] = f([ϵ[n_dim][gridᵢ] for (n_dim, gridᵢ) in enumerate(Tuple(CI))] .* Σ)
        end
    else
        @simd for CI in allCI
            feval[CI] = f([ϵ[n_dim][gridᵢ] for (n_dim, gridᵢ) in enumerate(Tuple(CI))] .* Σ + μ)
        end
    end

    # Iteratively integrate out each dimension, i.e. law of iterated expectations
    for n_dim in 1:(d - 1)
        iter = CartesianIndices(tuple(Tuple(1:ns[i] for i in 1:(d - n_dim))...,
                                      Tuple(1:1 for i in 1:n_dim)...))
        feval[iter, 1] .= dropdims(sum(mapslices(fᵢ -> fᵢ .* w[d - n_dim + 1], (@view feval[((1:ns[i] for i in 1:(d - n_dim + 1))...,
                                                                                (1 for i in 1:(n_dim - 1))...)...]),
                                                 dims = (d - n_dim) + 1), dims = (d - n_dim) + 1), dims = (d - n_dim) + 1)
    end

    # Handle final integration on its own
    return sum(w[1] .* (@view feval[:, (1 for i in 1:(d - 1))...])) / π^(d / 2)
end


text:
const QUAD_DATA = [
[(0.0,),
(2.0,)],
[(-0.5773502691896258, 0.5773502691896258,),
(1.0, 1.0,)],
[(-0.7745966692414834, 0.0, 0.7745966692414834,),
(0.5555555555555556, 0.8888888888888888, 0.5555555555555556,)],
[(-0.8611363115940526, -0.3399810435848563, 0.3399810435848563, 0.8611363115940526,),
(0.34785484513745385, 0.6521451548625462, 0.6521451548625462, 0.34785484513745385,)],
[(-0.906179845938664, -0.5384693101056831, 0.0, 0.5384693101056831, 0.906179845938664,),
(0.23692688505618908, 0.47862867049936647, 0.5688888888888889, 0.47862867049936647, 0.23692688505618908,)],
]


text:
function qti(qs, id)
    xdoc = XMLDocument()
    root = create_root(xdoc, "questestinterop")
    set_attributes(
        root,
        Dict(
            "xmlns"=>"http://www.imsglobal.org/xsd/ims_qtiasiv1p2",
            "xmlns:xsi"=>"http://www.w3.org/2001/XMLSchema-instance",
            "xsi:schemaLocation"=>"http://www.imsglobal.org/xsd/ims_qtiasiv1p2 http://www.imsglobal.org/xsd/ims_qtiasiv1p2p1.xsd"
        )
    )

    assessment = new_child(root, "assessment")
    set_attribute(assessment, "id", id)

    section = new_child(assessment, "section")
    set_attribute(section, "ident", "root_section")

    for q in qs
        item = new_child(section, "item")
        set_attributes(
            item,
            Dict(
                "ident" => q.id
            )
        )

        ## metadata
        itemmetadata = new_child(item, "itemmetadata")
        qtimetadata = new_child(itemmetadata, "qtimetadata")
        fields = [
            "question_type" => "multiple_choice_question"
        ]
        for field in fields
            xfield = new_child(qtimetadata, "qtimetadatafield")
            lab = new_child(xfield, "fieldlabel")
            add_text(lab, first(field))
            entry = new_child(xfield, "fieldentry")
            add_text(entry, last(field))
        end

        ## presentation
        pres = new_child(item, "presentation")
        mat = new_child(pres, "material")
        mattext = new_child(mat, "mattext")
        set_attributes(
            mattext,
            Dict(
                "texttype" => "text/html"
            )
        )
        for c in q.content
            add_text(mattext, string(c))
        end

        respid = "response1"
        response_lid = new_child(pres, "response_lid")
        set_attributes(
            response_lid,
            Dict(
                "ident" => respid,
                "rcardinality" => "Single"
            )
        )

        render_choice = new_child(response_lid, "render_choice")

        for r in q.responses
            response_label = new_child(render_choice, "response_label")
            set_attribute(response_label, "ident", r.id)

            mat = new_child(response_label, "material")
            mattext = new_child(mat, "mattext")
            set_attributes(
                mattext,
                Dict(
                    "texttype" => "text/html"
                )
            )
            add_text(mattext, string(r.content))
        end

        # correct response...
        correct = q.responses[something(findfirst(x -> x.correct, q.responses), 1)]

        resproc = new_child(item, "resproessing")
        outcomes = new_child(resproc, "outcomes")
        decvar = new_child(outcomes, "decvar")
        set_attributes(
            decvar,
            Dict(
                "maxvalue" => "100",
                "minvalue" => "0",
                "varname" => "SCORE",
                "vartype" => "Decimal"
            )
        )

        respcond = new_child(resproc, "respcondition")
        set_attribute(respcond, "continue", "No")
        condvar = new_child(respcond, "conditionvar")
        vareq = new_child(condvar, "varequal")
        set_attribute(vareq, "respident", "response1")
        add_text(vareq, correct.id)

        setvar = new_child(respcond, "setvar")
        set_attributes(
            setvar,
            Dict(
                "action" => "Set",
                "varname" => "SCORE"
            )
        )
        add_text(setvar, "100")

    end
    return(xdoc)
end

function package(name, qs::Pair...; media=Dict())

    package_dir = joinpath("imports", string(name))
    if !isdir(package_dir)
        try
            mkpath(package_dir)
        catch
            @warn "couldn't create package dir $package_dir"
        end
    end

    manifest_xdoc = XMLDocument()
    manifest = create_root(manifest_xdoc, "manifest")
    set_attributes(
        manifest,
        Dict(
            "identifier"=>"gb89c8d93f022f423d433a803104bde1a",
            "xmlns" => "http://www.imsglobal.org/xsd/imsccv1p1/imscp_v1p1",
            "xmlns:lom" => "http://ltsc.ieee.org/xsd/imsccv1p1/LOM/resource",
            "xmlns:imsmd" => "http://www.imsglobal.org/xsd/imsmd_v1p2",
            "xmlns:xsi" => "http://www.w3.org/2001/XMLSchema-instance",
            "xsi:schemaLocation" => "http://www.imsglobal.org/xsd/imsccv1p1/imscp_v1p1 http://www.imsglobal.org/xsd/imscp_v1p1.xsd http://ltsc.ieee.org/xsd/imsccv1p1/LOM/resource http://www.imsglobal.org/profile/cc/ccv1p1/LOM/ccv1p1_lomresource_v1p0.xsd http://www.imsglobal.org/xsd/imsmd_v1p2 http://www.imsglobal.org/xsd/imsmd_v1p2p2.xsd"
        )
    )

    println("Packaging resources:")

    new_child(manifest, "organizations")
    resources = new_child(manifest, "resources")
    res = new_child(resources, "resource")
    set_attributes(
        res,
        Dict(
            "identifier" => name,
            "type" => "imsqti_xmlv1p2"
        )
    )

    for (name, q) in qs
        q_xml = qti(q, name)
        filename = "$name.xml"
        save_file(q_xml, joinpath(package_dir, filename))

        println("✔ $(length(q)) questions → $filename")

        file = new_child(res, "file")
        set_attribute(file, "href", filename)
    end

    for (name, src) in media
        println("✔ $name")
        dest = joinpath(package_dir, name)
        isdir(dirname(dest)) || mkpath(dirname(dest))
        write(dest, src)

        file = new_child(res, "file")
        set_attribute(file, "href", name)
    end

    save_file(manifest_xdoc, joinpath(package_dir, "imsmanifest.xml"))
    run(`zip -r imports/$name.qti.zip $package_dir`)

end


text:
module OpCode

# this has to be a fixed size (single type) for efficiency (needs to
# pack into an array).  for simplicity(?) we go with the smallest
# format possible, which seems to be a type, an integer and a
# character

immutable Instruction{C}
    code::Int
    n::Int
    c::C
end

const EQ = 1  # fail if not eq (ignore n)
const GE = 2  # fail if not ge (ignore n)
const LE = 3  # fail if not le (ignore n)
const ALT = 4 # catch failure and jump to n (ignore c)

end


text:
using Distributed
import Nanosoldier, GitHub

nodes = Dict(Any => addprocs(1; exeflags="--project"))
@everywhere import Nanosoldier

auth = GitHub.authenticate(ENV["GITHUB_AUTH"])
secret = ENV["GITHUB_SECRET"]

config = Nanosoldier.Config(ENV["USER"], nodes, auth, secret;
                            workdir = joinpath(dirname(@__DIR__), "workdir"),
                            trackrepo = "vtjnash/julia",
                            reportrepo = "vtjnash/NanosoldierReports",
                            trigger = r"\@nanosoldier\s*`runtests\(.*?\)`",
                            admin = "vtjnash",
                            bucket = "julialang-reports/nanosoldier")

server = Nanosoldier.Server(config)

using Sockets
run(server, IPv4(0), 8888)


text:
# Randomization option slows down processing
# but improves quality of solution considerably
# Would be better to do randomization in place
function cddual(X::SparseOrFullMat,
	            Y::Vector;
	            C::Real = 1.0,
	            norm::Integer = 2,
	            randomized::Bool = true,
	            maxpasses::Integer = 1000,
	            eps::Real = 0.1)
	const PG_EPS = 1.0e-12

	# l: # of samples
	# n: # of features
	n, l = size(X)
	alpha = zeros(l)
	w = zeros(n)

	# Set U and D
	#  * L1-SVM: U = C, D[i] = 0
	#  * L2-SVM: U = Inf, D[i] = 1 / (2C)
	U = 0.0
	D = Array(Float64, l)
	if norm == 1
		U = C
		fill!(D, 0.0)
	elseif norm	== 2
		U = Inf
		fill!(D, 1.0 / (2.0 * C))
	else
		DomainError("Only L1-SVM and L2-SVM are supported")
	end

	# Set Qbar
	Qbar = Array(Float64, l)
	for i in 1:l
		Qbar[i] = D[i] + vecnorm(X[:, i], 2)^2
	end

	# Loop over examples
	converged = false
	pass = 0

	indices = randomized ? randperm(l) : [1:l]

	while !converged
		pass += 1
		if pass > maxpasses
			break
		end
		inv_ins = 0
		pg_max = -Inf
		pg_min = Inf

		if randomized
			shuffle!(indices)
		end

		# Process all observations
		for i in indices
			g = Y[i] * inner_prod(w, X[:, i]) - 1.0 + D[i] * alpha[i]

			if alpha[i] == 0.0
				pg = min(g, 0.0)
				if g <= 0.0
					inv_ins += 1
				end
			elseif alpha[i] == U
				pg = max(g, 0.0)
				if g >= 0.0
					inv_ins += 1
				end
			else
				pg = g
			end
			pg_max = max(pg_max, pg)
			pg_min = min(pg_min, pg)

			if abs(pg) > PG_EPS
				alphabar = alpha[i]
				alpha[i] = min(max(alpha[i] - g / Qbar[i], 0.0), U)
				add_vec!(w, X[:, i], e -> (alpha[i] - alphabar) * Y[i] * e)
			end
		end

		if (pg_max - pg_min <= eps) && (inv_ins == 0)
			converged = true
			break
		end
	end

	return SVMFit(w, pass, converged)
end


text:
import RIOPA

function julia_main()::Cint
    return RIOPA.main(ARGS)
end

if !isdefined(Base, :active_repl)
    julia_main()
end


text:
struct TextPrinter
end

function dprint(dict::AbstractDict; reorder::Bool=false)
    if reorder
        dict = sort(dict, by = x -> x[1])
    end
    for key in keys(dict)
        print(key, " => ", dict[key], "\n")
    end
end

max_stars_repo_name: oilneck/InverseIsing

max_issues_repo_name: oilneck/InverseIsing

max_forks_repo_name: oilneck/InverseIsing


text:
# ---
# title: 717. 1-bit and 2-bit Characters
# id: problem717
# author: zhwang
# date: 2022-02-20
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/1-bit-and-2-bit-characters/description/>
# hidden: true
# ---
#
# We have two special characters. The first character can be represented by one
# bit `0`. The second character can be represented by two bits (`10` or `11`).
#
# Now given a string represented by several bits. Return whether the last
# character must be a one-bit character or not. The given string will always end
# with a zero.
#
# **Example 1:**
#
#
#
#     Input:
#     bits = [1, 0, 0]
#     Output: True
#     Explanation:
#     The only way to decode it is two-bit character and one-bit character. So the last character is one-bit character.
#
#
# **Example 2:**
#
#
#
#     Input:
#     bits = [1, 1, 1, 0]
#     Output: False
#     Explanation:
#     The only way to decode it is two-bit character and two-bit character. So the last character is NOT one-bit character.
#
#
# **Note:**
#
# * `1 <= len(bits) <= 1000`.
# * `bits[i]` is always `0` or `1`.
#
#
## @lc code=start
using LeetCode

function is_one_bit_character(bits::Vector{Int})::Bool
    pos = findlast(x -> x == 0, @view bits[1:(end - 1)])
    pos === nothing && return length(bits) % 2 == 1
    return (length(bits) - pos) % 2 == 1
end
## @lc code=end


text:
import ImageView
using Colors, Gtk.ShortNames
using TestImages

grid, frames, c = ImageView.canvasgrid((2,2))
Gtk.showall(Window("canvasgrid", 800, 600) |> grid)
imshow(c[1,1], testimage("lighthouse"))
imshow(c[1,2], testimage("mountainstream"))
imshow(c[2,1], testimage("moonsurface"))
imshow(c[2,2], testimage("mandrill"))
sleep(0.01)


text:
@testset "_corresponding_type" begin

    for AS in [AbstractContinuousSystem, AbstractDiscreteSystem]
        for S in subtypes(AS)
            @test MathematicalSystems._corresponding_type(AS, fieldnames.(S)) == S
        end
    end
end

# ========
# Setup
# ========

n = 3; m = 2; l = 3
A = A1 = rand(n, n)
B = B1 = rand(n, m)
c = c1 = rand(n)
D = D1 = rand(n, l)
E = E1 = rand(l, l)
X = X1 = BallInf(zeros(n), 50.)
U = U1 = BallInf(zeros(m), 70.)
W = W1 = BallInf(zeros(l), 100.)

f1(x) = x'*x
f1(x, u) = x'*x + u'*u
f1(x, u, w) = x'*x + u'*u + w'*w

# ===================
# Error Handling
# ===================
@testset "@system miscellaneous" begin
    @testset "@system for arbitrary order of set specification" begin
        sys = @system(x' = Ax + Bu, u∈U, x∈X)
        @test sys == ConstrainedLinearControlContinuousSystem(A,B,X,U)
        sys = @system(x' = Ax + Bu + Dw, u∈U, w∈W, x∈X)
        @test sys == NoisyConstrainedLinearControlContinuousSystem(A,B,D,X,U,W)
    end
    @testset "@system for arbitrary order of rhs terms" begin
        sys = @system(x' = Bu + Ax, u∈U, x∈X)
        @test sys == ConstrainedLinearControlContinuousSystem(A,B,X,U)
        sys = @system(x' = Dw + Ax + Bu , u∈U, w∈W, x∈X)
        @test sys == NoisyConstrainedLinearControlContinuousSystem(A,B,D,X,U,W)
        sys = @system(x' = c + Ax + Bu , u∈U, x∈X)
        @test sys == ConstrainedAffineControlContinuousSystem(A,B,c,X,U)
    end

    @testset "@system for incorrect variable handling" begin
        @test_throws ArgumentError("state and input variables have the same name `a`") @system(a⁺= Aa + Bw, input:a)
        @test_throws ArgumentError("state and input variables have the same name `x`") @system(x⁺= Ax + Bx, input:x)
        @test_throws ArgumentError("input and noise variables have the same name `u`") @system(x⁺= Ax + Bu + Bw, input:u, noise:u)
        @test_throws ArgumentError("there is more than one input term `w`") @system(x⁺= Ax + Bw + Dw, input:w)
        @test_throws ArgumentError("there is more than one input term `u`") @system(x⁺= Ax + Bu + Bu, noise:u)
        @test_throws ArgumentError("there is more than one state term `x`") @system(x⁺= Ax + Bx)
        @test_throws ArgumentError("there is more than one constant term")  @system(x⁺= Ax + c + c)
        # @test_throws ArgumentError("input and noise variables have the same name `u`") @system(x⁺= Ax + Bu + Bw,  noise:u)
        # @test_throws ArgumentError("input and noise variables have the same name `w`") @system(x⁺= Ax + Bu + Bw,  input:w)
    end
end

# ===================
# Continuous systems
# ===================
@testset "@system for continous identity systems" begin
    @test @system(x' = 0, dim: 2) == ContinuousIdentitySystem(2)
    sys = @system x' = 0 dim: 2
    @test sys == ContinuousIdentitySystem(2)

    @test @system(x₁' = 0, dim=2) == ContinuousIdentitySystem(2)

    @test @system(x' = 0, dim: 2, x ∈ X) == ConstrainedContinuousIdentitySystem(2, X)
    @test @system(x' = 0, dim=2, x ∈ X1) == ConstrainedContinuousIdentitySystem(2, X)
    sys = @system x' = 0 dim=2 x ∈ X1
    @test sys == ConstrainedContinuousIdentitySystem(2, X)
end

@testset "@system for linear continous systems" begin
    @test @system(x' = A*x) == LinearContinuousSystem(A)
    @test @system(x1' = A1x1) == LinearContinuousSystem(A1)
    sys = @system x1' = A1x1
    @test sys == LinearContinuousSystem(A1)

    # automatic identification of rhs linearity
    @test @system(x' = -x) == LinearContinuousSystem(-1.0*IdentityMultiple(I, 1))
    @test @system(x' = x, dim=3) == LinearContinuousSystem(1.0*IdentityMultiple(I, 3))
    @test @system(x' = 2x, dim=3) == LinearContinuousSystem(2.0*IdentityMultiple(I, 3))

    @test @system(x' = A*x, x ∈ X) == ConstrainedLinearContinuousSystem(A,X)
    @test @system(x1' = A1x1, x1 ∈ X1) == ConstrainedLinearContinuousSystem(A1,X1)
    sys = @system x1' =A1x1     x1 ∈ X1
    @test sys == ConstrainedLinearContinuousSystem(A1,X1)
end

@testset "@system for linear control continuous systems" begin
    # if the state should be named `w`
    @test @system(w' = Aw + Bu) == LinearControlContinuousSystem(A, B)
    # but if the input should be named `w`
    # without specification, a noisy system is returned
    @test  @system(x' = Ax + Bw, x∈X, w∈W) == NoisyConstrainedLinearContinuousSystem(A, B, X, W)
    # but if we use the `input:w`, a controlled system is returned
    @test @system(x' = Ax + Bw, input:w) == LinearControlContinuousSystem(A, B)
    # and in general, if the input name is different from `u`
    @test @system(x' = Ax + Bu_1, input:u_1) == LinearControlContinuousSystem(A, B)
    @test @system(x' = Ax + B*u_1, input:u_1) == LinearControlContinuousSystem(A, B)

    @test @system(x' = A1x + B1u) == LinearControlContinuousSystem(A1, B1)
    @test @system(z_1' = A*z_1 + B*u_1, input:u_1) == LinearControlContinuousSystem(A, B)

    @test @system(x' = Ax + Bu, x ∈ X, u ∈ U) == ConstrainedLinearControlContinuousSystem(A, B, X, U)

    # if variable in front of input is ommitted add identity matrix
    @system(x' = Ax + u) == LinearControlContinuousSystem(A, IdentityMultiple(1.0*I,n))

    # if * are used x_ = A_*x_ + B_*u_, u_ is interpreted as input variable,
    # independent of the name used for u_
    @test @system(w' = A*w + B*u_1) == LinearControlContinuousSystem(A, B)
    # similarily for x_ = A_*x_ + B_*u_ + c_
    @test @system(w' = A*w + B*u_1 + c, w∈X, u_1∈U) == ConstrainedAffineControlContinuousSystem(A, B, c, X, U)

    # scalar cases
    @test @system(x' = 0.5x + u) == LinearControlContinuousSystem(hcat(0.5), I(1.0, 1))
    @test @system(x' = 0.5x + 1.) == AffineContinuousSystem(hcat(0.5), vcat(1.))
    @test @system(x' = x + u) == LinearControlContinuousSystem(I(1.0, 1), I(1.0, 1))
    @test @system(x' = x + 0.1u) == LinearControlContinuousSystem(I(1.0, 1), hcat(0.1))
    @test @system(x' = x + 0.1*u) == LinearControlContinuousSystem(I(1.0, 1), hcat(0.1))
    sys = @system(x' = 0.3x + 0.1u + 0.2, x∈X, u∈U)
    @test sys == ConstrainedAffineControlContinuousSystem(hcat(0.3), hcat(0.1), vcat(0.2), X, U)
end

@testset "@system for linear algebraic continous systems" begin
    # lhs needs a *
    @test_throws ArgumentError @system(Ex' = Ax)
    @test@system(E*x' = Ax) == LinearAlgebraicContinuousSystem(A, E)

    @test @system(E*x' = A*x) == LinearAlgebraicContinuousSystem(A, E)
    @test @system(E*x' = A*x, x∈X) == ConstrainedLinearAlgebraicContinuousSystem(A, E, X)
end

@testset "@system for affine continuous systems" begin
    b = [1.0]
    @test @system(x' = x + b) == AffineContinuousSystem(I(1), b)
    @test @system(x' = A*x + c) == AffineContinuousSystem(A, c)
    sys =  @system(z_1' = A*z_1 + B*v_1 + c1, z_1 ∈ X, v_1 ∈ U1, input:v_1)
    @test sys == ConstrainedAffineControlContinuousSystem(A, B, c1, X, U1)
   sys =  @system(x' = Ax + Bu + c)
    @test sys == AffineControlContinuousSystem(A, B, c)
end

@testset "@system for noisy continous systems" begin
    sys = @system(x' = f1(x, u, w), x ∈ X, u ∈ U, w ∈ W, dims=(1, 2,3))
    @test sys == NoisyConstrainedBlackBoxControlContinuousSystem(f1, 1, 2, 3, X, U, W)

    sys = @system(x' = Ax + Bu + Dw, x ∈ X, u ∈ U1, w ∈ W1)
    @test sys == NoisyConstrainedLinearControlContinuousSystem(A, B, D, X, U1, W1)
    sys = @system(z_1' = Az_1 + B*v_1 + c + Dd_1, z_1 ∈ X, v_1 ∈ U1, d_1 ∈ W1, input:v_1, d_1:noise)
    @test sys == NoisyConstrainedAffineControlContinuousSystem(A, B, c, D, X, U1, W1)

    # if variable in front of input or noise  is omitted, add identity matrix
    sys = @system(x' = Ax + u + c + w, x∈X, u∈U, w∈W)
    sys == NoisyConstrainedAffineControlContinuousSystem(A, IdentityMultiple(1.0*I,n), c, IdentityMultiple(1.0*I,n), X, U, W)
end

@testset "@system for black-box continous systems" begin
    @test_throws ArgumentError @system(x' = f1(x))
    @test_throws ArgumentError @system(x' = f1(x, u))
    sys =  @system(x' = f1(x), dim:3)
    @test sys == BlackBoxContinuousSystem(f1, 3)
    sys =  @system(x' = f1(x), x ∈ X, dim:2)
    @test sys == ConstrainedBlackBoxContinuousSystem(f1, 2, X)
    sys = @system(x' = f1(x, u), x ∈ X, u ∈ U, dims=(1, 2))
    @test sys == ConstrainedBlackBoxControlContinuousSystem(f1, 1, 2, X, U)
    # allow for arbitrary input with definition for rhs of the form f_(x_, u_)
    sys = @system(x' = f1(x, u123), x ∈ X, u123 ∈ U, dims=(2, 3))
    @test sys == ConstrainedBlackBoxControlContinuousSystem(f1, 2, 3, X, U)
end

# ==================
# Discrete systems
# ==================

@testset "@system for discrete identity systems" begin
    @test @system(x⁺ = x, dim=2) == DiscreteIdentitySystem(2)
    @test @system(x₁⁺ = x₁, dim: 2) == DiscreteIdentitySystem(2)

    @test @system(u⁺ = u, dim: 3, u ∈ U) == ConstrainedDiscreteIdentitySystem(3, U)
    @test @system(x1⁺ = x1, dim = 3, x1 ∈ X1) == ConstrainedDiscreteIdentitySystem(3, X1)

    # emoij support 😉
    🚈 = X
    sys = @system(👨⁺ = 👨, dim: 2, 👨∈🚈)
    @test sys == ConstrainedDiscreteIdentitySystem(2, X)
end

@testset "@system for linear discrete systems" begin
    @test @system(x⁺ = A*x) == LinearDiscreteSystem(A)
    @test @system(x⁺ = Ax) == LinearDiscreteSystem(A)
    @test @system(z⁺ = A*z) == LinearDiscreteSystem(A)
    @test @system(z⁺ = Az) == LinearDiscreteSystem(A)
    @test_throws ArgumentError @system(z⁺ = Ax)
    @test_throws ArgumentError @system(x⁺ = Az)
    @test @system(x1⁺ = A*x1) == LinearDiscreteSystem(A)
    @test @system(x1⁺ = Ax1) == LinearDiscreteSystem(A)

    @test @system(x⁺ = 2x, dim=3) == LinearDiscreteSystem(2.0*IdentityMultiple(I,3))

    @test @system(x⁺ = A1x, x ∈ X1) == ConstrainedLinearDiscreteSystem(A1, X1)
end

@testset "@system for affine discrete systems" begin
    @test @system(x⁺ = A*x + c) == AffineDiscreteSystem(A, c)
    @test @system(x⁺ = A1*x + c1) == AffineDiscreteSystem(A1, c1)
end

@testset "@system for linear algebraic discrete systems" begin
    sys = @system(E*x⁺ = Ax)
    @test sys == LinearAlgebraicDiscreteSystem(A, E)
    @test sys == @system(E1*x⁺ = A1*x)

    sys = @system(E*x⁺ = A*x, x ∈ X)
    @test sys == ConstrainedLinearAlgebraicDiscreteSystem(A, E, X)
end

@testset "@system for linear control discrete systems" begin
    sys = @system(x⁺ = A*x + B1*u)
    @test sys == LinearControlDiscreteSystem(A, B)
    @test sys == @system(x⁺ = A1*x + Bu)

    @test @system(x⁺ = A1x + B1u) == LinearControlDiscreteSystem(A1, B1)
    @test @system(z_1⁺ = A*z_1 + B*v_1, input:v_1) == LinearControlDiscreteSystem(A, B)

    # u or other symbols are used for inputs; here u is interpreted as input
    sys = @system(x⁺ = Ax + Bu, x ∈ X, u ∈ U1)
    @test sys == ConstrainedLinearControlDiscreteSystem(A, B, X, U1)

    # here v is interpreted as input (if the input has more than one letter, use a *)
    @system(x⁺ = Ax + B*u123, x ∈ X, u123 ∈ U1, input: u123)
    @test sys == ConstrainedLinearControlDiscreteSystem(A, B, X, U1)

    # if variable in front of input is omitted, add identity matrix
    @system(x⁺ = Ax + u) == LinearControlDiscreteSystem(A, IdentityMultiple(1.0*I,n))

    sys = @system(x⁺ = Ax + Bu123, x ∈ X, u123 ∈ U1, input: u123)
    @test sys == ConstrainedLinearControlDiscreteSystem(A, B, X, U1)
end

@testset "@system for noisy discrete systems" begin
    sys = @system(x⁺ = A1x + Dw, x ∈ X1, w ∈ W)
    @test sys == NoisyConstrainedLinearDiscreteSystem(A1, D, X1, W1)

    # by default w is reserved for noise; here w is interpreted as noise
    sys = @system(x⁺ = A1x + Dw, x ∈ X1, w ∈ W)
    @test sys == NoisyConstrainedLinearDiscreteSystem(A1, D, X1, W)

    sys = @system(x⁺ = A1x + Dw, x ∈ X1, w ∈ W, w:noise)
    @test sys == NoisyConstrainedLinearDiscreteSystem(A1, D, X1, W1)

    sys = @system(x⁺ = Ax + Bu + Dw, x ∈ X, u ∈ U1, w ∈ W1, w:noise, u:input)
    @test sys == NoisyConstrainedLinearControlDiscreteSystem(A, B, D, X, U1, W)

    sys = @system(x⁺ = Ax + Bu + c + Dw, x ∈ X, u ∈ U1, w ∈ W1, w:noise, u:input)
    @test sys == NoisyConstrainedAffineControlDiscreteSystem(A, B, c, D, X, U1, W1)

    # check error handling for repeated variable names
    @test_throws ArgumentError @system(x⁺ = Ax + Bw + Dw, x ∈ X, w ∈ U1, w ∈ W1, noise=w, input=w)
end

@testset "@system for black-box discrete systems" begin
    sys = @system(x⁺ = f1(x), dim=2)
    @test sys == BlackBoxDiscreteSystem(f1, 2)

    sys = @system(x⁺ = f1(x), dim: 2, x ∈ X)
    @test sys == ConstrainedBlackBoxDiscreteSystem(f1, 2, X)

    sys = @system(x⁺ = f1(x, u), x ∈ X, u ∈ U, dims = (2, 2))
    @test sys == ConstrainedBlackBoxControlDiscreteSystem(f1, 2, 2, X, U)

    sys = @system(x⁺ = f1(x, u, w), x ∈ X, u ∈ U, w ∈ W, dims = (2, 2, 2))
    @test sys == NoisyConstrainedBlackBoxControlDiscreteSystem(f1, 2, 2, 2, X, U, W)
end

# =======================
# Initial value problems
# =======================
@testset "@system with for an initial-value problem" begin
    # continuous ivp in floating-point
    ivp = @system(x' = -1.0x, x(0) ∈ Interval(-1, 1))
    @test ivp == IVP(LinearContinuousSystem(I(-1.0, 1)), Interval(-1, 1)) &&
          eltype(ivp.s.A) == Float64

    # discrete ivp in floating-point
    ivp = @system(x⁺ = -x, x(0) ∈ [1])
    @test ivp == IVP(LinearDiscreteSystem(I(-1.0, 1)), [1]) &&
          eltype(ivp.s.A) == Float64

    # initial state assignment doesn't match state variable
    @test_throws ArgumentError @system(x' = -x, t(0) ∈ Interval(-1.0, 1.0))
end


text:
export feature_importance

using ..JLBoostTrees: has_children

using DataFrames: DataFrame, leftjoin

"""
    feature_importance(jlt::JLBoostTree, df, loss, target)
    feature_importance(jlt::JLBoostTreeModel, df)
    feature_importance(jlt::JLBoostTreeModel, X, y::AbstractVector)

Return the feature of the tree computed on df
"""
feature_importance(jlt::JLBoostTreeModel, X, y::AbstractVector) = begin
    feature_importance(jlt, hcat(X, DataFrame(jlt.target => y)))
end

feature_importance(jlt::JLBoostTreeModel, df) = begin
    if (typeof(df[!, jlt.target]) <: CategoricalVector) & (jlt.loss isa LogitLogLoss)
        dfc = copy(df)
        dfc[!, jlt.target] = 2 .- categorical(dfc[!, jlt.target]).refs
        return feature_importance(trees(jlt), dfc, jlt.loss, jlt.target)
    end
    feature_importance(trees(jlt), df, jlt.loss, jlt.target)
end

add_dict!(d1, d2) = begin
    for (key, val) in d2
        if haskey(d1, key)
            d1[key] += val
        else
            d1[key] = val
        end
    end
    d1
end

dict_to_df(d) = begin
    sv = sum(values(d.freq_dict))
    d2 = Dict(key => (val/sv) for (key, val) in d.freq_dict)
    dfreq = DataFrame(feature = keys(d2) |> collect, Frequency = values(d2) |> collect)

    sv3 = sum(values(d.gain_dict))
    d3 = Dict(key => (val/sv3) for (key, val) in d.gain_dict)
    dgain = DataFrame(feature = keys(d3) |> collect, Quality_Gain = values(d3) |> collect)

    sv4 = sum(values(d.coverage_dict))
    d4 = Dict(key => (val/sv4) for (key, val) in d.coverage_dict)
    dcoverage = DataFrame(feature = keys(d4) |> collect, Coverage = values(d4) |> collect)

    tmp1 = leftjoin(dgain, dcoverage, on = :feature)
    tmp2 = leftjoin(tmp1,  dfreq, on = :feature)

    sort!(tmp2, [:Quality_Gain, :Coverage, :Frequency], rev = true)
end

feature_importance(jlt::AbstractVector{<:AbstractJLBoostTree}, df, loss, target) = begin
    res = feature_importance!.(jlt, Ref(df), Ref(loss), Ref(target))

    # combine
    freq_dict = res[1].freq_dict
    gain_dict = res[1].gain_dict
    coverage_dict = res[1].coverage_dict

    for r in @view(res[2:end])
        add_dict!(freq_dict, r.freq_dict)
        add_dict!(gain_dict, r.gain_dict)
        add_dict!(coverage_dict, r.coverage_dict)
    end
    d = (freq_dict = freq_dict, gain_dict = gain_dict, coverage_dict = coverage_dict)

    dict_to_df(d)
end

feature_importance(jlt::AbstractJLBoostTree, df, loss, target) = begin
    d = feature_importance!(jlt, df, loss, target)
    dict_to_df(d)
end

feature_importance!(jlt::AbstractJLBoostTree, df, loss, target, rows_bool = fill(true, nrow(df)), freq_dict = Dict{Symbol, Int}(), gain_dict = Dict{Symbol, Float64}(), coverage_dict = Dict{Symbol, Float64}(), Gs = JLBoost.g.(loss, getproperty(Tables.columns(df), target), jlt.weight), Hs = JLBoost.h.(loss, getproperty(Tables.columns(df), target), jlt.weight)) = begin
    if has_children(jlt)
        # compute the Quality/Gain. Coverage
        rows_bool_left = rows_bool .& (getproperty(Tables.columns(df), jlt.splitfeature) .<= jlt.split)

        rows_bool_right = rows_bool .& (.!rows_bool_left)

        G_left = sum(@view(Gs[rows_bool_left]))
        H_left = sum(@view(Hs[rows_bool_left]))
        G_right = sum(@view(Gs[rows_bool_right]))
        H_right = sum(@view(Hs[rows_bool_right]))

        # note that hyper parameters are not used to compute the gain
        #gain = G_left^2/H_left + G_right^2/H_right - (G_left + G_right)^2/(H_left + H_right)
        gain = (H_left == 0 ? 0 : G_left^2/H_left) + (H_right == 0 ? 0 : G_right^2/H_right) - (G_left + G_right)^2/(H_left + H_right)
        coverage = H_left + H_right

        if haskey(freq_dict, jlt.splitfeature)
            freq_dict[jlt.splitfeature] += 1
            gain_dict[jlt.splitfeature] += gain
            coverage_dict[jlt.splitfeature] += coverage
        else
            freq_dict[jlt.splitfeature] = 1
            gain_dict[jlt.splitfeature] = gain
            coverage_dict[jlt.splitfeature] = coverage
        end

        feature_importance!(jlt.children[1], df, loss, target, rows_bool_left,  freq_dict, gain_dict, coverage_dict, Gs, Hs)
        feature_importance!(jlt.children[2], df, loss, target, rows_bool_right, freq_dict, gain_dict, coverage_dict, Gs, Hs)
    end

    (freq_dict = freq_dict, gain_dict = gain_dict, coverage_dict = coverage_dict)
end

if false
    feature_importance(tree[1], gmsc, LogitLogLoss(), :SeriousDlqin2yrs)
end


text:
# runtest.jl
push!(LOAD_PATH, joinpath(Pkg.dir("PumiInterface"), "src"))
push!(LOAD_PATH, joinpath(Pkg.dir("SummationByParts"), "src"))
push!(LOAD_PATH, joinpath(Pkg.dir("MeshMovement"), "src"))
push!(LOAD_PATH, joinpath(Pkg.dir("PDESolver"), "src/Utils"))

using FFD  # have FFD initialize MPI
#using PDESolver
using PdePumiInterface
using SummationByParts
using ODLCommonTools
using ArrayViews
#using Utils
using MPI
# using MeshMovement
using Base.Test

# Include the actual tests
include("./test_b-splines.jl")
include("./test_serial_FFD_2D_mesh.jl")
include("test_2d.jl")
include("./test_contractWithdGdB.jl")
include("./test_constraints.jl")
#include("./test_serial_FFD_3D_mesh.jl")

# Debug

#MPI.Finalize()

Base.Test.exitstatus()


text:
# using Test
using Distributed
# import Distributed: RRID, workerpool

include("arraychannels_test.jl")
test_serialise()
test_synchronisation()
test_target_other_channel()
include("interleave_test.jl")
test_interleave()
include("reduce_test.jl")
test_reduce_two()
test_reduce_five()


text:
include("quadrature.jl")
include("recurrence_coefficients.jl")
include("discretization.jl")
include("logistic.jl")
include("comparison_gaussquadrature.jl")
include("quadrature_rules.jl")
include("show.jl")
include("types.jl")
include("auxfuns.jl")
include("densities.jl")
include("polynomial_chaos.jl")
include("multi_indices.jl")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "LLVMOpenMP"
version = v"12.0.1"

sources = [
    ArchiveSource(
        "https://github.com/llvm/llvm-project/releases/download/llvmorg-$(version)/openmp-$(version).src.tar.xz",
        "60fe79440eaa9ebf583a6ea7f81501310388c02754dbe7dc210776014d06b091"
    ),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/openmp-*/
# https://github.com/msys2/MINGW-packages/blob/d440dcb738/mingw-w64-clang/0901-cast-to-make-gcc-happy.patch
atomic_patch -p1 ../patches/0901-cast-to-make-gcc-happy.patch

platform_config=()
if [[ "${target}" == *-freebsd* ]]; then
    platform_config+=(-DCMAKE_SHARED_LINKER_FLAGS="-Wl,--version-script=$(pwd)/runtime/src/exports_so.txt")
elif [[ "${target}" == *-mingw* ]]; then
    # backport https://gitlab.kitware.com/cmake/cmake/-/commit/78f758a463516a78a9ec8d472080c6e61cb89c7f
    sed -i "s@/c  */Fo@-c -Fo@" /usr/share/cmake/Modules/CMakeASM_MASMInformation.cmake
    sed -i "s@libomp_append(asmflags_local /@libomp_append(asmflags_local -@" runtime/cmake/LibompHandleFlags.cmake
    if [[ "${target}" == *x86_64* ]]; then
        platform_config+=(-DLIBOMP_ASMFLAGS="-win64")
    fi
    platform_config+=(-DCMAKE_ASM_MASM_COMPILER="uasm")
elif [[ "${target}" == aarch64-apple-* ]]; then
    # Linking libomp requires the function `__divdc3`, which is implemented in
    # `libclang_rt.osx.a` from LLVM compiler-rt.
    platform_config+=(-DCMAKE_SHARED_LINKER_FLAGS="-L${libdir}/darwin -lclang_rt.osx")
fi
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" \
    -DLIBOMP_INSTALL_ALIASES=OFF \
    "${platform_config[@]}" \
    ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

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

# Dependencies that must be installed before this package can be built
dependencies = [
    # We use UASM only for Windows
    HostBuildDependency(PackageSpec(name="UASM_jll", uuid="bbf38c07-751d-5a2b-a7fc-5c0acd9bd57e")),
    # We need libclang_rt.osx.a for linking libomp, because this library provides the
    # implementation of `__divdc3`.
    BuildDependency(PackageSpec(name="LLVMCompilerRT_jll", uuid="4e17d02c-6bf5-513e-be62-445f41c75a11", version=v"12.0.0"); platforms=[Platform("aarch64", "macos")]),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"8")

max_stars_repo_name: cmhamel/Yggdrasil

max_issues_repo_name: cmhamel/Yggdrasil

max_forks_repo_name: cmhamel/Yggdrasil


text:
# %%
if !("." in LOAD_PATH) push!(LOAD_PATH, ".") end
using Revise
@time using CellularBase
using BenchmarkTools
using Random
# %%
mutable struct SandPile <: AbstractGrid{Int,2}
    state::Array{Int,2}
    otherstate::Array{Int,2}
    h_c::Int
    neighborhood::AbstractNeighborhood
    bc::BoundaryCondition
    function SandPile(state::Array{Int,2}, h_c::Int)
        new(state, similar(state), h_c, VonNeumannNeighborhood{2}(1), FixedMin())
    end
end
@inline iscritical(c::Int, g::SandPile) = c >= g.h_c
@inline iscritical(g::SandPile) = i -> iscritical(i, g)

CellularBase.newstate(s::SandPile) = s.otherstate
function CellularBase.state!(s::SandPile, newstate)
    temp = s.state
    s.state = s.otherstate
    s.otherstate = temp
end

function (_g::SandPile)(current::Int, neighbors::Array{Int, 1})
    h = iscritical(_g)
    return current + h(current)*-4 + count(h, neighbors)
end

CellularBase.possible_states(g::SandPile) = 0:(g.h_c + 3)
# %%

function getindexbc(grid, i)
    if checkbounds(Bool, grid, i)
        @inbounds return grid[i]
    else
        return 0
    end
end
getindexbc(grid) = i -> getindexbc(grid)
function evolveat(grid, hood, index, hc)
    res = grid[index]
    if res >= hc
        res -= 4
    end
    c = count(h -> hc <= getindexbc(grid, h + index), hood)
    res = res + c
    return res
end
function evolve(grid::Array{Int}, hood::Array{CartesianIndex{2},1}, hc::Int; kwargs...)
    simgrid = similar(grid)
    for index in CartesianIndices(grid)
        simgrid[index] = evolveat(grid, hood, index, hc)
    end
    return simgrid
end

function baresimulate(grid::Array{Int}, hc::Int, steps::Int)
    hood = CartesianIndex.([(-1, 0), (1, 0), (0, -1), (0, 1)])
    for i in 1:steps
        grid = evolve(grid, hood, hc)
    end
    return grid
end
function baresimulate!(cbgrid::SandPile, steps)
    for i in 1:steps
        evolve!(cbgrid)
    end
end

function gridtest(cbgrid::SandPile)
    for i = 1:400
        cbgrid(1, [1, 1, 1, 1])
    end
end
function gridtest(ngrid::Array{Int})
    hood = cartesians(VonNeumannNeighborhood{2}(1));
    for i in 1:400
        evolveat(ngrid, hood, CartesianIndex(1,1), 4)
    end
end

CellularBase.boundary_conditions(p::SandPile)::FixedMin = FixedMin()

# %%
println("===============")
Random.seed!(0);
hc = 4;
ngrid = rand(1:hc, 20, 20);
cbgrid = SandPile(copy(ngrid), 4);
hood = cartesians(VonNeumannNeighborhood{2}(1));

@btime baresimulate(ngrid, 4, 5)
# @btime evolve(ngrid, hood, 4)
# @btime evolve(cbgrid)
# @btime CellularBase.evolve!(cbgrid)
@btime simulate!(cbgrid, 5, store_results=false)
# @btime baresimulate!(cbgrid, 1)

# %%


text:
module ViscFoil

include("VorticityStream.jl")
# function integral_momentum_equation(dx, x, p, s)
#     θ, dθ = x                       # State vector
#     ρ_e, u_e, a_e, τ_w, δ_st = p    # Parameters
#     H = δ_st/θ                      # Shape parameter
#     c_f = τ_w/(1/2 * ρ_e * u_e^2)   # Skin-friction coefficient
#     M_e, H = u_e/a_e                # Mach number
#     dx[1] = dθ                      # LHS: dθ/ds
#     dx[2] = c_f / 2  - (H + 2 - M_e^2) * θ/u_e * ??? # RHS: c_f/2  - (H + 2 - Mₑ²) * θ/U_e * dU_e/ds
# end

function edge_velocity(σs, panels, wake_panels, u)
    all_panels          = [panels; wake_panels]
    foil_vortices       = vortex_matrix(panels, panels)
    all_sources         = source_matrix(panels, all_panels)
    wake_wake_vortices  = vortex_matrix(wake_panels, wake_panels)
    wake_foil_vortices  = vortex_matrix(wake_panels, panels)
    wake_boundary       = boundary_vector(colpoints.(wake_panels), u)

    P       = -foil_vortices^(-1) * all_sources
    P_s     =
    # D1  = visc_visc_vortices # PLACEHOLDER
    D1      = wake_foil_vortices * P + wake_wake_vortices
    U_inv   = [ (zeros ∘ length)(panels) ;
                    wake_boundary  ]
    D       = [ P
                D1 ]

    U_e     = U_inv + D * σs
end

end


text:
#
# Implements example 4.2 in Sarkka book, but due to technical reasons with
# regard to handling univariate and multivariate gaussians consistently,
# we simulate two random walks by making A, H, Q, R matrices as opposed to
# scalars specified in book.
#

function Example_4_2(; T = 400, seed = 1)

    # Set random number generator

    rg = MersenneTwister(seed)

    # Matrices in the dynamic model are:

    A = H = diagm(ones(2))
    Q = R = diagm(ones(2))

    # initialise

    m₀, P₀ = zeros(2), diagm(ones(2))

    x₀ = rand(rg, MvNormal(m₀, P₀))

    # store measurements and states

    yclean = Array{Vector{Float64},1}(undef, T)

    y = Array{Vector{Float64},1}(undef, T)

    x = Array{Vector{Float64},1}(undef, T)

    # first state

    x[1]      = A * x₀

    y[1]      = rand(rg, MvNormal(H * x[1], R))

    yclean[1] = H * x₀

    for k in 2:T

        x[k] = rand(rg, MvNormal(A * x[k-1], Q)) # (4.32)

        yclean[k] = H * x[k]

        y[k] = rand(rg, MvNormal(H * (x[k]), R)) # (4.33)

    end

    return yclean, y, x, A, H, Q, R, m₀, P₀

end

function runExample_4_2(seed = 1)

    # Simulate data

    yclean, y, x, A, H, Q, R, m₀, P₀ = Example_4_2(seed=seed)

    # Inference

    μfilter, Σfilter = filteringrecursion(y; A = A, H = H, Q = Q, R = R, m₀ = m₀, P₀ = P₀)

    μsmooth, Σsmooth = smoothingrecursion(y; A = A, H = H, Q = Q, R = R, m = μfilter, P = Σfilter)

    # collect positions

    truepositions = reduce(hcat, [xₖ[1:2] for xₖ in x])

    filtpositions = reduce(hcat, [xₖ[1:2] for xₖ in μfilter])

    smoopositions = reduce(hcat, [xₖ[1:2] for xₖ in μsmooth])

    # plot positions

    figure(1); cla()

    plot(truepositions[1,:], "-", label = "true")

    plot(filtpositions[1,:], "-", label = "filter")

    plot(smoopositions[1,:], "-", label = "smoother")

    legend()

    nothing

end


text:
# using PyPlot

flname = "errors_Elas_Dirichlet2D.dat"

Ps = 1:3         # Range of polynomial order
Ns = [9, 17, 33] # Range of grid size

data   = readdlm(flname; header=true)

data = data[1]

Err_uh1 = fill( 0.0, length(Ps), length(Ns) )
Err_uh2 = fill( 0.0, length(Ps), length(Ns) )
Err_σh1 = fill( 0.0, length(Ps), length(Ns) )
Err_σh2 = fill( 0.0, length(Ps), length(Ns) )
Err_σh4 = fill( 0.0, length(Ps), length(Ns) )
Err_ϵh1 = fill( 0.0, length(Ps), length(Ns) )
Err_ϵh2 = fill( 0.0, length(Ps), length(Ns) )
Err_ϵh4 = fill( 0.0, length(Ps), length(Ns) )

kk = 1
for ii in 1:length(Ps), jj in 1:length(Ns)
  Err_uh1[ii,jj] = data[kk,3]
  Err_uh2[ii,jj] = data[kk,4]
  Err_σh1[ii,jj] = data[kk,5]
  Err_σh2[ii,jj] = data[kk,6]
  Err_σh4[ii,jj] = data[kk,7]
  Err_ϵh1[ii,jj] = data[kk,8]
  Err_ϵh2[ii,jj] = data[kk,9]
  Err_ϵh4[ii,jj] = data[kk,10]
  kk += 1
end

h = 1 ./ ( Ns - 1 )
conv_uh1 = (log.( Err_uh1[:,end-2]) - log.( Err_uh1[:,end] ) ) / (log.( h[end-2]) - log.( h[end] ));
conv_uh2 = (log.( Err_uh2[:,end-2]) - log.( Err_uh2[:,end] ) ) / (log.( h[end-2]) - log.( h[end] ));

conv_σh1 = (log.( Err_σh1[:,end-2]) - log.( Err_σh1[:,end] ) ) / (log.( h[end-2]) - log.( h[end] ));
conv_σh2 = (log.( Err_σh2[:,end-2]) - log.( Err_σh2[:,end] ) ) / (log.( h[end-2]) - log.( h[end] ));
conv_σh4 = (log.( Err_σh4[:,end-2]) - log.( Err_σh4[:,end] ) ) / (log.( h[end-2]) - log.( h[end] ));

conv_ϵh1 = (log.( Err_ϵh1[:,end-2]) - log.( Err_ϵh1[:,end] ) ) / (log.( h[end-2]) - log.( h[end] ));
conv_ϵh2 = (log.( Err_ϵh2[:,end-2]) - log.( Err_ϵh2[:,end] ) ) / (log.( h[end-2]) - log.( h[end] ));
conv_ϵh4 = (log.( Err_ϵh4[:,end-2]) - log.( Err_ϵh4[:,end] ) ) / (log.( h[end-2]) - log.( h[end] ));

clf()
hu1 = loglog(h,Err_uh1)
hu2 = loglog(h,Err_uh2)

hσ1 = loglog(h,Err_σh1)
hσ2 = loglog(h,Err_σh2)
hσ4 = loglog(h,Err_σh4)


text:
struct FileCovariate{T<:AbstractFloat} <: AbstractCovariate{T}
    name::String
    length::Int64
    datapath::String
end

Base.length(covariate::FileCovariate{T}) where {T<:AbstractFloat} = covariate.length

function slice(covariate::FileCovariate{T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:AbstractFloat}
    #len = toobs - fromobs + 1
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    sizeT = sizeof(T)
    iostream = open(covariate.datapath)
    seek(iostream, sizeT * (fromobs - 1))
    buffer = Vector{T}(slicelength)
    map(Seq(Vector{T}, (buffer, iostream, fromobs, toobs, slicelength), nextdatachunk), SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}) do slice
        view(slice, 1:length(slice))
    end
    #data = Vector{UInt8}(sizeof(T) * len)
    # open(covariate.datapath) do f
    #     seek(f, sizeof(T) * (fromobs - 1))
    #     readbytes!(f, data)
    # end
    # covdata = reinterpret(T, data)
    # slice(covdata, 1, len, slicelength)
end

function Base.map(covariate::FileCovariate, dataframe::AbstractDataFrame)
    dataframe[getname(covariate)]
end

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
#!/usr/bin/env julia

#=
# data should be sorted and works by divide and conquer
# run-time complexity O(log n)
=#

function binary_search(search_for, data_array, low=1, high=-1)
  if high == -1
    high = length(data_array)
  end
  if high < low
    return -1
  end
  mid = Int(floor(low + (high - low) / 2))
  if data_array[mid] == search_for
    return mid
  end
  if data_array[mid] > search_for
    return binary_search(search_for, data_array, low, mid-1)
  else
    return binary_search(search_for, data_array, mid+1, high)
  end
end

     ## 1   2   3   4   5   6   7   8   9   10
data = [10, 11, 12, 22, 24, 27, 34, 47, 56, 67]

@time println("at index: ", binary_search(47, data))
@time println("at index: ", binary_search(12, data))
@time println("at index: ", binary_search(7, data))
@time println("at index: ", binary_search(67, data))
@time println("at index: ", binary_search(34, data))


text:
using NonconvexPavito
function convert_common_kwargs(opt::NonconvexPavito.PavitoIpoptCbcAlg, opt_kwargs;
    callback=nothing,
    maxiters=nothing,
    maxtime=nothing,
    abstol=nothing,
    reltol=nothing)

    conv_opt_kwargs = (; opt_kwargs...)

    if !isnothing(callback)
        @warn "common callback argument is currently not used by $(opt)"
    end

    if !isnothing(maxiters)
        @warn "common maxiters argument is currently not used by $(opt)"
    end

    if !isnothing(maxtime)
        conv_opt_kwargs = (; conv_opt_kwargs..., timeout = maxtime)
    end

    if !isnothing(abstol)
        @warn "common abstol argument is currently not used by $(opt)"
    end

    if !isnothing(reltol)
        conv_opt_kwargs = (; conv_opt_kwargs..., rel_gap =reltol)
    end

    return conv_opt_kwargs
end

function _create_options(opt::NonconvexPavito.PavitoIpoptCbcAlg;
    opt_kwargs=nothing,
    sub_options=nothing,
    convergence_criteria=nothing)

    if !isnothing(sub_options)
        options = (; options = !isnothing(opt_kwargs) ? NonconvexPavito.PavitoIpoptCbcOptions(; subsolver_options = IpoptOptions(sub_options...), opt_kwargs...) : NonconvexPavito.PavitoIpoptCbcOptions(subsolver_options = IpoptOptions(sub_options...)))
    else
        options = (; options = !isnothing(opt_kwargs) ? NonconvexPavito.PavitoIpoptCbcOptions(; subsolver_options = IpoptOptions(), opt_kwargs...) : NonconvexPavito.PavitoIpoptCbcOptions(subsolver_options = IpoptOptions()))
    end
    return options
end

check_optimizer_backend(opt::NonconvexPavito.PavitoIpoptCbcAlg) = false