text:
using DrWatson
@quickactivate "ExercisesRepo"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random

# %% phase space plots
fig = figure(figsize = (figx, figx/2))

ax = subplot(231)

ds = Systems.standardmap(k=1)
ic = [[0.1, 0.1], [2.5, 0.4], [1.88, 3.25]]
xs = range(0, stop = 2π, length = 7);
ys = range(0, stop = 2π, length = 3);
ys = copy(xs)
iters = 2000
dataset = Dataset{2, Float64}()
for (i, c) in enumerate(ic)
    tr = trajectory(ds, iters, c)
    ax.scatter(columns(tr)..., s = 5)
end
# ax.set_xticks([0, 2π])
ax.set_xticklabels([])
ax.set_yticklabels([])
# ax.set_yticks([0, 2π])
# ax.set_yticks([0, 2π])
ax.set_yticklabels([])
ax.set_xlabel("\$\\theta\$", labelpad = -5)
ax.set_ylabel("\$p\$",labelpad=0)
ax.set_title("Standard map")

ax = subplot(234)
n = 100
tr = trajectory(ds, n, ic[1])
ax.plot(0:n, tr[:, 2], c = COLORS[1], marker = "o", lw = 1.0)
ax.set_yticks([0, 2π])
ax.set_yticklabels(["0", "2\$\\pi\$"])
ax.set_ylabel("\$p\$",labelpad=-15)
ax.set_xticks([0, n])
ax.set_xlabel("\$n\$", labelpad = -15)

ds = Systems.lorenz()
using3D()
ax = subplot(232, projection = "3d")

ax.set_title("Lorenz-63", pad = 24)

llw = 2.0
tr = trajectory(ds, 100; Ttr = 100, dt = 0.01)
ax.plot3D(columns(tr)..., color = COLORS[1], lw = llw)
tr = trajectory(ds, 10; Ttr = 1000, dt = 0.01)
ax.plot3D(columns(tr)..., color = COLORS[2], lw = llw)
tr = trajectory(ds, 10, [10,20,40.0]; Ttr = 10, dt = 0.01)
ax.plot3D(columns(tr)..., color = COLORS[3], lw = llw)
ax.set_xlabel("\$x\$", labelpad = 0)
ax.set_ylabel("\$y\$", labelpad = 0)
ax.set_zlabel("\$z\$", labelpad = 0)
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.set_zticklabels([])
# ax.set_xticks([-20,20])
# ax.set_yticks([-20,20])
# ax.set_zticks([10,40])

tn = 30
ax = subplot(235)
tr = trajectory(ds, tn; Ttr = 40, dt = 0.01)
ax.plot(0:0.01:tn, tr[:, 1], c = COLORS[1])

ax.set_yticks([-15,15])
ax.set_ylabel("\$x\$", labelpad = -30)
ax.set_xticks([0, tn])
ax.set_xlabel("\$t\$", labelpad = -15)

ds = Systems.henonheiles()

u0s = [[0.0, -0.25, 0.42081, 0.0],
[0.0, 0.1, 0.5, 0.0],
[0.0, -0.31596, 0.354461, 0.0591255]]

ax = subplot(233, projection = "3d")
ax.clear()
ax.set_title("Hénon–Heiles", pad = 24)
for (i, u0) in enumerate(u0s)
    tr = trajectory(ds, 50, u0; Ttr = 50)
    ax.plot3D(tr[:, 1], tr[:, 2], tr[:, 3], c = COLORS[i], alpha = 1 -(i-1)*0.2)
end
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.set_zticklabels([])
ax.set_xlabel("\$x\$")
ax.set_ylabel("\$y\$")
ax.set_zlabel("\$p_x\$",labelpad=-5)

ax = subplot(236)
n = 150
tr = trajectory(ds, n, [0.0, -0.25, 0.48, 0.0]; Ttr = 100)
ax.plot(0:0.01:n, tr[:, 2], c = COLORS[1])
ax.set_xticks([0, n])
ax.set_yticks([-0.5, 0.5])
ax.set_ylabel("\$y\$", labelpad = -30)
ax.set_xlabel("\$t\$", labelpad = -15)

fig.tight_layout()
fig.subplots_adjust(top = 0.93, wspace = 0.4, hspace = 0.25, bottom = 0.1, left = 0.05)
wsave(plotsdir("trajectories"), fig)

# %% Magnetic pendulum
using LinearAlgebra
ma = Systems.magnetic_pendulum(α=0.2, ω=1.0, d=0.3)
∫ = integrator(ma)

g = range(-5, 5; length = 1500)
c = fill(0, length(g), length(g))
t = similar(c, Float64)
@show length(c)

for (i, x) ∈ enumerate(g)
    @show x
    for (j, y) in enumerate(g)
    reinit!(∫, SVector(x, y, 0, 0))
    t0 = ∫.t0
    step!(∫, 100.0)
    while ∫.u[3]^2 + ∫.u[4]^2 > 1e-3
        step!(∫)
    end
    s = SVector(∫.u[1], ∫.u[2])
    k = findmin([(s-m)⋅(s-m) for m in ma.f.magnets])[2]
    # scatter(x, y, color = "C$(k-1)", s = 1)
    t[i,j] = ∫.t - ∫.t0
    c[i,j] = k
    # break
end
end

cd(@__DIR__)
using FileIO
save("magneticpendulum.bson", Dict(:c => c, :g => g, :t => t))

# %% Magnetic pendulum zoomed
using LinearAlgebra
ma = Systems.magnetic_pendulum(α=0.2, ω=1.0, d=0.3)
∫ = integrator(ma)

gx = range(1.80, 1.95; length = 1000)
gy = range(0, 0.12; length = 1000)
c = fill(0, length(gx), length(gy))
t = similar(c, Float64)
@show length(c)

for (i, x) ∈ enumerate(gx)
    @show x
    for (j, y) in enumerate(gy)
        reinit!(∫, SVector(x, y, 0, 0))
        t0 = ∫.t0
        step!(∫, 100.0)
        while ∫.u[3]^2 + ∫.u[4]^2 > 1e-3
            step!(∫)
        end
        s = SVector(∫.u[1], ∫.u[2])
        k = findmin([(s-m)⋅(s-m) for m in ma.f.magnets])[2]
        # scatter(x, y, color = "C$(k-1)", s = 1)
        t[i,j] = ∫.t - ∫.t0
        c[i,j] = k
        # break
    end
end

cd(@__DIR__)
using FileIO
save("magneticpendulum_zoom.bson", Dict(:c => c, :gx => gx, :gy => gy, :t => t))

# %% Plot this
ma = Systems.magnetic_pendulum(α=0.2, ω=1.0, d=0.3)

cd(@__DIR__)
using FileIO
a = load("magneticpendulum.bson")
g = a[:g]; c = a[:c]

LC =  matplotlib.colors.ListedColormap
cmap = LC([matplotlib.colors.to_rgb("C$k") for k in 0:2])

