text:
# This file includes the simulation of four node network simulation

import Cairo, Fontconfig
using GraphPlot
using Graphs
using SimpleWeightedGraphs
using LinearAlgebra
using JLD2
using DynamicalNetworks
using DifferentialEquations
using PyPlot
pygui(true)
close("all")

# Helper functions
getwave(x, i, d) = getindex.(x, (i - 1) * d + 1)

# Node settings
nodetype = Lorenz           # Node type
P = diagm([1.0, 0.0, 0])    # State coupling matrix
ηt = 10.0                   # Threshold for maximum Lyapuno exponents

# Topology
n = 8
graph = SimpleWeightedGraph(cycle_graph(n))
add_edge!(graph, 2, 7, 1.0)
add_edge!(graph, 3, 6, 1.0)
gplot(graph, nodelabel = 1:nv(graph), layout = circular_layout)

# Coupling strength for the first symmetry
W1 = laplacian_matrix(graph)
Q1 = [
    0 0 0 0 0 0 0 1
    0 0 0 0 0 0 1 0
    0 0 0 0 0 1 0 0
    0 0 0 0 1 0 0 0
    0 0 0 1 0 0 0 0
    0 0 1 0 0 0 0 0
    0 1 0 0 0 0 0 0
    1 0 0 0 0 0 0 0
]
M1 = eigvecs(Q1)
G1 = inv(M1) * (W1 ./ diag(W1)) * M1
B1 = G1[1:4, 1:4]
D1 = G1[5:8, 5:8]
λt1 = eigvals(B1)
λs1 = eigvals(D1)
λt1[1] > λs1[2] || error("Synchronization is not possible")
ϵr1 = ηt ./ [λt1[1], λs1[2]]

# Coupling strength for the second symmetry
graph = SimpleWeightedGraph(cycle_graph(n))
add_edge!(graph, 2, 7, 1.0)
add_edge!(graph, 3, 6, 1.0)
add_edge!(graph, 5, 8, 5.0)
W2 = laplacian_matrix(graph)
Q2 = [
    0 0 0 1 0 0 0 0
    0 0 1 0 0 0 0 0
    0 1 0 0 0 0 0 0
    1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1
    0 0 0 0 0 0 1 0
    0 0 0 0 0 1 0 0
    0 0 0 0 1 0 0 0
]
M2 = eigvecs(Q2)
G2 = inv(M2) * (W2 ./ diag(W2)) * M2
B2 = G2[1:4, 1:4]
D2 = G2[5:8, 5:8]
λt2 = eigvals(B2)
λs2 = eigvals(D2)
λt2[1] > λs2[2] || error("Synchronization is not possible")
ϵr2 = ηt ./ [λt2[1], λs2[2]]

ϵl = maximum([ϵr1[1], ϵr2[1]])
ϵh = minimum([ϵr2[2], ϵr2[2]])
ϵ = ϵl + abs(ϵh - ϵl) * 0.5

# Callback
nbits = 100
tbit = 100.0
bits = rand(Bool, nbits)
ti = 0.0
tf = nbits * tbit
tswitch = ti:tbit:tf-tbit
srcnode1, dstnode1 = 5, 8     # Control link nodes
high = 5.0          # Control coupling strength
low = eps()         # Control coupling strength
condition(u, t, integrator) = integrator.t ∈ tswitch
function action!(integrator)
    net, bits = integrator.sol.prob.p
    if popfirst!(bits)
        add_edge!(net.graph, srcnode1, dstnode1, high)
    else
        add_edge!(net.graph, srcnode1, dstnode1, low)
    end
end
callback = DiscreteCallback(condition, action!)

# Network
σ = 0.0  # Noise strength
H = getnoisemat(graph) * σ
net = SDENetwork(nodetype, ϵ, graph, H, P)
sol = solvenet!(net, (ti, tf), (copy(bits),), solverkwargs = (saveat = 0.01,), probkwargs = (callback = callback, tstops = tswitch))
t = sol.t
x = sol.u

# Save simulation
datapath = joinpath(@__DIR__, "data")
isdir(datapath) || mkpath(datapath)
filename = joinpath(datapath, "data.jld2")
@save filename t x bits tbit


text:
cd(@__DIR__); import Pkg; Pkg.activate(".")
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
Pkg.develop(path="$(homedir())/github/SimpleSDMLayers.jl")
using Revise
using SimpleSDMLayers
using Test

# Import QGIS-cropped layer as example
ref = geotiff(SimpleSDMPredictor, "../../rasters/wc1_cropped.tif")

## Step 4: fix for geotiff subset dimensions bug ####

# Load things
layer = SimpleSDMPredictor(WorldClim, BioClim, 1)
layer_subset = layer[coords]
size(layer_subset) # should be 330, 570
left, right, bottom, top = coords
using Test
using ArchGDAL

# Preliminary tests
test = SimpleSDMPredictor(WorldClim, BioClim, 1; coords...)
ref
layer_subset
@test isequal(ref.grid, layer_subset.grid) # true
@test isequal(ref.grid, test.grid) # false
@test size(ref) == size(test) # false
@test stride(ref) == stride(layer_subset) # true
@test stride(ref) == stride(test)

## Investigate function
# SimpleSDMPredictor(...)
resolution = 10.0
file = SimpleSDMLayers._get_raster(WorldClim, BioClim, 1, resolution)
geotiff(SimpleSDMPredictor, file)
geotiff(SimpleSDMPredictor, file; coords...) # wrong dimensions

# geotiff(SimpleSDMPredictor, file; coords...)
import SimpleSDMLayers: _find_span
LT = SimpleSDMPredictor
bandnumber=1
dataset = ArchGDAL.read(file)
# function geotiff(
#     ::Type{LT},
#     file::AbstractString,
#     bandnumber::Integer=1;
#     left = -180.0,
#     right = 180.0,
#     bottom = -90.0,
#     top = 90.0
# ) where {LT<:SimpleSDMLayer}

    # This next block is reading the geotiff file, but also making sure that we
    # clip the file correctly to avoid reading more than we need.
    # This next block is reading the geotiff file, but also making sure that we
    # clip the file correctly to avoid reading more than we need.
    # layer = ArchGDAL.read(file) do dataset

        transform = ArchGDAL.getgeotransform(dataset)
        wkt = ArchGDAL.getproj(dataset)

        # The data we need is pretty much always going to be stored in the first
        # band, but this is not the case for the future WorldClim data.
        band = ArchGDAL.getband(dataset, bandnumber)
        T = ArchGDAL.pixeltype(band)

        # The nodata is not always correclty identified, so if it is not found, we assumed it is the smallest value in the band
        nodata = isnothing(ArchGDAL.getnodatavalue(band)) ? convert(T, ArchGDAL.minimum(band)) : convert(T, ArchGDAL.getnodatavalue(band))

        # Get the correct latitudes
        minlon = transform[1]
        maxlat = transform[4]
        maxlon = minlon + size(band,1)*transform[2]
        minlat = maxlat - abs(size(band,2)*transform[6])

        left = isnothing(left) ? minlon : max(left, minlon)
        right = isnothing(right) ? maxlon : min(right, maxlon)
        bottom = isnothing(bottom) ? minlat : max(bottom, minlat)
        top = isnothing(top) ? maxlat : min(top, maxlat)

        lon_stride, lat_stride = transform[2], transform[6]

        width = ArchGDAL.width(dataset)
        height = ArchGDAL.height(dataset)

        #global lon_stride, lat_stride
        #global left_pos, right_pos
        #global bottom_pos, top_pos

        lon_stride, left_pos, min_width = _find_span(width, minlon, maxlon, left) # wrong
        _, right_pos, max_width = _find_span(width, minlon, maxlon, right) # ok
        lat_stride, top_pos, max_height = _find_span(height, minlat, maxlat, top) # ok
        _, bottom_pos, min_height = _find_span(height, minlat, maxlat, bottom) # wrong

        max_height, min_height = height .- (min_height, max_height) .+ 1

        # We are now ready to initialize a matrix of the correct type.
        buffer = Matrix{T}(undef, length(min_width:max_width), length(min_height:max_height))
        ArchGDAL.read!(dataset, buffer, bandnumber, min_height:max_height, min_width:max_width)
        buffer = convert(Matrix{Union{Nothing,eltype(buffer)}}, rotl90(buffer))
        replace!(buffer, nodata => nothing)
        LT(buffer, left_pos-0.5lon_stride, right_pos+0.5lon_stride, bottom_pos-0.5lat_stride, top_pos+0.5lat_stride)
    # end

    # return layer

# end

# _find_span
_find_span(width, minlon, maxlon, left) # wrong
_find_span(width, minlon, maxlon, right) # ok
_find_span(height, minlat, maxlat, top) # ok
_find_span(height, minlat, maxlat, bottom) # wrong
n = width
m = minlon
M = maxlon
pos = left
# function _find_span(n, m, M, pos)
    pos > M && return nothing
    pos < m && return nothing
    # stride = (M - m) / n # cannot assign variable to function 🤦
    stride1 = (M - m) / n
    centers = (m + 0.5stride1):stride1:(M-0.5stride1)
    span_pos = last(findmin(abs.(pos .- centers)))
    return (stride, centers[span_pos], span_pos)
# end
# Stupid findmin thing 😡
last(findmin(abs.(pos .- centers)))
pos_diff = abs.(pos .- centers)
last(findmin(pos_diff))
filter(x -> x <= stride1, pos_diff)
findall(x -> x <= stride1, pos_diff)
pos_equal = isapprox.(pos_diff, 0.5stride1)
# Test time
@time _find_span(width, minlon, maxlon, left) # 3 alloc
# With fix
@time _find_span(width, minlon, maxlon, left) # error
@time _find_span(width, minlon, maxlon, left, :none) # error
@time _find_span(width, minlon, maxlon, left, :left) # 10 alloc, ok
@time _find_span(width, minlon, maxlon, right, :right) # ok
@time _find_span(height, minlat, maxlat, top, :top) # ok
@time _find_span(height, minlat, maxlat, bottom, :bottom) # ok
# Test fix
test = SimpleSDMPredictor(WorldClim, BioClim, 1; coords...)
@test size(test) == size(ref)
@test test.grid == ref.grid
@test stride(test) == stride(ref)

ext: jl


text:
module CombinatorialSpaces

using Reexport
using Requires

include("ArrayUtils.jl")
include("CombinatorialMaps.jl")
include("SimplicialSets.jl")
include("DualSimplicialSets.jl")
include("MeshInterop.jl")

function __init__()
  @require AbstractPlotting="537997a7-5e4e-5d89-9595-2241ea00577e" include("MeshGraphics.jl")
end

@reexport using .SimplicialSets
@reexport using .DualSimplicialSets
@reexport using .MeshInterop

end


text:
################################################################################
#
#   (10,3)c
#
################################################################################

