text:
import Base: isempty

export AffineMap,
       an_element,
       isempty,
       isbounded,
       ∈,
       vertices_list,
       constraints_list,
       linear_map

"""
    AffineMap{N<:Real, S<:LazySet{N}, NM, MAT<:AbstractMatrix{NM},
              VN<:AbstractVector{NM}} <: LazySet{N}

Type that represents an affine transformation ``M⋅X ⊕ v`` of a convex set ``X``.

### Fields

- `M` -- matrix/linear map
- `X` -- convex set
- `v` -- translation vector

### Notes

An affine map is the composition of a linear map and a translation. This type is
parametric in the coefficients of the linear map, `NM`, which may be different from
the numeric type of the wrapped set (`N`). However, the numeric type of the
translation vector should be `NM`.

### Examples

For the examples we create a ``3×2`` matrix, a two-dimensional unit square, and
a three-dimensional vector.
Then we combine them in an `AffineMap`.

```jldoctest constructors
julia> A = [1 2; 1 3; 1 4]; X = BallInf([0, 0], 1); b2 = [1, 2]; b3 = [1, 2, 3];

julia> AffineMap(A, X, b3)
AffineMap{Int64,BallInf{Int64},Int64,Array{Int64,2},Array{Int64,1}}([1 2; 1 3; 1 4], BallInf{Int64}([0, 0], 1), [1, 2, 3])
```

For convenience, `A` does not need to be a matrix but we also allow to use
`UniformScaling`s resp. scalars (interpreted as a scaling, i.e., a scaled
identity matrix).
Scaling by ``1`` is ignored and simplified to a pure `Translation`.

```jldoctest constructors
julia> using LinearAlgebra

julia> am = AffineMap(2I, X, b2)
AffineMap{Int64,BallInf{Int64},Int64,Diagonal{Int64,Array{Int64,1}},Array{Int64,1}}([2 0; 0 2], BallInf{Int64}([0, 0], 1), [1, 2])

julia> AffineMap(2, X, b2) == am
true

julia> AffineMap(1, X, b2)
Translation{Int64,Array{Int64,1},BallInf{Int64}}(BallInf{Int64}([0, 0], 1), [1, 2])
```

Applying a linear map to an `AffineMap` object combines the two maps into a new
`AffineMap` instance.
Again we can make use of the conversion for convenience.

```jldoctest constructors
julia> B = [2 0; 0 2]; am2 = B * am
AffineMap{Int64,BallInf{Int64},Int64,Array{Int64,2},Array{Int64,1}}([4 0; 0 4], BallInf{Int64}([0, 0], 1), [2, 4])

julia> 2 * am == am2
true
```

The application of an `AffineMap` to a `ZeroSet` or an `EmptySet` is simplified
automatically.

```jldoctest constructors
julia> AffineMap(A, ZeroSet{Int}(2), b3)
Singleton{Int64,Array{Int64,1}}([1, 2, 3])

julia> AffineMap(A, EmptySet{Int}(), b3)
EmptySet{Int64}()
```
"""
struct AffineMap{N<:Real, S<:LazySet{N}, NM, MAT<:AbstractMatrix{NM},
                 VN<:AbstractVector{NM}} <: LazySet{N}
    M::MAT
    X::S
    v::VN

    # default constructor with dimension match check
    function AffineMap(M::MAT, X::S, v::VN) where {N<:Real, S<:LazySet{N}, NM,
                                                   MAT<:AbstractMatrix{NM},
                                                   VN<:AbstractVector{NM}}

        @assert dim(X) == size(M, 2) "a matrix of size $(size(M)) cannot be " *
            "applied to a set of dimension $(dim(X))"

        @assert size(M, 1) == length(v) "a map with output dimension " *
            "$(size(M, 1)) is incompatible with a translation vector of " *
            "dimension $(length(v))"

        return new{N, S, NM, MAT, VN}(M, X, v)
    end
end

isoperationtype(::Type{<:AffineMap}) = true

# convenience constructor from a UniformScaling
function AffineMap(M::UniformScaling, X::LazySet, v::AbstractVector)
    return AffineMap(M.λ, X, v)
end

# convenience constructor from a scalar
function AffineMap(α::N, X::LazySet, v::AbstractVector) where {N<:Real}
    if α == one(N)
        return Translation(X, v)
    end
    return AffineMap(Diagonal(fill(α, length(v))), X, v)
end

# simplification for a LinearMap
for MAP in (:AbstractMatrix, :Real)
    @eval begin
        function LinearMap(map::$MAP, am::AffineMap)
            return AffineMap(map * am.M, am.X, map * am.v)
        end
    end
end

# ZeroSet is "almost absorbing" for the linear map (only the dimension changes)
# such that only the translation vector remains
function AffineMap(M::AbstractMatrix{N}, Z::ZeroSet{N}, v::AbstractVector{N}
                  ) where {N<:Real}
    @assert dim(Z) == size(M, 2) "a matrix of size $(size(M)) cannot be " *
        "applied to a set of dimension $(dim(Z))"

    @assert size(M, 1) == length(v) "a map with output dimension " *
        "$(size(M, 1)) is incompatible with a translation vector of " *
        "dimension $(length(v))"

    return Singleton(v)
end

# EmptySet is absorbing for AffineMap
function AffineMap(M::AbstractMatrix{N}, ∅::EmptySet{N}, v::AbstractVector{N}
                  ) where {N<:Real}
    return ∅
end

# ============================
# LazySet interface functions
# ============================

"""
    dim(am::AffineMap)::Int

Return the dimension of an affine map.

### Input

- `am` -- affine map

### Output

The dimension of an affine map.
"""
function dim(am::AffineMap)::Int
    return length(am.v)
end

"""
    σ(d::AbstractVector{N}, am::AffineMap{N}) where {N<:Real}

Return the support vector of an affine map.

### Input

- `d`  -- direction
- `am` -- affine map

### Output

The support vector in the given direction.
"""
function σ(d::AbstractVector{N}, am::AffineMap{N}) where {N<:Real}
    return am.M * σ(_At_mul_B(am.M, d), am.X) + am.v
end

"""
    ρ(d::AbstractVector{N}, am::AffineMap{N}) where {N<:Real}

Return the support function of an affine map.

### Input

- `d`  -- direction
- `am` -- affine map

### Output

The support function in the given direction.
"""
function ρ(d::AbstractVector{N}, am::AffineMap{N}) where {N<:Real}
    return ρ(_At_mul_B(am.M, d), am.X) + dot(d, am.v)
end

"""
    an_element(am::AffineMap)

Return some element of an affine map.

### Input

- `am` -- affine map

### Output

An element of the affine map. It relies on the `an_element` function of the
wrapped set.
"""
function an_element(am::AffineMap)
    return am.M * an_element(am.X) + am.v
end

"""
    isempty(am::AffineMap)::Bool

Return whether an affine map is empty or not.

### Input

- `am` -- affine map

### Output

`true` iff the wrapped set is empty and the affine vector is empty.
"""
function isempty(am::AffineMap)::Bool
    return isempty(am.X)
end

"""
    isbounded(am::AffineMap; cond_tol::Number=DEFAULT_COND_TOL)::Bool

Determine whether an affine map is bounded.

### Input

- `am`       -- affine map
- `cond_tol` -- (optional) tolerance of matrix condition (used to check whether
                the matrix is invertible)

### Output

`true` iff the affine map is bounded.

### Algorithm

We first check if the matrix is zero or the wrapped set is bounded.
If not, we perform a sufficient check whether the matrix is invertible.
If the matrix is invertible, then the map being bounded is equivalent to the
wrapped set being bounded, and hence the map is unbounded.
Otherwise, we check boundedness via [`isbounded_unit_dimensions`](@ref).
"""
function isbounded(am::AffineMap; cond_tol::Number=DEFAULT_COND_TOL)::Bool
    if iszero(am.M) || isbounded(am.X)
        return true
    end
    if isinvertible(am.M; cond_tol=cond_tol)
        return false
    end
    return isbounded_unit_dimensions(am)
end

"""
    ∈(x::AbstractVector{N}, am::AffineMap{N})::Bool where {N<:Real}

Check whether a given point is contained in the affine map of a convex set.

### Input

- `x`  -- point/vector
- `am` -- affine map of a convex set

### Output

`true` iff ``x ∈ am``.

### Algorithm

Note that ``x ∈ M⋅S ⊕ v`` iff ``M^{-1}⋅(x - v) ∈ S``.
This implementation does not explicitly invert the matrix, which is why it also
works for non-square matrices.

### Examples

```jldoctest
julia> am = AffineMap([2.0 0.0; 0.0 1.0], BallInf([1., 1.], 1.), [-1.0, -1.0]);

julia> [5.0, 1.0] ∈ am
false

julia> [3.0, 1.0] ∈ am
true
```

An example with a non-square matrix:

```jldoctest
julia> B = BallInf(zeros(4), 1.);

julia> M = [1. 0 0 0; 0 1 0 0]/2;

julia> [0.5, 0.5] ∈ M*B
true
```
"""
function ∈(x::AbstractVector{N}, am::AffineMap{N})::Bool where {N<:Real}
    return ∈(am.M \ (x - am.v), am.X)
end

"""
    vertices_list(am::AffineMap{N};
                  [apply_convex_hull]::Bool)::Vector{Vector{N}} where {N<:Real}

Return the list of vertices of a (polyhedral) affine map.

### Input

- `am`                -- affine map
- `apply_convex_hull` -- (optional, default: `true`) if `true`, apply the convex
                         hull operation to the list of vertices transformed by the
                         affine map

### Output

A list of vertices.

### Algorithm

This implementation computes all vertices of `X`, then transforms them through
the affine map, i.e. `x ↦ M*x + v` for each vertex `x` of `X`. By default, the
convex hull operation is taken before returning this list. For dimensions three
or higher, this operation relies on the functionality through the concrete
polyhedra library `Polyhedra.jl`.

If you are not interested in taking the convex hull of the resulting vertices under
the affine map, pass `apply_convex_hull=false` as a keyword argument.

Note that we assume that the underlying set `X` is polyhedral, either concretely
or lazily, i.e. there the function `vertices_list` should be applicable.
"""
function vertices_list(am::AffineMap{N};
                       apply_convex_hull::Bool=true)::Vector{Vector{N}} where {N<:Real}

    # collect vertices list of the wrapped set
    vlist_X = vertices_list(am.X)

    # create resulting vertices list
    vlist = [am.M * x + am.v for x in vlist_X]

    return apply_convex_hull ? convex_hull!(vlist) : vlist