fig = figure(figsize=(figx/2, figx/2))
pcolormesh(g, g, c'; cmap = cmap, shading = "gouraud")
gca().set_aspect("equal")
xticks([-5, 5])
yticks([-5, 5])
xlabel("\$x\$", labelpad=-30)
ylabel("\$y\$", labelpad=-30)
for m in ma.f.magnets
    scatter(m[1], m[2]; color = "white", edgecolor = "black", zorder = 99, s = 100)
end
tight_layout()
subplots_adjust(bottom = 0.1, top = 0.95)

# plot rando initial condition
# tr = trajectory(ma, 10000, [-1.91228, -3.59712, 0, 0]; dt = 0.1)
# x, y = columns(tr)
# plot(x, y, "C3", lw = 1, alpha = 0.5, zorder = 1)
fsave(joinpath(figdir, "magneticpendulum"), fig)

az = load("magneticpendulum_zoom.bson")
gx = az[:gx]; gy = az[:gy]; c = az[:c]

LC =  matplotlib.colors.ListedColormap
cmap = LC([matplotlib.colors.to_rgb("C$k") for k in 0:2])

fig = figure(figsize=(figx/2, figx/2))
pcolormesh(gx, gy, c'; cmap = cmap, shading = "gouraud")
gca().set_aspect("equal")

xticks([1.8, 1.95], size = 20)
yticks([0, 0.12], size = 20)

fig.savefig(joinpath(figdir, "magneticpendulum_zoom.png"))

# xticks([-5, 5])
# yticks([-5, 5])

max_stars_repo_name: moritz-g/ExercisesRepo

max_issues_repo_name: moritz-g/ExercisesRepo

max_forks_repo_name: moritz-g/ExercisesRepo


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/optimize/optimize_convex.jl
# Contains the solve_local_nlp! routine which computes the optimal value
# of a convex function. This is used to compute the upper bound in the
# branch and bound routine. A number of utility function required for
# solve_local_nlp! are also included.
#############################################################################

"""

Shifts the resulting local nlp objective value `f*` by `(1.0 + relative_tolerance/100.0)*f* + absolute_tolerance/100.0`.
This assumes that the local solvers relative tolerance and absolute tolerance is significantly lower than the global
tolerance (local problem is minimum).
"""
function stored_adjusted_upper_bound!(d::GlobalOptimizer, v::Float64)
    adj_atol = d._parameters.absolute_tolerance/100.0
    adj_rtol = d._parameters.relative_tolerance/100.0
    if v > 0.0
        d._upper_objective_value = v*(1.0 + adj_rtol) + adj_atol
    else
        d._upper_objective_value = v*(1.0 - adj_rtol) + adj_atol
    end
    return nothing
end

function _update_upper_variables!(d, m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType}
    for i = 1:_variable_num(FullVar(), m)
        v = m._upper_variables[i]
        l  = _lower_bound(FullVar(), m, i)
        u  = _upper_bound(FullVar(), m, i)
        if is_integer(FullVar(), m, i)
            l = ceil(l)
            u = floor(u)
        end
        is_fixed_int = l == u
        vi = _working_variable_info(m,i)
        if is_fixed(vi) || is_fixed_int
            MOI.add_constraint(d, v, ET(l))
        elseif is_less_than(vi)
            MOI.add_constraint(d, v, LT(u))
        elseif is_greater_than(vi)
            MOI.add_constraint(d, v, GT(l))
        elseif is_real_interval(vi)
            MOI.add_constraint(d, v, LT(u))
            MOI.add_constraint(d, v, GT(l))
        end
    end
    return
end

function _finite_mid(l::T, u::T) where T
    (isfinite(l) && isfinite(u)) && return 0.5*(l + u)
    isfinite(l) ? l : (isfinite(u) ? u : zero(T))
end
function _set_starting_point!(d, m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType}
    for i = 1:_variable_num(FullVar(), m)
        l  = _lower_bound(FullVar(), m, i)
        u  = _upper_bound(FullVar(), m, i)
        v = m._upper_variables[i]
        MOI.set(d, MOI.VariablePrimalStart(), v, _finite_mid(l, u))
    end
    return
end

"""
    LocalResultStatus

Status code used internally to determine how to interpret the results from the
solution of a local problem solve.
"""
@enum(LocalResultStatus, LRS_FEASIBLE, LRS_OTHER)

"""
$(SIGNATURES)

Takes an `MOI.TerminationStatusCode` and a `MOI.ResultStatusCode` and returns `true`
if this corresponds to a solution that is proven to be feasible.
Returns `false` otherwise.
"""
function local_problem_status(t::MOI.TerminationStatusCode, r::MOI.ResultStatusCode)

    if (t == MOI.OPTIMAL) && (r == MOI.FEASIBLE_POINT)
        return LRS_FEASIBLE
    elseif (t == MOI.LOCALLY_SOLVED) && (r == MOI.FEASIBLE_POINT)
        return LRS_FEASIBLE
    # This is default solver specific... the acceptable constraint tolerances
    # are set to the same values as the basic tolerance. As a result, an
    # acceptably solved solution is feasible but non necessarily optimal
    # so it should be treated as a feasible point
    elseif (t == MOI.ALMOST_LOCALLY_SOLVED) && (r == MOI.NEARLY_FEASIBLE_POINT)
        return LRS_FEASIBLE
    end
    return LRS_OTHER
end

function _unpack_local_nlp_solve!(m::GlobalOptimizer, d::T) where T

    tstatus = MOI.get(d, MOI.TerminationStatus())
    pstatus = MOI.get(d, MOI.PrimalStatus())
    m._upper_termination_status = tstatus
    m._upper_result_status = pstatus

    if local_problem_status(tstatus, pstatus) == LRS_FEASIBLE
        if is_integer_feasible(m)

            m._upper_feasibility = true
            obj_val = MOI.get(d, MOI.ObjectiveValue())
            stored_adjusted_upper_bound!(m, obj_val)
            m._best_upper_value = min(obj_val, m._best_upper_value)
            m._upper_solution .= MOI.get(d, MOI.VariablePrimal(), m._upper_variables)

            ip = m._input_problem
            _extract_primal_linear!(d, ip)
            _extract_primal_quadratic!(d, ip)
        end
    else
        m._upper_feasibility = false
        m._upper_objective_value = Inf
    end
    return
end

"""

Constructs and solves the problem locally on on node `y` updated the upper
solution informaton in the optimizer.
"""
function solve_local_nlp!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType}

    upper_optimizer = _upper_optimizer(m)
    MOI.empty!(upper_optimizer)

    for i = 1:m._working_problem._variable_count
        m._upper_variables[i] = MOI.add_variable(upper_optimizer)
    end
    _update_upper_variables!(upper_optimizer, m)
    _set_starting_point!(upper_optimizer, m)

    # add constraints
    ip = m._input_problem
    _add_constraint_store_ci_linear!(upper_optimizer, ip)
    _add_constraint_store_ci_quadratic!(upper_optimizer, ip)
     #add_soc_constraints!(m, upper_optimizer)

    # Add nonlinear evaluation block
    MOI.set(upper_optimizer, MOI.NLPBlock(), m._working_problem._nlp_data)
    MOI.set(upper_optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(upper_optimizer, MOI.ObjectiveFunction{SAF}(), m._working_problem._objective_saf)

    # Optimizes the object
    MOI.optimize!(upper_optimizer)

    #println(" ")
    #@show m._working_problem._objective_saf
    #@show upper_optimizer.optimizer.model.quadratic_le_constraints[1]
    #@show typeof(upper_optimizer.optimizer.model.quadratic_le_constraints[1])
    #@show upper_optimizer.optimizer.model.quadratic_ge_constraints
    #@show upper_optimizer.optimizer.model.quadratic_eq_constraints

    #@show upper_optimizer.optimizer.model.nlp_data
    #@show fieldnames(typeof(upper_optimizer.optimizer.model.nlp_data))

    _unpack_local_nlp_solve!(m, upper_optimizer)
end

function optimize!(::DIFF_CVX, m::GlobalOptimizer)
    solve_local_nlp!(m)
    return
end

ext: jl

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Plots
import .Frames # The dot is due to the fact that Frames is a submodule

default(markersize=1)

@userplot TrajPlot
"""`trajplot(T,args...)` Plots a trajectory of T-matrices in a single plot"""
trajplot
@recipe function trajplot(t::TrajPlot)
    title --> "Trajectory"
    layout --> 1
    T = t.args[1]
    label --> ["x" "y" "z"]
    squeeze(T[1:3,4,:],2)'
end

@userplot TrajPlot3
"""`plot_traj3(T, ls=".", plotFrame = 0)` Plots a trajectory of T-matrices in a 3D plot, with otional frames drawn of length `plotFrame`"""
trajplot3
@recipe function trajplot3(t::TrajPlot3, plotFrame = 0)
    T = t.args[1]
    xguide := "x"
    yguide := "y"
    zguide := "z"
    @series begin
        seriestype := :scatter3d
        (squeeze(T[1,4,:],(1,2)),squeeze(T[2,4,:],(1,2)),squeeze(T[3,4,:],(1,2)))
    end
    if plotFrame > 0
        for i = 1:size(T,3)
            @series begin
                plot!(T[:,:,i],plotFrame)
            end
        end
    end
end

@userplot Plot3Smart
"""`plot3smart(x,args...)` Makes a 3d plot of a matrix"""
plot3smart
@recipe function f(t::Plot3Smart)
    data = t.args[1]
    seriestype := :path3d
    (data[:,1][:],data[:,2][:],data[:,3][:])
end

# @recipe function f(::Type{Val{:mat3d}}, data)
#     @show typeof(data)
#     @show typeof(d[:x])
#     @show typeof(d[:y])
#     @show typeof(d[:z])
#     x := d[:y][:,1][:]
#     y := d[:y][:,2][:]
#     z := d[:y][:,3][:]
#     seriestype := :path3d
#     ()
# end


text:
function Pad5D(a,N1,N2,N3,N4,N5)

	n1,n2,n3,n4,n5 = size(a)
	b = zeros(Float32,N1,N2,N3,N4,N5)
	b[1:n1,1:n2,1:n3,1:n4,1:n5] = a

	return b
end


text:
function reindex(trian::Triangulation,indices)
  TriangulationPortion(trian,collect(Int,indices))
end

function reindex(trian::Triangulation,indices::Vector{Int})
  TriangulationPortion(trian,indices)
end

function reindex(trian::Triangulation,indices::IdentityVector)
  trian
end

"""
"""
struct TriangulationPortion{Dc,Dp,G} <: Triangulation{Dc,Dp}
  oldtrian::G
  cell_to_oldcell::Vector{Int}
  @doc """
  """
  function TriangulationPortion(oldtrian::Triangulation{Dc,Dp},cell_to_oldcell::Vector{Int}) where {Dc,Dp}
    new{Dc,Dp,typeof(oldtrian)}(oldtrian,cell_to_oldcell)
  end
end

function get_reffes(trian::TriangulationPortion)
  get_reffes(trian.oldtrian)
end

function get_cell_type(trian::TriangulationPortion)
  reindex(get_cell_type(trian.oldtrian),trian.cell_to_oldcell)
end

function get_cell_coordinates(trian::TriangulationPortion)
  reindex(get_cell_coordinates(trian.oldtrian),trian.cell_to_oldcell)
end

function get_cell_map(trian::TriangulationPortion)
  cell_map = get_cell_map(trian.oldtrian)
  reindex(cell_map,trian.cell_to_oldcell)
end

function get_cell_id(trian::TriangulationPortion)
  reindex(get_cell_id(trian.oldtrian),trian.cell_to_oldcell)
end

function restrict(f::AbstractArray,trian::TriangulationPortion)
  reindex(restrict(f,trian.oldtrian),trian.cell_to_oldcell)
end

function get_normal_vector(trian::TriangulationPortion)
  reindex(get_normal_vector(trian.oldtrian),trian.cell_to_oldcell)
end


text:
export Gridded

immutable Gridded{D<:Degree} <: InterpolationType end
Gridded{D<:Degree}(::D) = Gridded{D}()

griddedtype{D<:Degree}(::Type{Gridded{D}}) = D

typealias GridIndex{T} Union{AbstractVector{T}, Tuple}

# Because Ranges check bounds on getindex, it's actually faster to convert the
# knots to Vectors. It's also good to take a copy, so it doesn't get modified later.
immutable GriddedInterpolation{T,N,TCoefs,IT<:DimSpec{Gridded},K<:Tuple{Vararg{Vector}},pad} <: AbstractInterpolation{T,N,IT,OnGrid}
    knots::K
    coefs::Array{TCoefs,N}
end
function GriddedInterpolation{N,TCoefs,TWeights<:Real,IT<:DimSpec{Gridded},pad}(::Type{TWeights}, knots::NTuple{N,GridIndex}, A::AbstractArray{TCoefs,N}, ::IT, ::Val{pad})
    isleaftype(IT) || error("The b-spline type must be a leaf type (was $IT)")
    isleaftype(TCoefs) || warn("For performance reasons, consider using an array of a concrete type (eltype(A) == $(eltype(A)))")

    knts = mapcollect(knots...)
    for (d,k) in enumerate(knts)
        length(k) == size(A, d) || throw(DimensionMismatch("knot vectors must have the same number of elements as the corresponding dimension of the array"))
        length(k) == 1 && error("dimensions of length 1 not yet supported")  # FIXME
        issorted(k) || error("knot-vectors must be sorted in increasing order")
        iextract(IT, d) != NoInterp || k == collect(1:size(A, d)) || error("knot-vector should be the range 1:$(size(A,d)) for the method Gridded{NoInterp}")
    end
    c = one(TWeights)
    for _ in 2:N
        c *= c
    end
    T = typeof(c * one(TCoefs))

    GriddedInterpolation{T,N,TCoefs,IT,typeof(knts),pad}(knts, A)
end

# A type-stable version of map(collect, knots)
mapcollect() = ()
@inline mapcollect(k::AbstractVector) = (collect(k),)
@inline mapcollect(k1::AbstractVector, k2::AbstractVector...) = (collect(k1), mapcollect(k2...)...)

# Utilities for working either with scalars or tuples/tuple-types
iextract{T<:Gridded}(::Type{T}, d) = T
iextract{T<:GridType}(::Type{T}, d) = T

@generated function size{T,N,TCoefs,IT,K,pad}(itp::GriddedInterpolation{T,N,TCoefs,IT,K,pad}, d)
    quote
        d <= $N ? size(itp.coefs, d) - 2*padextract($pad, d) : 1
    end
end

function interpolate{TWeights,TCoefs,Tel,N,IT<:DimSpec{Gridded}}(::Type{TWeights}, ::Type{TCoefs}, knots::NTuple{N,GridIndex}, A::AbstractArray{Tel,N}, it::IT)
    GriddedInterpolation(TWeights, knots, A, it, Val{0}())
end
function interpolate{Tel,N,IT<:DimSpec{Gridded}}(knots::NTuple{N,GridIndex}, A::AbstractArray{Tel,N}, it::IT)
    interpolate(tweight(A), tcoef(A), knots, A, it)
end

interpolate!{TWeights,Tel,N,IT<:DimSpec{Gridded}}(::Type{TWeights}, knots::NTuple{N,GridIndex}, A::AbstractArray{Tel,N}, it::IT) = GriddedInterpolation(TWeights, knots, A, it, Val{0}())
function interpolate!{Tel,N,IT<:DimSpec{Gridded}}(knots::NTuple{N,GridIndex}, A::AbstractArray{Tel,N}, it::IT)
    interpolate!(tweight(A), tcoef(A), knots, A, it)
end

lbound(itp::GriddedInterpolation, d) = itp.knots[d][1]
ubound(itp::GriddedInterpolation, d) = itp.knots[d][end]

include("constant.jl")
include("linear.jl")
include("indexing.jl")


text:
using Images
using ImageView
using LinearAlgebra
using Plots

t = -3:.01:3

Utrue = [cos.(17*t).*exp.(-t.^2) sin.(11*t)]
Strue = [2.0 0.0; 0.0 0.5]
Vtrue = [sin.(5*t).*exp.(-t.^2) cos.(13*t)]

X = Utrue*Strue*Vtrue'
imshow(X)

sigma = 1
Xnoisy = X+sigma*randn(size(X))
imshow(Xnoisy)

U, S, V = svd(Xnoisy)

N = size(Xnoisy)[1]
cutoff = (4/sqrt(3))*sqrt(N)*sigma # Hard threshold
r = sum(S .> cutoff) # Keep modes w/ sig > cutoff
Xclean = U[:,1:r]*diagm(S[1:r])*V[:,1:r]'
imshow(Xclean)

cdS = cumsum(S)/sum(S) # Cumulative energy
r90 = sum(cdS .> 0.90) # Find r to capture 90% energy

X90 = U[:,1:r90]*diagm(S[1:r90])*V[:,1:r90]'
imshow(X90)

# plot singular values

p1 = plot(S, yaxis = :log, legend = false)
p2 = plot(S[1:r], yaxis = :log, legend = false)

p3 = plot(cdS, legend = false)
p4 = plot(cdS[1:r90], legend = false)
p5 = plot(cdS[1:r], legend = false)


text:
function observe(events::Events{T},
                 delay_infection::UnivariateDistribution,
                 delay_removal::UnivariateDistribution;
                 force::Bool = false) where T <: Union{SEIR, SIR}
  infection = fill(NaN, events.individuals)
  removal = fill(NaN, events.individuals)
  if force
    for i in findall(.!isnan.(events.infection))
      if events.infection[i] == -Inf
        infection[i] = -Inf
        if events.removal[i] == -Inf
          removal[i] = -Inf
        end
      elseif isnan(events.removal[i])
        infection[i] = events.infection[i] + rand(delay_infection)
      elseif events.removal[i] > -Inf
        infection_delay_ub = events.removal[i] - events.infection[i]
        infection[i] = events.infection[i] +
                       rand(Truncated(delay_infection, 0.0, infection_delay_ub))
        removal[i] = events.removal[i] + rand(delay_removal)
      end
      @debug "Infection observation of i = $i at t = $(round(infection[i], digits=3)) (actual infection onset at t = $(round(events.infection[i], digits=3)))"
      @debug "Removal observation of i = $i at t = $(round(removal[i], digits=3)) (actual removal at t = $(round(events.removal[i], digits=3)))"
    end
  else
    for i in findall(.!isnan.(events.infection))
      if events.infection[i] == -Inf
        infection[i] = -Inf
        if events.removal == -Inf
          removal[i] = -Inf
        end
      else
        infection_delay = rand(delay_infection)
        if isnan(events.removal[i])
          infection[i] = events.infection[i] + infection_delay
        else
          if infection_delay + events.infection[i] < events.removal[i]
            infection[i] = events.infection[i] + infection_delay
            removal[i] = events.removal[i] + rand(delay_removal)
          end
        end
      end
      @debug "Infection observation of i = $i t = $(round(infection[i], digits=3)) (actual infection onset at t = $(round(events.infection[i], digits=3)))"
      @debug "Removal observation of i = $i at t = $(round(removal[i], digits=3)) (actual removal at t = $(round(events.removal[i], digits=3)))"
    end
  end
  return EventObservations{T}(infection, removal)
end

function observe(sim::Simulation{T},
                 delay_infection::UnivariateDistribution,
                 delay_removal::UnivariateDistribution;
                 force::Bool = false) where T <: Union{SEIR, SIR}
  return observe(sim.events, delay_infection, delay_removal, force=force)
end

function observe(events::Events{T},
                 delay_infection::UnivariateDistribution) where T <: Union{SEI, SI}
  infection = fill(NaN, events.individuals)
  @simd for i in findall(.!isnan.(events.infection))
    if events.infection[i] == -Inf
      infection[i] = -Inf
    else
      infection[i] = events.infection[i] + rand(delay_infection)
    end
    @debug "Infection observation of i = $i at t = $(round(infection[i], digits=3))) (actual infection onset at t = $(round(events.infection[i], digits=3)))"
  end
  return EventObservations{T}(infection)
end

function observe(sim::Simulation{T},
                 delay_infection::UnivariateDistribution) where T <: Union{SEI, SI}
  return observe(sim.events, delay_infection)
end


text:
include(joinpath(@__DIR__, "pedagogy_graph.jl"))
using Random

dims = [5, 5] # if > (14, 14) the display becomes all black, don't know why !
periodic = false
H, J = 0.0, 0.01

ising = PRS.Ising(dims, J, H; periodic=periodic)

using Random
rng = Random.MersenneTwister(123)
@time sample = PRS.generate_sample_gibbs_perfect(rng, ising)

path = joinpath(@__DIR__, "output", "ising", "ising_grid_prs.pdf")
p = plot(ising, sample, dims, path; nodelabel=nothing)


text:
__precompile__()

module Sellmeier

const codefile = joinpath(dirname(@__FILE__), "literate_org_tangled_code.jl")
if isfile(codefile)
    include(codefile)
else
    error("Sellmeier not properly installed. Please run Pkg.build(\"Sellmeier\") then restart Julia.")
end

end


text:
path = "docs/src/assets/"

using ISA, Plots
𝐶₀ = AMFMtriplet(t->10t+100,t->25cos(t)+50,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(ψ₀; timeaxis=0.0:0.001:10.0)
png(path*"IS_exAMFM1.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->10t+100,t->25cos(t)+50,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=0.0:0.001:10.0)
png(path*"IS_exAMFM2.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->10t+100,t->25cos(t)+50,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=0.0:0.001:10.0, view="TF",
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"IS_exAMFM3.png")

#--------------------------

using ISA, Plots
𝐶₀ = AMFMtriplet(t->10sin(t),t->t^2+2t+50,1.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(ψ₀; timeaxis=0.0:0.001:10.0)
png(path*"IS_exAMFM4.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->10sin(t),t->t^2+2t+50,1.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=0.0:0.001:10.0)
png(path*"IS_exAMFM5.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->10sin(t),t->t^2+2t+50,1.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=0.0:0.001:10.0, view="TF",
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"IS_exAMFM6.png")

#----------------------------

using ISA, Plots
𝐶₀ = AMFMtriplet(t->5exp(0.1t)+15,t->t^2+t+25,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(ψ₀; timeaxis=0.0:0.001:10.0)
png(path*"IS_exAMFM7.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->5exp(0.1t)+15,t->t^2+t+25,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=0.0:0.001:10.0)
png(path*"IS_exAMFM8.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->5exp(0.1t)+15,t->t^2+t+25,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=0.0:0.001:10.0, view="TF",
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"IS_exAMFM9.png")


text:
using FEASTsolver
using Test

@testset "FEASTsolver.jl" begin
    # Write your tests here.
end


text:
using Test, JuCLI

@testset "JuCLI" begin

cli = CLI("mycli")

@test cli.name == "mycli"
@test CLI().name == "cli"

end # JuCLI


text:
# Node that implements a controller for a DC motor.
# Calling syntax: nodectrl.jl [MQTT server address]
using openbuildnet

# Process input arguments
mqttserver = length(ARGS)>0?String(ARGS[1]):"tcp://localhost:1883"

# Controller's states
x = zeros(3)
command = 0.0

# The parameters
A = [-0.82 1.0 0.82; 1.0 0.0 0.0; 0.0 1.0 0.0]
C = [12.62 -19.75 7.625]

# Create the node
workspace = "test2"
node = OBNNode("ctrl", workspace, mqttserver)

# Create inputs
u_v = create_input(node, "v", Float64)  # Create velocity input
u_sp = create_input(node, "sp", Float64) # Create setpoint input

# Create output
y_u = create_output(node, "u", Float64) # Create command output

# Create the CSV writer
dumpfile = open("controller.txt", "w")

# Create Callbacks
on_init(node) do
    set(y_u, command) # Initialize the output (not necessary)
    println("At ", sim_time(node), " simulation started.")
end

on_term(node) do
    close(dumpfile)  # Close the CSV file
    println("At ", sim_time(node), " simulation terminated.")
end

const MAINBLOCK = 0

on_block_output(node, MAINBLOCK) do
    global command
    command = C * x
    set(y_u, command[1])
end

on_block_state(node, MAINBLOCK) do
    # Get the inputs
    sp = get(u_sp)
    v = get(u_v)
    global x
    x = A * x
    x[1] += 32.0 * (sp - v)

    # Write to the CSV file
    write(dumpfile, join((sim_time(node), sp, v), "\t"), "\t",
          join(command, "\t"), "\t", join(x, "\t"), "\n")
end

println("Ready to run the controller node; please start all other nodes ...")
status = runsafe(node, 60)
println("Simulation stopped with status = ", status)


text:
Fun(::typeof(identity), d::Circle) = Fun(Laurent(d),[d.center,0.,d.radius])

Space(d::PeriodicDomain) = Fourier(d)
Space(d::Circle) = Laurent(d)

## Evaluation

Evaluation(d::PeriodicDomain,x::Number,n...) = Evaluation(Laurent(d),complex(x),n...)

## Definite Integral

DefiniteIntegral(d::PeriodicDomain) = DefiniteIntegral(Laurent(d))
DefiniteLineIntegral(d::PeriodicDomain) = DefiniteLineIntegral(Laurent(d))

## Toeplitz
union_rule(A::Space{<:PeriodicSegment}, B::Space{<:IntervalOrSegment}) =
    union(Space(Interval(domain(A))), B)


text:
c0 = chain(:standard,1,10,Float64,Float64)
c1 = chain(:standard,3,4,Int16,Float64)
c2 = chain(:gradient,2,3,Float32,Float32)
@test numparas(c0) == 1 && numsamples(c0) == 10
@test numparas(c1) == 3 && numsamples(c1) == 4
@test numparas(c2) == 2 && numsamples(c2) == 3

s0 = samples(:base,1,4,Float64,Float64)
s1 = samples(:base,3,2,Int16,Float64)
s2 = samples(:tensor,2,2,Float32,Float32)

#test that you can't add a base sample to a gradient chain
@test_throws MethodError store!(c2,s0,1)

#test that the number of parameters is the same for sample and chain
@test_throws AssertionError store!(c0,s1,1)
@test_throws AssertionError store!(c1,s0,1)

copy!(s0.values,1.0:4.0)
store!(c0,s0,4)
store!(c0,s0,3)
store!(c0,s0,2)
@test c0.values[:,1:3] == s0.values[:,4:-1:2]
@test isequal(c0.loglikelihood[1:3],s0.loglikelihood[4:-1:2])

copy!(s1.values,[1.0 4.0;2.0 5.0;3.0 6.0])
store!(c1,s1,2)
store!(c1,s1,1)
store!(c1,s1,2)
store!(c1,s1,2)
@test c1.values == [4.0 1.0 4.0 4.0;5.0 2.0 5.0 5.0;6.0 3.0 6.0 6.0]
@test isequal(c1.loglikelihood,s1.loglikelihood[[2,1,2,2]])

copy!(s2.values,[1.0 3.0;2.0 4.0])
store!(c2,s2,1)
store!(c2,s2,1)
store!(c2,s2,2)
@test c2.values == [1.0 1.0 3.0;2.0 2.0 4.0]
@test isequal(c2.loglikelihood,s2.loglikelihood[[1,1,2]])
@test isequal(c2.gradloglikelihood,s2.gradloglikelihood[:,[1,1,2]])

#these should give warnings that the chain is full
println("Two warnings should be printed below this line")
store!(c1,s1,1)
store!(c2,s2,1)

@test c0.proposed == 3 && accepted!(c0,IndicatorStationary{Float64}(zeros(3),[1,4,3,2])) == 3
@test c1.proposed == 4 && accepted!(c1,IndicatorStationary{Float32}(zeros(Float32,4),[1,2,1,2,2])) == 3
@test c2.proposed == 3 && accepted!(c2,IndicatorStationary{Float64}(zeros(3),[1,1,1,2])) == 1

immutable ChainTestModel <: AbstractModel
    parameters
end

p0 = parameters([0.0],[2.0])
p1 = parameters([0.0,1.0,2.0],[10.0,5.0,3.0])
p2 = parameters([2.0,2.0],[10.0,10.0])
m0 = ChainTestModel(p0)
m1 = ChainTestModel(p1)
m2 = ChainTestModel(p2)

@test_approx_eq logprior(c0,m0) logprior(p0,c0.values,Float64)
@test_approx_eq logprior(c1,m1) logprior(p1,c1.values,Float64)
@test_approx_eq logprior(c2,m2) logprior(p2,c2.values,Float64)

@test_approx_eq logposterior(c0,m0) loglikelihood(c0) + logprior(c0,m0)
@test_approx_eq logposterior(c1,m1) loglikelihood(c1) + logprior(c1,m1)
@test_approx_eq logposterior(c2,m2) loglikelihood(c2) + logprior(c2,m2)

###Test show()
println()
println("====================")
println("Test show() function")
show(c0)
show(c1)
show(c2)
println("End  show() function")
println("====================")
println()
nothing


text:
function gen_lists!(data, p, first_atom, next_atom)
    for iat in 1:length(p)
        ic = trunc(Int64,p[iat][1]/data.cutoff) + 1
        jc = trunc(Int64,p[iat][2]/data.cutoff) + 1
        next_atom[iat] = first_atom[ic, jc]
        first_atom[ic, jc] = iat
    end
    return first_atom, next_atom
end


text:
module CellDofBasesTests

using Test
using FillArrays
using Gridap.Arrays
using Gridap.Fields
using Gridap.ReferenceFEs
using Gridap.CellData
using Gridap.Fields: MockField, OtherMockBasis
using LinearAlgebra

reffe = QUAD4
ndofs = 3
b = get_shapefuns(reffe)
d = get_dof_basis(reffe)
c = fill(1.0,ndofs)
f = OtherMockBasis{2}(ndofs)

l = 10
bl = fill(b,l)
dl = Fill(d,l)
cl = fill(c,l)
fl = Fill(f,l)
ϕl = lincomb(fl,cl)

ref_style = Val(true)
axs = (Base.OneTo(num_dofs(reffe)),)
cell_axes = Fill(axs,l)
cb = GenericCellField(bl,ϕl,ref_style,cell_axes)
cd = GenericCellDofBasis(ref_style,dl)
test_cell_dof_basis(cd,cb)

@test evaluate(cd,cb) ≈ Fill(Matrix(I,num_dofs(reffe),num_dofs(reffe)),l)

end # module


text:
using DataFrames
import CSV, SQLite

include("paths.jl")

function readdata()
    beers = DataFrame(CSV.File(path_beers))[!, [:brewery_id, :name]]
    breweries = DataFrame(CSV.File(path_breweries))[!, [:id, :name]]
    geocodes = DataFrame(CSV.File(path_geocodes))[!, [:brewery_id, :latitude, :longitude]]

    locations = innerjoin(breweries, geocodes, on = :id => :brewery_id)
    return beers, locations
end

function createDB()
    beers, locations = readdata()
    db = SQLite.DB(path_database)
    beers |> SQLite.load!(db, path_db_beers)
    locations |> SQLite.load!(db, path_db_locations)
    return
end


text:
module Soundex
    export soundex!, soundex

    function ispseudovowel(chr::Char)
        return contains(['a', 'e', 'i', 'o', 'u', 'y', 'h', 'w'], chr)
    end

    isvowel(chr::Char) = contains(['a', 'e', 'i', 'o', 'u', 'y'], chr)

    function encode(chr::Char)
        if contains(['b', 'f', 'p', 'v'], chr)
            return '1'
        elseif contains(['c', 'g', 'j', 'k', 'q', 's', 'x', 'z'], chr)
            return '2'
        elseif contains(['d', 't'], chr)
            return '3'
        elseif 'l' == chr
            return '4'
        elseif contains(['m', 'n'], chr)
            return '5'
        elseif 'r' == chr
            return '6'
        elseif ispseudovowel(chr)
            return '9'
        else
            msg = @sprintf "Unknown character encountered %s" chr
            throw(ArgumentError(msg))
        end
    end

    function soundex!(output::Array{Char}, s::String)
        index = 0
        previous_chr, previous_encoding = char(0), '9'
        for chr in s
            chr = lowercase(chr)
            encoding = encode(chr)
            if index == 0
                index += 1
                output[index] = uppercase(chr)
            elseif index <= 3
                # Vowels along with 'h' and 'w' are not encoded
                if !ispseudovowel(chr)
                    # If two or more letters with the same number are adjacent in
                    #  the original name, only retain the first
                    if encoding != previous_encoding
                        if encoding == output[index]
                            # Two letters with the same number separated by a
                            # vowel are coded as two numbers
                            #
                            # Two letters with the same number separated by
                            # 'h' or 'w' are coded as a single number
                            if isvowel(previous_chr)
                                index += 1
                                output[index] = encoding
                            end
                        else
                            index += 1
                            output[index] = encoding
                        end
                    end
                end
            else
                # Only retain the first four encoded letters
                break
            end
            previous_chr, previous_encoding = chr, encoding
        end
        # Append with zeros if needed
        if index < 4
            for i in (index + 1):4
                output[i] = '0'
            end
        end
    end

    function soundex(s::String)
        output = Array(Char, 4)
        soundex!(output, s)
        return CharString(output)
    end
end # module


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -9.75097970893739 <= q <= 8.436779502208893)

                     add_NL_constraint(m, :(gelu(0.5456914178120993 + -0.994137029231358*gelu(0.4043679177252084 + -0.9159920740102909*$(x[1]) + 0.630569400327222*$(x[2]) + 0.15548796895771444*$(x[3])) + 0.38451384631825425*gelu(-0.7218630242098185 + -0.7615826451198555*$(x[1]) + 0.556650647741431*$(x[2]) + -0.7074053991868854*$(x[3])) + 0.9673084611668079*gelu(-0.4849043710003591 + 0.9895364493832828*$(x[1]) + -0.17212002157434236*$(x[2]) + -0.9444649970468406*$(x[3]))) + gelu(0.6143810270772536 + 0.32603207961941205*gelu(0.4043679177252084 + -0.9159920740102909*$(x[1]) + 0.630569400327222*$(x[2]) + 0.15548796895771444*$(x[3])) + 0.7660966505536861*gelu(-0.7218630242098185 + -0.7615826451198555*$(x[1]) + 0.556650647741431*$(x[2]) + -0.7074053991868854*$(x[3])) + 0.30918160231413605*gelu(-0.4849043710003591 + 0.9895364493832828*$(x[1]) + -0.17212002157434236*$(x[2]) + -0.9444649970468406*$(x[3]))) + gelu(-0.4079478457750749 + 0.06443443960646311*gelu(0.4043679177252084 + -0.9159920740102909*$(x[1]) + 0.630569400327222*$(x[2]) + 0.15548796895771444*$(x[3])) + 0.09351478982242645*gelu(-0.7218630242098185 + -0.7615826451198555*$(x[1]) + 0.556650647741431*$(x[2]) + -0.7074053991868854*$(x[3])) + -0.7258007123919756*gelu(-0.4849043710003591 + 0.9895364493832828*$(x[1]) + -0.17212002157434236*$(x[2]) + -0.9444649970468406*$(x[3]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m


text:
include("../src/SocpSolver.jl")
using .SocpSolver
import Convex
using Test

@time @testset "Convex Problem Depot tests" begin
    Convex.ProblemDepot.run_tests(;  exclude=[r"mip", r"exp", r"sdp"]) do problem
        Convex.solve!(problem, () -> SocpSolver.Optimizer(verbose=false))
    end
end;


text:
using HighDimPDE
using Random
using Test
using Flux
using Statistics
using CUDA
CUDA.functional() ? use_cuda = true : use_cuda = false
# import DiffEqFlux: FastChain, FastDense

#relative error l2
function rel_error_l2(u, uanal)
    if abs(uanal) >= 10 * eps(eltype(uanal))
        sqrt((u - uanal)^2 / u^2)
    else # overflow
        abs(u-uanal)
    end
end

# https://en.wikipedia.org/wiki/Heat_equation#Fundamental_solutions
@testset "DeepSplitting - heat equation - single point" begin # solving at one unique point
    for d in [1, 10, 50]
        println("test for d = ", d)
        tspan = (0f0, 5f-1)
        dt = 5f-1  # time step
        μ(x, p, t) = 0f0 # advection coefficients
        σ(x, p, t) = 1f-1 # diffusion coefficients

        u_anal(x, t) = sum(x.^2) + d * σ(0., 0., 0.)^2 * t
        g(x) = sum(x.^2, dims=1)

        # d = 10
        x0 = fill(2f0,d)
        hls = d + 10 #hidden layer size

        nn = Flux.Chain(Dense(d,hls,relu),
                        Dense(hls,hls,relu),
                        Dense(hls,1)) # Neural network used by the scheme

        opt = ADAM(0.01) #optimiser
        alg = DeepSplitting(nn, opt = opt)

        f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = 0f0 .* v_y

        # defining the problem
        prob = PIDEProblem(g, f, μ, σ, tspan,
                            x = x0
                            )
        # solving
        xs,ts,sol = solve(prob, alg, dt,
                        verbose = false,
                        use_cuda = use_cuda,
                        maxiters = 1000,
                        batch_size = 10000)
        u1 = sol[end]
        u1_anal = u_anal(x0, tspan[end])
        e_l2 = rel_error_l2(u1, u1_anal)
        println("rel_error_l2 = ", e_l2, "\n")
        @test e_l2 < 0.1
    end
end

@testset "DeepSplitting - heat equation - interval" begin # solving on interval [-5f-1, 5f-1]^d
    for d in [1, 10, 50]
        println("test for d = ", d)
        tspan = (0f0, 5f-1)
        dt = 5f-1  # time step
        μ(x, p, t) = 0f0 # advection coefficients
        σ(x, p, t) = 1f-1 # diffusion coefficients

        u_anal(x, t) = sum(x.^2) + d * σ(0., 0., 0.)^2 * t
        g(x) = sum(x.^2, dims=1)

        # d = 10
        u_domain = (fill(-5f-1, d), fill(5f-1, d))
        hls = d + 10 #hidden layer size

        nn = Flux.Chain(Dense(d,hls,relu),
                        Dense(hls,hls,relu),
                        Dense(hls,1)) # Neural network used by the scheme

        opt = ADAM(0.01) #optimiser
        alg = DeepSplitting(nn, opt = opt)

        f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = 0f0 .* v_y #TODO: this fix is not nice

        # defining the problem
        prob = PIDEProblem(g, f, μ, σ, tspan, u_domain = u_domain)
        # solving
        xs,ts,sol = solve(prob, alg, dt,
                        verbose = false,
                        use_cuda = use_cuda,
                        maxiters = 1000,
                        batch_size = 20000)
        u1 = [sol[end](x)[] for x in xs]
        u1_anal = u_anal.(xs, tspan[end])
        e_l2 = mean(rel_error_l2.(u1, u1_anal))
        println("rel_error_l2 = ", e_l2, "\n")
        @test e_l2 < 0.1
    end
end

@testset "DeepSplitting - heat equation - interval - Neumann BC" begin
    # no known analytical solutions, so we compare 2 runs
    sols = []
    xs = []
    d = 5
    for _ in 1:2
        tspan = (0f0, 5f-1)
        dt = 5f-1  # time step
        μ(x, p, t) = 0f0 # advection coefficients
        σ(x, p, t) = 1f-1 # diffusion coefficients

        u_anal(x, t) = sum(x.^2) + d * σ(0., 0., 0.)^2 * t
        g(x) = sum(x.^2, dims=1)

        # d = 10
        u_domain = (fill(-5f-1, d), fill(5f-1, d))
        hls = d + 10 #hidden layer size

        nn = Flux.Chain(Dense(d,hls,relu),
                        Dense(hls,hls,relu),
                        Dense(hls,1)) # Neural network used by the scheme

        opt = ADAM(0.01) #optimiser
        alg = DeepSplitting(nn, opt = opt)

        f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = 0f0 .* v_y #TODO: this fix is not nice

        # defining the problem
        prob = PIDEProblem(g, f, μ, σ, tspan,
                            u_domain = u_domain,
                            neumann = u_domain
                            )
        # solving
        xs,ts,sol = solve(prob, alg, dt,
                        verbose = false,
                        use_cuda = use_cuda,
                        maxiters = 1000,
                        batch_size = 20000,
                        )
        push!(sols, sol[end])
    end
    u1 = [sols[1](x)[] for x in xs]
    u1_anal = [sols[2](x)[] for x in xs]
    e_l2 = mean(rel_error_l2.(u1, u1_anal))
    println("rel_error_l2 = ", e_l2, "\n")
    @test e_l2 < 0.1

    # for fun: plotting
    # close("all")
    # fig, ax = plt.subplots()
    # xgrid = hcat([x[1:2] for x in xs]...)
    # #Deepsplitting sol
    # ax.scatter(xgrid[1,:], xgrid[2,:], c = reduce(vcat, sols[end].(xs))[:], s=.2)
    # ax.set_title("DeepSplitting")
    # gcf()
end

@testset "DeepSplitting - exponential - interval" begin
    for d in [1, 3, 10]
        println("test for d = ", d)
        tspan = (0f0, 5f-1)
        dt = 5f-2  # time step
        μ(x, p, t) = 0f0 # advection coefficients
        σ(x, p, t) = 0f0 # diffusion coefficients
        r = 1f-1

        u_anal(x, t) = (g(x) * exp(r * t))[]
        g(x) = sum(x.^2, dims=1) .+ 2f0

        # d = 10
        u_domain = (fill(-5f-1, d), fill(5f-1, d))
        hls = d + 50 #hidden layer size

        nn = Flux.Chain(Dense(d,hls,relu),
                        Dense(hls,hls,relu),
                        Dense(hls,1)) # Neural network used by the scheme

        opt = ADAM(0.01) #optimiser
        alg = DeepSplitting(nn, opt = opt)

        f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = r * v_y #TODO: this fix is not nice

        # defining the problem
        prob = PIDEProblem(g, f, μ, σ, tspan,
                            u_domain = u_domain,
                            )
        # solving
        xs,ts,sol = solve(prob, alg, dt,
                        verbose = false,
                        use_cuda = use_cuda,
                        maxiters = 1000,
                        batch_size = 100)
        u1 = [sol[end](x)[] for x in xs]
        u1_anal = u_anal.(xs, tspan[end])
        e_l2 = mean(rel_error_l2.(u1, u1_anal))
        println("rel_error_l2 = ", e_l2, "\n")
        @test e_l2 < 0.1
    end
end

@testset "DeepSplitting algorithm - allen cahn" begin
    batch_size = 2000
    train_steps = 500

    tspan = (3f-1, 6f-1)
    dt = 2f-2  # time step

    μ(x, p, t) = 0f0 # advection coefficients
    σ(x, p, t) = 1f0 # diffusion coefficients

    for d in [10]

        hls = d + 50 #hidden layer size

        nn = Flux.Chain(Dense(d,hls,tanh),
                        Dense(hls,hls,tanh),
                        Dense(hls,1)) # Neural network used by the scheme

        opt = ADAM(1e-3) #optimiser
        alg = DeepSplitting(nn, opt = opt )

        X0 = fill(0f0,d)  # initial point
        g(X) =  1f0 ./ (2f0 .+ 4f-1 * sum(X.^2, dims=1))   # initial condition
        a(u) = u - u^3
        f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = - a.(v_y) # nonlocal nonlinear part of the

        # defining the problem
        prob = PIDEProblem(g, f, μ, σ, tspan, x = X0 )
        # solving
        @time xs,ts,sol = solve(prob,
                        alg,
                        dt,
                        verbose = false,
                        # abstol = 1e-5,
                        use_cuda = use_cuda,
                        maxiters = train_steps,
                        batch_size=batch_size)
        u1 = sol[end]
        # value coming from \cite{Beck2017a}
        e_l2 = rel_error_l2(u1, 0.30879)
        @test e_l2 < 0.5 # this is quite high as a relative error.
        println("d = $d, rel_error_l2 = $e_l2")
    end
end

@testset "DeepSplitting - allen cahn local-  Neumann BC" begin
    batch_size = 1000
    train_steps = 1000

    tspan = (0f0, 5f-1)
    dt = 5f-2  # time step

    μ(x, p, t) = 0f0 # advection coefficients
    σ(x, p, t) = 1f0 # diffusion coefficients

    for d in [1,2,5]
        u_domain = (fill(-5f-1, d), fill(5f-1, d))

        hls = d + 50 #hidden layer size

        nn = Flux.Chain(Dense(d,hls,tanh),
                Dense(hls,hls,tanh),
                Dense(hls,1)) # Neural network used by the scheme

        opt = ADAM(1e-2) #optimiser
        alg = DeepSplitting(nn, opt = opt )

        X0 = fill(0f0,d)  # initial point
        g(X) = exp.(-0.25f0 * sum(X.^2,dims=1))   # initial condition
        a(u) = u - u^3
        f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = a.(v_y) # nonlocal nonlinear part of the

        # defining the problem
        prob = PIDEProblem(g, f, μ, σ, tspan, x = X0, neumann = u_domain )
        # solving
        @time xs,ts,sol = solve(prob,
                        alg,
                        dt,
                        verbose = false,
                        abstol = 1e-5,
                        use_cuda = use_cuda,
                        maxiters = train_steps,
                        batch_size=batch_size)
        u1 = sol[end]
        @test !isnan(u1)
        println("d = $d, u1 = $(u1)")
    end
end

# TODO: not working because of a bug in Zygote differentiation rules for adjoints
if false
    @testset "DeepSplitting - Black-Scholes-Barenblatt equation" begin
        batch_size = 20
        train_steps = 500

        tspan = (0f0, 1f0)
        dt = 25f-2  # time step

        r = 5f-2
        sigma = 4f-1

        μ(x, p, t) = 0f0 # advection coefficients
        σ(x, p, t) =  sigma * X.^2 # diffusion coefficients

        for d in [30]

            hls = d + 50 #hidden layer size

            nn = Flux.Chain(Dense(d,hls,tanh),
                            Dense(hls,hls,tanh),
                            Dense(hls,1)) # Neural network used by the scheme

            opt = ADAM(1e-3) #optimiser
            alg = DeepSplitting(nn, opt = opt )

            X0 = repeat([1.0f0, 0.5f0], div(d,2))  # initial point
            g(X) =  sum(X.^2, dims=1) # initial condition
            f(y,z,v_y,v_z,∇v_y,∇v_z, p, t) = r * (v_y .- sum(y .* ∇v_y, dims=1))

            # defining the problem
            prob = PIDEProblem(g, f, μ, σ, tspan, x = X0 )
            # solving
            @time xs,ts,sol = solve(prob,
                            alg,
                            dt,
                            verbose = false,
                            # abstol = 1e-5,
                            use_cuda = use_cuda,
                            maxiters = train_steps,
                            batch_size=batch_size)
            u1 = sol[end]
            @test !isnan(u1)
            println("d = $d, u1 = $(u1)")
        end
    end

    @testset "DeepSplitting - Hamilton Jacobi Bellman Equation" begin

        batch_size = 20
        train_steps = 500

        tspan = (0f0, 1f0)
        dt = 2f-1  # time step

        λ = 1.0f0

        μ(x, p, t) = 0f0 # advection coefficients
        σ(x, p, t) =  sqrt(2f0) # diffusion coefficients

        d = 20

        T = tspan[2]
        MC = 10^5
        W() = randn(d,1)
        u_analytical(x, t) = -(1/λ)*log(mean(exp(-λ*g(x .+ sqrt(2f0)*abs.(T-t).*W())) for _ = 1:MC))

        hls = d + 50 #hidden layer size

        nn = Flux.Chain(Dense(d,hls,tanh),
                        Dense(hls,hls,tanh),
                        Dense(hls,1)) # Neural network used by the scheme

        opt = ADAM(1e-3) #optimiser
        alg = DeepSplitting(nn, opt = opt )

        X0 = fill(0.0f0,d)  # initial point
        g(X) =  log.(5f-1 .+ 5f-1 * sum(X.^2, dims=1)) # initial condition
        f(y,z,v_y,v_z,∇v_y,∇v_z, p, t) = λ * sum(∇v_y.^2, dims=1)

        # defining the problem
        prob = PIDEProblem(g, f, μ, σ, tspan, x = X0 )
        # solving
        @time xs,ts,sol = solve(prob,
                        alg,
                        dt,
                        verbose = false,
                        # abstol = 1e-5,
                        use_cuda = use_cuda,
                        maxiters = train_steps,
                        batch_size=batch_size)

        u1 = sol[end]

        analytical_ans = u_analytical(x0, tspan[1])
        error_l2 = rel_error_l2(u1, analytical_ans)

        @test error_l2 < 2f0

        println("d = $d, u1 = $(u1)")

    end
end

@testset "DeepSplitting - Black-Scholes Equation with Default Risk" begin

    batch_size = 1000
    train_steps = 5000

    tspan = (0f0, 1f0)
    dt = 1.25f-1  # time step

    μ(x, p, t) = 0.02f0 * x # advection coefficients
    σ(x, p, t) = 0.2f0 * x # diffusion coefficients

    d = 20

    hls = d + 50 #hidden layer size

    nn = Flux.Chain(Dense(d,hls,tanh),
                    Dense(hls,hls,tanh),
                    Dense(hls,1)) # Neural network used by the scheme

    opt = ADAM(1e-3) #optimiser
    alg = DeepSplitting(nn, opt = opt )

    X0 = fill(100f0,d)  # initial point
    g(X) =  minimum(X, dims=1) # initial condition

    δ = 2.0f0/3f0
    R = 0.02f0

    vh = 50.0f0
    vl = 70.0f0
    γh = 0.2f0
    γl = 0.02f0

    Q(u) = (u .< vh) .* γh .+ (u .>= vl) .*  γl .+ ( (u .>= vh) .* (u .< vl)) .* (((γh - γl) / (vh - vl)) * (u .- vh) .+ γh)

    µc = 0.02f0
    σc = 0.2f0

    f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = -(1f0 - δ) * Q.(v_y) .* v_y .- R * v_y

    # defining the problem
    prob = PIDEProblem(g, f, μ, σ, tspan, x = X0 )
    # solving
    @time xs,ts,sol = solve(prob,
                    alg,
                    dt,
                    # verbose = true,
                    # abstol = 1e-5,
                    use_cuda = use_cuda,
                    maxiters = train_steps,
                    batch_size=batch_size)

    u1 = sol[end]
    analytical_ans = 60.781
    error_l2 = rel_error_l2(u1, analytical_ans)
    @test error_l2 < 0.1
    println("d = $d, error_l2 = $(error_l2)")
end

###################################################
########### NON LOCAL #############################
###################################################

@testset "DeepSplitting - Hamel example - udomain" begin
    tspan = (0f0,15f-2)
    dt = 5f-2 # time step
    μ(x, p, t) = 0f0 # advection coefficients
    σ(x, p, t) = 1f-1 #1f-1 # diffusion coefficients
    ss0 = 5f-2#std g0

    # Analytic sol
    function _SS(x, t, p)
        d = length(x)
        MM = σ(x, p, t) * ones(d)
        SSt = MM .* ((MM .* sinh.(MM *t) .+ ss0 .* cosh.( MM * t)) ./ (MM .* cosh.(MM * t ) .+ ss0 .* sinh.(MM * t)))
        return SSt
    end

    function uanal(x, t, p)
            d = length(x)
            return (2*π)^(-d/2) * prod(_SS(x, t, p) .^(-1/2)) * exp(-0.5 *sum(x .^2 ./ _SS(x, t, p)) )
    end

    sols = []
    xs = []
    for d in [5]
        U = 25f-2
        u_domain = (fill(-U, d), fill(U, d))

        batch_size = 10000
        train_steps = 2000
        K = 1

        hls = d + 50 #hidden layer size

        nn_batch = Flux.Chain(
                            # BatchNorm(d, affine = true, dim = 1),
                            Dense(d, hls, tanh),
                            # BatchNorm(hls, affine = true, dim = 1),
                            Dense(hls, hls, tanh),
                            # BatchNorm(hls, affine = true, dim = 1),
                            Dense(hls, 1, x->x^2)) # positive function

        opt = ADAM(1e-2)#optimiser
        alg = DeepSplitting(nn_batch, K=K, opt = opt, mc_sample = UniformSampling(u_domain[1], u_domain[2]) )

        g(x) = Float32((2*π)^(-d/2)) * ss0^(- Float32(d) * 5f-1) * exp.(-5f-1 *sum(x .^2f0 / ss0, dims = 1)) # initial condition
        m(x) = - 5f-1 * sum(x.^2, dims=1)
        vol = prod(u_domain[2] - u_domain[1])
        f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) =  max.(v_y, 0f0) .* (m(y) .- vol *  max.(v_z, 0f0) .* m(z)) # nonlocal nonlinear part of the

        # defining the problem
        prob = PIDEProblem(g, f, μ, σ, tspan,
                            u_domain = u_domain
                            )
        # solving
        xgrid,ts,sol = solve(prob,
                        alg,
                        dt,
                        verbose = false,
                        abstol = 1f-3,
                        maxiters = train_steps,
                        batch_size = batch_size,
                        use_cuda = use_cuda
                        )
        u1 = [sol[end](x)[] for x in xgrid]
        u1_anal = uanal.(xgrid, tspan[end], nothing)
        e_l2 = mean(rel_error_l2.(u1, u1_anal))
        println("rel_error_l2 = ", e_l2, "\n")
        @test e_l2 < 0.1
        push!(sols, sol[end])
    end
end

@testset "DeepSplitting algorithm - allen cahn non local - Neumann BC" begin
    batch_size = 2000
    train_steps = 1000
    K = 1
    tspan = (0f0, 5f-1)
    dt = 5f-2  # time step

    μ(x, p, t) = 0f0 # advection coefficients
    σ(x, p, t) = 1f-1 #1f-1 # diffusion coefficients

    for d in [1,2,5]
        u1s = []
        for _ in 1:2
            u_domain = (fill(-5f-1, d), fill(5f-1, d))

            hls = d + 50 #hidden layer size

            nn = Flux.Chain(Dense(d,hls,tanh),
                            Dense(hls,hls,tanh),
                            Dense(hls,1)) # Neural network used by the scheme

            opt = ADAM(1e-2) #optimiser
            alg = DeepSplitting(nn, K=K, opt = opt, mc_sample = UniformSampling(u_domain[1],u_domain[2]) )

            x = fill(0f0,d)  # initial point
            g(X) = exp.(-0.25f0 * sum(X.^2,dims=1))   # initial condition
            a(u) = u - u^3
            f(y,z,v_y,v_z,∇v_y,∇v_z, p, t) = a.(v_y) .- a.(v_z) #.* Float32(π^(d/2)) * σ_sampling^d .* exp.(sum(z.^2, dims = 1) / σ_sampling^2) # nonlocal nonlinear part of the

            # defining the problem
            prob = PIDEProblem(g, f, μ, σ, tspan, x = x, neumann = u_domain)
            # solving
            @time xs,ts,sol = solve(prob,
                            alg,
                            dt,
                            # verbose = true,
                            # abstol=1e-5,
                            use_cuda = use_cuda,
                            maxiters = train_steps,
                            batch_size=batch_size)
            push!(u1s, sol[end])
            println("d = $d, u1 = $(sol[end])")

        end
        e_l2 = mean(rel_error_l2.(u1s[1], u1s[2]))
        println("rel_error_l2 = ", e_l2, "\n")
        @test e_l2 < 0.1
    end
end


text:
"""
Initialise all the storage arrays that will be used in a GRAPE optimisation
"""
function init_GRAPE(X, n_timeslices, n_ensemble, A, n_controls)
    # states
    states = [similar(X) .* 0.0 for i = 1:n_timeslices+1, j = 1:n_ensemble]
    # costates
    costates = [similar(X) .* 0.0 for i = 1:n_timeslices+1, j = 1:n_ensemble]
    # exp(gens)
    propagators = [similar(A) .* 0.0 for i = 1:n_timeslices, j = 1:n_ensemble]

    fom = 0.0
    gradient = zeros(n_ensemble, n_controls, n_timeslices)
    # end
    return (states, costates, propagators, fom, gradient)
end

"""
Compute the eigen decomposition of a generator matrix G. Based on work by both Shai Machnes and Ville Bergholm.

Types are an issue here.
"""
function eig_factors(G; antihermitian = false, tol = 1e-10)
    if antihermitian
        d, v = eigen(1im * G)
        d = -1im .* d
    else
        d, v = eigen(G)
    end

    ooo = ones(size(d)[1])
    temp = d * ooo'
    diff = temp - temp'
    degenerate_mask = abs.(diff) .< tol
    diff[degenerate_mask] .= 1
    exp_d = exp.(d)

    temp = exp_d * ooo'
    exp_diff = temp - temp'
    zeta = exp_diff / diff
    zeta[degenerate_mask] .= temp[degenerate_mask]
    return (v, zeta, exp_d)
end

"""
Function to compute the matrix exponential using the eigenvalue decomposition of the matrix, this is slower than exp if your matrix can be a StaticArray
"""
function expm_exact_gradient(H::T, dt)::T where {T}
    dt_H = dt * H
    v, zeta, exp_d = eig_factors(dt_H, antihermitian = true)

    v * diagm(exp_d) * v'
end

"""
Copied from Ville's code, seems to be slower in Julia than simply taking trace but is faster in Matlab/Python
Function to compute trace(A @ B) efficiently.

Utilizes the identity trace(A @ B) == sum(transpose(A) * B).
Left side is O(n^3) to compute, right side is O(n^2).
"""
function trace_matmul(A, B)
    sum(transpose(A) .* B)
end


text:
function forward(backend::GPUBackend, state::MultinomialLogisticLossLayerState, inputs::Vector{Blob})
  pred      = inputs[1]
  label     = inputs[2]
  data_type = eltype(pred)

  spatial_dim, channels, num = split_dims(pred, state.op_dim)
  prob_dim = channels

  x_block = round(Int, ceil(convert(Float64, num)/CUDA.THREADS_PER_BLOCK_X))
  y_block = spatial_dim

  loss_blob = make_zero_blob(backend, Float32, 1, 1, 1, 1)

  if data_type == Float32
    kernel = backend.mocha.logistic_loss_forward_float
  elseif data_type == Float64
    kernel = backend.mocha.logistic_loss_forward_double
  else
    error("Unsupported data type $data_type")
  end

  if isa(state.weights_blob, NullBlob)
    weights = convert(Ptr{data_type}, 0)
  else
    weights = state.weights_blob.ptr.p
  end

  CUDA.launch(kernel, (x_block, y_block), (CUDA.THREADS_PER_BLOCK_X, 1),
      (pred.ptr.p, label.ptr.p, weights, num, spatial_dim, prob_dim, loss_blob.ptr.p))

  loss = Float32[0]
  copy!(loss, loss_blob)
  state.loss = state.layer.weight * loss[1] / (spatial_dim * num)
  destroy(loss_blob)
end


text:
export linear,∂linear

function linear(rᵢ)
  rₒ=rᵢ
end
function ∂linear(rᵢ)
  rₒ=constant(Float32(1),size(rᵢ))
end


text:
include("equations/index.jl")
include("reactors/index.jl")
include("scans/index.jl")


text:
"Holds the tableau of an Variational Partitioned Additive Runge-Kutta method."
const TableauVPARK = AbstractTableauSPARK{:vpark}

"Parameters for right-hand side function of Variational Partitioned Additive Runge-Kutta methods."
const ParametersVPARK = AbstractParametersSPARK{:vpark}

@doc raw"""
Variational partitioned additive Runge-Kutta integrator.

This integrator solves the following system of equations for the internal stages,
```math
\begin{aligned}
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} U_{n,j} , & i &= 1, ..., s , \\
P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} a_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} G_{n,j} , & i &= 1, ..., s , \\
\tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} U_{n,j} , & i &= 1, ..., r , \\
\tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} G_{n,j} , & i &= 1, ..., r , \\
\tilde{\Phi}_{n,i} &= 0 , & i &= 1, ..., r ,
\end{aligned}
```
with definitions
```math
\begin{aligned}
P_{n,i} &= \frac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\
F_{n,i} &= \frac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\
U_{n,i} &= \hphantom{-} \frac{\partial \phi}{\partial p} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\
G_{n,i} &=           -  \frac{\partial \phi}{\partial q} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\
\tilde{\Phi}_{n,i} &= \phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) , & i &= 1, ..., r ,
\end{aligned}
```
and update rule
```math
\begin{aligned}
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} V_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} U_{n,i} , \\
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b_{i} F_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} G_{n,i} .
\end{aligned}
```
"""
struct IntegratorVPARK{DT, TT, D, S, R, PT <: ParametersVPARK{DT,TT,D,S,R},
                                        ST <: NonlinearSolver{DT},
                                        IT <: InitialGuessIODE{TT}} <: AbstractIntegratorVSPARK{DT,TT,D,S,R}
    params::PT
    solver::ST
    iguess::IT
    caches::CacheDict{PT}

    function IntegratorVPARK(params::ParametersVPARK{DT,TT,D,S,R}, solver::ST, iguess::IT, caches) where {DT,TT,D,S,R,ST,IT}
        new{DT, TT, D, S, R, typeof(params), ST, IT}(params, solver, iguess, caches)
    end

    function IntegratorVPARK{DT,D}(equations::NamedTuple, tableau::TableauVPARK{TT}, Δt::TT) where {DT,TT,D}
        # get number of stages
        S = tableau.s
        R = tableau.r

        N = 3*D*S + 3*D*R

        if isdefined(tableau, :d) && length(tableau.d) > 0
            N += D
        end

        # create params
        params = ParametersVPARK{DT,D}(equations, tableau, Δt)

        # create cache dict
        caches = CacheDict(params)

        # create solver
        solver = create_nonlinear_solver(DT, N, params, caches)

        # create initial guess
        iguess = InitialGuessIODE(get_config(:ig_extrapolation), equations[:v̄], equations[:f̄], Δt)

        # create integrator
        IntegratorVPARK(params, solver, iguess, caches)
    end

    function IntegratorVPARK(equation::Union{IDAE{DT}, LDAE{DT}}, tableau::TableauVPARK, Δt; kwargs...) where {DT}
        IntegratorVPARK{DT, ndims(equation)}(get_functions(equation), tableau, Δt; kwargs...)
    end
end

GeometricBase.nconstraints(::IntegratorVPARK{DT,TT,D}) where {DT,TT,D} = D

function compute_stages!(x::Vector{ST}, cache::IntegratorCacheSPARK{ST,D,S,R},
                                        params::ParametersVPARK{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    local tpᵢ::TT
    local tλᵢ::TT

    for i in 1:S
        for k in 1:D
            # copy x to Y, Z
            cache.Yi[i][k] = x[3*(D*(i-1)+k-1)+1]
            cache.Zi[i][k] = x[3*(D*(i-1)+k-1)+2]
            cache.Vi[i][k] = x[3*(D*(i-1)+k-1)+3]

            # compute Q and P
            cache.Qi[i][k] = params.q[k] + params.Δt * cache.Yi[i][k]
            cache.Pi[i][k] = params.p[k] + params.Δt * cache.Zi[i][k]
        end

        # compute f(X)
        tpᵢ = params.t + params.Δt * params.tab.p.c[i]
        params.equs[:ϑ](tpᵢ, cache.Qi[i], cache.Vi[i], cache.Φi[i])
        params.equs[:f](tpᵢ, cache.Qi[i], cache.Vi[i], cache.Fi[i])

        cache.Φi[i] .-= cache.Pi[i]
    end

    for i in 1:R
        for k in 1:D
            # copy y to Y, Z and Λ
            cache.Yp[i][k] = x[3*D*S+3*(D*(i-1)+k-1)+1]
            cache.Zp[i][k] = x[3*D*S+3*(D*(i-1)+k-1)+2]
            cache.Λp[i][k] = x[3*D*S+3*(D*(i-1)+k-1)+3]

            # compute Q and V
            cache.Qp[i][k] = params.q[k] + params.Δt * cache.Yp[i][k]
            cache.Pp[i][k] = params.p[k] + params.Δt * cache.Zp[i][k]
        end

        # compute f(X)
        tλᵢ = params.t + params.Δt * params.tab.λ.c[i]
        params.equs[:u](tλᵢ, cache.Qp[i], cache.Pp[i], cache.Λp[i], cache.Up[i])
        params.equs[:g](tλᵢ, cache.Qp[i], cache.Pp[i], cache.Λp[i], cache.Gp[i])
        params.equs[:ϕ](tλᵢ, cache.Qp[i], cache.Pp[i], cache.Φp[i])
    end

    if isdefined(params.tab, :d) && length(params.tab.d) > 0
        for k in 1:D
            cache.μ[k] = x[3*D*S+3*D*R+k]
        end
    end
end

"Compute stages of variational partitioned additive Runge-Kutta methods."
function Integrators.function_stages!(y::Vector{ST}, b::Vector{ST}, params::ParametersVPARK{DT,TT,D,S,R},
                                      caches::CacheDict) where {ST,DT,TT,D,S,R}

    # get cache for internal stages
    cache = caches[ST]

    compute_stages!(y, cache, params)

    # compute b = - [(Y-AV-AU), (Z-AF-AG), Φ]
    for i in 1:S
        for k in 1:D
            b[3*(D*(i-1)+k-1)+1] = - cache.Yi[i][k]
            b[3*(D*(i-1)+k-1)+2] = - cache.Zi[i][k]
            b[3*(D*(i-1)+k-1)+3] = - cache.Φi[i][k]
            for j in 1:S
                b[3*(D*(i-1)+k-1)+1] += params.tab.q.a[i,j] * cache.Vi[j][k]
                b[3*(D*(i-1)+k-1)+2] += params.tab.p.a[i,j] * cache.Fi[j][k]
            end
            for j in 1:R
                b[3*(D*(i-1)+k-1)+1] += params.tab.q.α[i,j] * cache.Up[j][k]
                b[3*(D*(i-1)+k-1)+2] += params.tab.p.α[i,j] * cache.Gp[j][k]
            end
        end
    end

    # compute b = - [(Y-AV-AU), (Z-AF-AG), Φ]
    for i in 1:R
        for k in 1:D
            b[3*D*S+3*(D*(i-1)+k-1)+1] = - cache.Yp[i][k]
            b[3*D*S+3*(D*(i-1)+k-1)+2] = - cache.Zp[i][k]
            b[3*D*S+3*(D*(i-1)+k-1)+3] = - cache.Φp[i][k]
            for j in 1:S
                b[3*D*S+3*(D*(i-1)+k-1)+1] += params.tab.q̃.a[i,j] * cache.Vi[j][k]
                b[3*D*S+3*(D*(i-1)+k-1)+2] += params.tab.p̃.a[i,j] * cache.Fi[j][k]
            end
            for j in 1:R
                b[3*D*S+3*(D*(i-1)+k-1)+1] += params.tab.q̃.α[i,j] * cache.Up[j][k]
                b[3*D*S+3*(D*(i-1)+k-1)+2] += params.tab.p̃.α[i,j] * cache.Gp[j][k]
            end
        end
    end

    # compute b = - [Λ₁-λ]
    if params.tab.λ.c[1] == 0
        for k in 1:D
            b[3*D*S+3*(k-1)+3] = - cache.Λp[1][k] + params.λ[k]
        end
    end

    if isdefined(params.tab, :d) && length(params.tab.d) > 0
        for i in 1:S
            for k in 1:D
                b[3*(D*(i-1)+k-1)+3] -= cache.μ[k] * params.tab.d[i] / params.tab.p.b[i]
            end
        end

        for k in 1:D
            b[3*D*S+3*D*R+k] = 0
            for i in 1:S
                b[3*D*S+3*D*R+k] -= cache.Vi[i][k] * params.tab.d[i]
            end
        end
    end
end

function update_solution!(int::IntegratorVPARK{DT,TT}, sol::AtomicSolutionPDAE{DT,TT},
                          cache::IntegratorCacheSPARK{DT}=int.caches[DT]) where {DT,TT}
    # compute final update
    update_solution!(sol.q, sol.q̃, cache.Vi, int.params.tab.q.b, timestep(int))
    update_solution!(sol.p, sol.p̃, cache.Fi, int.params.tab.p.b, timestep(int))

    # compute projection
    update_solution!(sol.q, sol.q̃, cache.Up, int.params.tab.q.β, timestep(int))
    update_solution!(sol.p, sol.p̃, cache.Gp, int.params.tab.p.β, timestep(int))
    update_multiplier!(sol.λ, cache.Λp, int.params.tab.λ.b)
end


text:
function sortfile(infname::AbstractString, outfname::AbstractString, offset::Int=0, len::Int=0)
    dfs = HDFSClient(TeraSort.DFSHOST, TeraSort.DFSPORT, UserGroupInformation(TeraSort.DFSUSER))
    exists(dfs, infname) || (return 0)
    fptr = open(dfs, infname, "r")
    wptr = open(dfs, outfname, "w")
    try
        (len == 0) && (len = Int(filesize(fptr)))
        sortfile(fptr, wptr, 0, len)
    finally
        close(fptr)
        close(wptr)
    end
    len
end

function sortfile(fptr::IO, wptr::IO, offset::Int, len::Int)
    seek(fptr, offset)
    data = Array(UInt8, len)
    read!(fptr, data)
    datastr = bytestring(data)
    dataarr = split(datastr)
    sort!(dataarr)
    for s in dataarr
        write(wptr, s, '\n')
    end
    nothing
end

function sortrandfile(fname::AbstractString)
    println("sortrandfile starting at $(myid())")
    locfname = "$(fname)-$(myid())"
    d,f = splitdir(locfname)
    sortedfname = joinpath(TeraSort.OUTDIR, f)
    nbytes = sortfile(locfname, sortedfname)
    println("sortrandfile sorted $nbytes bytes at $(myid())")
    nbytes
end

function distsortfile()
    wrkrs = workers()
    nwrkrs = length(wrkrs)
    refs = Array(RemoteRef, nwrkrs)
    idx = 1
    fname = joinpath(INDIR, "randstr")
    for nid in wrkrs
        refs[idx] = @spawnat nid TeraSort.sortrandfile(fname)
        idx += 1
    end
    nbytes = 0
    for idx in 1:nwrkrs
        nbytes += fetch(refs[idx])
    end
    nbytes
end

function terasort()
    nbytes, t, b, g = @timed distsortfile()
    println("sorted $nbytes bytes")
    println("time:$t, memused:$b, gctime:$g")
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dashiconpicker

"""
    dashiconpicker(;kwargs...)

A DashIconPicker component.

Keyword arguments:
- `id` (String; required): The ID used to identify this component in Dash callbacks.
- `value` (String; optional): The value displayed in the input.
- `hideSearch` (Bool; optional): If true, the search input is not displayed. Default is false.
- `containerStyles` (Dict with Strings as keys and values of type String; optional): Styles for the picker container
- `buttonStyles` (Dict with Strings as keys and values of type String; optional): Styles for the picker button
- `buttonIconStyles` (Dict with Strings as keys and values of type String; optional): Styles for the picker button icon
- `pickerIconStyles` (Dict with Strings as keys and values of type String; optional): Styles for the icons inside of the picker
- `searchInputStyles` (Dict with Strings as keys and values of type String; optional): Styles for the search input inside of the picker
- `className` (String; optional): CSS class for div
"""
function dashiconpicker(; kwargs...)
        available_props = Symbol[:id, :value, :hideSearch, :containerStyles, :buttonStyles, :buttonIconStyles, :pickerIconStyles, :searchInputStyles, :className]
        wild_props = Symbol[]
        return Component("dashiconpicker", "DashIconPicker", "dash_react_font_awesome_icon_picker", available_props, wild_props; kwargs...)
end


text:
mutable struct RECMarket # Assumes completely inelastic demand

    # Parameters
    demand::Float64   # annual renewable energy (MWh)
    pricecap::Float64 # alternate compliance cost ($/MWh)

    # Variables
    shortfall::VariableRef # RPS shortfall variable (MWh)

    # Expressions
    shortfallcost::ExpressionRef # Compliance penalty costs ($)

    # Constraints
    minshortfall::GreaterThanConstraintRef # Minimum RE shortfall
    marketclearing::GreaterThanConstraintRef # REC requirement

    function RECMarket(demand::Float64, pricecap::Float64)

        @assert demand >= 0
        @assert pricecap >= 0

        new(demand, pricecap)

    end

end

function RECMarket(marketpath::String)

    rulesdata = DataFrame!(CSV.File(joinpath(marketpath, "rules.csv"),
                                    types=scenarios_market_param_types))

    demand = first(rulesdata.demand)
    pricecap = first(rulesdata.pricecap)

    return RECMarket(demand, pricecap)

end

function setup!(
    market::RECMarket, m::Model,
    ops::Operations{R,G1,G2,G3,I,T,P}, periodweights::Vector{Float64},
    s::AbstractScenario
) where {R,G1,G2,G3,I,T,P}

    # Variables

    market.shortfall = @variable(m)
    set_name(market.shortfall, "shortfall_rec_$(s.name)")

    # Expressions

    market.shortfallcost = @expression(m, market.shortfall * market.pricecap)

    # Constraints

    market.minshortfall = @constraint(m, market.shortfall >= 0)

    market.marketclearing = @constraint(m,
        sum(periodweights[p] * sum(ops.variable.recenergy[:,:,p])
            for p in 1:P) + market.shortfall >= market.demand)

    return

end

welfare(x::RECMarket) = -x.shortfallcost


text:
#=
  https://www.youtube.com/watch?v=e1Ykk_CqKTY&t=458s

  Probabilistic Programming: What It Is and How It Works - Noel Welsh

  We can see either 1, 2, or 3 cats.
  There are 3 different enticements:

  - Milkshake
  - Fish
  - Nothing

  And there are different probabilities how many cats there are given
  an enticement, see below.

  Now: We see 3 cats, what is the probability that it's a milkshake?

  The video got the following (for 3 cats):
   - milkshake: 0.42
   - fish: 0.04
   - nothing: 0.03

  Normalized to percentage (from the video):

  0.42/(0.42 + 0.04 + 0.03) milkshake
      0.85714285714285714286
  0.04/(0.42 + 0.04 + 0.03) fish
      0.081632653061224489796
  0.03/(0.42 + 0.04 + 0.03)  nothing
      0.061224489795918367347

  Here are two models.
  Model2 seems to be the correct one (i.e. corresponds to the video),
  i.e. where the enticements are mutual exclusive.
  In model1 we accept that fish and milkshake can coexist.

  See ~/cplint/cat.pl
      ~/blog/cat.blog
      ~/psi/cat.psi
=#

using Turing # , StatsPlots, DataFrames
include("jl_utils.jl")
# using ReverseDiff, Zygote, Tracker
# Turing.setadbackend(:reversediff)
# Turing.setadbackend(:zygote)
# Turing.setadbackend(:tracker)

# This is not correct
@model function cat1()
    NoCat = 0
    OneCat = 1
    TwoCats = 2
    ThreeCats = 3

    # probability of an enticement in the garden
    milkshake      ~ Bernoulli(0.6)
    fish           ~ Bernoulli(0.1)
    no_enticement  ~ Bernoulli(0.3)

    # Number of cats per enticement
    vs = [NoCat,OneCat,TwoCats,ThreeCats]
    vs_cat ~ DiscreteUniform(NoCat,ThreeCats)
    if !no_enticement && milkshake
        cat ~ Categorical(simplex([0.0,0.1,0.2,0.7]))
    elseif !no_enticement && fish
        cat ~ Categorical(simplex([0.0,0.2,0.3,0.4]))
    elseif no_enticement && !fish && !milkshake
        cat ~ Categorical(simplex([0.0,0.6,0.3,0.1]))
    else
        cat ~ Categorical(simplex([1.0,0.0,0.0,0.0]))
    end
    true ~ Dirac(vs_cat == vs[cat])
    true ~ Dirac(vs_cat == ThreeCats)

end

#=
println("Model 1:")
model = cat1()
num_chns = 4
# chns = sample(model, Prior(), MCMCThreads(), 1000, num_chns)
chns = sample(model, MH(), MCMCThreads(), 10_000, num_chns)
# chns = sample(model, MH(), 1000)
display(chns)
=#

# This is correct (according to the video)
#=
   Summary Statistics
            parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
                Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

            enticement    1.2108    0.5429     0.0014    0.0093   2253.8077    1.0013       70.4381
  enticement_milkshake    0.8532    0.3539     0.0009    0.0061   2150.1791    1.0013       67.1994
       enticement_fish    0.0824    0.2750     0.0007    0.0046   2344.6207    1.0014       73.2763
    enticement_nothing    0.0642    0.2452     0.0006    0.0041   2630.1451    1.0013       82.1997
                  cats    2.9992    0.0403     0.0001    0.0008   2814.4694    1.0013       87.9604

   Distributions of variable enticement
   milkshake  =>  136546  (0.853413)
   fish       =>   13181  (0.082381)
   no_enticement =>   10273  (0.064206)

=#
enticements = []
enticements_dict = Dict(1=>"milkshake",2=>"fish",3=>"no_enticement")
@model cat2() = begin

    milkshake = 1
    fish = 2
    no_enticement = 3

    # Priors
    # probability of an enticement in the garden
    enticement ~ Categorical(simplex([0.6,0.1,0.3]))
    # Note: This seems a little too complex...
    enticement_milkshake ~ Bernoulli(0.5)
    enticement_fish      ~ Bernoulli(0.5)
    enticement_nothing   ~ Bernoulli(0.5)
    cats ~ DiscreteUniform(1,3)
    if enticement == milkshake
        cats ~ Categorical(simplex([0.1,0.2,0.7]))
        enticement_milkshake ~ Bernoulli(1.0)
        enticement_fish      ~ Bernoulli(0.0)
        enticement_nothing   ~ Bernoulli(0.0)
    elseif enticement == fish
        cats ~ Categorical(simplex([0.2,0.4,0.4]))
        enticement_milkshake ~ Bernoulli(0.0)
        enticement_fish      ~ Bernoulli(1.0)
        enticement_nothing   ~ Bernoulli(0.0)
    else
        cats ~ Categorical(simplex([0.6,0.3,0.1]))
        enticement_milkshake ~ Bernoulli(0.0)
        enticement_fish      ~ Bernoulli(0.0)
        enticement_nothing   ~ Bernoulli(1.0)
    end

    # We observe 3 cats
    true ~ Dirac(cats == 3)

    return enticement
end

println("\n\nModel 2:")
model = cat2()
num_chns = 4
# chns = sample(model, Prior(), MCMCThreads(), 1000, num_chns)
chns = sample(model, MH(), MCMCThreads(), 40_000, num_chns)
# chns = sample(model, MH(), 40_000)
# chns = sample(model, PG(20), MCMCThreads(), 1000, num_chns)
# chns = sample(model, SMC(), MCMCThreads(), 1000, num_chns)
# chns = sample(model, IS(), MCMCThreads(), 1000, num_chns)

display(chns)

show_var_dist_pct(chns, :enticement,["milkshake","fish","no_enticement"])

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# strings are defined with double quotes
# like variables, strings can contain any unicode character
s1 = "The quick brown fox jumps over the lazy dog α,β,γ"
println(s1)
#> The quick brown fox jumps over the lazy dog α,β,γ

# [println](http://julia.readthedocs.org/en/latest/stdlib/base/#Base.println) adds a new line to the end of output
# [print](http://julia.readthedocs.org/en/latest/stdlib/base/#Base.print) can be used if you dont want that:
print("this")
#> this
print(" and")
#> and
print(" that.\n")
#> that.

# chars are defined with single quotes
c1 = 'a'
println(c1)
#> a
# the ascii value of a char can be found with Int():
println(c1, " ascii value = ", Int(c1))
#> a ascii value = 97
println("Int('α') == ", Int('α'))
#> Int('α') == 945

# so be aware that
println(Int('1') == 1)
#> false

# strings can be converted to upper case or lower case:
s1_caps = uppercase(s1)
s1_lower = lowercase(s1)
println(s1_caps, "\n", s1_lower)
#> THE QUICK BROWN FOX JUMPS OVER THE LAZY DOG Α,Β,Γ
#> the quick brown fox jumps over the lazy dog α,β,γ

# sub strings can be indexed like arrays:
# ([show](http://julia.readthedocs.org/en/latest/stdlib/base/#Base.show) prints the raw value)
show(s1[11]); println()
#> 'b'

# or sub strings can be created:
show(s1[1:10]); println()
#> "The quick "

# end is used for the end of the array or string
show(s1[end-10:end]); println()
#> "dog α,β,γ"

# julia allows string [Interpolation](http://julia.readthedocs.org/en/latest/manual/strings/#interpolation):
a = "wolcome"
b = "julia"
println("$a to $b.")
#> wolcome to julia.

# this can extend to evaluate statements:
println("1 + 2 = $(1 + 2)")
#> 1 + 2 = 3

# strings can also be concatenated using the * operator
# using * instead of + isn't intuitive when you start with Julia,
# however [people think it makes more sense](https://groups.google.com/forum/#!msg/julia-users/nQg_d_n0t1Q/9PSt5aya5TsJ)
s2 = "this" * " and" * " that"
println(s2)
#> this and that

# as well as the string function
s3 = string("this", " and", " that")
println(s3)
#> this and that

max_stars_repo_name: inkydragon/JuliaByExample

max_issues_repo_name: inkydragon/JuliaByExample

max_forks_repo_name: inkydragon/JuliaByExample


text:
using DiffResults, ForwardDiff, Distributions, StatsBase, StatsFuns

export ScoreParticleContainer

struct ScoreParticleContainer <: ParticleContainer
    X # Particles
    w # weights
    α # Mean of the score vector
    ℓ # Estimate of the log likelihood
    S # Score vector
    dg # Preallocated GradientResult
end

Base.size(ps::ScoreParticleContainer) = (size(ps.X,2),)
Base.length(ps::ScoreParticleContainer) = prod(size(ps))
Base.getindex(ps::ScoreParticleContainer,i::Int) = ps.X[:,i]
Base.setindex!(ps::ScoreParticleContainer,v,i::Int) = ps.X[:,i] = v
Base.IndexStyle(ps::ScoreParticleContainer) = IndexLinear()

resample!(ps,m,y,t,θ,sample_function=resample_stratified) = resample!(ps.X,ps.w,ps.α,ps.w[:],m.ξ,y,t,θ,length(ps.w),sample_function)

function resample!(X,w,α,ξs,ξ,y,t,θ,N,sample_function)
    for i in 1:N
        ξs[i] += ξ(X[:,i],y,t,θ)
    end
    Σ = logsumexp(ξs)
    idx = sample_function(exp.(ξs .- Σ),N)

    X[:,:] .= X[:,idx]
    w[:] .= w[idx]
    α[:,:] .= α[:,idx]
    ξs[idx],Σ
end

propagate!(ps,ξ,Σ,m,y,t,θ,λ) = propagate!(ps.X,ps.w,ps.α,ps.dg,ps.ℓ,ps.S,ξ,Σ,m.f,m.g,m.q,m.dp,y,t,θ,length(ps.w),λ)

function propagate!(X,w,α,dg,ℓ,S,ξ,Σ,f,g,q,dp,y,t,θ,N,λ)
    for i in 1:N
        dq = q(X[:,i],y,t,θ)
        Xn = rand(dq)

        dg = dp(dg,X[:,i],Xn,y,t,θ)

        w[i] += DiffResults.value(dg)
        #w[i] *= pdf(f(X[:,i],t,θ),Xn)
        #w[i] *= pdf(g(Xn,t,θ),y)
        w[i] -= ξ[i]
        w[i] -= logpdf(dq,Xn)

        α[:,i] *= λ
        α[:,i] .+= (1-λ)*S[:] .+ DiffResults.gradient(dg)
        X[:,i] .= Xn[:]
    end
    ℓn = logsumexp(w)
    ℓ[t+1] = ℓ[t] + ℓn - log(N) + Σ
    w .-= ℓn
    S[:] .= (α*exp.(w))[:]
end

function step(ps::ScoreParticleContainer,m::ProposalStateSpaceModel,y,t,θ,λ=0.95;Nthreshold=0.5,sample_function=resample_stratified)

    # Resample

    #rs = Resampler(ps,m,y,t,θ,sample_function)

    ξs,Σ = resample!(ps,m,y,t,θ,sample_function)

    # Propagate

    #ps2 = ScoreParticleContainer([propose(t,x,w,ξ,α,ps.S,y,θ,m.f,m.g,m.q,m.dp,ps.dg,λ) for (x,w,ξ,α) in rs],ps.ℓ,ps.dg,sum(rs.ξ),λ)

    propagate!(ps,ξs,Σ,m,y,t,θ,λ)

    #normalize_weights!(ps2)

    ps
end

function score_filter(m::ProposalStateSpaceModel,Y,θ,N;λ=0.95,Nthreshold=0.5,sample_function=resample_stratified,verbose=false)

    T = length(Y)

    μ0 = m.μ0
    f = m.f
    g = m.g

    if verbose
        p = Progress(T,1)
    end

    ps = initialize_particles(m,θ,Y,N,λ)

    for t in 1:T
        ps = step(ps,m,Y[t],t,θ,λ,Nthreshold=Nthreshold,sample_function=sample_function)
        if verbose
            next!(p)
        end
    end

    ps
end

function initialize_particles(μ0,f,g,θ,Y,N,T,λ;save_filter=false)

    # We should probably think of μ0 as a distribution for p(x_0)
    # not one for p(x_1)
    # Initialization gives X0, w = 1/N
    #x0 = rand(μ0(θ),N)

    X0 = rand(μ0(θ),N)
    # This used to be: cat(rand(μ0(θ),N)...,dims=3)
    # I think this was to make it work for matrix-valued random variables

    w = log.(ones(N)/N)
    ℓ = zeros(T+1)
    #w = [pdf(g(X0[:,end,i],1,θ),y)/N for i in 1:N]
    #ℓ = sum(w)
    #w ./= ℓ

    dg = DiffResults.GradientResult(θ)

    α = Array{Float64}(undef,length(θ),N)

    for i in 1:N
        dg = ForwardDiff.gradient!(dg,θ->logpdf(μ0(θ),X0[:,i]),θ)

        α[:,i] = DiffResults.gradient(dg)
    end

    S = α*exp.(w)

    ScoreParticleContainer(X0,w,α,ℓ,S,dg)
end
initialize_particles(m::ProposalStateSpaceModel,θ,Y,N,λ) = initialize_particles(m.μ0,m.f,m.g,θ,Y,N,length(Y),λ)


text:
module AF

include("AFDefs.jl")
include("AFError.jl")

export
    ArrayFire,
    getBackend,
    getSupportedBackends,
    scope!,
    @scope,
    register!

abstract AFImpl

include("AFLib.jl")
include("Results.jl")

type ArrayFire{T<:Backend}
    ptr::Ptr{Void}
    device::AFImpl
    create::AFImpl
    binary::AFImpl
    unary::AFImpl
    index::AFImpl
    vectorAlgos::AFImpl
    modify::AFImpl
    reduction::AFImpl
    math::AFImpl
    linAlg::AFImpl
    moveAndReorder::AFImpl
    statistics::AFImpl

    createHandle::Ptr{Void}
    createArray::Ptr{Void}
    retainArray::Ptr{Void}
    releaseArray::Ptr{Void}
    getDims::Ptr{Void}
    getDataPtr::Ptr{Void}
    getType::Ptr{Void}
    getNumDims::Ptr{Void}
    getElements::Ptr{Void}
    isEmpty::Ptr{Void}
    eval::Ptr{Void}
    freeList::AFImpl
    batch::Bool
    results::Results

    function ArrayFire()
        ptr = getLibPtr(T)
        new(
            ptr,
            AFDevice(ptr),
            Create(ptr),
            Binary(ptr),
            Unary(ptr),
            Index(ptr),
            VectorAlgos(ptr),
            Modify(ptr),
            Reduction(ptr),
            Math(ptr),
            LinAlg(ptr),
            MoveAndReorder(ptr),
            Statistics(ptr),
            Libdl.dlsym(ptr, :af_create_handle),
            Libdl.dlsym(ptr, :af_create_array),
            Libdl.dlsym(ptr, :af_retain_array),
            Libdl.dlsym(ptr, :af_release_array),
            Libdl.dlsym(ptr, :af_get_dims),
            Libdl.dlsym(ptr, :af_get_data_ptr),
            Libdl.dlsym(ptr, :af_get_type),
            Libdl.dlsym(ptr, :af_get_numdims),
            Libdl.dlsym(ptr, :af_get_elements),
            Libdl.dlsym(ptr, :af_is_empty),
            Libdl.dlsym(ptr, :af_eval),
            FreeList(),
            false,
            Results())
    end
end

getBackend{T<:Backend}(af::ArrayFire{T}) = T

function getSupportedBackends()
    backends = Vector{Any}()
    for b in [CPU, CUDA, OpenCL]
        try
            ArrayFire{b}()
            push!(backends, b)
        catch
        end
    end
    backends
end

include("AFDevice.jl")
include("AFArray.jl")
include("macros.jl")
include("Create.jl")
include("Binary.jl")
include("Unary.jl")
include("Index.jl")
include("VectorAlgos.jl")
include("Modify.jl")
include("Reduction.jl")
include("Math.jl")
include("LinAlg.jl")
include("MoveAndReorder.jl")
include("Statistics.jl")
include("FreeList.jl")

immutable ScopeHandle
    result::Function
    register::Function
end

function scope!(pred, af::ArrayFire)
    newScope!(af.freeList)
    try
        h = ScopeHandle(arr -> markResult!(af.freeList, arr), arr -> register!(af.freeList, arr))
        pred(h)
    finally
        endScope!(af.freeList)
    end
end

function register!(af::ArrayFire, arr::AFArray)
    register!(af.freeList, arr)
end

macro scope(af, expr)
    :( scope!(this -> $expr, $af) )
end

end # module

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
using Bokeh

plotfile("simple_glyphs.html")

x = linspace(0, 2pi)
y1 = sin(x)
y2 = cos(x)
y3 = tan(x)
plot([y1 y2 y3]) # since we haven't passed x this would use the range 1:50

plot(x, [y1 y2 y3], "rs|bo|g*")
showplot()


text:
struct SimpleBifurcationInterval{uType}
    i::Int
    u0::uType
    u1::uType
    h::Float64  # hType?
    direction::Int
    sweep::WeakRef
end

struct ContinuationSweep{uType <: AbstractVector}
    i::Int
    direction::Int
    u::Vector{uType}
    simple_bifurcation::Vector{SimpleBifurcationInterval{uType}}
    sol::WeakRef
end

ContinuationSweep(uType::Type{<: AbstractArray}, i, direction, sol) =
    ContinuationSweep(i,
                      direction,
                      Vector{uType}(),
                      Vector{SimpleBifurcationInterval{uType}}(),
                      sol)

Base.eltype(::Type{<: ContinuationSweep{uType}}) where uType = uType
Base.length(sweep::ContinuationSweep) = length(sweep.u)

struct ContinuationSolution{uType <: AbstractVector}
    sweeps::Vector{ContinuationSweep{uType}}
    prob
end

ContinuationSolution(uType::Type{<: AbstractArray}, prob) =
    ContinuationSolution(ContinuationSweep{uType}[], prob)

function push_point!(sweep::ContinuationSweep, cache)
    push_point!(sweep, cache.u)
    if cache.simple_bifurcation
        u0 = sweep.u[end-1]
        u1 = sweep.u[end]
        push!(sweep.simple_bifurcation, SimpleBifurcationInterval(
            length(sweep.u) - 1,
            u0,
            u1,
            cache.h,
            -cache.direction,  # re-flip to save direction at u0
            WeakRef(sweep),
        ))
    end
end

function push_point!(sweep::ContinuationSweep, u::AbstractArray)
    push!(sweep.u, copy(u))
end

function push_point!(sol::ContinuationSolution, args...)
    push_point!(sol.sweeps[end], args...)
end

function new_sweep!(sol::ContinuationSolution{uType},
                    direction) where uType
    push!(sol.sweeps, ContinuationSweep(
        uType,
        length(sol.sweeps) + 1,
        direction,
        WeakRef(sol)))
end

function sweeps_as_vectors(sol::ContinuationSolution{uType}, i) where uType
    if i == 0
        i = length(sol.sweeps[1].u[1])
    end
    x1 = [x[i] for x in sol.sweeps[1].u]
    x2 = [x[i] for x in sol.sweeps[2].u]
    curves = [vcat(reverse!(x2), x1)]
    for sweep in sol.sweeps[3:end]
        push!(curves, [x[i] for x in sweep.u])
    end
    return curves
end

abstract type LastPoint end

struct ContinuationLastPoint{S <: ContinuationSweep} <: LastPoint
    sweep::S
    is_simple_bifurcation::Bool
end

Base.propertynames(::ContinuationLastPoint) =
    (:i_sweep, :i_point, :u, :direction, :simple_bifurcation)

function Base.getproperty(point::ContinuationLastPoint, name::Symbol)
    if name in (:sweep, :is_simple_bifurcation)
        return getfield(point, name)
    end
    sweep = getfield(point, :sweep)
    is_simple_bifurcation = getfield(point, :is_simple_bifurcation)
    if name === :i_sweep
        return sweep.i
    elseif name === :i_point
        return length(sweep.u)
    elseif name === :u
        return sweep.u[end]
    elseif name === :direction
        return sweep.direction
    elseif name === :simple_bifurcation
        return is_simple_bifurcation ? sweep.simple_bifurcation[end] : nothing
    end
    throw(KeyError(name))
end


text:
using CUDAnative
using CuArrays
CuArrays.allowscalar(false)
try
  isdefined(NumericalFluxes, :rusanov!)
catch
  include("../../src/DGmethods/NumericalFluxes.jl")
end

include("DGBalanceLawProfileHarness.jl")

const _nstate = 5
const _ρ, _U, _V, _W, _E = 1:_nstate
const stateid = (ρid = _ρ, Uid = _U, Vid = _V, Wid = _W, Eid = _E)
const statenames = ("ρ", "U", "V", "W", "E")

const _nviscstates = 6
const _τ11, _τ22, _τ33, _τ12, _τ13, _τ23 = 1:_nviscstates

const _ngradstates = 3
const _states_for_gradient_transform = (_ρ, _U, _V, _W)

const _nauxstate = 3
const _a_x, _a_y, _a_z = 1:_nauxstate

const γ_exact = 7//5
const μ_exact = 1//100

@inline function preflux(Q, _...)
  γ::eltype(Q) = γ_exact
  @inbounds ρ, U, V, W, E = Q[_ρ], Q[_U], Q[_V], Q[_W], Q[_E]
  ρinv = 1 / ρ
  u, v, w = ρinv * U, ρinv * V, ρinv * W
  ((γ-1)*(E - ρinv * (U^2 + V^2 + W^2) / 2), u, v, w, ρinv)
end

# max eigenvalue
@inline function wavespeed(n, Q, aux, t, P, u, v, w, ρinv)
  γ::eltype(Q) = γ_exact
  @inbounds abs(n[1] * u + n[2] * v + n[3] * w) + sqrt(ρinv * γ * P)
end

# flux function
cns_flux!(F, Q, VF, aux, t) = cns_flux!(F, Q, VF, aux, t, preflux(Q)...)

@inline function cns_flux!(F, Q, VF, aux, t, P, u, v, w, ρinv)
  @inbounds begin
    ρ, U, V, W, E = Q[_ρ], Q[_U], Q[_V], Q[_W], Q[_E]

    τ11, τ22, τ33 = VF[_τ11], VF[_τ22], VF[_τ33]
    τ12 = τ21 = VF[_τ12]
    τ13 = τ31 = VF[_τ13]
    τ23 = τ32 = VF[_τ23]

    # inviscid terms
    F[1, _ρ], F[2, _ρ], F[3, _ρ] = U          , V          , W
    F[1, _U], F[2, _U], F[3, _U] = u * U  + P , v * U      , w * U
    F[1, _V], F[2, _V], F[3, _V] = u * V      , v * V + P  , w * V
    F[1, _W], F[2, _W], F[3, _W] = u * W      , v * W      , w * W + P
    F[1, _E], F[2, _E], F[3, _E] = u * (E + P), v * (E + P), w * (E + P)

    # viscous terms
    F[1, _U] -= τ11; F[2, _U] -= τ12; F[3, _U] -= τ13
    F[1, _V] -= τ21; F[2, _V] -= τ22; F[3, _V] -= τ23
    F[1, _W] -= τ31; F[2, _W] -= τ32; F[3, _W] -= τ33

    F[1, _E] -= u * τ11 + v * τ12 + w * τ13
    F[2, _E] -= u * τ21 + v * τ22 + w * τ23
    F[3, _E] -= u * τ31 + v * τ32 + w * τ33
  end
end

# Compute the velocity from the state
@inline function velocities!(vel, Q, _...)
  @inbounds begin
    # ordering should match states_for_gradient_transform
    ρ, U, V, W = Q[1], Q[2], Q[3], Q[4]
    ρinv = 1 / ρ
    vel[1], vel[2], vel[3] = ρinv * U, ρinv * V, ρinv * W
  end
end

# Visous flux
@inline function compute_stresses!(VF, grad_vel, _...)
  μ::eltype(VF) = μ_exact
  @inbounds begin
    dudx, dudy, dudz = grad_vel[1, 1], grad_vel[2, 1], grad_vel[3, 1]
    dvdx, dvdy, dvdz = grad_vel[1, 2], grad_vel[2, 2], grad_vel[3, 2]
    dwdx, dwdy, dwdz = grad_vel[1, 3], grad_vel[2, 3], grad_vel[3, 3]

    # strains
    ϵ11 = dudx
    ϵ22 = dvdy
    ϵ33 = dwdz
    ϵ12 = (dudy + dvdx) / 2
    ϵ13 = (dudz + dwdx) / 2
    ϵ23 = (dvdz + dwdy) / 2

    # deviatoric stresses
    VF[_τ11] = 2μ * (ϵ11 - (ϵ11 + ϵ22 + ϵ33) / 3)
    VF[_τ22] = 2μ * (ϵ22 - (ϵ11 + ϵ22 + ϵ33) / 3)
    VF[_τ33] = 2μ * (ϵ33 - (ϵ11 + ϵ22 + ϵ33) / 3)
    VF[_τ12] = 2μ * ϵ12
    VF[_τ13] = 2μ * ϵ13
    VF[_τ23] = 2μ * ϵ23
  end
end

@inline function stresses_penalty!(VF, nM, velM, QM, aM, velP, QP, aP, t)
  @inbounds begin
    n_Δvel = similar(VF, Size(3, 3))
    for j = 1:3, i = 1:3
      n_Δvel[i, j] = nM[i] * (velP[j] - velM[j]) / 2
    end
    compute_stresses!(VF, n_Δvel)
  end
end

numerical_flux!(x...) = NumericalFluxes.rusanov!(x..., cns_flux!, wavespeed,
                                                 preflux)
let
  DFloat = Float64
  dim = 2
  nelem = 100
  N = 4
  cpu_dg = DGProfiler(Array, DFloat, dim, nelem, N, _nstate, cns_flux!,
                      numerical_flux!;
                      stateoffset = ((_E, 20), (_ρ, 1)),
                      ngradstate = _ngradstates,
                      states_grad = _states_for_gradient_transform,
                      nviscstate = _nviscstates,
                      gradient_transform! = velocities!,
                      viscous_transform! = compute_stresses!,
                      viscous_penalty! = stresses_penalty!,
                      nauxstate = _nauxstate)

  gpu_dg = DGProfiler(CuArray, DFloat, dim, nelem, N, _nstate, cns_flux!,
                      numerical_flux!;
                      stateoffset = ((_E, 20), (_ρ, 1)),
                      ngradstate = _ngradstates,
                      states_grad = _states_for_gradient_transform,
                      nviscstate = _nviscstates,
                      gradient_transform! = velocities!,
                      viscous_transform! = compute_stresses!,
                      viscous_penalty! = stresses_penalty!,
                      nauxstate = _nauxstate)

  volumerhs!(cpu_dg)
  facerhs!(cpu_dg)
  volumeviscterms!(cpu_dg)
  faceviscterms!(cpu_dg)

  volumerhs!(gpu_dg)
  facerhs!(gpu_dg)
  volumeviscterms!(gpu_dg)
  faceviscterms!(gpu_dg)

  gpu_rhs = Array(gpu_dg.rhs)
  cpu_rhs = Array(cpu_dg.rhs)
  @show gpu_rhs ≈ cpu_rhs
  nothing
end

max_stars_repo_name: jkozdon/CLIMA

max_issues_repo_name: jkozdon/CLIMA

max_forks_repo_name: jkozdon/CLIMA