# Implementation
# - implementation 1
# - labels <: Any
function getUnitcell_10_3_c(
			unitcell_type  :: Type{U},
			implementation :: Val{1}
		) :: U where {LS,LB,S<:AbstractSite{LS,3},B<:AbstractBond{LB,3},U<:AbstractUnitcell{S,B}}

	# return a new Unitcell
	return newUnitcell(
		# Type of the unitcell
		U,
		# lattice vectors
		Vector{Float64}[
			Float64[1, 0, 0],
			Float64[-0.5, sqrt(3)/2., 0.0],
			Float64[0.0, 0.0, (3*sqrt(3))/2.],
		],
		# sites
		S[
			newSite(S, Float64[0.25, 1/(4*sqrt(3)), 1/(2*sqrt(3))], getDefaultLabelN(LS,1)),
			newSite(S, Float64[0.75, 1/(4*sqrt(3)), 2/sqrt(3)], getDefaultLabelN(LS,2)),
			newSite(S, Float64[0.5, 1/sqrt(3), 7/(2*sqrt(3))], getDefaultLabelN(LS,3)),
			newSite(S, Float64[0.75, 1/(4*sqrt(3)), 1/sqrt(3)], getDefaultLabelN(LS,4)),
			newSite(S, Float64[0.5, 1/sqrt(3), 5/(2*sqrt(3))], getDefaultLabelN(LS,5)),
			newSite(S, Float64[0.25, 1/(4*sqrt(3)), 4/sqrt(3)], getDefaultLabelN(LS,6)),
		],
		# bonds
		B[
			newBond(B, 1, 4, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 4, 2, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 2, 5, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 5, 3, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 3, 6, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 4, 1, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 2, 4, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 5, 2, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 3, 5, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 6, 3, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 4, 1, getDefaultLabel(LB), (1, 0, 0)),
			newBond(B, 1, 4, getDefaultLabel(LB), (-1, 0, 0)),
			newBond(B, 5, 2, getDefaultLabel(LB), (0, 1, 0)),
			newBond(B, 2, 5, getDefaultLabel(LB), (0, -1, 0)),
			newBond(B, 3, 6, getDefaultLabel(LB), (1, 1, 0)),
			newBond(B, 6, 3, getDefaultLabel(LB), (-1, -1, 0)),
			newBond(B, 1, 6, getDefaultLabel(LB), (0, 0, -1)),
			newBond(B, 6, 1, getDefaultLabel(LB), (0, 0, 1)),
		]
	)
end

# Implementation
# - implementation 4 (Kitaev)
# - labels <: Any
function getUnitcell_10_3_c(
			unitcell_type  :: Type{U},
			implementation :: Val{4}
		) :: U where {LS,LB,S<:AbstractSite{LS,3},B<:AbstractBond{LB,3},U<:AbstractUnitcell{S,B}}

	# return a new Unitcell
	return newUnitcell(
		# Type of the unitcell
		U,
		# lattice vectors
		Vector{Float64}[
			Float64[1, 0, 0],
			Float64[-0.5, sqrt(3)/2., 0.0],
			Float64[0.0, 0.0, (3*sqrt(3))/2.],
		],
		# sites
		S[
			newSite(S, Float64[0.25, 1/(4*sqrt(3)), 1/(2*sqrt(3))], getDefaultLabelN(LS,1)),
			newSite(S, Float64[0.75, 1/(4*sqrt(3)), 2/sqrt(3)], getDefaultLabelN(LS,2)),
			newSite(S, Float64[0.5, 1/sqrt(3), 7/(2*sqrt(3))], getDefaultLabelN(LS,3)),
			newSite(S, Float64[0.75, 1/(4*sqrt(3)), 1/sqrt(3)], getDefaultLabelN(LS,4)),
			newSite(S, Float64[0.5, 1/sqrt(3), 5/(2*sqrt(3))], getDefaultLabelN(LS,5)),
			newSite(S, Float64[0.25, 1/(4*sqrt(3)), 4/sqrt(3)], getDefaultLabelN(LS,6)),
		],
		# bonds
		B[
			newBond(B, 1, 4, getDefaultLabelX(LB), (0, 0, 0)),
			newBond(B, 4, 2, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 2, 5, getDefaultLabelX(LB), (0, 0, 0)),
			newBond(B, 5, 3, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 3, 6, getDefaultLabelX(LB), (0, 0, 0)),
			newBond(B, 4, 1, getDefaultLabelX(LB), (0, 0, 0)),
			newBond(B, 2, 4, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 5, 2, getDefaultLabelX(LB), (0, 0, 0)),
			newBond(B, 3, 5, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 6, 3, getDefaultLabelX(LB), (0, 0, 0)),
			newBond(B, 4, 1, getDefaultLabelY(LB), (1, 0, 0)),
			newBond(B, 1, 4, getDefaultLabelY(LB), (-1, 0, 0)),
			newBond(B, 5, 2, getDefaultLabelY(LB), (0, 1, 0)),
			newBond(B, 2, 5, getDefaultLabelY(LB), (0, -1, 0)),
			newBond(B, 3, 6, getDefaultLabelY(LB), (1, 1, 0)),
			newBond(B, 6, 3, getDefaultLabelY(LB), (-1, -1, 0)),
			newBond(B, 1, 6, getDefaultLabelZ(LB), (0, 0, -1)),
			newBond(B, 6, 1, getDefaultLabelZ(LB), (0, 0, 1)),
		]
	)
end


text:
"""
    LVmodel(N::Int64 = 20, tSpan::Vector{<:Real} = [0, 50]; β = 1, h = 0.1, random = false)

Generates a Lotka-Volterra model for N species over the given time span.
The output is a time series of abundance values that can be stored as two separate objects:
`t, Xapp = LVmodel()`

# Arguments
- `N::Int64`: the size of the model community, that is, the number of species in the system.
- `tSpan::Vector{<:Real}`: the time span along which computation is performed.
- `β::Real`: the order of derivation for all the differential equations of the LV model.
- `h::Real`: the step size for correction.
- `random::Bool`: whether the function should return a randomised (`true`) or constant (`false`) model.
"""
function LVmodel(N::Int64 = 20, tSpan::Vector{<:Real} = [0, 50]; β = 1, h = 0.1, random = false)

    # order of derivatives
    β = ones(N) * β

    if random == false

        Random.seed!(1234)

    end

    # initial abundances
    X0 = 2 * rand(N)

    if random == false

        Random.seed!(1234)

    end

    # define system of differential equations
    par = [2, 2 * rand(N), rand(N), 4 * rand(N, N), N]

    # evaluate numerical solution
    t, Xapp = FDEsolver(lv, tSpan, X0, β, par, h = h)

    return t, Xapp'

end

function lv(t, x, par)

    l = par[1] # Hill coefficient
    b = par[2] # growth rates
    k = par[3] # death rates
    K = par[4] # inhibition matrix
    N = par[5] # number of species

    F = zeros(N)

    for i in 1:N

        # inhibition functions
        f = prod(K[i, 1:end .!= i].^l ./ (K[i, 1:end .!= i] .^ l .+ x[1:end .!= i].^l))

        # ode
        F[i] = x[i] .* (b[i] .* f .- k[i] .* x[i])

    end

    return F

end

"""
    abundance_plot(se::SummarizedExperiment, assay_name::String)

Plots a time series of the abundance for each feature (row) of an assay in
a SummarizedExperiment object. The corresponding should include a column
named 'time' containing an array of sampling times

# Arguments
- `se::SummarizedExperiment`: the experiment object of interest.
- `assay_name::String`: the name of the assay to base the visualisation on.
"""
function abundance_plot(se::SummarizedExperiment, assay_name::String)

    if sum("time" .== names(coldata(se))) == 0

        error("coldata(se) should include a column named 'time' containing an array of sampling times; for example, add it with coldata(se).time = 1:10 if there are 10 samples.")

    end

    if sum(ismissing.(assay(se, assay_name))) != 0

        error("the input assay contains rows with missing values; they should be either removed or replaced with reasonable numerical values")

    end

    labels = reshape(se.rowdata.name, (1, length(se.rowdata.name)))

    p = plot(se.coldata.time, Number.(assay(se, assay_name))',
             label = labels, legend_position = :outerleft,
             xaxis = "Time", yaxis = "Abundance")

    return(p)

end

function size(se::SummarizedExperiment)

    nrow(se.rowdata), nrow(se.coldata)

end

nrow(se::SummarizedExperiment) = nrow(se.rowdata)
ncol(se::SummarizedExperiment) = nrow(se.coldata)

size(se::SummarizedExperiment, dim::Int64) = ifelse(dim == 1, nrow(se), ncol(se))


text:
using HNSW
using BenchmarkTools
using FileIO

const SUITE = BenchmarkGroup()
SUITE["build hnsw"] = BenchmarkGroup()
SUITE["knn"] = BenchmarkGroup()

for dimension ∈ (1,10,100)
    for points ∈ (1000,10000)
        data = [rand(dimension) for i=1:points]
	SUITE["build hnsw"]["init dim=$dimension, points=$points"] = @benchmarkable HierarchicalNSW($data)
        SUITE["build hnsw"]["dim=$dimension, points=$points"] = @benchmarkable add_to_graph!(HierarchicalNSW($data))
        hnsw = HierarchicalNSW(data)
        add_to_graph!(hnsw)
        for ef = (10, 100)
	        set_ef!(hnsw, ef)
            for K = (1,10)
                q = rand(dimension)
                SUITE["knn"]["K=$K nearest, ef=$ef"] = @benchmarkable knn_search($hnsw, $q, $K)
            end
        end
    end
end


text:
using Pkg; Pkg.activate(".")
include("../src/ARS.jl")
using Plots

μ, σ = 0.0, 5.0
support = (-10.0, 1.0)

f(x) = exp(-0.5(x - μ)^2 / σ^2) / sqrt(2pi * σ^2)

# Build the sampler and simulate 10,000 samples
sampler = RejectionSampler(
    f,
    support,
    (-3.0, -2.0),
    max_segments = 25,
    use_secants=true
)
@time sim = run_sampler!(sampler, 100_000);

x = range(-10.0, 10.0, length=100)
env = sampler.envelop
majorant = eval_envelop(sampler, x)
minorant = eval_secants(sampler, x)
target = f.(x)

histogram(sim, normalize = true, label = "Histogram", fillalpha=0.5)
plot!(x, [target majorant minorant], width = 5, label = ["Target" "Majorant" "Minorant"])

##

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

f2(x) = begin
    δ = x - μ
    fout = - 0.5δ^2 / σ^2
    gout = - δ / σ^2
    fout, gout
end

# Build the sampler and simulate 10,000 samples
sampler = RejectionSampler(
    f2,
    support,
    (-3.0, -2.0),
    autograd=false,
    apply_log=false,
    max_segments = 25,
    use_secants=true
)
@time sim = run_sampler!(sampler, 100_000);

x = range(-10.0, 10.0, length=100)
env = sampler.envelop
majorant = eval_envelop(sampler, x)
minorant = eval_secants(sampler, x)
target = f.(x)

histogram(sim, normalize = true, label = "Histogram", fillalpha=0.5)
plot!(x, [target majorant minorant], width = 5, label = ["Target" "Majorant" "Minorant"])


text:
using aoc2019.computer: load_program, io, run
using Match, LightGraphs

input = joinpath(@__DIR__, "input")
p = load_program(input)

struct Coord
    x::Int
    y::Int
end

mutable struct Maze
    pathing::SimpleDiGraph
    origin::Union{Coord, Nothing}
    oxygen::Union{Coord, Nothing}
    corridors::Set{Int}
    walls::Set{Int}
    ## ordered so recently discovered unknowns are visited first
    unknown::Array{Int, 1}
    vertices::Dict{Coord, Int}
end

NORTH = 1
SOUTH = 2
WEST = 3
EAST = 4

WALL = 0
OK = 1
OXY = 2
UNKNOWN = 3

function explore!(location, maze, stdin, stdout)
    loc = location
    while true
        ## add current location and neighbors to the maze
        update_maze!(maze, loc, OK)
        ## stop when there are no remaining unknown nodes
        if length(maze.unknown) == 0
            break
        end
        ## find the most recently discovered unknown location
        nearest_unknown = last(maze.unknown)
        edge = first(a_star(maze.pathing, maze.vertices[loc], nearest_unknown))
        next = first(collect(keys(filter(n -> n.second == edge.dst, maze.vertices))))
        ## try to move closer to nearest_unknown
        put!(stdin, dir(loc, next))
        res = take!(stdout)
        @match res begin
            0 => begin
                ## wall, move vertex from unknown to wall and loop
                t!(maze, edge.dst, WALL)
            end
            1 => begin
                ## corridor, move to location and loop
                loc = next
            end
            2 => begin
                ## oxygen!! move to location and update maze
                t!(maze, edge.dst, OXY)
                maze.oxygen = next
                loc = next
            end
        end
    end
end