end

"""
    constraints_list(am::AffineMap{N}) where {N<:Real}

Return the list of constraints of a (polyhedral) affine map.

### Input

- `am` -- affine map

### Output

The list of constraints of the affine map.

### Notes

We assume that the underlying set `X` is polyhedral, i.e., offers a method
`constraints_list(X)`.

### Algorithm

Falls back to the list of constraints of the translation of a lazy linear map.
"""
function constraints_list(am::AffineMap{N}) where {N<:Real}
    return _constraints_list_translation(LinearMap(am.M, am.X), am.v)
end

"""
    linear_map(M::AbstractMatrix{N}, am::AffineMap{N}) where {N<:Real}

Return the linear map of a lazy affine map.

### Input

- `M`  -- matrix
- `am` -- affine map

### Output

A set corresponding to the linear map of the lazy affine map of a set.
"""
function linear_map(M::AbstractMatrix{N}, am::AffineMap{N}) where {N<:Real}
     return translate(linear_map(M * am.M, am.X), M * am.v)
end


text:
module CompatTests

using FillArrays, ReverseDiff, Test

@test ReverseDiff.gradient(fill(2.0, 3)) do x
  sum(abs2.(x .- Zeros(3)))
end == fill(4.0, 3)

@test ReverseDiff.gradient(fill(2.0, 3)) do x
  sum(abs2.(x .- (1:3)))
end == [2, 0, -2]

end

ext: jl


text:
module Utils

import ..Lib
export invoke_and_check
export @invoke_and_check

const MAPPING = Lib.CEnum.namemap(Lib.helics_error_types)

function snakecase_to_camelcase(s::String)
    s = replace(s, r"['`]"=>"")
    s = split(s, r"[a-z]\K(?=[A-Z][a-z])|\W|_")
    s = join(map(uppercasefirst, s))
end

function snakecase_to_camelcase(s::Symbol)
    s = String(s)
    s = snakecase_to_camelcase(s)
    s = replace(s, "Helics" => "HELICS")
    s = Symbol(s)
end

abstract type HELICSException <: Exception end

for (sym, ans) in Lib.CEnum.name_value_pairs(Lib.helics_error_types)
    sym = snakecase_to_camelcase(sym)
    eval(:(struct $sym <: HELICSException
               msg::String
           end))
    eval(:(export $sym))
end

function invoke_and_check(f::Function, args...)
    err = Lib.helicsErrorInitialize()
    ref_err = Ref(err)
    r = f(args..., ref_err)
    err = ref_err[]
    if err.error_code == 0
        return r
    elseif err.error_code == -8
        func = String(nameof(f))
        c = Lib.helics_error_types(err.error_code)
        @warn "`$func` returned warning `$(unsafe_string(err.message))`"
        return r
    else
        func = String(nameof(f))
        s = snakecase_to_camelcase(MAPPING[err.error_code])
        S = getfield(Utils, Symbol(s))
        throw(S("`$func` errored with message `$(unsafe_string(err.message))`"))
    end
end

macro invoke_and_check(expr)
    pushfirst!(expr.args, :(invoke_and_check))
    return esc(expr)
end

end


text:
using  IntegratedOptics
λ = 1.55;
Lx = 8 * λ;
Ly = 8 * λ;
sim = Simulation(λ₀ = λ ;  grid = Grid(extent = (Lx,Ly) , spacing =  (λ/15,λ/15) ))
sim(setpml!,λ)
sim(setsrc!, PlaneWave(k̂ = (0,2pi/λ), ê = (0,0,1)))
sim(setTFSF! ,1.1λ)
Si = Material(ϵᵣ = 2.0 )
Rect(x) =  abs(x[2]-Ly/2) <=  λ/2
sim(setmaterial!,Si,Rect);
sim(solve!)


text:
myParse(x) = parse.(Int, x)

# hcat(x...) converts x::Vector{Vector} to Matrix
# hcat to get columns as rows instead

bits = hcat(myParse.(collect.(readlines("input.txt")))...)

#! Part 1

mcb = x -> count(isone, x) >= length(x)/2 ? 1 : 0;
lcb = x -> count(isone, x) < length(x)/2 ? 1 : 0;

γVec = [mcb(row) for row in eachrow(bits)]
εVec = replace(γVec, 0 => 1, 1 => 0)

γ = parse(Int, string(γVec...); base=2)
ε = parse(Int, string(εVec...); base=2)

power = γ * ε

#! Part 2

nRows = length(bits[:, 1])

function findRatingRecursive(vec, measure, rowIndex=1)

    if(length(vec[1, :]) == 1)
        return vec
    end

    if(rowIndex > nRows)
        rowIndex = 1
    end

    requiredBit = measure(vec[rowIndex, :])

    filteredVec::Vector{Vector{Int}} = []

    for col in eachcol(vec)
        if col[rowIndex] == requiredBit
            push!(filteredVec, col)
        end
    end

    findRatingRecursive(hcat(filteredVec...), measure, rowIndex+1)
end

oxygenRating = parse(Int, string(findRatingRecursive(bits, mcb)...); base=2)

co2Rating = parse(Int, string(findRatingRecursive(bits, lcb)...); base=2)

lifeSupportRating = oxygenRating * co2Rating

max_stars_repo_path: Day 3/main.jl

max_issues_repo_path: Day 3/main.jl

max_forks_repo_path: Day 3/main.jl


text:
using LinearAlgebra: diagind
export heuristic_loss
"""
Pick the edge with maximum flow
"""
function eFlow(values, flows, candidates::Vector{Tuple{Node, Node}})
    edge2flows = map(candidates) do (or, and)
        count_downflow(values, flows, nothing, or, and)
    end
    (max_flow, max_edge_id) = findmax(edge2flows)
    candidates[max_edge_id], max_flow
end

"""
Pick the variable with maximum sum of mutual information
"""
function vMI(values, flows, edge, vars::Vector{Var}, train_x)
    examples_id = downflow_all(values, flows, num_examples(train_x), edge...)
    sub_matrix = train_x[examples_id, vars]
    (_, mi) = mutual_information(sub_matrix; α=1.0)
    mi[diagind(mi)] .= 0
    scores = dropdims(sum(mi, dims = 1), dims = 1)
    var = vars[argmax(scores)]
    score = maximum(scores)
    var, score
end

"""
Pick the edge randomly
"""
function eRand(candidates::Vector{Tuple{Node, Node}})
    return rand(candidates)
end

"""
Pick the variable randomly
"""
function vRand(vars::Vector{Var})
    return Var(rand(vars))
end

function heuristic_loss(circuit::LogicCircuit, train_x; pick_edge="eFlow", pick_var="vMI")
    if isweighted(train_x)
        train_x, weights = split_sample_weights(train_x)
    else
        weights = nothing
    end

    candidates, variable_scope = split_candidates(circuit)
    values, flows = satisfies_flows(circuit, train_x; weights = nothing) # Do not use samples weights here
    if pick_edge == "eFlow"
        edge, flow = eFlow(values, flows, candidates)
    elseif pick_edge == "eRand"
        edge = eRand(candidates)
    else
        error("Heuristics $pick_edge to pick edge is undefined.")
    end

    or, and = edge
    vars = Var.(collect(variable_scope[and]))

    if pick_var == "vMI"
        var, score = vMI(values, flows, edge, vars, train_x)
    elseif pick_var == "vRand"
        var = vRand(vars)
    else
        error("Heuristics $pick_var to pick variable is undefined.")
    end

    return (or, and), var
end


text:
#Una función que verifica que haya tantos registros como se espera en el archivo:
function check_linenumber(archivo)

    #Usamos wc para contar el número de líneas. Guardamos la salida en un string:
    count_lines = read(`wc -l $archivo`, String)

    #Capturamos el resultado y operamos en él.
    number_lines_string = match(r"^(\d+) ", count_lines).captures[1]
    number_lines = parse(Int64, number_lines_string)
    number_registries = number_lines - 1 #Tenemos tantos registros como líneas menos uno (por el header)

    #Cargamos la última línea del archivo:
    final_line = read(`tail -n 1 $archivo`, String)

    #Capturamos el último número de caso del registro:
    expected_number_registries_string = match(r"^(\d+),", final_line).captures[1]
    expected_number_registries = parse(Int64, expected_number_registries_string)

    if number_registries == expected_number_registries

        return nothing

    else

        error("Mismatch in the number of registries in $archivo")
    end
end


text:
# # How to write transducers
#
# `Transducers` don't export public interface for implementing
# transducers (and reducible collections).  Let's import some handy
# ones:

using Transducers
using Transducers: Transducer, R_, next, inner, xform

# ## Stateless transducer

# Let's write manually what `Filter(x -> x isa Int) |> Map(x -> x +
# 1)` would do:

struct AddOneIfInt <: Transducer end

function Transducers.next(rf::R_{AddOneIfInt}, result, input)
    if input isa Int
#+
# Output `input + 1` is passed to the "inner" reducing step:
        next(inner(rf), result, input + 1)
    else
#+
# Filtering out is done by "doing nothing"; return `result`-so-far
# as-is:
        result
    end
end

# Ineed, for integer collection, it increments input by one:

addone_out1 = begin  # hide
collect(AddOneIfInt(), 1:5)
end  # hide

# Non integer elements are filtered out:

collect(AddOneIfInt(), Any[3, nothing, 2.0, missing, 5])

# ## Stateful transducer
#
# `AddOneIfInt` is a stateless transducer which is very easy to
# implement.  A stateful transducer needs a bit more code.

using Transducers: start, complete, wrap, unwrap, wrapping
using Random

# Let's define a transducer that spits out a random past element from
# the buffer:

struct RandomRecall <: Transducer
    history::Int
    seed::Int
end
RandomRecall() = RandomRecall(3, 0)
nothing  # hide

# A stateful transducer needs to implement [`Transducers.start`](@ref)
# to "allocate" its private state.  Here, the private state is a
# `buffer` and a random number generator state `rng`:

