text:
import Stripeline
const Sl = Stripeline

db = Sl.InstrumentDB(joinpath(dirname(pathof(Stripeline)), "..", "instrumentdb"))

@test length(db.focalplane) ≥ 49
@test length(db.detectors) ≥ 55
@test "I0" in keys(db.focalplane)

# Now test that everything works without specifying the paths

defaultdb = Sl.InstrumentDB()

@test sort(collect(keys(defaultdb.focalplane))) == sort(collect(keys(db.focalplane)))
@test sort(collect(keys(defaultdb.detectors))) == sort(collect(keys(db.detectors)))

# Check the high-level API

@test Sl.detector(defaultdb, "I0") != nothing
@test Sl.detector(defaultdb, 2) != nothing

@test Sl.spectrum(defaultdb, 2) != nothing
@test Sl.spectrum(defaultdb, "I0") != nothing

@test Sl.tnoise(defaultdb, 2) != nothing
@test Sl.tnoise(defaultdb, "I0") != nothing

@test Sl.bandpass(defaultdb, 2) != nothing
@test Sl.bandpass(defaultdb, "I0") != nothing

x, y = Sl.bandshape(defaultdb, "I0")
@test length(x) > 0
@test length(x) == length(y)

f1_i, f1_q, f1_u = Sl.fknee_hz(defaultdb, 10)
f2_i, f2_q, f2_u = Sl.fknee_hz(defaultdb, 10, tsys_k = 1000.0)

@test f2_q < f1_q
@test f2_u < f1_u

# Just a few basic checks, to verify that the function does not crash
(sensitivity, num) = Sl.sensitivity_tant(defaultdb, 10.0)
@test sensitivity > 0.0
@test num == 49

(sensitivity, num) = Sl.sensitivity_tant(defaultdb, 10.0, modules = Set([0, 1]))
@test sensitivity > 0.0
@test num == 14

@test Sl.t_to_trj(2.72548, 43e9) ≈ 1.82263012068
@test Sl.trj_to_t(Sl.t_to_trj(2.72548, 43e9), 43e9) ≈ 2.72548


text:
# ---
# cover: nonexistence.jpg
# ---

# ![valid_cover](../logo.png)

# front matter overrides image link

ext: jl


text:
function pairwiseAlloc(Z, W, i)
  J, K = size(Z)
  A = zeros(J, J)
  for r in 1:J
    for c in 1:J
      for k in 1:K
        if Z[r, k] == 1 && Z[c, k] == 1
          A[r, c] += W[i, k]
          A[c, r] += W[i, k]
        end
      end
    end
  end

  return A
end # pairwiseAlloc

function estimate_ZWi_index(Zs, Ws, i)
  As = [pairwiseAlloc(z, w, i) for (z, w) in zip(Zs, Ws)]

  Amean = mean(As)
  mse = [ mean((A - Amean) .^ 2) for A in As]

  return argmin(mse)
end


text:
export BrickletRS232Read
struct BrickletRS232Read
    message::Vector{Char}
    length::Integer
end

export BrickletRS232Configuration
struct BrickletRS232Configuration
    baudrate::Integer
    parity::Integer
    stopbits::Integer
    wordlength::Integer
    hardware_flowcontrol::Integer
    software_flowcontrol::Integer
end

export BrickletRS232ReadFrame
struct BrickletRS232ReadFrame
    message::Vector{Char}
    length::Integer
end