function update_maze!(maze, location, t)
    # add or update vertex and return vertex id
    id = add_or_update_vertex!(maze, location, OK)
    # update neighbors
    for direction in [NORTH, SOUTH, WEST, EAST]
        destination = coords(direction, location)
        if !exists(maze, destination)
            # add unvisited neighbor to the maze as unknown
            add_or_update_vertex!(maze, destination, UNKNOWN)
        end
        # add edge from location to neighbor if none exists
        nid = maze.vertices[destination]
        if !has_edge(maze.pathing, id, nid)
            add_edge!(maze.pathing, id, nid)
        end
    end
end

function add_or_update_vertex!(maze, location, t)
    if !exists(maze, location)
        add_vertex!(maze.pathing)
        id = nv(maze.pathing)
        maze.vertices[location] = id
        t!(maze, id, t)
        return id
    else
        id = maze.vertices[location]
        t!(maze, id, t)
        return id
    end
end

function exists(maze, location)
    return haskey(maze.vertices, location) ? true : false
end

function t!(maze, id, t)
    @match t begin
        0 => begin
            push!(maze.walls, id)
            filter!(n -> n != id, maze.unknown)
        end
        1 => begin
            push!(maze.corridors, id)
            filter!(n -> n != id, maze.unknown)
        end
        2 => begin
            ## mark oxygen as a corridor so it's navigable
            push!(maze.corridors, id)
            filter!(n -> n != id, maze.unknown)
        end
        3 => begin
            push!(maze.unknown, id)
        end
    end
end

function coords(dir::Int, location::Coord)
    x, y = location.x, location.y
    @match dir begin
        1 => return Coord(x, y + 1)
        2 => return Coord(x, y - 1)
        3 => return Coord(x + 1, y)
        4 => return Coord(x - 1, y)
        _ => throw("unreachable")
    end
end

function dir(location::Coord, destination::Coord)
    x, y = location.x, location.y
    dx, dy = destination.x, destination.y
    @match Pair(dx - x, dy - y) begin
        Pair(0, 1) => return NORTH
        Pair(0, -1) => return SOUTH
        Pair(-1, 0) => return WEST
        Pair(1, 0) => return EAST
        _ => throw("unreachable")
    end
end

## robot location, also the origin
location = Coord(0, 0)
## the maze
maze = Maze(
    SimpleDiGraph(),
    location,
    nothing,
    Set{Int}(),
    Set{Int}(),
    Array{Int, 1}(),
    Dict{Coord, Int}()
)

# start computer
stdin, stdout = io()
@async run(copy(p), stdin, stdout)

explore!(location, maze, stdin, stdout)

pathing = maze.pathing
vertices = maze.vertices
origin = vertices[maze.origin]
oxygen = vertices[maze.oxygen]

p1 = length(a_star(pathing, origin, oxygen))
@assert p1 == 236

corridors = collect(maze.corridors)
p2 = maximum(map(v -> length(a_star(pathing, oxygen, v)), corridors))
@assert p2 == 368

print("-----------------------------------------------------------------------\n")
print("oxygen system -- part one\n    distance to leak           : $p1\n")
print("oxygen system -- part two\n    distance to furthest point : $p2\n")
print("-----------------------------------------------------------------------\n")


text:
export Likelihood

@doc raw"""
    Likelihood(M<:ParameterizedMeasure, x)

"Observe" a value `x`, yielding a function from the parameters to ℝ.

Likelihoods are most commonly used in conjunction with an existing _prior_
measure to yield a new measure, the _posterior_. In Bayes's Law, we have

``P(θ|x) ∝ P(θ) P(x|θ)``

Here ``P(θ)`` is the prior. If we consider ``P(x|θ)`` as a function on ``θ``,
then it is called a likelihood.

Since measures are most commonly manipulated using `density` and `logdensity`,
it's awkward to commit a (log-)likelihood to using one or the other. To evaluate
a `Likelihood`, we therefore use `density` or `logdensity`, depending on the
circumstances. In the latter case, it is of course acting as a log-density.

For example,

    julia> ℓ = Likelihood(Normal{(:μ,)}, 2.0)
    Likelihood(Normal{(:μ,), T} where T, 2.0)

    julia> density(ℓ, (μ=2.0,))
    1.0

    julia> logdensity(ℓ, (μ=2.0,))
    -0.0

If, as above, the measure includes the parameter information, we can optionally
leave it out of the second argument in the call to `density` or `logdensity`.

    julia> density(ℓ, 2.0)
    1.0

    julia> logdensity(ℓ, 2.0)
    -0.0

With several parameters, things work as expected:

    julia> ℓ = Likelihood(Normal{(:μ,:σ)}, 2.0)
    Likelihood(Normal{(:μ, :σ), T} where T, 2.0)

    julia> logdensity(ℓ, (μ=2, σ=3))
    -1.0986122886681098

    julia> logdensity(ℓ, (2,3))
    -1.0986122886681098

    julia> logdensity(ℓ, [2, 3])
    -1.0986122886681098

---------

    Likelihood(M<:ParameterizedMeasure, constraint::NamedTuple, x)

In some cases the measure might have several parameters, and we may want the
(log-)likelihood with respect to some subset of them. In this case, we can use
the three-argument form, where the second argument is a constraint. For example,

    julia> ℓ = Likelihood(Normal{(:μ,:σ)}, (σ=3.0,), 2.0)
    Likelihood(Normal{(:μ, :σ), T} where T, (σ = 3.0,), 2.0)

Similarly to the above, we have

    julia> density(ℓ, (μ=2.0,))
    0.3333333333333333

    julia> logdensity(ℓ, (μ=2.0,))
    -1.0986122886681098

    julia> density(ℓ, 2.0)
    0.3333333333333333

    julia> logdensity(ℓ, 2.0)
    -1.0986122886681098

-----------------------

Finally, let's return to the expression for Bayes's Law,

``P(θ|x) ∝ P(θ) P(x|θ)``

The product on the right side is computed pointwise. To work with this in
MeasureTheory, we have a "pointwise product" `⊙`, which takes a measure and a
likelihood, and returns a new measure, that is, the unnormalized posterior that has density ``P(θ) P(x|θ)`` with respect to the base measure of the prior.

For example, say we have

    μ ~ Normal()
    x ~ Normal(μ,σ)
    σ = 1

and we observe `x=3`. We can compute the posterior measure on `μ` as

    julia> post = Normal() ⊙ Likelihood(Normal{(:μ, :σ)}, (σ=1,), 3)
    Normal() ⊙ Likelihood(Normal{(:μ, :σ), T} where T, (σ = 1,), 3)

    julia> logdensity(post, 2)
    -2.5
"""
struct Likelihood{F,X}
    f::F
    x::X
end

function Base.show(io::IO, ℓ::Likelihood{Tuple{M,NamedTuple{N,T}}}) where {M<:Type, N,T}
    (m,c) = ℓ.f
    println(io, "Likelihood(",m,", ",c,", ", ℓ.x, ")")
end

function Base.show(io::IO, ℓ::Likelihood)
    println(io, "Likelihood(",ℓ.f, ", ", ℓ.x, ")")
end

Likelihood(μ::M, x) where {M<:AbstractMeasure} = Likelihood(M, x)

function Likelihood(::Type{M}, constraint::NamedTuple, x) where {M <: ParameterizedMeasure}
    Likelihood((M, constraint), x)
end

function Likelihood(μ::M, constraint::NamedTuple, x) where {M<:AbstractMeasure}
    Likelihood((M, constraint), x)
end

function logdensity(ℓ::Likelihood{Tuple{M,NT}}, p::NamedTuple) where {M<:Type, NT <: NamedTuple}
    (D, constraint) = ℓ.f
    return logdensity(D(merge(p, constraint)), ℓ.x)
end

function logdensity(ℓ::Likelihood{Tuple{M,NT}}, p) where {M<:Type, NT <: NamedTuple}
    freevar = params(ℓ.f...)
    logdensity(ℓ, NamedTuple{freevar}(p))
end

logdensity(ℓ::Likelihood, p) = logdensity(ℓ.f(p), ℓ.x)


text:
const DEBUG_FLAG_NAME = "debug"
const DEBUG_FLAG_ENABLED_VALUE = "enabled"

function get_debug_flag(value = @load_preference(DEBUG_FLAG_NAME))
    value in (DEBUG_FLAG_ENABLED_VALUE, nothing) && return value == DEBUG_FLAG_ENABLED_VALUE
    @error "Invalid value `$( repr( value ) )` for preference `$DEBUG_FLAG_NAME`."
    return false
end

enable_debug() = @set_preferences!(DEBUG_FLAG_NAME => DEBUG_FLAG_ENABLED_VALUE)
disable_debug() = @delete_preferences!(DEBUG_FLAG_NAME)

const DEBUG_ENABLED = get_debug_flag()

macro ifdebug(ex)
    if DEBUG_ENABLED
        esc(ex)
    end
end


text:
# defines structures for holding ASWING Solution data

struct StabilityDerivatives
    CXu::Float64
    CYu::Float64
    CLu::Float64
    CXb::Float64
    CYb::Float64
    CLb::Float64
    CXa::Float64
    CYa::Float64
    CLa::Float64
    Clu::Float64
    Cmu::Float64
    Cnu::Float64
    Clb::Float64
    Cmb::Float64
    Cnb::Float64
    Cla::Float64
    Cma::Float64
    Cna::Float64
    CXp::Float64
    CYp::Float64
    CLp::Float64
    CXq::Float64
    CYq::Float64
    CLq::Float64
    CXr::Float64
    CYr::Float64
    CLr::Float64
    Clp::Float64
    Cmp::Float64
    Cnp::Float64
    Clq::Float64
    Cmq::Float64
    Cnq::Float64
    Clr::Float64
    Cmr::Float64
    Cnr::Float64
    CXdF::Array{Float64,1}
    CYdF::Array{Float64,1}
    CLdF::Array{Float64,1}
    CldF::Array{Float64,1}
    CmdF::Array{Float64,1}
    CndF::Array{Float64,1}
    CXdP::Array{Float64,1}
    CYdP::Array{Float64,1}
    CLdP::Array{Float64,1}
    CldP::Array{Float64,1}
    CmdP::Array{Float64,1}
    CndP::Array{Float64,1}
    function StabilityDerivatives(f_vba::AbstractArray{<:Real,2},
        m_vba::AbstractArray{<:Real,2}, f_rot::AbstractArray{<:Real,2},
        m_rot::AbstractArray{<:Real,2}, f_flap::AbstractArray{<:Real,2},
        m_flap::AbstractArray{<:Real,2}, f_peng::AbstractArray{<:Real,2},
        m_peng::AbstractArray{<:Real,2})

        new(f_vba[1,1], f_vba[2,1], -f_vba[3,1],
        f_vba[1,2], f_vba[2,2], -f_vba[3,2],
        f_vba[1,3], f_vba[2,3], -f_vba[3,3],
        m_vba[1,1], m_vba[2,1], m_vba[3,1],
        m_vba[1,2], m_vba[2,2], m_vba[3,2],
        m_vba[1,3], m_vba[2,3], m_vba[3,3],
        f_rot[1,1], f_rot[2,1], -f_rot[3,1],
        f_rot[1,2], f_rot[2,2], -f_rot[3,2],
        f_rot[1,3], f_rot[2,3], -f_rot[3,3],
        m_rot[1,1], m_rot[2,1], m_rot[3,1],
        m_rot[1,2], m_rot[2,2], m_rot[3,2],
        m_rot[1,3], m_rot[2,3], m_rot[3,3],
        f_flap[1,:], f_flap[2,:], -f_flap[3,:],
        m_flap[1,:], m_flap[2,:], m_flap[3,:],
        f_peng[1,:], f_peng[2,:], -f_peng[3,:],
        m_peng[1,:], m_peng[2,:], m_peng[3,:])
    end
end
eltype(::Type{StabilityDerivatives}) = Float64

