text:
#solves the EulerLagrane equations setup in DiscreteShape and DiscreteShapeBound.
function solve_EulerLagrange(ρ, I, iterations, p, q, R, λ, E, Z)
    #initilize common variables
    M, N = size(E)
    δp = zeros(Complex{Float64},axes(E))
    δq = zeros(Complex{Float64},axes(E))
    pq = zeros(Complex{Float64}, size(p))
    dRδp = zeros(Complex{Float64}, size(R))
    dRδq = zeros(Complex{Float64}, size(R))
    C = zeros(Complex{Float64}, size(R))
    w = centered(0.25*[0.0 1.0 0.0;1.0 0.0 1.0;0.0 1.0 0.0])
    wx, wy = setup_transform_values(M, N)
    for i = 1:iterations
        #calculate second derivatives of p and q
        δp = imfilter(p, reflect(w), "replicate")
        δq = imfilter(q, reflect(w), "replicate")

        for i in CartesianIndices(R)
            #calculate reflectance map
            R[i] = (ρ * (-I[1] * p[i] - I[2] * q[i] + I[3])) / sqrt(1 + p[i]^2
                + q[i]^2)
            #calculate partial derivatives of p and q in respect to p and q
            pq[i] = (1 + p[i]^2 + q[i]^2)
            dRδp[i] = ((-ρ * I[1]) / sqrt(pq[i])) + ((-I[1] * ρ) * p[i] - I[2]
                * ρ * q[i] + I[3] * ρ) * (-1 * p[i] * (pq[i]^(-3/2)))

            dRδq[i] = ((-ρ * I[2]) / sqrt(pq[i])) + (-I[1] * ρ * p[i] - I[2]
                * ρ * q[i] + I[3] * ρ) * (-1 * q[i] * (pq[i]^(-3/2)))
            #calculate new surface normals
            p[i] = δp[i] + (1 / (4 * λ)) * (E[i] - R[i]) * dRδp[i]
            q[i] = δq[i] + (1 / (4 * λ)) * (E[i] - R[i]) * dRδq[i]
        end

        #calculate Transform of surface normals
        fft!(p)
        fft!(q)

        for i in CartesianIndices(C)
            #compute C and prepare to update p and q
            C[i] = (-1im * wx[i] * p[i] + 1im * wy[i] * q[i]) / (2*π*(wx[i]^2 + wy[i]^2 + eps()))
            p[i] = 1im * wx[i] * C[i]
            q[i] = 1im * wy[i] * C[i]
        end

        #update p and q
        ifft!(p)
        ifft!(q)

        #recover Z by taking inverse transform of C
        ifft!(C)

        #take abs of C to get Z
        for i in CartesianIndices(Z)
            Z[i] = abs(C[i])
        end
    end
    return Z, p, q
end

# Set up values for fourier transform for Frankot method
function setup_transform_values(N, M)
    wx = [j for i in 1:M, j in -0.5:(1/(N-1)):0.5]
    wy = [j for j in -0.5:(1/(M-1)):0.5, i in 1:N]
    wx = ifftshift(wx)
    wy = ifftshift(wy)
    return wx, wy
end

# Set up values for fourier transform for Pentland method
function setup_transform_values_pentland(M, N)
    wx = zeros(M, N)
    wy = zeros(M, N)
    for i in CartesianIndices(wx)
        wx[i] = (2 * π * i[2]) / M
        wy[i] = (2 * π * (N - i[1] + 1)) / N
    end
    return wx, wy
end

# Sets up xy grid as array
function setup_xy(dim)
    x = zeros(dim , dim)
    y = zeros(dim , dim)
    xyrange = ceil(Int, -dim/2):floor(Int, (dim/2))
    for i in CartesianIndices(x)
        x[i] = xyrange[i[2]]
        y[i] = -xyrange[i[1]]
    end
    return x,y
end

function gen_mask(p::AbstractArray, q::AbstractArray,σ::Real)
    pq, pp = imgradients(p, KernelFactors.prewitt)
    qq, qp = imgradients(q, KernelFactors.prewitt)
    # Violation of Schwarz integrability condition + large jump in one derivative without a corrisponding jump in the other
    mask = abs.(pq - qp) + abs.(pp - qq)
    t = mean(mask)+σ*std(mask)
    for i in CartesianIndices(mask)
          mask[i] = mask[i] <= t ? 1 : 0
    end
    stop = false
    while !stop
        stop = true
        for i in CartesianIndices(mask)
            if mask[i] != 0 && i[1] != 1 && i[2] != 1 && i[1] != first(size(p)) && i[2] != last(size(p))
                count = 0
                count += mask[i[1]+1,i[2]] == 0 ? 1 : 0
                count += mask[i[1]-1,i[2]] == 0 ? 1 : 0
                count += mask[i[1],i[2]+1] == 0 ? 1 : 0
                count += mask[i[1],i[2]-1] == 0 ? 1 : 0
                if count > 2
                    mask[i] = 0
                    stop = false
                end
            end
        end
    end
    splitMask = copy(mask)
    while findfirst(x -> x==1 , splitMask) != nothing
        # @show count(x -> x == 1, splitMask)
        seed = [findfirst(x -> x==1 , splitMask)]
        val = maximum(splitMask) + 1
        while length(seed) > 0
            i = 1
            # @show length(seed)
            if splitMask[seed[i]] == 1
                splitMask[seed[i]] = val
                if seed[i][1] != 1
                    if splitMask[seed[i][1]-1, seed[i][2]] == 1
                        push!(seed, CartesianIndex(seed[i][1]-1, seed[i][2]))
                    end
                end
                if seed[i][2] != last(size(mask))
                    if splitMask[seed[i][1], seed[i][2]+1] == 1
                        push!(seed, CartesianIndex(seed[i][1], seed[i][2]+1))
                    end
                end
                if seed[i][1] != first(size(mask))
                    if splitMask[seed[i][1]+1, seed[i][2]] == 1
                        push!(seed, CartesianIndex(seed[i][1]+1, seed[i][2]))
                    end
                end
                if seed[i][2] != 1
                    if splitMask[seed[i][1], seed[i][2]-1] == 1
                        push!(seed, CartesianIndex(seed[i][1], seed[i][2]-1))
                    end
                end
            end
            deleteat!(seed,1)
            # @show length(seed)
        end
    end
    layerMask = zeros(Float64, size(mask)..., Int(maximum(splitMask)))
    layerMask[:,:,1] = mask
    for i in CartesianIndices(splitMask)
        if splitMask[i] != 0
            layerMask[i,Int(splitMask[i])] = 1
        end
    end
    return layerMask
end


text:
using PhysicsCodesLabXD
using Documenter

DocMeta.setdocmeta!(PhysicsCodesLabXD, :DocTestSetup, :(using PhysicsCodesLabXD); recursive=true)

makedocs(;
    modules=[PhysicsCodesLabXD],
    authors="PhysicsCodesLab",
    repo="https://github.com/PhysicsCodesLab/PhysicsCodesLabXD.jl/blob/{commit}{path}#{line}",
    sitename="PhysicsCodesLabXD.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://PhysicsCodesLab.github.io/PhysicsCodesLabXD.jl",
        assets=String[],
        mathengine = MathJax()
    ),
    pages=[
        "Home" => "index.md",
        "Manual" => ["Manual/introduction.md","Manual/category.md"],
        "Library" => [],
        "Index" => ["Index/index.md"]
    ],
)

deploydocs(;
    repo="github.com/PhysicsCodesLab/PhysicsCodesLabXD.jl",
    devbranch="main",
)

ext: jl


text:
using Graphs
using AtomicLevels

type LevelVertex
    config::Config
    term::Term
    calc_props::Dict
    exp_eng
    wfn
end
import Graphs.attributes
function attributes(l::LevelVertex, g::AbstractGraph)
    Dict{UTF8String,Any}("label" => "$(l.config) $(l.term)")
end

LevelVertex(config::Config, term::Term,
      exp_eng = Inf;
      kwargs...) = LevelVertex(config, term, Dict(kwargs), exp_eng, nothing)

level_graph_init = () -> graph(Vector{LevelVertex}(),Vector{Edge{LevelVertex}}())
level_graph = level_graph_init()

import Base.print, Base.show
function print(io::IO, lg::typeof(level_graph_init()))
    map(reverse(topological_sort_by_dfs(level_graph))) do l
        println(l.config, " ", l.term)
    end
end

function add_level!(dep::Union{LevelVertex,Void}, args...; kwargs...)
    s = LevelVertex(args...; kwargs...)
    global level_graph
    add_vertex!(level_graph, s)
    dep != nothing && add_edge!(level_graph, s, dep)
    s
end
add_level!(args...; kwargs...) = add_level!(nothing, args...; kwargs...)

function add_levels!(dep::Union{LevelVertex,Void}, config::Config, args...; kwargs...)
    ts = terms(config)
    for i in eachindex(ts)
        a = copy(args)
        if length(a) > 0 && length(a[1]) > 1
            a = (a[1][i],a[2:end]...)
        end
        dep = add_level!(dep, config, ts[i], a...; kwargs...)
    end
    dep
end
add_levels!(args...; kwargs...) = add_levels!(nothing, args...; kwargs...)

function clear_levels!()
    global level_graph
    level_graph = level_graph_init()
end

# Used to perform calculation using MPI codes if that is
# requested. However, certain level calculations may be explicitly
# marked as unsuitable for MPI codes (e.g. gst of He II), so these are
# masked.
function maybe_mpi_run(g::Function, d::Dict)
    if get(d, :nompi, false)
        delete!(d, :nompi)
        delete!(d, :mpi_np)
        g(d)
    elseif (np = get(d, :mpi_np, 1)) > 1
        delete!(d, :mpi_np)
        init_mpi(np, joinpath(homedir(), "tmp_mpi")) do
            g(d)
        end
    else
        g(d)
    end
end

function get_and_delete!(d::Dict, sym)
    v = d[sym]
    delete!(d, sym)
    v
end

function get_and_delete!(d::Dict, args...)
    map(args) do a
        get_and_delete!(d, a)
    end
end

get_wfn(wfn::Union{AbstractString,Void}) = wfn
get_wfn(wfn::Vector) = first(wfn)

function calc_level(f::Function, level::LevelVertex, wfn; kwargs...)
    wfn = maybe_mpi_run((d) -> f(level, d, get_wfn(wfn)),
                        merge(Dict(kwargs), level.calc_props))
    println(wfn)
    assert(isfile(wfn))
    wfn
end

# Loop over all levels in topological order such that the output of
# one calculation can be used as input for the next, that has
# requested that.
function calc_levels(f::Function; kwargs...)
    global level_graph

    map(reverse(topological_sort_by_dfs(level_graph))) do l
        e = out_edges(l, level_graph)
        if length(e) > 0
            wfn = first(e).target.wfn
        else
            wfn = nothing
        end

        l.wfn = calc_level(f, l, wfn; kwargs...)
    end
end

function calc_levels_atsp(Z; kwargs...)
    calc_levels(;kwargs...) do l,d,wfn0
        ncorr = get_and_delete!(d, :ncorr)
        hf_mchf_bp(l.config, l.term, Z, ncorr, wfn0; overwrite = false, d...)[2]
    end
end

function plot_level_graph(filename)
    stdin, proc = open(`neato -T$(splitext(filename)[2][2:end]) -o "$filename"`, "w")
    to_dot(AST.level_graph, stdin)
    close(stdin)
    filename
end

export add_level!, add_levels!, clear_levels!, calc_levels, calc_levels_atsp, plot_level_graph


text:
abstract type SquareRootFilter end

struct CholeskyFactorization <: SquareRootFilter end
CholeskySqrt() = CholeskyFactorization()

struct SVDFactorization <: SquareRootFilter end
SVDSqrt() = SVDFactorization()

struct ModifiedCholeskyFactorization <: SquareRootFilter end
CholeskyModSqrt() = ModifiedCholeskyFactorization()


text:
using Compose, Gadfly

set_default_plot_size(6inch, 3inch)

plot(x=[0.768448, 0.940515, 0.673959, 0.395453, 0.313244, 0.662555, 0.586022, 0.0521332, 0.26864, 0.108871, 0.163666, 0.473017, 0.865412, 0.617492, 0.285698, 0.463847, 0.275819, 0.446568, 0.582318, 0.255981],
     y=[0.70586, 0.291978, 0.281066, 0.792931, 0.20923, 0.918165, 0.614255, 0.802665, 0.555668, 0.940782, 0.48, 0.790201, 0.356221, 0.900925, 0.529253, 0.031831, 0.900681, 0.940299, 0.621379, 0.348173],
     color=repeat(["A","B"], inner=10),
     Guide.colorkey(title="Species", labels=["Name1","Name2"], pos=[0.85w,-0.3h]))