function Transducers.start(rf::R_{RandomRecall}, result)
    buffer = []
    rng = MersenneTwister(xform(rf).seed)
    private_state = (buffer, rng)
    return wrap(rf, private_state, start(inner(rf), result))
end

# Stateful transducer needs to unwrap its private state inside
# [`Transducers.next`](@ref) and then re-wrap it.  There is a helper
# function [`Transducers.wrapping`](@ref) does that with the `do`
# block:

function Transducers.next(rf::R_{RandomRecall}, result, input)
    wrapping(rf, result) do (buffer, rng), iresult
#+
# Pickup a random element to be passed to the inner reducing function.
# Replace it with the new incoming one in the buffer:
        if length(buffer) < xform(rf).history
            push!(buffer, input)
            iinput = rand(rng, buffer)
        else
            i = rand(rng, 1:length(buffer))
            iinput = buffer[i]
            buffer[i] = input
        end
#+
# Call the inner reducing function.  Note that `iresult` unwrapped by
# [`Transducers.wrapping`](@ref) must be passed to `next`:
        iresult = next(inner(rf), iresult, iinput)
        return (buffer, rng), iresult
    end
end

# Any transducer with custom [`Transducers.start`](@ref) must have a
# corresponding [`Transducers.complete`](@ref).  It is responsible for
# unwrapping the `result` and call the `complete` for the inner
# reducing function.

function Transducers.complete(rf::R_{RandomRecall}, result)
    _private_state, inner_result = unwrap(rf, result)
    return complete(inner(rf), inner_result)
end

# Here is how it works:

recall_out1 = begin  # hide
collect(RandomRecall(), 1:5)
end  # hide

# Indeed, it picks up some random elements from the past elements.
# With slightly more complex transducer:

collect(Filter(isodd) |> RandomRecall() |> Filter(x -> x > 10) |> Take(5), 1:100)

# Note that [`Transducers.complete`](@ref) can do more than `unwrap`
# and `complete`.  It is useful for, e.g., flushing the buffer.

function Transducers.complete(rf::R_{RandomRecall}, result)
    (buffer, _), iresult = unwrap(rf, result)
    for x in buffer
#+
# Note that inner `next` can be called more than one time inside
# `next` and `complete`:
        iresult = next(inner(rf), iresult, x)
    end
#+
# `complete` for inner reducing function must be called exactly once:
    return complete(inner(rf), iresult)
end

# This then adds 3 (`= RandomRecall().history`) more elements to the
# output:

recall_out2 = begin  # hide
collect(RandomRecall(), 1:5)
end  # hide


text:
mutable struct Node
    value::Int64
    parent::Union{Node, Nothing}
    left::Node
    right::Node

    Node(value::Int64) = new(value, nothing)
    Node(value::Int64, parent::Union{Node, Nothing}) = new(value, parent)
end

function insert(value::Int64, root::Node)::Node
    if value ≤ root.value
        if !isdefined(root, :left)
            root.left = Node(value, root)
        else
            root.left = insert(value, root.left)
        end
    else
        if !isdefined(root, :right)
            root.right = Node(value, root)
        else
            root.right = insert(value, root.right)
        end
    end
    return root
end

function find(value::Int64, root::Node)::Union{Node, Nothing}
    if value == root.value
        return root
    elseif value < root.value && isdefined(root, :left)
        return find(value, root.left)
    elseif value > root.value && isdefined(root, :right)
        return find(value, root.right)
    else
        return nothing
    end
end

function delete(value::Int64, root::Node)::Union{Node, Nothing}

end

root = Node(9)
root = insert(10, root)
show(find(10, root))
println("")
show(find(11, root))

max_stars_repo_name: beer-and-sanma/Karli

max_issues_repo_name: beer-and-sanma/Karli

max_forks_repo_name: beer-and-sanma/Karli


text:
#*****************************************************************************#
# Initialization functions for stochastic integrators                         #
#*****************************************************************************#

"Create integrator for stochastic explicit Runge-Kutta tableau."
function Integrators.Integrator(equation::SDE, tableau::TableauSERK, Δt; kwargs...)
    IntegratorSERK(equation, tableau, Δt; kwargs...)
end

"Create integrator for weak explicit Runge-Kutta tableau."
function Integrators.Integrator(equation::SDE, tableau::TableauWERK, Δt; kwargs...)
    IntegratorWERK(equation, tableau, Δt; kwargs...)
end

"Create integrator for stochastic fully implicit Runge-Kutta tableau."
function Integrators.Integrator(equation::SDE, tableau::TableauSIRK, Δt; kwargs...)
    IntegratorSIRK(equation, tableau, Δt; kwargs...)
end

"Create integrator for stochastic fully implicit partitioned Runge-Kutta tableau."
function Integrators.Integrator(equation::PSDE, tableau::TableauSIPRK, Δt; kwargs...)
    IntegratorSIPRK(equation, tableau, Δt; kwargs...)
end

"Create integrator for stochastic fully implicit split partitioned Runge-Kutta tableau."
function Integrators.Integrator(equation::SPSDE, tableau::TableauSISPRK, Δt; kwargs...)
    IntegratorSISPRK(equation, tableau, Δt; kwargs...)
end

"Create integrator for weak fully implicit Runge-Kutta tableau."
function Integrators.Integrator(equation::SDE, tableau::TableauWIRK, Δt; kwargs...)
    IntegratorWIRK(equation, tableau, Δt; kwargs...)
end

#*****************************************************************************#
# Integration functions for stochastic integrators                         #
#*****************************************************************************#

function Integrators.integrate!(int::StochasticIntegrator{DT,TT}, sol::Solution{AT,TT}, asol::AtomicSolution{DT,TT}, m::Int, n::Int) where {DT, TT, AT <: AbstractArray{DT}}
    # copy the increments of the Brownian Process
    get_increments!(sol, asol, n, m)

    integrate_common!(int, sol, asol, m, n)
end


text:
# Simulate suprcontinuum generation for parameters similar to Fig. 3 of
# Dudley et al, RMP 78 1135 (2006)
# Parameters from www.scgbook.info
T_T0 = 0.0284ps
T_FWHM = T_T0 * 2*log(1+sqrt(2)) # pulse time paramater passed as FWHM
betacoeff = [0.0/m, 0.0ps/m, -11.83e-3ps^2/m, 8.1038e-5ps^3/m, -9.5205e-8ps^4/m,
             2.0737e-10ps^5/m, -5.3943e-13ps^6/m,
             1.13486e-15ps^7/m, -2.5495e-18ps^8/m, 3.0524e-21ps^9/m, -1.7140e-24ps^10/m]
γnl = 0.11W/m
ω0 = 2π*c0/(835nm) |> THz# c is m/ps, exported from NonlinearOptics
P_pulse = 100.0W
fiber_length = 0.15m
dudley = ToyModel(;ω0=ω0, betacoeff=betacoeff,
                   nonlinearcoeff=γnl,
                   linearloss=0.0/m,
                   power_in=P_pulse,
                   pulsetime=T_FWHM,
                   length=fiber_length,
                   has_shock=true,
                   has_raman=true)

sol_dudley = solve(dudley, DynamicNLSE(); time_pts=2^13, time_window=12.5ps)
# solves in 1.3 seconds, compared to 27 seconds for scgbook code!


text:
using Distributions, StatsBase, ProgressMeter

function double_iex_multinomial_edit_accept_reject!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    posterior::SimPosterior,
    γ_curr::Float64,
    mcmc::SimIexInsertDeleteWithKick,
    P::CumCondProbMatrix,
    aux_data::InteractionSequenceSample{Int},
    suff_stat_curr::Float64
    )

    N = length(S_curr)
    dist = posterior.dist
    V = posterior.V
    K_inner, K_outer = (posterior.K_inner, posterior.K_outer)
    K_in_lb, K_in_ub = (K_inner.l, K_inner.u)
    data = posterior.data
    γ_prior = posterior.S_prior.γ
    mode_prior = posterior.S_prior.mode

    aux_mcmc = mcmc.aux_mcmc

    δ = rand(1:mcmc.ν_ed)  # Number of edits to enact
    rem_edits = δ # Remaining edits to allocate
    j = 0 # Keeps track how many interaction have been edited
    log_ratio = 0.0

    for i = 1:N

        # If at end we just assign all remaining edits to final interaction
        if i == N
            δ_tmp = rem_edits
        # Otherwise we sample the number of edits via rescaled Binomial
        else
            p = 1/(N-i+1)
            δ_tmp = rand(Binomial(rem_edits, p)) # Number of edits to perform on ith interaction
        end

        # println("   Index $i getting $δ_tmp edits...")
        # If we sampled zero edits we skip to next iteration
        if δ_tmp == 0
            continue
        else
            j += 1 # Increment j
            # Make edits ....
            @inbounds n = length(S_curr[i])
            d = rand(0:min(n-K_in_lb, δ_tmp))
            m = n + δ_tmp - 2*d

            # Catch invalid proposals
            if (m > K_in_ub)
                # Here we just reject the proposal
                for i in 1:N
                    copy!(S_prop[i], S_curr[i])
                end
                return 0, suff_stat_curr
            end

            ind_del = view(mcmc.ind_del, 1:d)
            ind_add = view(mcmc.ind_add, 1:(δ_tmp-d))
            vals_del = view(mcmc.vals, 1:d)

            # Sample indexing info and new entries (all in-place)
            StatsBase.seqsample_a!(1:n, ind_del)
            StatsBase.seqsample_a!(1:m, ind_add)

            # *** HERE IS DIFFERENT FROM MODEL SAMPLER ***
            # The delete_insert_informed() function does the sampling + editing
            log_ratio += delete_insert_informed!(
                S_prop[i],
                ind_del, ind_add, vals_del,
                P)

            mcmc.ind_update[j] = i # Store which interaction was updated

            # Add to log_ratio
            # log_prod_term += log(b - a + 1) - log(ub(m, δ_tmp) - lb(m, δ_tmp, model) +1)
            log_ratio += log(min(n-K_in_lb, δ_tmp)+1) - log(min(m-K_in_lb, δ_tmp)+1)

        end

        # Update rem_edits
        rem_edits -= δ_tmp

        # If no more left terminate
        if rem_edits == 0
            break
        end

    end

    aux_model = SIM(
        S_prop, γ_curr,
        dist,
        V,
        K_inner,
        K_outer
        )

    draw_sample!(aux_data, aux_mcmc, aux_model)

    aux_log_lik_ratio = -γ_curr * (
        mapreduce(x -> dist(x, S_curr), + , aux_data)
        - mapreduce(x -> dist(x, S_prop), +, aux_data)
    )

    suff_stat_prop = mapreduce(x -> dist(x, S_prop), + , data)
    log_lik_ratio = -γ_curr * (
        suff_stat_prop - suff_stat_curr
    )

    log_prior_ratio = -γ_prior * (
        dist(S_prop, mode_prior) - dist(S_curr, mode_prior)
    )

    log_multinom_term = log_multinomial_ratio(S_curr, S_prop)

    # Log acceptance probability
    log_α = log_lik_ratio + log_prior_ratio + aux_log_lik_ratio + log_ratio + log_multinom_term

    # Accept-reject step. Use info in mcmc.ind_update to know which interaction are to be copied over
    if log(rand()) < log_α
        for i in view(mcmc.ind_update, 1:j)
            copy!(S_curr[i], S_prop[i])
        end
        return 1, suff_stat_prop
    else
        for i in view(mcmc.ind_update, 1:j)
            copy!(S_prop[i], S_curr[i])
        end
        return 0, suff_stat_curr
    end