# contains Solution information for an operating point
struct SolutionParameters
    linear_acceleration::Array{Float64,1}
    angular_acceleration::Array{Float64,1}
    velocity::Array{Float64,1}
    rotation_rate::Array{Float64,1}
    position_aircraft_origin::Array{Float64,1}
    phi::Float64
    theta::Float64
    psi::Float64
    flap_defl_ctrl_var::Array{Float64,1}
    eng_pwr_ctrl_var::Array{Float64,1}
    err_int_Vinf::Float64
    err_int_beta::Float64
    err_int_alpha::Float64
    err_int_phi::Float64
    err_int_theta::Float64
    err_int_psi::Float64
    err_int_ROTxyz::Array{Float64,1}
    err_int_VACxyz::Array{Float64,1}
    position_aircraft_cg::Array{Float64,1}
    freestream_speed::Float64
    freestream_sideslip_angle::Float64
    freestream_incidence_angle::Float64
    indicated_V::Float64
    indicated_beta::Float64
    indicated_alpha::Float64
    drag_force::Float64
    side_force::Float64
    lift_force::Float64
    induced_drag_coefficient::Float64
    drag_coefficient::Float64
    sideforce_coefficient::Float64
    lift_coefficient::Float64
    rolling_moment_coefficient::Float64
    pitching_moment_coefficient::Float64
    yawing_moment_coefficient::Float64
    span_efficiency::Float64
    aero_force::Array{Float64,1}
    aero_moment::Array{Float64,1}
    accel_reaction_force::Array{Float64,1}
    accel_reaction_moment::Array{Float64,1}
    ground_reaction_force::Array{Float64,1}
    ground_reaction_moment::Array{Float64,1}
    engine_force::Array{Float64,1}
    engine_moment::Array{Float64,1}
    trefftz_plane_aero_force::Array{Float64,1}
    trefftz_plane_aero_moment::Array{Float64,1}
    mach::Float64
    altitude::Float64
    density::Float64
    speed_of_sound::Float64
    viscosity::Float64
    gravitational_acceleration::Float64
    load_factor::Float64
    mass::Float64
end
eltype(::Type{SolutionParameters}) = Float64

# contains spanwise input variable information interpolated to Solution nodes
struct SolutionGeometry
    s::Array{Float64,1}
    x::Array{Float64,1}
    y::Array{Float64,1}
    z::Array{Float64,1}
    twist::Array{Float64,1}
    EIcc::Array{Float64,1}
    EInn::Array{Float64,1}
    EIcn::Array{Float64,1}
    EIcs::Array{Float64,1}
    EIsn::Array{Float64,1}
    GJ::Array{Float64,1}
    EA::Array{Float64,1}
    GKc::Array{Float64,1}
    GKn::Array{Float64,1}
    mgcc::Array{Float64,1}
    mgnn::Array{Float64,1}
    mg::Array{Float64,1}
    Ccg::Array{Float64,1}
    Ncg::Array{Float64,1}
    Dmgcc::Array{Float64,1}
    Dmgnn::Array{Float64,1}
    Dmg::Array{Float64,1}
    DCcg::Array{Float64,1}
    DNcg::Array{Float64,1}
    Cea::Array{Float64,1}
    Nea::Array{Float64,1}
    Cta::Array{Float64,1}
    Nta::Array{Float64,1}
    tdeps::Array{Float64,1}
    tdgam::Array{Float64,1}
    Cshell::Array{Float64,1}
    Nshell::Array{Float64,1}
    Atshell::Array{Float64,1}
    radius::Array{Float64,1}
    Cdf::Array{Float64,1}
    Cdp::Array{Float64,1}
    chord::Array{Float64,1}
    Xax::Array{Float64,1}
    alpha::Array{Float64,1}
    Cm::Array{Float64,1}
    CLmax::Array{Float64,1}
    CLmin::Array{Float64,1}
    dCLda::Array{Float64,1}
    dCLdF::Array{Float64,2}
    dCMdF::Array{Float64,2}
    dCDdF::Array{Float64,2}
end
eltype(::Type{SolutionGeometry}) = Float64

# contains spanwise output variable information at Solution nodes
struct SolutionSpanwise
    x::Array{Float64,1}
    y::Array{Float64,1}
    z::Array{Float64,1}
    phi::Array{Float64,1}
    theta::Array{Float64,1}
    psi::Array{Float64,1}
    Mc::Array{Float64,1}
    Ms::Array{Float64,1}
    Mn::Array{Float64,1}
    Fc::Array{Float64,1}
    Fs::Array{Float64,1}
    Fn::Array{Float64,1}
    ux::Array{Float64,1}
    uy::Array{Float64,1}
    uz::Array{Float64,1}
    wx::Array{Float64,1}
    wy::Array{Float64,1}
    wz::Array{Float64,1}
    max_extensional_strain::Array{Float64,1}
    shear_stress::Array{Float64,1}
    chord::Array{Float64,1}
    cl::Array{Float64,1}
    cm::Array{Float64,1}
    f_aero::Array{Float64,1}
    f_reac::Array{Float64,1}
    da_eff::Array{Float64,1} # aeroelastic twist
    da_flp::Array{Float64,1} # additional twist due to flaps
end
eltype(::Type{SolutionSpanwise}) = Float64

# contains aswing solution data
struct Solution
    converged::Bool
    pnt_param::SolutionParameters
    span_geom::Array{SolutionGeometry,1}
    span_param::Array{SolutionSpanwise,1}
end
eltype(::Type{Solution}) = Float64

# converts ASWING output arrays to Solution struct
function getsolution(ipnt1::Integer=1, ipnt2::Integer=ASWING.NPOINT[1])

    sol = Array{Solution,1}(undef, max(ipnt2-ipnt1+1,0))
    for ipnt = ipnt1:ipnt2
        # Process converged boolean
        converged = ASWING.LCONV[ipnt]
        # Process operating point parameters
        linear_acceleration = ASWING.PARAM[KPVACX:KPVACZ,ipnt]
        angular_acceleration = ASWING.PARAM[KPRACX:KPRACZ,ipnt]
        velocity = ASWING.PARAM[KPVELX:KPVELZ,ipnt]
        rotation_rate = ASWING.PARAM[KPROTX:KPROTZ,ipnt]
        position_aircraft_origin = ASWING.PARAM[KPPOSX:KPPOSZ,ipnt]
        phi = ASWING.PARAM[KPBANK,ipnt]
        theta = ASWING.PARAM[KPELEV,ipnt]
        psi = ASWING.PARAM[KPHEAD,ipnt]
        flap_defl_ctrl_var = ASWING.PARAM[KPFLP1:KPFLP1+NFLPX-1,ipnt]
        eng_pwr_ctrl_var = ASWING.PARAM[KPENG1:KPENG1+NENGX-1,ipnt]
        err_int_Vinf = ASWING.PARAM[KPVIDT,ipnt]
        err_int_beta = ASWING.PARAM[KPBEDT,ipnt]
        err_int_alpha =ASWING.PARAM[KPALDT,ipnt]
        err_int_phi = ASWING.PARAM[KPBADT,ipnt]
        err_int_theta = ASWING.PARAM[KPELDT,ipnt]
        err_int_psi = ASWING.PARAM[KPHEDT,ipnt]
        err_int_ROTxyz = ASWING.PARAM[KPWDTX:KPWDTZ,ipnt]
        err_int_VACxyz = ASWING.PARAM[KPADTX:KPADTZ,ipnt]
        position_aircraft_cg = ASWING.PARAM[KPRCGX:KPRCGZ,ipnt]
        freestream_speed = ASWING.PARAM[KPVINF,ipnt]
        freestream_sideslip_angle = ASWING.PARAM[KPBETA,ipnt]
        freestream_incidence_angle = ASWING.PARAM[KPALFA,ipnt]
        indicated_V = ASWING.PARAM[KPVIAS,ipnt]
        indicated_beta = ASWING.PARAM[KPBREF,ipnt]
        indicated_alpha = ASWING.PARAM[KPAREF,ipnt]
        drag_force = ASWING.PARAM[KPDRAG,ipnt]
        side_force = ASWING.PARAM[KPSIDE,ipnt]
        lift_force = ASWING.PARAM[KPLIFT,ipnt]
        induced_drag_coefficient = ASWING.PARAM[KPCDIN,ipnt]
        drag_coefficient = ASWING.PARAM[KPCDRG,ipnt]
        sideforce_coefficient = ASWING.PARAM[KPCSID,ipnt]
        lift_coefficient = ASWING.PARAM[KPCLIF,ipnt]
        rolling_moment_coefficient = ASWING.PARAM[KPCROL,ipnt]
        pitching_moment_coefficient = ASWING.PARAM[KPCMOM,ipnt]
        yawing_moment_coefficient = ASWING.PARAM[KPCYAW,ipnt]
        span_efficiency = ASWING.PARAM[KPSEFF,ipnt]
        aero_force = ASWING.PARAM[KPAFOX:KPAFOZ,ipnt]
        aero_moment = ASWING.PARAM[KPAMOX:KPAMOZ,ipnt]
        accel_reaction_force = ASWING.PARAM[KPRFOX:KPRFOZ,ipnt]
        accel_reaction_moment = ASWING.PARAM[KPRMOX:KPRMOZ,ipnt]
        ground_reaction_force = ASWING.PARAM[KPGFOX:KPGFOZ,ipnt]
        ground_reaction_moment = ASWING.PARAM[KPGMOX:KPGMOZ,ipnt]
        engine_force = ASWING.PARAM[KPEFOX:KPEFOZ,ipnt]
        engine_moment = ASWING.PARAM[KPEMOX:KPEMOZ,ipnt]
        trefftz_plane_aero_force = ASWING.PARAM[KPTFOX:KPTFOZ,ipnt]
        trefftz_plane_aero_moment = ASWING.PARAM[KPTMOX:KPTMOZ,ipnt]
        mach = ASWING.PARAM[KPMACH,ipnt]
        altitude = ASWING.PARAM[KPALTK,ipnt]
        density = ASWING.PARAM[KPDENS,ipnt]
        speed_of_sound = ASWING.PARAM[KPVSOU,ipnt]
        viscosity = ASWING.PARAM[KPVISC,ipnt]
        gravitational_acceleration = ASWING.PARAM[KPGRAV,ipnt]
        load_factor = ASWING.PARAM[KPLOAD,ipnt]
        pnt_param = SolutionParameters(linear_acceleration,angular_acceleration,velocity,rotation_rate,
        position_aircraft_origin,phi,theta,psi,flap_defl_ctrl_var,eng_pwr_ctrl_var,
        err_int_Vinf,err_int_beta,err_int_alpha,err_int_phi,err_int_theta,err_int_psi,
        err_int_ROTxyz,err_int_VACxyz,position_aircraft_cg,freestream_speed,
        freestream_sideslip_angle,freestream_incidence_angle,indicated_V,
        indicated_beta,indicated_alpha,drag_force,side_force,lift_force,
        induced_drag_coefficient,drag_coefficient,sideforce_coefficient,
        lift_coefficient,rolling_moment_coefficient,pitching_moment_coefficient,
        yawing_moment_coefficient,span_efficiency,aero_force,aero_moment,
        accel_reaction_force,accel_reaction_moment,ground_reaction_force,
        ground_reaction_moment,engine_force,engine_moment,trefftz_plane_aero_force,
        trefftz_plane_aero_moment,mach,altitude,density,speed_of_sound,viscosity,
        gravitational_acceleration,load_factor,ASWING.MASS[1])
        #Process spanwise geometry and Solution
        span_geom = Array{SolutionGeometry,1}(undef, ASWING.NBEAM[1])
        span_sol = Array{SolutionSpanwise,1}(undef, ASWING.NBEAM[1])
        for ibeam = 1:ASWING.NBEAM[1]
            i1 = convert(Int,ASWING.IFRST[ibeam])
            i2 = convert(Int,ASWING.ILAST[ibeam])
            # spanwise geometry
            sgeom = ASWING.S[i1:i2]
            x = ASWING.QPAR[i1:i2,JXA]
            y = ASWING.QPAR[i1:i2,JYA]
            z = ASWING.QPAR[i1:i2,JZA]
            twist = ASWING.QPAR[i1:i2,JTW]
            EIcc = ASWING.QPAR[i1:i2,JECC]
            EInn = ASWING.QPAR[i1:i2,JENN]
            EIcn = ASWING.QPAR[i1:i2,JECN]
            EIcs = ASWING.QPAR[i1:i2,JECS]
            EIsn = ASWING.QPAR[i1:i2,JESN]
            GJ = ASWING.QPAR[i1:i2,JGJ]
            EA = ASWING.QPAR[i1:i2,JEA]
            GKc = ASWING.QPAR[i1:i2,JGKC]
            GKn = ASWING.QPAR[i1:i2,JGKN]
            mgcc = ASWING.QPAR[i1:i2,JMCC1]
            mgnn = ASWING.QPAR[i1:i2,JMNN1]
            mg = ASWING.QPAR[i1:i2,JMG1]
            Ccg = ASWING.QPAR[i1:i2,JCCG1]
            Ncg = ASWING.QPAR[i1:i2,JNCG1]
            Dmgcc = ASWING.QPAR[i1:i2,JMCC2]
            Dmgnn = ASWING.QPAR[i1:i2,JMNN2]
            Dmg = ASWING.QPAR[i1:i2,JMG2]
            DCcg = ASWING.QPAR[i1:i2,JCCG2]
            DNcg = ASWING.QPAR[i1:i2,JNCG2]
            Cea = ASWING.QPAR[i1:i2,JCEA]
            Nea = ASWING.QPAR[i1:i2,JNEA]
            Cta = ASWING.QPAR[i1:i2,JCTA]
            Nta = ASWING.QPAR[i1:i2,JNTA]
            tdeps = ASWING.QPAR[i1:i2,JTDE]
            tdgam = ASWING.QPAR[i1:i2,JTDG]
            Cshell = ASWING.QPAR[i1:i2,JCSH]
            Nshell = ASWING.QPAR[i1:i2,JNSH]
            Atshell = ASWING.QPAR[i1:i2,JASH]
            radius = ASWING.QPAR[i1:i2,JRAD]
            Cdf = ASWING.QPAR[i1:i2,JCDF]
            Cdp = ASWING.QPAR[i1:i2,JCDP]
            chord = ASWING.QPAR[i1:i2,JCH]
            Xax = ASWING.QPAR[i1:i2,JXAX]
            alpha = ASWING.QPAR[i1:i2,JAL]
            Cm = ASWING.QPAR[i1:i2,JCM]
            CLmax = ASWING.QPAR[i1:i2,JCMX]
            CLmin = ASWING.QPAR[i1:i2,JCMN]
            dCLda = ASWING.QPAR[i1:i2,JDCL]
            dCLdF = ASWING.QPAR[i1:i2,JCLF1:JCLF1+NFLPX-1]
            dCMdF = ASWING.QPAR[i1:i2,JCMF1:JCMF1+NFLPX-1]
            dCDdF = ASWING.QPAR[i1:i2,JCDF1:JCDF1+NFLPX-1]
            span_geom[ibeam] = SolutionGeometry(sgeom,x,y,z,twist,EIcc,EInn,EIcn,EIcs,EIsn,GJ,
            EA,GKc,GKn,mgcc,mgnn,mg,Ccg,Ncg,Dmgcc,Dmgnn,Dmg,DCcg,DNcg,Cea,Nea,Cta,Nta,
            tdeps,tdgam,Cshell,Nshell,Atshell,radius,Cdf,Cdp,chord,Xax,alpha,Cm,CLmax,
            CLmin,dCLda,dCLdF, dCMdF, dCDdF)
            # spanwise Solution
            x=ASWING.QPNT[i1:i2,1,ipnt]
            y=ASWING.QPNT[i1:i2,2,ipnt]
            z=ASWING.QPNT[i1:i2,3,ipnt]
            psi=ASWING.QPNT[i1:i2,4,ipnt]
            theta=ASWING.QPNT[i1:i2,5,ipnt]
            phi=ASWING.QPNT[i1:i2,6,ipnt]
            Mc=ASWING.QPNT[i1:i2,7,ipnt]
            Ms=ASWING.QPNT[i1:i2,8,ipnt]
            Mn=ASWING.QPNT[i1:i2,9,ipnt]
            Fc=ASWING.QPNT[i1:i2,10,ipnt]
            Fs=ASWING.QPNT[i1:i2,11,ipnt]
            Fn=ASWING.QPNT[i1:i2,12,ipnt]
            ux=ASWING.QPNT[i1:i2,13,ipnt]
            uy=ASWING.QPNT[i1:i2,14,ipnt]
            uz=ASWING.QPNT[i1:i2,15,ipnt]
            wx=ASWING.QPNT[i1:i2,16,ipnt]
            wy=ASWING.QPNT[i1:i2,17,ipnt]
            wz=ASWING.QPNT[i1:i2,18,ipnt]
            max_extensional_strain=ASWING.QPNT[i1:i2,19,ipnt]
            shear_stress=ASWING.QPNT[i1:i2,20,ipnt]
            chord=ASWING.QPNT[i1:(i2-1),21,ipnt]
            cl=ASWING.QPNT[i1:(i2-1),22,ipnt]
            cm=ASWING.QPNT[i1:(i2-1),23,ipnt]
            f_aero=ASWING.QPNT[i1:(i2-1),24,ipnt]
            f_reac=ASWING.QPNT[i1:(i2-1),25,ipnt]
            da_eff=ASWING.QPNT[i1:(i2-1),26,ipnt]
            da_flp=ASWING.QPNT[i1:(i2-1),27,ipnt]
            span_sol[ibeam] = SolutionSpanwise(x,y,z,psi,theta,phi,Mc,Ms,Mn,Fc,Fs,Fn,
            ux,uy,uz,wx,wy,wz,max_extensional_strain,shear_stress,chord,cl,cm,f_aero,
            f_reac,da_eff,da_flp)
        end

        # Create combined Solution
        sol[ipnt] = Solution(converged,pnt_param,span_geom,span_sol)
    end
    return sol