text:
module PseudoArcLengthContinuation
	using Plots, JLD2, Printf, Dates, LinearMaps, BlockArrays
	using Setfield: setproperties, @set, Lens, get, set, @lens
	using Parameters: @with_kw, @unpack
	using DocStringExtensions

	include("BorderedArrays.jl")
	include("LinearSolver.jl")
	include("EigSolver.jl")
	include("LinearBorderSolver.jl")
	include("Preconditioner.jl")
	include("Newton.jl")
	include("Continuation.jl")
	include("Bifurcations.jl")
	include("Predictor.jl")
	include("NormalForms.jl")
	include("DeflationOperator.jl")
	include("BorderedProblem.jl")
	include("Plotting.jl")
	include("Utils.jl")
	include("codim2/codim2.jl")
	include("codim2/FoldCont.jl")
	include("codim2/HopfCont.jl")

	include("periodicorbit/PeriodicOrbits.jl")
	include("periodicorbit/PeriodicOrbitUtils.jl")
	include("periodicorbit/Flow.jl")
	include("periodicorbit/StandardShooting.jl")
	include("periodicorbit/PoincareShooting.jl")
	include("periodicorbit/PeriodicOrbitFD.jl")
	include("periodicorbit/FloquetQaD.jl")

	# linear solvers
	export	DefaultLS, GMRESIterativeSolvers, GMRESIterativeSolvers!, GMRESKrylovKit,
			DefaultEig, EigArpack, EigIterativeSolvers, EigKrylovKit, EigArnoldiMethod, geteigenvector, AbstractEigenSolver

	# bordered nonlinear problems
	export	BorderedProblem, JacobianBorderedProblem, LinearSolverBorderedProblem

	# preconditioner based on deflation
	export	PrecPartialSchurKrylovKit, PrecPartialSchurArnoldiMethod

	# bordered linear problems
	export	MatrixBLS, BorderingBLS, MatrixFreeBLS, BorderedArray

	# nonlinear deflation
	export	DeflationOperator, DeflatedProblem, DeflatedLinearSolver, scalardM

	# predictors for continuation
	export	SecantPred, BorderedPred, NaturalPred

	# newton methods
	export	NewtonPar, newton, newtonDeflated, newtonPALC, newtonFold, newtonHopf, newtonBordered

	# continuation methods
	export	ContinuationPar, ContResult, continuation, continuation!, continuationFold, continuationHopf, continuationPOTrap, continuationBordered

	# iterators for continuation
	export	PALCIterable, iterate, PALCStateVariables, solution, getx, getp

	# codim2 Fold continuation
	export	FoldPoint, FoldProblemMinimallyAugmented, FoldLinearSolveMinAug, foldPoint

	# codim2 Hopf continuation
	export	HopfPoint, HopfProblemMinimallyAugmented, HopfLinearSolveMinAug

	# normal form
	export computeNormalForm

	# Periodic orbit computation based on Trapeze method
	export	PeriodicOrbitTrapProblem, continuationPOTrap, continuationPOTrapBPFromPO

	# Periodic orbit computation based on Shooting
	export	Flow, ShootingProblem, PoincareShootingProblem, continuationPOShooting, getPeriod, AbstractShootingProblem, extractPeriodShooting

	# Floquet multipliers computation
	export	FloquetQaDTrap, FloquetQaDShooting

	# guess for periodic orbit from Hopf bifurcation point
	export	guessFromHopf
end


text:
using Plots
order = 2 # 0 for local constant, 1 for local linear, 2 for local quadratic
k = 3 # number of regressors
srand(1) # set seed to enable testing by fact check
bandwidth = 0.07
n = 100000
neval = 100
x = rand(n,k)*pi*2.0
xeval = [pi*ones(neval,k-1) linspace(pi/2.,pi*1.5,neval)]
y = cos.(sum(x,2)) + cos.(2. * sum(x,2)) + 0.5*randn(n,1)
ytrue = cos.(sum(xeval,2)) + cos.(2. * sum(xeval,2))
weights = kernelweights(x, xeval, bandwidth, true, "gaussian", 200)
yhat = npreg(y, x, xeval, weights, order)
plot(xeval[:,k], [yhat ytrue])


text:
module Rand

export cryptographic_rand, create_block_cipher_key

include("gen/libcrypto_common.jl")
include("gen/rand.h.jl")

cryptographic_rand(len) = cryptographic_rand(zeros(Uint8, len))
function cryptographic_rand(a :: Array{Uint8})
    RAND_bytes(a, length(a))
    a
end

create_block_cipher_key(size) = cryptographic_rand(size)

end # module Rand


text:
#=

    This is a port of the SPPL model disease-infection.pynb
    (The arg max value of the posterior is 330 in the SPPL model, but this is based
     on a rather crude interval. Using a less crude, but much slower, interval give 333.)

    This model:
    parameters       mean       std   naive_se      mcse         ess      rhat   ess_per_sec
        Symbol    Float64   Float64    Float64   Float64     Float64   Float64       Float64

    num_people   671.4048   61.1521     0.6115    1.3344   2135.0113    1.0012      511.6250
       num_met   335.7222   27.9166     0.2792    0.5434   2402.8246    1.0005      575.8027

    Quantiles
    parameters       2.5%      25.0%      50.0%      75.0%      97.5%
        Symbol    Float64    Float64    Float64    Float64    Float64

    num_people   560.0000   628.0000   671.0000   713.0000   793.0000
       num_met   284.0000   316.0000   335.0000   354.0000   391.0250

=#

using Turing, StatsPlots, Distributions, StatsBase
using CSV
include("jl_utils.jl")

@model function disease_infection(num_infected=100)
    num_people ~ DiscreteUniform(500,1500) # values=range(500, 1500, 20))

    # switch (num_people) cases (n in range(500, 1500, 20)):
    num_met ~ Binomial(num_people, 0.5)
    # condition ((low[n] <= num_met) < high[n])
    # switch (num_met) cases (m in range(low[n], high[n])):
    num_infected ~ Binomial(num_met, 0.3)
end

# low  = {n : int(scipy.stats.binom.ppf(.01, n=n, p=.5)) for n in range(500, 1500)}
# high = {n : int(scipy.stats.binom.ppf(.99, n=n, p=.5)) for n in range(500, 1500)}

num_infected = 100
model = disease_infection(num_infected)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, HMC(0.1,6), 1_000)
# chns = sample(model, NUTS(), 1_000)

display(chns)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
mutable struct IterationChannel{D} <: AbstractChannel{D}
    data::Dict{Int,D}
    cond_take::Condition
    IterationChannel(data::Dict{Int,D}) where D = new{D}(data, Condition())
end

function put!(channel::IterationChannel, t, x)
    channel.data[t] = copy(x)
    notify(channel.cond_take)
    return channel
end

function take!(channel::IterationChannel, t)
    x = fetch(channel, t)
    delete!(channel.data, t)
    return x
end

isready(channel::IterationChannel) = length(channel.data) > 1
isready(channel::IterationChannel, t) = haskey(channel.data, t)

function fetch(channel::IterationChannel, t)
    wait(channel, t)
    return channel.data[t]
end

function wait(channel::IterationChannel, t)
    while !isready(channel, t)
        wait(channel.cond_take)
    end
end

mutable struct RunningAverageChannel{D} <: AbstractChannel{D}
    average::D
    data::Vector{D}
    buffer::Dict{Int,D}
    cond_put::Condition
    RunningAverageChannel(average::D, data::Vector{D}) where D = new{D}(average, data, Dict{Int,D}(), Condition())
end

function take!(channel::RunningAverageChannel, i::Integer)
    channel.buffer[i] = copy(channel.data[i])
end

function put!(channel::RunningAverageChannel, i::Integer, π::AbstractFloat)
    channel.average -= π*channel.buffer[i]
    channel.average += π*channel.data[i]
    delete!(channel.buffer, i)
    notify(channel.cond_put)
    return channel
end

function put!(channel::RunningAverageChannel{D}, i::Integer, x::D, π::AbstractFloat) where D
    channel.average -= π*channel.buffer[i]
    channel.average += π*x
    channel.data[i] = copy(x)
    delete!(channel.buffer, i)
    notify(channel.cond_put)
    return channel
end

isready(channel::RunningAverageChannel) = length(channel.buffer) == 0

function fetch(channel::RunningAverageChannel)
    wait(channel)
    return channel.average
end

function fetch(channel::RunningAverageChannel, i::Integer)
    return channel.data[i]
end

function wait(channel::RunningAverageChannel)
    while !isready(channel)
        wait(channel.cond_put)
    end
end

IteratedValue{T <: AbstractFloat} = RemoteChannel{IterationChannel{T}}
RunningAverage{D} = RemoteChannel{RunningAverageChannel{D}}
Decisions{A <: AbstractArray} = RemoteChannel{IterationChannel{A}}


text:
@testset "Capture rate" begin

    Q = range(-20, stop=35, length=4250)

    # potential
    # pot1
    pot1 = Potential(); pot1.name = "pot1"; pot1.nev = 30
    pot1.Q0 = 0; pot1.E0 = 1.5
    pot1.func_type = "harmonic"
    params = Dict()
    params["hw"] = 0.03
    pot1.params = params
    pot1.Q = Q
    fit_pot!(pot1)

    # pot2
    pot2 = Potential(); pot2.name = "pot2"; pot2.nev = 50
    pot2.Q0 = 2; pot2.E0 = 0.6
    pot2.func_type = "harmonic"
    params = Dict()
    params["hw"] = 0.03
    pot2.params = params
    pot2.Q = Q
    fit_pot!(pot2)

    solve_pot!(pot1)
    solve_pot!(pot2)

    # conf_coord
    Tmin = 10; Tmax = 800; NT = 100
    Volume = 1.28463E-21
    cut_off = 0.25; σ = 0.0075
    temperature = range(Tmin, stop=Tmax, length=NT)

    W = 0.205 # e-ph coupling
    g = 1 # degeneracy

    # build a configuration coordinate for the electron capture
    cc = conf_coord(pot1, pot2)
    cc.W = W; cc.g = g

    # calc_overlap!
    calc_overlap!(cc; cut_off = cut_off, σ = σ)

    using LinearAlgebra
    # zero-overlap between misaligned eigenvalues
    @test all(diag(cc.δ_matrix) .== 0)

    # calc_capt_coeff!
    calc_capt_coeff!(cc, Volume, temperature)
    @test all(diag(cc.partial_capt_coeff[:, :, 1]) .== 0)
end


text:
# This file was generated, do not modify it. # hide
pred_rfr_tm = MLJ.predict(rfr_tm, rows=test);
rms(pred_rfr_tm, y[test])


text:
# ============================================================
# Period Functions
# ============================================================

function bc_period!(ctr::AbstractVector{T}, ng = 1) where {T<:AbstractControlVolume}
    nx = length(ctr) - 2 * ng

    for i = 1:ng
        copy_ctr!(ctr[1-i], ctr[nx+1-i])
        copy_ctr!(ctr[nx+i], ctr[i])
    end

    return nothing
end

function bc_period!(ctr::AbstractMatrix{T}, ng; dirc) where {T<:AbstractControlVolume}
    if dirc == :x
        nx = size(ctr, 1) - 2 * ng
        for j in axes(ctr, 2), i = 1:ng
            copy_ctr!(ctr[1-i, j], ctr[nx-i+1, j])
        end
    elseif dirc == :y
        ny = size(ctr, 2) - 2 * ng
        for j = 1:ng, i in axes(ctr, 1)
            copy_ctr!(ctr[i, 1-j], ctr[i, ny-j+1])
        end
    end

    return nothing
end


text:
module IntegralVerification

using Statistics
using Random

export testIntegral, testExpressions

function doStatTest(testvalues, precision; verbose = true)
    output = Statistics.mean(testvalues),Statistics.std(testvalues)
    if verbose
        println("mean = ",output[1])
        println("std  = ", output[2])
    end
    if abs(output[1]) < precision && abs(output[2]) < precision
        return true
    end
    return false
end

function generateRandomArgs(randomVariables)
    vars = zeros(ComplexF64, length(randomVariables))
    for (i,f) in enumerate(randomVariables)
        vars[i] = f()
    end
    return vars
end

function testIntegral(integrator, integrand, analytic, randomVariables, Ntests,  precision = 10. ^ (-7))

    numeric(args...) = integrator(x -> integrand(x, args...))

    return testExpressions(numeric, analytic, randomVariables, Ntests,precision)

end

function testExpressions(expr1, expr2, Ntests)
    return testExpressions(expr1,expr2,[randn],Ntests)
end

function testExpressions(expr1,expr2,randomVariables,Ntests,precision = 10. ^ (-7))

    if Ntests < 2
        throw(ArgumentError)
    end

    outvalues = zeros(ComplexF64, Ntests)
    for i in 1:Ntests
        args = generateRandomArgs(randomVariables)
        outvalues[i] = expr1(args...) - expr2(args...)
    end

    passed = doStatTest(outvalues, precision)

    if passed
        println("passed")
    else
        println("failed")
    end
    return passed
end

end # module


text:
# Try to open at least all downloaded datasets
@testset "Reading rasters" begin
    for f in remotefiles
        @testset "Reading $f" begin
            ga = GeoArrays.read(joinpath(testdatadir, f))
        end
    end
end

@testset "Reading rasters first band" begin
    for f in remotefiles
        @testset "Reading $f" begin
            ga = GeoArrays.read(joinpath(testdatadir, f), band=1)
            @test last(size(ga)) == 1
        end
    end
end

@testset "Reading rasters streaming" begin
    for f in remotefiles
        @testset "Reading $f streaming" begin
            ga = GeoArrays.read(joinpath(testdatadir, f), masked=false)
            ga[1, 1, 1]
        end
    end
end

@testset "Reading rasters streaming first band" begin
    for f in remotefiles
        @testset "Reading $f streaming" begin
            ga = GeoArrays.read(joinpath(testdatadir, f), masked=false, band=1)
            @test last(size(ga)) == 1
            ga[1, 1, 1]
        end
    end
end

@testset "Reading rasters and writing" begin
    for (i, f) in enumerate(remotefiles)
        @testset "Reading and writing $f" begin
            ga = GeoArrays.read(joinpath(testdatadir, f))
            GeoArrays.write!(joinpath(tempdir(), "test_$i.tif"), ga)
            ga_copy = GeoArrays.read(joinpath(tempdir(), "test_$i.tif"))
            @test ga[1, 1, 1] === ga_copy[1, 1, 1]
            @test ga[end, end, end] === ga_copy[end, end, end]
        end
    end
