text:
export
    AbstractTreeNode,
    get_parent,
    get_children,
    rem_parent!,
    set_parent!,
    has_child,
    has_parent,
    validate_tree,
    validate_embedded_tree

"""
    abstract type AbstractTreeNode{E,ID}

E is element type, ID is id type
"""
abstract type AbstractTreeNode{ID} end
_id_type(n::AbstractTreeNode{ID}) where {ID} = ID

node_id(node::AbstractTreeNode) = node.id
get_parent(n::AbstractTreeNode) = n.parent
get_children(n::AbstractTreeNode) = n.children
has_child(parent::AbstractTreeNode,child::AbstractTreeNode) = haskey(get_children(parent),node_id(child))
has_parent(child::AbstractTreeNode,parent::AbstractTreeNode) = get_parent(child) == parent
function rem_parent!(child::AbstractTreeNode)
    delete!(get_children(get_parent(child)), node_id(child))
    child.parent = child
end
function set_parent!(child::AbstractTreeNode,parent::AbstractTreeNode)
    @assert !(child === parent)
    rem_parent!(child)
    get_children(parent)[child.id] = child
    child.parent = parent
    return true
end

"""
    get_root_node(n::AbstractTreeNode{E,ID}) where {E,ID}

Return the root node of a tree
"""
function get_root_node(n::AbstractTreeNode)
    # identify the root of the tree
    ids = Set{AbstractID}()
    node = n
    parent_id = nothing
    while !(parent_id === node_id(node))
        parent_id = node_id(node)
        if parent_id in ids
            throw(ErrorException("Tree is cyclic!"))
        end
        push!(ids,parent_id)
        node = get_parent(node)
    end
    return node
end

"""
    validate_tree(n::AbstractTreeNode)

Ensure that the subtree of n is in fact a tree--no cycles, and no duplicate ids
"""
function validate_sub_tree(n::N) where {ID,N<:AbstractTreeNode{ID}}
    # Breadth-first search
    node = n
    frontier = Dict{AbstractID,AbstractTreeNode}(node_id(node)=>node)
    explored = Set{AbstractID}()
    while !isempty(frontier)
        id,node = pop!(frontier)
        if node_id(node) in explored
            @warn "$node already explored--is there a cycle?"
            return false
        end
        push!(explored,node_id(node))
        for (child_id,child) in get_children(node)
            push!(frontier,child_id=>child)
        end
    end
    return true
end

"""
    validate_tree(n::AbstractTreeNode)

Ensure that the transform tree is in fact a tree--no cycles, and no duplicate
ids
"""
function validate_tree(n::AbstractTreeNode)
    node = get_root_node(n)
    validate_sub_tree(node)
end

"""
    validate_embedded_tree(graph,f=v->get_node(graph,v))

Verify that all graph edges are mirrored by the parent-child structure stored in
the nodes.
"""
function validate_embedded_tree(graph,f=v->get_node(graph,v),early_stop=false)
    valid = true
    for e in edges(graph)
        src = f(edge_source(e))
        dst = f(edge_target(e))
        if !(has_child(src,dst) && has_parent(dst,src))
            @warn "Property \"has_child(src,dst) && has_parent(dst,src)\" does not hold for edge $e"
            valid = false
        end
        early_stop && !valid ? break : nothing
    end
    for u in vertices(graph)
        parent = f(u)
        for v in vertices(graph)
            child = f(v)
            if has_child(parent,child) || has_parent(child,parent)
                if !(has_child(parent,child) && has_parent(child,parent))
                    if !(parent === child)
                        @warn "has_child($u,$v) = $(has_child(parent,child)) but has_parent($v,$u) = $(has_parent(child,parent))"
                        valid = false
                    end
                end
                if !has_edge(graph,u,v)
                    if !(parent === child)
                        @warn "has_child($u,$v) && has_parent($v,$u) but !has_edge(graph,$u,$v)"
                        valid = false
                    end
                end
            end
            early_stop && !valid ? break : nothing
        end
    end
    return valid
end

"""
	depth_first_search(node::AbstractTreeNode,goal_function,expand_function,
		neighbor_function=outneighbors)

Returns the first vertex satisfying goal_function(graph,v). Only expands v if
expand_function(graph,v) == true.
"""
function depth_first_search(node::AbstractTreeNode,
		goal_function,
		expand_function=v->true,
		explored=Set{_id_type(node)}(),
		skip_first=false,
		)
	if goal_function(node)
		if !(skip_first && isempty(explored))
			return node
		end
    end
    push!(explored,node_id(node))
	if expand_function(node)
		for (_,child) in get_children(node)
			if !(node_id(child) in explored)
				retval = depth_first_search(child,
					goal_function,
					expand_function,
                    explored)
                if !(retval === nothing)
					return retval
				end
			end
		end
	end
	return nothing
end

"""
    has_descendant(node::AbstractTreeNode,other::AbstractTreeNode)

Check if `other` is a descendant of `node`.
"""
function has_descendant(node::AbstractTreeNode,other::AbstractTreeNode)
    val = depth_first_search(node, v->( v === other))
    return !(val === nothing)
end

"""
    has_ancestor(node::AbstractTreeNode,other::AbstractTreeNode)

Check if `other` is an ancestor of `node`.
"""
function has_ancestor(node::AbstractTreeNode,other::AbstractTreeNode)
    has_descendant(other,node)
end

export
    CachedTreeNode,
    cached_element

"""
    CachedTreeNode{ID} <: AbstractTreeNode{ID}

Abstract type representing a node with a cached value. Concrete subtypes have a
cached element (accessed via `cached_element(n)`).
"""
abstract type CachedTreeNode{ID} <: AbstractTreeNode{ID} end

mutable struct TreeNode{E,ID} <: CachedTreeNode{ID}
    id::ID
    element::CachedElement{E}
    parent::TreeNode
    children::Dict{ID,TreeNode{E,ID}}
    function TreeNode{E,ID}(e::E) where {E,ID}
        t = new{E,ID}()
        t.id = get_unique_id(ID)
        t.element = CachedElement(e)
        t.parent = t
        t.children = Dict{ID,TreeNode}()
        return t
    end
end

"""
    cached_element(n::CachedTreeNode)   = n.element

Default method for retrieving the cached element of n.
"""
cached_element(n::CachedTreeNode)   = n.element
time_stamp(n::CachedTreeNode)       = time_stamp(cached_element(n))

"""
    cached_node_up_to_date(n::CachedTreeNode)

Check if n is up to date.
"""
function cached_node_up_to_date(n::CachedTreeNode)
    if is_up_to_date(cached_element(n))
        if !(n === get_parent(n))
            return time_stamp(n) > time_stamp(get_parent(n))
        else
            return true
        end
    end
    return false
end

"""
    set_cached_node_up_to_date!(n::CachedTreeNode,val=true)

Set "up to date" status of n to val.
"""
function set_cached_node_up_to_date!(n::CachedTreeNode,val=true,propagate=true)
    old_val = is_up_to_date(cached_element(n))
    set_up_to_date!(cached_element(n),val)
    # propagate "outdated" signal up the tree
    if propagate && old_val # reduce the amount of wasted forward propagation of "out-of-date" flag
        for (id,child) in get_children(n)
            set_cached_node_up_to_date!(child,false,propagate)
        end
    end
    return n
end

"""
    propagate_backward!(n::CachedTreeNode,args...)

Propagate information up the tree from n.
"""
function propagate_backward!(n::CachedTreeNode,args...)
    propagate_backward!(n,get_parent(n),args...)
end
"""
    propagate_backward!(child::CachedTreeNode,parent::CachedTreeNode,args...)

Propagate information from child to parent.
"""
propagate_backward!(child::CachedTreeNode,parent::CachedTreeNode,args...) = nothing
"""
    propagate_forward!(n::CachedTreeNode,args...)

Propagate information down the tree from n.
"""
function propagate_forward!(n::CachedTreeNode,args...)
    for (id,child) in get_children(n)
        propagate_forward!(n,child,args...)
    end
    return n
end
"""
    propagate_forward!(parent::CachedTreeNode,child::CachedTreeNode,args...)

Propagate information from parent to child.
"""
propagate_forward!(parent::CachedTreeNode,child::CachedTreeNode,args...) = nothing

"""
    update_element!(n::CachedTreeNode,element,args...)

Update cached element of n, and propagate relevant information forward and
backward via `propagate_forward!` and `propagate_backward!`
"""
function update_element!(n::CachedTreeNode,element,propagate=true,args...)
    update_element!(cached_element(n),element)
    set_cached_node_up_to_date!(n,true)
    if propagate
        propagate_backward!(n,args...)
        propagate_forward!(n,args...)
    end
    return n
end

"""
    get_cached_value!(n::CachedTreeNode)

Return the up to date cached value of n. This triggers a reach back to parent
nodes if necessary.
"""
function get_cached_value!(n::CachedTreeNode)
    if !cached_node_up_to_date(n)
        propagate_forward!(get_parent(n),n) # retrieve relevant info from parent
    end
    get_element(cached_element(n))
end


text:
using DrWatson
@quickactivate "ExercisesRepo"
using DelimitedFiles
cd(@__DIR__)

# %% 1 = Mackie Glass timeseries
using DelayDiffEq
function mackey_glass(du,u,h,p,t)
    beta,n,gamma,tau = p
    hist = h(p, t-tau)[1]
    du[1] = (beta*hist)/(1+hist^n) - gamma * u[1]
end
h(p,t) = 0
tau_d = 44
n = 10
β = 0.2
γ = 0.1
δt = 0.5
p = (β,n,γ,tau_d)

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

prob = DDEProblem(mackey_glass,u0,h,tspan,p; constant_lags=tau_d)
alg = MethodOfSteps(Tsit5())
sol = solve(prob,alg; adaptive=false, dt=δt)