end


text:
module DualOptimalFiltering

export dirichletkernel, filter_CIR, filter_WF, generate_CIR_trajectory

include("mcmc_sampler.jl")
include("common_utility_functions.jl")
include("pruning_functions.jl")
include("general_smoothing_functions.jl")
include("kde_for_pf_samples.jl")
include("dirichlet_kde.jl")
include("exact_L2_distances.jl")
include("post_process_Dirichlet_mixture_posterior.jl")

include("functions_CIR.jl")
include("filtering_CIR.jl")
include("approx_filtering_CIR.jl")
include("likelihood_CIR.jl")
include("smoothing_CIR.jl")
include("approx_smoothing_CIR.jl")
include("joint_smoothing_CIR.jl")
include("reparam_CIR.jl")
include("full_inference_CIR.jl")

include("functions_WF.jl")
include("filtering_WF.jl")
include("approx_filtering_WF.jl")
include("likelihood_WF.jl")
include("particle_filtering_WF.jl")
end # module

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
function Integrators.get_internal_variables(int::AbstractIntegratorSPARK{DT,TT,D,S,R}) where {DT, TT, D, S, R}
    Qi = create_internal_stage_vector(DT, D, S)
    Pi = create_internal_stage_vector(DT, D, S)
    Vi = create_internal_stage_vector(DT, D, S)
    Φi = create_internal_stage_vector(DT, D, S)

    Qp = create_internal_stage_vector(DT, D, R)
    Pp = create_internal_stage_vector(DT, D, R)
    Λp = create_internal_stage_vector(DT, D, R)
    Φp = create_internal_stage_vector(DT, D, R)

    solver = get_solver_status(int.solver)

    (Qi=Qi, Pi=Pi, Vi=Vi, Φi=Φi, Qp=Qp, Pp=Pp, Λp=Λp, Φp=Φp, solver=solver)
end

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

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

function Integrators.integrate_step!(int::AbstractIntegratorSPARK{DT,TT}, sol::AtomicSolutionPDAE{DT,TT},
                                     cache::IntegratorCacheSPARK{DT}=int.caches[DT]) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int.params, sol)

    # compute initial guess
    initial_guess!(int, sol, cache)

    # reset cache
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # check_jacobian(int.solver)
    # print_jacobian(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector fields at internal stages
    compute_stages!(int.solver.x, cache, int.params)

    # compute final update
    update_solution!(int, sol, cache)

    # copy internal stage variables
    sol.internal.Qi .= cache.Qi
    sol.internal.Pi .= cache.Pi
    sol.internal.Vi .= cache.Vi
    sol.internal.Φi .= cache.Φi

    sol.internal.Qp .= cache.Qp
    sol.internal.Pp .= cache.Pp
    sol.internal.Λp .= cache.Λp
    sol.internal.Φp .= cache.Φp

    # copy solution to initial guess
    update_vector_fields!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f)
end