end

@testset "Read second band" begin
    fn = joinpath(testdatadir, remotefiles[end-1])
    GeoArrays.read(fn, band=2)
end

@testset "Writing rasters" begin

    @testset "Simplest version" begin
        ga = GeoArray(rand(100, 200, 3))
        fn = GeoArrays.write!(joinpath(tempdir(), "test.tif"), ga)
        GeoArrays.read(fn)
        fn = GeoArrays.write!(joinpath(tempdir(), "test.img"), ga)
        GeoArrays.read(fn)
    end
    @testset "Nodata" begin
        ga = GeoArray(Array{Union{Missing,Int32}}(rand(1:10, 100, 200, 3)))
        fn = GeoArrays.write!(joinpath(tempdir(), "test_nodata.tif"), ga, 1)
        GeoArrays.read(fn)
        fn = GeoArrays.write!(joinpath(tempdir(), "test_nodata.img"), ga, 1)
        GeoArrays.read(fn)
        ga = GeoArray(fill(0, (100, 200, 3)))
        fn = GeoArrays.write(joinpath(tempdir(), "test_nodata_hardcoded.tiff"), ga; nodata=0)
        ga = GeoArrays.read(fn)
        @test all(ismissing.(ga))
    end
    @testset "COG" begin
        ga = GeoArray(Array{Union{Missing,Int32}}(rand(1:10, 2048, 2048, 3)))
        fn = GeoArrays.write(joinpath(testdatadir, "test_cog.tif"), ga; nodata=-1, shortname="COG", options=Dict("COMPRESSION" => "ZSTD"))
        GeoArrays.read(fn)
        ga = GeoArray(Array{Union{Missing,Float32}}(rand(1:10, 2048, 2048, 3)))
        fn = GeoArrays.write(joinpath(testdatadir, "test_cogf.tif"), ga; nodata=Inf, shortname="COG", options=Dict("COMPRESSION" => "ZSTD"))
        GeoArrays.read(fn)
    end
    @testset "Kwargs" begin
        ga = GeoArray(rand(100, 200, 3))
        fn = GeoArrays.write(joinpath(tempdir(), "test.tif"), ga; shortname="COG", nodata=1.0, options=Dict("compression" => "deflate"))
        GeoArrays.read(fn)
    end

end


text:
function get_stop_model(env::UrbanEnv, r::SVector{2, LaneTag})
    route = [env.roadway[tag] for tag in find_route(env, r)]
    intersection_entrances = get_start_lanes(env.roadway)
    if !(route[1] ∈ intersection_entrances)
        intersection = Lane{Float64}[]
        intersection_exits = Lane{Float64}[]
    else
        intersection_exits = get_exit_lanes(env.roadway)
        intersection=Lane{Float64}[route[1], route[2]]
    end
    navigator = RouteFollowingIDM(route=route, a_max=3., σ=1.)
    intersection_driver = StopIntersectionDriver(navigator= navigator,
                                                intersection=intersection,
                                                intersection_entrances = intersection_entrances,
                                                intersection_exits = intersection_exits,
                                                stop_delta=maximum(env.params.crosswalk_width),
                                                accel_tol=0.,
                                                priorities = env.priorities)
    crosswalk_drivers = Vector{CrosswalkDriver}(length(env.crosswalks))
    for i=1:length(env.crosswalks)
        cw_conflict_lanes = get_conflict_lanes(env.crosswalks[i], env.roadway)
        crosswalk_drivers[i] = CrosswalkDriver(navigator = navigator,
                                crosswalk = env.crosswalks[i],
                                conflict_lanes = cw_conflict_lanes,
                                intersection_entrances = intersection_entrances,
                                yield=!isempty(intersect(cw_conflict_lanes, route))
                                )
        # println(" yield to cw ", i, " ", crosswalk_drivers[i].yield)
    end
    return UrbanDriver(navigator=navigator,
                    intersection_driver=intersection_driver,
                    crosswalk_drivers=crosswalk_drivers
                    )
end

function get_ttc_model(env::UrbanEnv, r::SVector{2, LaneTag})
    route = [env.roadway[tag] for tag in find_route(env, r)]
    intersection_entrances = get_start_lanes(env.roadway)
    if !(route[1] ∈ intersection_entrances)
        intersection = Lane[]
        intersection_exits = Lane[]
    else
        intersection_exits = get_exit_lanes(env.roadway)
        intersection=Lane[route[1], route[2]]
    end
    navigator = RouteFollowingIDM(route=route, a_max=1., σ=1., v_des=6.0)
    intersection_driver = TTCIntersectionDriver(navigator = navigator,
                                            intersection = intersection,
                                            intersection_pos = VecSE2(env.params.inter_x,
                                                                      env.params.inter_y),
                                            stop_delta = maximum(env.params.crosswalk_width),
                                            accel_tol = 0.,
                                            priorities = env.priorities,
                                            ttc_threshold = (env.params.x_max - env.params.inter_x)/env.params.speed_limit
                                            )
    crosswalk_drivers = Vector{CrosswalkDriver}(undef, length(env.crosswalks))
    for i=1:length(env.crosswalks)
        cw_conflict_lanes = get_conflict_lanes(env.crosswalks[i], env.roadway)
        crosswalk_drivers[i] = CrosswalkDriver(navigator = navigator,
                                crosswalk = env.crosswalks[i],
                                conflict_lanes = cw_conflict_lanes,
                                intersection_entrances = intersection_entrances,
                                yield=!isempty(intersect(cw_conflict_lanes, route))
                                )
        # println(" yield to cw ", i, " ", crosswalk_drivers[i].yield)
    end
    return UrbanDriver(navigator=navigator,
                    intersection_driver=intersection_driver,
                    crosswalk_drivers=crosswalk_drivers
                    )
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
# This is unfortunately 10x slower than libknet8, 2x slower than CUDA.jl

"""
    cudnnReduceTensor(x; dims, op, compType, nanOpt, indices, alpha)
    cudnnReduceTensor(x, d::cudnnReduceTensorDescriptor; dims, indices, alpha)
    cudnnReduceTensor!(y, x; op, compType, nanOpt, indices, alpha, beta)
    cudnnReduceTensor!(y, x, d::cudnnReduceTensorDescriptor; indices, alpha, beta)

Return the result of the specified reduction operation applied to `x`.  Optionally `y` holds
the result and `d` specifies the operation.  Each dimension of the output tensor `y` must
match the corresponding dimension of the input tensor `x` or must be equal to 1. The
dimensions equal to 1 indicate the dimensions of `x` to be reduced.  Keyword arguments:

* `dims = ntuple(i->1,ndims(x))`: specifies the shape of the output when `y` is not given
* `indices = nothing`: previously allocated space for writing indices which can be generated for min and max ops only, can be a `CuArray` of `UInt8`, `UInt16`, `UInt32` or `UInt64`
* `alpha=1, beta=0` are used for scaling, i.e. `y .= alpha*op.(x1) .+ beta*y`

Keyword arguments that can be used when `reduceTensorDesc` is not specified:
* `op = CUDNN_REDUCE_TENSOR_ADD`: Reduction operation, ADD can be replaced with MUL, MIN, MAX, AMAX, AVG, NORM1, NORM2, MUL_NO_ZEROS
* `compType = (eltype(x) <: Float64 ? Float64 : Float32)`: Computation datatype
* `nanOpt = CUDNN_NOT_PROPAGATE_NAN`: NAN propagation policy, the other option is `CUDNN_PROPAGATE_NAN`
"""
cudnnReduceTensor, cudnnReduceTensor!

# Public methods
cudnnReduceTensor(x; o...)     = cudnnReduceTensorWithDefaults(x; o...)
cudnnReduceTensor!(y, x; o...) = cudnnReduceTensorWithDefaults(x; y, o...)
cudnnReduceTensor(x, d::cudnnReduceTensorDescriptor; o...)     = cudnnReduceTensorWithDefaults(x; reduceTensorDesc=d, o...)
cudnnReduceTensor!(y, x, d::cudnnReduceTensorDescriptor; o...) = cudnnReduceTensorWithDefaults(x; y, reduceTensorDesc=d, o...)

# Private method
function cudnnReduceTensorWithDefaults(
    x;
    op::cudnnReduceTensorOp_t = CUDNN_REDUCE_TENSOR_ADD,
    compType::DataType = (eltype(x) <: Float64 ? Float64 : Float32),
    nanOpt::cudnnNanPropagation_t = CUDNN_NOT_PROPAGATE_NAN,
    indices::Union{Vector{<:Unsigned},Nothing} = nothing,
    reduceTensorDesc::cudnnReduceTensorDescriptor = cudnnReduceTensorDescriptor(op, cudnnDataType(compType), nanOpt, cudnnReduceTensorIndices(op, indices), cudnnIndicesType(indices)),
    dims::Dims = ntuple(i->1,ndims(x)),
    y = similar(x, dims),
    alpha::Real = 1,
    beta::Real = 0,
    xDesc::cudnnTensorDescriptor = cudnnTensorDescriptor(x),
    yDesc::cudnnTensorDescriptor = cudnnTensorDescriptor(y),
)
    T = eltype(x)
    alpha, beta = scalingParameter(T,alpha), scalingParameter(T,beta)
    cudnnReduceTensorAD(x; reduceTensorDesc, alpha, xDesc, beta, yDesc, y, indices)
end

function cudnnReduceTensorIndices(op, indices)
    if indices !== nothing && op in (CUDNN_REDUCE_TENSOR_MIN, CUDNN_REDUCE_TENSOR_MAX)
        CUDNN_REDUCE_TENSOR_FLATTENED_INDICES
    else
        CUDNN_REDUCE_TENSOR_NO_INDICES
    end
end

cudnnIndicesType(::Nothing)=CUDNN_32BIT_INDICES
cudnnIndicesType(::Vector{UInt8})=CUDNN_8BIT_INDICES
cudnnIndicesType(::Vector{UInt16})=CUDNN_16BIT_INDICES
cudnnIndicesType(::Vector{UInt32})=CUDNN_32BIT_INDICES
cudnnIndicesType(::Vector{UInt64})=CUDNN_64BIT_INDICES
cudnnIndicesType(x)=error("Bad type $x for cudnnIndices, use Vector{UInt8, 16, 32 or 64}.")

# AD method
function cudnnReduceTensorAD(x; reduceTensorDesc, alpha, xDesc, beta, yDesc, y, indices)
    @workspace size=@argout(
        cudnnGetReductionWorkspaceSize(handle(), reduceTensorDesc, xDesc, yDesc, out(Ref{Csize_t}()))
    )[] workspace->begin
        cudnnReduceTensor(handle(), reduceTensorDesc, something(indices, C_NULL), sizeof(indices), workspace, sizeof(workspace), alpha, xDesc, x, beta, yDesc, y)
    end
    return y
end

# Deprecated
function cudnnReduceTensor(op::cudnnReduceTensorOp_t,
                           A::DenseCuArray{T,N}, C::DenseCuArray{T,N};
                           alpha=true, beta=false) where {T,N}
    @warn "cudnnReduceTensor(op,A,C) is deprecated, please use one of the methods in `@doc cudnnReduceTensor`." maxlog=1
    cudnnReduceTensor(A; y=C, op, alpha, beta)
end


text:
function get_lane_length(lane)
    ind = curveindex_end(lane.curve)
    s = lane.curve[ind].s
end

function convert_JuliaLaneId2VTDLaneId(laneTag,roadway)
    seg = roadway.segments[laneTag.segment]
    id_Julia = laneTag.lane
    lane_num = length(seg.lanes)
    id_VTD = -(lane_num-id_Julia)-1
    return id_VTD
end

function initialize_XML()
    doc = XMLDocument()
    r=ElementNode("OpenDRIVE")
    setroot!(doc,r)
    header = addelement!(r,"header","")
    header["revMajor"]=1
    header["revMinor"]=4
    header["name"]=""
    header["version"]="1.00"
    header["data"]= "Tue Mar 21 15:00:43 2017"
    header["north"]=0
    header["south"]=0
    header["east"]=0
    header["west"]=0
    return doc,r
end

function add_road_mark!(lane,markType="broken")
    roadMark=addelement!(lane,"roadMark")
    roadMark["sOffset"]=0.0
    roadMark["type"]=markType
    roadMark["weight"]="standard"
    roadMark["color"]="standard"
    roadMark["width"]=0.5
    roadMark["laneChange"]="both"
    return roadMark
end

function convert_lane!(laneTag,right,roadway)
    laneJulia = roadway.segments[laneTag.segment].lanes[laneTag.lane]
    laneWidth = laneJulia.width
    id = convert_JuliaLaneId2VTDLaneId(laneTag,roadway)
    lane=addelement!(right,"lane")
    lane["id"]=id
    lane["type"]="driving"
    lane["level"]="false"
    width=addelement!(lane,"width")
    width["sOffset"]=0
    width["a"]=laneWidth
    width["b"]=0.0
    width["c"]=0.0
    width["d"]=0.0
    roadMark=add_road_mark!(lane)
    link=addelement!(lane,"link","")
    if !isempty(laneJulia.entrances)
        prevId = convert_JuliaLaneId2VTDLaneId(laneJulia.entrances[1].target.tag,roadway)
        predecessor=addelement!(link,"predecessor")
        predecessor["id"]=prevId
    end
    if !isempty(laneJulia.exits)
        nextId = convert_JuliaLaneId2VTDLaneId(laneJulia.exits[1].target.tag,roadway)
        successor=addelement!(link,"successor")
        successor["id"]=nextId
    end
end