end

function double_iex_flip_accept_reject!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    posterior::SimPosterior,
    γ_curr::Float64,
    mcmc::SimIexInsertDeleteWithKick,
    P::CumCondProbMatrix,
    aux_data::InteractionSequenceSample{Int},
    suff_stat_curr::Float64
    )

    dist = posterior.dist
    V = posterior.V
    data = posterior.data
    γ_prior = posterior.S_prior.γ
    mode_prior = posterior.S_prior.mode
    K_inner, K_outer = (posterior.K_inner, posterior.K_outer)

    aux_mcmc = mcmc.aux_mcmc
    lengths = length.(S_curr)

    δ = rand(1:min(mcmc.ν_ed,sum(lengths)))
    # @show δ, lengths
    alloc = rand_restr_bins(lengths, δ)
    ind = mcmc.ind_add
    log_ratio = 0.0

    for (key,val) in pairs(alloc)
        ind_flip = view(ind, 1:val)
        StatsBase.seqsample_a!(1:lengths[key], ind_flip)
        log_ratio += flip_informed_excl!(
            S_prop[key],
            ind_flip,
            P
        )
    end

    aux_model = SIM(
        S_prop, γ_curr,
        dist,
        V,
        K_inner,
        K_outer
        )

    draw_sample!(aux_data, aux_mcmc, aux_model)

    aux_log_lik_ratio = -γ_curr * (
        mapreduce(x -> dist(x, S_curr), + , aux_data)
        - mapreduce(x -> dist(x, S_prop), +, aux_data)
    )

    suff_stat_prop = mapreduce(x -> dist(x, S_prop), + , data)
    log_lik_ratio = -γ_curr * (
        suff_stat_prop - suff_stat_curr
    )

    log_prior_ratio = -γ_prior * (
        dist(S_prop, mode_prior) - dist(S_curr, mode_prior)
    )

    log_multinom_term = log_multinomial_ratio(S_curr, S_prop)
    # Log acceptance probability
    log_α = log_lik_ratio + log_prior_ratio + aux_log_lik_ratio + log_ratio + log_multinom_term

    # Accept-reject step. Use info in mcmc.ind_update to know which interaction are to be copied over
    if log(rand()) < log_α
        for i in keys(alloc)
            copy!(S_curr[i], S_prop[i])
        end
        return 1, suff_stat_prop
    else
        for i in keys(alloc)
            copy!(S_prop[i], S_curr[i])
        end
        return 0, suff_stat_curr
    end

end

function double_iex_trans_dim_block_informed_accept_reject!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    γ_curr::Float64,
    posterior::SimPosterior,
    mcmc::SimIexInsertDeleteWithKick,
    p_ins::Categorical,
    aux_data::InteractionSequenceSample{Int},
    suff_stat_curr::Float64
    )

    K_inner, K_outer = (posterior.K_inner, posterior.K_outer)
    K_out_lb, K_out_ub = (K_outer.l, K_outer.u)
    K_in_ub = K_inner.u
    data = posterior.data
    dist = posterior.dist
    V = posterior.V
    γ_prior = posterior.γ_prior
    S_prior = posterior.S_prior

    ν_td = mcmc.ν_td
    ε_td = mcmc.ε_td
    ε_kick = mcmc.ε_kick
    curr_pointers = mcmc.curr_pointers
    prop_pointers = mcmc.prop_pointers
    aux_mcmc = mcmc.aux_mcmc

    log_ratio = 0.0

    N = length(S_curr)

    # Generate gamma proposal
    γ_tmp = rand_reflect(γ_curr, ε_td, 0.0, Inf)
    is_insert = rand(Bernoulli(0.5)) # Its relative value to current determines whether we insert or delete

    if is_insert
        γ_prop = γ_tmp + ε_kick
        ε = rand(1:ν_td) # How many to insert
        # Catch invalid proposal (ones which have zero probability)
        if (N + ε) > K_out_ub
            # Make no changes and imediately reject
            return γ_curr, 0, suff_stat_curr
        end
        ind_tr_dim = view(mcmc.ind_td, 1:ε) # Storage for where to insert
        StatsBase.seqsample_a!(1:(N+ε), ind_tr_dim) # Sample where to insert
        log_ratio += imcmc_multi_insert_prop_sample!(
            S_curr, S_prop,
            mcmc,
            ind_tr_dim,
            V,
            K_in_ub
            ) # Enact move and catch log ratio term
    else
        γ_prop = γ_tmp - ε_kick
        ε = rand(1:min(ν_td, N)) # How many to delete
        # Catch invalid proposal (would go to empty inter seq)
        if (N - ε) < K_out_lb
            return γ_curr, 0, suff_stat_curr
        end
        ind_tr_dim = view(mcmc.ind_td, 1:ε) # Storage
        StatsBase.seqsample_a!(1:N, ind_tr_dim) # Sample which to delete
        log_ratio += imcmc_multi_delete_prop_sample!(
            S_curr, S_prop,
            mcmc,
            ind_tr_dim,
            V
            ) # Enact move and catch log ratio
    end

    # Now do accept-reject step (**THIS IS WHERE WE DIFFER FROM MODEL SAMPLER***)
    aux_model = SIM(
        S_prop, γ_prop,
        dist,
        V,
        K_inner,
        K_outer
        )

    aux_data = draw_sample(aux_mcmc, aux_model)

    aux_suff_stat_curr = mapreduce(x -> dist(x, S_curr), + , aux_data)
    aux_suff_stat_prop = mapreduce(x -> dist(x, S_prop), +, aux_data)
    aux_log_lik_ratio = (-γ_curr * aux_suff_stat_curr) + (γ_prop * aux_suff_stat_prop)

    suff_stat_prop = mapreduce(x -> dist(x, S_prop), + , data)
    log_lik_ratio = (-γ_prop * suff_stat_prop) + (γ_curr * suff_stat_curr)

    log_prior_ratio = -S_prior.γ *(dist(S_prop, S_prior.mode) - dist(S_curr, S_prior.mode))
    log_prior_ratio += logpdf(γ_prior, γ_prop) - logpdf(γ_prior, γ_curr)

    log_multinom_term = log_multinomial_ratio(S_curr, S_prop)

    # Log acceptance probability
    log_α = log_lik_ratio + log_prior_ratio + aux_log_lik_ratio + log_ratio + log_multinom_term

    # println("Acceptance prob: $(exp(log_α))")
    # @show is_insert, suff_stat_curr, suff_stat_prop, log_lik_ratio, aux_suff_stat_curr, aux_suff_stat_prop, aux_log_lik_ratio, log_ratio, log_multinom_term, log_prior_ratio, log_α, exp(log_α)
    # Note that we copy interactions between S_prop (resp. S_curr) and prop_pointers (resp .curr_pointers) by hand.
    println("Currently at $(S_curr) and $γ_curr")
    println("Proposing $(S_prop) and $γ_prop  ")
    @show log_ratio
    println("Log acc prob: $log_α \n")
    if log(rand()) < log_α
        if is_insert
            for i in ind_tr_dim
                migrate!(S_curr, curr_pointers, i, 1)
                copy!(S_curr[i], S_prop[i])
            end
        else
            for i in Iterators.reverse(ind_tr_dim)
            migrate!(curr_pointers , S_curr, 1, i)
            end
        end
        return γ_prop, 1, suff_stat_prop
    else
        # Here we must delete the interactions which were added to S_prop

        if is_insert
            for i in Iterators.reverse(ind_tr_dim)
                migrate!(prop_pointers, S_prop, 1, i)
            end
        # Or reinsert the interactions which were deleted
        else
            for i in ind_tr_dim
                migrate!(S_prop, prop_pointers, i, 1)
                copy!(S_prop[i], S_curr[i])
            end
        end
        return γ_curr, 0, suff_stat_curr
    end

end

# Joint Distribution
# ------------------

# Gamma accept-reject

function accept_reject_gamma!(
    γ_curr::Float64,
    S_curr::InteractionSequence{Int},
    posterior::SimPosterior,
    mcmc::SimIexInsertDeleteWithKick,
    aux_data::InteractionSequenceSample{Int},
    suff_stat_curr::Float64
    )

    ε = mcmc.ε_ed
    aux_mcmc = mcmc.aux_mcmc

    γ_prop = rand_reflect(γ_curr, ε, 0.0, Inf)

    aux_model = SIM(
        S_curr, γ_prop,
        posterior.dist,
        posterior.V,
        posterior.K_inner, posterior.K_outer
        )
    draw_sample!(aux_data, aux_mcmc, aux_model)

    # Accept reject

    log_lik_ratio = (γ_curr - γ_prop) * suff_stat_curr
    aux_log_lik_ratio = (γ_prop - γ_curr) * sum_of_dists(aux_data, S_curr, posterior.dist)

    log_α = (
        logpdf(posterior.γ_prior, γ_prop)
        - logpdf(posterior.γ_prior, γ_curr)
        + log_lik_ratio + aux_log_lik_ratio
    )
    if log(rand()) < log_α
        return γ_prop, 1
    else
        return γ_curr, 0
    end