text:
"""
	grdcontour(cmd0::String="", arg1=nothing; kwargs...)

Reads a 2-D grid file or a GMTgrid type and produces a contour map by tracing each
contour through the grid.

Full option list at [`grdcontour`]($(GMTdoc)grdcontour.html)

Parameters
----------

- $(GMT.opt_J)
- **A** | **annot** :: [Type => Str or Number]       ``Arg = [-|[+]annot_int][labelinfo]``

    *annot_int* is annotation interval in data units; it is ignored if contour levels are given in a file.
    ($(GMTdoc)grdcontour.html#a)
- $(GMT.opt_B)
- **C** | **cont** | **contour** | **contours** | **levels** :: [Type => Str | Number | GMTcpt]  ``Arg = [+]cont_int``

    Contours to be drawn may be specified in one of three possible ways.
    ($(GMTdoc)grdcontour.html#c)
- **D** | **dump** :: [Type => Str]

    Dump contours as data line segments; no plotting takes place.
    ($(GMTdoc)grdcontour.html#d)
- **F** | **force** :: [Type => Str | []]

    Force dumped contours to be oriented so that higher z-values are to the left (-Fl [Default]) or right.
    ($(GMTdoc)grdcontour.html#f)
- **G** | **labels** :: [Type => Str]

    Controls the placement of labels along the quoted lines.
    ($(GMTdoc)grdcontour.html#g)
- $(GMT.opt_Jz)
- **L** | **range** :: [Type => Str]

    Limit range: Do not draw contours for data values below low or above high.
    ($(GMTdoc)grdcontour.html#l)
- **N** | **fill** | **colorize** :: [Type => Bool]

    Fill the area between contours using the discrete color table given by cpt.
    ($(GMTdoc)grdcontour.html#n)
- $(GMT.opt_P)
- **Q** | **cut** :: [Type => Str | Number]

    Do not draw contours with less than cut number of points.
    ($(GMTdoc)grdcontour.html#q)
- **S** | **smooth** :: [Type => Number]

    Used to resample the contour lines at roughly every (gridbox_size/smoothfactor) interval.
    ($(GMTdoc)grdcontour.html#s)
- **T** | **ticks** :: [Type => Str]

    Draw tick marks pointing in the downward direction every *gap* along the innermost closed contours.
    ($(GMTdoc)grdcontour.html#t)
- $(GMT.opt_R)
- $(GMT.opt_U)
- $(GMT.opt_V)
- **W** | **pen** :: [Type => Str | Number]

    Sets the attributes for the particular line.
    ($(GMTdoc)grdcontour.html#w)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- **Z** | **scale** :: [Type => Str]

    Use to subtract shift from the data and multiply the results by factor before contouring starts.
    ($(GMTdoc)grdcontour.html#z)
- $(GMT.opt_bo)
- $(GMT.opt_do)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_h)
- $(GMT.opt_p)
- $(GMT.opt_t)
"""
function grdcontour(cmd0::String="", arg1=nothing; first=true, kwargs...)

	length(kwargs) == 0 && cmd0 != "" && return monolitic("grdcontour", cmd0, arg1)
	arg2 = arg3 = nothing

	d, K, O = init_module(first, kwargs...)		# Also checks if the user wants ONLY the HELP mode
	dict_auto_add!(d)					# The ternary module may send options via another channel

	cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, " -JX" * split(def_fig_size, '/')[1] * "/0")
	cmd, = parse_common_opts(d, cmd, [:UVXY :params :bo :c :e :f :h :p :t], first)
	cmd  = parse_these_opts(cmd, d, [[:D :dump], [:F :force], [:L :range], [:Q :cut], [:S :smooth]])
	cmd  = parse_contour_AGTW(d::Dict, cmd::String)[1]
	cmd  = add_opt(d, cmd, 'Z', [:Z :scale], (factor = "+s", shift = "+o", periodic = "_+p"))

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, arg1)	# Find how data was transmitted
	if (isa(arg1, Array{<:Real}))		arg1 = mat2grid(arg1)	end

	# cmd, N_used, arg1, arg2, = get_cpt_set_R(d, cmd0, cmd, opt_R, got_fname, arg1, arg2, nothing, "grdcontour")
	cmd, N_used, arg1, arg2, = common_get_R_cpt(d, cmd0, cmd, opt_R, got_fname, arg1, arg2, nothing, "grdcontour")

	got_N_cpt = false		# Shits because 6.1 still cannot process N=cpt (6.1.1 can)
	if ((val = find_in_dict(d, [:N :fill :colorize], false)[1]) !== nothing)
		if (isa(val, GMTcpt))
			N = (N_used > 1) ? 1 : N_used		# Trickery because add_opt_cpt() is not able to deal with 3 argX
			if (isa(arg1, GMTgrid))
				cmd, arg2, arg3, = add_opt_cpt(d, cmd, [:N :fill :colorize], 'N', N, arg2, arg3)
    			else
				cmd, arg1, arg2, = add_opt_cpt(d, cmd, [:N :fill :colorize], 'N', N, arg1, arg2)
			end
			got_N_cpt = true
		else
			cmd *= " -N"
			del_from_dict(d, [:N :fill :colorize])
		end
	end

	if (!occursin(" -C", cmd))			# Otherwise ignore an eventual :cont because we already have it
		cmd, args, n, = add_opt(d, cmd, 'C', [:C :cont :contour :contours :levels], :data, Array{Any,1}([arg1, arg2, arg3]), (x = "",))
		if (n > 0)
			for k = 3:-1:1
				(args[k] === nothing) && continue
				if (isa(args[k], Array{<:Real}))
					cmd *= arg2str(args[k], ',')
					if (length(args[k]) == 1)  cmd *= ","  end		# A single contour needs to end with a ","
					break
				elseif (isa(args[k], GMTcpt))
					arg1, arg2, arg3 = args[:]
					break
				end
			end
		end
	end

	# -N option is bugged up to 6.1.0 because it lacked the apropriate KEY, so trickery is needed.
	if (occursin(" -N", cmd))
		if (occursin(" -C", cmd) && (isa(arg1, GMTcpt) || isa(arg2, GMTcpt)))
			if (!got_N_cpt)		# C=cpt, N=true. Must replicate the CPT into N
				isa(arg1, GMTcpt) ? arg2 = arg1 : arg3 = arg2
			end
		end
	elseif (GMTver <= v"6.1" && got_N_cpt && !occursin(" -C", cmd))	# N=cpt and no C. Work around the bug
		d[:C] = isa(arg1, GMTcpt) ? arg1 : arg2
	end

	opt_extra = "";		do_finish = true
	if (occursin("-D", cmd))
		opt_extra = "-D";		do_finish = false;	cmd = replace(cmd, opt_J => "")
	end

	_cmd, K = finish_PS_nested(d, ["grdcontour " * cmd], K)
	finish_PS_module(d, _cmd, opt_extra, K, O, do_finish, arg1, arg2, arg3)
end

# ---------------------------------------------------------------------------------------------------
function parse_contour_AGTW(d::Dict, cmd::String)
	# Common to both grd and ps contour
	if ((val = find_in_dict(d, [:A :annot], false)[1]) !== nothing && isa(val, Array{<:Real}))
		cmd *= " -A" * arg2str(val, ',')
		if (!occursin(",", cmd))  cmd *= ","  end
		del_from_dict(d, [:A :annot])
	elseif (isa(val, String) || isa(val, Symbol))
		arg = string(val)
		cmd *= (arg == "none") ? " -A-" : " -A" * arg
		del_from_dict(d, [:A :annot])
	else
		cmd = add_opt(d, cmd, 'A', [:A :annot],
		              (disable = ("_-", nothing, 1), none = ("_-", nothing, 1), single = ("+", nothing, 1), int = "", interval = "", labels = ("", parse_quoted)) )
	end
	cmd = add_opt(d, cmd, 'G', [:G :labels], ("", helper_decorated))
	cmd = add_opt(d, cmd, 'T', [:T :ticks], (local_high = ("h", nothing, 1), local_low = ("l", nothing, 1),
	                                         labels = "+l", closed = "_+a", gap = "+d") )
	opt_W = add_opt_pen(d, [:W :pen], "W", true)    # TRUE to also seek (lw,lc,ls)
	return cmd * opt_W, opt_W
end

# ---------------------------------------------------------------------------------------------------
grdcontour!(cmd0::String="", arg1=nothing; kw...) = grdcontour(cmd0, arg1; first=false, kw...)
grdcontour(arg1, cmd0::String=""; kw...) = grdcontour(cmd0, arg1; first=true, kw...)
grdcontour!(arg1, cmd0::String=""; kw...) = grdcontour(cmd0, arg1; first=false, kw...)


text:
using LSODA
using PiecewiseDeterministicMarkovProcesses, Plots, LinearAlgebra, Random

function n∞(v,v₃,v₄)
  ξ=(v-v₃)/v₄
  (1+tanh(ξ))/2
end

function τ(v,v₃,v₄,ϕ)
  ξ=(v-v₃)/v₄
  1/(ϕ*cosh(ξ/2))
end;

function m∞(v,v₁,v₂)
  (1+tanh((v-v₁)/v₂))/2
end

function α(v,v₃,v₄,ϕ)
    ξ=(v-v₃)/v₄
    ϕ*cosh(ξ/2)/(1+exp(-2*ξ))
end

function β(v,v₃,v₄,ϕ)
    ξ=(v-v₃)/v₄
    ϕ*cosh(ξ/2)/(1+exp(2*ξ))
end;

function f_ml_chv!(tt,x,xdot,data)
    (I,C,gL,gCa,gK,vL,vCa,vK,v₁,v₂,v₃,v₄,ϕ,Ntot,N) = data
    (v,t) = x
    tr = α(v,v₃,v₄,ϕ)*(Ntot-N)+β(v,v₃,v₄,ϕ)*N
    xdot[1] = (I-gL*(v-vL)-gCa*m∞(v,v₁,v₂)*(v-vCa)-gK*(N/Ntot)*(v-vK))/C/tr
    xdot[2] = 1.0/tr
    nothing
end;

function ml_chv(x0,parms,tf;n_jumps=500)
    (v,N) = x0
    (I,C,gL,gCa,gK,vL,vCa,vK,v₁,v₂,v₃,v₄,ϕ,Ntot) = parms
    eparms = [parms;N]
    t=0.0
    ta = Vector{Float64}()
    va = Vector{Float64}()
    Na = Vector{Float64}()
    push!(ta,t)
    push!(va,v)
    push!(Na,N)
    Flow(v_,t_,s_,eparms_) = LSODA.lsoda((tt,x,xdot,data)->f_ml_chv!(tt,x,xdot,eparms_),
        [v_;t_],
        [0.0,s_],
        abstol=1e-9,#Vector([1.e-10,1.e-6]),
        reltol=1e-7,
        nbsteps=10000)
    n = 1 #number of jumps, fairer to compare with PDMP
    while t<tf && n<n_jumps
        s = -log(rand())
        # res = LSODA.lsoda((tt,x,xdot,data)->f_ml_chv!(tt,x,xdot,eparms),
        #     [v;t],
        #     [0.0,s],
        #     abstol=1e-9,#Vector([1.e-10,1.e-6]),
        #     reltol=1e-7,
        #     nbsteps=10000)
        res = Flow(v,t,s,eparms)
        v = res[end,1]
        t = res[end,end]
        # Update N
        opn = α(v,v₃,v₄,ϕ)*(Ntot-N)
        cls = β(v,v₃,v₄,ϕ)*N
        p=opn/(opn+cls)
        if rand()<p
            N=N+1
            eparms[end]=N
        else
            N=N-1
            eparms[end]=N
        end
        push!(ta,t)
        push!(va,v)
        push!(Na,N)
        n += 1
    end
    return(ta,va,Na)
end

parms_chv=[100.0,20.0,2.0,4.4,8.0,-60.0,120.0,-84.0,-1.2,18.0,2.0,30,0.04,40]
x0_chv=[-50.0;20.0]
tf_chv=100000.;

Random.seed!(123)
sol_chv=ml_chv(x0_chv,parms_chv,tf_chv,n_jumps=660)
plot(sol_chv[1],sol_chv[2])

println("="^70)
Random.seed!(123)
@time begin
    for i in 1:100
        out=ml_chv(x0_chv,parms_chv,tf_chv,n_jumps=660)
    end
end

################################################################################
################################################################################
################################################################################

function f_ml_pdmp!(xcdot, xc, xd, t, parms)
  (I,C,gL,gCa,gK,vL,vCa,vK,v₁,v₂,v₃,v₄,ϕ,Ntot) = parms
  (v,) = xc
  (N,) = xd
  xcdot[1] = (I-gL*(v-vL)-gCa*m∞(v,v₁,v₂)*(v-vCa)-gK*(N/Ntot)*(v-vK))/C
  nothing
end

function r_ml_pdmp!(rate, xc, xd, t, parms, sum_rate)
  (I,C,gL,gCa,gK,vL,vCa,vK,v₁,v₂,v₃,v₄,ϕ,Ntot) = parms
  (v,) = xc
  (N,) = xd
  if sum_rate==false
      rate[1] = α(v,v₃,v₄,ϕ)*(Ntot-N)
      rate[2] = β(v,v₃,v₄,ϕ)*N
      return 0.
  else
      return α(v,v₃,v₄,ϕ)*(Ntot-N)+β(v,v₃,v₄,ϕ)*N
  end
end

rate_ = zeros(2)
xc0 = [-50.0]
xd0 = [20]
xd0 |> typeof |> println
nu_ml = reshape([[1];[-1]],2,1)
parms_chv_pdmp = [100.0,20.0,2.0,4.4,8.0,-60.0,120.0,-84.0,-1.2,18.0,2.0,30,0.04,40]
tf_pdmp = 100000.;

Random.seed!(123)
# sol_chv_pdmp=PDMP.pdmp(xc0, xd0, f_ml_pdmp!, r_ml_pdmp!, nu_ml, parms_chv_pdmp, 0.0, tf_pdmp, false, ode=:lsoda, n_jumps = 500);
sol_chv_pdmp=PDMP.chv!(660, xc0, xd0, f_ml_pdmp!, r_ml_pdmp!,PDMP.Delta_dummy, nu_ml, parms_chv_pdmp, 0.0, tf_pdmp, ode=:lsoda)
plot(sol_chv_pdmp.time,sol_chv_pdmp.xc[1,:]-sol_chv[2])

# 1.022143 seconds (8.14 M allocations: 443.023 MiB, 12.90% gc time)
# 1.072882 seconds (8.35 M allocations: 445.167 MiB, 12.10% gc time)
# v0.7
# 0.832310 seconds (11.66 M allocations: 458.525 MiB, 11.06% gc time)
# 0.832129 seconds (12.07 M allocations: 464.681 MiB, 11.09% gc time)
Random.seed!(123)
@time begin
    for i in 1:100
        # PDMP.pdmp(xc0, xd0, f_ml_pdmp!, r_ml_pdmp!, nu_ml, parms_chv_pdmp, 0.0, tf_pdmp, false, ode=:lsoda, n_jumps = 500);
        res_pdmp = PDMP.chv!(660, xc0, xd0, f_ml_pdmp!, r_ml_pdmp!,PDMP.Delta_dummy, nu_ml, parms_chv_pdmp, 0.0, tf_pdmp, ode=:lsoda);
    end