s = [u[1] for u in sol.u]
s = s[4001:end]
writedlm("1.csv", s)

# %% 2 = Roessler trajectory (1st parameter set)
using DynamicalSystemsBase
roe = Systems.roessler([1.0, 0, 0]; a=0.2, b=0.2, c=5.7)
sroe = trajectory(roe, 500; dt = 0.05, Ttr = 100.0)
writedlm("2.csv", sroe)

# %% 3 = Lorenz trajectory (1st parameter set, dense time sampling)
lo = Systems.lorenz([0, 10.0, 0.0]; σ=10, ρ=28, β=8/3)
tr = trajectory(lo, 100; dt = 0.01, Ttr = 100)
writedlm("3.csv", tr)

# %% 4 = Roessler trajectory, low sampling, 2nd parameter set
ro = Systems.roessler([1.0, 1.0, 1.0], a=0.15, b = 0.2, c=10)
tr = trajectory(ro, 1250; dt = 0.2, Ttr = 100)
writedlm("4.csv", tr)

# %% 5 = collection of periodic and quasiperiodic orbits
# some from Henon Heiles
using DynamicalSystems, OrdinaryDiffEq, PyPlot, DelimitedFiles
diffeq = (alg = Vern9(), reltol = 1e-9, abstol = 1e-9)
hh = Systems.henonheiles()

xs = Vector{Float64}[]
for u in (
        [0.0, 0.1, 0.5, 0.0], # quasiperiodic
        [0.0, -0.0910355, 0.459522, -0.173339], # higher order quasiperiodic
        [0.0, 0.30266568064689636, 0.42056546332015626, 0.0] # periodic
    )

    tr = trajectory(hh, 500.0, u; Ttr = 10, dt = 0.1, diffeq...)
    x = regularize(tr[:, 1])
    x .+= 0.05randn(length(x))
    push!(xs, x)
end

# Also add a quasi-periodic orbit from Hose Noover for good measure
nh = Systems.nosehoover([0, 1.0, 0])
tr = trajectory(nh, 500.0; Ttr = 10, dt = 0.1, diffeq...)
x = regularize(tr[:, 1])
x .+= 0.05randn(length(x))
push!(xs, x)

# chaotic & period-4 roessler
ro = Systems.roessler(ones(3))
for c in (4.0, 5.7)
    set_parameter!(ro, 3, c)
    tr = trajectory(ro, 500.0; Ttr = 1000, dt = 0.1, diffeq...)
    x = regularize(tr[:, 1])
    x .+= 0.05randn(length(x))
    push!(xs, x)
end

# from standard map (quasiperiod around period 3)
sm = Systems.standardmap()

for (i, u) in enumerate((
        SVector(0.8121, 1.6243),
        SVector(0.877, 1.565),
    ))
    tr = trajectory(sm, 5000, u; Ttr = 10)
    x = regularize(tr[:, 1])
    i == 1 && (x .+= 0.02randn(length(x)))
    push!(xs, x)
end

# Fourth part is just oscillatory ARMA, a-la chapter 7
using Random, Statistics
Random.seed!(77163)
η = randn(15000)
s = ones(15000)
for n in 5:15000
    s[n] = 1.625s[n-1] - 0.284s[n-2] - 0.355s[n-3] + η[n] - 0.26η[n-1]
end
x = s[end-5000:end]
x ./= std(x)
push!(xs, x)

# Last part 3 cosines, to confuse for periodic → quasiperiodic
t = range(0, 100π; length = 5001)
x = @. cos(t) + 0.5cos(2t + π/3) + 0.2cos(13t + π/5)
x ./= std(x)
x .+= 0.05randn(length(x))
push!(xs, x)

writedlm("5.csv", hcat(xs...))

# %% 6 = towel map trajectory
using DynamicalSystems, DelimitedFiles

ds = Systems.towel()

tr = trajectory(ds, 10000; Ttr = 100)

writedlm("6.csv", tr)

# %% 7 = Non-stationary dynamics of logistic map
ds = Systems.logistic()
rs = 3.60:0.0001:3.9
integ = integrator(ds, π/4)
x = [integ.u]; sizehint!(x, length(rs)*3)

for r in rs
    set_parameter!(ds, 1, r)
    for i in 1:3
        DynamicalSystems.step!(integ)
        push!(x, integ.u)
    end
end

x .+= 0.05randn(length(x))

writedlm("7.csv", x)

# %% 8 = oscillatory ARMA, embedded in 4d
using ARFIMA, Random

rng = MersenneTwister(5132786517)
x = arma(rng, 10000, 0.4, SVector(1.624, -0.284, -0.355), SVector(-0.96))
τ = estimate_delay(x, "mi_min")
e = embed(x, 4, τ)

# This looks awesome!
# using3D()
# cols = columns(e)
# scatter3D(cols[1], cols[2], cols[3]; c = cols[4])
writedlm("8.csv", x)

# %% 9 = 6dimensional Lorenz96 with 2% AR noise - Periodic or Quasiperiodic (don't know!)
using DynamicalSystems
D = 6
u = fill(0.1, D)
u[1] = 0.2
u[3] = 0.3
lo = Systems.lorenz96(D, u; F = 5.5)

tr = regularize(trajectory(lo, 100.0; Ttr = 10.0))
co = columns(tr)
for x in co
    noise = 0.02randn(length(x))
    noise = arma(length(x), 0.4, SVector(0.624, -0.284, -0.355, +0.355))
    noise ./= std(noise) / (0.02)
    x .+= noise
end

writedlm("9.csv", hcat(co...))

# scatter3D(co[1], co[2], co[3]; c = co[4])

# %% 10 = 5dimensional Lorenz96 with 2% AR noise - Chaotic trajectory
using DynamicalSystems
D = 5
u = fill(0.1, D)
u[1] = 0.2
u[3] = 0.3
lo = Systems.lorenz96(D, u; F = 8.0)

tr = regularize(trajectory(lo, 100.0; Ttr = 10.0))
co = columns(tr)
for x in co
    noise = 0.02randn(length(x))
    noise = arma(length(x), 0.4, SVector(0.624, -0.284, -0.355, +0.355))
    noise ./= std(noise) / (0.02)
    x .+= noise
end

# scatter3D(co[1], co[2], co[3]; c = co[4])
writedlm("10.csv", hcat(co...))

ext: jl

max_stars_repo_name: moritz-g/ExercisesRepo

max_issues_repo_name: moritz-g/ExercisesRepo

max_forks_repo_name: moritz-g/ExercisesRepo


text:
# spectplan.jl

export SPECTplan

"""
    SPECTplan
Struct for storing key factors for a SPECT system model
- `T` datatype of work arrays
- `imgsize` size of image
- `px,pz` psf dimension
- `imgr [nx, ny, nz]` 3D rotated version of image
- `add_img [nx, ny, nz]` 3D image for adding views and backprojection
- `mumap [nx,ny,nz]` attenuation map, must be 3D, possibly zeros()
- `mumapr [nx, ny, nz]` 3D rotated mumap
- `exp_mumapr [nx, nz]` 2D exponential rotated mumap
- `psfs [px,pz,ny,nview]` point spread function, must be 4D, with `px and `pz` odd, and symmetric for each slice
- `nview` number of views, must be integer
- `viewangle` set of view angles, must be from 0 to 2π
- `interpmeth` interpolation method, :one means 1d, :two means 2d
- `mode` pre-allcoation method, :fast means faster, :mem means use less memory
- `dy` voxel size in y direction (dx is the same value)
- `nthread` number of CPU threads used to process data, must be integer
- `planrot` Vector of struct `PlanRotate`
- `planpsf` Vector of struct `PlanPSF`
Currently code assumes the following:
* each of the `nview` projection views is `[nx,nz]`
* `nx = ny`
* uniform angular sampling
* `psf` is symmetric
* multiprocessing using # of threads specified by `Threads.nthreads()`
"""
struct SPECTplan{T}
    T::DataType # default type for work arrays etc.
    imgsize::NTuple{3, Int}
    px::Int
    pz::Int
    imgr::Union{Array{T, 3}, Vector{Array{T, 3}}} # 3D rotated image, (nx, ny, nz)
    add_img::Union{Array{T, 3}, Vector{Array{T, 3}}}
    mumap::Array{T, 3} # [nx,ny,nz] attenuation map, must be 3D, possibly zeros()
    mumapr::Union{Array{T, 3}, Vector{Array{T, 3}}} # 3D rotated mumap, (nx, ny, nz)
    exp_mumapr::Vector{Matrix{T}} # 2D exponential rotated mumap, (nx, ny)
    psfs::Array{T, 4} # PSFs must be 4D, [px, pz, ny, nview], finally be centered psf
    nview::Int # number of views
    viewangle::StepRangeLen{T}
    interpmeth::Symbol
    mode::Symbol
    dy::T
    nthread::Int # number of threads
    planrot::Vector{PlanRotate}
    planpsf::Vector{PlanPSF}

    """
        SPECTplan(mumap, psfs, dy; T, viewangle, interpmeth, nthread, mode)
    """
    function SPECTplan(
        mumap::Array{<:RealU, 3},
        psfs::Array{<:RealU, 4},
        dy::RealU;
        T::DataType = promote_type(eltype(mumap), Float32),
        viewangle::StepRangeLen{<:RealU} = (0:size(psfs, 4) - 1) / size(psfs, 4) * T(2π), # set of view angles
        interpmeth::Symbol = :two, # :one is for 1d interpolation, :two is for 2d interpolation
        nthread::Int = Threads.nthreads(),
        mode::Symbol = :fast,
    )

        # convert to the same type
        dy = convert(T, dy)
        mumap .= T.(mumap)
        psfs .= T.(psfs)

        (nx, ny, nz) = size(mumap) # typically 128 x 128 x 81

        isequal(nx, ny) || throw("nx != ny")
        (iseven(nx) && iseven(ny)) || throw("nx odd")

        imgsize = (nx, ny, nz)
        # check psf
        px, pz, _, nview = size(psfs)
        (isodd(px) && isodd(pz)) || throw("non-odd size psfs")
        all(mapslices(x -> x == reverse(x, dims=:), psfs, dims = [1, 2])) ||
            throw("asym. psf")

        # check interpidx
        (interpmeth === :one || interpmeth === :two) || throw("bad interpmeth")

        # remember to check if nthread == Threads.nthreads()
        (nthread == Threads.nthreads()) || throw("bad nthread")

        # check mode
        (mode === :fast || mode === :mem) || throw("bad mode")

        if mode === :fast
            # imgr stores 3D image in different view angles
            imgr = [Array{T, 3}(undef, nx, ny, nz) for id = 1:nthread]
            # add_img stores 3d image for backprojection
            add_img = [Array{T, 3}(undef, nx, ny, nz) for id = 1:nthread]
            # mumapr stores 3D mumap in different view angles
            mumapr = [Array{T, 3}(undef, nx, ny, nz) for id = 1:nthread]
        else
            # imgr stores 3D image in different view angles
            imgr = Array{T, 3}(undef, nx, ny, nz)
            # add_img stores 3d image for backprojection
            add_img = Array{T ,3}(undef, nx, ny, nz)
            # mumapr stores 3D mumap in different view angles
            mumapr = Array{T, 3}(undef, nx, ny, nz)
        end

        exp_mumapr = [Matrix{T}(undef, nx, nz) for id = 1:nthread]

        planrot = plan_rotate(nx; T, method = interpmeth)

        planpsf = plan_psf(; nx, nz, px, pz, nthread, T)

        new{T}(T, # default type for work arrays etc.
            imgsize,
            px,
            pz,
            imgr, # 3D rotated image, (nx, ny, nz)
            add_img,
            mumap, # [nx,ny,nz] attenuation map, must be 3D, possibly zeros()
            mumapr, # 3D rotated mumap, (nx, ny, nz)
            exp_mumapr,
            psfs, # PSFs must be 4D, [px, pz, ny, nview], finally be centered psf
            nview, # number of views
            viewangle,
            interpmeth,
            mode,
            dy,
            nthread, # number of threads
            planrot,
            planpsf,
        )
    end