end

function accept_reject!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    γ_curr::Float64,
    posterior::SimPosterior,
    mcmc::SimIexInsertDeleteWithKick,
    P::CumCondProbMatrix,
    p_ins::Categorical,
    aux_data::InteractionSequenceSample{Int},
    acc_count::Vector{Int},
    count::Vector{Int},
    suff_stat_curr::Float64
    )

    β = mcmc.β
    if rand() < β
        # Here we (i) do update move on mode conditonal, and (ii) do update in gamma conditional
        if rand() < 0.5
            was_accepted, suff_stat_curr = double_iex_multinomial_edit_accept_reject!(
                S_curr, S_prop,
                posterior, γ_curr,
                mcmc, P,
                aux_data,
                suff_stat_curr
            )
            acc_count[1] += was_accepted
            count[1] += 1
        else
            was_accepted, suff_stat_curr = double_iex_flip_accept_reject!(
                S_curr, S_prop,
                posterior, γ_curr,
                mcmc, P,
                aux_data,
                suff_stat_curr
            )
            acc_count[2] += was_accepted
            count[2] += 1
        end

        # And do accept-reject on dispersion
        γ_curr, was_accepted =  accept_reject_gamma!(
            γ_curr,
            S_curr,
            posterior,
            mcmc,
            aux_data,
            suff_stat_curr,
        )
        acc_count[3] += was_accepted
        count[3] += 1

    else
        γ_curr, was_accepted, suff_stat_curr = double_iex_trans_dim_block_informed_accept_reject!(
                S_curr, S_prop,
                γ_curr, posterior,
                mcmc, p_ins,
                aux_data,
                suff_stat_curr
        )

        acc_count[4] += was_accepted
        count[4] += 1
    end
    return γ_curr, suff_stat_curr
end

function draw_sample!(
    sample_out_S::Union{InteractionSequenceSample{Int},SubArray},
    sample_out_gamma::Union{Vector{Float64},SubArray},
    mcmc::SimIexInsertDeleteWithKick,
    posterior::SimPosterior;
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    S_init::InteractionSequence{Int}=sample_frechet_mean(posterior.data, posterior.dist),
    γ_init::Float64=5.0,
    loading_bar::Bool=true
    )

    if loading_bar
        iter = Progress(
            length(sample_out_S) * lag + burn_in, # How many iters
            1,  # At which granularity to update loading bar
            "Chain for n = $(posterior.sample_size) (joint)....")  # Loading bar. Minimum update interval: 1 second
    end

    # Define aliases for pointers to the storage of current vals and proposals
    curr_pointers = mcmc.curr_pointers
    prop_pointers = mcmc.prop_pointers
    aux_mcmc = mcmc.aux_mcmc

    S_curr = InteractionSequence{Int}()
    S_prop = InteractionSequence{Int}()
    for i in 1:length(S_init)
        migrate!(S_curr, curr_pointers, i, 1)
        migrate!(S_prop, prop_pointers, i, 1)
        copy!(S_curr[i], S_init[i])
        copy!(S_prop[i], S_init[i])
    end
    γ_curr = copy(γ_init)

    sample_count = 1 # Keeps which sample to be stored we are working to get
    i = 1 # Keeps track all samples (included lags and burn_ins)

    acc_count = [0,0,0,0]  #[edit,flip,gamma,block_td]
    count = [0,0,0,0]

    aux_data = [[Int[]] for i in 1:posterior.sample_size]
    # Initialise the aux_data
    aux_model = SIM(
        S_curr, γ_curr,
        posterior.dist,
        posterior.V,
        posterior.K_inner,
        posterior.K_outer)
    draw_sample!(aux_data, aux_mcmc, aux_model)
    # Initialise sufficient statistic
    suff_stat_curr = mapreduce(
        x -> posterior.dist(S_curr, x),
        +,
        posterior.data
    )
    suff_stats = Float64[suff_stat_curr]
    # Get informed proposal matrix
    P, vmap, vmap_inv = get_informed_proposal_matrix(posterior, mcmc.α)
    p_ins = get_informed_insertion_dist(posterior, mcmc.α)
    while sample_count ≤ length(sample_out_S)
        # Store values
        if (i > burn_in) & (((i-1) % lag)==0)
            sample_out_S[sample_count] = deepcopy(S_curr)
            sample_out_gamma[sample_count] = copy(γ_curr)
            sample_count += 1
        end

        γ_curr, suff_stat_curr = accept_reject!(
            S_curr, S_prop, γ_curr,
            posterior,
            mcmc, P, p_ins,
            aux_data,
            acc_count, count,
            suff_stat_curr
        )

        if loading_bar
            next!(iter)
        end
        i += 1
    end

    for i in 1:length(S_curr)
        migrate!(curr_pointers, S_curr, 1, 1)
        migrate!(prop_pointers, S_prop, 1, 1)
    end

    ed_acc_prob = acc_count[1] / count[1]
    flip_acc_prob = acc_count[2] / count[2]
    γ_acc_prob = acc_count[3] / count[3]
    block_td_acc_prob = acc_count[4] / count[4]
    return ed_acc_prob, flip_acc_prob, γ_acc_prob, block_td_acc_prob, suff_stats
end

function draw_sample(
    mcmc::SimIexInsertDeleteWithKick,
    posterior::SimPosterior;
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    S_init::InteractionSequence{Int}=sample_frechet_mean(posterior.data, posterior.dist),
    γ_init::Float64=5.0,
    loading_bar::Bool=true
    )

    sample_out_S = Vector{InteractionSequence{Int}}(undef, desired_samples)
    sample_out_gamma = Vector{Float64}(undef, desired_samples)

    draw_sample!(
        sample_out_S,
        sample_out_gamma,
        mcmc,
        posterior,
        burn_in=burn_in, lag=lag,
        S_init=S_init, γ_init=γ_init,
        loading_bar=loading_bar
    )

    return (S=sample_out_S, gamma=sample_out_gamma)
end

function (mcmc::SimIexInsertDeleteWithKick)(
    posterior::SimPosterior;
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    S_init::InteractionSequence{Int}=sample_frechet_mean(posterior.data, posterior.dist),
    γ_init::Float64=5.0,
    loading_bar::Bool=true
    )

    sample_out_S = Vector{InteractionSequence{Int}}(undef, desired_samples)
    sample_out_gamma = Vector{Float64}(undef, desired_samples)

    ed_acc_prob, flip_acc_prob, γ_acc_prob, block_td_acc_prob, suff_stats = draw_sample!(
        sample_out_S,
        sample_out_gamma,
        mcmc,
        posterior,
        burn_in=burn_in, lag=lag,
        S_init=S_init, γ_init=γ_init,
        loading_bar=loading_bar
    )

    p_measures = Dict(
        "Dipsersion acceptance probability" => γ_acc_prob,
        "Edit allocation acceptance probability" => ed_acc_prob,
        "Flip allocation acceptance probability" => flip_acc_prob,
        "Trans-dimensional (block) acceptance probability" => block_td_acc_prob
        )

    return SimPosteriorMcmcOutput(
        sample_out_S,
        sample_out_gamma,
        posterior,
        suff_stats,
        p_measures
    )

end


text:
"Return the Fourier coefficients of the bimodal wrapped Gaussian."
function GaussCoef(θ)
    φ₁ = θ[1]
    φ₂ = θ[2]
    σ = θ[3]
    a1 = exp(-(σ^2) / 2) * (cos(φ₁) + cos(φ₂))
    b1 = exp(-(σ^2) / 2) * (sin(φ₁) + sin(φ₂))
    a2 = exp(-2(σ^2)) * (cos(2φ₁) + cos(2φ₂))
    b2 = exp(-2(σ^2)) * (sin(2φ₁) + sin(2φ₂))
    return a1, b1, a2, b2
end
"Return the Fourier coefficients from the matrix of cross spectra."
function freqdir_fouriercoef2(S)
    denom = sqrt(abs((S[2, 2] + S[3, 3]) * S[1, 1]))
    a1 = imag(S[2, 1]) / denom
    b1 = imag(S[3, 1]) / denom
    a2 = real((S[2, 2] - S[3, 3]) / (S[2, 2] + S[3, 3]))
    b2 = 2 * real((S[2, 3]) / (S[2, 2] + S[3, 3]))
    return a1, b1, a2, b2
end
"Function to compute the least squares objective for method of moments from Ewans 1998."
function momentLS(θ, S)
    a1h, b1h, a2h, b2h = freqdir_fouriercoef2(S)
    a1, b1, a2, b2 = GaussCoef(θ)
    return (a1 - a1h)^2 + (b1 - b1h)^2 + (a2 - a2h)^2 + (b2 - b2h)^2
end
"Function to optimise moments objective at a given frequency."
function getPars(θ₀, S)
    res = optimize(x -> momentLS(x, S), θ₀)
    if Optim.converged(res)
        return [
            min(res.minimizer[1], res.minimizer[2]),
            max(res.minimizer[1], res.minimizer[2]),
            res.minimizer[3],
        ]
    else
        return fill(NaN, size(θ₀))
    end
end
"Function to compute the peak seperation under a given model at a frequency."
function truepeakseperation(θ, ω, ωₚ)
    β = θ[2]
    ν = θ[3]
    return β / exp(ν * ((ω > ωₚ) * ((ωₚ / ω) - 1.0) + 1.0))
end

"Function to compute the means and angular width of the bimodal spreading at a frequency."
function truemeansigma(θ, ω, ωₚ)
    σₗ = θ[4]
    σᵣ = θ[5]
    φₘ = θ[1]
    φₛ = truepeakseperation(θ, ω, ωₚ) / 2
    return φₘ - φₛ, φₘ + φₛ, (σₗ - (σᵣ / 3) * (4 * (ωₚ / ω)^2 - (ωₚ / ω)^8))
end

"Function to compute the means and angular width of the bimodal spreading at a frequency returning a vector."
function truemeansigmavec(θ, ω, ωₚ)
    φₘ₁, φₘ₂, σ = truemeansigma(θ, ω, ωₚ)
    return [φₘ₁, φₘ₂, σ]