end

################################################################################
################################################################################
################################################################################
xd0chv = copy(x0_chv)
eparms = [parms_chv;20]
@time(f_ml_chv!(1.0,x0_chv,xd0chv,eparms))
    println("--> CHV result xdot = $xd0chv for x = $x0_chv")
    xc1 = [xc0;1.]
    xdotpdmp=copy(xc1)
    @time (PDMP.f_CHV!(f_ml_pdmp!, r_ml_pdmp!,0.,xc1,xdotpdmp,xd0,parms_chv_pdmp))
    println("-->PDMP result xdot = $xdotpdmp for x = $xc1, xd = $xd0")
    f_pdmp = @time  (tt,x,xdot)->PDMP.f_CHV!(f_ml_pdmp!, r_ml_pdmp!,tt,x,xdot,xd0,parms_chv_pdmp)
    @time f_pdmp(0.,xc1,xdotpdmp)
    f_chv = @time (tt,x,xdot,data)->f_ml_chv!(tt,x,xdot,eparms)
    @time f_chv(0.,xc1,xdotpdmp,eparms)


text:
const INFECTION_HORIZON = 14

const INFECTION_DATA_PATH = joinpath(@__DIR__, "data", "Sample50.csv")

const INF_DATAFRAME = begin
    dlm = readdlm(INFECTION_DATA_PATH, ',')
    DataFrame(dlm[2:end, :], string.(Int.(dlm[1,:])))
end

const INF_DIST = FitInfectionDistributions(INF_DATAFRAME, INFECTION_HORIZON, 50)

const VIRAL_LOADS_DATA_PATH = joinpath(@__DIR__, "data", "raw_viral_load.csv")

const VIRAL_LOADS = begin
    dlm = readdlm(VIRAL_LOADS_DATA_PATH, ',')
    DataFrame(dlm[2:end, :], string.(Int.(dlm[1,:])))
end

const LIMIT_OF_DETECTION = 6

const POS_TEST_PROBS = [
    prop_above_LOD(VIRAL_LOADS, day, LIMIT_OF_DETECTION)
    for day in 1:INFECTION_HORIZON
]

const SYMPTOMATIC_ISOLATION_PROB = 0.30

const ASYMPTOMATIC_PROB = 0.10

const SYMPTOM_DIST = LogNormal(1.644, 0.363)

const SYMPTOM_PROBS = [
    cdf(SYMPTOM_DIST, i) - cdf(SYMPTOM_DIST, i-1) for i in 1:INFECTION_HORIZON
]


text:
using CUDAdrv, CUDAnative, CuArrays

using Test

"Dummy kernel doing 100 FMAs."
function kernel_100fma(a, b, c, out)
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x
    @inbounds a_val = a[i]
    @inbounds b_val = b[i]
    @inbounds c_val = c[i]

    for j in 1:33
        a_val = CUDAnative.fma(a_val, b_val, c_val)
        b_val = CUDAnative.fma(a_val, b_val, c_val)
        c_val = CUDAnative.fma(a_val, b_val, c_val)
    end

    @inbounds out[i] = CUDAnative.fma(a_val, b_val, c_val)

    return
end

function peakflops(n::Integer=5000, dev::CuDevice=CuDevice(0))
    ctx = CuContext(dev)

    dims = (n, n)
    a = round.(rand(Float32, dims) * 100)
    b = round.(rand(Float32, dims) * 100)
    c = round.(rand(Float32, dims) * 100)
    out = similar(a)

    d_a = CuArray(a)
    d_b = CuArray(b)
    d_c = CuArray(c)
    d_out = CuArray(out)

    len = prod(dims)
    threads = min(len, 1024)
    blocks = len ÷ threads

    # warm-up
    @cuda kernel_100fma(d_a, d_b, d_c, d_out)
    synchronize(ctx)

    secs = CUDAdrv.@elapsed begin
        @cuda blocks=blocks threads=threads kernel_100fma(d_a, d_b, d_c, d_out)
    end
    flopcount = 200*len
    flops = flopcount / secs

    destroy!(ctx)
    return flops
end