end

"""
    show(io::IO, ::MIME"text/plain", plan::SPECTplan)
"""
function Base.show(io::IO, ::MIME"text/plain", plan::SPECTplan{T}) where {T}
    t = typeof(plan)
    println(io, t)
    for f in (:imgsize, :px, :pz, :nview, :viewangle, :interpmeth, :mode, :dy, :nthread)
        p = getproperty(plan, f)
        t = typeof(p)
        println(io, " ", f, "::", t, " ", p)
    end
    for f in (:mumap, )
        p = getproperty(plan, f)
        println(io, " ", f, ":", " ", summary(p))
    end
    println(io, " (", sizeof(plan), " bytes)")
end

"""
    sizeof(::SPECTplan)
Show size in bytes of `SPECTplan` object.
"""
function Base.sizeof(ob::T) where {T <: Union{SPECTplan}}
    sum(f -> sizeof(getfield(ob, f)), fieldnames(typeof(ob)))
end


text:
# 1D diffusion problem

# Packages and inclusions
using ModelingToolkit,DiffEqOperators,LinearAlgebra,Test,OrdinaryDiffEq, DomainSets
using ModelingToolkit: Differential

# Tests
@testset "Test 00: Dt(u(t,x)) ~ Dxx(u(t,x))" begin
    # Method of Manufactured Solutions
    u_exact = (x,t) -> exp.(-t) * cos.(x)

    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    Dt = Differential(t)
    Dxx = Differential(x)^2

    # 1D PDE and boundary conditions
    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
    bcs = [u(0,x) ~ cos(x),
           u(t,0) ~ exp(-t),
           u(t,Float64(π)) ~ -exp(-t)]

    # Space and time domains
    domains = [t ∈ Interval(0.0,1.0),
               x ∈ Interval(0.0,Float64(π))]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])

    # Method of lines discretization
    dx = range(0.0,Float64(π),length=30)
    order = 2
    discretization = MOLFiniteDifference([x=>dx],t)
    discretization_edge = MOLFiniteDifference([x=>dx],t;grid_align=edge_align)
    # Explicitly specify order of centered difference
    discretization_centered = MOLFiniteDifference([x=>dx],t;centered_order=order)
    # Higher order centered difference
    discretization_centered_order4 = MOLFiniteDifference([x=>dx],t;centered_order=4)

    for disc in [discretization, discretization_edge, discretization_centered, discretization_centered_order4]
        # Convert the PDE problem into an ODE problem
        prob = discretize(pdesys,disc)

        # Solve ODE problem
        sol = solve(prob,Tsit5(),saveat=0.1)

        if disc.grid_align == center_align
            x = dx[2:end-1]
        else
            x = (dx[1:end-1]+dx[2:end])/2
        end
        t = sol.t

        # Test against exact solution
        for i in 1:length(sol)
            exact = u_exact(x, t[i])
            u_approx = sol.u[i]
            @test all(isapprox.(u_approx, exact, atol=0.01))
        end
    end
end

@testset "Test 01: Dt(u(t,x)) ~ D*Dxx(u(t,x))" begin
    # Parameters, variables, and derivatives
    @parameters t x D
    @variables u(..)
    Dt = Differential(t)
    Dxx = Differential(x)^2

    # 1D PDE and boundary conditions
    eq  = Dt(u(t,x)) ~ D*Dxx(u(t,x))
    bcs = [u(0,x) ~ -x*(x-1)*sin(x),
           u(t,0) ~ 0.0,
           u(t,1) ~ 0.0]

    # Space and time domains
    domains = [t ∈ Interval(0.0,1.0),
               x ∈ Interval(0.0,1.0)]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)],[D=>10.0])

    # Method of lines discretization
    dx = 0.1
    order = 2
    discretization = MOLFiniteDifference([x=>dx],t)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    sol = solve(prob,Tsit5(),saveat=0.1)

    # Test
    n = size(sol,1)
    t_f = size(sol,3)
    @test sol[end] ≈ zeros(n) atol = 0.001;
end

# @test_set "Test 02: Dt(u(t,x)) ~ Dx(D(t,x))*Dx(u(t,x))+D(t,x)*Dxx(u(t,x))" begin
@test_broken begin
    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..) D(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dxx = Differential(x)^2

    D = (t,x) -> 0.999 + 0.001 * t * x
    DxD = expand_derivatives(Dx(D(t,x)))

    # 1D PDE and boundary conditions

    eq  = [ Dt(u(t,x)) ~ DxD*Dx(u(t,x))+D(t,x)*Dxx(u(t,x)),]

    bcs = [u(0,x) ~ -x*(x-1)*sin(x),
           u(t,0) ~ 0.0,
           u(t,1) ~ 0.0]

    # Space and time domains
    domains = [t ∈ Interval(0.0,1.0),
               x ∈ Interval(0.0,1.0)]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])

    # Method of lines discretization
    dx = 0.1
    order = 2
    discretization = MOLFiniteDifference([x=>dx],t)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    sol = solve(prob,Tsit5(),saveat=0.1)

    # Test
    n = size(sol,1)
    t_f = size(sol,3)
    @test_broken sol[:,1,t_f] ≈ zeros(n) atol=0.01;
end

@testset "Test 03: Dt(u(t,x)) ~ Dxx(u(t,x)), homogeneous Neumann BCs" begin
    # Method of Manufactured Solutions
    u_exact = (x,t) -> exp.(-t) * cos.(x)

    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dxx = Differential(x)^2

    # 1D PDE and boundary conditions
    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
    bcs = [u(0,x) ~ cos(x),
           Dx(u(t,0)) ~ 0,
           Dx(u(t,Float64(pi))) ~ 0]

    # Space and time domains
    domains = [t ∈ Interval(0.0,1.0),
               x ∈ Interval(0.0,Float64(pi))]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])

    # Method of lines discretization
    dx = range(0.0,Float64(π),length=300)
    order = 2
    discretization = MOLFiniteDifference([x=>dx],t)
    discretization_edge = MOLFiniteDifference([x=>dx],t;grid_align=center_align)

    # Convert the PDE problem into an ODE problem
    for disc in [discretization, discretization_edge]
        prob = discretize(pdesys,disc)

        # Solve ODE problem
        sol = solve(prob,Tsit5(),saveat=0.1)

        if disc.grid_align == center_align
            x_sol = dx[2:end-1]
        else
            x_sol = (dx[1:end-1]+dx[2:end])/2
        end
        t_sol = sol.t

        # Test against exact solution
        for i in 1:length(sol)
            exact = u_exact(x_sol, t_sol[i])
            u_approx = sol.u[i]
            @test all(isapprox.(u_approx, exact, atol=0.01))
            @test sum(u_approx) ≈ 0 atol=1e-10
        end
    end
end