end

##
"Objective function struct for moments. Requires known peak frequency from some other fit."
struct MomentObjective{T}
    moments::Vector{Vector{Float64}}
    Ω::T
    ωₚ::Float64
    function MomentObjective(ts, Δ, θ₀, ωₚ; ωₗ = 0, ωᵤ = Inf)
        B = BartlettPeriodogram(ts, Δ)[end÷2+1:end]
        S = B[ωₗ.<=B.Ω.<=ωᵤ]
        momentest = [getPars(truemeansigmavec(θ₀, ω, ωₚ), s) for (ω, s) ∈ zip(S.Ω, S.ordinates)]
        goodestind = [!isnan(m[1]) for m ∈ momentest]
        new{typeof(S.Ω)}(momentest[goodestind], S.Ω[goodestind], ωₚ)
    end
end
function momentobjective_summand(m, ω, θ, ωₚ)
    φₘ₁, φₘ₂, σ = truemeansigma(θ, ω, ωₚ)
    return (min(φₘ₁, φₘ₂) - m[1])^2 + (max(φₘ₁, φₘ₂) - m[2])^2 + (σ - m[3])^2
end
function (f::MomentObjective)(θ)
    return sum(momentobjective_summand(m, ω, θ, f.ωₚ) for (m, ω) ∈ zip(f.moments, f.Ω))
end

function momentsfit(model, ts, Δ, θ₀; ωₗ = 0, ωᵤ = Inf)
    ## fit LS
    B = BartlettPeriodogram(ts[:, 1], Δ, 128)[end÷2+1:end]
    Bz = B[ωₗ.<=B.Ω.<=ωᵤ]
    res1 = optimize(θ -> marginalbounds(θ,ωₗ,ωᵤ) ? squarediff(model(θ), Bz) : Inf, θ₀[1:4])
    ## fit moments
    if Optim.converged(res1)
        momobj = MomentObjective(ts, Δ, θ₀[5:end], res1.minimizer[2]; ωₗ = ωₗ, ωᵤ = ωᵤ)
        res2 = optimize(θ -> spreadingbounds(θ) ? momobj(θ) : Inf, θ₀[5:end])
        if Optim.converged(res2)
            return (true, [res1.minimizer; res2.minimizer])
        end
    end
    return (false, θ₀)
end


text:
using Agents, Random

mutable struct Tree <: AbstractAgent
    id::Int
    pos::Tuple{Int, Int}
    status::Bool  # true is green and false is burning
end

function model_initiation(; f, p, griddims)
    space = GridSpace(griddims, moore = true)
    properties = Dict(:f => f, :p => p)
    forest = AgentBasedModel(Tree, space; properties=properties)
    for node in nodes(forest)
        add_agent!(node, forest, true)
    end
    return forest
end

# f is combustion prob., p is regrowth prob.
println("time to initialize model")
@time forest = model_initiation(f=0.0001, p=0.01, griddims=(400, 400))

# ## Defining the step!
# Because of the way the forest fire model is defined, we only need a
# stepping function for the model

function forest_step!(forest)
  for node in nodes(forest, by = :random)
    nc = get_node_contents(node, forest)
    ## the cell is empty, maybe a tree grows here
    if length(nc) == 0
        rand() ≤ forest.properties[:p] && add_agent!(node, forest, true)
    else
      tree = forest[nc[1]] # by definition only 1 agent per node
      if tree.status == false  # if it is has been burning, remove it.
        kill_agent!(tree, forest)
      else
        if rand() ≤ forest.properties[:f]  # the tree ignites spntaneously
          tree.status = false
        else  # if any neighbor is on fire, set this tree on fire too
          for cell in node_neighbors(node, forest)
            neighbors = get_node_contents(cell, forest)
            length(neighbors) == 0 && continue
            if any(n -> !forest.agents[n].status, neighbors)
              tree.status = false
              break
            end
          end
        end
      end
    end
  end
end

println("time to do a single step!")
step!(forest, dummystep, forest_step!) # compile
@time for i in 1:50
           step!(forest, dummystep, forest_step!)
       end


text:
@testset "datatip" begin
    @testset "local datatips" begin
        let str = """
            function localdatatip(word, column, row, startRow, context) # L0
              position = row - startRow                                 # L1
              ls = locals(context, position, column)                    # L2
              filter!(ls) do l                                          # L3
                l[:name] == word &&                                     # L4
                l[:line] < position                                     # L5
              end                                                       # L6
              # there should be zero or one element in `ls`             # L7
              map(l -> localdatatip(l, word, startRow), ls)             # L8
            end                                                         # L9
            """,
            localdatatip(word, line) = Atom.localdatatip(word, Inf, line + 1, 0, str)[1]

            @test localdatatip("row", 1) == 0 # line
            @test localdatatip("position", 2) == Dict(:type => :snippet, :value => "position = row - startRow") # binding string
            @test localdatatip("l", 4) == 3 # line
        end

        # ignore dot accessors
        let str = """
            function withdots(expr::CSTParser.EXPR)
                bind = CSTParser.bindingof(expr.args[1])
                val = bind.val
                return val
            end
            """,
            localdatatip(word, line) = Atom.localdatatip(word, Inf, line + 1, 0, str)[1]

            @test localdatatip("expr.args", 1)[:value] == "expr::CSTParser.EXPR"
            @test localdatatip("bind.val", 2)[:value] == "bind = CSTParser.bindingof(expr.args[1])"
        end

        # don't error on fallback case
        @test_nowarn @test Atom.localdatatip("word", 1, 1, 0, "") == []
    end

    @testset "code block search" begin
        using Base.Docs
        using Atom: searchcodeblocks

        @eval Main begin
            @doc """
                codeblocktest()

            this doc should n't be captured

            ```julia
            julia> codeblocktest()
            you should move away from atom-ide-ui asap
            ```
            """
            codeblocktest() = print("you should move away from atom-ide-ui asap")
        end

        let codes = searchcodeblocks(@doc codeblocktest)
            @test filter(codes) do code
                occursin(r"codeblocktest()", code)
            end |> !isempty
            @test filter(codes) do code
                occursin(r"this doc should n't be captured", code)
            end |> isempty
            @test filter(codes) do code
                occursin(r"julia> codeblocktest()", code)
            end |> !isempty
        end
    end

    @testset "global datatips" begin
        using Atom: globaldatatip

        ## method datatip
        @eval Main begin
            @doc """
                datatipmethodtest()

            this doc should be shown in datatip
            """
            datatipmethodtest() = nothing
        end

        let datatip = globaldatatip("Main", "datatipmethodtest")
            @test datatip isa Vector
            firsttip = datatip[1]
            secondtip = datatip[2]
            @test firsttip[:type] == :snippet
            @test occursin(r"datatipmethodtest()", firsttip[:value])
            @test secondtip[:type] == :markdown
            @test occursin(r"this doc should be shown in datatip", secondtip[:value])
        end

        ## variable datatip
        @eval Main datatipvariabletest = "this string should be shown in datatip"

        let datatip = globaldatatip("Main", "datatipvariabletest")
            @test datatip isa Vector
            firsttip = datatip[1]
            @test firsttip[:type] == :snippet
            @test occursin(r"this string should be shown in datatip", firsttip[:value])
        end
    end
end


text:
"""
	module ViewerGL

3D `interactive Viewer` for geometric and topological data.

Helper module for Julia's native ``OpenGL visualization``, forked from [Plasm.jl](https://github.com/plasm-language/pyplasm/tree/master/src/plasm.jl). To be used with geometric models and geometric expressions from [LinearAlgebraicRepresentation.jl](https://github.com/cvdlab/LinearAlgebraicRepresentation.jl), the simplest data structures for geometric and solid modeling :-)
"""
module ViewerGL

	using LinearAlgebra,StaticArrays
	using ModernGL
	using GLFW
	import Base:*

	include("Point.jl")
	include("Box.jl")
	include("Matrix.jl")
	include("Quaternion.jl")
	include("Frustum.jl")

	include("GLUtils.jl")
	include("GLVertexBuffer.jl")
	include("GLVertexArray.jl")
	include("GLMesh.jl")
	include("GLShader.jl")
	include("GLPhongShader.jl")

	include("Viewer.jl")
	include("Geometry.jl")
	include("GLText.jl")
	include("GLColorBuffer.jl")

	const GLFrame2 = GLAxis(Point3d(0,0,0),Point3d(1,1,0))
	const GLFrame = GLAxis(Point3d(0,0,0),Point3d(1,1,1))

end # module


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct OrderedGroup <: Edit
        on::Array{Symbol,1}
        var::Symbol
        val::Array{Any,1}
    end

maybe, if on and var are the same, we can just fill in groups? i'm thinking SCTG group.

# Arguments
- `on::Array{Symbol,1}`: name of columns containing information specific to a particular level
- `var::Symbol`: name of column containing information of what we will unstack on
- `val::Array{Any,1}`: ordered list of values to unstack on. If empty, unstack in order of appearance.
"""
mutable struct OrderedGroup <: Edit
    "name of columns containing information specific to a particular level"
    on::Array{Symbol,1}
    "name of column containing information of what we will unstack on"
    var::Symbol
    "ordered list of values to unstack on. If empty, unstack in order of appearance."
    val::Array{Any,1}
end

function OrderedGroup(; on, var, val, )
    OrderedGroup(on, var, val, )
end

"""Get [`OrderedGroup`](@ref) `on`."""
get_on(value::OrderedGroup) = value.on
"""Get [`OrderedGroup`](@ref) `var`."""
get_var(value::OrderedGroup) = value.var
"""Get [`OrderedGroup`](@ref) `val`."""
get_val(value::OrderedGroup) = value.val

"""Set [`OrderedGroup`](@ref) `on`."""
set_on!(value::OrderedGroup, val) = value.on = val
"""Set [`OrderedGroup`](@ref) `var`."""
set_var!(value::OrderedGroup, val) = value.var = val
"""Set [`OrderedGroup`](@ref) `val`."""
set_val!(value::OrderedGroup, val) = value.val = val

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
"""
    expect(σ,ψ), expect(σ,ρ)

