text:
using SHA

@testset "Server" begin

    host = "127.0.0.1"
    port = 8080
    server_url = "http://" * host * ":" * string(port)

    num_comm_rounds = 1
    fraction_clients = 0.1f0
    num_total_clients = 10

    chunksize = 256
    is_client = false

    # dummy values for testing purpose
    weights = rand(Float32, 199210)
    strategy = Fed.Server.federated_averaging
    eval_hook = (weights::Vector{Float32}) -> (0.0f0, 0.0f0)

    base_config = Fed.Config.BaseConfig(
        server_url,
        num_comm_rounds,
        fraction_clients,
        num_total_clients
    )

    function start_stop_server(config)
        central_node = Fed.Server.CentralNode(
            host,
            port,
            weights,
            strategy,
            eval_hook,
            config
        )

        # start and stop server
        task = @async Fed.Server.start(central_node)

        @test length(central_node.client_manager) == 0

        @async Base.throwto(task,  InterruptException())
    end

    @testset "Start vanilla config" begin
        config = Fed.Config.VanillaConfig(base_config)
        start_stop_server(config)
    end

    @testset "Start quantized config" begin
        config = Fed.Config.QuantizedConfig{UInt8}(base_config)
        start_stop_server(config)
    end

    @testset "Start QDiff config" begin
        config = Fed.Config.QDiffConfig{UInt8}(
            base_config,
            chunksize,
            is_client
        )
        start_stop_server(config)
    end

    @testset "Start GD config" begin
        config = Fed.Config.GDConfig{UInt8}(
            base_config,
            chunksize,
            sha1,
            0x05,
            host,
            9090,
            is_client
        )
        start_stop_server(config)
    end

end


text:
module Gimbal

using ModelingToolkit, DifferentialEquations

include("components.jl")
include("plant.jl")
include("full_system.jl")
include("uncertainty.jl")
include("optimize.jl")

end

ext: jl


text:
# ===== VSB Verifications: RBF Gaussian Plots =====

println("=== BEGINNING GaussPlotting.jl ===")

# == Imports ==
using VSB
using Plots
pyplot()
using LaTeXStrings

cd("/Users/Damyn/Documents/BYU/FLOW Lab/VSB/verifications")

println("   PLOTTING...")
Xg = -1:0.01:1
Yg1 = [VSB.RBF_gauss(X,A=1.0,deriv = 0,sigma=0.2) for X in Xg]
Yg2 = [VSB.RBF_gauss(X,A=0.33 ,deriv = 1,sigma=0.2) for X in Xg]
Yg3 = [VSB.RBF_gauss(X,A=0.1153,deriv = 2,sigma=0.2) for X in Xg]
plt = plot(Xg, Yg1, line = (:blue), label = L"\phi(x)")
plot!(Xg,Yg2, line = (:red), label = L"\nabla \phi(x)")
plot!(Xg,Yg3, line = (:green), label = L"∇^2 ϕ(x)")
savefig(plt,"VSB_Verification_RBFGauss.pdf")

println("=== END OF GaussPlotting.jl ===")
# L"\nabla^2 \phi(|\textbf{x} - \textbf{x}_i|), A = 0.5"


text:
module DashEcharts
using Dash

const resources_path = realpath(joinpath( @__DIR__, "..", "deps"))
const version = "0.0.12.1"

include("jl/''_dashecharts.jl")

function __init__()
    DashBase.register_package(
        DashBase.ResourcePkg(
            "dash_echarts",
            resources_path,
            version = version,
            [
                DashBase.Resource(
    relative_package_path = "dash_echarts.min.js",
    external_url = "https://unpkg.com/dash_echarts@0.0.12.1/dash_echarts/dash_echarts.min.js",
    dynamic = nothing,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "dash_echarts.min.js.map",
    external_url = "https://unpkg.com/dash_echarts@0.0.12.1/dash_echarts/dash_echarts.min.js.map",
    dynamic = true,
    async = nothing,
    type = :js
)
            ]
        )

    )
end
end

max_stars_repo_name: Covarians/dash-echarts

max_issues_repo_name: Covarians/dash-echarts

max_forks_repo_name: Covarians/dash-echarts


text:
using SparseArrays

# Read GSET files
function graph_from_file(filename; dimacs=false)
    file_vec = readlines(filename)
    if !dimacs
        n, m = parse.(Int, split(strip(file_vec[1]), ' '))
    else
        n, _, m = parse.(Int, split(strip(file_vec[1]), ' '))
    end

    G = spzeros(Float64, n, n)
    if !dimacs
        for line in 2:length(file_vec)
            i, j, v = parse.(Int, split(strip(file_vec[line]), ' '))
            G[i, j] = G[j, i] = v
        end
    else
        for line in 2:length(file_vec)
            i, j = parse.(Int, split(strip(file_vec[line]), ' '))
            G[i, j] = G[j, i] = 1.0
        end
    end

    return G
end

# MAXCUT cost function (true -- not used in optimization)
function cost(G, y)
    return 0.25*(sum(G) - dot(y, G*y))
end

# COSMO native solve
function solve_with_COSMO(C; settings = nothing)
    n = size(C, 1)

    sdp = COSMO.Model()
    Cvec = vec_symm(C)
    M = spzeros(Int(n*(n+1)/2), n)
    q = -ones(n)
    for i in 1:n
        M[COSMO.mat_to_svec_ind(i,i), i] = -1.0
    end
    constraints = [COSMO.Constraint(M, Cvec, COSMO.PsdConeTriangle)]
    if isnothing(settings)
        settings = COSMO.Settings(
            verbose = true,
            eps_abs = 1e-5,
            eps_rel = 1e-5,
            decompose = true,
            merge_strategy = COSMO.CliqueGraphMerge,
            max_iter = 1000,
            # rho = 1e-5,
            # alpha = 1.6,
            kkt_solver = COSMO.QdldlKKTSolver,
            # kkt_solver = with_options(COSMO.MKLPardisoKKTSolver, msg_level_on = true),
            verbose_timing = true,
        )
    end
    BLAS.set_num_threads(1)
    COSMO.assemble!(
        sdp,
        spzeros(n, n),
        q,
        constraints,
        settings = settings,
        # x0 = ones(n),
    )
    results = COSMO.optimize!(sdp)
    return results
end