function handle_road_connection!(seg,road,r,roadway)
    link = addelement!(road,"link","")
    road["junction"] = -1
    lane = seg.lanes[1]
    if !isempty(lane.entrances)
        prevSegId = lane.entrances[1].target.tag.segment
        predecessor = addelement!(link,"predecessor")
        predecessor["elementType"]="road"
        predecessor["elementId"]=prevSegId
        predecessor["contactPoint"]="end"
    end
    if !isempty(lane.exits)
        nextSegId = lane.exits[1].target.tag.segment
        successor = addelement!(link,"successor")
        successor["elementType"]="road"
        successor["elementId"]=nextSegId
        successor["contactPoint"]="start"
    end
end

function convert_seg!(seg,r,roadway)
    road = addelement!(r,"road")
    road["name"]=""
    road["id"]=seg.id
    lane=seg.lanes[1]
    s = get_lane_length(lane)
    road["length"]= s
    thetype = addelement!(road,"type")
    thetype["s"]=0
    thetype["type"]="rural"

    handle_road_connection!(seg,road,r,roadway)

    #deal with planView
    planView=addelement!(road,"planView")
    width = lane.width
    geometry = addelement!(planView,"geometry")
    geometry["s"]=0.0
    geometry["x"]=lane.curve[1].pos.x
    geometry["y"]=lane.curve[1].pos.y
    geometry["hdg"]=mod2pi(lane.curve[1].pos.θ)
    geometry["length"]=s
    if lane.curve[1].k == 0 || isnan(lane.curve[1].k)
        line=addelement!(geometry,"line")
    else
        arc=addelement!(geometry,"arc")
        arc["curvature"]=lane.curve[1].k
    end

    #deal with lanes
    lanes=addelement!(road,"lanes")
    laneOffset=addelement!(lanes,"laneOffset")
    a= 0.0
    for i = 1:length(seg.lanes)
        a += seg.lanes[i].width
    end
    a -= seg.lanes[1].width/2
    laneOffset["s"]=0.0
    laneOffset["a"]=a
    laneOffset["b"]=0
    laneOffset["c"]=0
    laneOffset["d"]=0
    laneSection=addelement!(lanes,"laneSection")
    laneSection["s"]=0
    left=addelement!(laneSection,"left")
    center=addelement!(laneSection,"center")
    right=addelement!(laneSection,"right")

    lane=addelement!(center,"lane","")
    lane["id"]=0
    lane["type"]="driving"
    lane["level"]="false"
    add_road_mark!(lane,"solid")

    for lane in seg.lanes
        convert_lane!(lane.tag,right,roadway)
    end
end

function handle_junction!(r,id,junction,roadway)
    junctionId = id
    junctionVTD = addelement!(r,"junction")
    junctionVTD["name"] = ""
    junctionVTD["id"] = junctionId
    connectionID = 0
    for connection in junction.connections
        # Get an attribute value by name.
        connectionVTD = addelement!(junctionVTD,"connection")
        connectionVTD["id"] = connectionID
        connectionID += 1
        connectionVTD["incomingRoad"] = connection.source
        connectionVTD["connectingRoad"] = connection.path
        connectionVTD["contactPoint"] = "start"
        for i = 1:length(connection.laneConnections)
            laneLink = addelement!(connectionVTD,"laneLink")
            laneLink["from"] = convert_JuliaLaneId2VTDLaneId(LaneTag(connection.source,connection.laneConnections[i][1]),roadway)
            laneLink["to"] = convert_JuliaLaneId2VTDLaneId(LaneTag(connection.path,i),roadway)
        end
        for road in findall("road", r)
            if parse(Int,road["id"]) == connection.source
                for successor in findall("link/successor", road)
                    successor["elementType"] = "junction"
                    successor["elementId"] = junctionId
                    delete!(successor, "contactPoint")
                end
                for successor in findall("lanes/laneSection/right/lane/link/successor", road)
                    unlink!(successor)
                end
            end
            if parse(Int,road["id"]) == connection.path
                road["junction"] = junctionId
                for roadi in findall("road", r)
                    if parse(Int,roadi["id"]) == roadway.segments[connection.path].lanes[1].exits[1].target.tag.segment
                        for predecessor in findall("link/predecessor", roadi)
                            predecessor["elementType"] = "junction"
                            predecessor["elementId"] = junctionId
                            delete!(predecessor, "contactPoint")
                        end
                        for predecessor in findall("lanes/laneSection/right/lane/link/predecessor", roadi)
                            unlink!(predecessor)
                        end
                    end
                end
            end
        end
    end
end

function handle_junctions(r,junctions,roadway)
    for i = 1:length(junctions)
        handle_junction!(r,i,junctions[i],roadway)
    end
end

function convert_roadway!(r,roadway)
    for seg in roadway.segments
        convert_seg!(seg,r,roadway)
    end
    #resolveRoadJunction!(r,roadway)
end

max_stars_repo_name: sisl/UrbanDrivingSimulation

max_issues_repo_name: sisl/UrbanDrivingSimulation

max_forks_repo_name: sisl/UrbanDrivingSimulation


text:
function make_test_worker_utility()
    # `true` is debugging switch
    uFct = UtilityFunctionsLH.CRRA(2.0, true);
    return WorkerUtility(uFct, 0.97, 0.0)
end

function validate(u :: WorkerUtility{F1}) where F1
    return  (u.discFactor > 0.5)  &&  !isinf(u.utilFixed)
end

function Base.show(io :: IO, u :: WorkerUtility)
    betaStr = round(u.discFactor, digits = 3);
    uStr = round(u.utilFixed, digits = 2);
    print(io,  "WorkerUtility:  beta: $betaStr,  uFixed: uStr")
end

discount_factor(u :: WorkerUtility{F1}) where F1 = u.discFactor;
util_fixed(u :: WorkerUtility{F1}) where F1 = u.utilFixed;

## ---------------  Methods

"""
	$(SIGNATURES)

Lifetime utility for given lifetime income `ltIncome`, life-span `T`.
Includes fixed utility `utilFixed`.
"""
function lifetime_utility(u :: WorkerUtility{F1}, T :: Integer,  R :: F1, ltIncome) where F1
    pBeta = discount_factor(u);

    utilFixed = UtilityFunctionsLH.lifetime_utility(util_fixed(u), pBeta, T);
    utilFlow  = UtilityFunctionsLH.lifetime_utility(u.util,
        pBeta, R, T, max.(0.0001, ltIncome));
    return utilFixed + utilFlow
end

"""
	$(SIGNATURES)

Marginal utility of wealth.
"""
function mu_wealth(u :: WorkerUtility{F1}, T :: Integer, R :: F1, ltIncome) where F1
    return UtilityFunctionsLH.mu_wealth(u.util, u.discFactor, R, T, ltIncome)
end

"""
	$(SIGNATURES)

Consumption at work start.
"""
function cons_age1(u :: WorkerUtility{F1}, T :: Integer, R :: F1, ltIncome) where F1
    return UtilityFunctionsLH.cons_age1(u.util, u.discFactor, R, T, ltIncome)
end

function cons_path(u :: WorkerUtility{F1}, T :: Integer, R :: F1, ltIncome) where F1
    return UtilityFunctionsLH.cons_path(u.util, u.discFactor, R, T, ltIncome)
end

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


text:
# GraffSDK integration file

using GraffSDK
using DocStringExtensions # temporary while $(SIGNATURES) is in use in this file

const Graff = GraffSDK

"""
    GraffSLAM

An object definition containing the require variables to leverage the server side SLAM solution per user, robot, and session.
"""
mutable struct GraffSLAM
  robotId::AbstractString
  sessionId::AbstractString
  graffconf::Union{Nothing, GraffSDK.GraffConfig}
  robot
  session

  # Constructor
  GraffSLAM(
    robotId::AbstractString,
    sessionId::AbstractString,
    graffconf::Union{Nothing, GraffSDK.GraffConfig};
    robot=nothing,
    session=nothing
   ) = new(
        robotId,
        sessionId,
        graffconf,
        robot,
        session
      )
end

"""
$(SIGNATURES)

Get  Graff configuration, default filepath location assumed as `~/Documents/graffConfig.json`.
"""
function loadGraffConfig(;
            console::Vector{Symbol}=[],
            filepath::AS=joinpath(ENV["HOME"],"Documents","graffConfig.json")
         ) where {AS <: AbstractString}
  #
  println(" - Retrieving Graff Configuration...")
  configFile = open(filepath)
  configData = JSON.parse(readstring(configFile))
  close(configFile)
  cfg = Unmarshal.unmarshal(GraffConfig, configData) # graffConfig =

  :userId in console ? (println("userId: "); cfg.userId = readline(stdin);) : nothing
  :robotId in console ? (println("robotId: "); cfg.robotId = readline(stdin);) : nothing
  :sessionId in console ? (println("sessionId: "); cfg.sessionId = readline(stdin);) : nothing

  return cfg
end

"""
$(SIGNATURES)

Confirm that the robot already exists, create if it doesn't.
"""
function graffRobot(
           graffConfig::GraffConfig,
           robotId::AS
         ) where {AS <: AbstractString}
  #
  println(" - Creating or retrieving robot '$robotId'...")
  robot = nothing
  if(GraffSDK.isRobotExisting(graffConfig, robotId))
      println(" -- Robot '$robotId' already exists, retrieving it...")
      robot = getRobot(graffConfig, robotId)
  else
      # Create a new one
      println(" -- Robot '$robotId' doesn't exist, creating it...")
      newRobot = RobotRequest(robotId, "My New Bot", "Description of my neat robot", "Active")
      robot = addRobot(graffConfig, newRobot)
  end
  println(robot)
  robot
end

"""
$(SIGNATURES)

Get sessions, if it already exists, add to it.
"""
function graffSession(
            graffConfig::GraffConfig,
            robotId::AS,
            sessionId::AS
         ) where {AS<: AbstractString}
  #
  println(" - Creating or retrieving data session '$sessionId' for robot...")
  session = nothing
  if(GraffSDK.isSessionExisting(graffConfig, robotId, sessionId))
      println(" -- Session '$sessionId' already exists for robot '$robotId', retrieving it...")
      session = getSession(graffConfig, robotId, sessionId)
  else
      # Create a new one
      println(" -- Session '$sessionId' doesn't exist for robot '$robotId', creating it...")
      newSessionRequest = SessionDetailsRequest(sessionId, "Submersible vehicle.", "Pose3", false)
      session = addSession(graffConfig, robotId, newSessionRequest)
  end
  println(session)
  session
end

"""
$(SIGNATURES)

Intialize the `sslaml` object using configuration file defined in `graffconfpath`.
"""
function initialize!(sslaml::GraffSLAM;
            graffconfpath::AS=joinpath(ENV["HOME"],"Documents","graffConfig.json")
         ) where {AS <: AbstractString}
  # 1. Get a Graff configuration
  sslaml.graffconf = loadGraffConfig(filepath=graffconfpath)

  # 2. Confirm that the robot already exists, create if it doesn't.
  sslaml.robot = graffRobot(sslaml.graffconf, sslaml.robotId)

  # Make sure that this session is not already populated
  if isSessionExisting(slam_client.graffconf, sslaml.robotId, sslaml.sessionId)
      @warn "There is already a session named '$sessionId' for robot '$robotId'. This example will fail if it tries to add duplicate nodes. We strongly recommend providing a new session name."
      print(" Should we delete it? [Y/N] - ")
      if lowercase(strip(readline(stdin))) == "y"
          deleteSession(sslaml.graffconf, sslaml.robotId, sslaml.sessionId)
          println(" -- Session '$sessionId' deleted!")
      else
          @warn "Okay, but this may break the example! Caveat Userus!"
      end
  end

  # 3. Create or retrieve the session.
  sslaml.session = graffSession(sslaml.graffconf, sslaml.robotId, sslaml.sessionId)

  nothing
end


text:
abstract type VSControl <: InverterComponent end

@def vscontrol_ports begin
    state_input = [:id_o, :iq_o, :id_c, :iq_c, :vd_cap, :vq_cap]
    inner_input = [Vdo_var, Vqo_var, v_control_var, ω_control_var]
end

"""
Parameters of an inner loop controller using virtual impedance, voltage controller and current controller.
# Conmutable structor
```julia
CombinedVIwithVZ(kpv, kiv, kffv, rv, lv, kpc, kic, kffi, ωad, kad)
```

# Arguments
*  `kpv`::Float64 : voltage controller proportional gain
*  `kiv`::Float64 : voltage controller integral gain
*  `kffv`::Float64 : Binary variable to enable feed-forward gain of voltage.
*  `rv`::Float64 : virtual resistance
*  `lv`::Float64 : virtual inductance
*  `kpc`::Float64 : current controller proportional gain
*  `kic`::Float64 : current controller integral gain
*  `kffi`::Float64 : Binary variable to enable feed-forward gain of current
*  `ωad`::Float64 : active damping filter cutoff frequency (rad/sec)
*  `kad`::Float64 : active damping gain
"""
mutable struct CombinedVIwithVZ <: VSControl #(D'Arco EPSR122 Model)
    kpv::Float64
    kiv::Float64
    kffv::Float64
    rv::Float64
    lv::Float64
    kpc::Float64
    kic::Float64
    kffi::Float64
    ωad::Float64
    kad::Float64
    n_states::Int64
    states::Vector{Symbol}
    ports::Ports

        function CombinedVIwithVZ(kpv::Float64,
                                  kiv::Float64,
                                  kffv::Float64,
                                  rv::Float64,
                                  lv::Float64,
                                  kpc::Float64,
                                  kic::Float64,
                                  kffi::Float64,
                                  ωad::Float64,
                                  kad::Float64)

            @vscontrol_ports

            new(kpv,
                kiv,
                kffv,
                rv,
                lv,
                kpc,
                kic,
                kffi,
                ωad,
                kad,
                6,
                [:ξ_d, :ξ_q, :γ_d, :γ_q, :ϕ_d, :ϕ_q],
                Ports(state_input, inner_input))

        end