@testset "Test 03a: Dt(u(t,x)) ~ Dxx(u(t,x)), Neumann BCs" begin
    # Method of Manufactured Solutions
    u_exact = (x,t) -> exp.(-t) * sin.(x)

    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dxx = Differential(x)^2

    # 1D PDE and boundary conditions
    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
    bcs = [u(0,x) ~ sin(x),
           Dx(u(t,0)) ~ exp(-t),
           Dx(u(t,Float64(pi))) ~ -exp(-t)]

    # Space and time domains
    domains = [t ∈ Interval(0.0,1.0),
               x ∈ Interval(0.0,Float64(pi))]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])

    # Method of lines discretization
    dx = range(0.0,Float64(π),length=30)
    order = 2
    discretization = MOLFiniteDifference([x=>dx],t)
    discretization_edge = MOLFiniteDifference([x=>dx],t;grid_align=edge_align)

    # Convert the PDE problem into an ODE problem
    for disc ∈ [discretization, discretization_edge]
        prob = discretize(pdesys,disc)

        # Solve ODE problem
        sol = solve(prob,Tsit5(),saveat=0.1)

        if disc.grid_align == center_align
            x = dx[2:end-1]
        else
            x = (dx[1:end-1]+dx[2:end])/2
        end
        t = sol.t

        # Test against exact solution
        # exact integral based on Neumann BCs
        integral_u_exact = t -> sum(sol.u[1] * dx[2]) + 2 * (exp(-t) - 1)
        for i in 1:length(sol)
            exact = u_exact(x, t[i])
            u_approx = sol.u[i]
            @test all(isapprox.(u_approx, exact, atol=0.01))
            # test mass conservation
            integral_u_approx = sum(u_approx * dx[2])
            @test integral_u_exact(t[i]) ≈ integral_u_approx atol=1e-13
        end
    end
end

@testset "Test 04: Dt(u(t,x)) ~ Dxx(u(t,x)), Neumann + Dirichlet BCs" begin
    # Method of Manufactured Solutions
    u_exact = (x,t) -> exp.(-t) * sin.(x)

    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dxx = Differential(x)^2

    # 1D PDE and boundary conditions
    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
    bcs = [u(0,x) ~ sin(x),
           u(t,0) ~ 0.0,
           Dx(u(t,Float64(pi))) ~ -exp(-t)]

    # Space and time domains
    domains = [t ∈ Interval(0.0,1.0),
               x ∈ Interval(0.0,Float64(pi))]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])

    # Method of lines discretization
    dx = range(0.0,Float64(π),length=30)
    order = 2
    discretization = MOLFiniteDifference([x=>dx],t)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    sol = solve(prob,Tsit5(),saveat=0.1)
    x = dx[2:end-1]
    t = sol.t

    # Test against exact solution
    for i in 1:length(sol)
        exact = u_exact(x, t[i])
        u_approx = sol.u[i]
        @test all(isapprox.(u_approx, exact, atol=0.01))
    end
end

@testset "Test 05: Dt(u(t,x)) ~ Dxx(u(t,x)), Robin BCs" begin
    # Method of Manufactured Solutions
    u_exact = (x,t) -> exp.(-t) * sin.(x)

    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dxx = Differential(x)^2

    # 1D PDE and boundary conditions
    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
    bcs = [u(0,x) ~ sin(x),
           u(t,-1.0) + 3Dx(u(t,-1.0)) ~ exp(-t) * (sin(-1.0) + 3cos(-1.0)),
           4u(t,1.0) + Dx(u(t,1.0)) ~ exp(-t) * (4sin(1.0) + cos(1.0))]

    # Space and time domains
    domains = [t ∈ Interval(0.0,1.0),
               x ∈ Interval(-1.0,1.0)]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])

    # Method of lines discretization
    dx = 0.01
    order = 2
    discretization = MOLFiniteDifference([x=>dx],t)
    discretization_edge = MOLFiniteDifference([x=>dx],t;grid_align=edge_align)

    for disc ∈ [discretization, discretization_edge]
        # Convert the PDE problem into an ODE problem
        prob = discretize(pdesys,disc)

        # Solve ODE problem
        sol = solve(prob,Tsit5(),saveat=0.1)
        x = (-1:dx:1)
        if disc.grid_align == center_align
            x = x[2:end-1]
        else
            x = (x[1:end-1].+x[2:end])/2
        end
        t = sol.t

        # Test against exact solution
        for i in 1:length(sol)
            exact = u_exact(x, t[i])
            u_approx = sol.u[i]
            @test all(isapprox.(u_approx, exact, atol=0.1))
        end
    end
end

@testset "Test 06: Dt(u(t,x)) ~ Dxx(u(t,x)), time-dependent Robin BCs" begin
    # Method of Manufactured Solutions
    u_exact = (x,t) -> exp.(-t) * sin.(x)

    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dxx = Differential(x)^2

    # 1D PDE and boundary conditions
    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
    bcs = [u(0,x) ~ sin(x),
          t^2 * u(t,-1.0) + 3Dx(u(t,-1.0)) ~ exp(-t) * (t^2 * sin(-1.0) + 3cos(-1.0)),
          4u(t,1.0) + t * Dx(u(t,1.0)) ~ exp(-t) * (4sin(1.0) + t * cos(1.0))]

    # Space and time domains
    domains = [t ∈ Interval(0.0,1.0),
               x ∈ Interval(-1.0,1.0)]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])

    # Method of lines discretization
    dx = 0.01
    order = 2
    discretization = MOLFiniteDifference([x=>dx],t)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    sol = solve(prob,Rodas4(),reltol=1e-6,saveat=0.1)

    x = (-1:dx:1)
    t = sol.t

    # Test against exact solution
    for i in 1:length(sol)
       exact = u_exact(x, t[i])
       # Due to structural simplification
       # [u2 -> u(n-1), u(1), u(n)]
       # Will be fixed by sol[u]
       u_approx = [sol.u[i][end-1];sol.u[i][1:end-2];sol.u[i][end]]
       @test all(isapprox.(u_approx, exact, atol=0.01))
    end
end

@testset "Test 07: Dt(u(t,r)) ~ 1/r^2 * Dr(r^2 * Dr(u(t,r))) (Spherical Laplacian)" begin
    # Method of Manufactured Solutions
    # general solution of the spherical Laplacian equation
    # satisfies Dr(u(t,0)) = 0
    u_exact = (r,t) -> exp.(-t) * sin.(r) ./ r

    # Parameters, variables, and derivatives
    @parameters t r
    @variables u(..)
    Dt = Differential(t)
    Dr = Differential(r)

    # 1D PDE and boundary conditions

    eq  = Dt(u(t,r)) ~ 1/r^2 * Dr(r^2 * Dr(u(t,r)))
    bcs = [u(0,r) ~ sin(r)/r,
           Dr(u(t,0)) ~ 0,
           u(t,1) ~ exp(-t) * sin(1)]
        #    Dr(u(t,1)) ~ -exp(-t) * sin(1)]

    # Space and time domains
    domains = [t ∈ Interval(0.0,1.0),
               r ∈ Interval(0.0,1.0)]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,r],[u(t,r)])

    # Method of lines discretization
    dr = 0.1
    order = 2
    discretization = MOLFiniteDifference([r=>dr],t)
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    sol = solve(prob,Tsit5(),saveat=0.1)

    r = (0:dr:1)[2:end-1]
    t = sol.t

    # Test against exact solution
    for i in 1:length(sol)
        exact = u_exact(r, t[i])
        u_approx = sol.u[i]
        @test all(isapprox.(u_approx, exact, atol=0.01))
    end
end

@testset "Test 08: Dt(u(t,r)) ~ 4/r^2 * Dr(r^2 * Dr(u(t,r))) (Spherical Laplacian)" begin
    # Method of Manufactured Solutions
    # general solution of the spherical Laplacian equation
    # satisfies Dr(u(t,0)) = 0
    u_exact = (r,t) -> exp.(-4t) * sin.(r) ./ r

    # Parameters, variables, and derivatives
    @parameters t r
    @variables u(..)
    Dt = Differential(t)
    Dr = Differential(r)

    # 1D PDE and boundary conditions

    eq  = Dt(u(t,r)) ~ 4/r^2 * Dr(r^2 * Dr(u(t,r)))
    bcs = [u(0,r) ~ sin(r)/r,
           Dr(u(t,0)) ~ 0,
           u(t,1) ~ exp(-4t) * sin(1)]

    # Space and time domains
    domains = [t ∈ Interval(0.0,1.0),
               r ∈ Interval(0.0,1.0)]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,r],[u(t,r)])

    # Method of lines discretization
    dr = 0.1
    order = 2
    discretization = MOLFiniteDifference([r=>dr],t)
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    sol = solve(prob,Tsit5(),saveat=0.1)

    r = (0:dr:1)[2:end-1]
    t = sol.t

    # Test against exact solution
    for i in 1:length(sol)
        exact = u_exact(r, t[i])
        u_approx = sol.u[i]
        @test all(isapprox.(u_approx, exact, atol=0.01))
    end
end