# Random
unzip(a) = map(x->getfield.(a, x), fieldnames(eltype(a)))
vec_symm(X) = X[LinearAlgebra.tril(trues(size(X)))']


text:
@testset "fsm_stream_reassembler_seq.jl" begin
    @testset begin
        rs = StreamReassembler(65000)
        push_substring!(rs, "abcd", 0)
        @test rs.cur_idx == 4
        read!(rs.output) == "abcd"
        @test !eof(rs.output)

        push_substring!(rs, "efgh", 4)
        @test rs.cur_idx == 8
        read!(rs.output) == "efgh"
        @test !eof(rs.output)
    end

    @testset begin
        rs = StreamReassembler(65000)
        push_substring!(rs, "abcd", 0)
        @test rs.cur_idx == 4
        @test !eof(rs.output)

        push_substring!(rs, "efgh", 4)
        @test rs.cur_idx == 8
        read!(rs.output) == "abcdefgh"
        @test !eof(rs.output)
    end

    @testset begin
        rs = StreamReassembler(65000)
        ss = IOBuffer()
        for i in 0:99
            @test rs.cur_idx == 4i
            push_substring!(rs, "abcd", 4i)
            @test !eof(rs.output)
            write(ss, "abcd")
        end
        @test read!(rs.output) == take!(ss) |> String
        @test !eof(rs.output)
    end

    @testset begin
        rs = StreamReassembler(65000)

        for i in 0:99
            @test rs.cur_idx == 4i
            push_substring!(rs, "abcd", 4i)
            @test !eof(rs.output)
            @test read!(rs.output) == "abcd"
        end
        @test !eof(rs.output)
    end
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
module stochastic_semiclassical

export schroedinger_semiclassical, master_semiclassical

using ...bases, ...states, ...operators
using ...operators_dense, ...operators_sparse
using ...semiclassical
import ...semiclassical: recast!, State, dmaster_h_dynamic
using ...timeevolution
import ...timeevolution: integrate_stoch
import ...timeevolution.timeevolution_schroedinger: dschroedinger, dschroedinger_dynamic
using ...stochastic
import ...stochastic.stochastic_master: dmaster_stochastic, dmaster_stoch_dynamic, dlindblad

const DecayRates = Union{Vector{Float64}, Matrix{Float64}, Void}

"""
    semiclassical.schroedinger_stochastic(tspan, state0, fquantum, fclassical[; fout, ...])

Integrate time-dependent Schrödinger equation coupled to a classical system.

# Arguments
* `tspan`: Vector specifying the points of time for which the output should
        be displayed.
* `state0`: Initial semi-classical state [`semiclassical.State`](@ref).
* `fquantum`: Function `f(t, psi, u) -> H` returning the time and or state
        dependent Hamiltonian.
* `fclassical`: Function `f(t, psi, u, du)` calculating the possibly time and
        state dependent derivative of the classical equations and storing it
        in the vector `du`.
* `fstoch_quantum=nothing`: Function `f(t, psi, u) -> Hs` that returns a vector
        of operators corresponding to the stochastic terms of the Hamiltonian.
        NOTE: Either this function or `fstoch_classical` has to be defined.
* `fstoch_classical=nothing`: Function `f(t, psi, u, du)` that calculates the
        stochastic terms of the derivative `du`.
        NOTE: Either this function or `fstoch_quantum` has to be defined.
* `fout=nothing`: If given, this function `fout(t, state)` is called every time
        an output should be displayed. ATTENTION: The given state is neither
        normalized nor permanent!
* `noise_processes=0`: Number of distinct white-noise processes in the equation.
        This number has to be equal to the total number of noise operators
        returned by `fstoch_quantum`. Add 1 if `fstoch_classical` is specificed.
        If unset, the number is automatically calculated from function outputs.
        NOTE: Set this number if you want to avoid an initial calculation of
        function outputs!
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function schroedinger_semiclassical(tspan, state0::State{Ket}, fquantum::Function,
                fclassical::Function; fstoch_quantum::Union{Void, Function}=nothing,
                fstoch_classical::Union{Void, Function}=nothing,
                fout::Union{Function,Void}=nothing,
                noise_processes::Int=0,
                kwargs...)
    tspan_ = convert(Vector{Float64}, tspan)
    dschroedinger_det(t::Float64, state::State{Ket}, dstate::State{Ket}) = semiclassical.dschroedinger_dynamic(t, state, fquantum, fclassical, dstate)

    if isa(fstoch_quantum, Void) && isa(fstoch_classical, Void)
        throw(ArgumentError("No stochastic functions provided!"))
    end

    x0 = Vector{Complex128}(length(state0))
    recast!(state0, x0)
    state = copy(state0)
    dstate = copy(state0)

    if noise_processes == 0
        n = 0
        if isa(fstoch_quantum, Function)
            fs_out = fstoch_quantum(0.0, state0.quantum, state0.classical)
            n += length(fs_out)
        end
        if isa(fstoch_classical, Function)
            n += 1
        end
    else
        n = noise_processes
    end

    dschroedinger_stoch(t::Float64, state::State{Ket}, dstate::State{Ket}, index::Int) =
        dschroedinger_stochastic(t, state, fstoch_quantum, fstoch_classical, dstate, index)
    integrate_stoch(tspan_, dschroedinger_det, dschroedinger_stoch, x0, state, dstate, fout, n; kwargs...)
end

"""
    stochastic.master_semiclassical(tspan, rho0, H, Hs, J; <keyword arguments>)

Time-evolution according to a stochastic master equation.

For dense arguments the `master` function calculates the
non-hermitian Hamiltonian and then calls master_nh which is slightly faster.

# Arguments
* `tspan`: Vector specifying the points of time for which output should
        be displayed.
* `rho0`: Initial density operator. Can also be a state vector which is
        automatically converted into a density operator.
* `fquantum`: Function `f(t, rho, u) -> (H, J, Jdagger)` or
        `f(t, rho, u) -> (H, J, Jdagger, rates)` giving the deterministic
        part of the master equation.
* `fclassical`: Function `f(t, rho, u, du)` that calculates the classical
        derivatives `du`.
* `fstoch_quantum=nothing`: Function `f(t, rho, u) -> Js, Jsdagger` or
        `f(t, psi, u) -> Js, Jsdagger, rates_s` that returns the stochastic
        operator for the superoperator of the form `Js[i]*rho + rho*Jsdagger[i]`.
* `fstoch_classical=nothing`: Function `f(t, rho, u, du)` that calculates the
        stochastic terms of the derivative `du`.
* `fstoch_H=nothing`: Function `f(t, rho, u) -> Hs` providing a vector of operators
        that correspond to stochastic terms of the Hamiltonian.
* `fstoch_J=nothing`: Function `f(t, rho, u) -> (J, Jdagger)` or
        `f(t, rho, u) -> (J, Jdagger, rates)` giving a stochastic
        Lindblad term.
* `rates=nothing`: Vector or matrix specifying the coefficients (decay rates)
        for the jump operators. If nothing is specified all rates are assumed
        to be 1.
* `rates_s=nothing`: Vector or matrix specifying the coefficients (decay rates)
        for the stochastic jump operators. If nothing is specified all rates are assumed
        to be 1.
* `fout=nothing`: If given, this function `fout(t, rho)` is called every time
        an output should be displayed. ATTENTION: The given state rho is not
        permanent! It is still in use by the ode solver and therefore must not
        be changed.
* `noise_processes=0`: Number of distinct white-noise processes in the equation.
        This number has to be equal to the total number of noise operators
        returned by all stochastic functions. Add 1 for classical noise.
        If unset, the number is calculated automatically from the function outputs.
        NOTE: Set this number if you want to avoid an initial calculation of
        function outputs!
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function master_semiclassical(tspan::Vector{Float64}, rho0::State{DenseOperator},
                fquantum::Function, fclassical::Function;
                fstoch_quantum::Union{Function, Void}=nothing,
                fstoch_classical::Union{Function, Void}=nothing,
                fstoch_H::Union{Function, Void}=nothing, fstoch_J::Union{Function, Void}=nothing,
                rates::DecayRates=nothing, rates_s::DecayRates=nothing,
                fout::Union{Function,Void}=nothing,
                noise_processes::Int=0,
                kwargs...)

    tmp = copy(rho0.quantum)

    if isa(rates_s, Matrix{Float64})
        throw(ArgumentError("A matrix of stochastic rates is ambiguous! Please provide a vector of stochastic rates.
        You may want to set them as ones or use diagonaljumps."))
    end
    if isa(fstoch_quantum, Void) && isa(fstoch_classical, Void) && isa(fstoch_H, Void) && isa(fstoch_J, Void)
        throw(ArgumentError("No stochastic functions provided!"))
    end

    if noise_processes == 0
        n = 0
        if isa(fstoch_quantum, Function)
            fq_out = fstoch_quantum(0, rho0.quantum, rho0.classical)
            n += length(fq_out[1])
        end
        if isa(fstoch_classical, Function)
            n += 1
        end
        if isa(fstoch_H, Function)
            n += length(fstoch_H(0, rho0.quantum, rho0.classical))
        end
        if isa(fstoch_J, Function)
            n += length(fstoch_J(0, rho0.quantum, rho0.classical)[1])
        end
    else
        n = noise_processes
    end

    dmaster_determ(t::Float64, rho::State{DenseOperator}, drho::State{DenseOperator}) =
            dmaster_h_dynamic(t, rho, fquantum, fclassical, rates, drho, tmp)
    if isa(fstoch_H, Void) && isa(fstoch_J, Void)
        dmaster_stoch_std(t::Float64, rho::State{DenseOperator},
                        drho::State{DenseOperator}, index::Int) =
            dmaster_stoch_dynamic(t, rho, fstoch_quantum, fstoch_classical,
                        rates_s, drho, tmp, index)
        integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch_std, rho0, fout, n; kwargs...)
    else
        dmaster_stoch_gen(t::Float64, rho::State{DenseOperator},
                        drho::State{DenseOperator}, index::Int) =
            dmaster_stoch_dynamic_general(t, rho, fstoch_quantum,
                        fstoch_classical, fstoch_H, fstoch_J, rates, rates_s,
                        drho, tmp, index)
        integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch_gen, rho0, fout, n; kwargs...)
    end
end
master_semiclassical(tspan::Vector{Float64}, psi0::State{Ket}, args...; kwargs...) =
        master_semiclassical(tspan, dm(psi0), args...; kwargs...)

function dschroedinger_stochastic(t::Float64, state::State{Ket}, fstoch_quantum::Function,
            fstoch_classical::Function, dstate::State{Ket}, index::Int)
    H = fstoch_quantum(t, state.quantum, state.classical)
    if index <= length(H)
        dschroedinger(state.quantum, H[index], dstate.quantum)
    else
        fstoch_classical(t, state.quantum, state.classical, dstate.classical)
    end
end
function dschroedinger_stochastic(t::Float64, state::State{Ket}, fstoch_quantum::Function,
            fstoch_classical::Void, dstate::State{Ket}, index::Int)
    H = fstoch_quantum(t, state.quantum, state.classical)
    dschroedinger(state.quantum, H[index], dstate.quantum)
end
function dschroedinger_stochastic(t::Float64, state::State{Ket}, fstoch_quantum::Void,
            fstoch_classical::Function, dstate::State{Ket}, index::Int)
    fstoch_classical(t, state.quantum, state.classical, dstate.classical)
end

function dmaster_stoch_dynamic(t::Float64, state::State{DenseOperator}, fstoch_quantum::Function,
            fstoch_classical::Void,
            rates_s::DecayRates, dstate::State{DenseOperator}, tmp::DenseOperator,
            index::Int)
    fstoch_quantum_(t, rho) = fstoch_quantum(t, state.quantum, state.classical)
    dmaster_stoch_dynamic(t, state.quantum, fstoch_quantum_, rates_s,
                dstate.quantum, tmp, index)
end
function dmaster_stoch_dynamic(t::Float64, state::State{DenseOperator}, fstoch_quantum::Void,
            fstoch_classical::Function,
            rates_s::DecayRates, dstate::State{DenseOperator}, tmp::DenseOperator,
            index::Int)
    fstoch_classical(t, state.quantum, state.classical, dstate.classical)
end
function dmaster_stoch_dynamic(t::Float64, state::State{DenseOperator}, fstoch_quantum::Function,
            fstoch_classical::Function,
            rates_s::DecayRates, dstate::State{DenseOperator}, tmp::DenseOperator,
            index::Int)
    result = fstoch_quantum(t, state.quantum, state.classical)
    if index <= length(result[1])
        @assert 2 <= length(result) <= 3
        if length(result) == 2
            Js, Jsdagger = result
            rates_s_ = rates_s
        else
            Js, Jsdagger, rates_s_ = result
        end
        dmaster_stochastic(state.quantum, nothing, rates_s_, Js, Jsdagger, dstate.quantum, tmp, index)
    else
        fstoch_classical(t, state.quantum, state.classical, dstate.classical)
    end
end

function dmaster_stoch_dynamic_general(t::Float64, state::State{DenseOperator},
            fstoch_quantum::Union{Function, Void}, fstoch_classical::Union{Function, Void},
            fstoch_H::Function, fstoch_J::Void, rates::DecayRates, rates_s::DecayRates,
            dstate::State{DenseOperator}, tmp::DenseOperator, index::Int)
    H = fstoch_H(t, state.quantum, state.classical)
    if index <= length(H)
        operators.gemm!(-1.0im, H[index], state.quantum, 0.0, dstate.quantum)
        operators.gemm!(1.0im, state.quantum, H[index], 1.0, dstate.quantum)
    else
        dmaster_stoch_dynamic(t, state, fstoch_quantum, fstoch_classical,
                rates_s, dstate, tmp, index-length(H))
    end
end
function dmaster_stoch_dynamic_general(t::Float64, state::State{DenseOperator},
            fstoch_quantum::Union{Function, Void}, fstoch_classical::Union{Function, Void},
            fstoch_H::Void, fstoch_J::Function, rates::DecayRates, rates_s::DecayRates,
            dstate::State{DenseOperator}, tmp::DenseOperator, index::Int)
    result_J = fstoch_J(t, state.quantum, state.classical)
    if index <= length(result_J[1])
        @assert 2 <= length(result_J) <= 3
        if length(result_J) == 2
            J, Jdagger = result_J
            rates_ = rates
        else
            J, Jdagger, rates_ = result_J
        end
        dlindblad(state.quantum, rates_, J, Jdagger, dstate.quantum, tmp, index)
    else
        dmaster_stoch_dynamic(t, state, fstoch_quantum, fstoch_classical,
                rates_s, dstate, tmp, index-length(result_J[1]))
    end
end
function dmaster_stoch_dynamic_general(t::Float64, state::State{DenseOperator},
            fstoch_quantum::Union{Function, Void}, fstoch_classical::Union{Function, Void},
            fstoch_H::Function, fstoch_J::Function, rates::DecayRates, rates_s::DecayRates,
            dstate::State{DenseOperator}, tmp::DenseOperator, index::Int)
    H = fstoch_H(t, state.quantum, state.classical)
    if index <= length(H)
        operators.gemm!(-1.0im, H[index], state.quantum, 0.0, dstate.quantum)
        operators.gemm!(1.0im, state.quantum, H[index], 1.0, dstate.quantum)
    else
        dmaster_stoch_dynamic_general(t, state, fstoch_quantum, fstoch_classical,
                nothing, fstoch_J, rates, rates_s, dstate, tmp, index-length(H))
    end
end

function integrate_master_stoch(tspan, df::Function, dg::Function,
                        rho0::State{DenseOperator}, fout::Union{Void, Function},
                        n::Int;
                        kwargs...)
    tspan_ = convert(Vector{Float64}, tspan)
    x0 = Vector{Complex128}(length(rho0))
    recast!(rho0, x0)
    state = copy(rho0)
    dstate = copy(rho0)
    integrate_stoch(tspan_, df, dg, x0, state, dstate, fout, n; kwargs...)
end

function recast!(state::State, x::SubArray{Complex128, 1})
    N = length(state.quantum)
    copy!(x, 1, state.quantum.data, 1, N)
    copy!(x, N+1, state.classical, 1, length(state.classical))
    x
end
function recast!(x::SubArray{Complex128, 1}, state::State)
    N = length(state.quantum)
    copy!(state.quantum.data, 1, x, 1, N)
    copy!(state.classical, 1, x, N+1, length(state.classical))
end

end # module


text:
struct Preserve{V} <: State{V}
    value::V
end

# Preserve is the only State that can store value `nothing`
Preserve(; unit, optional, _value, _type, _...) = begin
    U = value(unit)
    V = valuetype(_type, U)
    v = _value
    #V = promote_type(V, typeof(v))
    V = optional ? Union{V,Nothing} : V
    Preserve{V}(v)
end

constructortags(::Val{:Preserve}) = (:unit, :optional)

updatetags!(d, ::Val{:Preserve}; _...) = begin
    !haskey(d, :optional) && (d[:optional] = false)
end

genvartype(v::VarInfo, ::Val{:Preserve}; V, _...) = begin
    if istag(v, :optional)
        V = @q Union{$V,Nothing}
    end
    @q Preserve{$V}
end

gendefault(v::VarInfo, ::Val{:Preserve}) = gendefaultvalue(v, parameter=true)

genupdate(v::VarInfo, ::Val{:Preserve}, ::MainStep; kw...) = nothing


text:
#####
##### Convinience functions
#####

"""
    total_extent(topology, H, Δ, L)

Returns the total extent, including halo regions, of constant-spaced
`Periodic` and `Flat` dimensions with number of halo points `H`,
constant grid spacing `Δ`, and interior extent `L`.
"""
total_extent(topology, H, Δ, L) = L + (2H - 1) * Δ

"""
    total_extent(::Type{Bounded}, H, Δ, L)

Returns the total extent of, including halo regions, of constant-spaced
`Bounded` and `Flat` dimensions with number of halo points `H`,
constant grid spacing `Δ`, and interior extent `L`.
"""
total_extent(::Type{Bounded}, H, Δ, L) = L + 2H * Δ

"""
    total_length(loc, topo, N, H=0)

Returns the total length (number of nodes), including halo points, of a field
located at `Cell` centers along a grid dimension of length `N` and with halo points `H`.
"""
total_length(loc, topo, N, H=0) = N + 2H

"""
    total_length(::Type{Face}, ::Type{Bounded}, N, H=0)

Returns the total length, including halo points, of a field located at
cell `Face`s along a grid dimension of length `N` and with halo points `H`.
"""
total_length(::Type{Face}, ::Type{Bounded}, N, H=0) = N + 1 + 2H

#####
##### << Nodes >>
#####

@inline interior_indices(loc, topo, N) = 1:N
@inline interior_indices(::Type{Face}, ::Type{Bounded}, N) = 1:N+1

@inline interior_x_indices(loc, grid) = interior_indices(loc, topology(grid, 1), grid.Nx)
@inline interior_y_indices(loc, grid) = interior_indices(loc, topology(grid, 2), grid.Ny)
@inline interior_z_indices(loc, grid) = interior_indices(loc, topology(grid, 3), grid.Nz)

@inline interior_parent_indices(loc, topo, N, H) = 1+H:N+H
@inline interior_parent_indices(::Type{Face}, ::Type{Bounded}, N, H) = 1+H:N+1+H

# Node by node
@inline xnode(::Type{Cell}, i, grid) = @inbounds grid.xC[i]
@inline xnode(::Type{Face}, i, grid) = @inbounds grid.xF[i]

@inline ynode(::Type{Cell}, j, grid) = @inbounds grid.yC[j]
@inline ynode(::Type{Face}, j, grid) = @inbounds grid.yF[j]

@inline znode(::Type{Cell}, k, grid) = @inbounds grid.zC[k]
@inline znode(::Type{Face}, k, grid) = @inbounds grid.zF[k]

# Convenience is king
@inline xC(i, grid) = xnode(Cell, i, grid)
@inline xF(i, grid) = xnode(Face, i, grid)

@inline yC(j, grid) = ynode(Cell, j, grid)
@inline yF(j, grid) = ynode(Face, j, grid)

@inline zC(k, grid) = znode(Cell, k, grid)
@inline zF(k, grid) = znode(Face, k, grid)

"""
    xnodes(loc, grid)

Returns a view over the interior `loc=Cell` or `loc=Face` nodes
on `grid` in the x-direction. For `Bounded` directions,
`Face` nodes include the boundary points.

See `znodes` for examples.
"""
xnodes(::Type{Cell}, grid; reshape=false) = reshape ? Base.reshape(view(grid.xC, 1:grid.Nx), grid.Nx, 1, 1) : view(grid.xC, 1:grid.Nx)

"""
    ynodes(loc, grid)

Returns a view over the interior `loc=Cell` or `loc=Face` nodes
on `grid` in the y-direction. For `Bounded` directions,
`Face` nodes include the boundary points.

See `znodes` for examples.
"""
ynodes(::Type{Cell}, grid; reshape=false) = reshape ? Base.reshape(view(grid.yC, 1:grid.Ny), 1, grid.Ny, 1) : view(grid.yC, 1:grid.Ny)

"""
    znodes(loc, grid)

Returns a view over the interior `loc=Cell` or `loc=Face` nodes
on `grid` in the z-direction. For `Bounded` directions,
`Face` nodes include the boundary points.

Examples
========

```jldoctest
julia> using Oceananigans, Oceananigans.Grids

julia> horz_periodic_grid = RegularCartesianGrid(size=(3, 3, 3), extent=(2π, 2π, 1),
                                                 topology=(Periodic, Periodic, Bounded));

julia> zC = znodes(Cell, horz_periodic_grid)
3-element view(OffsetArray(::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}, 0:4), 1:3) with eltype Float64:
 -0.8333333333333331
 -0.4999999999999999
 -0.16666666666666652

 julia> zF = znodes(Face, horz_periodic_grid)
4-element view(OffsetArray(::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}, 0:5), 1:4) with eltype Float64:
 -1.0
 -0.6666666666666666
 -0.33333333333333337
 -4.44089209850063e-17
```
"""
znodes(::Type{Cell}, grid; reshape=false) = reshape ? Base.reshape(view(grid.zC, 1:grid.Nz), 1, 1, grid.Nz) : view(grid.zC, 1:grid.Nz)

function xnodes(::Type{Face}, grid; reshape=false)
    xF = view(grid.xF, interior_indices(Face, topology(grid, 1), grid.Nx))
    return reshape ? Base.reshape(xF, length(xF), 1, 1) : xF
end

function ynodes(::Type{Face}, grid; reshape=false)
    yF = view(grid.yF, interior_indices(Face, topology(grid, 2), grid.Ny))
    return reshape ? Base.reshape(yF, 1, length(yF), 1) : yF
end

function znodes(::Type{Face}, grid; reshape=false)
    zF = view(grid.zF, interior_indices(Face, topology(grid, 3), grid.Nz))
    return reshape ? Base.reshape(zF, 1, 1, length(zF)) : zF
end

"""
    nodes(loc, grid; reshape=false)

Returns a 3-tuple of views over the interior nodes
at the locations in `loc` in `x, y, z`.

If `reshape=true`, the views are reshaped to 3D arrays
with non-singleton dimensions 1, 2, 3 for `x, y, z`, respectively.
These reshaped arrays can then be used in broadcast operations with 3D fields
or arrays.

See `xnodes`, `ynodes`, and `znodes`.
"""
function nodes(loc, grid::AbstractGrid; reshape=false)
    if reshape
        x, y, z = nodes(loc, grid; reshape=false)

        N = (length(x), length(y), length(z))

        x = Base.reshape(x, N[1], 1, 1)
        y = Base.reshape(y, 1, N[2], 1)
        z = Base.reshape(z, 1, 1, N[3])

        return (x, y, z)
    else
        return (xnodes(loc[1], grid),
                ynodes(loc[2], grid),
                znodes(loc[3], grid))
    end
end

#####
##### Convinience functions
#####

unpack_grid(grid) = grid.Nx, grid.Ny, grid.Nz, grid.Lx, grid.Ly, grid.Lz

#####
##### Input validation
#####

function validate_topology(topology)
    for T in topology
        if !isa(T(), AbstractTopology)
            e = "$T is not a valid topology! " *
                "Valid topologies are: Periodic, Bounded, Flat."
            throw(ArgumentError(e))
        end
    end

    return topology
end

"""Validate that an argument tuple is the right length and has elements of type `argtype`."""
function validate_tupled_argument(arg, argtype, argname)
    length(arg) == 3        || throw(ArgumentError("length($argname) must be 3."))
    all(isa.(arg, argtype)) || throw(ArgumentError("$argname=$arg must contain $argtype s."))
    all(arg .> 0)           || throw(ArgumentError("Elements of $argname=$arg must be > 0!"))
    return nothing
end

coordinate_name(i) = i == 1 ? "x" : i == 2 ? "y" : "z"

function validate_dimension_specification(i, c)
    name = coordinate_name(i)
    length(c) == 2       || throw(ArgumentError("$name length($c) must be 2."))
    all(isa.(c, Number)) || throw(ArgumentError("$name=$c should contain numbers."))
    c[2] >= c[1]         || throw(ArgumentError("$name=$c should be an increasing interval."))
    return nothing
end

function validate_regular_grid_size_and_extent(FT, size, extent, halo, x, y, z)
    validate_tupled_argument(size, Integer, "size")
    validate_tupled_argument(halo, Integer, "halo")

    # Find domain endpoints or domain extent, depending on user input:
    if !isnothing(extent) # the user has specified an extent!

        (!isnothing(x) || !isnothing(y) || !isnothing(z)) &&
            throw(ArgumentError("Cannot specify both length and x, y, z keyword arguments."))

        validate_tupled_argument(extent, Number, "extent")

        Lx, Ly, Lz = extent

        # An "oceanic" default domain:
        x = (  0, Lx)
        y = (  0, Ly)
        z = (-Lz,  0)

    else # isnothing(extent) === true implies that user has not specified a length

        (isnothing(x) || isnothing(y) || isnothing(z)) &&
            throw(ArgumentError("Must supply length or x, y, z keyword arguments."))

        for (i, c) in enumerate((x, y, z))
            validate_dimension_specification(i, c)
        end

        Lx = x[2] - x[1]
        Ly = y[2] - y[1]
        Lz = z[2] - z[1]
    end

    return FT(Lx), FT(Ly), FT(Lz), FT.(x), FT.(y), FT.(z)
end

function validate_vertically_stretched_grid_size_and_xy(FT, size, halo, x, y)
    validate_tupled_argument(size, Integer, "size")
    validate_tupled_argument(halo, Integer, "halo")

    for (i, c) in enumerate((x, y))
            validate_dimension_specification(i, c)
        end

        Lx = x[2] - x[1]
        Ly = y[2] - y[1]

    return FT(Lx), FT(Ly), FT.(x), FT.(y)
end


text:
#!/usr/bin/env julia

using Pkg
Pkg.activate(normpath(joinpath(@__DIR__, ".")))

using ColorSchemes
using Images
using Colors
using EchogramColorSchemes
using TestImages
using ImageQualityIndexes

zx = ColorScheme([colorant"#000000",
                  colorant"#0000ee",
                  colorant"#ee0000",
                  colorant"#ee00ee",
                  colorant"#00ee00",
                  colorant"#00eeee",
                  colorant"#eeee00",
                  colorant"#eeeeee",
                  colorant"#0000ff",
                  colorant"#ff0000",
                  colorant"#ff00ff",
                  colorant"#00ff00",
                  colorant"#00ffff",
                  colorant"#ffff00",
                  colorant"#ffffff" ],
                 "ZX Spectrum",
                 "The colors used on the Sinclair ZX Spectrum computer.")

spectra = ColorScheme(
    [RGB{N0f8}(r, g, b) for r in [0, 1/3 , 2/3 , 1] for g in [0, 1/3 , 2/3 , 1] for b in [0, 1/3 , 2/3 , 1]],
    "Spectra",
    "The colours used on the Sinclair ZX Spectrum SPECTRA interface.")

function show_colorscheme(cs)
    cs = Lab.(cs)
    cs = sort(cs, by=x->x.l)
    img = Array{RGB{N0f8}}(undef, 10, length(cs) *10)
    for i in 1:length(cs)
        for j in (i-1) * 10 + 1:(i-1) * 10 + 10
            for k in 1:10
                img[k,j] = cs[i]
            end
        end
    end
    return img
end

save("zx.png", show_colorscheme(zx))
save("spectra.png", show_colorscheme(spectra))

"""
    rgb_dist(c1, c2)

Return the squared euclidean distance between two colours, c1 and c2
in RGB space.
"""
function rgb_dist(c1, c2)
    rdiff = float(red(c1)) - float(red(c2))
    gdiff = float(green(c1)) - float(green(c2))
    bdiff = float(blue(c1)) - float(blue(c2))
    return rdiff * rdiff + gdiff * gdiff + bdiff * bdiff
end

"""
    nearest(c, cs, f=colordiff)

Return the nearest colour in colorscheme cs to colour c using the
given colour distance metric f (which defaults to CIEDE2000).
"""
function nearest(c, cs, f=colordiff)
    d = f.(c, cs.colors)
    return cs.colors[argmin(d)]
end

"""
    recolor(img, colorscheme::ColorScheme)

Recolour image using colours from colorscheme, returning the
recoloured image. Use the given colour distance metric f (which
defaults to CIEDE2000).
"""
function recolor(img, colorscheme::ColorScheme, f=colordiff)
    map(x -> nearest(x, colorscheme, f), img)
end

function test(name, img, cs=spectra)
    save("$(name).png", img)
    M1 =  colorfulness(img)

    recoloured =  recolor(img, cs, rgb_dist)
    M2 =  colorfulness(recoloured)
    save("$(name)1.png", recoloured)

    recoloured = recolor(img, cs)
    M3 =  colorfulness(recoloured)
    save("$(name)2.png", recoloured )

    println("Colourfulness $name $M1 $M2 $M3")
end

test("sblair", load("blair.png"))

test("smandril", testimage("mandril_color.tif"))

test("slighthouse", testimage("lighthouse.png"))

test("sbarbara", testimage("barbara_color.png"))

test("schelsea", testimage("chelsea.png"))

test("zxblair", load("blair.png"), zx)

test("zxmandril", testimage("mandril_color.tif"), zx)

test("zxlighthouse", testimage("lighthouse.png"), zx)

test("zxbarbara", testimage("barbara_color.png"), zx)

test("zxchelsea", testimage("chelsea.png"), zx)


text:
######### Stan batch program example  ###########

using StanSample, Test
ProjDir = @__DIR__

dyes ="
data {
  int BATCHES;
  int SAMPLES;
  real y[BATCHES, SAMPLES];
  // vector[SAMPLES] y[BATCHES];
}

parameters {
  real<lower=0> tau_between;
  real<lower=0> tau_within;
  real theta;
  real mu[BATCHES];
}

transformed parameters {
  real sigma_between;
  real sigma_within;
  sigma_between <- 1/sqrt(tau_between);
  sigma_within <- 1/sqrt(tau_within);
}

model {
  theta ~ normal(0.0, 1E5);
  tau_between ~ gamma(.001, .001);
  tau_within ~ gamma(.001, .001);

  mu ~ normal(theta, sigma_between);
  for (n in 1:BATCHES)
    y[n] ~ normal(mu[n], sigma_within);
}

generated quantities {
  real sigmasq_between;
  real sigmasq_within;

  sigmasq_between <- 1 / tau_between;
  sigmasq_within <- 1 / tau_within;
}
"

data = joinpath(ProjDir, "dyes.json")
sm = SampleModel("dyes", dyes);
rc = stan_sample(sm; data)

if success(rc)
  samples = read_samples(sm)

  # Fetch cmdstan summary data frame.
  df = read_summary(sm)
  @test df[df.parameters .== :theta, :mean][1] ≈ 1527.5 atol=50.0
  @test df[df.parameters .== Symbol("mu[6]"), :mean][1] ≈ 1487.3 rtol=0.1

end


text:
abstract type ReserveDirection end
abstract type ReserveUp <: ReserveDirection end
abstract type ReserveDown <: ReserveDirection end
abstract type Reserve{T <: ReserveDirection} <: Service end

# Note that src/models/services_struct_types.jl and the testset
# "Test struct type collections" in test/test_services.jl must be modified if any types are
# added here.


text:
using Test
using MPI
using CLIMA.Topologies
using CLIMA.Grids
using CLIMA.MPIStateArrays

function main()
  T = Float64
  Nhorz = 3
  Nstack = 5
  N = 4
  DA = Array

  MPI.Initialized() || MPI.Init()
  Sys.iswindows() || (isinteractive() && MPI.finalize_atexit())

  comm = MPI.COMM_WORLD
  crank = MPI.Comm_rank(comm)
  csize = MPI.Comm_size(comm)

  Rrange = T.(accumulate(+,1:Nstack+1))
  topology = StackedCubedSphereTopology(MPI.COMM_SELF, Nhorz, Rrange;
                                        boundary=(1,2))
  grid = DiscontinuousSpectralElementGrid(topology; FloatType=T,
                                          DeviceArray=DA, polynomialorder=N,
                                          meshwarp=Topologies.cubedshellwarp)

  #=
  @show elems       = topology.elems
  @show realelems   = topology.realelems
  @show ghostelems  = topology.ghostelems
  @show sendelems   = topology.sendelems
  @show elemtocoord = topology.elemtocoord
  @show elemtoelem  = topology.elemtoelem
  @show elemtoface  = topology.elemtoface
  @show elemtoordr  = topology.elemtoordr
  @show elemtobndy  = topology.elemtobndy
  @show nabrtorank  = topology.nabrtorank
  @show nabrtorecv  = topology.nabrtorecv
  @show nabrtosend  = topology.nabrtosend
  =#

  # Check xyz matches before comm
  x = @view grid.vgeo[:, Grids._x, :]
  y = @view grid.vgeo[:, Grids._y, :]
  z = @view grid.vgeo[:, Grids._z, :]

  @test x[grid.vmapM] ≈ x[grid.vmapP]
  @test y[grid.vmapM] ≈ y[grid.vmapP]
  @test z[grid.vmapM] ≈ z[grid.vmapP]

  Np = (N+1)^3
  xyz = MPIStateArray{Tuple{Np, 3}, T, DA}(topology.mpicomm,
                                           length(topology.elems),
                                           realelems=topology.realelems,
                                           ghostelems=topology.ghostelems,
                                           sendelems=topology.sendelems,
                                           nabrtorank=topology.nabrtorank,
                                           nabrtorecv=topology.nabrtorecv,
                                           nabrtosend=topology.nabrtosend)
  xyz.Q[:,:,topology.realelems] .=
        @view grid.vgeo[:, [Grids._x, Grids._y, Grids._z], topology.realelems]
  MPIStateArrays.start_ghost_exchange!(xyz)
  MPIStateArrays.finish_ghost_exchange!(xyz)

  # Check xyz matches after
  x = @view xyz.Q[:, 1, :]
  y = @view xyz.Q[:, 2, :]
  z = @view xyz.Q[:, 3, :]

  @test x[grid.vmapM] ≈ x[grid.vmapP]
  @test y[grid.vmapM] ≈ y[grid.vmapP]
  @test z[grid.vmapM] ≈ z[grid.vmapP]

  nothing
end
isinteractive() || main()

isinteractive() || MPI.Finalize()

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
# Implement some extensions to Euclidean and Squared-Euclidean distances.
for (d, D) in [(:sqeuclidean, :SqEuclidean), (:euclidean, :Euclidean)]
    @eval begin
        function ew(d::$D, x::AV{<:Real}, x′::AV{<:Real})
            return colwise($D(), reshape(x, 1, :), reshape(x′, 1, :))
        end
        function pw(d::$D, x::AV{<:Real}, x′::AV{<:Real})
            return pw($D(), reshape(x, 1, :), reshape(x′, 1, :); dims=2)
        end

        ew(::$D, x::AV{T}) where {T<:Real} = zeros(T, length(x))
        pw(::$D, x::AV{<:Real}) = pairwise($D(), x, x)

        ew(::$D, x::ColVecs{<:Real}, x′::ColVecs{<:Real}) = colwise($D(), x.X, x′.X)
        pw(::$D, x::ColVecs{<:Real}, x′::ColVecs{<:Real}) = pw($D(), x.X, x′.X; dims=2)

        ew(::$D, x::ColVecs{T}) where {T<:Real} = zeros(T, length(x))
        pw(::$D, x::ColVecs{<:Real}) = pairwise($D(), x.X; dims=2)
    end
end

@adjoint function pairwise(::Euclidean, X::AV{<:Real})
    D, back = Zygote.pullback(X->pairwise(SqEuclidean(), X), X)
    D .= sqrt.(D)
    return D, function(Δ)
        Δ = Δ ./ (2 .* D)
        Δ[diagind(Δ)] .= 0
        return (nothing, first(back(Δ)))
    end
end


text:
module CameraModels

using StaticArrays, CoordinateTransformations

export Ray, CameraModel, PixelCoordinate, Vector3, Vector2, Point3
export pixel2ray, point2pixel, canreproject, origin, direction, sensorsize, columns, rows, lookdirection, updirection

const PixelCoordinate = SVector{2, Float64}
const Vector2 = SVector{2, Float64}

const Point3 = SVector{3, Float64}
const Vector3 = SVector{3, Float64}

abstract type CameraModel end

origin3d = zeros(Point3)

struct Ray
    origin::Point3
    direction::Vector3
end

"""
    origin(ray)

Return the direction of the ray as a (Vector3)
"""
function origin end

origin(vector::Vector3) = origin3d
origin(ray::Ray) = vector.origin

"""
    direction(ray)

Return the origin of ray, typically just a zero $(Point3) for normal cameras.
"""
function direction end

direction(vector::Vector3) = vector
direction(ray::Ray) = vector.direction

"""
    columns(model::CameraModel)

Returns the width of the camera sensor.
"""
function columns end

"""
    rows(model::CameraModel)

Returns the height of the camera sensor.
"""
function rows end

"""
    sensorsize(model::CameraModel)

Return the size of the camera sensor. By default calling out to columns(model) and rows(model) to build an SVector{2}

`sensorsize(cameramodel::CameraModel) = SVector{2}(columns(cameramodel), rows(cameramodel))`
"""
function sensorsize end

sensorsize(cameramodel::CameraModel) = SVector{2}(columns(cameramodel), rows(cameramodel))

"""
    pixel2ray(cameramodel::CameraModel, pixelcoordinate::$(PixelCoordinate))

Returns the ray in space (origin + direction) corresponding to this `pixelcoordinate`.
"""
function pixel2ray end

"""
    point2pixel(camera::CameraModel, pointincamera::$(Point3))

Returns the pixel location onto which the 3D coordinate `pointincamera` is projected.
"""
function point2pixel end

"""
    lookdirection(camera::CameraModel)

Return the lookdirection of this camera model.
"""
function lookdirection end

"""
    updirection(camera::CameraModel)

Return the updirection of this camera model.
"""
function updirection end

"""
    canreproject(camera::CameraModel)

Confirms if point2pixel is implemented for this camera model.
"""
canreproject(camera::CameraModel) = true

include("Pinhole.jl")

end # module


text:
# # Results of all experiments

using Rhone2020  #hide
const R = Rhone2020  #hide

#
## Run the calculations #hide
include("../read_measurements.jl"); #hide
include("../derived_quantities.jl"); #hide

# Filled symbols: experiments picked for results in the paper, criteria..
#
# - all sensors were below the water table in the borehole (if a sensor was not in the water column it was recording atmospheric pressure)
# - discharge values at the two sensors are consistent
#
# Relative position of the CTDs, depth increasing: CTD-145, CTD-309, CTD-265 (see `raw_data.pdf` for more details)
#
#
# Hydraulic gradients are computed from pressure smoothed with moving average (time window 2 minutes for AM15 and 30 minutes for AM13).
# Uncertainties of the hydraulic gradient are not shown in the plots, they would hardly be visible.
#
# ### AM13, 8 August to 19 August
#+
# In AM13 only CTD-309 and CTD-265 were operating.
#+ #hide
R.plot_data_unselected("0808", ctd309, ctd265, ctd145, mid_surf_309, mid_309_265, mid_145_309, pick, e_p) #hide
#+ #hide
R.plot_data_unselected("0908", ctd309, ctd265, ctd145, mid_surf_309, mid_309_265, mid_145_309, pick, e_p) #hide
#+ #hide
R.plot_data_unselected("1008", ctd309, ctd265, ctd145, mid_surf_309, mid_309_265, mid_145_309, pick, e_p) #hide
#+ #hide
R.plot_data_unselected("1108", ctd309, ctd265, ctd145, mid_surf_309, mid_309_265, mid_145_309, pick, e_p) #hide
#+ #hide
R.plot_data_unselected("1208", ctd309, ctd265, ctd145, mid_surf_309, mid_309_265, mid_145_309, pick, e_p) #hide
#+ #hide
R.plot_data_unselected("1308", ctd309, ctd265, ctd145, mid_surf_309, mid_309_265, mid_145_309, pick, e_p) #hide
#+ #hide
R.plot_data_unselected("1408", ctd309, ctd265, ctd145, mid_surf_309, mid_309_265, mid_145_309, pick, e_p) #hide
#+ #hide
R.plot_data_unselected("1708", ctd309, ctd265, ctd145, mid_surf_309, mid_309_265, mid_145_309, pick, e_p) #hide
#+ #hide
R.plot_data_unselected("1808", ctd309, ctd265, ctd145, mid_surf_309, mid_309_265, mid_145_309, pick, e_p) #hide
#+ #hide
R.plot_data_unselected("1908", ctd309, ctd265, ctd145, mid_surf_309, mid_309_265, mid_145_309, pick, e_p) #hide

# ### AM15, 20 and 21 August
#+
# In theory the third CTD in AM15 allowed the characterisation of a second channel section between the upper two CTDs.
# However, the discharges between the upper two sensors did not agree to within their uncertainties.
# Either there was an inflow in between or the discharge at the uppermost sensor was erroneous, potentially because the salt tracer had not mixed across the cross-section or was not completely dissolved yet.
# In both cases our results would not be meaningful, hence we decided to focus on the measurements between the two lower sensors.
#
# On 20 August there was hardly any discharge, which was already visible in the field.
#+ #hide
R.plot_data_unselected("2008",ctd309, ctd265, ctd145, mid_surf_309, mid_309_265, mid_145_309, pick, e_p) #hide
#+ #hide
R.plot_data_unselected("2108",ctd309, ctd265, ctd145, mid_surf_309, mid_309_265, mid_145_309, pick, e_p) #hide


text:
using ParamLevelSet
using jInv.Mesh
using LinearAlgebra
using SparseArrays
using Test
plotting = false;
if plotting
	using jInvVis
	using PyPlot
end

println("~~~~~~~~~~~~~~~~Testing extended RBFs.~~~~~~~~~~~~~~~~~~")

n = [32,32,32];

Mesh = getRegularMesh([0.0;3.0;0.0;3.0;0.0;3.0],n);
midmesh = [(Mesh.domain[2]+Mesh.domain[1]/2) ; (Mesh.domain[4]+Mesh.domain[3]/2); (Mesh.domain[5]+Mesh.domain[6]/2)];
lenmesh = [(Mesh.domain[2]-Mesh.domain[1]) ; (Mesh.domain[4]-Mesh.domain[3]); (Mesh.domain[6]-Mesh.domain[5])];
numParamOfRBF = 10;
nRBF = 5
m5 = zeros(nRBF*5);
m10 = zeros(nRBF*numParamOfRBF);

for k=1:nRBF
	offset10 = (k-1)*numParamOfRBF + 1;
	m10[offset10] = randn();
	beta = max(abs(randn()),0.25);
	# println(beta)
	# m10[offset10+1] = beta; m10[offset10+4] = beta; m10[offset10+6] = beta; ## for the L version.
	m10[offset10+1] = beta^2; m10[offset10+4] = beta^2; m10[offset10+6] = beta^2; ## for the A version.
	m10[(offset10+7):(offset10+9)] = midmesh + 0.1*lenmesh.*randn(size(lenmesh));
	offset5 = (k-1)*5 + 1;
	m5[offset5] = m10[offset10];
	m5[offset5+1] = beta;
	m5[(offset5+2):(offset5+4)] = m10[(offset10+7):(offset10+9)];
end

u5, = ParamLevelSetModelFunc(Mesh,m5;computeJacobian=0,numParamOfRBF = 5);
u10, = ParamLevelSetModelFunc(Mesh,m10;computeJacobian=0,numParamOfRBF = 10);

@test norm(u5-u10) < 1e-5

n = [32,32,32];
Mesh = getRegularMesh([0.0;5.0;0.0;5.0;0.0;5.0],n);
midmesh = [(Mesh.domain[2]+Mesh.domain[1])/2 ; (Mesh.domain[4]+Mesh.domain[3])/2; (Mesh.domain[5]+Mesh.domain[6])/2];
lenmesh = [(Mesh.domain[2]-Mesh.domain[1]) ; (Mesh.domain[4]-Mesh.domain[3]); (Mesh.domain[6]-Mesh.domain[5])];
numParamOfRBF = 10;
nRBF = 5
m = zeros(nRBF*numParamOfRBF);

for k=1:nRBF
	offset = (k-1)*numParamOfRBF + 1;
	m[offset] = randn();
	A = randn(3,3); A = 0.3*A'*A + 0.5*Matrix(1.0I,3,3); #A = eye(3) + 1e-5;
	# L = chol(A)';
	# m[(offset+1):(offset+6)] = L[L.!=0.0];
	m[(offset+1):(offset+6)] = A[tril(A).!=0.0];

	# m[offset+1] = 1.0; m[offset+4] = 1.0; m[offset+6] = 1.0;
	m[(offset+7):(offset+9)] = midmesh + 0.05*lenmesh.*randn(size(lenmesh));
	#println(midmesh + 0.1*lenmesh.*randn(size(lenmesh)))
end
u, = ParamLevelSetModelFunc(Mesh,m;computeJacobian=0,numParamOfRBF = 10);
if plotting
	figure()
	u = reshape(u,tuple(Mesh.n...));
	plotModel(u)
end

sigmaH = getDefaultHeaviside();
us, = ParamLevelSetModelFunc(Mesh,m;computeJacobian=0,sigma = sigmaH,numParamOfRBF = 10);
if plotting
	figure()
	plotModel(us)
end

bf = 1;
dm = 0.01*zeros(length(m));

for k=1:nRBF
	offset = (k-1)*numParamOfRBF + 1;
	dm[offset] = 0.01*randn();
	dm[(offset+1):(offset+6)] = 0.01*randn(6);# 0.01*[randn(); randn(); 0.0; randn(); 0.0; randn()];#
	dm[(offset+7):(offset+9)] = 0.01*randn(3)
end

hh = 1.0;
# u0,II,JJ,VV = ParamLevelSetModelFunc_old(Mesh,m;computeJacobian=1,bf = bf,numParamOfRBF = 10);
# J0 = sparse(II,JJ,VV,prod(Mesh.n),length(m))

u0,JBuilder = ParamLevelSetModelFunc(Mesh,m;computeJacobian=1,bf = bf,numParamOfRBF = 10);
J0 = getSparseMatrix(JBuilder);

for k=0:8
	hhh = (0.5^k)*hh;
	ut = ParamLevelSetModelFunc(Mesh,m+ hhh*dm;computeJacobian=0,bf = bf,numParamOfRBF = 10)[1];
	println("norm(ut-u0): ",norm(ut[:]-u0[:]),", norm(ut - u0 - J0*dm): ",norm(ut[:] - u0[:] - J0*hhh*dm));
end

println("With Heaviside func");

hh = 1.0;

# u0,II,JJ,VV = ParamLevelSetModelFunc_old(Mesh,m;computeJacobian=1,sigma=sigmaH,bf = bf,numParamOfRBF = 10);
# J0 = sparse(II,JJ,VV,prod(Mesh.n),length(m));

u0,JBuilder = ParamLevelSetModelFunc(Mesh,m;computeJacobian=1,sigma=sigmaH,bf = bf,numParamOfRBF = 10);
J0 = getSparseMatrix(JBuilder);
for k=0:8
	hhh = (0.5^k)*hh;
    ut = ParamLevelSetModelFunc(Mesh,m + hhh*dm;computeJacobian=0,sigma=sigmaH,bf = bf,numParamOfRBF = 10)[1];
	println("norm(ut-u0): ",norm(ut[:]-u0[:]),", norm(ut - u0 - J0*dm): ",norm(ut[:] - u0[:] - hhh*(J0*dm)));
end


text:
function ellipsoid(x)
  res = 0.0
  cumsum = 0.0
  for xx in x
    cumsum += xx
    res += cumsum^2
  end
  res
end

function sort_population(population, mu)
 for i=1:length(pop); fitnesses(i) = pop{i}.F; end;
 [sorted_fitnesses, index] = sort(fitnesses);
 for i=1:mu; sorted_pop{i} = pop{index(i)}; end
end

function cma_es(p)
  lambda = 12
  mu = 3

end

function cma_es(p;
  trace = true,

  # Stopping criteria related
  max_seconds = 4*numdims(p), max_evals_per_dim = 1e7,
  ftol = 1e-7, xtol = 1e-10, stol = 1e-10,
  max_rounds_without_improvement = 500,

  # Starting points, will be random unless specified
  xmean = false,

  # Algorithm specific params:
  covarMatrixSampler = CholeskyCovarSampler,
  utilitiesFunc = log_utilities,

  lambda = 12,
  mu = 3,
  decompose_covar_prob = 0.4,
  xmean = false,
  sigma = 0.05*rand(1:8)*minimum(diameters(search_space(p))),
  sigmaMin = 1e-10)

  )

  N = numdims(p)
  ss = search_space(p)
  max_evals = max_evals_per_dim * N

  tau = sqrt(N)
  tau_c = N^2
  tau_sigma = sqrt(N)

  C = covarMatrixSampler(N)
  utilities = utilitiesFunc(mu, lambda)

  s = zeros(N,1)
  s_sigma = zeros(N,1)

  xbest = xmean = rand_individual(ss)   # Current best mean value.
  fbest = eval1(xbest, p)
  num_fevals = 1
  archive = BlackBoxOptim.TopListArchive(N, 10)
  add_candidate!(archive, fbest, xbest[:], num_fevals)

  fevals_last_best = num_fevals
  next_print_covar = 100
  termination_reason = "?"

  start_time = time()

  # Now lets optimize! Ensure we run at least one iteration.
  while(true)

    if (time() - start_time) > max_seconds
      termination_reason = "Exceeded time budget"
      break
    end

    if num_fevals > max_evals
      termination_reason = "Exceeded function eval budget"
      break
    end

    # Decompose only with a given probability => saves time
    if rand() <= decompose_covar_prob
      decompose!(C)
    end

    # Generate new population
    ss = multivariate_normal_sample(C, N, lambda) # std
    z = sigma * ss
    xs = repmat(xmean, 1, lambda) + z               # n*lambda

    #
    # Evaluate fitness
    fitnesses = eval_fitnesses(p, xs, lambda)
    num_fevals += lambda

    # Check if best new fitness is best ever and print some info if tracing.
    indbest = indmin(fitnesses)
    fbest_new = fitnesses[indbest]

    if fbest_new < fbest
      xbest = xs[:, indbest]
      fbest = fbest_new
      add_candidate!(archive, fbest, xbest, num_fevals)
      fevals_last_best = num_fevals

      if fitness_is_within_ftol(p, ftol, fbest)
        termination_reason = "Within ftol"
        break
      end

      if trace
        println("$(num_fevals): Best fitness = $(fbest)")
        if num_fevals > next_print_covar
          next_print_covar = num_fevals + 100
          println("covar summary: ", sumstats(C.C, (x) -> @sprintf("%.2e", x)))
          println("sigma: ", sigma)
        end
      end
    else
      if (num_fevals - fevals_last_best) > max_rounds_without_improvement * lambda
        termination_reason = "Max rounds without improvement reached"
        break
      end
    end

    # Assign weights to the best individuals according to the utilities vector.
    weights = assign_weights(lambda, fitnesses, utilities)

    # Recombination
    rw = z * weights
    rstd =
    xmean += rw