end

"""
*  `kpc`::Float64 : current controller proportional gain
*  `kic`::Float64 : current controller integral gain
*  `kpv`::Float64 : voltage controller proportional gain
*  `kiv`::Float64 : voltage controller integral gain
*  `ω_ad`::Float64 : active damping filter cutoff frequency (rad/sec)
*  `kad`::Float64 : active damping gain
*  `lv`::Float64 : virtual inductance
*  `rv`::Float64 : virtual resistance
"""
mutable struct DirectControl <: VSControl
    val::Float64
    n_states::Int64
    states::Vector{Symbol}
    ports::Ports

        function DirectControl(val::Float64)

            @vscontrol_ports

            new(0,
                Vector{Symbol}(),
                Ports(state_input, inner_input))

        end

end


text:
"""
    AbstractPreconditioner

Preconditioners for the iterative solvers mostly focused on GPUs

"""
abstract type AbstractPreconditioner end

"""
    BlockJacobiPreconditioner

Creates an object for the block-Jacobi preconditioner

* `nblocks::Int64`: Number of partitions or blocks.
* `blocksize::Int64`: Size of each block.
* `nJs::Int64`: Size of the blocks. For the GPUs these all have to be of equal size.
* `partitions::Vector{Vector{Int64}}``: `npart` partitions stored as lists
* `cupartitions`: `partitions` transfered to the GPU
* `lpartitions::Vector{Int64}``: Length of each partitions.
* `culpartitions::Vector{Int64}``: Length of each partitions, on the GPU.
* `blocks`: Dense blocks of the block-Jacobi
* `cublocks`: `Js` transfered to the GPU
* `map`: The partitions as a mapping to construct views
* `cumap`: `cumap` transferred to the GPU`
* `part`: Partitioning as output by Metis
* `cupart`: `part` transferred to the GPU
* `P`: The sparse precondition matrix whose values are updated at each iteration
"""
struct BlockJacobiPreconditioner{AT,GAT,VI,GVI,MT,GMT,MI,GMI,SMT} <: AbstractPreconditioner
    nblocks::Int64
    blocksize::Int64
    partitions::MI
    cupartitions::Union{GMI,Nothing}
    lpartitions::VI
    culpartitions::Union{GVI,Nothing}
    blocks::AT
    cublocks::Union{GAT,Nothing}
    map::VI
    cumap::Union{GVI,Nothing}
    part::VI
    cupart::Union{GVI,Nothing}
    P::SMT
    id::Union{GMT,MT}
    function BlockJacobiPreconditioner(J, npart, device=CPU()) where {}
        if device == CPU()
            AT  = Array{Float64,3}
            GAT = Nothing
            VI  = Vector{Int64}
            GVI = Nothing
            MT = Matrix{Float64}
            GMT = Nothing
            MI  = Matrix{Int64}
            GMI  = Nothing
            SMT = SparseMatrixCSC{Float64,Int64}
        elseif device == CUDADevice()
            AT  = Array{Float64,3}
            GAT = CuArray{Float64,3}
            VI  = Vector{Int64}
            GVI = CuVector{Int64}
            MT = Matrix{Float64}
            GMT = CuMatrix{Float64}
            MI  = Matrix{Int64}
            GMI = CuMatrix{Int64}
            SMT = CUDA.CUSPARSE.CuSparseMatrixCSR{Float64}
            J = SparseMatrixCSC(J)
        else
            error("Unknown device type")
        end
        m, n = size(J)
        if npart < 2
            error("Number of partitions `npart` should be at" *
                  "least 2 for partitioning in Metis")
        end
        adj = build_adjmatrix(J)
        g = LightGraphs.Graph(adj)
        part = Metis.partition(g, npart)
        partitions = Vector{Vector{Int64}}()
        for i in 1:npart
            push!(partitions, [])
        end
        for (i,v) in enumerate(part)
            push!(partitions[v], i)
        end
        lpartitions = VI(undef, npart)
        lpartitions = length.(partitions)
        blocksize = maximum(length.(partitions))
        blocks = AT(undef, blocksize, blocksize, npart)
        # Get partitions into bit typed structure
        bpartitions = MI(undef, blocksize, npart)
        bpartitions .= 0.0
        for i in 1:npart
            bpartitions[1:length(partitions[i]),i] .= VI(partitions[i])
        end
        id = MT(I, blocksize, blocksize)
        for i in 1:npart
            blocks[:,:,i] .= id
        end
        nmap = 0
        for b in partitions
            nmap += length(b)
        end
        map = VI(undef, nmap)
        part = VI(undef, nmap)
        for b in 1:npart
            for (i,el) in enumerate(partitions[b])
                map[el] = i
                part[el] = b
            end
        end
        row = Vector{Float64}()
        col = Vector{Float64}()
        nzval = Vector{Float64}()

        for b in 1:npart
            for x in partitions[b]
                for y in partitions[b]
                    push!(row, x)
                    push!(col, y)
                    push!(nzval, 1.0)
                end
            end
        end
        P = sparse(row, col, nzval)
        if isa(device, CUDADevice)
            id = GMT(I, blocksize, blocksize)
            cubpartitions = GMI(bpartitions)
            culpartitions = GVI(lpartitions)
            cublocks = GAT(blocks)
            cumap = CUDA.cu(map)
            cupart = CUDA.cu(part)
            P = SMT(P)
        else
            cublocks = nothing
            cubpartitions = nothing
            cumap = nothing
            cupart = nothing
            id = MT(I, blocksize, blocksize)
            culpartitions = nothing
        end
        return new{AT,GAT,VI,GVI,MT,GMT,MI,GMI,SMT}(npart, blocksize, bpartitions, cubpartitions, lpartitions, culpartitions, blocks, cublocks, map, cumap, part, cupart, P, id)
    end
end

"""
    build_adjmatrix

Build the adjacency matrix of a matrix A corresponding to the undirected graph

"""
function build_adjmatrix(A)
    rows = Int64[]
    cols = Int64[]
    vals = Float64[]
    rowsA = rowvals(A)
    m, n = size(A)
    for i = 1:n
        for j in nzrange(A, i)
            push!(rows, rowsA[j])
            push!(cols, i)
            push!(vals, 1.0)

            push!(rows, i)
            push!(cols, rowsA[j])
            push!(vals, 1.0)
        end
    end
    return sparse(rows,cols,vals,size(A,1),size(A,2))
end

"""
    fillblock_gpu

Fill the dense blocks of the preconditioner from the sparse CSR matrix arrays

"""
@kernel function fillblock_gpu!(blocks, blocksize, partition, map, rowPtr, colVal, nzVal, part, lpartitions, id)
    b = @index(Global, Linear)
    for i in 1:blocksize
        for j in 1:blocksize
            blocks[i,j,b] = id[i,j]
        end
    end

    @inbounds for i in 1:lpartitions[b]
        @inbounds for j in rowPtr[partition[i,b]]:rowPtr[partition[i,b]+1]-1
            if b == part[colVal[j]]
                @inbounds blocks[map[partition[i,b]], map[colVal[j]], b] = nzVal[j]
            end
        end
    end
end

"""
    fillP_gpu

Update the values of the preconditioner matrix from the dense Jacobi blocks

"""
@kernel function fillP_gpu!(blocks, partition, map, rowPtr, colVal, nzVal, part, lpartitions)
    b = @index(Global, Linear)
    @inbounds for i in 1:lpartitions[b]
        @inbounds for j in rowPtr[partition[i,b]]:rowPtr[partition[i,b]+1]-1
            if b == part[colVal[j]]
                @inbounds nzVal[j] += blocks[map[partition[i,b]], map[colVal[j]],b]
            end
        end
    end
end

function _update_gpu(j_rowptr, j_colval, j_nzval, p)
    nblocks = p.nblocks
    fillblock_gpu_kernel! = fillblock_gpu!(CUDADevice())
    fillP_gpu_kernel! = fillP_gpu!(CUDADevice())
    # Fill Block Jacobi" begin
    ev = fillblock_gpu_kernel!(p.cublocks, size(p.id,1), p.cupartitions, p.cumap, j_rowptr, j_colval, j_nzval, p.cupart, p.culpartitions, p.id, ndrange=nblocks, dependencies=Event(CUDADevice()))
    wait(ev)
    # Invert blocks begin
    blocklist = Array{CuArray{Float64,2}}(undef, nblocks)
    for b in 1:nblocks
        blocklist[b] = p.cublocks[:,:,b]
    end
    CUDA.@sync pivot, info = CUDA.CUBLAS.getrf_batched!(blocklist, true)
    CUDA.@sync pivot, info, blocklist = CUDA.CUBLAS.getri_batched(blocklist, pivot)
    for b in 1:nblocks
        p.cublocks[:,:,b] .= blocklist[b]
    end
    p.P.nzVal .= 0.0
    # Move blocks to P" begin
    ev = fillP_gpu_kernel!(p.cublocks, p.cupartitions, p.cumap, p.P.rowPtr, p.P.colVal, p.P.nzVal, p.cupart, p.culpartitions, ndrange=nblocks, dependencies=Event(CUDADevice()))
    wait(ev)
    return p.P
end

"""
    function update(J::CuSparseMatrixCSR, p)

Update the preconditioner `p` from the sparse Jacobian `J` in CSR format for the GPU

1) The dense blocks `cuJs` are filled from the sparse Jacobian `J`
2) To a batch inversion of the dense blocks using CUBLAS
3) Extract the preconditioner matrix `p.P` from the dense blocks `cuJs`

"""
function update(J::CUSPARSE.CuSparseMatrixCSR, p)
    _update_gpu(J.rowPtr, J.colVal, J.nzVal, p)
end
function update(J::Transpose{T, CUSPARSE.CuSparseMatrixCSR{T}}, p) where T
    Jt = CUSPARSE.CuSparseMatrixCSC(J.parent)
    _update_gpu(Jt.colPtr, Jt.rowVal, Jt.nzVal, p)
end

@kernel function update_cpu_kernel!(colptr, rowval, nzval, p, lpartitions)
    nblocks = length(p.partitions)
    # Fill Block Jacobi
    b = @index(Global, Linear)
    blocksize = size(p.id,1)
    for i in 1:blocksize
        for j in 1:blocksize
            p.blocks[i,j,b] = p.id[i,j]
        end
    end
    for k in 1:lpartitions[b]
        i = p.partitions[k,b]
        for j in colptr[i]:colptr[i+1]-1
            if b == p.part[rowval[j]]
                @inbounds p.blocks[p.map[rowval[j]], p.map[i], b] = nzval[j]
            end
        end
    end
    # Invert blocks
    p.blocks[:,:,b] .= inv(p.blocks[:,:,b])
    # Move blocks to P
    for k in 1:lpartitions[b]
        i = p.partitions[k,b]
        for j in p.P.colptr[i]:p.P.colptr[i+1]-1
            if b == p.part[p.P.rowval[j]]
                @inbounds p.P.nzval[j] += p.blocks[p.map[p.P.rowval[j]], p.map[i], b]
            end
        end
    end
end

"""
    function update(J::SparseMatrixCSC, p)

Update the preconditioner `p` from the sparse Jacobian `J` in CSC format for the CPU

Note that this implements the same algorithm as for the GPU and becomes very slow on CPU with growing number of blocks.

"""
function update(J::SparseMatrixCSC, p)
    kernel! = update_cpu_kernel!(CPU())
    p.P.nzval .= 0.0
    ev = kernel!(J.colptr, J.rowval, J.nzval, p, p.lpartitions, ndrange=p.nblocks, dependencies=Event(CPU()))
    wait(ev)
    return p.P
end
function update(J::Transpose{T, SparseMatrixCSC{T, I}}, p) where {T, I}
    ix, jx, zx = findnz(J.parent)
    update(sparse(jx, ix, zx), p)
end

is_valid(precond::BlockJacobiPreconditioner) = _check_nan(precond.P)
_check_nan(P::SparseMatrixCSC) = !any(isnan.(P.nzval))
_check_nan(P::CUSPARSE.CuSparseMatrixCSR) = !any(isnan.(P.nzVal))

function Base.show(precond::BlockJacobiPreconditioner)
    npartitions = precond.npart
    nblock = div(size(precond.P, 1), npartitions)
    println("#partitions: $npartitions, Blocksize: n = ", nblock,
            " Mbytes = ", (nblock*nblock*npartitions*8.0)/(1024.0*1024.0))
    println("Block Jacobi block size: $(precond.nJs)")
end


text:
function generate_sor(s::Int, a::Int, rng::AbstractRNG)
    sp = s+a
    o = sp+randn(rng)
    r = sp^2
    return sp, o, r
end

function iter(s, a)
    sp, o, r = generate_sor(s, a, Base.GLOBAL_RNG)
    return sp+o+r
end

N = 1_000_000

s = 1
a = 1
r_sum = 0.0
o_sum = 0.0
@code_warntype iter(s,a)
for i in 1:N
    sp, o, r = generate_sor(s, a, Base.GLOBAL_RNG)
    r_sum += r
    o_sum += o
    s = sp
end

rng = MersenneTwister(1)
s = 1
a = 1
r_sum = 0.0
@time for i in 1:N
    sp, o, r = generate_sor(s, a, rng)
    r_sum += r
    o_sum += o
    s = sp
end
@show r_sum