Compute the expectation value of an operator \$σ\$ given a state ket \$|ψ⟩\$ or a density matrix \$ρ\$. The expectation value is defined as
```math
\\begin{align*}
E(σ,|ψ⟩) &= ⟨ψ|σ|ψ⟩, \\\\
E(σ,ρ) &= \\textrm{tr}(σρ).
\\end{align*}
```

A specialized method exists for vector of `Ket` or `Operator` inputs.
"""
expect(σ::Operator,ψ::Ket) = (dimsmatch(σ,ψ); dot(data(ψ),data(σ)*data(ψ)))
expect(ψ::Ket,σ::Operator) = expect(σ,ψ)
expect(σ::Operator,ρ::Operator) = (dimsmatch(σ,ρ); trace(data(σ)*data(ρ)))

# Faster expectation value methods when passing in many inputs
function expect!(res,σ::Operator,states::Vector{Ket{T,D}}) where {T,D}
    dimsmatch(σ,states[1])
    tmp = Vector{ComplexF64}(undef,prod(dims(σ)))
    for (i,ψ) in enumerate(states)
        mul!(tmp,data(σ),data(ψ))
        res[i] = dot(data(ψ),tmp)
    end
    return res
end

function expect!(res,σ::Operator,states::Vector{Operator{T,D}}) where {T,D}
    dimsmatch(σ,states[1])
    superσ = super(data(σ))
    tmp = Vector{ComplexF64}(undef,prod(dims(σ))^2)
    N = length(tmp); sqrtNp1 = isqrt(N)+1
    fill!(res, zero(eltype(tmp)))
    for (i,ρ) in enumerate(states)
        mul!(tmp,superσ,vec(data(ρ)))
        for n = 1:sqrtNp1:N
            res[i] += tmp[n]
        end
    end
    return res
end

"""
    levelprobs(ψ), levelprobs(ψ,s)

Compute the level occupation probabilities. For a `Ket`, this simply corresponds to the absolute square of the amplitude of each level. For an `Operator`, the function returns the diagonal.

A system index, or vector of indices, can be passed as a second argument. In that case, the full system will first be partial traced to keep only the desired index. Level occupation probabilities are then calculated from the resulting reduced density matrix. If a vector of indices is passed, occupation probabilities are calculated for a fully reduced density matrix for each index.

Specialized methods exists for vectors of `Ket` or `Operator` inputs.
"""
levelprobs(ψ::Ket{T,N}) where {T,N} = abs2(ψ)
levelprobs(ψ::Ket{T,N},s::Integer) where {T,N} = real(diag(ptrace(ψ,ntuple_sans_m(s,Val(N)))))
levelprobs(ρ::Operator{T,N}) where {T,N} = real(diag(ρ))
levelprobs(ρ::Operator{T,N},s::Integer) where {T,N} = real(diag(ptrace(ρ,ntuple_sans_m(s,Val(N)))))
levelprobs(ψ::QuObject,S::AbstractVector) = map(s->levelprobs(ψ,s),S)

# Faster levelprobs methods when passing in many inputs
function levelprobs(states::Vector{<:QuObject},S::Union{Integer,AbstractVector})
    N = length(states)
    D = dims(states[1])
    probs = map(S) do s
        s > length(D) && throw(ArgumentError("System index $s is larger than the number of systems, $M."))
        P = Matrix{Float64}(undef,N,D[s])
        for n = 1:N
            P[n,:] = levelprobs(states[n],s)
        end
        return P
    end
    return probs
end

function levelprobs(states::Vector{<:QuObject})
    N = length(states)
    P = Matrix{Float64}(undef,N,prod(dims(states[1])))
    for n = 1:N
        P[n,:] = levelprobs(states[n])
    end
    return P
end

"""
    purity(ρ)

Compute the purity of a quantum state `ρ`, defined by
```math
\\begin{align*}
γ = \\textrm{tr}(ρ^†ρ).
\\end{align*}
```
This is the same as the square of the Frobenius norm of a matrix.

The purity of a quantum state is bounded below by `1/d`, for a maximally mixed state, and above by `1`, for a pure state. Note that the purity of a `Ket` is `1` by definition.

The purity is related to the linear entropy ``S_L`` by ``S_L = 1-γ``.
"""
function purity(ρ::Operator)
    checkdensityop(ρ)
    return norm(ρ)^2
end

purity(ψ::QuVector) = real(one(eltype(ψ)))

"""
    fidelity(ρ,σ), fidelity(ρ,ψ), fidelity(ψ,ϕ)

Compute the fidelity between density matrices \$ρ\$ and \$σ\$, a density matrix \$ρ\$ and a ket \$|ψ⟩\$, or two kets \$|ψ⟩\$ and \$|ϕ⟩\$. The fidelity in those three cases is defined as
```math
\\begin{align*}
F(ρ,σ) &= \\textrm{tr}\\sqrt{ρ^{1/2}σρ^{1/2}}, \\\\
F(ρ,|ψ⟩) &= \\sqrt{⟨ψ|ρ|ψ⟩}, \\\\
F(|ψ⟩,|ϕ⟩) &= \\left|⟨ψ|ϕ⟩\\right|.
\\end{align*}
```
See also [`fidelity2`](@ref), which is the square of the state fidelity.
"""
function fidelity(ρ::Operator,σ::Operator)
    # ref: https://quantumcomputing.stackexchange.com/q/9891
    dimsmatch(ρ,σ)
    checkdensityop(ρ)
    checkdensityop(σ)
    D = real.(eigvals(ρ*σ))
    return sum(d -> d>0 ? √(d) : zero(d), D)
end
fidelity(ρ::Operator,ψ::Ket) = sqrt(fidelity2(ρ,ψ))
fidelity(ψ::Ket,ρ::Operator) = fidelity(ρ,ψ)
fidelity(ψ::Ket,ϕ::Ket) = abs(dot(ψ,ϕ))

"""
    fidelity2(ρ,ψ)

Compute the Uhlmann state fidelity between density matrices \$ρ\$ and \$σ\$, a density matrix \$ρ\$ and a ket \$|ψ⟩\$, or two kets \$|ψ⟩\$ and \$|ϕ⟩\$. The Uhlmann state fidelity is simply defined as the square of the "regular" state [`fidelity`](@ref).
"""
fidelity2(ρ::Operator,σ::Operator) = abs2(fidelity(ρ,σ))
function fidelity2(ρ::Operator,ψ::Ket)
    checkdensityop(ρ)
    return abs(expect(ρ,ψ))
end
fidelity2(ψ::Ket,ρ::Operator) = fidelity2(ρ,ψ)
fidelity2(ψ::Ket,ϕ::Ket) = abs2(dot(ψ,ϕ))

function fidelity2(A::NTuple{N,T},ϕ::Ket) where {N,T<:Ket}
    # If sum(A) = ψ is a valid Ket, then this function will calculate the state fidelity between ψ and ϕ, ignoring all relative phases between the different parts of ψ (as well as the global phase, of course)
    # F² = ∑ᵢ(|⟨ψᵢ,ϕ⟩|² + ∑ⱼ₌₁ⁱ⁻¹2*|⟨ψᵢ,ϕ⟩|*|⟨ψⱼ,ϕ⟩|)
    f2 = fidelity2.(A,ϕ) # (|⟨ψ₁,ϕ⟩|², |⟨ψ₂,ϕ⟩|², ...)
    res = 0.0
    for i = 1:length(f2)
        res += f2[i]
        for j = 1:i-1
            res += 2*sqrt(f2[i])*sqrt(f2[j])
        end
    end
    return res
end

"""
    entanglement_fidelity(U,V)

Compute the entanglement fidelity of an imperfect quantum operation `V` with respect to an ideal operation `U`. The entanglement fidelity is given by the following formula:
```math
F_\\text{e}(\\mathcal{E}) = ⟨ϕ|(I⊗\\mathcal{E})(ϕ)|ϕ⟩,
```
where \$|ϕ⟩ = \\frac{1}{\\sqrt{d}}\\sum_{i=1}^d|i⟩|i⟩\$ is the maximally entangled state and \$\\mathcal{E}\$ is a trace preserving quantum operation.

When trying to find the fidelity of a quantum gate \$U\$, implemented imperfectly as \$V\$, then \$\\mathcal{E} = U^†∘V\$. In the case where \$U\$ and \$V\$ are a simple operators,
```math
\\begin{align*}
(I⊗\\mathcal{E})(ϕ) = (I⊗U^†V)|ϕ⟩⟨ϕ|(I⊗U^†V)^† \\quad\\text{and so}
\\end{align*}
```
```math
\\begin{align*}
F_\\text{e}(\\mathcal{E}) &= ⟨ϕ|(I⊗U^†V)|ϕ⟩⟨ϕ|(I⊗U^†V)^†|ϕ⟩ \\\\
  &= |⟨ϕ|(I⊗U^†V)|ϕ⟩|^2 \\\\
  &= |⟨U,V⟩|^2/d^2 \\\\
\\end{align*}
```

The entanglement fidelity is related to the [`average gate fidelity`](@ref gate_fidelity) by
```math
\\begin{align*}
\\overline{F}(\\mathcal{E}) =  \\frac{d F_\\text{e}(\\mathcal{E}) + 1}{d + 1}.
\\end{align*}
```
"""
function entanglement_fidelity(U::Operator,V::Operator,d::Integer=prod(dims(U)))
    dimsmatch(U,V)
    return abs2(inner(U,V))/d^2
end

"""
    gate_fidelity(U,V)

Compute the average gate fidelity of an imperfect quantum operation `V` with respect to an ideal operation `U`.