@testset "Test 10: linear diffusion, two variables, mixed BCs" begin
    # Method of Manufactured Solutions
    u_exact = (x,t) -> exp.(-t) * cos.(x)
    v_exact = (x,t) -> exp.(-t) * sin.(x)

    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..) v(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dxx = Dx^2

    # 1D PDE and boundary conditions
    eqs = [Dt(u(t,x)) ~ Dxx(u(t,x)),
           Dt(v(t,x)) ~ Dxx(v(t,x))]
    bcs = [u(0,x) ~ cos(x),
           v(0,x) ~ sin(x),
           u(t,0) ~ exp(-t),
           Dx(u(t,1)) ~ -exp(-t) * sin(1),
           Dx(v(t,0)) ~ exp(-t),
           v(t,1) ~ exp(-t) * sin(1)]

    # Space and time domains
    domains = [t ∈ Interval(0.0,1.0),
               x ∈ Interval(0.0,1.0)]

    # PDE system
    pdesys = PDESystem(eqs,bcs,domains,[t,x],[u(t,x),v(t,x)])

    # Method of lines discretization
    l = 100
    dx = range(0.0,1.0,length=l)
    order = 2
    discretization = MOLFiniteDifference([x=>dx],t)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    sol = solve(prob,Tsit5(),saveat=0.1)

    x_sol = dx[2:end-1]
    t_sol = sol.t

    # Test against exact solution
    for i in 1:length(sol)
        @test all(isapprox.(u_exact(x_sol, t_sol[i]), sol.u[i][1:l-2], atol=0.01))
        @test all(isapprox.(v_exact(x_sol, t_sol[i]), sol.u[i][l-1:end], atol=0.01))
    end
end

@testset "Test 11: linear diffusion, two variables, mixed BCs, with parameters" begin
    @parameters t x
    @parameters Dn, Dp
    @variables u(..) v(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dxx = Differential(x)^2

    eqs  = [Dt(u(t,x)) ~ Dn * Dxx(u(t,x)) + u(t,x)*v(t,x),
            Dt(v(t,x)) ~ Dp * Dxx(v(t,x)) - u(t,x)*v(t,x)]
    bcs = [u(0,x) ~ sin(pi*x/2),
        v(0,x) ~ sin(pi*x/2),
        u(t,0) ~ 0.0, Dx(u(t,1)) ~ 0.0,
        v(t,0) ~ 0.0, Dx(v(t,1)) ~ 0.0]

    domains = [t ∈ Interval(0.0,1.0),
            x ∈ Interval(0.0,1.0)]

    pdesys = PDESystem(eqs,bcs,domains,[t,x],[u(t,x),v(t,x)],[Dn=>0.5, Dp=>2])
    discretization = MOLFiniteDifference([x=>0.1],t)
    prob = discretize(pdesys,discretization)
    @test prob.p == [0.5,2]
    # Make sure it can be solved
    sol = solve(prob,Tsit5())
end

@testset "Test 12: linear diffusion, two variables, mixed BCs, different independent variables" begin
    # Method of Manufactured Solutions
    u_exact = (x,t) -> exp.(-t) * cos.(x)
    v_exact = (y,t) -> exp.(-t) * sin.(y)

    # Parameters, variables, and derivatives
    @parameters t x y
    @variables u(..) v(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dxx = Dx^2
    Dy = Differential(y)
    Dyy = Dy^2

    # 1D PDE and boundary conditions
    eqs = [Dt(u(t,x)) ~ Dxx(u(t,x)),
           Dt(v(t,y)) ~ Dyy(v(t,y))]
    bcs = [u(0,x) ~ cos(x),
           v(0,y) ~ sin(y),
           u(t,0) ~ exp(-t),
           Dx(u(t,1)) ~ -exp(-t) * sin(1),
           Dy(v(t,0)) ~ exp(-t),
           v(t,2) ~ exp(-t) * sin(2)]

    # Space and time domains
    domains = [t ∈ Interval(0.0,1.0),
               x ∈ Interval(0.0,1.0),
               y ∈ Interval(0.0,2.0)]

    # PDE system
    pdesys = PDESystem(eqs,bcs,domains,[t,x,y],[u(t,x),v(t,y)])

    # Method of lines discretization
    l = 100
    dx = range(0.0,1.0,length=l)
    dy = range(0.0,2.0,length=l)
    order = 2
    discretization = MOLFiniteDifference([x=>dx,y=>dy],t)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    sol = solve(prob,Tsit5(),saveat=0.1)

    x_sol = dx[2:end-1]
    y_sol = dy[2:end-1]
    t_sol = sol.t

    # Test against exact solution
    for i in 1:length(sol)
        @test all(isapprox.(u_exact(x_sol, t_sol[i]), sol.u[i][1:l-2], atol=0.01))
        @test all(isapprox.(v_exact(y_sol, t_sol[i]), sol.u[i][l-1:end], atol=0.01))
    end
end

@testset "Test 13: one linear diffusion with mixed BCs, one ODE" begin
    # Method of Manufactured Solutions
    u_exact = (x,t) -> exp.(-t) * sin.(x)
    v_exact = (t) -> exp.(-t)

    @parameters t x
    @variables u(..) v(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dxx = Dx^2

    # 1D PDE and boundary conditions
    eqs = [Dt(u(t,x)) ~ Dxx(u(t,x)),
           Dt(v(t)) ~ -v(t)]
    bcs = [u(0,x) ~ sin(x),
           v(0) ~ 1,
           u(t,0) ~ 0,
           Dx(u(t,1)) ~ exp(-t) * cos(1)]

    # Space and time domains
    domains = [t ∈ Interval(0.0,1.0),
               x ∈ Interval(0.0,1.0)]

    # PDE system
    pdesys = PDESystem(eqs,bcs,domains,[t,x],[u(t,x),v(t)])

    # Method of lines discretization
    l = 100
    dx = range(0.0,1.0,length=l)
    order = 2
    discretization = MOLFiniteDifference([x=>dx],t)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    sol = solve(prob,Tsit5(),saveat=0.1)

    x_sol = dx[2:end-1]
    t_sol = sol.t

    # Test against exact solution
    for i in 1:length(sol)
        @test all(isapprox.(u_exact(x_sol, t_sol[i]), sol.u[i][1:length(x_sol)], atol=0.01))
        @test v_exact(t_sol[i]) ≈ sol.u[i][end] atol=0.01
    end
end

@testset "Test error 01: Test Invalid Centered Order" begin
    # Method of Manufactured Solutions
    u_exact = (x,t) -> exp.(-t) * cos.(x)

    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    Dt = Differential(t)
    Dxx = Differential(x)^2

    # 1D PDE and boundary conditions
    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
    bcs = [u(0,x) ~ cos(x),
           u(t,0) ~ exp(-t),
           u(t,Float64(π)) ~ -exp(-t)]

    # Space and time domains
    domains = [t ∈ Interval(0.0,1.0),
               x ∈ Interval(0.0,Float64(π))]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])

    # Method of lines discretization
    dx = range(0.0,Float64(π),length=30)

    # Explicitly specify and invalid order of centered difference
    for order in 1:6
        discretization = MOLFiniteDifference([x=>dx],t;centered_order=order)
        if order % 2 != 0
            @test_throws ArgumentError discretize(pdesys,discretization)
        else
            discretize(pdesys,discretization)
        end
    end
end


text:
using Documenter, RandomNumbers, DocumenterMarkdown
import Random123

makedocs(
    modules = [RandomNumbers],
    format = :markdown
)

deploydocs(
    repo   = "github.com/sunoru/RandomNumbers.jl.git",
    deps   = Deps.pip("pygments", "mkdocs==0.17.5", "python-markdown-math", "mkdocs-material==2.9.4"),
    make = () -> run(`mkdocs build`),
    target = "site"
)


text:
urls = [
    "https://github.com/tawheeler/Vec.jl.git",
    "https://github.com/tawheeler/AutomotiveDrivingModels.jl.git",
    "https://github.com/wulfebw/AutoRisk.jl.git",
    "https://github.com/wulfebw/Coord.jl.git"
]

for url in urls
    try
        Pkg.clone(url)
    catch e
        println("Exception when cloning $(url): $(e)")
    end
end

Pkg.build("AutomotiveDrivingModels")
Pkg.build("AutoRisk")
Pkg.build("Coord")


text:
rot(θ) = [cos(θ) -sin(θ); sin(θ) cos(θ)]

convert(::Type{Matrix}, dgm::PersistenceDiagram; skipinf=true) =
    hcat(( [birth(i), death(i)] for i in dgm if !skipinf || !isinf(birthx(i)) )...)

"""
    wasserstein(dgm1, dgm2)

Calculate Wasserstein distance between persistent diagrams `dgm1` and `dgm2`.
"""
function wasserstein(dgm1::PersistenceDiagram{T}, dgm2::PersistenceDiagram{T}) where {T}
    m1 = convert(Matrix, dgm1)
    m2 = convert(Matrix, dgm2)
    m, n = min(size(m1)...), min(size(m2)...)
    if m == 0
        m = 1
        m1 = zeros(T, 2, m)
    end
    if n == 0
        n = 1
        m2 = zeros(T, 2, n)
    end

    CSM = fill(0.0, n+m, n+m)
    ul = view(CSM,1:m,1:n)
    pairwise!(ul, Euclidean(), m1, m2, dims=2)
    maxd = maximum(ul)

    ur = @view CSM[1:m,(n+1):m+n]
    ur[:,:] .= maxd.*ones(m,m)
    view(ur, diagind(ur,0)) .= (m1'*rot(π/4))[:,2]

    ll = @view CSM[(m+1):m+n,1:n]
    ll[:,:] .= maxd.*ones(n,n)
    view(ll, diagind(ll,0)) .= (m2'*rot(π/4))[:,2]

    # Run the hungarian algorithm
    _, dist = hungarian(CSM)

    return dist
end

"""
    bottleneck(dgm1, dgm2)

Calculate bottleneck distance between persistent diagrams `dgm1` and `dgm2`.
"""
function bottleneck(dgm1::PersistenceDiagram{T}, dgm2::PersistenceDiagram{T}) where {T}
    return zero(T)
end


text:
module TestFMDatasetML100K

using FactorizationMachines
using Base.Test

info("Testing reading from libsvm format...")

# read data in columnar sparse matrix
XFromFile, yFromFile = read_libsvm("data/small_train.libfm")

@test typeof(XFromFile) == SparseMatrixCSC{Float64,Int64}
I, J, V = findnz(XFromFile)

# check observations
@test size(XFromFile) == (10, 5)
@test I == [1,3,7,10,1,4,7,10,1,5,7,10,2,3,9,10,2,4,9,10]
@test J == [1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5]
@test V == [1.0,1.0,1.0,12.5,1.0,1.0,1.0,20.0,1.0,1.0,1.0,78.0,1.0,1.0,1.0,12.5,1.0,1.0,1.0,20.0]

#check predictions
@test typeof(yFromFile) == Array{Float64,1}
@test size(yFromFile, 1) == 5

info("Test reading from libsvm formate ended.")
end


text:
@testset "4.1.4.1 (a+b sin)^m (A+B sin+C sin^2)" begin
    (A, B, C, a, b, e, f, m, x, ) = @variables A B C a b e f m x

    #= ::Package:: =#

    #= ::Title::Closed:: =#
    #=Integrands*of*the*form*(b*sin(e+f*x))^m*(A+C*sin(e+f*x)^2)=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(b*sin(e+f*x))^m*(A+C*sin(e+f*x)^2)*when*A*(m+2)+C*(m+1)=0=#

    @test_int [sin(e + f*x)^m*(1 + m - (2 + m)*sin(e + f*x)^2), x, 1, (cos(e + f*x)*sin(e + f*x)^(1 + m))/f]

    @test_int [sin(e + f*x)^5*(6 - 7*sin(e + f*x)^2), x, 1, (cos(e + f*x)*sin(e + f*x)^6)/f]
    @test_int [sin(e + f*x)^4*(5 - 6*sin(e + f*x)^2), x, 1, (cos(e + f*x)*sin(e + f*x)^5)/f]
    @test_int [sin(e + f*x)^3*(4 - 5*sin(e + f*x)^2), x, 1, (cos(e + f*x)*sin(e + f*x)^4)/f]
    @test_int [sin(e + f*x)^2*(3 - 4*sin(e + f*x)^2), x, 1, (cos(e + f*x)*sin(e + f*x)^3)/f]
    @test_int [sin(e + f*x)^1*(2 - 3*sin(e + f*x)^2), x, 1, (cos(e + f*x)*sin(e + f*x)^2)/f]
    @test_int [sin(e + f*x)^0*(1 - 2*sin(e + f*x)^2), x, 3, (cos(e + f*x)*sin(e + f*x))/f]
    @test_int [csc(e + f*x)^1*(0 - 1*sin(e + f*x)^2), x, 1, cos(e + f*x)/f]
    @test_int [csc(e + f*x)^2*(-1 - 0*sin(e + f*x)^2), x, 2, cot(e + f*x)/f]
    @test_int [csc(e + f*x)^3*(-2 + 1*sin(e + f*x)^2), x, 1, (cot(e + f*x)*csc(e + f*x))/f]
    @test_int [csc(e + f*x)^4*(-3 + 2*sin(e + f*x)^2), x, 1, (cot(e + f*x)*csc(e + f*x)^2)/f]
    @test_int [csc(e + f*x)^5*(-4 + 3*sin(e + f*x)^2), x, 1, (cot(e + f*x)*csc(e + f*x)^3)/f]

    #= ::Title::Closed:: =#
    #=Integrands*of*the*form*(a+b*sin(e+f*x))^m*(A+C*sin(e+f*x)^2)=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(a+a*sin(e+f*x))^m*(A+C*sin(e+f*x)^2)=#

    @test_int [(a + a*sin(e + f*x))^m*(A + C*sin(e + f*x)^2), x, 4, If(13>=8, (C*cos(e + f*x)*(a + a*sin(e + f*x))^m)/(f*(2 + 3*m + m^2)) - (2^(1/2 + m)*(C*(1 + m + m^2) + A*(2 + 3*m + m^2))*cos(e + f*x)*HypergeometricFunctions._₂F₁(1/2, 1/2 - m, 3/2, (1/2)*(1 - sin(e + f*x)))*(1 + sin(e + f*x))^(-(1/2) - m)*(a + a*sin(e + f*x))^m)/(f*(1 + m)*(2 + m)) - (C*cos(e + f*x)*(a + a*sin(e + f*x))^(1 + m))/(a*f*(2 + m)), (C*cos(e + f*x)*(a + a*sin(e + f*x))^m)/(f*(2 + 3*m + m^2)) - (2^(1/2 + m)*(C*(1 + m + m^2) + A*(2 + 3*m + m^2))*cos(e + f*x)*HypergeometricFunctions._₂F₁(1/2, 1/2 - m, 3/2, (1/2)*(1 - sin(e + f*x)))*(1 + sin(e + f*x))^(-(1/2) - m)*(a + a*sin(e + f*x))^m)/(f*(2 + 3*m + m^2)) - (C*cos(e + f*x)*(a + a*sin(e + f*x))^(1 + m))/(a*f*(2 + m)))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(a+b*sin(e+f*x))^m*(A+C*sin(e+f*x)^2)=#

    @test_int [(a + b*sin(e + f*x))^m*(A - A*sin(e + f*x)^2), x, 7, (4*sqrt(2)*A*AppellF1(1/2, -(3/2), -m, 3/2, (1/2)*(1 - sin(e + f*x)), (b*(1 - sin(e + f*x)))/(a + b))*cos(e + f*x)*(a + b*sin(e + f*x))^m)/(((a + b*sin(e + f*x))/(a + b))^m*(f*sqrt(1 + sin(e + f*x)))) - (4*sqrt(2)*A*AppellF1(1/2, -(1/2), -m, 3/2, (1/2)*(1 - sin(e + f*x)), (b*(1 - sin(e + f*x)))/(a + b))*cos(e + f*x)*(a + b*sin(e + f*x))^m)/(((a + b*sin(e + f*x))/(a + b))^m*(f*sqrt(1 + sin(e + f*x))))]

    @test_int [(a + b*sin(e + f*x))^m*(A + C*sin(e + f*x)^2), x, 8, -((C*cos(e + f*x)*(a + b*sin(e + f*x))^(1 + m))/(b*f*(2 + m))) + (sqrt(2)*a*(a + b)*C*AppellF1(1/2, 1/2, -1 - m, 3/2, (1/2)*(1 - sin(e + f*x)), (b*(1 - sin(e + f*x)))/(a + b))*cos(e + f*x)*(a + b*sin(e + f*x))^m)/(((a + b*sin(e + f*x))/(a + b))^m*(b^2*f*(2 + m)*sqrt(1 + sin(e + f*x)))) - (sqrt(2)*(a^2*C + b^2*(C*(1 + m) + A*(2 + m)))*AppellF1(1/2, 1/2, -m, 3/2, (1/2)*(1 - sin(e + f*x)), (b*(1 - sin(e + f*x)))/(a + b))*cos(e + f*x)*(a + b*sin(e + f*x))^m)/(((a + b*sin(e + f*x))/(a + b))^m*(b^2*f*(2 + m)*sqrt(1 + sin(e + f*x))))]

    @test_int [sin(e + f*x)^5*(A + C*sin(e + f*x)^2), x, 3, -(((A + C)*cos(e + f*x))/f) + ((2*A + 3*C)*cos(e + f*x)^3)/(3*f) - ((A + 3*C)*cos(e + f*x)^5)/(5*f) + (C*cos(e + f*x)^7)/(7*f)]

    #= ::Title::Closed:: =#
    #=Integrands*of*the*form*(a+b*sin(e+f*x))^m*(A+B*sin(e+f*x)+C*sin(e+f*x)^2)=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(a+a*sin(e+f*x))^m*(A+B*sin(e+f*x)+C*sin(e+f*x)^2)=#

    @test_int [(a + a*sin(e + f*x))^m*(A + B*sin(e + f*x) + C*sin(e + f*x)^2), x, 4, If(13>=8, ((C - B*(2 + m))*cos(e + f*x)*(a + a*sin(e + f*x))^m)/(f*(1 + m)*(2 + m)) - (2^(1/2 + m)*(B*m*(2 + m) + C*(1 + m + m^2) + A*(2 + 3*m + m^2))*cos(e + f*x)*HypergeometricFunctions._₂F₁(1/2, 1/2 - m, 3/2, (1/2)*(1 - sin(e + f*x)))*(1 + sin(e + f*x))^(-(1/2) - m)*(a + a*sin(e + f*x))^m)/(f*(1 + m)*(2 + m)) - (C*cos(e + f*x)*(a + a*sin(e + f*x))^(1 + m))/(a*f*(2 + m)), ((C - B*(2 + m))*cos(e + f*x)*(a + a*sin(e + f*x))^m)/(f*(1 + m)*(2 + m)) - (2^(1/2 + m)*(B*m*(2 + m) + C*(1 + m + m^2) + A*(2 + 3*m + m^2))*cos(e + f*x)*HypergeometricFunctions._₂F₁(1/2, 1/2 - m, 3/2, (1/2)*(1 - sin(e + f*x)))*(1 + sin(e + f*x))^(-(1/2) - m)*(a + a*sin(e + f*x))^m)/(f*(2 + 3*m + m^2)) - (C*cos(e + f*x)*(a + a*sin(e + f*x))^(1 + m))/(a*f*(2 + m)))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(a+b*sin(e+f*x))^m*(A+B*sin(e+f*x)+C*sin(e+f*x)^2)=#

    @test_int [(a + b*sin(e + f*x))^m*(A + (A + C)*sin(e + f*x) + C*sin(e + f*x)^2), x, 7, -((4*sqrt(2)*C*AppellF1(1/2, -(3/2), -m, 3/2, (1/2)*(1 - sin(e + f*x)), (b*(1 - sin(e + f*x)))/(a + b))*cos(e + f*x)*(a + b*sin(e + f*x))^m)/(((a + b*sin(e + f*x))/(a + b))^m*(f*sqrt(1 + sin(e + f*x))))) - (2*sqrt(2)*(A - C)*AppellF1(1/2, -(1/2), -m, 3/2, (1/2)*(1 - sin(e + f*x)), (b*(1 - sin(e + f*x)))/(a + b))*cos(e + f*x)*(a + b*sin(e + f*x))^m)/(((a + b*sin(e + f*x))/(a + b))^m*(f*sqrt(1 + sin(e + f*x))))]

    @test_int [(a + b*sin(e + f*x))^m*(A + B*sin(e + f*x) + C*sin(e + f*x)^2), x, 8, -((C*cos(e + f*x)*(a + b*sin(e + f*x))^(1 + m))/(b*f*(2 + m))) + (sqrt(2)*(a + b)*(a*C - b*B*(2 + m))*AppellF1(1/2, 1/2, -1 - m, 3/2, (1/2)*(1 - sin(e + f*x)), (b*(1 - sin(e + f*x)))/(a + b))*cos(e + f*x)*(a + b*sin(e + f*x))^m)/(((a + b*sin(e + f*x))/(a + b))^m*(b^2*f*(2 + m)*sqrt(1 + sin(e + f*x)))) - (sqrt(2)*(a^2*C + b^2*C*(1 + m) + A*b^2*(2 + m) - a*b*B*(2 + m))*AppellF1(1/2, 1/2, -m, 3/2, (1/2)*(1 - sin(e + f*x)), (b*(1 - sin(e + f*x)))/(a + b))*cos(e + f*x)*(a + b*sin(e + f*x))^m)/(((a + b*sin(e + f*x))/(a + b))^m*(b^2*f*(2 + m)*sqrt(1 + sin(e + f*x))))]
    end

max_stars_repo_path: test/4 Trig functions/4.1 Sine/4.1.4.1 (a+b sin)^m (A+B sin+C sin^2).jl


text:
using YaoHIR: X, Z, S, T, SWAP, Rz, Rx, shift

ZXDiagram(bir::BlockIR) = convert_to_zxd(bir)
Chain(zxd::ZXDiagram) = convert_to_chain(zxd)

convert_to_gate(::Val{:X}, loc) = Gate(X, Locations(loc))
convert_to_gate(::Val{:Z}, loc) = Gate(Z, Locations(loc))
convert_to_gate(::Val{:H}, loc) = Gate(H, Locations(loc))
convert_to_gate(::Val{:S}, loc) = Gate(S, Locations(loc))
convert_to_gate(::Val{:Sdag}, loc) = Gate(AdjointOperation(S), Locations(loc))
convert_to_gate(::Val{:T}, loc) =  Gate(T, Locations(loc))
convert_to_gate(::Val{:Tdag}, loc) = Gate(AdjointOperation(S), Locations(loc))
convert_to_gate(::Val{:SWAP}, loc1, loc2) = Gate(SWAP, Locations((loc1, loc2)))
convert_to_gate(::Val{:CNOT}, loc, ctrl) = Ctrl(Gate(X, Locations(loc)), CtrlLocations(ctrl))
convert_to_gate(::Val{:CZ}, loc, ctrl) = Ctrl(Gate(Z, Locations(loc)), CtrlLocations(ctrl))
function convert_to_gate(::Val{:Rz}, loc, theta)
    if theta isa Phase
        theta = theta * π
        theta = theta.ex
    end
    return Gate(Rz(theta), Locations(loc))
end
function convert_to_gate(::Val{:Rx}, loc, theta)
    if theta isa Phase
        theta = theta * π
        theta = theta.ex
    end
    return Gate(Rx(theta), Locations(loc))
end
function convert_to_gate(::Val{:shift}, loc, theta)
    if theta isa Phase
        theta = theta * π
        theta = theta.ex
    end
    return Gate(shift(theta), Locations(loc))
end

function push_gate!(circ::Chain, gargs...)
    push!(circ.args, convert_to_gate(gargs...))
    return circ
end
function pushfirst_gate!(circ::Chain, gargs...)
    pushfirst!(circ.args, convert_to_gate(gargs...))
    return circ
end

function unwrap_ssa_phase(theta, ir::Core.Compiler.IRCode)
    if theta isa Core.SSAValue
        return Phase(theta, ir.stmts[theta.id][:type])
    elseif theta isa QuoteNode
        return theta.value
    elseif theta isa Core.Const
        return theta.val
    elseif theta isa Number
        return theta
    else
        error("expect SSAValue or Number")
    end
end

function convert_to_zxd(root::YaoHIR.BlockIR)
    circ = ZXDiagram(root.nqubits)

    for gate in YaoHIR.leaves(root.circuit)
        @switch gate begin
            @case Gate(&Z, loc::Locations{Int})
                push_gate!(circ, Val(:Z), plain(loc), 1//1)
            @case Gate(&X, loc::Locations{Int})
                push_gate!(circ, Val(:X), plain(loc), 1//1)
            @case Gate(&H, loc::Locations{Int})
                push_gate!(circ, Val(:H), plain(loc))
            @case Gate(&S, loc::Locations{Int})
                push_gate!(circ, Val(:Z), plain(loc), 1//2)
            @case Gate(&T, loc::Locations{Int})
                push_gate!(circ, Val(:Z), plain(loc), 1//4)
            @case Gate(shift(theta), loc::Locations{Int})
                theta = unwrap_ssa_phase(theta, root.parent)
                push_gate!(circ, Val(:Z), plain(loc), (1/π)*theta)
            @case Gate(Rx(theta), loc::Locations{Int})
                theta = unwrap_ssa_phase(theta, root.parent)
                push_gate!(circ, Val(:X), plain(loc), (1/π)*theta)
            @case Gate(Ry(theta), loc::Locations{Int})
                theta = unwrap_ssa_phase(theta, root.parent)
                push_gate!(circ, Val(:X), plain(loc),  1//2)
                push_gate!(circ, Val(:Z), plain(loc),  (1/π) * theta)
                push_gate!(circ, Val(:X), plain(loc), -1//2)
            @case Gate(Rz(theta), loc::Locations{Int})
                theta = unwrap_ssa_phase(theta, root.parent)
                push_gate!(circ, Val(:Z), plain(loc), (1/π)*theta)
            @case Gate(AdjointOperation(&S), loc::Locations{Int})
                push_gate!(circ, Val(:Z), plain(loc), 3//2)
            @case Gate(AdjointOperation(&T), loc::Locations{Int})
                push_gate!(circ, Val(:Z), plain(loc), 7//4)
            @case Ctrl(Gate(&X, loc::Locations{Int}), ctrl::CtrlLocations{Int}) # CNOT
                push_gate!(circ, Val(:CNOT), plain(loc), plain(ctrl))
            @case Ctrl(Gate(&Z, loc::Locations{Int}), ctrl::CtrlLocations{Int}) # CZ
                push_gate!(circ, Val(:CZ), plain(loc), plain(ctrl))
            @case _
                error("$gate is not supported")
        end
    end
    return circ
end

function convert_to_chain(circ::ZXDiagram{TT, P}) where {TT, P}
    spider_seq = spider_sequence(circ)
    vs = spiders(circ)
    locs = Dict()
    nqubit = nqubits(circ)
    qc = []
    frontier_v = ones(TT, nqubit)

    while sum([frontier_v[i] <= length(spider_seq[i]) for i = 1:nqubit]) > 0
        for q = 1:nqubit
            if frontier_v[q] <= length(spider_seq[q])
                v = spider_seq[q][frontier_v[q]]
                nb = ZXCalculus.neighbors(circ, v)
                if length(nb) <= 2
                    θ = phase(circ, v) * π
                    if spider_type(circ, v) == ZXCalculus.SpiderType.Z
                        if phase(circ, v) == 1
                            push!(qc, Gate(Z, Locations(q)))
                        elseif phase(circ, v) == 1//2
                            push!(qc, Gate(S, Locations(q)))
                        elseif phase(circ, v) == 3//2
                            push!(qc, Gate(AdjointOperation(S), Locations(q)))
                        elseif phase(circ, v) == 1//4
                            push!(qc, Gate(T, Locations(q)))
                        elseif phase(circ, v) == 7//4
                            push!(qc, Gate(AdjointOperation(T), Locations(q)))
                        elseif phase(circ, v) != 0
                            if θ isa Phase
                                θ = θ.ex
                            end
                            push!(qc, Gate(shift(θ), Locations(q)))
                        end
                    elseif spider_type(circ, v) == ZXCalculus.SpiderType.X
                        if phase(circ, v) == 1
                            push!(qc, Gate(X, Locations(q)))
                        else phase(circ, v) != 0
                            if θ isa Phase
                                θ = θ.ex
                            end
                            push!(qc, Gate(Rx(θ), Locations(q)))
                        end
                    elseif spider_type(circ, v) == ZXCalculus.SpiderType.H
                        push!(qc, Gate(H, Locations(q)))
                    end

                    frontier_v[q] += 1
                elseif length(nb) == 3
                    v1 = nb[[qubit_loc(circ, u) != q for u in nb]][1]
                    if spider_type(circ, v1) == SpiderType.H
                        v1 = setdiff(ZXCalculus.neighbors(circ, v1), [v])[1]
                    end
                    if sum([findfirst(isequal(u), spider_seq[qubit_loc(circ, u)]) != frontier_v[qubit_loc(circ, u)] for u in [v, v1]]) == 0
                        if phase(circ, v) != 0
                            if spider_type(circ, v) == ZXCalculus.SpiderType.Z
                                if phase(circ, v) == 1
                                    push!(qc, Gate(Z, Locations(qubit_loc(circ, v))))
                                elseif phase(circ, v) == 1//2
                                    push!(qc, Gate(S, Locations(qubit_loc(circ, v))))
                                elseif phase(circ, v) == 3//2
                                    push!(qc, Gate(S', Locations(qubit_loc(circ, v))))
                                elseif phase(circ, v) == 1//4
                                    push!(qc, Gate(T, Locations(qubit_loc(circ, v))))
                                elseif phase(circ, v) == 7//4
                                    push!(qc, Gate(T', Locations(qubit_loc(circ, v))))
                                else
                                    θ = phase(circ, v)*π
                                    if θ isa Phase
                                        θ = θ.ex
                                    end
                                    push!(qc, Gate(shift(θ), Locations(qubit_loc(circ, v))))
                                end
                            else
                                if phase(circ, v) == 1
                                    push!(qc, Gate(X, Locations(qubit_loc(circ, v))))
                                else
                                    θ = phase(circ, v)*π
                                    if θ isa Phase
                                        θ = θ.ex
                                    end
                                    push!(qc, Gate(Rx(θ), Locations(qubit_loc(circ, v))))
                                end
                            end
                        end
                        if phase(circ, v1) != 0
                            if spider_type(circ, v1) == ZXCalculus.SpiderType.Z
                                if phase(circ, v1) == 1
                                    push!(qc, Gate(Z, Locations(qubit_loc(circ, v1))))
                                elseif phase(circ, v1) == 1//2
                                    push!(qc, Gate(S, Locations(qubit_loc(circ, v1))))
                                elseif phase(circ, v1) == 3//2
                                    push!(qc, Gate(S', Locations(qubit_loc(circ, v1))))
                                elseif phase(circ, v1) == 1//4
                                    push!(qc, Gate(T, Locations(qubit_loc(circ, v1))))
                                elseif phase(circ, v1) == 7//4
                                    push!(qc, Gate(T', Locations(qubit_loc(circ, v1))))
                                else
                                    θ = phase(circ, v1)*π
                                    if θ isa Phase
                                        θ = θ.ex
                                    end
                                    push!(qc, Gate(shift(θ), Locations(qubit_loc(circ, v1))))
                                end
                            else
                                if phase(circ, v1) == 1
                                    push!(qc, Gate(X, Locations(qubit_loc(circ, v1))))
                                else
                                    θ = phase(circ, v1)*π
                                    if θ isa Phase
                                        θ = θ.ex
                                    end
                                    push!(qc, Gate(Rx(θ), Locations(qubit_loc(circ, v1))))
                                end
                            end
                        end

                        if spider_type(circ, v) == spider_type(circ, v1) == ZXCalculus.SpiderType.Z
                            push!(qc, Ctrl(Gate(Z, Locations(qubit_loc(circ, v))), CtrlLocations(qubit_loc(circ, v1))))
                        elseif spider_type(circ, v) == ZXCalculus.SpiderType.Z
                            push!(qc, Ctrl(Gate(X, Locations(qubit_loc(circ, v1))), CtrlLocations(qubit_loc(circ, v))))
                        elseif spider_type(circ, v) == ZXCalculus.SpiderType.X
                            push!(qc, Ctrl(Gate(X, Locations(qubit_loc(circ, v))), CtrlLocations(qubit_loc(circ, v1))))
                        end
                        for u in [v, v1]
                            frontier_v[qubit_loc(circ, u)] += 1
                        end
                    end
                else
                    error("ZX-diagram without a circuit structure is not supported!")
                end
            end
        end
    end
    return Chain(qc...)
end


text:
function estimate_regressand_variance(obj::MultipleLinearRegression)
    return jcall(obj, "estimateRegressandVariance", jdouble, ())
end

function estimate_regression_parameters(obj::MultipleLinearRegression)
    return jcall(obj, "estimateRegressionParameters", Vector{jdouble}, ())
end

function estimate_regression_parameters_standard_errors(obj::MultipleLinearRegression)
    return jcall(obj, "estimateRegressionParametersStandardErrors", Vector{jdouble}, ())
end

function estimate_regression_parameters_variance(obj::MultipleLinearRegression)
    return jcall(obj, "estimateRegressionParametersVariance", Vector{Vector{jdouble}}, ())
end

function estimate_residuals(obj::MultipleLinearRegression)
    return jcall(obj, "estimateResiduals", Vector{jdouble}, ())
end


text:
module PolyOpt

using ClusteredLowRankSolver, BasesAndSamples, AbstractAlgebra

export min_f

function invariant_basis(x,y,z, d)
    # create a vector with a precise type
    v = [(x*y*z)^0]
    for deg=1:d, j=0:div(deg,3), i=0:div(deg-3j,2)
        push!(v, (x+y+z)^(deg-2i-3j) * (x*y+y*z+z*x)^i * (x*y*z)^j)
    end
    return v
end

function min_f(d)
    obj = Objective(0, Dict(), Dict(:M => 1))

    FF = RealField
    R, (x,y,z) = PolynomialRing(FF, ["x", "y", "z"])
    # The polynomial f:
    f =  x^4 + y^4 + z^4 - 4x*y*z + x + y + z

    # An invariant basis up to degree d:
    basis = invariant_basis(x, y, z, 2d)
    # For the sum-of-squares polynomials we have to
    # select elements of the basis based on the degree
    degrees = [total_degree(p) for p in basis]

    # generate samples and a good basis
    cheb_points = [vcat(sample_points_chebyshev(2d+k)...) for k=0:2]
    samples_grid = [[cheb_points[1][i+1], cheb_points[2][j+1], cheb_points[3][k+1]]
        for i=0:2d for j=0:2d+1 for k=0:2d+2]
    basis, samples = approximatefekete(basis, samples_grid)

    psd_dict = Dict()
    symmetry_weights = [[[R(1)]],
                        [[R((x-y)*(y-z)*(z-x))]],
                        [[1/sqrt(FF(2))*(2x-y-z),1/sqrt(FF(2))*(2y*z-x*z-x*y)],
                            [sqrt(FF(3)/FF(2))*(y-z),sqrt(FF(3)/FF(2))*(x*z-x*y)]]]
    for swi=1:length(symmetry_weights)
        rank = length(symmetry_weights[swi])
        # This has in general too many entries, so we remove the ones with too high degree.
        vecs = [kron(symmetry_weights[swi][r], basis) for r=1:rank]
        for r=1:rank
            # the length of the symmetric basis we tensored the symmetric weight with
            l = length(basis)
            # we keep the elements with degree at most d
            keep_idx = [i for i=1:length(vecs[r]) if total_degree(symmetry_weights[swi][r][div(i-1,l)+1]) + degrees[(i-1)%l+1] <= d]
            vecs[r] = vecs[r][keep_idx]
        end
        vecs = [v for v in vecs if length(v) > 0]
        if length(vecs) > 0
            psd_dict[Block((:trivariatesos,swi))] = LowRankMatPol([R(1) for r=1:length(vecs)], vecs)
        end
    end

    constr = Constraint(f, psd_dict, Dict(:M => R(1)), samples)
    pol_problem = LowRankPolProblem(true, obj, [constr])
    sdp = ClusteredLowRankSDP(pol_problem)

    solvesdp(sdp)
end

end


text:
using Pkg, DrWatson

@quickactivate "StatisticsWithJuliaPlutoNotebooks"

function lineSearch(inputFilename, outputFilename, keyword)
    infile  = open(inputFilename, "r")
    outfile = open(outputFilename,"w")

    for (index, line) in enumerate(split(read(infile, String), "\n"))
        if occursin(keyword, line)
            println(outfile, "$index: $line")
        end
    end
    close(infile)
    close(outfile)
end

lineSearch(datadir("earth.txt"), datadir("waterLines.txt"), "water") |> display


text:
using POMDPModels
using Test

let
    problem = InvertedPendulum()
    policy = RandomPolicy(problem)
    sim = RolloutSimulator(MersenneTwister(1), 20)

    simulate(sim, problem, policy)

    sv = convert_s(Array{Float64}, (0.5, 0.25), problem)
    @test sv == [0.5, 0.25]
    s = convert_s(Tuple{Float64,Float64}, sv, problem)
    @test s == (0.5, 0.25)
end


text:
abstract Layout

type Row <: Layout
    chidren::Vector{Layout}
end

type Col <: Layout
    chidren::Vector{Layout}
end

type Element <: Layout
    template::Symbol
    dataset::Int
end

abstract Token

type OpenRow <: Token end
type CloseRow <: Token end
type SplitCol <: Token end
type EofToken <: Token end
type StrToken <: Token
    value::String
end
type NumToken <: Token
    value::Int
end

type _tokenizer_next{T}
    this::T
end

type Tokenizer
    code::String
    i::Int

    next::_tokenizer_next{Tokenizer}
    Tokenizer(code) = begin
        this = new(code, 1)
        this.next = _tokenizer_next{Tokenizer}(this)
        this
    end
end

function (_this::_tokenizer_next)()
    code = _this.this.code

    while _this.this.i <= length(code)
        i = _this.this.i
        _this.this.i += 1
        if 'a' <= code[i] <= 'z'
            while _this.this.i <= length(code) && 'a' <= code[_this.this.i] <= 'z'
                _this.this.i += 1
            end
            return StrToken(code[i:_this.this.i-1])
        elseif '1' <= code[i] <= '9'
            return NumToken(Int(code[i]) - Int('0'))
        elseif code[i] == '['
            return OpenRow()
        elseif code[i] == ']'
            return CloseRow()
        elseif code[i] == '|'
            return SplitCol()
        end # ignore any other characters, including whitespaces, uppercase letters etc.
    end

    EofToken()
end

function parse_layout(tokens::Tokenizer)
    chidren = []

    while true
        next = tokens.next()

        if isa(next, CloseRow) || isa(next, EofToken)
            if length(chidren) == 1 && isa(chidren[], Element)
                return chidren[]
            elseif SplitCol() in chidren
                result, this = Layout[], Layout[]

                for c in chidren
                    if isa(c, SplitCol)
                        if isempty(this)
                            throw(ParseError("unexpected SplitCol()"))
                        elseif length(this) == 1
                            push!(result, this[])
                        else
                            push!(result, Col(this))
                        end

                        this = Layout[]
                    else
                        push!(this, c)
                    end
                end

                if length(this) == 1
                    push!(result, this[])
                elseif length(this) > 1
                    push!(result, Col(this))
                end

                return Row(result)
            else
                return Col(chidren)
            end
        elseif isa(next, StrToken)
            push!(chidren, Element(next.value, 0))
        elseif isa(next, NumToken)
            if isa(chidren[end], Element)
                chidren[end].dataset = next.value
            else
                throw(ParseError("unexpected $next"))
            end
        elseif isa(next, OpenRow)
            push!(chidren, parse_layout(tokens))
        elseif isa(next, SplitCol)
            push!(chidren, next)
        end
    end
end