rng = MersenneTwister(1)
s = 1
a = 1
r_sum = 0.0
@time for i in 1:N
    sp = s+a
    o = sp+randn(rng)
    r = sp^2
    r_sum += r
    o_sum += o
    s = sp
end
@show r_sum

rng = MersenneTwister(1)
s = 1
a = 1
r_sum = 0.0
@time for i in 1:N
    sp = s+a
    o = sp+randn(rng)
    r = sp^2
    r_sum += r
    o_sum += o
    s = sp
end
@show r_sum


text:
export MNIST
module MNIST
    using ImageCore
    using ColorTypes

    export

        traintensor,
        testtensor,

        trainlabels,
        testlabels,

        traindata,
        testdata,

        convert2image,
        convert2features,

        download_helper

    const DEFAULT_DIR = abspath(joinpath(dirname(@__FILE__), "..", "..", "datasets", "mnist"))

    include("Reader/Reader.jl")
    import .Reader.download_helper
    include("interface.jl")
    include("utils.jl")

    Reader.download_helper(; nargs...) = Reader.download_helper(DEFAULT_DIR; nargs...)
end


text:
module TMLEEpistasis

using Base.Sys

if occursin("Intel", cpu_info()[1].model)
    using MKL
end
using DataFrames
using MLJBase
using CSV
using TMLE
using TOML
using BGEN
using HighlyAdaptiveLasso
using EvoTrees
using MLJModels
using MLJLinearModels
using Serialization
using JLD2
using SnpArrays
using Mmap
using HypothesisTests

import MLJBase: fit, transform, target_scitype, input_scitype

###############################################################################
# INCLUDES

include("ukbb.jl")
include("models.jl")
include("estimators.jl")
include("utils.jl")
include("genotypes_and_queries.jl")
include("confounders.jl")
include("phenotypes.jl")
include("grm.jl")
include("sieve_plateau.jl")
include("summary.jl")

###############################################################################
# EXPORTS

export UKBBVariantRun
export build_genotypes_and_queries
export filter_chromosome, merge_beds, adapt_flashpca
export prepare_phenotypes, tmle_phenotypes_batches
export sieve_variance_plateau
export build_summary

end


text: isvalid(::Type{ASCII}, c::Chars) = c < 0x80


text:
function load_data(case::Int64)::Tuple{ms_type,mp_type,ps_type,pp_type,u_type}

    df_sets = CSV.read("$(pwd())/data/df_sets.csv" ;delim=",",use_mmap=false)
    P  = 1:df_sets[1,:P]                                             # set of technologies 𝒫
    G  = 1:df_sets[1,:G]                                             # set of conventional technologies 𝒢
    B  = 1:df_sets[1,:B]                                             # set of storage technologies ℬ
    R  = 1:df_sets[1,:R]                                             # set of renewable technologies ℛ
    S  = 1:df_sets[1,:S]                                             # set of seasons 𝒮
    Hs = 1:df_sets[1,:H]                                             # set of hours ℋ per seasons
    G0 = 1:df_sets[1,:G0]                                            # set of conventional technologies without nuclear
    gn = df_sets[1,:gn]                                              # index of nuclear generation
    b0 = df_sets[1,:b0]                                              # starting index of storage technologies
    r0 = df_sets[1,:r0]                                              # starting index of renewable technologies

    df_unc_sets = CSV.read("$(pwd())/data/df_unc_sets_case$(case).csv" ;delim=",",use_mmap=false)
    NI1 = df_unc_sets[1,:NΩ1]                                        # number of stochastic nodes at present time
    NI2 = df_unc_sets[1,:NΩ2]                                        # number of stochastic nodes in 5 years
    NI3 = df_unc_sets[1,:NΩ3]                                        # number of stochastic nodes in 10 years
    NI0 = NI1 + NI2                                                  # number of investiment nodes
    NI  = NI2 + NI3                                                  # number of operational nodes

    df_unc_params = CSV.read("$(pwd())/data/df_unc_params_case$(case).csv" ;delim=",",use_mmap=false)
    H      = zeros(NI,2)                                             # create empty matrix of rhs terms uncertain paramters H
    C      = zeros(NI,2)                                             # create empty matrix of cost coefficient uncertain paramters C
    H[:,1] = df_unc_params[:,:ρco2 ]                                 # set values of CO₂ budget (H)
    H[:,2] = -df_unc_params[:,:νD  ]                                 # set values of demand growth (H)
    C[:,1] = df_unc_params[:,:c_co2]                                 # set values of CO₂ emission cost (C)
    C[:,2] = df_unc_params[:,:c_ur ]                                 # set values of nuclear fuel cost (C)

    df_inv_params_P = CSV.read("$(pwd())/data/df_inv_params_P.csv" ;delim=",",use_mmap=false)
    c_inv_0   = df_inv_params_P[:,:c_inv_0].*exp10(3)                # investment cost of technologies at present time (£/GW)
    c_inv_5   = df_inv_params_P[:,:c_inv_5  ].*exp10(3)              # investment cost of technologies in 5 years (£/GW)
    c_fix     = df_inv_params_P[:,:c_fix    ].*exp10(3)              # fixed O&M cost of technologies (£/GWyr)
    x_hist_5  = df_inv_params_P[:,:x_hist_5 ].*exp10(-3)             # historical values of technologies in 5 years (GW)
    x_hist_10 = df_inv_params_P[:,:x_hist_10].*exp10(-3)             # historical values of technologies in 10 years (GW)
    x_max     = df_inv_params_P[:,:x_max    ].*exp10(-3)             # maximum capacity of technologies (GW)
    c_inv     = hcat(repeat(c_inv_0 ,1,NI1),repeat(c_inv_5  ,1,NI2)) # investment cost of technologies (£/GW)
    x_hist    = hcat(repeat(x_hist_5,1,NI2),repeat(x_hist_10,1,NI3)) # historical values of technologies (GW)
    πi0       = vcat(1/NI1*ones(NI1),1/NI2*ones(NI2))                # probabilities associated to investiment nodes
    πi        = vcat(1/NI2*ones(NI2),1/NI3*ones(NI3))                # probabilities associated to operational nodes
    itoi0     = Vector{Array{Int64,1}}()                             # map that, for given operational node i, it gives its ancestor investiment nodes
    for i in 1:NI2 push!(itoi0,[1]) end
    for i in (NI2+1):(NI2+NI3) push!(itoi0,[1,Int64(((i-NI2-1)-(i-NI2-1)%NI2)/NI2+1)+1]) end

    df_inv_params_O = CSV.read("$(pwd())/data/df_inv_params_O.csv" ;delim=",",use_mmap=false)
    κ = df_inv_params_O[1,:κ]                                        # years of operational problem (yr)

    df_oper_params_G = CSV.read("$(pwd())/data/df_oper_params_G.csv" ;delim=",",use_mmap=false)
    c_OMvarG = df_oper_params_G[:,:c_varOM]                          # set variable OM cost (conventional technologies) (£/MWh)
    c_fuelG  = df_oper_params_G[:,:c_fuel ]                          # set fuel cost (conventional technologies) (£/MWh)
    em_co2G  = df_oper_params_G[:,:em_co2 ]                          # set CO2 emissions (conventional technologies) (tCO₂/MWh)
    ηG       = df_oper_params_G[:,:η      ]                          # set efficiency (conventional technologies)
    rampG    = df_oper_params_G[:,:ramp   ]                          # ramping limitations (conventional technologies) (MW)

    df_oper_params_B = CSV.read("$(pwd())/data/df_oper_params_B.csv" ;delim=",",use_mmap=false)
    ηB = df_oper_params_B[:,:η]                                      # set efficiency (storage technologies)
    PB = df_oper_params_B[:,:P]                                      # set ch/disch power (storage technologies) (MW/MW)

    df_oper_params_R = CSV.read("$(pwd())/data/df_oper_params_R.csv" ;delim=",",use_mmap=false)
    PR = zeros(R[end],S[end],Hs[end])                                # create empty array of renewable energy production (MW/MW)
    PR[1,1,:] = df_oper_params_R[:,:P1_S1]                           # set values of onshore  wind, season 1
    PR[1,2,:] = df_oper_params_R[:,:P1_S2]                           # set values of onshore  wind, season 2
    PR[1,3,:] = df_oper_params_R[:,:P1_S3]                           # set values of onshore  wind, season 3
    PR[1,4,:] = df_oper_params_R[:,:P1_S4]                           # set values of onshore  wind, season 4
    PR[2,1,:] = df_oper_params_R[:,:P2_S1]                           # set values of offshore wind, season 1
    PR[2,2,:] = df_oper_params_R[:,:P2_S2]                           # set values of offshore wind, season 2
    PR[2,3,:] = df_oper_params_R[:,:P2_S3]                           # set values of offshore wind, season 3
    PR[2,4,:] = df_oper_params_R[:,:P2_S4]                           # set values of offshore wind, season 4
    PR[3,1,:] = df_oper_params_R[:,:P3_S1]                           # set values of pv solar     , season 1
    PR[3,2,:] = df_oper_params_R[:,:P3_S2]                           # set values of pv solar     , season 2
    PR[3,3,:] = df_oper_params_R[:,:P3_S3]                           # set values of pv solar     , season 3
    PR[3,4,:] = df_oper_params_R[:,:P3_S4]                           # set values of pv solar     , season 4

    df_oper_params_D = CSV.read("$(pwd())/data/df_oper_params_D.csv" ;delim=",",use_mmap=false)
    PD = zeros(S[end],Hs[end])                                       # create empty array of energy demand (MW)
    PD[1,:] = df_oper_params_D[:,:P_S1]                              # set values of energy demand, season 1
    PD[2,:] = df_oper_params_D[:,:P_S2]                              # set values of energy demand, season 2
    PD[3,:] = df_oper_params_D[:,:P_S3]                              # set values of energy demand, season 3
    PD[4,:] = df_oper_params_D[:,:P_S4]                              # set values of energy demand, season 4

    df_oper_params_O = CSV.read("$(pwd())/data/df_oper_params_O.csv" ;delim=",",use_mmap=false)
    νD      = df_oper_params_O[1,:νD     ]                           # scaling demand factor
    c_shed  = df_oper_params_O[1,:c_shed ]                           # shedding demand cost (£/MW)
    α       = df_oper_params_O[1,:α      ]                           # seasonal weight
    c_co2   = df_oper_params_O[1,:c_co2  ]                           # CO₂ cost (£/tCO₂)
    co2_lim = df_oper_params_O[1,:co2_lim]                           # CO₂ yearly budget (tCO₂)
    ρco2    = df_oper_params_O[1,:ρco2   ]                           # scaling factor of CO₂ budget

    ms  = ms_type(P,1:NI0,1:NI)
    mp  = mp_type(κ,x_hist,x_max,c_inv,c_fix,πi0,πi,itoi0)
    ps  = ps_type(P,G,B,R,S,Hs,G0,gn,b0,r0)
    pp  = pp_type(c_OMvarG,c_fuelG,em_co2G,ηG,rampG,ηB,PB,PR,PD,νD,c_shed,c_co2,co2_lim,ρco2,α)
    unc = u_type(NI,P[end]+size(H)[2],P[end],size(H)[2],size(C)[2],H,C)

    return ms,mp,ps,pp,unc
end


text:
# Used "Conjugate Bayesian analysis of the Gaussian distribution" by Murphy as
# a reference.  Note that there were some typos in that document so the code
# here may not correspond exactly.

struct NormalGamma{T<:Real} <: ContinuousUnivariateDistribution
    mu::T
    nu::T     # scales precision of Normal
    shape::T
    rate::T

    function NormalGamma{T}(mu::T, nu::T, sh::T, r::T) where T<:Real
    	nu > zero(nu) && sh > zero(sh) && r > zero(r) || error("Both shape and scale must be positive")
    	new(T(mu), T(nu), T(sh), T(r))
    end
end

function NormalGamma(mu::Real, nu::Real, sh::Real, r::Real)
    T = promote_type(typeof(mu), typeof(nu), typeof(sh), typeof(r))
    return NormalGamma{T}(T(mu),T(nu),T(sh),T(r))
end

mu(d::NormalGamma) = d.mu
nu(d::NormalGamma) = d.nu
shape(d::NormalGamma) = d.shape
scale(d::NormalGamma) = 1. / d.rate
rate(d::NormalGamma) = d.rate

insupport(::Type{NormalGamma}, x::T, tau2::T) where T<:Real =
    isfinite(x) && zero(tau2) <= tau2 < Inf

# Probably should guard agains dividing by and taking the log of 0.
function pdf(d::NormalGamma, x::T, tau2::T) where T<:Real
    Zinv = d.rate.^d.shape / gamma(d.shape) * sqrt(d.nu / (2.0*pi))
    return Zinv * tau2.^(d.shape-0.5) * exp(-0.5*tau2*(d.nu*(x-d.mu).^2 + 2.0*d.rate))
end
function logpdf(d::NormalGamma, x::T, tau2::T) where T<:Real
    lZinv = d.shape*log(d.rate) - lgamma(d.shape) + 0.5*(log(d.nu) - log(2.0*pi))
    return lZinv + (d.shape-0.5)*log(tau2) - 0.5*tau2*(d.nu*(x-d.mu).^2 + 2*d.rate)
end

function rand(d::NormalGamma)
    # Guard against invalid precisions
    tau2 = rand(Gamma(d.shape, scale(d)))
    if tau2 <= zero(Float64)
        tau2 = eps(Float64)
    end
    mu = rand(Normal(d.mu, sqrt(1.0/(tau2*d.nu))))
    return mu, tau2
end


text:
using Weave: parse_options, parse_markdown