The average gate fidelity is related to the [`entanglement fidelity`](@ref entanglement_fidelity) by
```math
\\begin{align*}
\\overline{F}(\\mathcal{E}) =  \\frac{d F_\\text{e}(\\mathcal{E}) + 1}{d + 1}.
\\end{align*}
```
"""
function gate_fidelity(U::Operator,V::Operator,d::Integer=prod(dims(U)))
    dimsmatch(U,V)
    return (abs2(inner(U,V)) + d)/(d^2 + d)
end


text:
# # Crane

# ## Model

module crane    #jl

using ReachabilityAnalysis

A = [ 0.0           1.0               0.0   0.0   0.0           0.0;
     -0.417533     -3.1931759963     39.24  0.0 -14.825331     11.123344;
      0.0           0.0               0.0   1.0   0.0           0.0;
      0.0417533     0.31931759963    -4.905 0.0   1.4825331    -1.1123344;
      0.0638407957 -0.32473339016573  0.0   0.0  -3.7332068901 -0.7007592976;
      0.0853437452 -0.72366802635628  0.0   0.0  -5.9714023436 -2.2736115136]

function model(X0)
    S = @system(x' = Ax)
    return IVP(S, X0)
end

end #jl

# ## References

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
using Iterators

function sumSubset(set::Array{Int64,1}, subsetNumber::Int64)
    result = 0
    count = 0
    for i = 0:length(set)-1
        if ((2^i & subsetNumber) > 0)
            count += 1
            result += set[i+1]
        end
    end
    return (result, count)
end

function checkSubsets(set::Array{Int64,1})
    num_subsets = 2^length(set)-1
    for i = 1:num_subsets
        for j = 1:num_subsets
            if (i & j == 0) # treating i,j as bit flags, this is the condition for disjoint
                (sum_i, count_i) = sumSubset(set, i)
                (sum_j, count_j) = sumSubset(set, j)
                #println("i: $count_i, $sum_i; j: $count_j, $sum_j")
                if (sum_i == sum_j)
                    return false
                end
                if (count_i > count_j && sum_i <= sum_j)
                    return false
                end
                if (count_i < count_j && sum_i >= sum_j)
                    return false
                end
            end
        end
    end
    return true
end

# Main
best = nothing
best_sum = Inf

for i in subsets(13:48, 7) # Using rough heuristics -- this takes ~2min to run
    if (checkSubsets(i) && sum(i) < best_sum)
        println("Found one!")
        best_sum = sum(i)
        best = i
    end
end

@show join(best)

#=
Note this is relatively "brute-force". There's a more elegant way to do this which could "learn"
from past failures. e.g. once you know 10 & 15 add to 25, no set with these three numbers should be
considered.
=#

max_stars_repo_name: mebrunet/project-euler

max_issues_repo_name: mebrunet/project-euler

max_forks_repo_name: mebrunet/project-euler


text:
## Norms
function space_trace_norm(u1::Matrix{T}, xgrid::Vector, Nspacep::Int, ind::Int) where {T<:Real}
    out = zero(T)
    for i in 2:Nspacep
        out += (xgrid[i] - xgrid[i - 1]) * abs2(u1[i, ind])
    end
    sqrt(out)
end
function space_trace_norm(u1::Matrix{T}, xgrid::AbstractRange, Nspacep::Int, ind::Int) where {T<:Real}
    out = zero(T)
    for i in 2:Nspacep
        out += abs2(u1[i, ind])
    end
    sqrt(step(xgrid) * out)
end

function space_trace_dist(u1::Matrix{T}, v::Vector{T}, xgrid::Vector, Nspacep::Int, ind::Int) where {T<:Real}
    out = zero(T)
    for i in 2:Nspacep
        out += (xgrid[i] - xgrid[i - 1]) * abs2(u1[i, ind] - v[i])
    end
    sqrt(out)
end
function space_trace_dist(u1::Matrix{T}, v::Vector{T}, xgrid::AbstractRange, Nspacep::Int, ind::Int) where {T<:Real}
    out = zero(T)
    for i in 2:Nspacep
        out += abs2(u1[i, ind] - v[i])
    end
    sqrt(step(xgrid) * out)
end

function space_dist(u1::Matrix{T}, u2::Matrix{T}, xgrid::Vector, k::Real) where {T<:Real}
    NX, NT = size(u1)
    dout = zero(T)
    for i in 1:NT
        for j in 2:NX
            dout += (xgrid[j] - xgrid[j - 1]) * abs2(u1[j, i] - u2[j, i])
        end
    end
    sqrt(k * dout)
end
function space_dist(u1::Matrix{T}, u2::Matrix{T}, xgrid::AbstractRange, k::Real) where {T<:Real}
    NX, NT = size(u1)
    dout = zero(T)
    for i in 1:NT
        for j in 2:NX
            dout += abs2(u1[j, i] - u2[j, i])
        end
    end
    sqrt(k * step(xgrid) * dout)
end


text:
using Base.Test

# add at least 3 worker processes
if nworkers() < 3
    n = max(3, min(8, Sys.CPU_CORES))
    addprocs(n; exeflags=`--check-bounds=yes`)
end
@assert nprocs() > 3
@assert nworkers() >= 3

using DistributedArrays

@everywhere srand(1234 + myid())

include("darray.jl")


text:
using Gadfly, DataFrames, Compat, CSV

set_default_plot_size(6inch, 3inch)

plot(CSV.read(joinpath(dirname(@__DIR__), "data","issue120.csv")), x=:x1, y=:x2)


text:
### Error functions.

function error_function_selector(error::ASCIIString)
	if error == "absolute_error"
		return "absolute_error", absolute_error!, absolute_error_prime!
	elseif error == "cross_entropy"
		return "cross_entropy", cross_entropy!, cross_entropy_prime!
	else
		return "squared_error", squared_error!, squared_error_prime!
	end
end

function absolute_error!{T<:FloatingPoint}(YH::AbstractVector{T}, Y::AbstractVector{T}, E::AbstractVector{T}; eta::T=1e-10)
	@inbounds begin
		for i = 1:length(E)
			if YH[i] > Y[i]
				E[i] = YH[i] - Y[i]
			else
				E[i] = Y[i] - YH[i]
			end
		end
	end
end

function absolute_error_prime!{T<:FloatingPoint}(YH::AbstractVector{T}, Y::AbstractVector{T}, DE_DYH::AbstractVector{T}; eta::T=1e-10)
	@inbounds begin
		for i = 1:length(Y)
			if YH[i] > Y[i]
				DE_DYH[i] = 1.0
			else
				DE_DYH[i] = -1.0
			end
		end
	end
end

function cross_entropy!{T<:FloatingPoint}(YH::AbstractVector{T}, Y::AbstractVector{T}, E::AbstractVector{T}; eta::T=1e-10)
	@inbounds begin
		for i = 1:length(E)
			E[i] = -((Y[i] * log(YH[i] + eta)) + ((1.0 - Y[i]) * log(1.0 - YH[i] + eta)))
		end
	end
end

function cross_entropy_prime!{T<:FloatingPoint}(YH::AbstractVector{T}, Y::AbstractVector{T}, DE_DYH::AbstractVector{T}; eta::T=1e-10)
	@inbounds begin
		for i = 1:length(Y)
			DE_DYH[i] = ((1.0 - Y[i]) / (1.0 - YH[i] + eta)) - (Y[i] / (YH[i] + eta))
		end
	end
end

function squared_error!{T<:FloatingPoint}(YH::AbstractVector{T}, Y::AbstractVector{T}, E::AbstractVector{T})
	@inbounds begin
		for i = 1:length(E)
			E[i] = 0.5 * abs2(YH[i] - Y[i])
		end
	end
end

function squared_error_prime!{T<:FloatingPoint}(YH::AbstractVector{T}, Y::AbstractVector{T}, DE_DYH::AbstractVector{T})
	@inbounds begin
		for i = 1:length(Y)
			DE_DYH[i] = YH[i] - Y[i]
		end
	end
end


text:
"""
    p1 / p2

I/ Acts as joinpath.
"""
(/)(s...) = joinpath(s...)

"""
    isdef(el)
    isanydef(obj)

I/ Check if something is defined (not nothing). Check if any field of an
object is defined.
"""
isdef(el)     = (el !== nothing)
isanydef(obj) = any(isdef, (getfield(obj, f) for f ∈ fieldnames(typeof(obj))))

"""
    reset!(obj; exclude, inits...)

I/ Reset a mutable object of type `T` by creating a "fresh" object calling the
constructor `T(; inits...)` and copying the fields of that fresh object into
the from those in `exclude`
"""
function reset!(
            obj::T;
            exclude=Symbol[],
            inits...
        )::T where T

    # create a new object of the same type, assumes there is
    # a constructor that accepts empty input
    fresh = T(; inits...)
    for fn ∈ fieldnames(T)
        fn ∈ exclude && continue
        # set all fields to the field value given by the default
        # constructor (but keeps the address of the parent object)
        setfield!(obj, fn, deepcopy(getfield(fresh, fn)))
    end
    return obj
end

"""
    nvec(n, T)

I/ Create a `Vector{T}` with `n`
"""
nvec(n::Int, T)    = [T() for _ in 1:n]
nvec(n::Int, T, s) = [
    T((
        get(GLE_DEFAULTS, Symbol("$(s)_$k"), nothing)
        for k in fieldnames(T)
      )...)
    for _ in 1:n
]

const GLE_INVALID_PAT = r"missing|NaN|Inf"

"""
    csv_writer(path, z)

I/ Write data `z` to `path`.
"""
function csv_writer(
            path::String,
            z,
            hasmissing::Bool
        )::Nothing

    dirpath = splitdir(path)[1]
    isdir(dirpath) || mkpath(dirpath)
    if hasmissing
        tempio = IOBuffer()
        writedlm(tempio, z)
        # NOTE assumes it's fine to materialize the buffer with huge arrays
        # it's probably not ideal but in general should be fine, huge arrays
        # are more likely to happen with 3D objects (mesh) which are somewhat
        # less likely to have missings
        temps = String(take!(tempio))
        write(path, replace(temps, GLE_INVALID_PAT => "?"))
    else
        writedlm(path, z)
    end
    return
end


text:
function list_scripts(dir, ext=".jl")
    paths = []
    for (root, _dirs, files) in walkdir(dir)
        for name in files
            if endswith(name, ext)
                push!(paths, joinpath(root, name))
            end
        end
    end
    return paths
end

plots_loaded = false
for path in list_scripts(joinpath(dirname(@__FILE__), "src/gallery"))
    pngpath = path[1:end-length(".jl")] * ".png"
    if mtime(path) < mtime(pngpath)
        info("Skip running: $path")
        continue
    end
    if ! plots_loaded
        # Load Plots.jl only if necessary.
        @info("using Plots...")
        @time using Plots
        gr()
        plots_loaded = true
    end
    info("Plotting: $path")
    plt = @time include(path)
    plt = plot(plt, dpi=30)  # plot!(plt, dpi=30) didn't work
    savefig(plt, pngpath)
end