export BrickletRS232Identity
struct BrickletRS232Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletRS232
"""
Communicates with RS232 devices
"""
mutable struct BrickletRS232 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletRS232(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_rs232")
        deviceInternal = package.BrickletRS232(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export write
"""
    $(SIGNATURES)

Writes a string of up to 60 characters to the RS232 interface. The string
can be binary data, ASCII or similar is not necessary.

The length of the string has to be given as an additional parameter.

The return value is the number of bytes that could be written.

See :func:`Set Configuration` for configuration possibilities
regarding baudrate, parity and so on.
"""
function write(device::BrickletRS232, message, length)
    return device.deviceInternal.write(message, length)
end

export read
"""
    $(SIGNATURES)

Returns the currently buffered message. The maximum length
of message is 60. If the returned length is 0, no new data was available.

Instead of polling with this function, you can also use
callbacks. See :func:`Enable Read Callback` and :cb:`Read` callback.
"""
function read(device::BrickletRS232)
    return device.deviceInternal.read()
end

export enable_read_callback
"""
    $(SIGNATURES)

Enables the :cb:`Read` callback. This will disable the :cb:`Frame Readable` callback.

By default the callback is disabled.
"""
function enable_read_callback(device::BrickletRS232)
    device.deviceInternal.enable_read_callback()
end

export disable_read_callback
"""
    $(SIGNATURES)

Disables the :cb:`Read` callback.

By default the callback is disabled.
"""
function disable_read_callback(device::BrickletRS232)
    device.deviceInternal.disable_read_callback()
end

export is_read_callback_enabled
"""
    $(SIGNATURES)

Returns *true* if the :cb:`Read` callback is enabled,
*false* otherwise.
"""
function is_read_callback_enabled(device::BrickletRS232)
    return device.deviceInternal.is_read_callback_enabled()
end

export set_configuration
"""
    $(SIGNATURES)

Sets the configuration for the RS232 communication.

Hard-/Software flow control can either be on or off but not both simultaneously on.
"""
function set_configuration(device::BrickletRS232, baudrate, parity, stopbits, wordlength, hardware_flowcontrol, software_flowcontrol)
    device.deviceInternal.set_configuration(baudrate, parity, stopbits, wordlength, hardware_flowcontrol, software_flowcontrol)
end

export get_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Configuration`.
"""
function get_configuration(device::BrickletRS232)
    return device.deviceInternal.get_configuration()
end

export set_break_condition
"""
    $(SIGNATURES)

Sets a break condition (the TX output is forced to a logic 0 state).
The parameter sets the hold-time of the break condition.

.. versionadded:: 2.0.2\$nbsp;(Plugin)
"""
function set_break_condition(device::BrickletRS232, break_time)
    device.deviceInternal.set_break_condition(break_time)
end

export set_frame_readable_callback_configuration
"""
    $(SIGNATURES)

Configures the :cb:`Frame Readable` callback. The frame size is the number of bytes, that have to be readable to trigger the callback.
A frame size of 0 disables the callback. A frame size greater than 0 enables the callback and disables the :cb:`Read` callback.

By default the callback is disabled.

.. versionadded:: 2.0.4\$nbsp;(Plugin)
"""
function set_frame_readable_callback_configuration(device::BrickletRS232, frame_size)
    device.deviceInternal.set_frame_readable_callback_configuration(frame_size)
end

export get_frame_readable_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Frame Readable Callback Configuration`.

.. versionadded:: 2.0.4\$nbsp;(Plugin)
"""
function get_frame_readable_callback_configuration(device::BrickletRS232)
    return device.deviceInternal.get_frame_readable_callback_configuration()
end

export read_frame
"""
    $(SIGNATURES)

Returns up to one frame of bytes from the read buffer.
The frame size is configured with :func:`Set Frame Readable Callback Configuration`.
If the returned length is 0, no new data was available.

.. versionadded:: 2.0.4\$nbsp;(Plugin)
"""
function read_frame(device::BrickletRS232)
    return device.deviceInternal.read_frame()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletRS232)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletRS232, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
# This script shows how to compute risk-adjusted linearizations using sparse arrays
# for caches and sparse differentiation for the Jacobians of μ, ξ, and 𝒱 .
using RiskAdjustedLinearizations, LinearAlgebra, SparseArrays
using BenchmarkTools, Test, SparseDiffTools

# Settings
define_functions = true
time_methods     = true
algorithm        = :relaxation # Note that while both methods work with sparse array caches, only relaxation works for sparse differentiation.
N_approx         = 10          # Number of periods ahead used for forward-difference equations

if define_functions
    include(joinpath(dirname(@__FILE__), "..", "nk_with_capital", "nk_with_capital.jl"))
end

# Set up
m_nk = NKCapital(; N_approx = N_approx) # create parameters

# Sparse arrays for caches
m = nk_capital(m_nk; sparse_arrays = true)

## The keyword sparse_arrays tells the nk_capital function
## to make Γ₅, Γ₆ sparse arrays and to add the keyword arguments
## `Λ_cache_init = dims -> spzeros(dims...)`
## `Σ_cache_init = dims -> spzeros(dims...)`
## when calling the constructor, i.e.
## RiskAdjustedLinearization(...; Λ_cache_init = ...)
solve!(m; algorithm = algorithm, verbose = :none)

# Risk-adjusted linearization with sparse differentiation
# for the Jacobians of μ, ξ, and 𝒱 in addition to sparse caches.

## The first approach directly uses
## the constructor of a RiskAdjustedLinearization
## to determine the sparsity pattern.
## Since the initial guess for Ψ is a matrix of zeros,
## we will only use sparse differentiation for μ and ξ at first.
## Within nk_capital, the keyword `sparse_jacobian = [:μ, :ξ]` is passed
## to the constructor, i.e.
## RiskAdjustedLinearization(...; sparse_jacobian = ...)
m_sparsejac = nk_capital(m_nk; sparse_arrays = true, sparse_jacobian = [:μ, :ξ])

## Check the caches for the Jacobians of μ and ξ are actually sparse
@test issparse(m_sparsejac[:Γ₁])
@test issparse(m_sparsejac[:Γ₂])
@test issparse(m_sparsejac[:Γ₃])
@test issparse(m_sparsejac[:Γ₄])

## Now solve the model! Note that sparse differentiation
## can be fragile sometimes and result in NaNs or undefined
## numbers appearing during calls to `nlsolve`. Re-running
## `solve!` repeatedly or reconstructing `m_sparsejac` again
## will usually lead to a successful run.
solve!(m_sparsejac; algorithm = algorithm)
@test norm(steady_state_errors(m_sparsejac), Inf) < 1e-8

## The second approach calls `update_sparsity_pattern!`
## on an existing `RiskAdjustedLinearization`

### Create dictionaries for specifying the sparsity pattern
### Here, we will tell `m_sparsejac`
### to now use sparse differentiation for J𝒱
### and to use a new sparsity pattern for the
### Jacobians of μ and ξ.
sparsity = Dict()
colorvec = Dict()
sparsity[:J𝒱] = sparse(m[:JV])
@test isempty(sparsity[:J𝒱].nzval) # However, in the case of the current NK model, the entropy actually does not depend on z.
colorvec[:J𝒱] = ones(Int64, size(sparsity[:J𝒱], 2)) # Instead, we just pass in a coloring vector of this form
jac_to_sym = (μz = :Γ₁, μy = :Γ₂, ξz = :Γ₃, ξy = :Γ₄)
for k in [:μz, :μy, :ξz, :ξy]
    sparsity[k] = m_sparsejac[jac_to_sym[k]]
    colorvec[k] = matrix_colors(sparsity[k])
end

@test !issparse(m_sparsejac[:JV]) # JV is not sparse currently

update_sparsity_pattern!(m_sparsejac, [:μ, :ξ]; sparsity = sparsity) # Don't have to provide the matrix coloring vector (note 𝒱 is not included
                                                                      # b/c calling matrix_color on its sparsity pattern will error)
update_sparsity_pattern!(m_sparsejac, [:μ, :ξ, :𝒱]; sparsity = sparsity, # But if the coloring vector already exists,
                         colorvec = colorvec)                            # then you may as well pass that information, too.

@test issparse(m_sparsejac[:JV]) # Now JV is sparse

### If you ever need to look at the sparsity pattern or coloring vector,
### you can call `linearized_system(m_sparsejac).sparse_jac_caches` or
### `m_sparsejac.linearization.sparse_jac_caches`, which is a
### `NamedTuple` whose values include the sparsity pattern and coloring vector,
### as well as a cache used by forwarddiff_color_jacobian!
for k in [:μz, :μy, :ξz, :ξy, :J𝒱]
    @test issparse(m_sparsejac.linearization.sparse_jac_caches[k][:sparsity])
end

### Now solve again with sparse Jacobian of 𝒱, too!
solve!(m_sparsejac, m_sparsejac.z .* 1.01, m_sparsejac.y .* 1.01, m_sparsejac.Ψ .* 1.01; algorithm = algorithm)
@test norm(steady_state_errors(m_sparsejac), Inf) < 1e-8

if time_methods
    m_dense = nk_capital(m_nk)
    zinit = deepcopy(m_dense.z)
    yinit = deepcopy(m_dense.y)
    Ψinit = deepcopy(m_dense.Ψ)

    solve!(m; algorithm = :deterministic, verbose = :none)
    zdet = deepcopy(m.z)
    ydet = deepcopy(m.y)
    Ψdet = deepcopy(m.Ψ)

    println("Deterministic steady state with dense caches for Γ₅, Γ₆, Λ, and Σ")
    @btime begin
        update!(m_dense, zinit, yinit, Ψinit)
        solve!(m_dense; algorithm = :deterministic, verbose = :none)
    end

    println("Deterministic steady state with sparse caches for Γ₅, Γ₆, Λ, and Σ")
    @btime begin
        update!(m, zinit, yinit, Ψinit)
        solve!(m; algorithm = :deterministic, verbose = :none)
    end

    println("Relaxation with dense caches for Γ₅, Γ₆, Λ, and Σ")
    @btime begin
        update!(m_dense, zdet, ydet, Ψdet)
        solve!(m_dense; algorithm = :relaxation, verbose = :none)
    end

    println("Relaxation with sparse caches for Γ₅, Γ₆, Λ, and Σ")
    @btime begin
        update!(m, zdet, ydet, Ψdet)
        solve!(m; algorithm = :relaxation, verbose = :none)
    end

    println("Homotopy with dense caches for Γ₅, Γ₆, Λ, and Σ")
    @btime begin
        update!(m_dense, zdet, ydet, Ψdet)
        solve!(m_dense; algorithm = :homotopy, verbose = :none)
    end

    println("Homotopy with sparse caches for Γ₅, Γ₆, Λ, and Σ")
    @btime begin
        update!(m, zdet, ydet, Ψdet)
        solve!(m; algorithm = :homotopy, verbose = :none)
    end

    println("Relaxation with sparse caches and differentiation")
    @btime begin
        update!(m_sparsejac, zdet, ydet, Ψdet)
        solve!(m; algorithm = :relaxation, verbose = :none)
    end
end


text:
# PROX_SEPARABLE   Evaluate the prox operator of a fully separable function.
#
# Arguments:
#
#  v is the point at which to evaluate the operator.
#  fp is a subgradient oracle for the function.
#  lambda (optional) is the proximal parameter; defaults to 1.
#  l (optional) is a lower bound for x; defaults to -Inf.
#  u (optional) is an upper bound for x; defaults to Inf.
#  x0 (optional) is a value at which to warm start the algorithm.
#  tol (optional) is a stopping tolerance.
#  MAX_ITER (optional) is the maximum number of iterations.
#
# Examples:
#
#  v = randn(n,1);
#  x = prox_separable(v, 1, @(w) sign(w));
#  [x iter] = prox_separable(v, 1, @(w) sign(w));
#
# This function can be called in vectorized form if fp is vectorized,
# i.e., if fp works elementwise, then v, l, u, and x0 can be vectors.
function prox_separable(
    v,
    lambda::Real = NaN,
    fp::Function,
    l,
    u,
    x0,
    tol::Real = NaN,
    maxiter::Integer = 0
)
    n = length(v)

    if isnan(lambda)
        lambda = 1.0
    end
    rho = 1 / lambda

    if any(isnan(l)) || isempty(l)
        l = -inf(n, 1)
    end
    if any(isnan(u)) || isempty(u)
        u = inf(n, 1)
    end
    if any(isnan(x0)) || isempty(x0)
        x0 = zeros(n, 1)
    end
    if isnan(tol)
        tol = 1e-8
    end
    if maxiter == 0
        maxiter = 500
    end

    iter = 0
    x = max(l, min(x0, u))

    while any(u .- l .> tol) && iter < MAX_ITER
        g = fp(x) .+ rho .* (x .- v)

        idx = g .> 0
        l[idx] = max(l[idx], x[idx] .- g[idx] ./ rho)
        u[idx] = x[idx]

        idx = !idx
        u[idx] = min(u[idx], x[idx] .- g[idx] ./ rho)
        l[idx] = x[idx]

        x = (l .+ u) ./ 2
        iter += 1
    end

    if any(u .- l .> tol)
        @printf(
            STDERR,
            "Warning: %d entries did not converge; max interval size = %f.",
            sum(u .- l .> tol),
            max(u .- l)
        )
    end

    return x, iter
end


text:
#===================================
Methods using Interval Optimisation
===================================#
using .IntervalOptimisation

function _enclose(mse::MooreSkelboeEnclosure, f::Function, dom::Interval_or_IntervalBox;
                  kwargs...)

    global_min, _ = minimise(X->f(X...), dom, structure=mse.structure, tol=mse.tol)
    global_max, _ = maximise(X->f(X...), dom, structure=mse.structure, tol=mse.tol)

    return Interval(inf(global_min), sup(global_max))
end


text:
# Create associative arrays from sparse matrices
using JLD

# Iterate through files
(@isdefined Nfile) || (Nfile = 8)

for i = 1:Nfile
    assoctime = @elapsed begin
        # Create filename
        fname = joinpath(Base.source_dir(),"data", string(i))
        println(fname)
        println("On file: "*string(i))

        # Open files, read data, and close files
        fidRow = open(fname * "r.txt","r")
        fidCol = open(fname * "c.txt","r")
        fidVal = open(fname * "v.txt","r")
        rowStr = readlines(fidRow)
        colStr = readlines(fidCol)
        valStr = readlines(fidVal)
        close(fidRow)
        close(fidCol)
        close(fidVal)

        # Construct adjacency associative array and sum duplicates
        A = Assoc(rowStr,colStr,1)
        save(fname*".A.jld","A",A)
        #saveassoc(fname*".A.jld",A)

        if i == 1; global testA = A; end

        # Get versions with duplicates summed together
        rowStr,colStr,valStr = find(A)
        Nedge = length(rowStr)
        edgeBit = ceil(Int,log10(i.*Nedge))

        # Create identifier for each edge
        edgeStr = reverse.(lpad.(collect(((i-1).*Nedge) .+ (1:Nedge)),edgeBit,"0"))

        # Format out and in edge strings
        outStr = CatStr(fill("Out",length(rowStr)),"|",rowStr)
        inStr = CatStr(fill("In",length(colStr)),"|",colStr)

        # Create directed incidence Assoc
        E = Assoc([edgeStr; edgeStr],[outStr; inStr],[valStr; valStr])
        save(fname*".E.jld","E",E)
        #saveassoc(fname*".E.jld",E)

        if i == 1; global testE = E; end

    end

    println("Time: "*string(assoctime)*", Edges/sec: "*string(length(rowStr)./assoctime))
end


text:
agent = RL.Agent(
    policy = RL.QBasedPolicy(
        learner = RL.DQNLearner(
            approximator = RL.NeuralNetworkApproximator(
                model = SeaPearl.VariableOutputCPNN(
                    graphChain = Flux.Chain(
                        SeaPearl.EdgeFtLayer(; v_dim=numInFeatures => 32, e_dim= 1 => 4),
                        SeaPearl.EdgeFtLayer(; v_dim=32 => 32, e_dim= 4 => 4),
                        SeaPearl.EdgeFtLayer(; v_dim=32 => 32, e_dim= 4 => 4),
                        SeaPearl.EdgeFtLayer(; v_dim=32 => 32, e_dim= 4 => 4),
                    ),
                    nodeChain = Flux.Chain(
                        Flux.Dense(32, 32, relu),
                        Flux.Dense(32, 32, relu),
                    ),
                    outputChain = Flux.Dense(64, 1),
                ),
                optimizer = ADAM(0.0001f0)
            ),
            target_approximator = RL.NeuralNetworkApproximator(
                model = SeaPearl.VariableOutputCPNN(
                    graphChain = Flux.Chain(
                        SeaPearl.EdgeFtLayer(; v_dim=numInFeatures => 32, e_dim= 1 => 4),
                        SeaPearl.EdgeFtLayer(; v_dim=32 => 32, e_dim= 4 => 4),
                        SeaPearl.EdgeFtLayer(; v_dim=32 => 32, e_dim= 4 => 4),
                        SeaPearl.EdgeFtLayer(; v_dim=32 => 32, e_dim= 4 => 4),
                    ),
                    nodeChain = Flux.Chain(
                        Flux.Dense(32, 32, relu),
                        Flux.Dense(32, 32, relu),
                    ),
                    outputChain = Flux.Dense(64, 1),
                ),
                optimizer = ADAM(0.0001f0)
            ),
            loss_func = Flux.Losses.huber_loss,
            stack_size = nothing,
            γ = 0.99f0,
            batch_size = 1,#32,
            update_horizon = 1,
            min_replay_history = 1,
            update_freq = 1,
            target_update_freq = 100,
        ),
        explorer = RL.EpsilonGreedyExplorer(
            ϵ_stable = 0.01,
            kind = :exp,
            ϵ_init = 1.0,
            warmup_steps = 0,
            decay_steps = 2000,
            step = 1,
            is_break_tie = false,
            #is_training = true,
            rng = MersenneTwister(33)
        )
    ),
    trajectory = RL.CircularArraySARTTrajectory(
        capacity = 4000,
        state = SeaPearl.TsptwTrajectoryState[] => ()
    )
)

max_stars_repo_name: corail-research/SeaPearlZoo

max_issues_repo_name: corail-research/SeaPearlZoo

max_forks_repo_name: corail-research/SeaPearlZoo


text:
using HTTP: HTTP, request, URI
using JSON3: JSON3
using LibPQ: LibPQ, Connection, execute, load!
using DataFrames: DataFrame

conn = Connection("dbname = sdad");
page = 7139
while true
    global page
    response = request("POST",
                       URI(scheme = "https",
                           host = "api.usaspending.gov",
                           path = "/api/v2/recipient/duns"),
                       ["Content-Type" => "application/json"],
                       JSON3.write(Dict("sort" => "name", "limit" => 100, "page" => page)))
    json = JSON3.read(response.body)
    data = DataFrame(json.results)
    data[!,:duns] = something.(data[!,:duns], missing)
    execute(conn, "BEGIN;")
    load!(data,
          conn,
          string("INSERT INTO us_spending.recipients VALUES(",
                 join(("\$$i" for i in 1:size(data, 2)), ','),
                 ") ON CONFLICT DO NOTHING;"))
    execute(conn, "COMMIT;")
    json.page_metadata.hasNext || break
    page = json.page_metadata.next
    println(page)
    # sleep(0.5)
end

execute(conn, "CREATE TABLE ")


text:
module Frame
import Base: ==, AbstractIOBuffer

@enum FRAME_TYPES DATA=0x0 HEADERS=0x1 PRIORITY=0x2 RST_STREAM=0x3 SETTINGS=0x4 PUSH_PROMISE=0x5 PING=0x6 GOAWAY=0x7 WINDOW_UPDATE=0x8 CONTINUATION=0x9

struct FrameHeader
    length::UInt32
    typ::FRAME_TYPES
    flags::UInt8
    stream_identifier::UInt32
end

function decode_header(buf)
    length_arr = read(buf, 3)
    length = UInt32(length_arr[1]) << 16 + UInt32(length_arr[2]) << 8 + UInt32(length_arr[3])

    typ = FRAME_TYPES(read(buf, 1)[1])
    flags = read(buf, 1)[1]
    stream_identifier_arr = read(buf, 4)
    stream_identifier = UInt32(stream_identifier_arr[1]) << 24 + UInt32(stream_identifier_arr[2]) << 16 +
        UInt32(stream_identifier_arr[3]) << 8 + UInt32(stream_identifier_arr[4])

    @assert stream_identifier & 0x8000000000 == 0

    return FrameHeader(length, typ, flags, stream_identifier)
end

function encode_header(header::FrameHeader)
    buf = IOBuffer()

    write(buf, UInt8(header.length >> 16), UInt8((header.length >> 8) & 0x000000ff), UInt8(header.length & 0x000000ff))
    write(buf, UInt8(header.typ))
    write(buf, header.flags)

    @assert header.stream_identifier & 0x8000000000 == 0

    write(buf, UInt8(header.stream_identifier >> 24), UInt8((header.stream_identifier >> 16) & 0x000000ff),
          UInt8((header.stream_identifier >> 8) & 0x000000ff), UInt8(header.stream_identifier & 0x000000ff))

    return take!(buf)
end

struct UnimplementedError <: Exception end

include("Frame/utils.jl")
include("Frame/data.jl")
include("Frame/headers.jl")
include("Frame/priority.jl")
include("Frame/rst_stream.jl")
include("Frame/settings.jl")
include("Frame/push_promise.jl")
include("Frame/ping.jl")
include("Frame/goaway.jl")
include("Frame/window_update.jl")
include("Frame/continuation.jl")

function decode(buf)
    header = decode_header(buf)
    payload = read(buf, header.length)
    @assert length(payload) == header.length

    if header.typ == DATA
        return decode_data(header, payload)
    elseif header.typ == HEADERS
        return decode_headers(header, payload)
    elseif header.typ == PRIORITY
        return decode_priority(header, payload)
    elseif header.typ == RST_STREAM
        return decode_rst_stream(header, payload)
    elseif header.typ == SETTINGS
        return decode_settings(header, payload)
    elseif header.typ == PUSH_PROMISE
        return decode_push_promise(header, payload)
    elseif header.typ == PING
        return decode_ping(header, payload)
    elseif header.typ == GOAWAY
        return decode_goaway(header, payload)
    elseif header.typ == WINDOW_UPDATE
        return decode_window_update(header, payload)
    elseif header.typ == CONTINUATION
        return decode_continuation(header, payload)
    else
        throw(ParseError())
    end
end

encode(frame::DataFrame) = encode_data(frame)
encode(frame::HeadersFrame) = encode_headers(frame)
encode(frame::PriorityFrame) = encode_priority(frame)
encode(frame::RstStreamFrame) = encode_rst_stream(frame)
encode(frame::SettingsFrame) = encode_settings(frame)
encode(frame::PushPromiseFrame) = encode_push_promise(frame)
encode(frame::PingFrame) = encode_ping(frame)
encode(frame::GoawayFrame) = encode_goaway(frame)
encode(frame::WindowUpdateFrame) = encode_window_update(frame)
encode(frame::ContinuationFrame) = encode_continuation(frame)

export encode, decode, DataFrame, HeadersFrame, PriorityFrame, RstStreamFrame, SettingsFrame, PushPromiseFrame, PingFrame, GoawayFrame, WindowUpdateFrame, ContinuationFrame

end


text:
"""
    is_in_cell(p::Coord, cell::Cell)

Determines if a point (such as a Ray origin) is inside a given cell
"""
function is_in_cell(p::Coord, cell::Cell)
    result = navigate_tree(p, cell.regions, cell.definition)
    return result
end

function navigate_tree(p::Coord, r::Array{Region}, ex::Expr)
    global _p = Coord(p.x, p.y, p.z)

	# Check if Complement
	if ex.args[1] == :~
		if typeof(ex.args[2]) == typeof(1)
			return ~ r[ex.args[2]]
		else
			return ~ navigate_tree(p, r, ex.args[2])
		end
	end

	if typeof(ex.args[2]) == typeof(1)
		# Case 1 - Both operands are leaves
		if typeof(ex.args[3]) == typeof(1)
			if ex.args[1] == :^
            	return r[ex.args[2]] ^ r[ex.args[3]]
			end
			if ex.args[1] == :|
            	return r[ex.args[2]] | r[ex.args[3]]
			end
		end
		# Case 2 - Left operand is leaf, right is not
		if typeof(ex.args[3]) != typeof(1)
			if ex.args[1] == :^
            	return r[ex.args[2]] ^ navigate_tree(p, r, ex.args[3])
			end
			if ex.args[1] == :|
            	return r[ex.args[2]] | navigate_tree(p, r, ex.args[3])
			end
		end
	end

	if typeof(ex.args[2]) != typeof(1)
		# Case 3 - left operand is not leaf, but right is
		if typeof(ex.args[3]) == typeof(1)
			if ex.args[1] == :^
            	return navigate_tree(p, r, ex.args[2]) ^ r[ex.args[3]]
			end
			if ex.args[1] == :|
            	return navigate_tree(p, r, ex.args[2]) | r[ex.args[3]]
			end
		end
		# Case 4 - Neither operand is a leaf
		if typeof(ex.args[3]) != typeof(1)
			if ex.args[1] == :^
            	return navigate_tree(p, r, ex.args[2]) ^ navigate_tree(p, r, ex.args[3])
			end
			if ex.args[1] == :|
            	return navigate_tree(p, r, ex.args[2]) | navigate_tree(p, r, ex.args[3])
			end
		end
	end
end

"""
    find_cell_id(p::Coord, geometry::Geometry)

Finds the cell id that a point resides within
"""
function find_cell_id(p::Coord, geometry::Geometry)
    for i = 1:length(geometry.cells)
        if is_in_cell(p, geometry.cells[i]) == true
            return i
        end
    end
    return -1
end


text:
mutable struct DoorKeyDirected{T, R} <: AbstractGridWorld
    world::GridWorldBase{Tuple{Agent, Wall, Goal, Door, Key}}
    agent_pos::CartesianIndex{2}
    agent_dir::AbstractDirection
    reward::T
    rng::R
    terminal_reward::T
    goal_pos::CartesianIndex{2}
    door_pos::CartesianIndex{2}
    key_pos::CartesianIndex{2}
    partition::CartesianIndices{2,Tuple{UnitRange{Int},UnitRange{Int}}}
    has_key::Bool
    done::Bool
end

@generate_getters(DoorKeyDirected)
@generate_setters(DoorKeyDirected)

mutable struct DoorKeyUndirected{T, R} <: AbstractGridWorld
    world::GridWorldBase{Tuple{Agent, Wall, Goal, Door, Key}}
    agent_pos::CartesianIndex{2}
    reward::T
    rng::R
    terminal_reward::T
    goal_pos::CartesianIndex{2}
    door_pos::CartesianIndex{2}
    key_pos::CartesianIndex{2}
    partition::CartesianIndices{2,Tuple{UnitRange{Int},UnitRange{Int}}}
    has_key::Bool
    done::Bool
end

@generate_getters(DoorKeyUndirected)
@generate_setters(DoorKeyUndirected)

#####
# Directed
#####

function DoorKeyDirected(; T = Float32, height = 8, width = 8, rng = Random.GLOBAL_RNG)
    objects = (AGENT, WALL, GOAL, DOOR, KEY)
    world = GridWorldBase(objects, height, width)

    room = Room(CartesianIndex(1, 1), height, width)
    place_room!(world, room)

    goal_pos = CartesianIndex(height - 1, width - 1)
    world[GOAL, goal_pos] = true

    agent_pos = CartesianIndex(2, 2)
    world[AGENT, agent_pos] = true

    door_pos = CartesianIndex(2, 3)
    world[DOOR, door_pos] = true

    partition = CartesianIndices((2:height-1, door_pos[2]:door_pos[2]))
    world[WALL, partition] .= true
    world[WALL, door_pos] = false

    key_pos = CartesianIndex(3, 2)
    world[KEY, key_pos] = true

    agent_dir = RIGHT
    reward = zero(T)
    terminal_reward = one(T)
    has_key = false
    done = false

    env = DoorKeyDirected(world, agent_pos, agent_dir, reward, rng, terminal_reward, goal_pos, door_pos, key_pos, partition, has_key, done)

    RLBase.reset!(env)

    return env
end

RLBase.state_space(env::DoorKeyDirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = nothing
const DOOR_KEY_DIRECTED_LAYERS = SA.SVector(2, 3, 4, 5)
RLBase.state(env::DoorKeyDirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = get_grid(get_world(env), get_agent_pos(env), get_agent_dir(env), get_half_size(env), DOOR_KEY_DIRECTED_LAYERS)

RLBase.state_space(env::DoorKeyDirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = nothing
RLBase.state(env::DoorKeyDirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = (copy(get_grid(env)), get_agent_dir(env))

RLBase.action_space(env::DoorKeyDirected, ::RLBase.DefaultPlayer) = (DIRECTED_NAVIGATION_ACTIONS..., PICK_UP)
RLBase.reward(env::DoorKeyDirected, ::RLBase.DefaultPlayer) = get_reward(env)
RLBase.is_terminated(env::DoorKeyDirected) = get_done(env)

function RLBase.reset!(env::DoorKeyDirected{T}) where {T}
    world = get_world(env)
    height = get_height(env)
    width = get_width(env)
    rng = get_rng(env)

    world[AGENT, get_agent_pos(env)] = false
    world[WALL, get_partition(env)] .= false
    world[GOAL, get_goal_pos(env)] = false
    world[DOOR, get_door_pos(env)] = false
    world[KEY, get_key_pos(env)] = false

    new_door_pos, new_partition, region1, region2 = generate_partition(rng, height, width)
    region1, region2 = Random.shuffle(SA.SVector(region1, region2))

    set_door_pos!(env, new_door_pos)
    world[DOOR, new_door_pos] = true

    set_partition!(env, new_partition)
    world[WALL, new_partition] .= true
    world[WALL, new_door_pos] = false

    new_goal_pos = rand(rng, pos -> !any(@view world[:, pos]), region2)
    set_goal_pos!(env, new_goal_pos)
    world[GOAL, new_goal_pos] = true

    new_key_pos = rand(rng, pos -> !any(@view world[:, pos]), region1)
    set_key_pos!(env, new_key_pos)
    world[KEY, new_key_pos] = true

    new_agent_pos = rand(rng, pos -> !any(@view world[:, pos]), region1)
    set_agent_pos!(env, new_agent_pos)
    world[AGENT, new_agent_pos] = true

    set_agent_dir!(env, rand(rng, DIRECTIONS))

    set_reward!(env, zero(T))
    set_has_key!(env, false)
    set_done!(env, false)

    return nothing
end

function (env::DoorKeyDirected{T})(action::AbstractTurnAction) where {T}
    new_dir = turn(action, get_agent_dir(env))
    set_agent_dir!(env, new_dir)
    world = get_world(env)

    if world[GOAL, get_agent_pos(env)]
        set_done!(env, true)
        set_reward!(env, get_terminal_reward(env))
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end

function (env::DoorKeyDirected{T})(action::AbstractMoveAction) where {T}
    world = get_world(env)

    agent_pos = get_agent_pos(env)
    dest = move(action, get_agent_dir(env), agent_pos)
    if !(world[DOOR, dest] || world[WALL, dest]) || (world[DOOR, dest] && get_has_key(env))
        world[AGENT, agent_pos] = false
        set_agent_pos!(env, dest)
        world[AGENT, dest] = true
    end

    if world[GOAL, get_agent_pos(env)]
        set_done!(env, true)
        set_reward!(env, get_terminal_reward(env))
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end

#####
# Undirected
#####

function DoorKeyUndirected(; T = Float32, height = 8, width = 8, rng = Random.GLOBAL_RNG)
    objects = (AGENT, WALL, GOAL, DOOR, KEY)
    world = GridWorldBase(objects, height, width)

    room = Room(CartesianIndex(1, 1), height, width)
    place_room!(world, room)

    goal_pos = CartesianIndex(height - 1, width - 1)
    world[GOAL, goal_pos] = true

    agent_pos = CartesianIndex(2, 2)
    world[AGENT, agent_pos] = true

    door_pos = CartesianIndex(2, 3)
    world[DOOR, door_pos] = true

    partition = CartesianIndices((2:height-1, door_pos[2]:door_pos[2]))
    world[WALL, partition] .= true
    world[WALL, door_pos] = false

    key_pos = CartesianIndex(3, 2)
    world[KEY, key_pos] = true

    reward = zero(T)
    terminal_reward = one(T)
    has_key = false
    done = false

    env = DoorKeyUndirected(world, agent_pos, reward, rng, terminal_reward, goal_pos, door_pos, key_pos, partition, has_key, done)

    RLBase.reset!(env)

    return env
end

RLBase.state_space(env::DoorKeyUndirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = nothing
const DOOR_KEY_UNDIRECTED_LAYERS = SA.SVector(2, 3, 4, 5)
RLBase.state(env::DoorKeyUndirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = get_grid(get_world(env), get_agent_pos(env), get_half_size(env), DOOR_KEY_UNDIRECTED_LAYERS)

RLBase.state_space(env::DoorKeyUndirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = nothing
RLBase.state(env::DoorKeyUndirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = copy(get_grid(env))

RLBase.action_space(env::DoorKeyUndirected, player::RLBase.DefaultPlayer) = (UNDIRECTED_NAVIGATION_ACTIONS..., PICK_UP)
RLBase.reward(env::DoorKeyUndirected, ::RLBase.DefaultPlayer) = get_reward(env)
RLBase.is_terminated(env::DoorKeyUndirected) = get_done(env)

function RLBase.reset!(env::DoorKeyUndirected{T}) where {T}
    world = get_world(env)
    height = get_height(env)
    width = get_width(env)
    rng = get_rng(env)

    world[AGENT, get_agent_pos(env)] = false
    world[WALL, get_partition(env)] .= false
    world[GOAL, get_goal_pos(env)] = false
    world[DOOR, get_door_pos(env)] = false
    world[KEY, get_key_pos(env)] = false

    new_door_pos, new_partition, region1, region2 = generate_partition(rng, height, width)
    region1, region2 = Random.shuffle(SA.SVector(region1, region2))

    set_door_pos!(env, new_door_pos)
    world[DOOR, new_door_pos] = true

    set_partition!(env, new_partition)
    world[WALL, new_partition] .= true
    world[WALL, new_door_pos] = false

    new_goal_pos = rand(rng, pos -> !any(@view world[:, pos]), region2)
    set_goal_pos!(env, new_goal_pos)
    world[GOAL, new_goal_pos] = true

    new_key_pos = rand(rng, pos -> !any(@view world[:, pos]), region1)
    set_key_pos!(env, new_key_pos)
    world[KEY, new_key_pos] = true

    new_agent_pos = rand(rng, pos -> !any(@view world[:, pos]), region1)
    set_agent_pos!(env, new_agent_pos)
    world[AGENT, new_agent_pos] = true

    set_reward!(env, zero(T))
    set_has_key!(env, false)
    set_done!(env, false)

    return nothing
end

function (env::DoorKeyUndirected{T})(action::AbstractMoveAction) where {T}
    world = get_world(env)

    agent_pos = get_agent_pos(env)
    dest = move(action, agent_pos)
    if !(world[DOOR, dest] || world[WALL, dest]) || (world[DOOR, dest] && get_has_key(env))
        world[AGENT, agent_pos] = false
        set_agent_pos!(env, dest)
        world[AGENT, dest] = true
    end

    if world[GOAL, get_agent_pos(env)]
        set_done!(env, true)
        set_reward!(env, get_terminal_reward(env))
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end

#####
# Common
#####

function generate_partition(rng, height, width)
    if rand(rng) < 0.5
        # vertical partition
        door_pos = rand(rng, CartesianIndices((2:height-1, 3:width-2)))
        partition = CartesianIndices((2:height-1, door_pos[2]:door_pos[2]))
        region1 = CartesianIndices((2:height-1, 2:door_pos[2]-1))
        region2 = CartesianIndices((2:height-1, door_pos[2]+1:width-1))
        return door_pos, partition, region1, region2
    else
        # horizontal partition
        door_pos = rand(rng, CartesianIndices((3:height-2, 2:width-1)))
        partition = CartesianIndices((door_pos[1]:door_pos[1], 2:width-1))
        region1 = CartesianIndices((2:door_pos[1]-1, 2:width-1))
        region2 = CartesianIndices((door_pos[1]+1:height-1, 2:width-1))
        return door_pos, partition, region1, region2
    end
end

function (env::Union{DoorKeyDirected{T}, DoorKeyUndirected{T}})(::PickUp) where {T}
    world = get_world(env)
    agent_pos = get_agent_pos(env)

    if world[KEY, agent_pos]
        set_has_key!(env, true)
        world[KEY, agent_pos] = false
    end

    if world[GOAL, get_agent_pos(env)]
        set_done!(env, true)
        set_reward!(env, get_terminal_reward(env))
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end


text:
# functions to generate distance fields for stroke types
function circle{T}(uv::Vec{2,T})
    T(1)-norm(uv)
end

function dots(gap, resolution)
    X = linspace(0, gap, resolution[1])
    Y = linspace(0, 1, resolution[2])
    Float16[
        circle(Vec(x-0.5,y-0.5))
    for x=X, y=Y]
end

dist(a, b) = abs(a-b)
mindist(x, a, b) = min(dist(a, x), dist(b, x))
function gappy(x, ps)
    n = length(ps)
    x <= first(ps) && return first(ps) - x
    for j=1:(n-1)
        p0 = ps[j]
        p1 = ps[min(j+1, n)]
        if p0 <= x && p1 >= x
            return mindist(x, p0, p1) * (isodd(j) ? 1 : -1)
        end
    end
    return last(ps) - x
end
function ticks(points, resolution)
    Float16[gappy(x, points) for x=linspace(first(points),last(points), resolution)]
end

points = [0., 0.01, 0.02, 0.03, 0.09, 0.6, 0.8, 0.9]

text = Texture(ticks(points, 100), x_repeat=:repeat)


text:
export discreteApprox!, discreteApprox, discreteNormalApprox, discreteNormalApprox!

# ----------------------- objective functions for max entropy calcs --------------------------

function expΔTx!(tmpvec::Vector, ΔT::AbstractMatrix, x::AbstractVector)
  mul!(tmpvec, ΔT, x)
  tmpvec .= exp.(tmpvec)
end

# objective
function entropyObjective_f!(tmpvec::Vector, x::Vector, q::Vector, ΔT::AbstractMatrix)
  expΔTx!(tmpvec, ΔT, x)
  return dot(q, tmpvec)
end

# gradient
function entropyObjective_g!(grad::Vector, tmpvec::Vector, x::Vector, q::Vector, ΔT::AbstractMatrix)
  n,L = size(ΔT)
  expΔTx!(tmpvec, ΔT, x)
  tmpvec .*= q
  for l = 1:L
    grad[l] = dot(tmpvec, @view(ΔT[:,l]))
  end
end

function entropyObjective_fg!(grad::Vector, tmpvec::Vector, x::Vector, q::Vector, ΔT::AbstractMatrix)
  entropyObjective_g!(grad, tmpvec, x, q, ΔT)
  return sum(tmpvec)
end

function entropyObjective_h!(hess::Matrix{T}, tmpvec::Vector, x::Vector, q::Vector, ΔT::AbstractMatrix) where {T}
  n,L = size(ΔT)
  expΔTx!(tmpvec, ΔT, x)
  tmpvec .*= q
  fill!(hess, zero(T))
  for k = 1:L
    for l = 1:L
      hess[l,k] = sum(@view(ΔT[:,l]) .* tmpvec .* @view(ΔT[:,k]))
    end
  end
end

# ----------------------- wrappers --------------------------

function discreteApprox!(p::AbstractVector, λfinal::AbstractVector, err::AbstractVector, tmp::Vector, q0::Vector, ΔT::AbstractMatrix{T}) where {T}

  l = length(λfinal)
  n = length(q0)
  (n,l) == size(ΔT)  || throw(DimensionMismatch())
  n == length(p)     || throw(DimensionMismatch())
  l == length(err)   || throw(DimensionMismatch())
  n == length(tmp)   || throw(DimensionMismatch())

  # test that initial value is finite
  λ0 = zeros(T,l)
  grad = zeros(T,l)
  f0 = entropyObjective_fg!(grad, tmp, λ0, q0, ΔT)
  !isfinite(f0)         && return Inf
  !all(isfinite.(grad)) && return Inf

  tdf = TwiceDifferentiable(
    (x::Vector)               -> entropyObjective_f!(       tmp, x, q0, ΔT),
    (grad::Vector, x::Vector) -> entropyObjective_g!( grad, tmp, x, q0, ΔT),
    (grad::Vector, x::Vector) -> entropyObjective_fg!(grad, tmp, x, q0, ΔT),
    (hess::Matrix, x::Vector) -> entropyObjective_h!( hess, tmp, x, q0, ΔT),
    λ0
  )

  opt = Optim.optimize(tdf, λ0, Newton())
  λ1 = opt.minimizer
  J = opt.minimum

  # update gradient
  entropyObjective_g!(grad, tmp, λ1, q0, ΔT)

  if norm(grad, Inf) < 1e-9 && all(isfinite.(grad)) && all(isfinite.(λ1)) && 0.0 < J < Inf && maximum(abs.(grad ./ J)) < 1e-5
    λfinal .= λ1
    expΔTx!(tmp, ΔT, λ1)
    p .= q0 .* tmp ./ J
    err .= grad ./ J
    return J
  else
    return Inf
  end

end

function ΔTmat!(ΔT::Matrix, dev::AbstractVector, Tbar::Vector)
  n,L = size(ΔT)
  length(dev) == n || throw(DimensionMismatch())
  length(Tbar) == L || throw(DimensionMismatch())

  ΔT[:,1] .= dev
  for j = 2:L
    ΔT[:,j] .= ΔT[:,j-1] .* dev
  end
  for j = 1:L
    ΔT[:,j] .-= Tbar[j]
  end
end

function discreteApprox!(P::AbstractMatrix, y::AbstractVector{T}, S::Union{AbstractVector, Base.Iterators.ProductIterator}, zval::Function, pdffun::Function, scaled_moments::Vector, scale_factor::Real, maxMoments::Integer, κ::Real) where {T<:Real}

  nS = length(S)
  n = length(y)
  0 < maxMoments < n || throw(error("Must use 1 to $n-1 moments or fewer"))
  (nS,n) == size(P)  || throw(DimensionMismatch())

  # Initialize elements that will be returned
  Λ          = zeros(T, nS, maxMoments)
  JN         = zeros(T, nS)
  approxErr  = zeros(T, nS, maxMoments)
  numMoments = zeros(Int, nS)

  # preallocate these, which will be updated each iteration
  ΔT  = Array{T}(undef,n, maxMoments)
  z   = Array{T}(undef,n)
  q   = Array{T}(undef,n)
  tmp = Array{T}(undef,n)

  for (i,st) in enumerate(S)
    z .= zval.(y, st)
    q .= max.(pdffun.(z), κ)
    z ./= scale_factor
    ΔTmat!(ΔT, z, scaled_moments)
    updated = false
    for l in maxMoments:-1:2
      J = discreteApprox!(@view(P[i,:]), @view(Λ[i,1:l]), @view(approxErr[i,1:l]), tmp, q, @view(ΔT[:,1:l]))
      if isfinite(J)
        JN[i], numMoments[i] = (J, l)
        updated = true
        break
      end
    end
    if !updated
      sumq = sum(q)
      P[i,:] .= q ./ sumq
    end

  end

  return JN, Λ, numMoments, approxErr
end

# ----------------------- wrappers --------------------------

function discreteApprox(y::AbstractVector{T}, S::Union{AbstractVector, Base.Iterators.ProductIterator}, zval::Function, pdffun::Function, scaled_moments::Vector, scale_factor::Real, maxMoments::Integer, κ::Real) where {T<:Real}
  n = length(y)
  nS = length(S)
  P = Array{T}(nS,n)
  return (P, discreteApprox!(P, y, S, zval, pdffun, scaled_moments, scale_factor, maxMoments, κ)...)
end

function discreteNormalApprox!(P::AbstractMatrix, y::AbstractVector, S::Union{AbstractVector, Base.Iterators.ProductIterator}, zval::Function, maxMoments::Integer, κ::Real)
  scale_factor = maximum(abs.(y))
  scaled_moments = [m for m in NormCentralMoment(maxMoments, 1.0/scale_factor)]
  discreteApprox!(P, y, S, zval, normpdf, scaled_moments, scale_factor, maxMoments, κ)
end

function discreteNormalApprox(y::AbstractVector{T}, S::Union{AbstractVector, Base.Iterators.ProductIterator}, zval::Function, maxMoments::Integer=2, κ::Real=1e-8) where {T}
  n = length(y)
  P = Matrix{T}(undef,n,n)
  out = discreteNormalApprox!(P, y, S, zval, maxMoments, κ)
  return (P, out...)
end


text:
"""
Run gibbs sampler.
"""
function gibbs_sample!(
        model::SeqModel,
        spikes::Vector{Spike},
        initial_assignments::Vector{Int64},
        num_samples::Int64,
        extra_split_merge_moves::Int64,
        split_merge_window::Float64,
        save_every::Int64;
        verbose::Bool=false
    )

    # Initialize spike assignments.
    assignments = initial_assignments
    recompute!(model, spikes, assignments)

    # ======== THINGS TO SAVE ======== #

    # number of samples to save.
    n_saved_samples = Int(round(num_samples / save_every))

    # log likelihoods.
    log_p_hist = zeros(n_saved_samples)

    # parent assignments.
    assignment_hist = zeros(
        Int64, length(spikes), n_saved_samples
    )

    # latent event assignment ids, times, types, amplitudes.
    latent_event_hist = Vector{EventSummaryInfo}[]

    # global variables (offsets, response amplitudes, etc.)
    globals_hist = SeqGlobals[]

    # Order to iterate over spikes.
    spike_order = collect(1:length(spikes))

    # ======== MAIN LOOP ======== #

    for s = 1:num_samples

        # Update spike assignments in random order.
        Random.shuffle!(spike_order)
        for i = spike_order
            remove_datapoint!(model, spikes[i], assignments[i])
            assignments[i] = gibbs_add_datapoint!(model, spikes[i])
        end

        # Add extra split merge moves.
        split_merge_sample!(
            model,
            spikes,
            extra_split_merge_moves,
            assignments,
            split_merge_window,
        )

        # Update latent events.
        gibbs_update_latents!(model)

        # Update globals
        gibbs_update_globals!(model, spikes, assignments)

        # Store results
        if (s % save_every) == 0

            # Index into hist vectors.
            j = Int(s / save_every)

            # Note that the likelihood should be computed after
            # gibbs_update_latents! and before spike parent assignments
            # are recomputed. Otherwise, newly added latent events are
            # not initialized.
            log_p_hist[j] = log_like(model, spikes)

            # Save assignments.
            assignment_hist[:, j] .= assignments

            # Save latent event information.
            push!(
                latent_event_hist,
                event_list_summary(model)
            )

            # Save global variables.
            push!(globals_hist, deepcopy(model.globals))

            # Display progress.
            verbose && print(s, "-")
        end
    end

    # Finish progress bar.
    verbose && (n_saved_samples > 0) && println("Done")

    return (
        assignments,
        assignment_hist,
        log_p_hist,
        latent_event_hist,
        globals_hist
    )
end

"""
Adds spikes `s` to an existing sequence event, to a new sequence event,
or to the background process.

For each sequence event k = 1 ... K, we compute

    prob[k] = p(x_i | z_i = k, x_{neq i}) * (N_k + alpha)

The probability of forming a new cluster is

    prob[K + 1] propto p(x_i | z_i = K + 1) * alpha * (V(K + 1) / V(K))

where p(x_i | z_i = K + 1) is the marginal probability of a singleton cluster.
See section 6 of Miller & Harrison (2018).

The probability of the background is

    prob[K + 2] propto p(x_i | bkgd) * lambda0 * (1 + beta)

where lambda0 = m.bkgd_rate and (alpha, beta) are the shape and
rate parameters of the gamma prior on sequence event amplitudes.
"""
function gibbs_add_datapoint!(model::SeqModel, x::Spike)

    # Create log-probability vector to sample assignments.
    #
    #  - We need to sample K + 2 possibilities. There are K existing clusters
    #    we could assign to. We could also form a new cluster (index K + 1),
    #    or assign the spike to the background (index K + 2).

    # Shape and rate parameters of gamma prior on latent event amplitude.
    α = model.priors.seq_event_amplitude.α
    β = model.priors.seq_event_amplitude.β

    K = num_sequence_events(model)

    # Grab vector without allocating new memory.
    log_probs = resize!(model._K_buffer, K + 2)

    # Iterate over model events, indexed by k = {1, 2, ..., K}.
    for (k, event) in enumerate(model.sequence_events)

        # Check if event is too far away to be considered. If sampled_type < 0,
        # then the event timestamp hasn't been sampled yet, so we can't give
        # up yet. Be aware that this previously led to a subtle bug, which was
        # very painful to fix.
        too_far = abs(x.timestamp - event.sampled_timestamp) > model.max_sequence_length
        if too_far && (event.sampled_type > 0)
            log_probs[k] = -Inf

        # Compute probability of adding spike to cluster k.
        else
            log_probs[k] = (
                log(event.spike_count + α)
                + log_posterior_predictive(event, x, model)
            )
        end
    end

    # New cluster probability.
    log_probs[K + 1] = (
        model.new_cluster_log_prob
        + log_posterior_predictive(x, model)
    )

    # Background probability
    log_probs[K + 2] = (
        model.bkgd_log_prob
        + model.globals.bkgd_log_proportions[x.neuron]
        - log(model.max_time)
    )

    # Sample new assignment for spike x.
    z = sample_logprobs!(log_probs)

    # New sample corresponds to background, do nothing.
    if z == (K + 2)
        return -1

    # New sample corresponds to new sequence event / cluster.
    elseif z == (K + 1)
        return add_event!(model, x)  # returns new assignment

    # Otherwise, add datapoint to existing sequence event. Note
    # that z is an integer in [1:K], while assignment indices
    # can be larger and non-contiguous.
    else
        k = model.sequence_events.indices[z]  # look up assignment index.
        return add_datapoint!(model, x, k)
    end

end

function gibbs_update_globals!(
        model::SeqModel,
        spikes::Vector{Spike},
        assignments::AbstractVector{Int64}
    )

    K = num_sequence_events(model)
    N = num_neurons(model)
    R = num_sequence_types(model)

    priors = model.priors
    globals = model.globals

    # === RESAMPLE BACKGROUND SPIKE PARAMETERS === #

    num_bkgd_spikes = 0
    bkgd_counts = zeros(N)

    for (i, x) in enumerate(spikes)
        if assignments[i] < 0
            num_bkgd_spikes += 1
            bkgd_counts[x.neuron] += 1
        end
    end

    # Sample proportions of background spikes across neurons
    # as a probability vector then map to log space for
    # future computations.
    rand!(
        posterior(
            bkgd_counts, priors.bkgd_proportions
        ),
        globals.bkgd_log_proportions
    ) # multinomial - symmetric dirichlet conjugate pair.
    map!(
        log,
        globals.bkgd_log_proportions,
        globals.bkgd_log_proportions
    )

    # Sample the rate of background spikes. Here, we need to adjust
    # for the length of the time interval, max_time. The scaling
    # property of the gamma distribution implies that:
    #
    #   bkgd_amp * max_time ~ RateGamma(α_bkgd, β_bkgd / max_time)
    #
    # If we observe num_bkgd, then the posterior (by Poisson-gamma
    # conjugacy) is:
    #
    #   bkgd_amp * T | num_bkgd_spikes ~ RateGamma(num_bkgd + α_bkgd, 1 + β_bkgd / T)
    #
    # Now apply the gamma scaling property again, dividing by T this
    # time, so we get:
    #
    #   bkgd_amp | num_bkgd_spikes ~ RateGamma(num_bkgd + α_bkgd, T + β_bkgd)
    #
    globals.bkgd_amplitude = rand(
        RateGamma(
            num_bkgd_spikes + priors.bkgd_amplitude.α,
            priors.bkgd_amplitude.β + model.max_time
        )
    )

    # === RESAMPLE SEQUENCE TYPE PROPORTIONS === #

    seq_type_counts = zeros(R)
    for event in model.sequence_events
        seq_type_counts[event.sampled_type] += 1
    end

    rand!(
        posterior(
            seq_type_counts, priors.seq_type_proportions
        ),
        model.globals.seq_type_log_proportions
    )  # multinomial - symmetric dirichlet conjugate pair.
    map!(
        log,
        globals.seq_type_log_proportions,
        globals.seq_type_log_proportions
    )

    # === RESAMPLE NEURON RESPONSE PROFILES === #

    # TODO -- preallocate arrays for this? Or is it not worth it?
    spk_count = zeros(Int64, N, R)
    spk_fm = zeros(N, R)
    spk_sm = zeros(N, R)

    for (i, x) = enumerate(spikes)

        # spike assignment to latent event.
        k = assignments[i]

        # skip if spike is assigned to background.
        (k < 0) && continue

        event = model.sequence_events[k]
        n = x.neuron
        r = event.sampled_type
        w = event.sampled_warp
        offset = (x.timestamp - event.sampled_timestamp) / model.priors.warp_values[w]

        # compute sufficient statistics
        spk_count[n, r] += 1
        spk_fm[n, r] += offset
        spk_sm[n, r] += offset * offset

    end

    for r = 1:R

        rand!(
            posterior(
                view(spk_count, :, r),
                priors.neuron_response_proportions
            ),
            view(globals.neuron_response_log_proportions, :, r)
        ) # multinomial - dirichlet conjugate pair.

        for n = 1:N

            (
                globals.neuron_response_offsets[n, r],
                globals.neuron_response_widths[n, r]
            ) =
            rand(
                posterior(
                    spk_count[n, r],
                    spk_fm[n, r],
                    spk_sm[n, r],
                    priors.neuron_response_profile
                )
            ) # normal - norm-inv-chi-squared conjugate pair.

        end

    end

    map!(
        log,
        globals.neuron_response_log_proportions,
        globals.neuron_response_log_proportions
    )

    # === RECOMPUTE NEW CLUSTER PROBABILITIES === #

    α = priors.seq_event_amplitude.α
    β = priors.seq_event_amplitude.β

    # TODO : make this a field of globals, not model.
    model.new_cluster_log_prob = (
        log(α)
        + log(model.priors.seq_event_rate)  # TODO: resample this as well.
        + log(model.max_time)
        + α * (log(β) - log(1 + β))
    )

    # TODO : make this a field of globals, not model.
    model.bkgd_log_prob = (
        log(globals.bkgd_amplitude)
        + log(model.max_time)
        + log(1 + β)
    )

    # === RECOMPUTE SUFFICIENT STATISTICS === #

    # Because the sequence types have been resampled, the sufficient statistics
    # computed for each sequence event must be re-computed (they depend on the
    # neuron offsets, relative amplitudes, etc.!)

    recompute!(model, spikes, assignments)

end

"""
Resamples sequence type, timestamp, and amplitude. For all
latent events.
"""
function gibbs_update_latents!(model::SeqModel)

    # Grab length-R vector (already pre-allocated).
    log_probs = model._RW_buffer

    for event in model.sequence_events

        # We should only be resampling non-empty events.
        @assert (event.spike_count > 0)

        # Sample sequence type.
        log_probs .= event.seq_type_posterior
        ind = sample_logprobs!(vec(log_probs))
        r, w = Tuple(CartesianIndices(size(log_probs))[ind])
        event.sampled_type = r
        event.sampled_warp = w

        # Sample event time, t ~ N(μ, σ2), given sequence type.
        σ2 = 1 / event.summed_precisions[r, w]
        μ = event.summed_potentials[r, w] * σ2
        event.sampled_timestamp = μ + randn() * sqrt(σ2)

        # Sample event amplitude.
        event.sampled_amplitude = rand(
            posterior(
                event.spike_count,
                model.priors.seq_event_amplitude
            )
        ) # Poisson - gamma conjugate pair.

    end

end


text:
using BetterFileWatching

using Test

function poll(query::Function, timeout::Real=Inf64, interval::Real=1/20)
    start = time()
    while time() < start + timeout
        if query()
            return true
        end
        sleep(interval)
    end
    return false
end

"Like @async except it prints errors to the terminal. 👶"
macro asynclog(expr)
	quote
		@async begin
			# because this is being run asynchronously, we need to catch exceptions manually
			try
				$(esc(expr))
			catch ex
				bt = stacktrace(catch_backtrace())
				showerror(stderr, ex, bt)
				rethrow(ex)
			end
		end
	end
end

@testset "Basic - $(method)" for method in ["new", "legacy"]
    test_dir = tempname(cleanup=false)
    mkpath(test_dir)

    j(args...) = joinpath(test_dir, args...)

    write(j("script.jl"), "nice(123)")

    mkdir(j("mapje"))
    write(j("mapje", "cool.txt"), "hello")

    events = []
    last_length = Ref(0)

    watch_task = if method == "new"
        @asynclog watch_folder(test_dir) do e
            @info "Event" e
            push!(events, e)
        end
    else
        @asynclog while true
            e = watch_folder(test_dir)
            @info "Event" e
            push!(events, e)
        end
    end

    sleep(2)

    @test length(events) == 0

    function somethingdetected()
        result = poll(3, 1/100) do
            length(events) > last_length[]
        end
        sleep(.5)
        last_length[] = length(events)
        result
    end

    write(j("mapje", "cool2.txt"), "hello")
    @test somethingdetected()

    write(j("mapje", "cool2.txt"), "hello again!")
    @test somethingdetected()

    write(j("asdf.txt"), "")
    @test somethingdetected()

    mv(j("mapje", "cool2.txt"), j("mapje", "cool3.txt"))
    @test somethingdetected()

    rm(j("mapje", "cool3.txt"))
    @test somethingdetected()

    sleep(2)
    @test_nowarn schedule(watch_task, InterruptException(); error=true)
end


text:
using SustainableFinance
using Test

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

max_stars_repo_name: tlorans/SustainableFinance

max_issues_repo_name: tlorans/SustainableFinance

max_forks_repo_name: tlorans/SustainableFinance


text:
include("discrete_laplacian.jl")

using LinearAlgebra
using Test

@testset "ALS linsolv" begin
    n = 10
    L = Lap(n,3)
    b = ones(n,n,n)
    x0 = randn(n,n,n)
    L_tt = tto_decomp(reshape(L,n,n,n,n,n,n))
    b_tt = ttv_decomp(b)
    x0_tt = ttv_decomp(x0)
    x = L\b[:]
    x_tt = als_linsolv(L_tt,b_tt,x0_tt)
    y = ttv_to_tensor(x_tt)
    @test isapprox(y[:],x)
end

@testset "ALS eigsolv" begin
    n = 10
    L = Lap(n,3)
    E = eigmin(L)
    x0 = randn(n,n,n)
    L_tt = tto_decomp(reshape(L,n,n,n,n,n,n))
    x0_tt = ttv_decomp(x0)
    E_tt,x_tt = als_eigsolv(L_tt,x0_tt)
    y = ttv_to_tensor(x_tt)
    @test isapprox(E,E_tt[end])
    #with the iterative solver
    E_tt,x_tt = als_eigsolv(L_tt,x0_tt,it_solver=true)
    @test isapprox(E,E_tt[end],atol=1e-8)
    #with scheduler
    x0_tt = ttv_decomp(ones(n,n,n))
    E_tt,x_tt = als_eigsolv(L_tt,x0_tt;sweep_schedule=[2,3,4],rmax_schedule=[5,7,10],noise_schedule=[1e-3,1e-4,0.0])
    y = ttv_to_tensor(x_tt)
    @test isapprox(E,E_tt[end])
end

@testset "ALS Generalized eigsolv" begin
    n = 5
    L = randn(n^3,n^3)
    L = L*L'+I
    S = randn(n^3,n^3)
    S = S*S'+I
    F = eigen(L,S)
    x0 = randn(n,n,n)
    L_tt = tto_decomp(reshape(L,n,n,n,n,n,n))
    S_tt = tto_decomp(reshape(S,n,n,n,n,n,n))
    x0_tt = ttv_decomp(x0)
    E_tt,x_tt = als_gen_eigsolv(L_tt,S_tt,x0_tt)
    y = ttv_to_tensor(x_tt)
    @test isapprox(real(F.values[1]),E_tt[end])
    E_tt,x_tt = als_gen_eigsolv(L_tt,S_tt,x0_tt,it_solver=true)
    @test isapprox(real(F.values[1]),E_tt[end],atol=1e-8)
end

max_stars_repo_name: msdupuy/Tensor-Train-Julia

max_issues_repo_name: msdupuy/Tensor-Train-Julia

max_forks_repo_name: msdupuy/Tensor-Train-Julia


text:
export ℍ
export 💞
export 🌐
export cartesian
export geographic
export sphere
export ℝ³
export ℂ

struct ℍ
    su2::Array{Complex,2}
    q::Array{Float64}
end

struct 🌐
    basemanifold::Array{Float64,3}
    basecolor::Array{Float64,3}
    basecenter::Array{Float64}
    markermanifold::Array{Float64,4}
    markercolor::Array{Float64,4}
    markercenter::Array{Float64,2}
    markerradius::Float64
    segments::Integer
end

ℍ(q::Array{Float64}) = begin
    z₁, z₂ = Complex(q[1], q[2]), Complex(q[3], q[4])
    s = [z₁ conj(z₂); -z₂ conj(z₁)]
    ℍ(s, q)
end

ℍ(s::Array{Complex,2}) = begin
    q = [real(s[1]); imag(s[1]); -real(s[3]); -imag(s[3])]
    ℍ(s, q)
end

Base.imag(h::ℍ) = [imag(h.su2[1]); -real(h.su2[3]); -imag(h.su2[3])]

Base.adjoint(h::ℍ) = ℍ(convert(Array{Complex,2}, adjoint(h.su2)))

function Base.:*(h₁::ℍ, h₂::ℍ)
    s = h₁.su2 * h₂.su2
    q = [real(s[1]); imag(s[1]); -real(s[3]); -imag(s[3])]
    ℍ(s, q)
end

💞(p::Array{Float64}, h::ℍ) = imag(adjoint(h) * ℍ([0; p]) * h)

function 💞(p::Array{Float64,2}, h::ℍ)
    r = similar(p)
    for i in 1:size(p, 1)
        r[i, :] = 💞(p[i, :], h)
    end
    r
end

function 💞(p::Array{Float64,3}, center::Array{Float64}, h::ℍ)
    r = similar(p)
    c = reshape(repeat(center', size(p, 2)), size(p, 2), size(p, 3))
    for i in 1:size(p, 1)
        oldposition = p[i, :, :] - c
        r[i, :, :] = 💞(oldposition, h) + c
    end
    r
end

ℝ³(b::Array{Complex}) = begin
    x, y = real.(b), imag.(b)
    d = x .^ 2 .+ y .^ 2 .+ 1
    [2x ./ d 2y ./ d (d .- 2) ./ d]
end

cartesian(x::Array{Float64,2}) = begin
    ϕ = x[:, 1]
    θ = x[:, 2]
    r = x[:, 3]
    y₁ = r .* cos.(θ) .* cos.(ϕ)
    y₂ = r .* cos.(θ) .* sin.(ϕ)
    y₃ = r .* sin.(θ)
    [y₁ y₂ y₃]
end

function geographic(y)
    samples = size(y, 1)
    g = Array{Float64,2}(undef, samples, 3)
    for i in 1:samples
        if y[i, 1] > 0
            ϕ = atan(y[i, 2] / y[i, 1])
        elseif y[i, 2] > 0
            ϕ = atan(y[i, 2] / y[i, 1]) + pi
        else
            ϕ = atan(y[i, 2] / y[i, 1]) - pi
        end
        r = LinearAlgebra.norm(y[i, :])
        θ = asin(y[i, 3] / r)
        g[i, :] = [ϕ; θ; r]
    end
    g
end

function geographic(b::Array{Complex})
    geographic(ℝ³(b))
end

function ℂ(p::Array{Float64,2})
    xyz = cartesian(p)
    d = 1 .- xyz[:, 3]
    x, y = xyz[:, 1] ./ d, xyz[:, 2] ./ d
    convert(Array{Complex}, Complex.(x, y))
end

function sphere(center, radius, segments)
    manifold = Array{Float64,3}(undef, segments, segments, 3)
    for i in 1:segments
        θ = (i - 1) / (segments - 1) * pi - pi / 2
        for j in 1:segments
            ϕ = (j - 1) / (segments - 1) * 2pi - pi
            manifold[i, j, :] = vec(cartesian([ϕ -θ radius])) + center
        end
    end
    manifold
end

function 🌐(basecenter::Array{Float64},
            basecolor::Array{Float64},
            baseradius::Float64,
            markercenter::Array{Float64,2},
            markercolor2::Array{Float64,2},
            markerradius::Float64,
            segments::Integer)
    basemanifold = sphere(basecenter, baseradius - markerradius, segments)
    basecolor = reshape(repeat(basecolor', segments^2), segments, segments, 3)
    samples = size(markercenter, 1)
    markermanifold = Array{Float64,4}(undef, samples, segments, segments, 3)
    markercolor4 = similar(markermanifold)
    for i in 1:samples
        markermanifold[i, :, :, :] = sphere((markercenter[i, :] .* baseradius) + basecenter,
                                            markerradius,
                                            segments)
        markercolor4[i, :, :, :] = reshape(repeat(markercolor2[i, :]', segments^2),
                                           segments,
                                           segments,
                                           3)
    end
    🌐(basemanifold,
       basecolor,
       basecenter,
       markermanifold,
       markercolor4,
       markercenter,
       markerradius,
       segments)
end


text:
# Test runner

include("unit.jl")


text:
module TrussTopOptProblems

using Ferrite, StaticArrays, LinearAlgebra
using SparseArrays
using ..TopOpt
using ..TopOpt.Utilities
using Setfield
# using ..TopOpt: PENALTY_BEFORE_INTERPOLATION
# using ..Utilities: @forward_property

abstract type AbstractFEAMaterial end
struct TrussFEAMaterial{T} <: AbstractFEAMaterial
    E::T # Young's modulus
    ν::T # Poisson's ratio
end

abstract type AbstractFEACrossSec end
struct TrussFEACrossSec{T} <: AbstractFEACrossSec
    A::T # cross section area
end

import Ferrite: assemble!

# TODO apply @params to all the problem struct

include("grids.jl")
include("problem_types.jl")
include("matrices_and_vectors.jl")
include("elementinfo.jl")
include("buckling.jl")
include(joinpath("TrussIO", "TrussIO.jl"))
using .TrussIO
include(joinpath("TrussVisualization", "TrussVisualization.jl"))
using .TrussVisualization

export TrussGrid, TrussProblem, TrussFEACrossSec, TrussFEAMaterial
export parse_truss_json

end # module


text:
function calclp_dic(d, EL, ER, SL, SR)
    logmu = logmean(d)
    ll = zero(EL[1])
    for i in eachindex(EL)
        if isfinite(EL[i])
            if ER[i] < SL[i]
                ll += logmu+log(eqcdf(d,SR[i]-ER[i])-eqcdf(d,SL[i]-ER[i])-(eqcdf(d,SR[i]-EL[i])-eqcdf(d,SL[i]-EL[i]))) -
                 log(ER[i] - EL[i]) - log(SR[i] - SL[i])
            elseif ER[i] < SR[i]
                ll += logmu+log(eqcdf(d,SR[i]-ER[i])-(eqcdf(d,SR[i]-EL[i])-eqcdf(d,SL[i]-EL[i]))) -
                 log(ER[i] - EL[i]) - log(SR[i] - SL[i])
            elseif SR[i] <= ER[i]
                ll += logmu+log(eqcdf(d,SR[i]-EL[i]) - eqcdf(d,SL[i]-EL[i])) -
                  log(ER[i] - EL[i]) - log(SR[i] - SL[i])
            end
        else
            if ER[i] < SL[i]
                ll += log(eqcdf(d,SR[i]-ER[i])-eqcdf(d,SL[i]-ER[i])) - log(SR[i] - SL[i])
            elseif ER[i] < SR[i]
                ll += log(eqcdf(d,SR[i]-ER[i])) - log(SR[i] - SL[i])
            end
        end
    end
    return -ll
end

# function calclp_dicrt(d, EL, ER, SL, SR, Tmax)
#     n = length(EL)
#     ll = 0
#     logmu = logmean(d)
#     for i in 1:n
#         if ER[i]<SL[i]
#             ll += logmu+logsubexp(log(eqcdf(d,SR[i]-ER[i])-eqcdf(d,SL[i]-ER[i])),log(eqcdf(d,SR[i]-EL[i])-eqcdf(d,SL[i]-EL[i])))
#         else
#             ll += logmu+logsubexp(log(eqcdf(d,SR[i]-ER[i])-(SL[i]-ER[i])),log(eqcdf(d,SR[i]-EL[i])-eqcdf(d,SL[i]-EL[i])))
#         end
#         ll += log(cdf2(d, Tmax-SR[i]))
#     end
#     return -ll
# end

#######
#interval censored
function calclp_ic(d, EL, ER, S)
    n = length(S)
    ll = 0.0
    mu = mean(d)
    for i in 1:n
        if isfinite(EL[i])
            ll += logdiffcdf(d, S[i]-EL[i], S[i]-ER[i]) - log(ER[i] - EL[i])
        else
            ll += logccdf(d, S[i]-ER[i]) - logmean(d)
        end
    end
    return -ll
end

function calclp_icrt(d, EL, ER, S, Tmax)
    n = length(S)
    ll = 0.0
    for i in 1:n
        if isfinite(EL[i])
            ll += log(ccdf(d,S[i]-ER[i]) - ccdf(d,S[i]-EL[i])) - (logmean(d)+log(eqcdf(d,Tmax-ER[i])-eqcdf(d,Tmax-EL[i])))
        else
            ll += log(ccdf(d,S[i]-ER[i])) - log(eqcdf(d,Tmax-ER[i]))
        end
    end
    return -ll
end

max_stars_repo_name: abikoushi/IntervalCensored

max_issues_repo_name: abikoushi/IntervalCensored

max_forks_repo_name: abikoushi/IntervalCensored


text:
abstract type DiagnosticsGroupParams end

"""
    DiagnosticsGroup

Holds a set of diagnostics that share a collection interval, a filename
prefix, an output writer, an interpolation, and any extra parameters.
"""
mutable struct DiagnosticsGroup{DGP <: Union{Nothing, DiagnosticsGroupParams}}
    name::String
    init::Function
    fini::Function
    collect::Function
    interval::String
    out_prefix::String
    writer::AbstractWriter
    interpol::Union{Nothing, InterpolationTopology}
    params::DGP

    DiagnosticsGroup(
        name,
        init,
        fini,
        collect,
        interval,
        out_prefix,
        writer,
        interpol,
        params = nothing,
    ) = new{typeof(params)}(
        name,
        init,
        fini,
        collect,
        interval,
        out_prefix,
        writer,
        interpol,
        params,
    )
end

function GenericCallbacks.init!(dgngrp::DiagnosticsGroup, solver, Q, param, t)
    @info @sprintf(
        """
    Diagnostics: %s
        initializing at %8.2f""",
        dgngrp.name,
        t,
    )
    dgngrp.init(dgngrp, t)
    dgngrp.collect(dgngrp, t)
    return nothing
end
function GenericCallbacks.call!(dgngrp::DiagnosticsGroup, solver, Q, param, t)
    @tic diagnostics
    @info @sprintf(
        """
    Diagnostics: %s
        collecting at %8.2f""",
        dgngrp.name,
        t,
    )
    dgngrp.collect(dgngrp, t)
    @toc diagnostics
    return nothing
end
function GenericCallbacks.fini!(dgngrp::DiagnosticsGroup, solver, Q, param, t)
    @info @sprintf(
        """
    Diagnostics: %s
        finishing at %8.2f""",
        dgngrp.name,
        t,
    )
    dgngrp.collect(dgngrp, t)
    dgngrp.fini(dgngrp, t)
    return nothing
end

include("atmos_les_default.jl")
include("atmos_gcm_default.jl")
include("atmos_les_core.jl")
include("atmos_les_default_perturbations.jl")
include("atmos_refstate_perturbations.jl")
include("atmos_turbulence_stats.jl")
include("atmos_mass_energy_loss.jl")
include("dump_init.jl")
include("dump_state.jl")
include("dump_aux.jl")
include("dump_spectra.jl")

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
struct UncertainConstraint <: RobustConstraintType
    name::String
end

struct UncertainConst{S}  <: RobustConstraint{S}
    name::String
    func::GenericAffExpr
    set::S
end
export UncertainConstraint,UncertainConst


text:
# The entire non-mutating rod sequence.
# For organization, the entire front half accepts arrays and
# returns tuples of arrays. Convenience overloads are at the
# end of the file

"""
        cross3(a,b)

Cross products of the columns of two 3 x N arrays. Output is also 3 x N
"""
function cross3(a,b)
    hcat(a[:,2].*b[:,3] .- a[:,3].*b[:,2],
         a[:,3].*b[:,1] .- a[:,1].*b[:,3],
         a[:,1].*b[:,2] .- a[:,2].*b[:,1])
end

"""
        dot3(a,b)

Dot products of the columns of two 3 x N arrays. Output is 1 x N
"""
function dot3(a,b)
    sum(a.*b, dims=2)
end

"""
        norm3(a,b)

Norms of the columns of two 3 x N arrays. Output is 1 x N
"""
function norm3(a)
    sqrt.(sum(abs2, a, dims=2))
end

"""
        edges(x)

Edges of a centerline `x`, stored as a 3 x N array. Output is 3 x N-1
"""
function edges(x)
    x1 = @view x[2:end,:]
    x2 = @view x[1:end-1,:]
    x1 .- x2
end

"""
        tangents(x)

Tangents of a centerline `x`, stored as a 3 x N array. Output is 3 x N-1
"""
function tangents(x)
    e = edges(x)
    e./norm3(e)
end

"""
        ptransport(v, t1, t2)

Parallel transport of vectors `v` from tangents `t1` to tangents 't2'.

All inputs are (3 x N). Each column of `t1` and `t2` must have unit norm.

Uses the Rodriguez rotation formula. There is a singularity in this definition
when t1 == t2, which is handled by switching to an small-θ expansion
(see ptransport_inner_coefficient).

"""
function ptransport(v, t1, t2)
    X = cross3(t1,t2);
    s2 = dot3(X,X);
    c = dot3(t1,t2);
    ptic = ptransport_inner_coefficient.(s2,c)
    c.*v .+ cross3(X,v) .+ ptic.*dot3(X,v).*X
end

function ptransport_inner_coefficient(s2,c)
    opt1 = (1-c)/s2
    opt2 = 1/2 + (1-c)/4 + (1-c)^2/8*(1 + s2/4)
    s2 > PTRANSPORT_TOL^2 ?  opt1 : opt2
end

"""
        rotate_orthogonal_unit(v, t, θ)

Rotate unit vectors `v` around axes `t1` by angles `θ`.

`v` and `t` are (3 x N), `θ` is (1 x N).

The columns `v[:,i]` and `t[:,i]` must be orthogonal unit vectors. The function
uses the Rodriguez rotation formula speciailized to this case.
"""
function rotate_orthogonal_unit(v, t, θ)
    cos.(θ).*v .+ sin.(θ).*cross3(t,v)
end

"""
        full_kinematics(x, d)

Returns the edge lengths, curvatures, and twists of a rod.

# Arguments
- `x`: centerline vertices, size (3 x N)
- `d`: directors, size (3 x N-1). Unit length and orthogonal to `edges(x)`

# Returns
- `l`: edge lengths (1 x N-1)
- `κ`: curvature normals (2 x N-2), expressed in the director frame
- 'τ': scalar twist at each vertex (1 x N-2)

Twists are calculated by parallel transporting the directors across each
vertex, and comparing with the original. This avoids the reference
twist calculations used in Lestringant, Audoly, and Kochmann (2020)
"""
function full_kinematics(x, d1)

    # Tangents and directors
    nv = size(x,1)
    #e =x[:, 2:nv] .- x[:, 1:nv-1]
    e = edges(x)
    l = norm3(e)
    t = e./l
    d2 = cross3(t, d1)

    # Split into halves of the rod
    t1 = @view t[1:nv-2,:]
    t2 = @view t[2:nv-1,:]
    d1l = @view d1[1:nv-2,:]
    d1r = @view d1[2:nv-1,:]
    d2l = @view d2[1:nv-2,:]
    d2r = @view d2[2:nv-1,:]

    # Curvatures
    k = 2 .*cross3(t1, t2)./(1.0 .+ dot3(t1, t2))
    κ  = hcat(dot3(k, (d2l + d2r)./2),
             -dot3(k, (d1l + d1r)./2))

    # Twists
    d1_transport = ptransport(d1l, t1, t2)
    cosτ = dot3(d1_transport, d1r)
    sinτ = dot3(d1_transport, d2r)
    τ = -1.0.*atan.(sinτ, cosτ)

    l, κ, τ
end
"""
        rod_update(x0, d0, Δx, Δθ)

Update the centerline and directors of a rod.

The vertices `x0` are translated `Δx`. The directors `d0` are parallel
transported to the new centerline, then rotated by `Δθ`

# Arguments
- `x0`: centerline vertices, size (3 x N)
- `d0`: directors, size (3 x N-1). Unit length and orthogonal to `edges(x0)`
- `Δx`: displacements, size (3 x N)
- `Δ0`: rotations, size (1 x N-1).

# Returns
- `x`, `d`: new centerline and directors.

"""
function rod_update(x0, d0, Δx, Δθ)

    x = x0 + Δx
    t0 = tangents(x0)
    t = tangents(x)
    d = ptransport(d0, t0, t)
    d = rotate_orthogonal_unit(d, t, Δθ)

    x, d
end

"""
        elastic_energy(l0,κ0,τ0, l,κ,τ, k,B,β)

Elastic energy, calculated from stretching, bending, and twisting strains.

# Arguments
- `l0`, `κ0`, `τ0`: reference lengths, curvatures, and twists
- `l`, `κ`, `τ`: current lengths, curvatures, and twists
- `k`, `B`, `β`: stretching, bending, and twisting stiffness.

Lengths, curvatures, and twists should be (1 x N-1), (2 x N-2), and (1 x N-2),
respectively.

 `k` and `β` may be scalar or similar to `l` and `τ`. The bending stiffness
 `B` may constant isotropic (scalar), constant anisotropic (2 x 1),
 or non-constant anisotropic (2 x N-2).
"""
function elastic_energy(l0,κ0,τ0, l,κ,τ, k,B,β)

    # Reference voronoi lengths
    vl0 = [l0[1]; l0[2:end-1]/2] + [l0[2:end-1]/2; l0[end]]

    # Energies
    Es = 1/2*sum(k./l0.*(l - l0).^2)
    Et = 1/2*sum(β./vl0.*(τ - τ0).^2)
    Eb = 1/2*sum(B./vl0.*(κ - κ0).^2)
    Es + Et + Eb
end

"""
        gravitational_energy(x, m, g)

Gravitational energy, calculated from centerline.

# Arguments
- `x`: centerline vertices, 3 x N
- `m`: edges masses, 1 x N-1. Typically `ρ*A*l0`
- `g`: magnitude of gravitational acceleration (positive).

Gravity is assumed to act in the negative `z` direction.
"""
function gravitational_energy(x, m, g)
    nv = size(x,1)
    z = (x[1:nv-1, 3] .+ x[2:nv,3])
    sum(m.*z)*g/2
end

# ------------------------------------------------------------------------------
# Convenience type overloads
# ------------------------------------------------------------------------------
edges(r::basic_rod) = edges(r.x)
tangents(r::basic_rod) = tangents(r.x)

function full_kinematics(r::basic_rod)
    l, κ, τ = full_kinematics(r.x, r.d)
    rod_strains(l, κ, τ)
end

function rod_update(r::basic_rod, dr::rod_delta)
    x, d = rod_update(r.x, r.d, dr.Δx, dr.Δθ)
    basic_rod(x,d)
end

function elastic_energy(ref_strains::rod_strains, strains::rod_strains, rod_props::rod_properties)
    l0, κ0, τ0 = ref_strains.l, ref_strains.κ, ref_strains.τ
    l, κ, τ = strains.l, strains.κ, strains.τ
    k, B, β = rod_props.k, rod_props.B, rod_props.β
    elastic_energy(l0,κ0,τ0, l,κ,τ, k,B,β)
end


text:
include("day7_solution.jl")
using Test

test_string_1 = "light red bags contain 1 bright white bag, 2 muted yellow bags."
expected_output_1 = [("light red", "bright white", 1), ("light red", "muted yellow", 2)]
@test edges_from_rule(test_string_1) == expected_output_1

test_string_2 = "faded blue bags contain no other bags."
expected_output_2 = [("faded blue", nothing, 0)]
@test edges_from_rule(test_string_2) == expected_output_2


text:
# TODO: store types other than QBool

struct QBuffer
  bits::Bits
end

QBuffer(n::Integer) = QBuffer(Bits{n}())

Base.show(io::IO, buf::QBuffer) = print(io, "QBuffer(", length(buf), ")")

Base.length(buf::QBuffer) = nbits(buf.bits)
Base.lastindex(buf::QBuffer) = length(buf)

function Base.getindex(buf::QBuffer, i::Integer)
  @assert 1 ≤ i ≤ length(buf)
  QBool(Bit(buf.bits.offset+i-1, buf.bits.state))
end

function Base.getindex(buf::QBuffer, i::UnitRange)
  @assert 1 ≤ i.start ≤ length(buf)
  @assert 1 ≤ i.stop  ≤ length(buf)
  QBuffer(Bits{length(i)}(buf.bits.offset+i.start-1, buf.bits.state))
end


text:
"""
    advection_pde!(duhat,uhat,p,t)

RHS for the advection equation ``u_t = - u_x`` for numerical integration in Fourier space.

"""
function advection_pde!(duhat,uhat,p,t)
    N, k, dL, nu, alpha = p;
    uhat[Int(N/2)+1] = 0; # Set the most negative mode to zero to prevent an asymmetry
    duhat .= -alpha*im*k.*uhat;
end

"""
    advection_diffusion_pde!(duhat,uhat,p,t)

RHS for the advection-diffusion equation \$u_t = - u_x + ν u_{xx}\$ for numerical integration in Fourier space where ν is the viscosity.

"""
function advection_diffusion_pde!(duhat,uhat,p,t)
    N, k, dL, nu, alpha = p;
    uhat[Int(N/2)+1] = 0; # Set the most negative mode to zero to prevent an asymmetry
    duhat .= -alpha*im*k.*uhat + nu*(im*k).^2 .*uhat;
end

"""
    viscous_burgers_pde!(duhat,uhat,p,t)

RHS for the viscous Burgers equation \$u_t = - u u_x + ν u_{xx}\$ for numerical integration in Fourier space where ν is the viscosity.

"""
function viscous_burgers_pde!(duhat,uhat,p,t)
    N, k, dL, nu = p;
    alpha = 1.0;
    uhat[Int(N/2)+1] = 0.0; # Set the most negative mode to zero to prevent an asymmetry
    uhat_nonlinear = quadratic_nonlinear(uhat,N,dL,alpha);
    duhat .= -nu*k.^2 .*uhat .- uhat_nonlinear;
end

"""
    inviscid_burgers_pde!(duhat,uhat,p,t)

RHS for the inviscid Burgers equation \$u_t = - u u_x\$ for numerical integration in Fourier space.

"""
function inviscid_burgers_pde!(duhat,uhat,p,t)
    N, k, dL, nu = p;
    alpha = 1.0;
    uhat[Int(N/2)+1] = 0; # Set the most negative mode to zero to prevent an asymmetry
    uhat_nonlinear = quadratic_nonlinear(uhat,N,dL,alpha);
    duhat .= -uhat_nonlinear;
end

"""
    quadratic_nonlinear!(uhat,N,dL,alpha)

Compute the convolution sum \$\\frac{ik}{2}\\sum_{p+q=k} u_p u_q\$ resulting from the quadratic nonlinearity of Burgers equation \$u u_x\$ in Fourier space. Convolution sum is padded with the 3/2 rule for dealiasing.

"""
function quadratic_nonlinear(uhat,N,dL,alpha)
    M = Int((3/2)*N); # For dealiasing -> For MZ models must be 3*N to account for unresolved modes + dealiasing
    k_M = reduce(vcat,(2*π/dL)*[0:M/2-1 -M/2:-1]); # Wavenumbers for convolution sum with padding
    u_M = zeros(Complex{Float64},Int(M));
    u_M[1:Int(N/2)] = uhat[1:Int(N/2)];
    u_M[end-Int(N/2)+2:end] = uhat[Int(N/2)+2:end];
    u_M_sum = fft_norm(ifft_norm(u_M).*ifft_norm(u_M));
    u_M_convolution = alpha*im.*k_M/2.0 .*u_M_sum;
    uhat_nonlinear  = zeros(Complex{Float64},Int(N)); # Extract N modes from convolution
    uhat_nonlinear[1:Int(N/2)] = u_M_convolution[1:Int(N/2)];
    uhat_nonlinear[end-Int(N/2)+2:end] = u_M_convolution[end-Int(N/2)+2:end];
    return uhat_nonlinear
end

"""
    generate_fourier_solution(L1,L2,tspan,N,initial_condition,pde_function;dt=1e-3,nu=0.1,rtol=1e-10,atol=1e-14)

Generate the solution for a given `pde_function` and `initial_condition` on a periodic domain using a `N` mode Fourier expansion.

"""
function generate_fourier_solution(L1,L2,tspan,N,initial_condition,pde_function;dt=1e-3,nu=0.1,rtol=1e-10,atol=1e-14,alpha=1.0)
    # Transform random initial condition to Fourier domain
    uhat0 = fft_norm(initial_condition);
    dL = abs(L2-L1);

    # Generate Fourier Galerkin solution for N
    k = reduce(vcat,(2*π/dL)*[0:N/2-1 -N/2:-1]);
    p = [N,k,dL,nu,alpha]
    t_length = Int(round(tspan[2]/dt)+1);

    # Solve the system of ODEs in Fourier domain
    prob = ODEProblem(pde_function,uhat0,tspan,p);
    sol = solve(prob,DP5(),reltol=rtol,abstol=atol,saveat = dt)

    u_sol = zeros(t_length,N);
    for j in 1:size(sol.t,1) # Reshape output and plot
        u_sol[j,:] = real.(ifft_norm(sol.u[j]));
    end
    return u_sol, sol.u, k
end

"""
    central_difference(u_j,u_jpos,u_jneg,mu)

Compute the second order central difference for the viscous term of the viscous Burgers equation \$u_t = - u u_x + ν u_{xx}\$. `mu` is equal to \$ \\frac{\\nu}{\\Delta x^2}\$.

"""
function central_difference(u_j,u_jpos,u_jneg,mu)
    ux = zeros(size(u_j,2));
    ux = mu.*(u_jpos-2*u_j+u_jneg);
    return ux
end

"""
    minmod(x,y)

"""
function minmod(x,y)
    return sign(x)*max(0,min(abs(x),y*sign(x)));
end

"""
    ub(ulv,urv)

"""
function ub(ulv,urv)
    if ulv != urv
        return ((1/2)*urv^2 - (1/2)*ulv^2) / (urv - ulv); # ̄u_{j+/-1/2}
    else
        return ulv
    end
end

"""
    fl(ulv,urv,ubv)

"""
function fl(ulv,urv,ubv)
    return (1/2).*((1/2).*ulv.^2 .+ (1/2).*urv.^2 .- abs.(ubv).*(urv .- ulv)); # f_{j+/-1/2}
# fln(ulv,urv,ubv) = (1/2)*(flux.(urv) - flux.(ulv) - abs.(ubv).*(ulv - urv)); # f_{j+/-1/2}
end

"""
    ulpl(ujp,uj,ujn,kappa,omega)

"""
function ulpl(ujp,uj,ujn,kappa,omega)
    return uj .+ ((1-kappa)/4).*minmod.((uj .- ujn),omega.*(ujp .- uj)) .+ ((1+kappa)/4).*minmod.((ujp .- uj),omega.*(uj .- ujn)); # u_{j+1/2}^L
end

"""
    urpl(ujpp,ujp,uj,kappa,omega)

"""
function urpl(ujpp,ujp,uj,kappa,omega)
    return ujp  .- ((1+kappa)/4).*minmod.((ujp .- uj),omega.*(ujpp .- ujp)) .- ((1-kappa)/4).*minmod.((ujpp .- ujp),omega.*(ujp .- uj)); # u_{j+1/2}^R
end

"""
    ulnl(uj,ujn,ujnn,kappa,omega)

"""
function ulnl(uj,ujn,ujnn,kappa,omega)
    return ujn .+ ((1-kappa)/4).*minmod.((ujn .- ujnn),omega.*(uj .- ujn)) .+ ((1+kappa)/4).*minmod.((uj .- ujn),omega.*(ujn .- ujnn)); # u_{j-1/2}^L
end

"""
    urnl(ujp,uj,ujn,kappa,omega)

"""
function urnl(ujp,uj,ujn,kappa,omega)
    return uj .- ((1+kappa)/4).*minmod.((uj .- ujn),omega.*(ujp .- uj)) .- ((1-kappa)/4).*minmod.((ujp .- uj),omega.*(uj .- ujn)); # u_{j-1/2}^R
end

"""
    muscl_minmod_RHS!(du,u,p,t)

"""
function muscl_minmod_RHS!(du,u,p,t)
    dx, kappa, omega, nu = p;

    ulpvp = zeros(size(u,1));
    urpvp = zeros(size(u,1));
    ulnvp = zeros(size(u,1));
    urnvp = zeros(size(u,1));

    # ulp(ujp,uj,ujn,omega)
    ulpvp[2:end-1] .= ulpl(u[3:end],u[2:end-1],u[1:end-2],kappa,omega);
    ulpvp[1] = ulpl(u[2],u[1],u[end],kappa,omega);
    ulpvp[end] = ulpl(u[1],u[end],u[end-1],kappa,omega);

    # urp(ujpp,ujp,uj,omega)
    urpvp[1:end-2] .= urpl(u[3:end],u[2:end-1],u[1:end-2],kappa,omega);
    urpvp[end-1] = urpl(u[1],u[end],u[end-1],kappa,omega);
    urpvp[end] = urpl(u[2],u[1],u[end],kappa,omega);

    ubpp = ub.(ulpvp,urpvp);
    fpp = fl(ulpvp,urpvp,ubpp);

    # uln(uj,ujn,ujnn,omega)
    ulnvp[3:end] .= ulnl(u[3:end],u[2:end-1],u[1:end-2],kappa,omega);
    ulnvp[2] = ulnl(u[2],u[1],u[end],kappa,omega);
    ulnvp[1] = ulnl(u[1],u[end],u[end-1],kappa,omega);

    # urn(ujp,uj,ujn,omega)
    urnvp[2:end-1] .= urnl(u[3:end],u[2:end-1],u[1:end-2],kappa,omega);
    urnvp[1] = urnl(u[2],u[1],u[end],kappa,omega);
    urnvp[end] = urnl(u[1],u[end],u[end-1],kappa,omega);

    ubnp = ub.(ulnvp,urnvp);
    fnp = fl(ulnvp,urnvp,ubnp);

    du .= -(1/dx).*(fpp .- fnp);
end

"""
    muscl_minmod_viscous_RHS!(du,u,p,t)

"""
function muscl_minmod_viscous_RHS!(du,u,p,t)
    dx, kappa, omega, nu = p;

    mu = nu/(dx)^2;

    ulpvp = zeros(size(u,1));
    urpvp = zeros(size(u,1));
    ulnvp = zeros(size(u,1));
    urnvp = zeros(size(u,1));
    ux_viscous = zeros(size(u,1));

    # ulp(ujp,uj,ujn,omega)
    ulpvp[2:end-1] .= ulpl(u[3:end],u[2:end-1],u[1:end-2],kappa,omega);
    ulpvp[1] = ulpl(u[2],u[1],u[end],kappa,omega);
    ulpvp[end] = ulpl(u[1],u[end],u[end-1],kappa,omega);

    # urp(ujpp,ujp,uj,omega)
    urpvp[1:end-2] .= urpl(u[3:end],u[2:end-1],u[1:end-2],kappa,omega);
    urpvp[end-1] = urpl(u[1],u[end],u[end-1],kappa,omega);
    urpvp[end] = urpl(u[2],u[1],u[end],kappa,omega);

    ubpp = ub.(ulpvp,urpvp);
    fpp = fl(ulpvp,urpvp,ubpp);

    # uln(uj,ujn,ujnn,omega)
    ulnvp[3:end] .= ulnl(u[3:end],u[2:end-1],u[1:end-2],kappa,omega);
    ulnvp[2] = ulnl(u[2],u[1],u[end],kappa,omega);
    ulnvp[1] = ulnl(u[1],u[end],u[end-1],kappa,omega);

    # urn(ujp,uj,ujn,omega)
    urnvp[2:end-1] .= urnl(u[3:end],u[2:end-1],u[1:end-2],kappa,omega);
    urnvp[1] = urnl(u[2],u[1],u[end],kappa,omega);
    urnvp[end] = urnl(u[1],u[end],u[end-1],kappa,omega);

    ux_viscous[1] = central_difference(u[1],u[2],u[end],mu);
    ux_viscous[2:end-1] = central_difference(u[2:end-1],u[3:end],u[1:end-2],mu);
    ux_viscous[end] = central_difference(u[end],u[1],u[end-1],mu);

    ubnp = ub.(ulnvp,urnvp);
    fnp = fl(ulnvp,urnvp,ubnp);

    du .= -(1/dx).*(fpp .- fnp) .+ ux_viscous;
end

"""
    generate_muscl_minmod_solution(L1,L2,t_end,N,u0,pde_function_handle;dt=1e-4,kappa=-1)

"""
function generate_muscl_minmod_solution(L1,L2,t_end,N,u0,pde_function_handle;dt=1e-4,kappa=-1,nu=0.1)

    omega = ((3-kappa)/(1-kappa))

    dL = abs(L2-L1);
    # Set up periodic domain
    j = reduce(vcat,[0:1:N-1]);
    x = (dL.*j)./N;
    dx = x[2]-x[1];
    t = reduce(vcat,[0:dt:t_end]);

    p = [dx,kappa,omega,nu]
    t_span = (0,t_end);

    prob = ODEProblem(pde_function_handle,u0,t_span,p);
    sol = solve(prob,BS3(),reltol=1e-6,abstol=1e-8,saveat = dt);

    u_sol = zeros(size(sol.t,1),size(x,1));
    for i in 1:size(sol.t,1)
        u_sol[i,:] = sol.u[i];
    end

    return u_sol

end

"""
    generate_muscl_reduced(L1,L2,t_end,dt,M,N,ic_func,pde_function_handle;kappa=-1,nu=0.1)

"""
# function generate_muscl_reduced(L1,L2,t_end,dt,M,N,ic_func,pde_function_handle;kappa=-1,nu=0.1)
function generate_muscl_reduced(L1,L2,t_end,dt,M,N,ic,pde_function_handle;kappa=-1,nu=0.1)
    dL = abs(L2 - L1);
    j = reduce(vcat,[0:1:N-1]);
    x = (dL.*j)./N;
    j_full = reduce(vcat,[0:1:M-1]);
    x_full = (dL.*j_full)./M;
    # ic = ic_func.(x_full);
    u_full = generate_muscl_minmod_solution(L1,L2,t_end,M,ic,pde_function_handle;dt=dt,kappa=kappa,nu=nu)
    u_reduced = solution_spatial_sampling(x,x_full,u_full);
    return u_full, u_reduced
end

"""
    get_1D_energy_fft(u_solution)

Compute the energy in the Fourier domain using the scaling of \$ \\frac{1}{N} \$. Note: this does not include the 2π multiplier found in Parseval's identity for Fourier series and computes \$ \\frac{1}{2} \\sum \\vert \\hat{u}_k \\vert^2 \$.

"""
function get_1D_energy_fft(u_solution) # To Do: Add and extraction of specific mode step
    energy = zeros(size(u_solution,1))
    for i in 1:size(u_solution,1)
        u_hat = fft_norm(u_solution[i,:]);
        energy[i] = (1/2)*real((u_hat'*u_hat));
    end
    return energy
end

"""
    get_1D_energy_custom(basis,u_solution,L1,L2,weights;multiplier=1/(4*pi))

Compute the energy in the custom basis domain: \$ \\frac{1}{2} \\sum \\vert \\a_k \\vert^2 \$. Multiplier defaults to \$ \\frac{1}{4\\pi} \$ to match the Fourier calculation.

"""
function get_1D_energy_custom(basis,u_solution,nodes,weights;multiplier=1/(4*pi))
    energy = zeros(size(u_solution,1));
    for i in 1:size(u_solution,1)
        a_hat = expansion_coefficients(basis,u_solution[i,:],nodes,weights);
        energy[i] = get_1D_energy_custom_coefficients(a_hat,zeta=multiplier);
    end
    return energy
end

"""
    get_1D_energy_custom_coefficients(u_coefficients;zeta=1/(4*pi))

Compute the energy in the custom basis domain: \$ \\zeta \\sum \\vert \\a_k \\vert^2 \$. \$ \\zeta \$ defaults to \$ \\frac{1}{4\\pi} \$ to match the Fourier calculation.

"""
function get_1D_energy_custom_coefficients(u_coefficients;zeta=1/(4*pi))
        # energy[i] = (1/2)*(a_hat'*a_hat);
        energy = zeta*(u_coefficients'*u_coefficients);
    return energy
end

"""
    mode_extractor(uhat,N)

"""
function mode_extractor(uhat,N)
    uhat_N = zeros(typeof(uhat[1]),N);
    uhat_N[1:Int(N/2)] = uhat[1:Int(N/2)];
    uhat_N[end-Int(N/2)+2:end] = uhat[end-Int(N/2)+2:end];
    return uhat_N
end

"""
    get_1D_energy_upwind(u_solution_full,u_solution,N)

"""
function get_1D_energy_upwind(u_solution_full,u_solution,N)
    energy = zeros(size(u_solution,1))
    for i in 1:size(u_solution,1)
        u_upwind_fft = mode_extractor(fft_norm(u_solution_full[i,:]),N);
        energy[i] = (1/2)*real(u_upwind_fft'*u_upwind_fft);
    end
    return energy
end

"""
    spectral_approximation_fourier(x_locations,k,coefficients)

"""
function spectral_approximation_fourier(x_locations,k,coefficients)
    approximation = zeros(size(x_locations,1));
    for i = 1:length(k)
        approximation += coefficients[i]*exp.(im*k[i]*(x_locations));
    end
    return approximation
end


text:
using Pkg; Pkg.add("Gadfly");

using Documenter, FocusedBlindDecon
using Gadfly

makedocs(
	 format = Documenter.HTML(
		  prettyurls = get(ENV, "CI", nothing) == "true"),
   sitename = "Focused Blind Deconvolution",
   pages = ["Home" => "index.md",
	    "Library" => "library.md",
	    "Tutorial I" => "tut1.md"
	    #"SeisForwExpt" => Any[
			#	  "LBSD" => "lsbd.md",
#				  "Basic usage" => "Fdtd/reuse_expt.md",
#				  "Generate snaps" => "Fdtd/create_snaps.md",
#				  ],
	    ]
#	    "Seismic Born Modeling" => "FWI/born_map.md",
#	    "Seismic Full Waveform Inversion" => "FWI/gradient_accuracy.md",
#    modules = []
)

deploydocs(
    repo   = "github.com/pawbz/FocusedBlindDecon.jl.git",
)


text:
"""
    clique_percolation(g, k=3)

Community detection using the clique percolation algorithm. Communities are potentionally overlapping.
Return a vector of vectors `c` such that `c[i]` is the set of vertices in community `i`.
The parameter `k` defines the size of the clique to use in percolation.

### References
- [Palla G, Derenyi I, Farkas I J, et al.] (https://www.nature.com/articles/nature03607)

# Examples
```jldoctest
julia> using LightGraphs

julia> clique_percolation(clique_graph(3, 2))
2-element Array{BitSet,1}:
 BitSet([4, 5, 6])
 BitSet([1, 2, 3])

julia> clique_percolation(clique_graph(3, 2), k=2)
1-element Array{BitSet,1}:
 BitSet([1, 2, 3, 4, 5, 6])

julia> clique_percolation(clique_graph(3, 2), k=4)
0-element Array{BitSet,1}
```
"""
function clique_percolation end

@traitfn function clique_percolation(g::::(!IsDirected); k = 3)
  kcliques = filter(x->length(x)>=k, maximal_cliques(g))
  nc = length(kcliques)
  # graph with nodes represent k-cliques
  h = Graph(nc)
  # vector for counting common nodes between two cliques efficiently
  x = falses(nv(g))
  for i = 1:nc
    x[kcliques[i]] .= true
    for j = i+1:nc
        sum(x[kcliques[j]]) >= k-1 && add_edge!(h, i, j)
    end
    # reset status
    x[kcliques[i]] .= false
  end
  components = connected_components(h)
  communities = [BitSet() for i=1:length(components)]
  for (i,component) in enumerate(components)
    push!(communities[i], vcat(kcliques[component]...)...)
  end
  return communities
end


text:
module ldtiSubTests

using LinearAlgebra
using Subsystems
using Test

@testset "ldti Subsystem" begin
    A = Matrix{Float64}(I, 2, 2);
    B = reshape([1., 0.], :, 1)
    C = I;

    @testset "defines correct system" begin
        s = LdtiSubsystem(1, A, B, C)
        @test s.A == A
        @test s.C == I
        @test s.B == B
    end

    @testset "implements correct conservation" begin
        s = LdtiSubsystem(1, A, B, C)
        s1 = LdtiSubsystem(2, A, B, C)
        k = [0.1 0; 0 0]

        addNeighbour(s, s1, k, true)
        @test s.A == A - k

        removeNeighbour(s, 2, true)
        @test s.A == A
        addNeighbour(s, s1, k, true)
        setNeighbour(s, s1, 2, 0.2*k, true)
        @test s.A == A - 0.2*k

        s = LdtiSubsystem(1, A, B, C)
        s2 = LdtiSubsystem(3, A, B, C)

        addNeighbour(s, s1, k, true)
        addNeighbour(s, s2, k, true)
        @test s.A == A - 2*k
    end

    @testset "updateE" begin
        s = LdtiSubsystem(1, A, B, C)
        s1 = LdtiSubsystem(2, A, B, C)
        k = [0.1 0; 0 0]

        addNeighbour(s, s1, k)
        @test s.E == k

        removeNeighbour(s, 2)
        @test s.E == zeros(s.nx, 0)
        addNeighbour(s, s1, k)
        setNeighbour(s, s1, 2, 0.2*k)
        @test s.E == 0.2*k

        s = LdtiSubsystem(1, A, B, C)
        s2 = LdtiSubsystem(3, A, B, C)

        addNeighbour(s, s1, k)
        addNeighbour(s, s2, k)
        @test s.E == repeat(k, 1, 2)
    end

    @testset "getNeighbourStates" begin
        @testset "gets the correct state" begin
            s1 = LdtiSubsystem(1, A, B, C)
            s2 = LdtiSubsystem(2, A, B, C)
            s3 = LdtiSubsystem(3, A, B, C)
            Aij = rand(2,2)
            addNeighbour(s1, s2, Aij)
            addNeighbour(s1, s3, Aij)
            x2 = rand(2)
            x3 = rand(2)
            setInitialState(s2, x2)
            setInitialState(s3, x3)
            x1bar = getNeighbourStates(s1)
            @test all(x1bar .≈ vcat(x2, x3))
        end
        @testset "gets correct zero on no neighbours" begin
            s1 = LdtiSubsystem(1, A, B, C)
            x1bar = getNeighbourStates(s1)
            @test x1bar == zero(s1.x)
        end
    end

    @testset "getOutboundWeights" begin
        s1 = LdtiSubsystem(1, A, B, C)
        s2 = LdtiSubsystem(2, A, B, C)
        s3 = LdtiSubsystem(3, A, B, C)
        A1x = rand(2,2)
        A31 = rand(2,2)
        A21 = rand(2,2)
        addNeighbour(s1, s2, A1x)
        addNeighbour(s1, s3, A1x)
        addNeighbour(s2, s1, A21)
        addNeighbour(s3, s1, A31)

        @test getOutboundWeights(s1) == [A21, A31]
    end

    @testset "updateState" begin
        @testset "implements correct state equation" begin
            s = LdtiSubsystem(1, A, B, C)
            updateState(s, 1.0)
            @test s.xnext == [1.,0.]
            @test s.x == [0., 0.]
        end
    end
end

end