@testset "`parse_options`" begin

# general
@test isempty(parse_options(""))
@test (:opt => nothing) in parse_options("opt = nothing")

# Weave style -- semicolon separated
let opts = parse_options("opt1 = 1; opt2 = \"2\"")
    @test (:opt1 => 1) in opts
    @test (:opt2 => "2") in opts
end
# robust parsing
@test let opts = parse_options("invalid; valid = nothing")
    @test (:valid => nothing) in opts
    true
end

# RMarkdown style -- comma separated
let opts = parse_options("opt1 = 1, opt2 = \"2\"")
    @test (:opt1 => 1) in opts
    @test (:opt2 => "2") in opts
end
# robust parsing
@test let opts = parse_options("invalid, valid = nothing")
    @test (:valid => nothing) in opts
    true
end

end

@testset "`parse_markdown` Julia markdown format" begin

get_options(str) =
    (chunk = first(last(parse_markdown(str))); @test hasproperty(chunk, :options); chunk.options)

# Julia markdown
@test get_options("```julia\n```") |> length === 1
@test get_options("```julia \n```") |> length === 1
@test get_options("```{julia}\n```") |> length === 1
@test get_options("```{julia }\n```") |> length === 1

# Weave style -- semicolon separated
@test get_options("```julia;\n```") |> length === 1
@test get_options("```julia ;\n```") |> length === 1
@test get_options("```julia; \n```") |> length === 1
@test (:opt => nothing) in get_options("```julia; opt = nothing\n```")
@test (:opt => nothing) in get_options("```{julia; opt = nothing}\n```")
let opts = get_options("```julia; opt1 = 1; opt2 = \"2\"\n```")
    @test (:opt1 => 1) in opts
    @test (:opt2 => "2") in opts
end
let opts = get_options("```{julia; opt1 = 1; opt2 = \"2\"}\n```")
    @test (:opt1 => 1) in opts
    @test (:opt2 => "2") in opts
end

# RMarkdown style -- comma separated
@test get_options("```julia,\n```") |> length === 1
@test get_options("```julia ,\n```") |> length === 1
@test get_options("```julia, \n```") |> length === 1
@test (:opt => nothing) in get_options("```julia, opt = nothing\n```")
@test (:opt => nothing) in get_options("```{julia, opt = nothing}\n```")
let opts = get_options("```{julia, opt1 = 1, opt2 = \"2\"}\n```")
    @test (:opt1 => 1) in opts
    @test (:opt2 => "2") in opts
end
let opts = get_options("```julia, opt1 = 1, opt2 = \"2\"\n```")
    @test (:opt1 => 1) in opts
    @test (:opt2 => "2") in opts
end

end

@testset "`parse_markdown` pandoc format" begin

get_options(str) =
    (chunk = first(last(parse_markdown(str; is_pandoc = true))); @test hasproperty(chunk, :options); chunk.options)

@test get_options("<<>>=\n@") |> length === 1
@test (:opt => nothing) in get_options("<<opt = nothing>>=\n@")
let opts = get_options("<<opt1 = 1; opt2 = \"2\">>=\n@")
    @test (:opt1 => 1) in opts
    @test (:opt2 => "2") in opts
end
let opts = get_options("<<opt1 = 1, opt2 = \"2\">>=\n@")
    @test (:opt1 => 1) in opts
    @test (:opt2 => "2") in opts
end

end

# TODO: tests for `"script"` format ?


text:
# Base
import Base.nextpow
function nextpow2(bins::T) where T<:AbstractRange
    eT = eltype(bins)
    newend = eT(bins.step)*(Base.nextpow(2,length(bins)))

    return bins[1]:eT(bins.step):newend
end
function nextpow(power::Int, bins::T) where T<:AbstractRange
    eT = eltype(bins)
    newend = eT(bins.step)*(Base.nextpow(power,length(bins)))
    return bins[1]:eT(bins.step):newend
end
@deprecate nextpow2 nextpow

## Base - Arrays & Iterators
function slicedimview(A::AbstractArray, d::Integer, i)
    d >= 1 || throw(ArgumentError("dimension must be ≥ 1"))
    nd = ndims(A)
    d > nd && (i == 1 || throw_boundserror(A, (ntuple(k->Colon(),nd)..., ntuple(k->1,d-1-nd)..., i)))
    view(A,Base.setindex(axes(A), i, d)...)
end

unzip(a) = map(x->getfield.(a, x), fieldnames(eltype(a)))
function flatten(X::Array{A,1}) where A <: Array{Z} where Z
    V = Array{Z}(undef, sum(length.(X)))
    st = 1
    for x in X
        lx = length(x)
        V[st:st+lx-1] .= x
        st += lx
    end
    return V
end

@inline function allequal(x::AbstractArray)
    isempty(x) && error("Equality Undefined for empty arrays")
	length(x)  < 2 && return true
    e1 = x[1]
    i = 2
    @inbounds for i=2:length(x)
        x[i] == e1 || return false
    end
    return true
end

using LinearAlgebra
"""
	offdiagind(M::AbstractArray, k::Int=0)
Get indicies not on k-th diagonal of an array.
"""
offdiagind(M::AbstractArray, k::Int=0) = setdiff(LinearIndices(M),diagind(M,k))

"""
	offdiag(x::AbstractArray)
Get all values in array not on the main diagonal.
"""
offdiag(x::AbstractArray) = Base.unsafe_getindex(x,utils.offdiagind(x))