println(peakflops())


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>20.0, :alg=>GLGM06(δ=0.1)),
    "X0"   => Hyperrectangle(low=fill(0.9, 15), high=fill(1.1, 15)),
    "info" => Dict(
        "name" => "vehicle_platoon_5",
        "dim" => 15,
        "linear" => true,
        "hybrid" => false
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text: Genie.secret_token!("e425d80f1e7d6ea38a1f2517ebc6104f553668bfd0c824813907a00116e820e0")


text:
"""
substring(source)

converts to a SubString{String}, trying to use an internal
string buffer. Works without content copying on Utf8String and IOShared objects

substring is intended to be used with (offset,size,source) tuples as argument on
byte buffers, not the String-oriented (source,firstIndex,lastIndex) notation introduced in
julia Base for AbstractString sources.

However, in all cases, UTF8 validity is checked (that is built into SubString constructor).
Only use substring for Utf8 string processing!!

"""
function substring end

"an empty string buffer"
const EMPTYSTRING = ""

usize(s::Union{String,SubString{String}}) = ncodeunits(s)%UInt64

"SubString construction using offset and size in bytes, with check of UTF8 content validity. "
Base.@propagate_inbounds function Base.SubString(ofs::UInt32, size::UInt64,s::String)
#Base.@propagate_inbounds function Base.SubString{String}(ofs::UInt32, size::UInt64,s::String)
    @boundscheck checkulimit(ofs+size,ncodeunits(s)%UInt64)
    return @inbounds SubString(s,ofs+1,thisind(s,(ofs+size)%Int))
end

substring(s::AbstractString) = substring(s,firstindex(s),lastindex(s))

substring(s::String,first::Integer,last::Integer) = SubString(s,first,last)

"expensive default (converts to String!)"
Base.@propagate_inbounds substring(s::AbstractString,first::Integer,last::Integer) = substring(string(s[first:last]))

Base.@propagate_inbounds substring(s::AbstractString,first::Integer) = substring(s,first, lastind(s))

Base.@propagate_inbounds substring(offset::UInt32, size::UInt64, s::AbstractString) = SubString(offset,size,string(s))

Base.@propagate_inbounds substring(offset::UInt32, size::UInt64, s::String) = SubString(offset,size,s)

"universal fallback implementation in (offset,size)-notation"
substring(offset::UInt32, size::UInt64, s) = substring(offset, size, string(s))

Base.@propagate_inbounds function substring(offset::UInt32, size::UInt64, s::SubString{String})
    @boundscheck checksize(offset+size,usize(s)) # necessary to ensure substring is in bounds of s, not only of s.string
    @inbounds return SubString(offset+s.offset,size,s.string)
end

## enhancements of write and read - with prefix t to avoid type piracy

"""
optimized 'substring writing' of code units from a string buffer segment.

Not overloading Base.write because of inconsistency with generic multi-parameter-write
"""
function sswrite(io::IO, ofs::UInt32, size::UInt64, s::String)
    @boundscheck inbounds(s,ofs+size)
    GC.@preserve s unsafe_write(io, pointer(s)+ofs, size) # maybe fails on 32 bit julia implementations? size is not type UInt
end

"""
optimized reading of bytes (aka code units) into a string buffer
"""
function Base.read(io::IO, size::UInt64, ::Type{String})
    ss = _string_n(size)
    GC.@preserve ss unsafe_read(io, pointer(ss), UInt(size))
end


text:
module LaserLab
include("dffunctions.jl")
include("setup.jl")
include("obectives.jl")
include("gaussianLaser.jl")
include("tfunctions.jl")
include("util.jl")
include("histos.jl")
include("fits.jl")
include("images.jl")
end

max_stars_repo_name: jjgomezcadenas/laserlab

max_issues_repo_name: jjgomezcadenas/laserlab

max_forks_repo_name: jjgomezcadenas/laserlab


text:
# DESCRIPTION
# Computes the linearized dynamics of a nonlinear system (acfun), given the
# trim state and control vectors (x0 and u0).
#
# INPUT
# - acfun: name of system dynamics function
# - x0: trim state vector
# - u0: trim control vector
# - xdot_in: state derivative
# - t: time
#
# OUTPUT
# - F, G, M, A, B: coefficient matrices of linearized dynamics
#
# ------------------------------------------------------------------------------

function linearize(acfun,x0,u0,xdot_in,t)
    # init state-space matrices
    F=zeros(NSTATES,NSTATES)
    G=zeros(NSTATES,NCTRLS)
    M=Matrix{Float64}(I, NSTATES, NSTATES)
    A=zeros(NSTATES,NSTATES)
    B=zeros(NSTATES,NCTRLS)
    # F matrix
    x_p=zeros(NSTATES)
    xd_p=zeros(NSTATES)
    u_p=zeros(NCTRLS)
    for k=1:NSTATES
        #x_p=x0
        x_p[:]=x0[:]
        x_p[k]=x_p[k]+DELXLIN[k]
        xdot_p1=acfun(xdot_in,x_p,u0,t)
        #xdot_p1=H60!(xdot_in,x_p,u0,t)
        x_p[k]=x_p[k]-2*DELXLIN[k]
        #xdot_p2=eval(acfun,x_p,u0,xdot_in)
        xdot_p2=acfun(xdot_in,x_p,u0,t)
        F[:,k]=(xdot_p1-xdot_p2)/(2*DELXLIN[k])
    end
    # M matrix
    xd_p=zeros(NSTATES)
    for k=1:NSTATES
        #xd_p=xdot_in
        xd_p[:]=xdot_in[:]
        xd_p[k]=xd_p[k]+DELXLIN[k]
        xdot_p1=acfun(xd_p,x0,u0,t)
        xd_p[k]=xd_p[k]-2*DELXLIN[k]
        xdot_p2=acfun(xd_p,x0,u0,t)
        M[:,k]=M[:,k]-(xdot_p1-xdot_p2)/(2*DELXLIN[k])
    end
    # G matrix
    u_p=zeros(NCTRLS)
    for k=1:NCTRLS
        #u_p=u0
        u_p[:]=u0[:]
        u_p[k]=u_p[k]+DELCLIN[k]
        xdot_p1=acfun(xdot_in,x0,u_p,t)
        u_p[k]=u_p[k]-2*DELCLIN[k]
        xdot_p2=acfun(xdot_in,x0,u_p,t)
        G[:,k]=(xdot_p1-xdot_p2)/(2*DELCLIN[k])
    end
    A=M\F
    B=M\G
    return F, G, M, A, B
end

max_stars_repo_name: umbertosaetti/J-GenHel

max_issues_repo_name: umbertosaetti/J-GenHel

max_forks_repo_name: umbertosaetti/J-GenHel


text:
using Turing, Test

N = 5
setchunksize(4*N)
alg = HMC(0.2, 4, 1000)

@model vdemo(::Type{T}=Float64) where {T} = begin
  v = Vector{Matrix{T}}(undef, N)
  v ~ [Wishart(7, [1 0.5; 0.5 1])]
end

t_vec = @elapsed res_vec = sample(vdemo(), alg)

@model vdemo() = begin
  v = Vector{Matrix{Real}}(undef, N)
  for i = 1:N
    v[i] ~ Wishart(7, [1 0.5; 0.5 1])
  end
end

t_loop = @elapsed res = sample(vdemo(), alg)


text:
function commander(c::Client, m::Message, ::Val{:game_master})
    startswith(m.content, COMMAND_PREFIX * "gm") || return nothing

    regex = Regex(COMMAND_PREFIX * raw"gm( help| in| out)? *(.*)$")

    matches = match(regex, m.content)

    if matches === nothing || matches.captures[1] ∈ (" help", nothing)
        help_commander(c, m, :game_master)
    elseif matches.captures[1] == " in"
        @info "opt-in was required" m.content m.author.username m.author.discriminator
        opt_in(c, m, :game_master)
    elseif matches.captures[1] == " out"
        @info "opt-out was required" m.content m.author.username m.author.discriminator
        opt_out(c, m, :game_master)
    else
        reply(c, m, "Sorry, are you playing me? Please check out `gm help`")
    end
    return nothing
end

function help_commander(c::Client, m::Message, ::Val{:game_master})
    return reply(
        c,
        m,
        """
        How to play with the `gm` command:
        ```
        gm help
        gm in
        gm out
        ```
        The commands `in` and `out` are to opt-in and opt-out of the game. Playing data of an `out` player are kept until the end of the game. An `out` player can come `in` anytime to resume its participation!
        """,
    )
end


text:
function encode(s)
    encoded = ""
    if !isempty(s)
        char_diff = append!([s[i-1]==s[i] for i=2:length(s)], false)
        count = 1
        for i=1:length(char_diff)
            if char_diff[i]
                count += 1
            else
                encoded = encoded * (count>1 ? string(count) : "") * s[i]
                count = 1
            end
        end
    end
    return encoded
end

function decode(s)
    n_buffer = ""
    decoded = ""
    for c in s
        if isletter(c) || isspace(c)
            if isempty(n_buffer)
                decoded = decoded*c
            else
                decoded = decoded*c^parse(Int, n_buffer)
            end
            n_buffer = ""
        else
            n_buffer = n_buffer*c
        end
    end
    return decoded
end

encoded = encode("   hsqq qww  ")

max_stars_repo_name: AbstractGeek/exercism-expeditions

max_issues_repo_name: AbstractGeek/exercism-expeditions

max_forks_repo_name: AbstractGeek/exercism-expeditions


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

const Chars = Union{Char,Tuple{Vararg{Char}},AbstractVector{Char},Set{Char}}

"""
    search(string::AbstractString, chars::Chars, [start::Integer])

Search for the first occurrence of the given characters within the given string. The second
argument may be a single character, a vector or a set of characters, a string, or a regular
expression (though regular expressions are only allowed on contiguous strings, such as ASCII
or UTF-8 strings). The third argument optionally specifies a starting index. The return
value is a range of indexes where the matching sequence is found, such that `s[search(s,x)] == x`:

`search(string, "substring")` = `start:end` such that `string[start:end] == "substring"`, or
`0:-1` if unmatched.

`search(string, 'c')` = `index` such that `string[index] == 'c'`, or `0` if unmatched.

```jldoctest
julia> search("Hello to the world", "z")
0:-1

julia> search("JuliaLang","Julia")
1:5
```
"""
function search(s::AbstractString, c::Chars, i::Integer)
    if isempty(c)
        return 1 <= i <= nextind(s,endof(s)) ? i :
               throw(BoundsError(s, i))
    end
    if i < 1 || i > nextind(s,endof(s))
        throw(BoundsError(s, i))
    end
    while !done(s,i)
        d, j = next(s,i)
        if d in c
            return i
        end
        i = j
    end
    return 0
end
search(s::AbstractString, c::Chars) = search(s,c,start(s))

in(c::Char, s::AbstractString) = (search(s,c)!=0)

function _searchindex(s, t, i)
    if isempty(t)
        return 1 <= i <= nextind(s,endof(s)) ? i :
               throw(BoundsError(s, i))
    end
    t1, j2 = next(t,start(t))
    while true
        i = search(s,t1,i)
        if i == 0 return 0 end
        c, ii = next(s,i)
        j = j2; k = ii
        matched = true
        while !done(t,j)
            if done(s,k)
                matched = false
                break
            end
            c, k = next(s,k)
            d, j = next(t,j)
            if c != d
                matched = false
                break
            end
        end
        if matched
            return i
        end
        i = ii
    end
end

function _search_bloom_mask(c)
    UInt64(1) << (c & 63)
end

_nthbyte(s::String, i) = codeunit(s, i)
_nthbyte(a::ByteArray, i) = a[i]

function _searchindex(s::Union{String,ByteArray}, t::Union{String,ByteArray}, i)
    n = sizeof(t)
    m = sizeof(s)

    if n == 0
        return 1 <= i <= m+1 ? max(1, i) : 0
    elseif m == 0
        return 0
    elseif n == 1
        return search(s, _nthbyte(t,1), i)
    end

    w = m - n
    if w < 0 || i - 1 > w
        return 0
    end

    bloom_mask = UInt64(0)
    skip = n - 1
    tlast = _nthbyte(t,n)
    for j in 1:n
        bloom_mask |= _search_bloom_mask(_nthbyte(t,j))
        if _nthbyte(t,j) == tlast && j < n
            skip = n - j - 1
        end
    end

    i -= 1
    while i <= w
        if _nthbyte(s,i+n) == tlast
            # check candidate
            j = 0
            while j < n - 1
                if _nthbyte(s,i+j+1) != _nthbyte(t,j+1)
                    break
                end
                j += 1
            end

            # match found
            if j == n - 1
                return i+1
            end

            # no match, try to rule out the next character
            if i < w && bloom_mask & _search_bloom_mask(_nthbyte(s,i+n+1)) == 0
                i += n
            else
                i += skip
            end
        elseif i < w
            if bloom_mask & _search_bloom_mask(_nthbyte(s,i+n+1)) == 0
                i += n
            end
        end
        i += 1
    end

    0
end

searchindex(s::ByteArray, t::ByteArray, i) = _searchindex(s,t,i)

"""
    searchindex(s::AbstractString, substring, [start::Integer])

Similar to [`search`](@ref), but return only the start index at which
the substring is found, or `0` if it is not.

```jldoctest
julia> searchindex("Hello to the world", "z")
0

julia> searchindex("JuliaLang","Julia")
1

julia> searchindex("JuliaLang","Lang")
6
```
"""
searchindex(s::AbstractString, t::AbstractString, i::Integer) = _searchindex(s,t,i)
searchindex(s::AbstractString, t::AbstractString) = searchindex(s,t,start(s))
searchindex(s::AbstractString, c::Char, i::Integer) = _searchindex(s,c,i)
searchindex(s::AbstractString, c::Char) = searchindex(s,c,start(s))

function searchindex(s::String, t::String, i::Integer=1)
    # Check for fast case of a single byte
    # (for multi-byte UTF-8 sequences, use searchindex on byte arrays instead)
    if endof(t) == 1
        search(s, t[1], i)
    else
        _searchindex(s, t, i)
    end
end

function _search(s, t, i::Integer)
    idx = searchindex(s,t,i)
    if isempty(t)
        idx:idx-1
    else
        idx:(idx > 0 ? idx + endof(t) - 1 : -1)
    end
end

search(s::AbstractString, t::AbstractString, i::Integer=start(s)) = _search(s, t, i)
search(s::ByteArray, t::ByteArray, i::Integer=start(s)) = _search(s, t, i)

function rsearch(s::AbstractString, c::Chars)
    j = search(RevString(s), c)
    j == 0 && return 0
    endof(s)-j+1
end

"""
    rsearch(s::AbstractString, chars::Chars, [start::Integer])

Similar to [`search`](@ref), but returning the last occurrence of the given characters within the
given string, searching in reverse from `start`.

```jldoctest
julia> rsearch("aaabbb","b")
6:6
```
"""
function rsearch(s::AbstractString, c::Chars, i::Integer)
    e = endof(s)
    j = search(RevString(s), c, e-i+1)
    j == 0 && return 0
    e-j+1
end

function _rsearchindex(s, t, i)
    if isempty(t)
        return 1 <= i <= nextind(s,endof(s)) ? i :
               throw(BoundsError(s, i))
    end
    t = RevString(t)
    rs = RevString(s)
    l = endof(s)
    t1, j2 = next(t,start(t))
    while true
        i = rsearch(s,t1,i)
        if i == 0 return 0 end
        c, ii = next(rs,l-i+1)
        j = j2; k = ii
        matched = true
        while !done(t,j)
            if done(rs,k)
                matched = false
                break
            end
            c, k = next(rs,k)
            d, j = next(t,j)
            if c != d
                matched = false
                break
            end
        end
        if matched
            return nextind(s,l-k+1)
        end
        i = l-ii+1
    end
end

function _rsearchindex(s::Union{String,ByteArray}, t::Union{String,ByteArray}, k)
    n = sizeof(t)
    m = sizeof(s)

    if n == 0
        return 0 <= k <= m ? max(k, 1) : 0
    elseif m == 0
        return 0
    elseif n == 1
        return rsearch(s, _nthbyte(t,1), k)
    end

    w = m - n
    if w < 0 || k <= 0
        return 0
    end

    bloom_mask = UInt64(0)
    skip = n - 1
    tfirst = _nthbyte(t,1)
    for j in n:-1:1
        bloom_mask |= _search_bloom_mask(_nthbyte(t,j))
        if _nthbyte(t,j) == tfirst && j > 1
            skip = j - 2
        end
    end

    i = min(k - n + 1, w + 1)
    while i > 0
        if _nthbyte(s,i) == tfirst
            # check candidate
            j = 1
            while j < n
                if _nthbyte(s,i+j) != _nthbyte(t,j+1)
                    break
                end
                j += 1
            end

            # match found
            if j == n
                return i
            end

            # no match, try to rule out the next character
            if i > 1 && bloom_mask & _search_bloom_mask(_nthbyte(s,i-1)) == 0
                i -= n
            else
                i -= skip
            end
        elseif i > 1
            if bloom_mask & _search_bloom_mask(_nthbyte(s,i-1)) == 0
                i -= n
            end
        end
        i -= 1
    end

    0
end

rsearchindex(s::ByteArray, t::ByteArray, i::Integer) = _rsearchindex(s,t,i)

"""
    rsearchindex(s::AbstractString, substring, [start::Integer])

Similar to [`rsearch`](@ref), but return only the start index at which the substring is found, or `0` if it is not.

```jldoctest
julia> rsearchindex("aaabbb","b")
6

julia> rsearchindex("aaabbb","a")
3
```
"""
rsearchindex(s::AbstractString, t::AbstractString, i::Integer) = _rsearchindex(s,t,i)
rsearchindex(s::AbstractString, t::AbstractString) = (isempty(s) && isempty(t)) ? 1 : rsearchindex(s,t,endof(s))

function rsearchindex(s::String, t::String)
    # Check for fast case of a single byte
    # (for multi-byte UTF-8 sequences, use rsearchindex instead)
    if endof(t) == 1
        rsearch(s, t[1])
    else
        _rsearchindex(s, t, sizeof(s))
    end
end

function rsearchindex(s::String, t::String, i::Integer)
    # Check for fast case of a single byte
    # (for multi-byte UTF-8 sequences, use rsearchindex instead)
    if endof(t) == 1
        rsearch(s, t[1], i)
    elseif endof(t) != 0
        _rsearchindex(s, t, nextind(s, i)-1)
    elseif i > sizeof(s)
        return 0
    elseif i == 0
        return 1
    else
        return i
    end
end

function _rsearch(s, t, i::Integer)
    idx = rsearchindex(s,t,i)
    if isempty(t)
        idx:idx-1
    else
        idx:(idx > 0 ? idx + endof(t) - 1 : -1)
    end
end

rsearch(s::AbstractString, t::AbstractString, i::Integer=endof(s)) = _rsearch(s, t, i)
rsearch(s::ByteArray, t::ByteArray, i::Integer=endof(s)) = _rsearch(s, t, i)

"""
    contains(haystack::AbstractString, needle::AbstractString)

Determine whether the second argument is a substring of the first.

```jldoctest
julia> contains("JuliaLang is pretty cool!", "Julia")
true
```
"""
contains(haystack::AbstractString, needle::AbstractString) = searchindex(haystack,needle)!=0

in(::AbstractString, ::AbstractString) = error("use contains(x,y) for string containment")