end

function cma_recombine(w, std)
  return sum(w)/length(w), sum(std)/length(std)
end

##############################################################################
# sorts population w.r.t. the individuals fitness in ascending order
function sorted_pop = SortPop(pop, mu);
 for i=1:length(pop); fitnesses(i) = pop{i}.F; end;
 [sorted_fitnesses, index] = sort(fitnesses);
 for i=1:mu; sorted_pop{i} = pop{index(i)}; end
end

############################################################################
# initialization:
############################################################################
#Cov = Matrix{Float64}(I, n, n);           # initial covariance matrix
# initializing individual population:
#Individual.y = yInit;
#Individual.w = 0;
#Individual.std = 0;
#Individual.F = fitness(Individual.y);
#for i=1:mu; ParentPop{i} = Individual; end;
#yParent = yInit;        # initial centroid parent

###############################################################################
# evolution loop of the CMA-ES
###############################################################################
while(1)
 SqrtCov = chol(Cov)';                    # "square root" of covariance matrix
 for l = 1:lambda;                        # generate lambda offspring
  OffspringIndividual.std = randn(n,1);   # line (L1a)
  OffspringIndividual.w = sigma*(SqrtCov*OffspringIndividual.std); # line (L1a)
  OffspringIndividual.y = yParent + OffspringIndividual.w;         # line (L1b)
  OffspringIndividual.F = fitness(OffspringIndividual.y);  # determine fitness (L1c)
  OffspringPop{l} = OffspringIndividual;                   # offspring complete
 end;
 ParentPop = SortPop(OffspringPop, mu);   # sort population and take mu best
 disp(ParentPop{1}.F);                    # display best fitness in population
 Recombinant = CMArecomb(ParentPop);      # (L2) perform recombination
 yParent = yParent + Recombinant.w;       # (L2) calculate new centroid parent
 s = (1-1/tau)*s + sqrt(mu/tau*(2-1/tau))*Recombinant.w/sigma;   # line (L3)
 Cov = (1-1/tau_c)*Cov + (s/tau_c)*s';                           # line (L4)
 Cov = (Cov + Cov')/2;                    # enforce symmetry of cov matrix
 s_sigma = (1-1/tau_sigma)*s_sigma + sqrt(mu/tau_sigma*(2-1/tau_sigma))*Recombinant.std; # line (L5)
 sigma = sigma*exp((s_sigma'*s_sigma - n)/(2*n*sqrt(n)));        # line (L6)
 if (sigma < sigmaMin ) break; end;       # termination condition
end
###############################################################################


text:
using BenchmarkTools
using Distances
using Printf

include("benchmarks.jl")

# BenchmarkTools stores things in a Dict so it loses ordering but we want to print the table
# in a special order. Therefore define an order here:

order = [
    :SqEuclidean,
    :Euclidean,
    :Cityblock,
    :Chebyshev,
    :Minkowski,
    :Hamming,
    :CosineDist,
    :CorrDist,
    :ChiSqDist,
    :KLDivergence,
    :RenyiDivergence,
    :RenyiDivergence,
    :RenyiDivergence,
    :RenyiDivergence,
    :JSDivergence,
    :BhattacharyyaDist,
    :HellingerDist,
    :WeightedSqEuclidean,
    :WeightedEuclidean,
    :WeightedCityblock,
    :WeightedMinkowski,
    :WeightedHamming,
    :SqMahalanobis,
    :Mahalanobis,
    :Haversine,
    :BrayCurtis,
]

BenchmarkTools.DEFAULT_PARAMETERS.seconds = 2.0 # Long enough

# Tuning
if !isfile(@__DIR__, "params.json")
    tuning = tune!(SUITE; verbose = true);
    BenchmarkTools.save("params.json", "SUITE", params(SUITE))
end
loadparams!(SUITE, BenchmarkTools.load("params.json")[1], :evals, :samples);

# Run and judge
results = run(SUITE; verbose = true)
judgement = minimum(results)

# Output the comparison table
getname(T::DataType) = T.name.name

function print_table(judgement)
    for typ in ("colwise", "pairwise")
        io = IOBuffer()
        println(io, "|  distance  |  loop  |  $typ  |  gain  |")
        println(io, "|----------- | -------| ----------| -------|")
        sorted_distances = sort(collect(judgement[typ]), by = y -> findfirst(x -> x == getname(y[1]), order))

        for (dist, result) in sorted_distances
            t_loop = BenchmarkTools.time(result["loop"])
            t_spec = BenchmarkTools.time(result["specialized"])
            print(io, "| ", getname(dist), " |")
            print(io, @sprintf("%9.6fs | %9.6fs | %7.4f |\n", t_loop / 1e9, t_spec / 1e9, (t_loop / t_spec)))
        end
        print(STDOUT, String(take!(io)))
        println()
    end
end

print_table(judgement)


text:
using Documenter
push!(LOAD_PATH, joinpath(@__DIR__, "..", "src"))
using Diary

DocMeta.setdocmeta!(Diary, :DocTestSetup, :(using Diary); recursive=true)
makedocs(
    sitename = "Diary.jl 📔",
    modules = [Diary],
    pages = [
        "Index" => "index.md",
        "introduction.md",
        "configuration.md",
        "how_it_works.md",
        "developer_reference.md",
    ],
)

deploydocs(
    repo = "github.com/dalum/Diary.jl.git",
)


text:
_debug_args(expr::LazyExpression) = expr.args
_debug_args(expr::LazyExpression{<:FunctionWrapper}) = expr.f.obj[].args

function findallocs(io::IO, x, depth = 0, argnum = nothing)
    depth > 0 && print(io, "  "^depth)
    argnum != nothing && print(io, "[$argnum]: ")
    if x isa LazyExpression || x isa Parameter
        x isa Parameter && setdirty!(x)
        allocs = @allocated x()
        print(io, "$x: ")
        color = allocs > 0 ? :light_red : :green
        printstyled(io, allocs, color=color)
        print(io, " bytes")
        println(io)
    else
        println(io, typeof(x))
    end
    if x isa LazyExpression
        for (argnum, arg) in enumerate(_debug_args(x))
            findallocs(io, arg, depth + 1, argnum)
        end
    end
end

"""
$(SIGNATURES)

Utility function that can be used to track down allocations in [`LazyExpression`](@ref)s.

# Examples

The following session shows the output of `findallocs` if the expression doesn't allocate:

```jldoctest; setup = :(using Parametron; model = Parametron.mock_model())
julia> x = [Variable(model) for i in 1 : 2];

julia> param = Parameter{Int}(() -> 2, model)
Parameter{Int64, …}(…)

julia> expr = @expression param * x
LazyExpression{FunctionWrapper{…}(LazyExpression{typeof(Parametron.Functions.scale!), …}(…))}(…)

julia> Parametron.findallocs(expr)
LazyExpression{FunctionWrapper{…}(LazyExpression{typeof(Parametron.Functions.scale!), …}(…))}(…): 0 bytes
  [1]: Array{LinearTerm{Int64},1}
  [2]: Parameter{Int64, …}(…): 0 bytes
  [3]: Array{Variable,1}
```

In this session, `param` allocates, and `findallocs` reports the allocation:

```jldoctest; setup = :(using Parametron; using LinearAlgebra; model = Parametron.mock_model())
julia> x = [Variable(model) for i in 1 : 2];

julia> param = Parameter(() -> zeros(2), model)
Parameter{Array{Float64,1}, …}(…)

julia> expr = @expression param ⋅ x
LazyExpression{FunctionWrapper{…}(LazyExpression{typeof(Parametron.Functions.vecdot!), …}(…))}(…)

julia> Parametron.findallocs(expr)
LazyExpression{FunctionWrapper{…}(LazyExpression{typeof(Parametron.Functions.vecdot!), …}(…))}(…): 0 bytes
  [1]: AffineFunction{Float64}
  [2]: Parameter{Array{Float64,1}, …}(…): 96 bytes
  [3]: Array{Variable,1}
```
"""
findallocs(x) = findallocs(stdout, x)


text:
#!/usr/bin/env julia

@testset "synapse recovery" begin

    om = sparse([0 3; 3 0]);

    tm = find_trunk_mapping(om, [1,2])
    @test tm[2] == 1
    @test tm[1] == 2

    om = sparse([3 0; 0 3]);

    tm = find_trunk_mapping(om, [1,2])
    @test tm[1] == 1
    @test tm[2] == 2

    uncorr_es = Array{Any}(2,4); corr_es = Array{Any}(2,4);

    uncorr_es[1,:] = [1, [1,2], [1,1,1], 1];
    uncorr_es[2,:] = [2, [3,4], [2,2,2], 1];

    corr_es[1,:] = [1, [1,2], [1,1,1], 1];
    corr_es[2,:] = [2, [5,6], [3,3,3], 1];

    tm1 = Dict( 2 => 2, 6 => 4 )
    tm2 = Dict( 2 => 2, 6 => 6 )
    @test synapse_recovery(uncorr_es, corr_es, tm1) == 1
    @test synapse_recovery(uncorr_es, corr_es, tm2) == (1/2)
end


text:
"""
    read_charmatrix(filename)

Reads the file `filename` into a matrix of characters.
"""
function read_charmatrix(filename::AbstractString)
    mat = mapreduce(collect, hcat, eachline(filename))::Matrix{Char}
    permutedims(mat)
end


text:
fourportstring = """
  # GHz S RI R 50.0
  1.0 1.11 1.11 1.12 1.12 1.13 1.13 1.14 1.14
  1.21 1.21 1.22 1.22 1.23 1.23 1.24 1.24
  1.31 1.31 1.32 1.32 1.33 1.33 1.34 1.34
  1.41 1.41 1.42 1.42 1.43 1.43 1.44 1.44
  2.0 2.11 2.11 2.12 2.12 2.13 2.13 2.14 2.14
  2.21 2.21 2.22 2.22 2.23 2.23 2.24 2.24
  2.31 2.31 2.32 2.32 2.33 2.33 2.34 2.34
  2.41 2.41 2.42 2.42 2.43 2.43 2.44 2.44
  """

fourportdata = TouchstoneData(
  [
    DataPoint( 1e9,  [
      1.11+1.11im 1.12+1.12im 1.13+1.13im 1.14+1.14im;
      1.21+1.21im 1.22+1.22im 1.23+1.23im 1.24+1.24im;
      1.31+1.31im 1.32+1.32im 1.33+1.33im 1.34+1.34im;
      1.41+1.41im 1.42+1.42im 1.43+1.43im 1.44+1.44im;
    ] ),
    DataPoint( 2e9,  [
      2.11+2.11im 2.12+2.12im 2.13+2.13im 2.14+2.14im;
      2.21+2.21im 2.22+2.22im 2.23+2.23im 2.24+2.24im;
      2.31+2.31im 2.32+2.32im 2.33+2.33im 2.34+2.34im;
      2.41+2.41im 2.42+2.42im 2.43+2.43im 2.44+2.44im;
    ] ),
  ],
  NoiseDataPoint[],
  Options( 1e9, :ScatteringParameters, :RealImaginary, 50.0 )
)

# String:
@test write_touchstone_string( parse_touchstone_string( fourportstring, 4 ) ) == fourportstring
@test parse_touchstone_string( write_touchstone_string( fourportdata ), 4 ) == fourportdata

# Stream:
buffer = PipeBuffer()
write_touchstone_stream( buffer, fourportdata )
@test parse_touchstone_stream( buffer, 4 ) == fourportdata
close( buffer )

# File:
filename = joinpath( tempdir(), "test.s4p" )
write_touchstone_file( filename, fourportdata )
@test parse_touchstone_file( filename ) == fourportdata


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    eltypes = (N0f8, N0f16, Float32, Float64)        # eltypes of parametric colors
    pctypes = (Gray, RGB, AGray, GrayA, ARGB, RGBA)  # parametric colors
    cctypes = (Gray24, AGray32, RGB24, ARGB32)       # non-parametric colors
    realtypes = (Float16, Float32, Float64, Int)     # types for mixed Normed/Real operations
    for R in realtypes, T in eltypes, C in pctypes
        precompile(Tuple{typeof(+),C{T},C{T}})
        precompile(Tuple{typeof(-),C{T},C{T}})
        precompile(Tuple{typeof(*),R,C{T}})
        precompile(Tuple{typeof(*),C{T},R})
        precompile(Tuple{typeof(/),C{T},R})
    end
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    isplanar,
    isouterplanar

DataStructure = Vector{Vector{Vector{Int}}}

struct DFSState
    graph::PlainGraph

    rank::Dict{Int,Int} # tree node n, dfsindex(n)
    inedge::Dict{Int,Int} # tree node n, inedge(n)
    isthick::Dict{Int,Bool} # edge e, isthick(e)

    source::Dict{Int,Int} # edge e, source node dfsindex
    loworder::Dict{Int,Vector{Int}} # tree edge e, list of outedge o in low(o) order
    cotree::Dict{Int,Int} # cotree edge e, low(e)
    treeedge::Vector{Int} # tree edges in dfs order

    function DFSState(graph)
        new(graph, Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), [])
    end
end

dfs!(state::DFSState) = dfs!(state, 1)

function dfs!(state::DFSState, u::Int)
    state.rank[u] = length(state.rank) + 1
    buckets = Dict{Int,Vector{Int}}() # low(e), edges
    lows = Dict{Int,Int}() # edge, lownode(e)
    for (inc, adj) in neighbors(state.graph, u)
        inc == get(state.inedge, u, nothing) && continue # Predecessor
        haskey(state.cotree, inc) && continue # Visited
        state.source[inc] = state.rank[u]
        if !haskey(state.rank, adj)
            # Tree node
            push!(state.treeedge, inc)
            state.inedge[adj] = inc
            low = dfs!(state, adj)
        else
            # Cotree node
            low = adj
            state.cotree[inc] = state.rank[adj]
            state.isthick[inc] = false
        end
        lows[inc] = low
        rank = state.rank[low]
        haskey(buckets, rank) || (buckets[rank] = Int[])
        push!(buckets[rank], inc)
    end
    u == 1 && return # Root
    # low(e) ordering
    inedge = state.inedge[u]
    fringes = [i for i in values(state.cotree) if i < state.rank[u]]
    state.isthick[inedge] = length(unique(fringes)) > 1
    sorted = Int[]
    for k in sort(collect(keys(buckets)))
        append!(sorted, sort(buckets[k], by=x->state.isthick[x]))
    end
    state.loworder[inedge] = sorted
    return lows[state.loworder[inedge][1]]
end

function merge!(ds1::DataStructure, ds2::DataStructure,
        cotree::Dict{Int,Int}; verbose=false)
    if isempty(ds1)
        append!(ds1, ds2)
        verbose && println("Trunk cells")
        return true
    end
    monoc(a, b) = cotree[a[end]] <= cotree[b[1]]
    branch = Int[]
    rootcells = DataStructure()
    # Collect branch cells
    for (alike, oppo) in ds2
        if !isempty(oppo)
            verbose && println("Not planer: dichromatic branch found")
            return false
        end
        if length(alike) == 1 && cotree[ds1[1][1][1]] == cotree[alike[1]]
            c = Vector{Int}[alike, []]
            verbose && println("Root cell $(c)")
            push!(rootcells, c)
        else
            append!(branch, alike)
        end
    end
    if isempty(branch)
        prepend!(ds1, rootcells)
        return true
    end
    # Merge branch into trunk
    for (i, cell) in enumerate(ds1)
        (alike, oppo) = cell
        if monoc(alike, branch)
            if isempty(oppo) || monoc(oppo, branch)
                verbose && println("Next cell")
                continue
            else
                append!(alike, branch)
                verbose && println("Merge alike")
                prepend!(ds1, rootcells)
                return true
            end
        else
            if isempty(oppo)
                # Collect higher trunk cells
                for (halike, hoppo) in ds1[i:end]
                    if !isempty(hoppo)
                        verbose && println("Not planer: dichromatic branch found")
                        return false
                    else
                        append!(oppo, halike)
                    end
                end
                newds = ds1[1:i-1]
                push!(newds, [branch, oppo])
                empty!(ds1)
                append!(ds1, newds)
                verbose && println("New opposite")
                prepend!(ds1, rootcells)
                return true
            elseif monoc(oppo, branch)
                append!(oppo, branch)
                verbose && println("Merge opposite")
                prepend!(ds1, rootcells)
                return true
            else
                verbose && println("Not planar: trichromatic")
                return false
            end
        end
    end
    verbose && println("Next cell")
    push!(ds1, Vector{Int}[branch, []])
    prepend!(ds1, rootcells)
    return true
end

function remove!(ds::DataStructure, edges)
    newds = DataStructure()
    for cell in ds
        (alike, oppo) = cell
        setdiff!(alike, edges)
        setdiff!(oppo, edges)
        if !isempty(alike)
            push!(newds, cell)
        end
    end
    empty!(ds)
    append!(ds, newds)
end

function planaritytest(graph::OrderedGraph; verbose=false)
    # Do DFS to determine treeedge, cotree, loworder, source
    state = DFSState(plaingraph(graph))
    dfs!(state)
    verbose && println(state)

    # L-R check
    ds = Dict(c => Vector{Vector{Int}}[[[c], []]] for c in keys(state.cotree))
    while !isempty(state.treeedge)
        e = pop!(state.treeedge)
        ds[e] = DataStructure()
        for i in state.loworder[e]
            verbose && println("stem: $(i), ds: $(ds[i])")
            if !merge!(ds[e], ds[i], state.cotree, verbose=verbose)
                empty!(state.treeedge)
                return false
            end
        end
        verbose && println("inedge: $(e), ds: $(ds[e])")
        finished = filter(
            x -> state.cotree[x] == state.source[e], keys(state.cotree))
        verbose && println("finished: $(finished)\n")
        remove!(ds[e], finished)
    end
    return true
end

function outerplanaritytest(graph::UndirectedGraph; verbose=false)
    # Add a node connects all other nodes
    newg = plaingraph(graph)
    newnode = addnode!(newg)
    for n in nodeset(graph)
        addedge!(newg, newnode, n)
    end
    return planaritytest(newg, verbose=verbose)
end

"""
    isplanar(graph::UndirectedGraph) -> Bool

Return whether the graph is planar.

# Reference

1. de Fraysseix, H., & Ossona de Mendez, P. (2012). Trémaux trees and planarity.
   European Journal of Combinatorics, 33(3), 279–293.
   https://doi.org/10.1016/j.ejc.2011.09.012
"""
@cache function isplanar(graph::UndirectedGraph)
    if isdefined(graph, :cache) && haskey(graph.cache, :isouterplanar)
        graph.cache[:isouterplanar] && return true
    end
    edgecount(graph) < 9 && return true
    edgecount(graph) > nodecount(graph) * 3 - 6 && return false
    return planaritytest(graph, verbose=verbose)
end

"""
    isouterplanar(graph::UndirectedGraph) -> Bool

Return whether the graph is outerplanar. The outerplanarity test is based on
a planarity test (see [`isplanar`](@ref)).
"""
@cache function isouterplanar(graph::UndirectedGraph)
    if isdefined(graph, :cache) && haskey(graph.cache, :isplanar)
        graph.cache[:isplanar] || return false
    end
    edgecount(graph) < 6 && return true
    edgecount(graph) > nodecount(graph) * 2 - 6 && return false
    return outerplanaritytest(graph, verbose=verbose)
end


text:
module ens_forecast

include("ssa.jl")
include("analog.jl")
include("da.jl")

using .SSA
using .Analog
using .DA

export forecast

using Statistics
using LinearAlgebra
using Random

using NearestNeighbors
using Distributions
using PyCall

struct Forecast_Info
    errs
    errs_uncorr
    crps
    crps_uncorr
    spread
    r_errs
    ens_errs
    ens
    x_trues
    errs_m
    r_forecasts
    stds
    errs_y_fcst
end

xskillscore = pyimport("xskillscore")
xarray = pyimport("xarray")

function is_valid(p, model, integrator, t, Δt, test_time, bounds)
    end_point = integrator(model, p, t, t + test_time, Δt)
    return all(bounds[1] .<= end_point .<= bounds[2])
end

function forecast(; E::Array{float_type, 2}, model, model_err, integrator,
                  m::Integer, Δt::float_type, window::Integer, cycles::Integer,
                  outfreq::Integer, D::Integer, k, k_r, r, tree, tree_r,
                  osc_vars=1:D, means, stds, err_pct::float_type, mp, check_bounds,
                  test_time, bounds, da, R, inflation, var_model, α) where {float_type<:AbstractFloat}
    if da
        R_inv = inv(R)
    end
    x_true = E[:, end]
    x0 = copy(x_true)

    errs = []
    errs_uncorr = []
    crps = []
    crps_uncorr = []
    spread = []
    r_errs_hist = []
    ens_errs = []
    ens = []
    x_trues = []
    errs_m = []
    r_forecasts = []
    errs_y_fcst = []

    t = 0.0
    r_forecast = nothing
    hybrid_fcsts = zeros(size(E))*NaN

    for cycle=1:cycles
        println(cycle)

        x_m = mean(E, dims=2)
        append!(errs_uncorr, sqrt(mean((x_m .- x_true)[osc_vars, :].^2)))

        if (r_forecast !== nothing)
            r_ens = vcat([find_point(r, tree, E[:, i], k, 0) for i=1:m]...)
            r_errs = sqrt.(mean((r_ens .- r_forecast).^2, dims=2))
            append!(r_errs_hist, r_errs)
            append!(ens_errs, sqrt.(mean((E .- x_true).^2, dims=1)'))
            append!(r_forecasts, r_forecast)

            if !da
                E_mp = E[:, (sortperm(r_errs[:]))[1:mp]]
                E_mp_array = xarray.DataArray(data=E_mp, dims=["dim", "member"])
                E_array = xarray.DataArray(data=E, dims=["dim", "member"])
                x_m = mean(E_mp, dims=2)
                append!(errs, sqrt(mean((x_m .- x_true)[osc_vars, :].^2)))
                if var_model !== nothing
                    append!(errs_y_fcst, sqrt(mean((mean(hybrid_fcsts, dims=2) .- x_true).^2)))
                end
                append!(crps, xskillscore.crps_ensemble(x_true, E_mp_array).values[1])
                append!(crps_uncorr, xskillscore.crps_ensemble(x_true, E_array).values[1])
                append!(ens, E)
            else
                E_array = xarray.DataArray(data=E, dims=["dim", "member"])
                E = etkf(E=E, R_inv=R_inv, inflation=inflation,
                         H=x->find_point(r, tree, x, k, 0), y=r_forecast)
                E_corr_array = xarray.DataArray(data=E, dims=["dim", "member"])
                x_m = mean(E, dims=2)
                append!(errs, sqrt(mean((x_m .- x_true)[osc_vars, :].^2)))
                append!(crps, xskillscore.crps_ensemble(x_true, E_corr_array).values[1])
                append!(crps_uncorr, xskillscore.crps_ensemble(x_true, E_array).values[1])
                append!(ens, E)
            end
            append!(x_trues, x_true)
            append!(errs_m, sqrt(mean((means .- x_true).^2)))
        end

        ens_spread = mean(std(E, dims=2))
        append!(spread, ens_spread)

        if !check_bounds
            perts = rand(MvNormal(zeros(float_type, D), diagm(0=>(err_pct*stds)[:].^2)), m)
        else
            perts = Array{float_type}(undef, D, m)
            i = 0
            while i < m
                pert = rand(MvNormal(zeros(float_type, D), diagm(0=>(err_pct*stds)[:].^2)),
                            1)
                if is_valid(x_true + pert[:], model_err, integrator, t, Δt, test_time,
                            bounds)
                    i += 1
                    perts[:, i] = pert
                end
            end
        end
        E = x_true .+ perts

        x_m = mean(E, dims=2)
        p2 = find_point(r, tree, x_m, k, 0)
        r_forecast = find_point(r, tree_r, p2, k_r, window)

        for i=1:m
            integration = integrator(model_err, E[:, i], t,
                                     t + window*outfreq*Δt, Δt, inplace=false)
            E[:, i] = integration[end, :]
            if var_model !== nothing
                if window >= var_model.k_ar
                    y_fcst = var_model.forecast(integration[end-var_model.k_ar+1:end, :],
                                                window)[end, :]
                    hybrid_fcsts[:, i] = α*E[:, i] + (1 - α)*y_fcst
                end
            end
        end

        x_true = integrator(model, x_true, t, t + window*outfreq*Δt, Δt)

        t += window*outfreq*Δt
    end

    ens_errs = reshape(ens_errs, m, :)'
    r_errs_hist = reshape(r_errs_hist, m, :)'
    ens = reshape(ens, D, m, :)
    x_trues = reshape(x_trues, D, :)
    r_forecasts = reshape(r_forecasts, length(osc_vars), :)
    return Forecast_Info(errs, errs_uncorr, crps, crps_uncorr, spread,
                         r_errs_hist, ens_errs, ens, x_trues, errs_m,
                         r_forecasts, stds, errs_y_fcst)
end

end


text:
function TransformGenerator(arg0::jint, arg1::TransformProvider, arg2::jdouble)
    return TransformGenerator((jint, TransformProvider, jdouble), arg0, arg1, arg2)
end

function generate(obj::TransformGenerator, arg0::AbsoluteDate, arg1::AbsoluteDate)
    return jcall(obj, "generate", List, (AbsoluteDate, AbsoluteDate), arg0, arg1)
end


text:
export Population

using Distributions
using Random
using QuadGK
import Random.seed!

if !@isdefined(MyInt)
  const MyInt = UInt8 # When called by a fuction in Main.CGP, may have to agree with MyInt defined in CGP.jl
end
const DIST_TYPE = Dict{Any,Float64}
const IDIST_TYPE = Dict{Int64,Float64}
const SDIST_TYPE = Dict{String,Float64}
const DIST2_TYPE = Dict{Tuple{Any,Any},Float64}
const IDIST2_TYPE = Dict{Tuple{Int64,Int64},Float64}
const SDIST2_TYPE = Dict{Tuple{String,String},Float64}
const Population = Vector
const IPopulation = Array{Int,1}
const MIPopulation = Array{MyInt,1}
const SPopulation = Array{String,1}
const FPopulation = Array{Float64,1}
const PopList = Array{Population,1}

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
### Abstract Interface

# No Time Solution : Forward to `A.u`
Base.getindex(A::AbstractNoTimeSolution,i::Int) = A.u[i]
Base.getindex(A::AbstractNoTimeSolution,I::Vararg{Int, N}) where {N} = A.u[I]
Base.getindex(A::AbstractTimeseriesSolution,I::AbstractArray{Int}) = solution_slice(A,I)
Base.setindex!(A::AbstractNoTimeSolution, v, i::Int) = (A.u[i] = v)
Base.setindex!(A::AbstractNoTimeSolution, v, I::Vararg{Int, N}) where {N} = (A.u[I] = v)
Base.size(A::AbstractNoTimeSolution) = size(A.u)

Base.summary(A::AbstractNoTimeSolution) = string(nameof(typeof(A))," with uType ",eltype(A.u))
Base.show(io::IO, A::AbstractNoTimeSolution) = (print(io,"u: ");show(io, A.u))
Base.show(io::IO, m::MIME"text/plain", A::AbstractNoTimeSolution) = (print(io,"u: ");show(io,m,A.u))

# Symbol Handling

# For handling ambiguities
for T in [Int, Colon]
    @eval Base.@propagate_inbounds Base.getindex(A::AbstractTimeseriesSolution, I::$T) = A.u[I]
end
Base.@propagate_inbounds Base.getindex(A::AbstractTimeseriesSolution, I::Union{Int,AbstractArray{Int},CartesianIndex,Colon}...) = A.u[I[end]][Base.front(I)...]
Base.@propagate_inbounds Base.getindex(A::AbstractTimeseriesSolution, i::Int,::Colon) = [A.u[j][i] for j in 1:length(A)]
Base.@propagate_inbounds Base.getindex(A::AbstractTimeseriesSolution, ::Colon,i::Int) = A.u[i]
Base.@propagate_inbounds Base.getindex(A::AbstractTimeseriesSolution, i::Int,II::AbstractArray{Int}) = [A.u[j][i] for j in II]
Base.@propagate_inbounds function Base.getindex(A::AbstractTimeseriesSolution, ii::CartesianIndex)
    ti = Tuple(ii)
    i = last(ti)
    jj = CartesianIndex(Base.front(ti))
    return A.u[i][jj]
end

Base.@propagate_inbounds function Base.getindex(A::AbstractTimeseriesSolution,sym)
  if issymbollike(sym)
      i = sym_to_index(sym,A)
  else
      i = sym
  end

  if i == nothing
      if issymbollike(i) && Symbol(i) == getindepsym(A)
        A.t
      else
        observed(A,sym,:)
      end
  else
      A[i,:]
  end
end

Base.@propagate_inbounds function Base.getindex(A::AbstractTimeseriesSolution,sym,args...)
    if issymbollike(sym)
        i = sym_to_index(sym,A)
    else
        i = sym
    end

    if i == nothing
      if issymbollike(i) && Symbol(i) == getindepsym(A)
        A.t[args...]
      else
        observed(A,sym,args...)
      end
    else
        A[i,args...]
    end
end

function observed(A::AbstractTimeseriesSolution,sym,i::Int)
  getobserved(A)(sym,A[i],A.prob.p,A.t[i])
end

function observed(A::AbstractTimeseriesSolution,sym,i::AbstractArray{Int})
  getobserved(A).((sym,),A.u[i],(A.prob.p,),A.t[i])
end

function observed(A::AbstractTimeseriesSolution,sym,i::Colon)
  getobserved(A).((sym,),A.u,(A.prob.p,),A.t)
end

## AbstractTimeseriesSolution Interface

Base.summary(A::AbstractTimeseriesSolution) = string(
                      TYPE_COLOR, nameof(typeof(A)),
                      NO_COLOR, " with uType ",
                      TYPE_COLOR, eltype(A.u),
                      NO_COLOR, " and tType ",
                      TYPE_COLOR, eltype(A.t), NO_COLOR)

function Base.show(io::IO, A::AbstractTimeseriesSolution)
  println(io,string("retcode: ",A.retcode))
  println(io,string("Interpolation: "),interp_summary(A.interp))
  print(io,"t: ")
  show(io, A.t)
  println(io)
  print(io,"u: ")
  show(io, A.u)
end
function Base.show(io::IO, m::MIME"text/plain", A::AbstractTimeseriesSolution)
  println(io,string("retcode: ",A.retcode))
  println(io,string("Interpolation: "),interp_summary(A.interp))
  print(io,"t: ")
  show(io,m,A.t)
  println(io)
  print(io,"u: ")
  show(io,m,A.u)
end
TreeViews.hastreeview(x::DESolution) = true
function TreeViews.treelabel(io::IO,x::DESolution,
                             mime::MIME"text/plain" = MIME"text/plain"())
  show(io,mime,Text(Base.summary(x)))
end

RecursiveArrayTools.tuples(sol::AbstractTimeseriesSolution) = tuple.(sol.u,sol.t)

function Base.iterate(sol::AbstractTimeseriesSolution,state=0)
  state >= length(sol) && return nothing
  state += 1
  return (solution_new_tslocation(sol,state),state)
end

DEFAULT_PLOT_FUNC(x...) = (x...,)
DEFAULT_PLOT_FUNC(x,y,z) = (x,y,z) # For v0.5.2 bug

@recipe function f(sol::AbstractTimeseriesSolution;
                   plot_analytic=false,
                   denseplot = (sol.dense ||
                                typeof(sol.prob) <: AbstractDiscreteProblem) &&
                                !(typeof(sol) <: AbstractRODESolution) &&
                                !(hasfield(typeof(sol),:interp) &&
                                  typeof(sol.interp) <: SensitivityInterpolation),
                   plotdensity = min(Int(1e5),sol.tslocation==0 ?
                                (typeof(sol.prob) <: AbstractDiscreteProblem ?
                                max(1000,100*length(sol)) :
                                max(1000,10*length(sol))) :
                                1000*sol.tslocation),
                   tspan = nothing, axis_safety = 0.1,
                   vars=nothing)

  syms = getsyms(sol)
  int_vars = interpret_vars(vars,sol,syms)
  strs = cleansyms(syms)

  tscale = get(plotattributes, :xscale, :identity)
  plot_vecs,labels = diffeq_to_arrays(sol,plot_analytic,denseplot,
                                      plotdensity,tspan,axis_safety,
                                      vars,int_vars,tscale,strs)

  tdir = sign(sol.t[end]-sol.t[1])
  xflip --> tdir < 0
  seriestype --> :path

  # Special case labels when vars = (:x,:y,:z) or (:x) or [:x,:y] ...
  if typeof(vars) <: Tuple && (issymbollike(vars[1]) && issymbollike(vars[2]))
    xguide --> issymbollike(int_vars[1][2]) ? Symbol(int_vars[1][2]) : strs[int_vars[1][2]]
    yguide --> issymbollike(int_vars[1][3]) ? Symbol(int_vars[1][3]) : strs[int_vars[1][3]]
    if length(vars) > 2
      zguide --> issymbollike(int_vars[1][4]) ? Symbol(int_vars[1][4]) : strs[int_vars[1][4]]
    end
  end

  if (!any(issymbollike,getindex.(int_vars,1)) && getindex.(int_vars,1) == zeros(length(int_vars))) ||
     (!any(issymbollike,getindex.(int_vars,2)) && getindex.(int_vars,2) == zeros(length(int_vars))) ||
     all(t->Symbol(t)==getindepsym_defaultt(sol),getindex.(int_vars,1)) || all(t->Symbol(t)==getindepsym_defaultt(sol),getindex.(int_vars,2))
    xguide --> "$(getindepsym_defaultt(sol))"
  end
  if length(int_vars[1]) >= 3 && ((!any(issymbollike,getindex.(int_vars,3)) && getindex.(int_vars,3) == zeros(length(int_vars))) ||
     all(t->Symbol(t)==getindepsym_defaultt(sol),getindex.(int_vars,3)))
    yguide --> "$(getindepsym_defaultt(sol))"
  end
  if length(int_vars[1]) >= 4 && ((!any(issymbollike,getindex.(int_vars,4)) && getindex.(int_vars,4) == zeros(length(int_vars))) ||
     all(t->Symbol(t)==getindepsym_defaultt(sol),getindex.(int_vars,4)))
    zguide --> "$(getindepsym_defaultt(sol))"
  end

  if (!any(issymbollike,getindex.(int_vars,2)) && getindex.(int_vars,2) == zeros(length(int_vars))) ||
      all(t->Symbol(t)==getindepsym_defaultt(sol),getindex.(int_vars,2))
    if tspan === nothing
      if tdir > 0
        xlims --> (sol.t[1],sol.t[end])
      else
        xlims --> (sol.t[end],sol.t[1])
      end
    else
      xlims --> (tspan[1],tspan[end])
    end
  else
    mins = minimum(sol[int_vars[1][2],:])
    maxs = maximum(sol[int_vars[1][2],:])
    for iv in int_vars
      mins = min(mins,minimum(sol[iv[2],:]))
      maxs = max(maxs,maximum(sol[iv[2],:]))
    end
    xlims --> ((1-sign(mins)*axis_safety)*mins,(1+sign(maxs)*axis_safety)*maxs)
  end

  # Analytical solutions do not save enough information to have a good idea
  # of the axis ahead of time
  # Only set axis for animations
  if sol.tslocation != 0 && !(typeof(sol) <: AbstractAnalyticalSolution)
    if all(getindex.(int_vars,1) .== DEFAULT_PLOT_FUNC)
      mins = minimum(sol[int_vars[1][3],:])
      maxs = maximum(sol[int_vars[1][3],:])
      for iv in int_vars
        mins = min(mins,minimum(sol[iv[3],:]))
        maxs = max(maxs,maximum(sol[iv[3],:]))
      end
      ylims --> ((1-sign(mins)*axis_safety)*mins,(1+sign(maxs)*axis_safety)*maxs)

      if length(int_vars[1]) >= 4
        mins = minimum(sol[int_vars[1][4],:])
        maxs = maximum(sol[int_vars[1][4],:])
        for iv in int_vars
          mins = min(mins,minimum(sol[iv[4],:]))
          maxs = max(mins,maximum(sol[iv[4],:]))
        end
        zlims --> ((1-sign(mins)*axis_safety)*mins,(1+sign(maxs)*axis_safety)*maxs)
      end
    end
  end

  label --> reshape(labels,1,length(labels))
  (plot_vecs...,)
end

function getsyms(sol)
  if has_syms(sol.prob.f)
    return sol.prob.f.syms
  else
    return keys(sol.u[1])
  end
end

function getindepsym(sol)
  if has_indepsym(sol.prob.f)
    return sol.prob.f.indepsym
  else
    return nothing
  end
end

# Only for compatibility!
function getindepsym_defaultt(sol)
  if has_indepsym(sol.prob.f)
    return sol.prob.f.indepsym
  else
    return :t
  end
end

function getobserved(sol)
  if has_syms(sol.prob.f)
    return sol.prob.f.observed
  else
    return DEFAULT_OBSERVED
  end
end

cleansyms(syms::Nothing) = nothing
cleansyms(syms::Tuple) = collect(cleansym(sym) for sym in syms)
cleansyms(syms::Vector{Symbol}) = cleansym.(syms)
cleansyms(syms::LinearIndices) = nothing
cleansyms(syms::CartesianIndices) = nothing

function cleansym(sym::Symbol)
  str = String(sym)
  # MTK generated names
  rules = ("₊"=>".", "⦗"=>"(", "⦘"=>")")
  for r in rules
    str = replace(str, r)
  end
  return str
end

sym_to_index(sym,sol::SciMLSolution) = sym_to_index(sym,getsyms(sol))
sym_to_index(sym,syms) = findfirst(isequal(Symbol(sym)),syms)
issymbollike(x) = typeof(x) <: Symbol ||
                  Symbol(parameterless_type(typeof(x))) == :Operation ||
                  Symbol(parameterless_type(typeof(x))) == :Variable ||
                  Symbol(parameterless_type(typeof(x))) == :Sym ||
                  Symbol(parameterless_type(typeof(x))) == :Num ||
                  Symbol(parameterless_type(typeof(x))) == :Term

function diffeq_to_arrays(sol,plot_analytic,denseplot,plotdensity,tspan,axis_safety,vars,int_vars,tscale,strs)
  if tspan === nothing
    if sol.tslocation == 0
      end_idx = length(sol)
    else
      end_idx = sol.tslocation
    end
    start_idx = 1
  else
    start_idx = something(findfirst(x -> x>=tspan[1], sol.t), 1)
    end_idx = something(findlast(x -> x<=tspan[end], sol.t), length(sol))
  end

  # determine type of spacing for plott
  densetspacer = if tscale in [:ln, :log10, :log2]
    (start, stop, n) -> 10.0.^range(log10(start), stop=log10(stop), length=n)
  else
    (start, stop, n) -> range(start;stop=stop,length=n)
  end

  if denseplot
    # Generate the points from the plot from dense function
    if tspan === nothing && !(typeof(sol) <: AbstractAnalyticalSolution)
      plott = collect(densetspacer(sol.t[start_idx],sol.t[end_idx],plotdensity))
    elseif typeof(sol) <: AbstractAnalyticalSolution
      tspan = sol.prob.tspan
      plott = collect(densetspacer(tspan[1],tspan[end],plotdensity))
    else
      plott = collect(densetspacer(tspan[1],tspan[end],plotdensity))
    end
    plot_timeseries = sol(plott)
    if plot_analytic
      if typeof(sol.prob.f) <: Tuple
        plot_analytic_timeseries = [sol.prob.f[1](Val{:analytic},t,sol.prob.u0) for t in plott]
      else
        plot_analytic_timeseries = [sol.prob.f(Val{:analytic},t,sol.prob.u0) for t in plott]
      end
    else
      plot_analytic_timeseries = nothing
    end
  else
    # Plot for sparse output: use the timeseries itself
    if sol.tslocation == 0
      plott = sol.t
      plot_timeseries = sol.u
      if plot_analytic
        plot_analytic_timeseries = sol.u_analytic
      else
        plot_analytic_timeseries = nothing
      end
    else
      if tspan === nothing
        plott = sol.t[start_idx:end_idx]
      else
        plott = collect(densetspacer(tspan[1],tspan[2],plotdensity))
      end

      plot_timeseries = sol.u[start_idx:end_idx]
      if plot_analytic
        plot_analytic_timeseries = sol.u_analytic[start_idx:end_idx]
      else
        plot_analytic_timeseries = nothing
      end
    end
  end

  dims = length(int_vars[1])
  for var in int_vars
    @assert length(var) == dims
  end
  # Should check that all have the same dims!
  plot_vecs,labels = solplot_vecs_and_labels(dims,int_vars,plot_timeseries,plott,sol,plot_analytic,plot_analytic_timeseries,strs)
end

function interpret_vars(vars,sol,syms)
  if vars !== nothing && syms !== nothing
    # Do syms conversion
    tmp_vars = []
    for var in vars
      if typeof(var) <: Union{Tuple,AbstractArray} #eltype(var) <: Symbol # Some kind of iterable
        tmp = []
        for x in var
          if issymbollike(x)
            found = sym_to_index(x,syms)
            push!(tmp,found == nothing && getindepsym_defaultt(sol) == x ? 0 : something(found,x))
          else
            push!(tmp,x)
          end
        end
        if typeof(var) <: Tuple
          var_int = tuple(tmp...)
        else
          var_int = tmp
        end
      elseif issymbollike(var)
        found = sym_to_index(var,syms)
        var_int = found == nothing && getindepsym_defaultt(sol) == var ? 0 : something(found,var)
      else
        var_int = var
      end
      push!(tmp_vars,var_int)
    end
    if typeof(vars) <: Tuple
      vars = tuple(tmp_vars...)
    else
      vars = tmp_vars
    end
  end

  if vars === nothing
    # Default: plot all timeseries
    if typeof(sol[1]) <: Union{Tuple,AbstractArray}
      vars = collect((DEFAULT_PLOT_FUNC,0, i) for i in plot_indices(sol[1]))
    else
      vars = [(DEFAULT_PLOT_FUNC,0, 1)]
    end
  end

  if typeof(vars) <: Integer
    vars = [(DEFAULT_PLOT_FUNC,0, vars)]
  end

  if typeof(vars) <: AbstractArray
    # If list given, its elements should be tuples, or we assume x = time
    tmp = Tuple[]
    for x in vars
      if typeof(x) <: Tuple
        if typeof(x[1]) <: Int
          push!(tmp,tuple(DEFAULT_PLOT_FUNC,x...))
        else
          push!(tmp,x)
        end
      else
        push!(tmp,(DEFAULT_PLOT_FUNC,0, x))
      end
    end
    vars = tmp
  end

  if typeof(vars) <: Tuple
    # If tuple given...
    if typeof(vars[end-1]) <: AbstractArray
      if typeof(vars[end]) <: AbstractArray
        # If both axes are lists we zip (will fail if different lengths)
        vars = collect(zip([DEFAULT_PLOT_FUNC for i in eachindex(vars[end-1])],vars[end-1], vars[end]))
      else
        # Just the x axis is a list
        vars = [(DEFAULT_PLOT_FUNC,x, vars[end]) for x in vars[end-1]]
      end
    else
      if typeof(vars[2]) <: AbstractArray
        # Just the y axis is a list
        vars = [(DEFAULT_PLOT_FUNC,vars[end-1], y) for y in vars[end]]
      else
        # Both axes are numbers
        if typeof(vars[1]) <: Int || issymbollike(vars[1])
          vars = [tuple(DEFAULT_PLOT_FUNC,vars...)]
        else
          vars = [vars]
        end
      end
    end
  end

  # Here `vars` should be a list of tuples (x, y).
  @assert(typeof(vars) <: AbstractArray)
  @assert(eltype(vars) <: Tuple)
  vars
end

function add_labels!(labels,x,dims,sol,strs)
  lys = []
  for j in 3:dims
    if !issymbollike(x[j]) && x[j] == 0
      push!(lys,"$(getindepsym_defaultt(sol)),")
    elseif issymbollike(x[j])
      push!(lys,"$(x[j]),")
    else
      if strs !== nothing
        push!(lys,"$(strs[x[j]]),")
      else
        push!(lys,"u$(x[j]),")
      end
    end
  end
  lys[end] = chop(lys[end]) # Take off the last comma
  if !issymbollike(x[2]) && x[2] == 0 && dims == 3
    # if there are no dependence in syms, then we add "(t)"
    if strs !== nothing && (typeof(x[3]) <: Int && endswith(strs[x[3]], r"(.*)")) ||
                           (issymbollike(x[3]) && endswith(string(x[3]), r"(.*)"))
      tmp_lab = "$(lys...)"
    else
      tmp_lab = "$(lys...)($(getindepsym_defaultt(sol)))"
    end
  else
    if strs !== nothing && !issymbollike(x[2]) && x[2] != 0
      tmp = strs[x[2]]
      tmp_lab = "($tmp,$(lys...))"
    else
      if !issymbollike(x[2]) && x[2] == 0
        tmp_lab = "($(getindepsym_defaultt(sol)),$(lys...))"
      elseif issymbollike(x[2])
        tmp_lab = "($(x[2]),$(lys...))"
      else
        tmp_lab = "(u$(x[2]),$(lys...))"
      end
    end
  end
  if x[1] != DEFAULT_PLOT_FUNC
    push!(labels,"$(x[1])$(tmp_lab)")
  else
    push!(labels,tmp_lab)
  end
  labels
end

function add_analytic_labels!(labels,x,dims,sol,strs)
  lys = []
  for j in 3:dims
    if x[j] == 0 && dims == 3
      push!(lys,"$(getindepsym_defaultt(sol)),")
    else
      if strs !== nothing
        push!(lys,string("True ",strs[x[j]],","))
      else
        push!(lys,"True u$(x[j]),")
      end
    end
  end
  lys[end] = lys[end][1:end-1] # Take off the last comma
  if x[2] == 0
    tmp_lab = "$(lys...)($(getindepsym_defaultt(sol)))"
  else
    if strs !== nothing
      tmp = string("True ",strs[x[2]])
      tmp_lab = "($tmp,$(lys...))"
    else
      tmp_lab = "(True u$(x[2]),$(lys...))"
    end
  end
  if x[1] != DEFAULT_PLOT_FUNC
    push!(labels,"$(x[1])$(tmp_lab)")
  else
    push!(labels,tmp_lab)
  end
end

function u_n(timeseries::AbstractArray, n::Int,sol,plott,plot_timeseries)
  # Returns the nth variable from the timeseries, t if n == 0
  if n == 0
    return plott
  elseif n == 1 && !(typeof(sol[1]) <: Union{AbstractArray,ArrayPartition})
    return timeseries
  else
    tmp = Vector{eltype(sol[1])}(undef, length(plot_timeseries))
    for j in 1:length(plot_timeseries)
      tmp[j] = plot_timeseries[j][n]
    end
    return tmp
  end
end

function u_n(timeseries::AbstractArray, sym,sol,plott,plot_timeseries)
  @assert issymbollike(sym)
  if getindepsym_defaultt(sol) == Symbol(sym)
    return plott
  else
    getobserved(sol).((sym,),eachcol(timeseries),(sol.prob.p,),plott)
  end
end

function solplot_vecs_and_labels(dims,vars,plot_timeseries,plott,sol,plot_analytic,plot_analytic_timeseries,strs)
  plot_vecs = []
  labels = String[]
  for x in vars
    tmp = []
    for j in 2:dims
      push!(tmp, u_n(plot_timeseries, x[j],sol,plott,plot_timeseries))
    end

    f = x[1]

    tmp = f.(tmp...)

    tmp = tuple((getindex.(tmp,i) for i in eachindex(tmp[1]))...)
    for i in eachindex(tmp)
      if length(plot_vecs) < i
        push!(plot_vecs,[])
      end
      push!(plot_vecs[i],tmp[i])
    end
    add_labels!(labels,x,dims,sol,strs)
  end

  if plot_analytic
    analytic_plot_vecs = []
    for x in vars
      tmp = []
      for j in 2:dims
        push!(tmp, u_n(plot_analytic_timeseries, x[j],sol,plott,plot_analytic_timeseries))
      end
      f = x[1]
      tmp = f.(tmp...)
      tmp = tuple((getindex.(tmp,i) for i in eachindex(tmp[1]))...)
      for i in eachindex(tmp)
        push!(plot_vecs[i],tmp[i])
      end
      add_analytic_labels!(labels,x,dims,sol,strs)
    end
  end
  plot_vecs = [hcat(x...) for x in plot_vecs]
  plot_vecs,labels
end

plot_indices(A::AbstractArray) = eachindex(A)
plot_indices(A::ArrayPartition) = eachindex(A)


text:
# Basic model stuff
# model types
struct AffineModelType <: ModelType end
struct TermStructureConsistentModelType <: ModelType end
struct CalibratedModelType <: ModelType end # general model type

# type aliases
const AffineModel = Model{AffineModelType}
const TermStructureConsistentModel = Model{TermStructureConsistentModelType}
const CalibratedModel = Model{CalibratedModelType}

mutable struct ModelCommon
  observe::ObserverMixin
end

ModelCommon() = ModelCommon(ObserverMixin())

get_observers(m::Model) = m.common.observe.observers


text:
module ExampleNeuralNetworks

# package code goes here
include("feedforwardunit.jl")
include("activation_functions.jl")

export
    AbstractFeedForwardUnit,
    RectifiedLinearUnit,
    SigmoidLinearUnit

end # module

max_stars_repo_name: helgith/ExampleNeuralNetworks

max_issues_repo_name: helgith/ExampleNeuralNetworks

max_forks_repo_name: helgith/ExampleNeuralNetworks


text:
# agnostic factors for internal use

"""
$(TYPEDEF)
"""
mutable struct GenericMarginal <: AbstractRelativeRoots
    Zij::Array{Float64,1}
    Cov::Array{Float64,1}
    W::Array{Float64,1}
    GenericMarginal() = new()
    GenericMarginal(a,b,c) = new(a,b,c)
end

getSample(::CalcFactor{<:GenericMarginal}) = [0]

Base.@kwdef mutable struct PackedGenericMarginal <: AbstractPackedFactor
    Zij::Array{Float64,1}
    Cov::Array{Float64,1}
    W::Array{Float64,1}
end
function convert(::Type{PackedGenericMarginal}, d::GenericMarginal)
  return PackedGenericMarginal(d.Zij, d.Cov, d.W)
end
function convert(::Type{GenericMarginal}, d::PackedGenericMarginal)
  return GenericMarginal(d.Zij, d.Cov, d.W)
end

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


text:
export physical
abstract type DoubledBasis <: AbstractHilbert end

@inline physical(h::DoubledBasis)                = h.basis
@inline nsites_physical(h::DoubledBasis)         =
    nsites(physical(h))
@inline spacedimension_physical(h::DoubledBasis) =
    spacedimension(physical(h))
@inline indexable(h::DoubledBasis) = spacedimension(h) != 0
@inline is_homogeneous(h::DoubledBasis) = is_homogeneous(physical(h))

@inline function checkbounds_doubled(σ::ADoubleState, i)
    @inbounds N = length(first(σ)) + length(last(σ))
    i>= 1 && i<= N || Base.throw_boundserror(σ, i)
    return nothing
end

# generalMethods
function flipat!(rng::AbstractRNG, σ::ADoubleState, h::DoubledBasis, i)
    @boundscheck checkbounds_doubled(σ, i)

    hp = physical(h)
    np = nsites_physical(h)

    @inbounds res = i > np ? flipat!(rng, col(σ), hp, i-np) : flipat!(rng, row(σ), hp, i)
    return res
end

function setat!(σ::ADoubleState, h::DoubledBasis, i::Int, val)
    @boundscheck checkbounds_doubled(σ, i)

    hp = physical(h)
    np = nsites_physical(h)

    @inbounds old = i > np ? setat!(row(σ), hp, i-np, val) : setat!(col(σ), hp, i, val)
    return old
end