"""
	function rotatexy(x::AbstractVector,y::AbstractVector,θ::Real)
Rotate a set of x and y co-ordinationes by θ degrees.
"""
function rotatexy(x::AbstractVector,y::AbstractVector,θ::Real)
    θ = deg2rad(θ)
    rot = [cos(θ) -sin(θ);
           sin(θ) cos(θ)]
    x2,y2 = (rot*[x y]') |> x->(x[1,:],x[2,:])
    return x2,y2
end

function rotatexy(xy::AbstractArray,θ::Real)
    θ = deg2rad(θ)
    rot = [cos(θ) -sin(θ);
           sin(θ) cos(θ)]
    return rot*xy
end

### find
function _check_length(n::Int, lX::Int)
	if n > lX
		error("n > length(X)")
	elseif n == lX
		return true
	end
	return false
end

"""
	findnmax(X::AbstractArray, n::Int)
Finds `n`-largest elements in array `X`.
"""
function findnmax(X::AbstractArray, n::Int)
    _check_length(n, length(X)) && return collect(LinearIndices(X))
	x = copy(X)
    inds = zeros(Int64,n)
    min = typemin(eltype(X))
    for i in 1:n
        _,inds[i] = findmax(x)
        x[inds[i]] = min
    end
    return inds
end

function findnmin(X::AbstractArray, n::Int)
    _check_length(n, length(X)) && return collect(LinearIndices(X))
	x = copy(X)
    inds = zeros(Int64,n)
    min = typemax(eltype(X))
    for i in 1:n
        _,inds[i] = findmin(x)
        x[inds[i]] = min
    end
    return inds
end


text:
abstract type AbstractPDE{N} end

Geometry.ambient_dimension(pde::AbstractPDE{N}) where {N} = N

"""
    abstract type AbstractKernel{T}

A kernel functions `K` with the signature `K(target,source)::T`.

# See also: [`GenericKernel`](@ref), [`SingleLayerKernel`](@ref), [`DoubleLayerKernel`](@ref), [`AdjointDoubleLayerKernel`](@ref), [`HyperSingularKernel`](@ref)
"""
abstract type AbstractKernel{T} end

return_type(::AbstractKernel{T}) where {T} = T

ambient_dimension(K::AbstractKernel) = ambient_dimension(pde(K))

"""
    pde(K::AbstractKernel)

Return the underlying `AbstractPDE` when `K` correspond to the kernel of an
integral operator derived from a partial differential equation.
"""
pde(k::AbstractKernel) = k.pde

"""
    struct GenericKernel{T,F} <: AbstractKernel{T}

An [`AbstractKernel`](@ref) with `kernel` of type `F`.
"""
struct GenericKernel{T,F} <: AbstractKernel{T}
    kernel::F
end

"""
    struct SingleLayerKernel{T,Op} <: AbstractKernel{T}

Given an operator `Op::Abstract`, construct its free-space single-layer kernel (i.e. the
fundamental solution).
"""
struct SingleLayerKernel{T,Op} <: AbstractKernel{T}
    pde::Op
end
SingleLayerKernel{T}(op) where {T} = SingleLayerKernel{T,typeof(op)}(op)
SingleLayerKernel(op)              = SingleLayerKernel{default_kernel_eltype(op)}(op)

"""
    struct DoubleLayerKernel{T,Op} <: AbstractKernel{T}

Given an operator `Op`, construct its free-space double-layer kernel. This
corresponds to the `γ₁` trace of the [`SingleLayerKernel`](@ref). For operators
such as [`Laplace`](@ref) or [`Helmholtz`](@ref), this is simply the normal
derivative of the fundamental solution respect to the source variable.
"""
struct DoubleLayerKernel{T,Op} <: AbstractKernel{T}
    pde::Op
end
DoubleLayerKernel{T}(op) where {T}     = DoubleLayerKernel{T,typeof(op)}(op)
DoubleLayerKernel(op)                  = DoubleLayerKernel{default_kernel_eltype(op)}(op)

"""
    struct AdjointDoubleLayerKernel{T,Op} <: AbstractKernel{T}

Given an operator `Op`, construct its free-space adjoint double-layer kernel. This
corresponds to the `transpose(γ₁,ₓ[G])`, where `G` is the
[`SingleLayerKernel`](@ref). For operators such as [`Laplace`](@ref) or
[`Helmholtz`](@ref), this is simply the normal derivative of the fundamental
solution respect to the target variable.
"""
struct AdjointDoubleLayerKernel{T,Op} <: AbstractKernel{T}
    pde::Op
end
AdjointDoubleLayerKernel{T}(op) where {T} = AdjointDoubleLayerKernel{T,typeof(op)}(op)
AdjointDoubleLayerKernel(op)              = AdjointDoubleLayerKernel{default_kernel_eltype(op)}(op)

"""
    struct HyperSingularKernel{T,Op} <: AbstractKernel{T}

Given an operator `Op`, construct its free-space hypersingular kernel. This
corresponds to the `transpose(γ₁,ₓγ₁[G])`, where `G` is the
[`SingleLayerKernel`](@ref). For operators such as [`Laplace`](@ref) or
[`Helmholtz`](@ref), this is simply the normal derivative of the fundamental
solution respect to the target variable of the `DoubleLayerKernel`.
"""
struct HyperSingularKernel{T,Op} <: AbstractKernel{T}
    pde::Op
end
HyperSingularKernel{T}(op) where {T} = HyperSingularKernel{T,typeof(op)}(op)
HyperSingularKernel(op)              = HyperSingularKernel{default_kernel_eltype(op)}(op)

# a trait for the kernel type
struct SingleLayer end
struct DoubleLayer end
struct AdjointDoubleLayer end
struct HyperSingular end

kernel_type(::SingleLayerKernel)        = SingleLayer()
kernel_type(::DoubleLayerKernel)        = DoubleLayer()
kernel_type(::AdjointDoubleLayerKernel) = AdjointDoubleLayer()
kernel_type(::HyperSingularKernel)      = HyperSingular()

combined_field_coefficients(::SingleLayerKernel)        = (0,-1)
combined_field_coefficients(::DoubleLayerKernel)        = (1,0)
combined_field_coefficients(::AdjointDoubleLayerKernel) = (0,-1)
combined_field_coefficients(::HyperSingularKernel)      = (1,0)

################################################################################
################################# LAPLACE ######################################
################################################################################

"""
    struct Laplace{N}

Laplace equation in `N` dimension: Δu = 0.
"""
struct Laplace{N} <: AbstractPDE{N} end

Laplace(;dim=3) = Laplace{dim}()

function Base.show(io::IO,pde::Laplace)
    print(io,"Δu = 0")
end

default_kernel_eltype(::Laplace)  = Float64
default_density_eltype(::Laplace) = Float64

function (SL::SingleLayerKernel{T,Laplace{N}})(target,source)::T  where {N,T}
    x = coords(target)
    y = coords(source)
    r = x - y
    d = norm(r)
    d == 0 && (return zero(T))
    if N==2
        return -1/(2π)*log(d)
    elseif N==3
        return 1/(4π)/d
    else
        notimplemented()
    end
end

function (DL::DoubleLayerKernel{T,Laplace{N}})(target,source)::T where {N,T}
    x = coords(target)
    y = coords(source)
    ny = normal(source)
    r = x - y
    d = norm(r)
    d == 0 && (return zero(T))
    if N == 2
        return 1/(2π)/(d^2) .* dot(r,ny)
    elseif N==3
        return 1/(4π)/(d^3) .* dot(r,ny)
    else
        notimplemented()
    end
end

function (ADL::AdjointDoubleLayerKernel{T,Laplace{N}})(target,source)::T where {N,T}
    x = coords(target)
    y = coords(source)
    nx = normal(target)
    r = x - y
    d = norm(r)
    d == 0 && (return zero(T))
    if N==2
        return -1/(2π)/(d^2) .* dot(r,nx)
    elseif N==3
        return -1/(4π)/(d^3) .* dot(r,nx)
    end
end

function (HS::HyperSingularKernel{T,Laplace{N}})(target,source)::T where {N,T}
    x = coords(target)
    y = coords(source)
    nx = normal(target)
    ny = normal(source)
    r = x - y
    d = norm(r)
    d == 0 && (return zero(T))
    if N==2
        ID = SMatrix{2,2,Float64,4}(1,0,0,1)
        RRT = r*transpose(r) # r ⊗ rᵗ
        return 1/(2π)/(d^2) * transpose(nx)*(( ID -2*RRT/d^2  )*ny)
    elseif N==3
        ID = SMatrix{3,3,Float64,9}(1,0,0,0,1,0,0,0,1)
        RRT = r*transpose(r) # r ⊗ rᵗ
        return 1/(4π)/(d^3) * transpose(nx)*(( ID -3*RRT/d^2  )*ny)
    end
end

################################################################################
################################# Helmholtz ####################################
################################################################################

"""
    struct Helmholtz{N,T}

Helmholtz equation in `N` dimensions: Δu + k²u = 0.
"""
struct Helmholtz{N,K} <: AbstractPDE{N}
    k::K
end

Helmholtz(;k,dim=3) = Helmholtz{dim,typeof(k)}(k)

function Base.show(io::IO,pde::Helmholtz)
    # k = parameters(pde)
    print(io,"Δu + k u = 0")
end

parameters(pde::Helmholtz) = pde.k

default_kernel_eltype(::Helmholtz)  = ComplexF64
default_density_eltype(::Helmholtz) = ComplexF64

function (SL::SingleLayerKernel{T,S})(target, source)::T where {T,S <: Helmholtz}
    x = coords(target)
    y = coords(source)
    N = ambient_dimension(pde(SL))
    k = parameters(pde(SL))
    r = x - y
    d = norm(r)
    d == 0 && (return zero(T))
    if N == 2
        return im / 4 * hankelh1(0, k * d)
    elseif N == 3
        return 1 / (4π) / d * exp(im * k * d)
    end
end

# Double Layer Kernel
function (DL::DoubleLayerKernel{T,S})(target, source)::T where {T,S <: Helmholtz}
    x, y, ny = coords(target), coords(source), normal(source)
    N = ambient_dimension(pde(DL))
    k = parameters(pde(DL))
    r = x - y
    d = norm(r)
    d == 0 && (return zero(T))
    if N == 2
        return im * k / 4 / d * hankelh1(1, k * d) .* dot(r, ny)
    elseif N == 3
        return 1 / (4π) / d^2 * exp(im * k * d) * ( -im * k + 1 / d ) * dot(r, ny)
    end
end

# Adjoint double Layer Kernel
function (ADL::AdjointDoubleLayerKernel{T,S})(target, source)::T where {T,S <: Helmholtz}
    x, y, nx = coords(target), coords(source), normal(target)
    N = ambient_dimension(pde(ADL))
    k = parameters(pde(ADL))
    r = x - y
    d = norm(r)
    d == 0 && (return zero(T))
    if N == 2
        return -im * k / 4 / d * hankelh1(1, k * d) .* dot(r, nx)
    elseif N == 3
        return -1 / (4π) / d^2 * exp(im * k * d) * ( -im * k + 1 / d ) * dot(r, nx)
    end
end

# Hypersingular kernel
function (HS::HyperSingularKernel{T,S})(target, source)::T where {T,S <: Helmholtz}
    x, y, nx, ny = coords(target), coords(source), normal(target), normal(source)
    N = ambient_dimension(pde(HS))
    k = parameters(pde(HS))
    r = x - y
    d = norm(r)
    d == 0 && (return zero(T))
    if N == 2
        RRT = r * transpose(r) # r ⊗ rᵗ
        # TODO: rewrite the operation below in a more clear/efficient way
        return transpose(nx) * ((-im * k^2 / 4 / d^2 * hankelh1(2, k * d) * RRT + im * k / 4 / d * hankelh1(1, k * d) * I) * ny)
    elseif N == 3
        RRT   = r * transpose(r) # r ⊗ rᵗ
        term1 = 1 / (4π) / d^2 * exp(im * k * d) * ( -im * k + 1 / d ) * I
        term2 = RRT / d * exp(im * k * d) / (4 * π * d^4) * (3 * (d * im * k - 1) + d^2 * k^2)
        return  transpose(nx) * (term1 + term2) * ny
    end
end

################################################################################
################################# Elastostatic #################################
################################################################################

"""
    struct Elastostatic{N,T} <: AbstractPDE{N}

Elastostatic equation in `N` dimensions: μΔu + (μ+λ)∇(∇⋅u) = 0. Note that the
displacement u is a vector of length `N` since this is a vectorial problem.
"""
struct Elastostatic{N,T} <: AbstractPDE{N}
    μ::T
    λ::T
end
Elastostatic(;μ,λ,dim=3)               = Elastostatic{dim}(promote(μ,λ)...)
Elastostatic{N}(μ::T,λ::T) where {N,T} = Elastostatic{N,T}(μ,λ)

function Base.show(io::IO,pde::Elastostatic)
    print(io,"μΔu + (μ+λ)∇(∇⋅u) = 0")
end

parameters(pde::Elastostatic) = pde.μ, pde.λ

default_kernel_eltype(::Elastostatic{N}) where {N}  = SMatrix{N,N,Float64,N*N}
default_density_eltype(::Elastostatic{N}) where {N} = SVector{N,Float64}

function (SL::SingleLayerKernel{T,S})(target,source)::T  where {T,S<:Elastostatic}
    N   = ambient_dimension(pde(SL))
    μ,λ = parameters(pde(SL))
    ν = λ/(2*(μ+λ))
    x = coords(target)
    y = coords(source)
    r = x .- y
    d = norm(r)
    d == 0 && return zero(T)
    RRT = r*transpose(r) # r ⊗ rᵗ
    if N==2
        ID = SMatrix{2,2,Float64,4}(1,0,0,1)
        return 1/(8π*μ*(1-ν))*(-(3-4*ν)*log(d)*ID + RRT/d^2)
        # return (λ + 3μ)/(4*π*(N-1)*μ*(λ+2μ))* (-log(d)*ID + (λ+μ)/(λ+3μ)*RRT/d^2)
    elseif N==3
        ID = SMatrix{3,3,Float64,9}(1,0,0,0,1,0,0,0,1)
        return 1/(16π*μ*(1-ν)*d)*((3-4*ν)*ID + RRT/d^2)
    end
end

function (DL::DoubleLayerKernel{T,S})(target,source)::T where {T,S<:Elastostatic}
    N = ambient_dimension(pde(DL))
    μ,λ = parameters(pde(DL))
    ν = λ/(2*(μ+λ))
    x = coords(target)
    y = coords(source)
    ny = normal(source)
    ν = λ/(2*(μ+λ))
    r = x .- y
    d = norm(r)
    d == 0 && return zero(T)
    RRT = r*transpose(r) # r ⊗ rᵗ
    drdn = -dot(r,ny)/d
    if N==2
        ID = SMatrix{2,2,Float64,4}(1,0,0,1)
        return -1/(4π*(1-ν)*d)*(drdn*((1-2ν)*ID + 2*RRT/d^2) + (1-2ν)/d*(r*transpose(ny) - ny*transpose(r)))
    elseif N==3
        ID = SMatrix{3,3,Float64,9}(1,0,0,0,1,0,0,0,1)
        return -1/(8π*(1-ν)*d^2)*(drdn * ((1-2*ν)*ID + 3*RRT/d^2) + (1-2*ν)/d*(r*transpose(ny) - ny*transpose(r)))
    end
end

function (ADL::AdjointDoubleLayerKernel{T,S})(target,source)::T where {T,S<:Elastostatic}
    # DL = DoubleLayerKernel{T}(pde(ADL))
    # return -DL(x,y,nx) |> transpose
    N   = ambient_dimension(pde(ADL))
    μ,λ = parameters(pde(ADL))
    ν = λ/(2*(μ+λ))
    x = coords(target)
    nx = normal(target)
    y = coords(source)
    ν = λ/(2*(μ+λ))
    r = x .- y
    d = norm(r)
    d == 0 && return zero(T)
    RRT = r*transpose(r) # r ⊗ rᵗ
    drdn = -dot(r,nx)/d
    if N==2
        ID = SMatrix{2,2,Float64,4}(1,0,0,1)
        out = -1/(4π*(1-ν)*d)*(drdn*((1-2ν)*ID + 2*RRT/d^2) + (1-2ν)/d*(r*transpose(nx) - nx*transpose(r)))
        return -transpose(out)
    elseif N==3
        ID = SMatrix{3,3,Float64,9}(1,0,0,0,1,0,0,0,1)
        out =  -1/(8π*(1-ν)*d^2)*(drdn * ((1-2*ν)*ID + 3*RRT/d^2) + (1-2*ν)/d*(r*transpose(nx) - nx*transpose(r)))
        return -transpose(out)
    end
end

function (HS::HyperSingularKernel{T,S})(target,source)::T where {T,S<:Elastostatic}
    N = ambient_dimension(pde(HS))
    μ,λ = parameters(pde(HS))
    ν = λ/(2*(μ+λ))
    x = coords(target)
    nx = normal(target)
    y = coords(source)
    ny = normal(source)
    r = x .- y
    d = norm(r)
    d == 0 && return zero(T)
    RRT   = r*transpose(r) # r ⊗ rᵗ
    drdn  = dot(r,ny)/d
    if N==2
        ID = SMatrix{2,2,Float64,4}(1,0,0,1)
        return μ/(2π*(1-ν)*d^2)* (2*drdn/d*( (1-2ν)*nx*transpose(r) + ν*(dot(r,nx)*ID + r*transpose(nx)) - 4*dot(r,nx)*RRT/d^2 ) +
                                  2*ν/d^2*(dot(r,nx)*ny*transpose(r) + dot(nx,ny)*RRT) +
                                  (1-2*ν)*(2/d^2*dot(r,nx)*r*transpose(ny) + dot(nx,ny)*ID + ny*transpose(nx)) -
                                  (1-4ν)*nx*transpose(ny)
                                  )
    elseif N==3
        ID = SMatrix{3,3,Float64,9}(1,0,0,0,1,0,0,0,1)
        return μ/(4π*(1-ν)*d^3)* (3*drdn/d*( (1-2ν)*nx*transpose(r) + ν*(dot(r,nx)*ID + r*transpose(nx)) - 5*dot(r,nx)*RRT/d^2 ) +
                                  3*ν/d^2*(dot(r,nx)*ny*transpose(r) + dot(nx,ny)*RRT) +
                                  (1-2*ν)*(3/d^2*dot(r,nx)*r*transpose(ny) + dot(nx,ny)*ID + ny*transpose(nx)) -
                                  (1-4ν)*nx*transpose(ny)
                                  )
    end
end

################################################################################
################################# Maxwell ######################################
################################################################################

"""
    Maxwell{T} <: AbstractPDE{3}

Normalized Maxwell's equation ∇ × ∇ × E - k² E = 0, where
k = ω √ϵμ.
"""
struct Maxwell{T} <: AbstractPDE{3}
    k::T
end

Maxwell(;dim=3,k::T) where {T}        = Maxwell{T}(k)

parameters(pde::Maxwell) = pde.k

function Base.show(io::IO,pde::Maxwell)
    # k = parameters(pde)
    print(io,"∇ × ∇ × E - k² E = 0")
end

default_kernel_eltype(::Maxwell)   = SMatrix{3,3,ComplexF64,9}
default_density_eltype(::Maxwell)  = SVector{3,ComplexF64}

# Single Layer kernel for Maxwell is the dyadic Greens function
function (SL::SingleLayerKernel{T,S})(target,source)::T  where {T,S<:Maxwell}
    k  = parameters(pde(SL))
    x  = coords(target)
    y  = coords(source)
    rvec = x - y
    r = norm(rvec)
    r == 0 && return zero(T)
    # helmholtz greens function
    g   = exp(im*k*r)/(4π*r)
    gp  = im*k*g - g/r
    gpp = im*k*gp - gp/r + g/r^2
    RRT = rvec*transpose(rvec) # rvec ⊗ rvecᵗ
    G   = g*I + 1/k^2*(gp/r*I + (gpp/r^2 - gp/r^3)*RRT)
    # TODO: when multiplying by a density, it is faster to exploit the outer
    # product format isntead of actually assemblying the matrices.
    return G
end

# Double Layer Kernel
function (DL::DoubleLayerKernel{T,S})(target,source)::T where {T,S<:Maxwell}
    k  = parameters(pde(DL))
    x  = coords(target)
    y  = coords(source)
    ny = normal(source)
    rvec = x - y
    r      = norm(rvec)
    r == 0 && return zero(T)
    g      = exp(im*k*r)/(4π*r)
    gp     = im*k*g - g/r
    rcross = cross_product_matrix(rvec)
    ncross = cross_product_matrix(ny)
    return ncross * rcross * gp/r |> transpose
end

"""
    _curl_y_green_tensor_maxwell(x, y, k)

Returns `∇ʸ × G(x, y)` where `G` is the Green tensor for Maxwell's equations
with wavenumber `k`.
"""
function _curl_y_green_tensor_maxwell(x, y, k)
    rvec = x - y
    r    = norm(rvec)
    g    = exp(im*k*r)/(4π*r)
    gp   = im*k*g - g/r
    rcross = cross_product_matrix(rvec)
    curl_G = -gp/r*rcross
    return curl_G
end

"""
    _curl_x_green_tensor_maxwell(x, y, k)

Returns `∇ˣ × G(x, y)` where `G` is the Green tensor for Maxwell's equations
with wavenumber `k`.
"""
function _curl_x_green_tensor_maxwell(x, y, k)
    rvec = x - y
    r = norm(rvec)
    g    = exp(im*k*r)/(4π*r)
    gp  = im*k*g - g/r
    rcross = cross_product_matrix(rvec)
    curl_G = gp/r*rcross
    return curl_G
end

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom