text:
"""
Forward Learner (core)
Glaescher et al. 2010

-- TSAS is initialized to Uniform (1/nr_states)
-- Q value at goal state is set to 0.
"""
mutable struct ForwardLearner
    nr_states::Int
    nr_actions::Int
    η::Float64
    γ::Float64
    initvalue::Float64
    R::Array{Float64, 1}
    QSA::Array{Float64, 2}
    TSAS::Array{Float64, 3}
    spe::Array{Float64, 1}
    do_debug_print::Bool
end
function ForwardLearner(; nr_states = 10, nr_actions = 4, η = 0.1,
    γ = .9, init_R = 0.,
    initvalue = init_R / (1. - γ),
    R = zeros(Float64, nr_states) .+ init_R,
    QSA = zeros(Float64, nr_states, nr_actions) .+ initvalue,
    TSAS = (1.0 / nr_states)*ones(nr_states, nr_actions, nr_states),
    do_debug_print = false)

    ForwardLearner(nr_states, nr_actions, η, γ, initvalue, R, QSA, TSAS,
        [0.], do_debug_print)
end
export ForwardLearner
function do_debug_print_final(learner::ForwardLearner)
    println("Final:")
    for s in 1:size(learner.QSA,1)
        @show learner.QSA[s, :]
    end
    println("--------------")
end
""" Update: main function """
function update!(learner::ForwardLearner, state, action, next_state, reward)
    learner.spe[1] = 1. - learner.TSAS[state, action, next_state]
    tsasprime = learner.TSAS[state, action, next_state] + learner.η * learner.spe[1]
    learner.TSAS[state, action, :] .*= (1. - learner.η)
    learner.TSAS[state, action, next_state] = copy(tsasprime)
    learner.R[next_state] = reward #rewards are observed. R is deterministic
    dovalueiteration!(learner)
    if abs(sum(learner.TSAS[state, action, :]) - 1.) > 0.001
        error("Tsas does not sum to 1 in row but to $(sum(TSAS[state,action,:]))")
    end
    if learner.do_debug_print
        do_debug_print_final(learner)
    end
end
export update!

function dovalueiteration!(learner::Union{ForwardLearner,ForwardLearnerBF}; maxIter = 100)

    epsilon = 0.0001 + 0.01 * maximum(learner.R) #stop iteration if max change is smaller than epsilon
    iter = 0
    maxChange = Inf
    while maxChange>epsilon && iter < maxIter
        maxChange = 0.
        for s in 1:learner.nr_states
            if (learner.R[s] > 0.) # in our paradigm, the reward state is final. Q_fwd[s,:] is undefined
                learner.QSA[s,:] .= 0.
            else
                for a in 1:learner.nr_actions
                    qsa = 0.
                    for sPrime in 1:learner.nr_states
                        qsa += learner.TSAS[s,a,sPrime] * (learner.R[sPrime] + learner.γ * maximum(learner.QSA[sPrime,:]))
                    end
                    maxChange = max(maxChange, abs(learner.QSA[s,a] - qsa))
                    learner.QSA[s,a] = copy(qsa)
                end
            end
        end
        iter +=1
    end
end


text:
function isallreal{T<:AbstractArray}(a::T)
    all(map(x->isa(x, Real), a))
end

function islexfirst{T<:AbstractArray,U<:AbstractArray}(a::T, b::U)
    isallreal(a) && isallreal(b) || throw(DomainError())
    for i in 1:min(length(a), length(b))
        x = a[i]
        y = b[i]
        x != y || continue
        return x < y
    end
    return length(a) < length(b)
end

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
Idea from https://discourse.julialang.org/t/i-have-vector-uint8-i-need-bitvector/2286/5
"""
function make_bitvector(v::AbstractVector{T}) where T<:Union{UInt8,UInt16,UInt32,UInt64}
    siz = sizeof(v)
    bv = falses(siz<<3)
    unsafe_copyto!(reinterpret(Ptr{UInt8}, pointer(bv.chunks)), pointer(v), siz)
    bv
end
make_bitvector(v, nb::Integer) = make_bitvector(v)[1:nb]

get_uint16(arr::Vector{UInt8}, k::Integer) = reinterpret(UInt16, view(arr, k + 1:-1:k))[]
function set_uint16!(arr::Vector{UInt8}, k::Integer, val::Integer)
    v_k = view(arr, k + 1:-1:k)
    v_k .= reinterpret(UInt8, [UInt16(val)])
end
function tcp_write_byte(sockfd::Integer, buf::AbstractVector{UInt8}, idx)
    r_buf = Ref(buf, idx)
    ret = ccall(:write, Csize_t, (Cint, Ref{UInt8}, Csize_t),
                sockfd, r_buf, 1)
    LM._strerror(ret, "tcp_write_byte()")

    return ret
end

function server_main_loop(ctx::TcpContext, mb_mapping::Ptr{ModbusMapping}; verbose=false)
    header_length = modbus_get_header_length(ctx)
    while true
        local ret_code,query
        while true
            ret_code,query = modbus_receive(ctx)
            # println("ret_code = $(ret_code)   query = $(query[1:ret_code])")
            # println("query_type $(query[hdr_length + 1])")
            ret_code == 0 || break
        end
        # do not close connection on bad CRC
        if ret_code < 0 && ret_code != EMBBADCRC
            break
        end
        if query[header_length + 1] == 0x03
            verbose && println("Read holding registers")
            if get_uint16(query, header_length + 4) == UT_REGISTERS_NB_SPECIAL
                verbose && println("Set an incorrect number of values")
                set_uint16!(query, header_length + 4, UT_REGISTERS_NB_SPECIAL - 1)
            elseif get_uint16(query, header_length + 2)  == UT_REGISTERS_ADDRESS_SPECIAL
                verbose && println("Reply to this special register address by an exception")
                modbus_reply_exception(ctx, query, LM.MODBUS_EXCEPTION_SLAVE_OR_SERVER_BUSY)
            elseif get_uint16(query, header_length + 2) ==
                UT_REGISTERS_ADDRESS_INVALID_TID_OR_SLAVE
                raw_req = [0xFF, 0x03, 0x02, 0x00, 0x00]
                verbose && println("Reply with an invalid TID or slave")
                modbus_send_raw_request(ctx, raw_req)
            elseif get_uint16(query, header_length + 2) ==
                UT_REGISTERS_ADDRESS_SLEEP_500_MS
                verbose && println("Sleep 0.5 s before replying")
                sleep(0.5)
            elseif get_uint16(query, header_length + 2) ==
                UT_REGISTERS_ADDRESS_BYTE_SLEEP_5_MS
                # Test low level only available in TCP mode
                # Catch the reply and send reply byte a byte
                req = [0x00, 0x1C, 0x00, 0x00, 0x00, 0x05 ,0xFF, 0x03, 0x02, 0x00, 0x00]
                w_s = modbus_get_socket(ctx)
                # Copy TID
                req[2] = query[2]
                for i = 1:length(req)
                    sleep(0.005)
                    rc = tcp_write_byte(w_s, req, i)
                    rc > -1 || break
                end
            end
        end
        rc = modbus_reply(ctx, query, mb_mapping)
        verbose && println("rc = $(rc)")
        rc >= 0 || break
    end

    return Libc.errno()
end

function init_server_mapping()
    bits_tab = Vector{UInt8}(make_bitvector(UT_BITS_TAB, UT_BITS_NB))
    input_bits_tab = Vector{UInt8}(make_bitvector(UT_INPUT_BITS_TAB, UT_INPUT_BITS_NB))
    mb_mapping = modbus_mapping_new_start_address(
        UT_BITS_ADDRESS, UT_BITS_NB,
        UT_INPUT_BITS_ADDRESS, UT_INPUT_BITS_NB,
        UT_REGISTERS_ADDRESS, UT_REGISTERS_NB_MAX,
        UT_INPUT_REGISTERS_ADDRESS, UT_INPUT_REGISTERS_NB
    )

    for (k, x) in enumerate(UT_BITS_TAB)
        unsafe_store!(unsafe_load(mb_mapping).tab_bits, x, k)
    end
    for (k, x) in enumerate(input_bits_tab)
        unsafe_store!(unsafe_load(mb_mapping).tab_input_bits, x, k)
    end
    for (k, x) in enumerate(UT_INPUT_REGISTERS_TAB)
        unsafe_store!(unsafe_load(mb_mapping).tab_input_registers, x, k)
    end
    for (k, x) in enumerate(UT_REGISTERS_TAB)
        unsafe_store!(unsafe_load(mb_mapping).tab_registers, x, k)
    end

    return mb_mapping
end

function unit_test_server(ctx::TcpContext, mb_mapping::Ptr{ModbusMapping}; verbose=false)
    modbus_set_debug(ctx, verbose)
    srv = modbus_tcp_listen(ctx, 1)
    srv > 0 || error("$(Libc.strerror(Libc.errno()))")
    sa = modbus_tcp_accept(ctx, srv)
    sa > 0 || error("$(Libc.strerror(Libc.errno()))")
    rc = server_main_loop(ctx, mb_mapping; verbose)
    verbose && println("Quit the loop: $(Libc.strerror(rc))")
    tcp_close(sa)
    tcp_close(srv)
end

function start_unit_test_server(ctx, mb_mapping)
    t = @task unit_test_server(ctx, mb_mapping)
    t.sticky = false;
    schedule(t)
end


text:
# ==========================================
# Convenience hybrid automaton constructors
# ==========================================

import HybridSystems: HybridSystem

# hybrid automaton with one location and a self-loop TODO comparare / use OneStateAutomaton
function HybridSystem(mode::AbstractContinuousSystem, reset_map::AbstractMap)
    automaton = LightAutomaton(1)
    add_transition!(automaton, 1, 1, 1)
    return HybridSystem(automaton, [mode], [reset_map], [AutonomousSwitching()])
end

# hybrid automaton constructors with default switchings
function HybridSystem(automaton, modes, resetmaps)
    m = nmodes(automaton)
    switchings = fill(AutonomousSwitching(), m)
    return HybridSystem(automaton, modes, resetmaps, switchings)
end

function HybridSystem(; automaton, modes, resetmaps)
    return HybridSystem(automaton, modes, resetmaps)
end

#=
# hybrid automaton with only one mode TODO compare with OneStateAutomaton
struct HA1{ST, IT, TT, GT} <: AbstractHybridSystem
    mode::ST
    invariant::IT
    transition::TT
    guard::GT
end

function HA1(sys::AbstractSystem, tr::AbstractMap)
    # normalize(sys) ? for cases without invariant
    HA1(sys, stateset(sys), tr, stateset(tr))
end
=#

# TODO refactor => MathematicalSystems (?)
HybridSystems.mode(ivp::InitialValueProblem{<:HybridSystem}, i::Integer) = mode(system(ivp), i)

# ===============================================
# Hybrid systems with time-triggered transitions
# ===============================================

import MathematicalSystems: system, statedim, initial_state

# particular HA where transitions are only time-triggered
abstract type AbstractHybridAutomatonwithClockedLinearDynamics <: AbstractHybridSystem end

# type alias
const AHACLD = AbstractHybridAutomatonwithClockedLinearDynamics

# abstarct supertype for the types of switching between discrete modes that are
# handled by hybrid systems with time-triggeerd transitions
abstract type AbstractSwitching end

# deterministic switchings occur at prescribed time stamps, i.e. there is no jitter
struct DeterministicSwitching <: AbstractSwitching end

#const deterministic_switching = DeterministicSwitching()

# non-deterministic switchings may occur at any time between a prescribed interval,
# i.e. there jitter is represented by an interval
struct NonDeterministicSwitching <: AbstractSwitching end

#const nondeterministic_switching() = DeterministicSwitching()

"""
    HACLD1{T<:AbstractSystem, MT, N, J} <: AHACLD

Single-mode hybrid automaton with clocked linear dynamics.

### Fields

- `sys`       -- system
- `rmap`      -- reset map
- `Tsample`   -- sampling time
- `ζ`         -- jitter
- `switching` -- (optional) value that

### Notes

This type is parametric in:

- `T`  -- system type
- `MT` -- type of the reset map
- `N`  -- numeric type, applies to the sampling time and jitter
- `J`  -- type associated to the jitter

The type associated to the jitter, `J`, can be one of the following:

- `Missing`     -- no jitter, i.e. switchings are deterministic
- `Number`      -- symetric jitter, i.e. non-deterministic switchings occur in the
                   intervals `[Tsample - ζ, Tsample + ζ]`
- `IA.Interval` -- nonsymetric jitter, i.e. non-deterministic switchings occur in the
                   intervals `[Tsample + inf(ζ), Tsample + sup(ζ)]`; note that
                   the infimum is expected to be negative for most use cases, i.e.
                   when the jitter interval is centered at zero

The following getter functions are available:

- `initial_state`  -- initial state of the continuous mode
- `jitter`         -- return the jitter
- `reset_map`      -- return the reset map
- `sampling_time`  -- return the sampling time
- `statedim`       -- dimension of the state-space
- `system`         -- return the continuous mode
- `switching`      -- return the type of switching

Non-deterministic switching:

 tstart       Ts-ζ⁻         tend
 [-------------|-------------]

In the following, suppose that the continuous post-operator has fixed step-size
`δ > 0`. If `F` denotes the flowpipe, then

  F[1]    F[2]    F[3]    F[4]    F[5]        F[k]
[------][------][------][------][------]  ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅  [------]

 `R = array(F)` denote the array of reach-sets time-span for reach reach-set is of the form:

Similarly we compute tstart and tend for the supremum part Ts-ζ⁺
"""
struct HACLD1{T<:AbstractSystem, MT, N, J, S<:AbstractSwitching} <: AHACLD
    sys::T
    rmap::MT
    Tsample::N
    ζ::J
    switching::S
end

# getter functions
system(hs::HACLD1) = hs.sys
reset_map(hs::HACLD1) = hs.rmap
sampling_time(hs::HACLD1) = hs.Tsample
statedim(hs::HACLD1) = statedim(hs.sys)
initial_state(hs::HACLD1) = initial_state(hs.sys)
jitter(hs::HACLD1{T, MT, N, Missing}) where {T, MT, N} = TimeInterval(zero(N), zero(N))
jitter(hs::HACLD1{T, MT, N, J}) where {T, MT, N, J}  = hs.ζ
# TODO cleanup
#jitter(hs::HACLD1{T, MT, N, J}) where {T, MT, N, J<:Real} = TimeInterval(-ζ, ζ)
#jitter(hs::HACLD1{T, MT, N, J}) where {T, MT, N, J}  = _promote_tspan(hs.ζ)
switching(::HACLD1{T, MT, N, J, S}) where {T, MT, N, J, S} = S

# default constructor without jitter
function HACLD1(sys::T, rmap::MT, Tsample::N) where {T, MT, N}
    return HACLD1(sys, rmap, Tsample, missing, DeterministicSwitching())
end

# constructor when jitter is a real number
function HACLD1(sys::T, rmap::MT, Tsample::N, ζ::J) where {T, MT, N, J<:Real}
    if iszero(ζ)
        switching = DeterministicSwitching()
        ζint = TimeInterval(ζ, ζ)
    else
        switching = NonDeterministicSwitching()
        ζint = TimeInterval(-ζ, ζ)
    end
    return HACLD1(sys, rmap, Tsample, ζint, switching)
end

# should be treated separately because IA.Interval <: Number
function HACLD1(sys::T, rmap::MT, Tsample::N, ζ::J) where {T, MT, N, J<:IA.Interval}
    return HACLD1(sys, rmap, Tsample, ζ, NonDeterministicSwitching())
end

# constructor when jitter is an interval: check whether its diameter is zero or not
# we promote ζ to accept tspan-like inputs (tuples, vectors, etc.)
function HACLD1(sys::T, rmap::MT, Tsample::N, ζ::J) where {T, MT, N, J}
    ζint = _promote_tspan(ζ)
    switching = diam(ζint) > zero(N) ? NonDeterministicSwitching() : DeterministicSwitching()
    return HACLD1(sys, rmap, Tsample, ζint, switching)
end

function _check_dim(ivp::InitialValueProblem{<:HACLD1}; throw_error::Bool=true)
    S = system(ivp)
    X0 = initial_state(ivp)
    _check_dim(S, X0, throw_error=throw_error)
end


text:
module ParametersTests

using ReTest
using BenchmarkExt
using BenchmarkExt: Parameters

@testset "Parameters test" begin
    @test BenchmarkExt.DEFAULT_PARAMETERS == Parameters()

    p = Parameters(seconds = 1, gctrial = false)
    oldseconds = BenchmarkExt.DEFAULT_PARAMETERS.seconds
    oldgctrial = BenchmarkExt.DEFAULT_PARAMETERS.gctrial
    BenchmarkExt.DEFAULT_PARAMETERS.seconds = p.seconds
    BenchmarkExt.DEFAULT_PARAMETERS.gctrial = p.gctrial

    @test p == Parameters()
    @test Parameters(p; evals = 3, time_tolerance = .32) == Parameters(evals = 3, time_tolerance = .32)

    BenchmarkExt.DEFAULT_PARAMETERS.seconds = oldseconds
    BenchmarkExt.DEFAULT_PARAMETERS.gctrial = oldgctrial

    p = Parameters(seconds = 1, gctrial = false, samples = 2, evals = 2, overhead = 42,
                   gcsample = false, time_tolerance = 0.043, memory_tolerance = 0.15)
    oldseconds = BenchmarkExt.DEFAULT_PARAMETERS.seconds
    oldgctrial = BenchmarkExt.DEFAULT_PARAMETERS.gctrial
    old_time_tolerance = BenchmarkExt.DEFAULT_PARAMETERS.time_tolerance
    old_memory_tolerance = BenchmarkExt.DEFAULT_PARAMETERS.memory_tolerance
    oldsamples = BenchmarkExt.DEFAULT_PARAMETERS.samples
    oldevals = BenchmarkExt.DEFAULT_PARAMETERS.evals
    oldoverhead = BenchmarkExt.DEFAULT_PARAMETERS.overhead
    oldgcsample = BenchmarkExt.DEFAULT_PARAMETERS.gcsample
    BenchmarkExt.DEFAULT_PARAMETERS.seconds = p.seconds
    BenchmarkExt.DEFAULT_PARAMETERS.gctrial = p.gctrial
    BenchmarkExt.DEFAULT_PARAMETERS.time_tolerance = p.time_tolerance
    BenchmarkExt.DEFAULT_PARAMETERS.memory_tolerance = p.memory_tolerance
    BenchmarkExt.DEFAULT_PARAMETERS.samples = p.samples
    BenchmarkExt.DEFAULT_PARAMETERS.evals = p.evals
    BenchmarkExt.DEFAULT_PARAMETERS.overhead = p.overhead
    BenchmarkExt.DEFAULT_PARAMETERS.gcsample = p.gcsample

    @test p == Parameters()
    @test p == Parameters(p)

    BenchmarkExt.DEFAULT_PARAMETERS.seconds = oldseconds
    BenchmarkExt.DEFAULT_PARAMETERS.gctrial = oldgctrial
    BenchmarkExt.DEFAULT_PARAMETERS.time_tolerance = old_time_tolerance
    BenchmarkExt.DEFAULT_PARAMETERS.memory_tolerance = old_memory_tolerance
    BenchmarkExt.DEFAULT_PARAMETERS.samples = oldsamples
    BenchmarkExt.DEFAULT_PARAMETERS.evals = oldevals
    BenchmarkExt.DEFAULT_PARAMETERS.overhead = oldoverhead
    BenchmarkExt.DEFAULT_PARAMETERS.gcsample = oldgcsample
end

end # module


text:
struct MetaboliteWithCoefficient
    coeff::Float64
    metabolite::Metabolite
    MetaboliteWithCoefficient(c, m) = new(Float64(c), m)
end

function Base.:*(coeff, met::Metabolite)
    return MetaboliteWithCoefficient(coeff, met)
end

function Base.:+(
    m1::Union{Metabolite,MetaboliteWithCoefficient},
    m2::Union{Metabolite,MetaboliteWithCoefficient},
)
    if typeof(m1) == Metabolite
        m1 = MetaboliteWithCoefficient(1.0, m1)
    end
    if typeof(m2) == Metabolite
        m2 = MetaboliteWithCoefficient(1.0, m2)
    end
    return MetaboliteWithCoefficient[m1, m2]
end

function Base.:+(
    m1::Vector{MetaboliteWithCoefficient},
    m2::Union{Metabolite,MetaboliteWithCoefficient},
)
    if typeof(m2) == Metabolite
        m2 = MetaboliteWithCoefficient(1.0, m2)
    end
    return push!(m1, m2)
end

function _mkrxn(substrates, products)
    metdict = Dict{String,Float64}()

    if typeof(substrates) == Metabolite
        metdict[substrates.id] = get(metdict, substrates.id, 0.0) - 1.0
    elseif typeof(substrates) == MetaboliteWithCoefficient
        metdict[substrates.metabolite.id] =
            get(metdict, substrates.metabolite.id, 0.0) - 1.0 * abs(substrates.coeff)
    elseif typeof(products) == Vector{MetaboliteWithCoefficient}
        for mwc in substrates
            metdict[mwc.metabolite.id] =
                get(metdict, mwc.metabolite.id, 0.0) - 1.0 * abs(mwc.coeff)
        end
    end

    if typeof(products) == Metabolite
        metdict[products.id] = get(metdict, products.id, 0.0) + 1.0
    elseif typeof(products) == MetaboliteWithCoefficient
        metdict[products.metabolite.id] =
            get(metdict, products.metabolite.id, 0.0) + abs(products.coeff)
    elseif typeof(products) == Vector{MetaboliteWithCoefficient}
        for mwc in products
            metdict[mwc.metabolite.id] =
                get(metdict, mwc.metabolite.id, 0.0) + 1.0 * abs(mwc.coeff)
        end
    end

    return metdict
end

"""
    ⟶(
        substrates::Union{
            Nothing,
            Metabolite,
            MetaboliteWithCoefficient,
            Vector{MetaboliteWithCoefficient},
        },
        products::Union{
            Nothing,
            Metabolite,
            MetaboliteWithCoefficient,
            Vector{MetaboliteWithCoefficient}
        },
    )

Make a forward-only [`Reaction`](@ref) from `substrates` and `products`.
An equivalent alternative is `→`.
"""
function ⟶(
    substrates::Union{
        Nothing,
        Metabolite,
        MetaboliteWithCoefficient,
        Vector{MetaboliteWithCoefficient},
    },
    products::Union{
        Nothing,
        Metabolite,
        MetaboliteWithCoefficient,
        Vector{MetaboliteWithCoefficient},
    },
)
    metdict = _mkrxn(substrates, products)
    return Reaction("", metdict, :forward)
end
const → = ⟶

"""
    ⟵(
        substrates::Union{
            Nothing,
            Metabolite,
            MetaboliteWithCoefficient,
            Vector{MetaboliteWithCoefficient},
        },
        products::Union{
            Nothing,
            Metabolite,
            MetaboliteWithCoefficient,
            Vector{MetaboliteWithCoefficient}
        },
    )

Make a reverse-only [`Reaction`](@ref) from `substrates` and `products`.
An equivalent alternative is `←`.
"""
function ⟵(
    substrates::Union{
        Nothing,
        Metabolite,
        MetaboliteWithCoefficient,
        Vector{MetaboliteWithCoefficient},
    },
    products::Union{
        Nothing,
        Metabolite,
        MetaboliteWithCoefficient,
        Vector{MetaboliteWithCoefficient},
    },
)
    metdict = _mkrxn(substrates, products)
    return Reaction("", metdict, :reverse)
end
const ← = ⟵

"""
    ⟷(
        substrates::Union{
            Nothing,
            Metabolite,
            MetaboliteWithCoefficient,
            Vector{MetaboliteWithCoefficient},
        },
        products::Union{
            Nothing,
            Metabolite,
            MetaboliteWithCoefficient,
            Vector{MetaboliteWithCoefficient}
        },
    )

Make a bidirectional (reversible) [`Reaction`](@ref) from `substrates` and `products`.
An equivalent alternative is `↔`.
"""
function ⟷(
    substrates::Union{
        Nothing,
        Metabolite,
        MetaboliteWithCoefficient,
        Vector{MetaboliteWithCoefficient},
    },
    products::Union{
        Nothing,
        Metabolite,
        MetaboliteWithCoefficient,
        Vector{MetaboliteWithCoefficient},
    },
)
    metdict = _mkrxn(substrates, products)
    return Reaction("", metdict, :bidirectional)
end
const ↔ = ⟷


text:
@recipe function plot(p::CapsuleParticle)
    @series shape(p.outer)
    @series shape(p.inner)
end

# Plot a vector of particles
@recipe function plot(particles::AbstractParticles)
    for particle in particles
        @series begin
            particle
        end
    end
end

# Plot the shape of a particle
@recipe function plot(particle::AbstractParticle)
    shape(particle)
end

@recipe function plot(shape::Shape{T,2}) where T

    grid --> false
    xguide --> "x"
    yguide --> "y"
    aspect_ratio := 1.0
    label --> ""
    linecolor --> :grey

    x, y = boundary_functions(shape)

    (x, y, 0.0, 1.0)

end


text:
#= Paper
BBVI:
http://proceedings.mlr.press/v33/ranganath14.pdf

Supplementary:
http://proceedings.mlr.press/v33/ranganath14-supp.pdf

# Another useful resource:
http://www.it.uu.se/research/systems_and_control/education/2018/pml/lectures/VILectuteNotesPart3.pdf

This gradient is for the score function estimator. I don't
understand why the gradioent of the sampled qauntity (x)
is not needed. But that's how it is.
=#

# implement for normal, gamma, invgamma, dirichlet, beta, lognormal

import SpecialFunctions: digamma, lgamma

"""
gradient of log of normal density with respect to μ and log(σ)
"""
function grad_logq_normal(x, m, log_s)
  s = exp(log_s)
  z = (x - m) / s
  # Gradient wrt μ
  grad_mu = z / s
  # Gradient wrt log(σ)
  grad_log_s = -1 + z ^ 2

  return grad_mu, grad_log_s
end

"""
gradient of log of gamma density with respect to log(shape) and log(scale)
"""
function grad_logq_gamma(x, log_shape, log_scale)
  a = exp(log_shape)
  b = exp(log_scale)
  grad_log_shape = -digamma(a) * a - a * log_scale + a * log(x)
  grad_log_scale = -a + (x / b)

  return grad_log_shape, grad_log_scale
end

"""
gradient of log of beta density with respect to log(a) and log(b)
"""
function grad_logq_beta(p, log_a, log_b)
  a = exp(log_a)
  b = exp(log_b)
  grad_log_a = digamma(a + b) * a - digamma(a) * a + a * log(p)
  grad_log_b = digamma(a + b) * b - digamma(b) * b + b * log1p(-p)

  return grad_log_a, grad_log_b
end

"""
gradient of log of dirichlet density with respect to log(alpha)
"""
function grad_logq_dirichlet(p, log_alpha)
  a = exp.(log_alpha)
  grad_log_a = digamma(sum(a)) * a - digamma.(a) .* a + a .* log.(p)
  return grad_log_a
end


text:
"""
    print_tree(tree; kwargs...)
    print_tree(io::IO, tree; kwargs...)
    print_tree(f::Function, io::IO, tree; kwargs...)

Print a text representation of `tree` to the given `io` object.

# Arguments

* `f::Function` - custom implementation of [`printnode`](@ref) to use. Should have the
  signature `f(io::IO, node)`.
* `io::IO` - IO stream to write to.
* `tree` - tree to print.
* `maxdepth::Integer = 5` - truncate printing of subtrees at this depth.
* `indicate_truncation::Bool = true` - print a vertical ellipsis character beneath
  truncated nodes.
* `charset::TreeCharSet` - [`TreeCharSet`](@ref) to use to print branches.
* `printkeys::Union{Bool, Nothing}` - Whether to print keys of child nodes (using
  `pairs(children(node))`). A value of `nothing` uses [`printkeys_default`](@ref) do decide the
  behavior on a node-by-node basis.

# Examples

```jldoctest; setup = :(using AbstractTrees)
julia> tree = [1:3, "foo", [[[4, 5], 6, 7], 8]];

julia> print_tree(tree)
Vector{Any}
├─ UnitRange{Int64}
│  ├─ 1
│  ├─ 2
│  └─ 3
├─ "foo"
└─ Vector{Any}
   ├─ Vector{Any}
   │  ├─ Vector{Int64}
   │  │  ├─ 4
   │  │  └─ 5
   │  ├─ 6
   │  └─ 7
   └─ 8

julia> print_tree(tree, maxdepth=2)
Vector{Any}
├─ UnitRange{Int64}
│  ├─ 1
│  ├─ 2
│  └─ 3
├─ "foo"
└─ Vector{Any}
   ├─ Vector{Any}
   │  ⋮
   │
   └─ 8

julia> print_tree(tree, charset=AbstractTrees.ASCII_CHARSET)
Vector{Any}
+-- UnitRange{Int64}
|   +-- 1
|   +-- 2
|   \\-- 3
+-- "foo"
\\-- Vector{Any}
    +-- Vector{Any}
    |   +-- Vector{Int64}
    |   |   +-- 4
    |   |   \\-- 5
    |   +-- 6
    |   \\-- 7
    \\-- 8
```

"""
print_tree

"""
    printnode(io::IO, node)

Print a compact representation of a single node.
"""
printnode(io::IO, node) = show(IOContext(io, :compact => true, :limit => true), node)

"""
    repr_node(node; context=nothing)

Get the string representation of a node using [`printnode`](@ref). This works
analagously to `Base.repr`.

`context` is an `IO` or `IOContext` object whose attributes are used for the
I/O stream passed to `printnode`.
"""
function repr_node(node; context=nothing)
    buf = IOBuffer()
    io = context === nothing ? buf : IOContext(buf, context)
    printnode(io, node)
    return String(take!(buf))
end

const _CharArg = Union{AbstractString, Char}

"""
    TreeCharSet(mid, terminator, skip, dash, trunc, pair)

Set of characters (or strings) used to pretty-print tree branches in [`print_tree`](@ref).

# Fields

* `mid::String` - "Forked" branch segment connecting to middle children.
* `terminator::String` - Final branch segment connecting to last child.
* `skip::String` - Vertical branch segment.
* `dash::String` - Horizontal branch segmentt printed to the right of `mid` and `terminator`.
* `trunc::String` - Used to indicate the subtree has been truncated at the maximum depth.
* `pair::String` - Printed between a child node and its key.
"""
struct TreeCharSet
    mid::String
    terminator::String
    skip::String
    dash::String
    trunc::String
    pair::String

    function TreeCharSet(mid::_CharArg, terminator::_CharArg, skip::_CharArg, dash::_CharArg, trunc::_CharArg, pair::_CharArg)
        return new(String(mid), String(terminator), String(skip), String(dash), String(trunc), String(pair))
    end
end

"""
    TreeCharSet(base::TreeCharSet; fields...)

Create a new `TreeCharSet` by modifying select fields of an existing instance.
"""
function TreeCharSet(base::TreeCharSet;
                     mid = base.mid,
                     terminator = base.terminator,
                     skip = base.skip,
                     dash = base.dash,
                     trunc = base.trunc,
                     pair = base.pair,
                     )
    return TreeCharSet(mid, terminator, skip, dash, trunc, pair)
end

"""Default `charset` argument used by [`print_tree`](@ref)."""
const DEFAULT_CHARSET = TreeCharSet("├", "└", "│", "─", "⋮", " => ")
"""Charset using only ASCII characters."""
const ASCII_CHARSET = TreeCharSet("+", "\\", "|", "--", "...", " => ")

function TreeCharSet()
    Base.depwarn("The 0-argument constructor of TreeCharSet is deprecated, use AbstractTrees.DEFAULT_CHARSET instead.", :TreeCharSet)
    return DEFAULT_CHARSET
end

"""
    printkeys_default(children)::Bool

Whether a collection of children should be printed with its keys by default.

The base behavior is to print keys for all collections for which `keys()` is defined, with the
exception of `AbstractVector`s and tuples.
"""
printkeys_default(children) = applicable(keys, children)
printkeys_default(children::AbstractVector) = false
printkeys_default(children::Tuple) = false

"""
    print_child_key(io::IO, key)

Print the key for a child node.
"""
print_child_key(io::IO, key) = show(io, key)
print_child_key(io::IO, key::CartesianIndex) = show(io, Tuple(key))

function _print_tree(printnode::Function,
                     io::IO,
                     tree;
                     maxdepth::Int,
                     indicate_truncation::Bool,
                     charset::TreeCharSet,
                     printkeys::Union{Bool, Nothing},
                     roottree = tree,
                     depth::Int = 0,
                     prefix::String = "",
                     )

    # Print node representation

    # Get node representation as string
    toprint = tree != roottree && isa(treekind(roottree), IndexedTree) ? roottree[tree] : tree
    str = repr_node(toprint, context=io)

    # Copy buffer to output, prepending prefix to each line
    for (i, line) in enumerate(split(str, '\n'))
        i != 1 && print(io, prefix)
        println(io, line)
    end

    # Node children
    c = isa(treekind(roottree), IndexedTree) ? childindices(roottree, tree) : children(roottree, tree)

    # No children?
    isempty(c) && return

    # Reached max depth?
    if depth >= maxdepth
        # Print truncation char(s)
        if indicate_truncation
            println(io, prefix, charset.trunc)
            println(io, prefix)
        end

        return
    end

    # Print keys?
    this_printkeys = applicable(keys, c) && (printkeys === nothing ? printkeys_default(c) : printkeys)

    # Print children
    s = Iterators.Stateful(this_printkeys ? pairs(c) : c)

    while !isempty(s)
        child_prefix = prefix

        if this_printkeys
            child_key, child = popfirst!(s)
        else
            child = popfirst!(s)
            child_key = nothing
        end

        print(io, prefix)

        # Last child?
        if isempty(s)
            print(io, charset.terminator)
            child_prefix *= " " ^ (textwidth(charset.skip) + textwidth(charset.dash) + 1)
        else
            print(io, charset.mid)
            child_prefix *= charset.skip * " " ^ (textwidth(charset.dash) + 1)
        end

        print(io, charset.dash, ' ')

        # Print key
        if this_printkeys
            buf = IOBuffer()
            print_child_key(IOContext(buf, io), child_key)
            key_str = String(take!(buf))

            print(io, key_str, charset.pair)

            child_prefix *= " " ^ (textwidth(key_str) + textwidth(charset.pair))
        end

        _print_tree(printnode, io, child;
            maxdepth=maxdepth, indicate_truncation=indicate_truncation, charset=charset,
            printkeys=printkeys, roottree=roottree, depth=depth + 1, prefix=child_prefix)
    end
end

function print_tree(f::Function,
                    io::IO,
                    tree;
                    maxdepth::Int = 5,
                    indicate_truncation::Bool = true,
                    charset::TreeCharSet = DEFAULT_CHARSET,
                    printkeys::Union{Bool, Nothing} = nothing,
                    )
    _print_tree(f, io, tree; maxdepth=maxdepth, indicate_truncation=indicate_truncation, charset=charset, printkeys=printkeys)
end

function print_tree(f::Function, io::IO, tree, maxdepth; kwargs...)
    Base.depwarn("Passing maxdepth as a positional argument is deprecated, use as a keyword argument instead.", :print_tree)
    print_tree(f, io, tree; maxdepth=maxdepth, kwargs...)
end

print_tree(io::IO, tree, args...; kwargs...) = print_tree(printnode, io, tree, args...; kwargs...)
print_tree(tree, args...; kwargs...) = print_tree(stdout::IO, tree, args...; kwargs...)

"""
    repr_tree(tree; context=nothing, kw...)

Get the string result of calling [`print_tree`](@ref) with the supplied arguments.

The `context` argument works as it does in `Base.repr`.
"""
function repr_tree(tree; context=nothing, kw...)
    buf = IOBuffer()
    io = context === nothing ? buf : IOContext(buf, context)
    print_tree(io, tree; kw...)
    return String(take!(buf))
end


text:
using Test
using MieSeries

@testset "Test 1" begin
    x = 0.3
    l_max = 3
    m = 0
    p, dp = MieSeries.associated_legendre(x; l_max, m)
    @test p[1] ≈ 0.3
    @test p[2] ≈ -0.365
    @test p[3] ≈ -0.3825
    @test dp[1] ≈ 1
    @test dp[2] ≈ 0.9
    @test dp[3] ≈ -0.825
end

@testset "Test 2" begin
    x = -0.3
    l_max = 3
    m = 0
    p, dp = MieSeries.associated_legendre(x; l_max, m)
    @test p[1] ≈ -0.3
    @test p[2] ≈ -0.365
    @test p[3] ≈ 0.3825
    @test dp[1] ≈ 1
    @test dp[2] ≈ -0.9
    @test dp[3] ≈ -0.825
end

@testset "Test 3" begin
    x = 0.3
    l_max = 3
    m = 1
    p, dp = MieSeries.associated_legendre(x; l_max, m)
    @test p[1] ≈ -0.953939201416945649152621586023226540254623425250545753908
    @test p[2] ≈ -0.858545281275251084237359427420903886229161082725491178517
    @test p[3] ≈ 0.7869998411689801605509128084691618957100643258317002469742
    @test dp[1] ≈ 0.3144854510165755
    @test dp[2] ≈ -2.578780698335919
    @test dp[3] ≈ -4.55217690346493
end

@testset "Test 4" begin
    x = -0.3
    l_max = 3
    m = 1
    p, dp = MieSeries.associated_legendre(x; l_max, m)
    @test p[1] ≈ -0.953939201416945649152621586023226540254623425250545753908
    @test p[2] ≈ 0.858545281275251084237359427420903886229161082725491178517
    @test p[3] ≈ 0.7869998411689801605509128084691618957100643258317002469742
    @test dp[1] ≈ -0.3144854510165755
    @test dp[2] ≈ -2.578780698335919
    @test dp[3] ≈ 4.55217690346493
end

@testset "Test 5" begin
    x = 0.8
    l_max = 50
    m = 1
    p, dp = MieSeries.associated_legendre(x; l_max, m)
    @test p[48] ≈ 5.96348895739654
    @test p[49] ≈ 2.40302362973485973633385371736914932219523828838
    @test p[50] ≈ -2.201121967241330605081500626464831230994512944496
    @test dp[48] ≈ 328.9545358407495
    @test dp[49] ≈ 566.5997821783876
    @test dp[50] ≈ 584.9974550170323
end

max_stars_repo_name: Riarrieta/MieSeries

max_issues_repo_name: Riarrieta/MieSeries

max_forks_repo_name: Riarrieta/MieSeries


text:
using ..Concepts
using ..ModelFitting
import LinearAlgebra

export IndexTracker

function has_same_dimension(data)
  if length(data) == 0
    return false
  end
  first_sz = size(data[1])
  for _data in data
    if size(_data) != first_sz
      return false
    end
  end
  return true
end

mutable struct IndexTracker{T<:Any} <: AbstractTracker
  indices::Optional{Dict{T, Array{CartesianIndex{N}, 1}}} where N<:Any
  dimension::Optional{Tuple{Vararg{Int64}}}

  function IndexTracker{T}() where {T<:Any}
    return new(nothing, nothing)
  end

  function IndexTracker{T}(data::Vararg{Array{T, N}}) where {T<:Any, N<:Any}
    if !has_same_dimension(data)
      @warn("Data stream of different dimension. Won't construct.")
      throw(DimensionMismatch())
    end
    new_object = new(Dict{T,Array{CartesianIndex{N}, 1}}(), size(data[1]))
    for _data in data
      disjoint_join(new_object, _data)
    end
    return new_object
  end

end

@overload
function Base.getindex(object::IndexTracker{T}, i::T) where T<:Any
  if isnothing(object.indices)
    @warn "indices are not constructed."
    return nothing
  end
  return object.indices[i]
end

@overload
function Base.getproperty(object::IndexTracker{T}, sym::Symbol) where T<:Any
  if sym == :keys
    return collect(keys(object.indices))
  elseif sym == :dimension
    return getfield(object, :dimension)
  elseif sym == :indices
    return getfield(object, :indices)
  elseif sym == :size
    return object.dimension
  elseif sym == :dim
    return object.dimension
  end
end

@overload
function Base.size(object::IndexTracker)
  return object.dimension
end

@overload
function Concepts.provide(object::IndexTracker{T}, data::Vararg{Array{T}}) where T<:Any
    return IndexTracker{T}(data);
end

@overload
function Concepts.groupby(obj::IndexTracker{T}, list::Array{T}) where T<:Any
    # result = Dict{T, Dict{T, Array{<:CartesianIndex}}}()
    result = Dict{T, Any}()
    for a_key in collect(unique(obj.keys))
        result[a_key] = Dict{T, Array{<:CartesianIndex}}()
        for sym in list
            result[a_key][sym] = intersect(obj[a_key], obj[sym])
        end
    end
    return result
end

# @overload
# function Base.convert(::Type{Array{<:CartesianIndex}}, x::Union{Array{Int64, 1}, Array{CartesianIndex}})
#   if typeof(x) <: Array{Int64,1}
#     return [CartesianIndex{1}(_x) for _x in x]
#   end
#   return x
# end

@overload
function Concepts.type_conversion(::Type{Array{<:CartesianIndex}}, x::Union{Array{Int64, 1}, Array{CartesianIndex{2}, 1}})
  if typeof(x) <: Array{Int64,1}
    return [CartesianIndex{1}(_x) for _x in x]
  end
  return x
end

# function disjoint_partition(a::IndexTrakcer{T}, b::Array{T, N}) where {T<:Any, N<:Any}
#     if isnothing(a.dimension)
#     a.dimension = size(b)
#     elseif size(a) != size(b)
#         @show(size(a))
#         @show(size(b))
#         throw(DimensionMismatch())
#     end
#     if isnothing(a.indices)
#         a.indices = Dict{T, Array{CartesianIndex{N}, 1}}()
#     end
#     if length(intersect(unique(a.keys), unique(b))) > 0
#         @warn("New View is not disjoint from the old")
#         throw(MethodError())
#     end
#     for sym in collect(unique(b))
#         # a.indices[sym] = convert(Array{<:CartesianIndex} ,findall(x -> x == sym, b))
#         a.indices[sym] = type_conversion(Array{<:CartesianIndex}, findall(x -> x == sym, b))
#     end
# end

@overload
function Concepts.disjoint_join(a::IndexTracker{T}, b::Array{T, N}) where {T<:Any, N<:Any}
  if isnothing(a.dimension)
    a.dimension = size(b)
  elseif size(a) != size(b)
    @show(size(a))
    @show(size(b))
    throw(DimensionMismatch())
  end

  if isnothing(a.indices)
    a.indices = Dict{T, Array{CartesianIndex{N}, 1}}()
  end
  if length(intersect(unique(a.keys), unique(b))) > 0
    @warn("New View is not disjoint from the old")
    throw(MethodError())
  end
  for sym in collect(unique(b))
      # a.indices[sym] = convert(Array{<:CartesianIndex} ,findall(x -> x == sym, b))
      a.indices[sym] = type_conversion(Array{<:CartesianIndex}, findall(x -> x == sym, b))
  end
end


text:
# This file was generated, do not modify it. # hide
using PyPlot

figure(figsize=(8,6))
plot(X.LStat, y, ls="none", marker="o")
Xnew = (LStat = collect(range(extrema(X.LStat)..., length=100)),)
plot(Xnew.LStat, predict(mach_uni, Xnew))

savefig("assets/literate/ISL-lab-3-lm1.svg") # hide


text:
module TestParser

using ConfigEnv
using Test

@testset "parsing" begin

    #comment
    @test ConfigEnv.parse("# Comment") == Dict()

    #spaces without quotes
    @test begin
        p = ConfigEnv.parse("TEST=hi  the  re")
        count(c -> c == ' ', collect(p["TEST"])) == 4
    end

    #single quotes
    @test ConfigEnv.parse("TEST=''")["TEST"] == ""
    @test ConfigEnv.parse("TEST='something'")["TEST"] == "something"

    #double quotes
    @test ConfigEnv.parse("TEST=\"\"")["TEST"] == ""
    @test ConfigEnv.parse("TEST=\"something\"")["TEST"] == "something"

    #inner quotes are mantained
    @test ConfigEnv.parse("TEST=\"\"json\"\"")["TEST"] == "\"json\""
    @test ConfigEnv.parse("TEST=\"'json'\"")["TEST"] == "'json'"
    @test ConfigEnv.parse("TEST=\"\"\"")["TEST"] == "\""
    @test ConfigEnv.parse("TEST=\"'\"")["TEST"] == "'"

    #line breaks
    @test ConfigEnv.parse("TEST=\"\\n\"")["TEST"] == "" #It's null because of final trim
    @test ConfigEnv.parse("TEST=\"\\n\\n\\nsomething\"")["TEST"] == "something"
    @test ConfigEnv.parse("TEST=\"something\\nsomething\"")["TEST"] == "something\nsomething"
    @test ConfigEnv.parse("TEST=\"something\\n\\nsomething\"")["TEST"] == "something\n\nsomething"
    @test ConfigEnv.parse("TEST='\\n'")["TEST"] == "\\n"
    @test ConfigEnv.parse("TEST=\\n")["TEST"] == "\\n"

    #empty vars
    @test ConfigEnv.parse("TEST=")["TEST"] == ""

    #trim spaces with and without quotes
    @test ConfigEnv.parse("TEST='  something  '")["TEST"] == "something"
    @test ConfigEnv.parse("TEST=\"  something  \"")["TEST"] == "something"
    @test ConfigEnv.parse("TEST=  something  ")["TEST"] == "something"
    @test ConfigEnv.parse("TEST=    ")["TEST"] == ""
end

end # module


text:
@testset "Identity" begin

    @testset "Float64" begin

        @test_throws ErrorException percentidentity(res"AH", res"AGH")

        @test percentidentity(res"AH", res"AH") == 100.
        @test percentidentity(res"AH", res"AG") == 50.
        @test percentidentity(res"AH", res"RG") == 0.
        @test percentidentity(res"AH-", res"AG-") == 50.
        @test percentidentity(res"A--", res"AG-") == 50.
        # Columns with XAA aren't counted
        @test percentidentity(res"AXA-", res"AG--") == 100 .* (1+0+0+0) / (1+0+1+0)
        @test percentidentity(res"AAX-", res"AG--") == 100 .* (1+0+0+0) / (1+1+0+0)
        @test percentidentity(res"AH-", res"AX-") == 100.
        @test percentidentity(res"AH-", res"XG-") == 0.
        @test percentidentity(res"AGG", res"AHX") == 50.

        @test isnan(percentidentity(res"---", res"---"))
        @test isnan(percentidentity(res"XX-", res"-XX"))
        @test isnan(percentidentity(res"XXX", res"XXX"))
    end

    @testset "MSA 2x3" begin
        αβ = (1,20,21,22)
        for i₁ in αβ, j₁ in αβ, k₁ in αβ
            seq₁ = Residue[i₁,j₁,k₁]
            for i₂ in αβ, j₂ in αβ, k₂ in αβ
                seq₂ = Residue[i₂,j₂,k₂]
                val = percentidentity(seq₁, seq₂)
                if  sum((seq₁ .== GAP) .& (seq₂ .== GAP)) +
                    sum((seq₁ .== XAA) .| (seq₂ .== XAA)) == 3

                    @test isnan(val)
                elseif seq₁ == seq₂
                    @test val == 100.0
                else
                    @test 100.0 >= val >= 0.0
                end
            end
        end
    end

    @testset "Bool" begin

        @test percentidentity(res"A--", res"AG-", 40.)
        @test !percentidentity(res"A--", res"AG-",60)
        # Columns with XAA aren't counted
        @test percentidentity(res"AXA-", res"AG--", 40.)
        @test !percentidentity(res"AAX-", res"AG--", 60.)
        @test percentidentity(res"AH-", res"AX-", 100.)
        @test percentidentity(res"AH-", res"XG-", 0.)
        @test !percentidentity(res"AGG", res"AHX", 62.) # 50% < 62%
    end

    @testset "MSA" begin

        fasta = read(joinpath(pwd(), "data", "Gaoetal2011.fasta"), FASTA)
        id = percentidentity(fasta)

        @test id[1,1] == 100.0

        @test isapprox(id[1,2], 83.33, atol=0.01)
        @test isapprox(id[1,3], 83.33, atol=0.01)
        @test isapprox(id[3,1], 83.33, atol=0.01)
        @test isapprox(id[1,6], 33.33, atol=0.01)
        @test isapprox(id[4,5], 83.33, atol=0.01)

        @test id[5,6] == 100.0

        @test maximum(id) == 100.0
        @test isapprox(minimum(id), 33.33, atol=0.01)

        @testset "SequenceIdentityMatrix" begin

            @test isapprox(sum(id[:,3]), 100.0 + 50.0 + 2*(200/6) + 2*(500/6))
            id[4,3] = 80
            @test isapprox(sum(id[:,3]), 100.0 + 80.0 + 2*(200/6) + 2*(500/6))
        end

        @testset "Gaps" begin

            aln = read(joinpath(pwd(), "data", "gaps.txt"), Raw)
            id = percentidentity(aln)

            @test id[1,1] == 100.0
            @test id[1,2] == 90.0
            @test id[1,3] == 80.0
            @test id[2,3] ≈ 800/9
        end

        @testset "Mean percent identity" begin

            aln = Residue[ '-' '-' 'G' 'G' 'G' '-'
                           '-' '-' '-' 'G' 'G' 'G' ]
            # identities    0   0   0   1   1   0 sum 2
            # aligned res   0   0   1   1   1   1 sum 4
            @test percentidentity(aln)[1, 2] == 50.0 # 2 / 4
            @test meanpercentidentity(aln)   == 50.0

            msa    = rand(Residue, 400, 2)
            msa300 = msa[1:300, :]

            @test mean(getlist(percentidentity(msa300))) == meanpercentidentity(msa300)
            @test isapprox( mean(getlist(percentidentity(msa))),
                            meanpercentidentity(msa), atol=0.5 )
            @test mean(getlist(percentidentity(msa))) == meanpercentidentity(msa,exact=true)
        end
    end

    @testset "Percent Similarity" begin

        fasta = read(joinpath(pwd(), "data", "Gaoetal2011.fasta"), FASTA)

        @testset "Gaps" begin

            @test percentsimilarity(res"AH", res"IM") == 50.0

            @test percentsimilarity(res"-AH", res"--H") == 50.0
            @test percentsimilarity(res"-AH", res"-AH") == 100.0
            @test percentsimilarity(res"-AH",res"-AH") ==
                percentsimilarity(res"AH", res"AH")
            @test percentsimilarity(res"-AH",res"-AH") ==
                percentsimilarity(res"--AH", res"--AH")
            # Residues outside the alphabet aren't counted (i.e.: XAA)
            @test percentsimilarity(res"-AH",res"-AH") ==
                percentsimilarity(res"XXAH", res"AAAH")
                @test percentsimilarity(res"AH",res"AH") ==
                    percentsimilarity(res"GGAH", res"AAAH",
                        ReducedAlphabet("AH"))
        end

        @testset "Using SMS's Ident and Sim residue groups" begin

            sim = percentsimilarity(fasta,
                  ReducedAlphabet("(GAVLI)(FYW)(ST)(KRH)(DENQ)P(CM)"))

            @test eltype(sim) == Float64
            @test sim[1,1] == 100.0

            @test isapprox(sim[1,2],  83.33, atol=0.01)
            @test isapprox(sim[1,3], 100.00, atol=0.01)
            @test isapprox(sim[1,4],  66.67, atol=0.01)
            @test isapprox(sim[1,5],  50.00, atol=0.01)
            @test isapprox(sim[1,6],  50.00, atol=0.01)
            @test isapprox(sim[2,3],  83.33, atol=0.01)
            @test isapprox(sim[2,4],  50.00, atol=0.01)
            @test isapprox(sim[2,5],  50.00, atol=0.01)
            @test isapprox(sim[2,6],  50.00, atol=0.01)
            @test isapprox(sim[3,4],  66.67, atol=0.01)
            @test isapprox(sim[3,5],  50.00, atol=0.01)
            @test isapprox(sim[3,6],  50.00, atol=0.01)
            @test isapprox(sim[4,5],  83.33, atol=0.01)
            @test isapprox(sim[4,6],  83.33, atol=0.01)
            @test isapprox(sim[5,6], 100.00, atol=0.01)
        end

        @testset "Using Bio3D's (2.2) seqidentity residue groups" begin

            sim = percentsimilarity(fasta,
                  ReducedAlphabet("(GA)(MVLI)(FYW)(ST)(KRH)(DE)(NQ)PC"), out=Float16)
            bio3d = [   1.000 0.833 1.000 0.667 0.500 0.500
                        0.833 1.000 0.833 0.500 0.500 0.500
                        1.000 0.833 1.000 0.667 0.500 0.500
                        0.667 0.500 0.667 1.000 0.833 0.833
                        0.500 0.500 0.500 0.833 1.000 1.000
                        0.500 0.500 0.500 0.833 1.000 1.000 ] .* 100.00

            @test eltype(sim) == Float16

            for i in 1:6
                for j in 1:6
                    @test isapprox(sim[i,j], bio3d[i,j], atol=0.1)
                end
            end
        end

        @testset "Percent identity" begin

            for αβ in (UngappedAlphabet(), GappedAlphabet())

                @test_throws ErrorException percentsimilarity(res"AH", res"AGH", αβ)

                @test percentsimilarity(res"AH", res"AH", αβ) == 100.
                @test percentsimilarity(res"AH", res"AG", αβ) == 50.
                @test percentsimilarity(res"AH", res"RG", αβ) == 0.
                @test percentsimilarity(res"AH-",res"AG-", αβ) == 50.
                @test percentsimilarity(res"A--",res"AG-", αβ) == 50.
                # Columns with XAA aren't counted
                @test percentsimilarity(res"AXA-",res"AG--", αβ) == 100 .* (1+0+0+0) / (1+0+1+0)
                @test percentsimilarity(res"AAX-",res"AG--", αβ) == 100 .* (1+0+0+0) / (1+1+0+0)
                @test percentsimilarity(res"AH-", res"AX-", αβ) == 100.
                @test percentsimilarity(res"AH-", res"XG-", αβ) == 0.
                @test percentsimilarity(res"AGG", res"AHX", αβ) == 50.

                @test isnan(percentsimilarity(res"---", res"---", αβ))
                @test isnan(percentsimilarity(res"XX-", res"-XX", αβ))
                @test isnan(percentsimilarity(res"XXX", res"XXX", αβ))
            end
        end
    end
end


text:
# EXCLUDE FROM TESTING
if Base.find_package("CUDA") !== nothing
    using CUDA
    using CUDAKernels
    const has_cuda = true
else
    const has_cuda = false
end
if Base.find_package("AMDGPU") !== nothing
    using AMDGPU
    using ROCKernels
    const has_rocm = true
else
    const has_rocm = false
end


text:
function subbac!(u::Matrix{Float64}, b::Vector{Float64})
# Back-substitution on an Upper Triangle
  n = size(u, 1)
  for i in n:-1:1
    total = b[i]
    if i < n
      for j in i+1:n
        total -= u[i, j] * b[j]
      end
    end
    b[i] = total / u[i, i]
  end
end

function subbac(u::Matrix{Float64}, b::Vector{Float64})
# Back-substitution on an Upper Triangle
  bt = deepcopy(b)
  n = size(u, 1)
  for i in n:-1:1
    total = bt[i]
    if i < n
      for j in i+1:n
        total -= u[i, j] * bt[j]
      end
    end
    bt[i] = total / u[i, i]
  end
  bt
end

export
  subbac!,
  subbac


text:
### Global surface temperature map

for i = 1:length(casenames)

    fig, ax = plt[:subplots](1, 1, figsize=(12,8))

    casename = casenames[i]

    if match(r"STD", casename) != nothing
        continue
    end

    cmap = plt[:get_cmap]("coolwarm")
    cmap[:set_over]("r")
    cmap[:set_under]("b")
    cntrs = range(-2, step=2, stop=30)

    Dataset("$nc_file_dir/$casename.ocn.h.ma.fv4x5.0001-0020.nc", "r") do ds

        println(ds.dim["time"])

        if match(r"SSM_NK", casename) != nothing
            T_ML = ds["T"][:, :, 1, t_rng]
        elseif match(r"SSM_SOM", casename) != nothing
            T_ML = ds["T_ML"][:, :, t_rng]
        end

        T_ML = mreplace(T_ML)
        T_ML = mean(T_ML, dims=(3,))[:, :, 1] .- 273.15

        cmapping = plt[:contourf](lon, lat, T_ML'[:,:], cntrs, cmap=cmap, antialised=true, extend="both")

        #ax[:contour](lon, lat, mask'[:,:], [.5], colors="#000000")

        cba = plt[:colorbar](cmapping, ax=ax)

        ax[:set_title](format("Mean SST ({:d} years): $casename", years))
        ax[:set_xlabel]("Lat [deg]")
        ax[:set_ylabel]("Lon [deg]")
        cba[:set_label]("Temperature [degC]")
    end

    fig[:savefig]("$(casename)_map_SST.png", dpi=200)

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
using Pkg
Pkg.test("DynamicBoundspODEsDiscrete"; coverage=true)


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# Read in / generate multiple deployments, and simulate each one.

using JEMSS
using Random

# parameters:
configFilename = "sim_config.xml"
outputFolder = "output"
conf = 0.95 # statistical confidence level, note that sim output files may have different conf value
doPrint = true
# parallel (multi-threading):
parallel = false
numThreads = 1
# additional parameters to use if config does not include deployments:
deploymentsOutputFilename = joinpath(outputFolder, "deployments.csv")
numDeployments = 2 # number of deployments to generate and test
deploymentRng = Random.MersenneTwister(0) # for reproducible results, if generating random deployments

# parameter checks
@assert(isfile(configFilename))
@assert(0 <= conf < 1) # should be between 0 and 1, but not equal to 1 otherwise confidence interval is infinite
@assert(isa(parallel, Bool))
@assert(numThreads >= 1)

isdir(outputFolder) || mkdir(outputFolder)

# load sim
doPrint && println("Loading sim.")
sim = initSim(configFilename, createBackup = true)
(numAmbs, numStations) = (sim.numAmbs, sim.numStations)
@assert(numAmbs > 0 && numStations > 1) # otherwise only one deployment is possible
@assert(sim.numReps >= 2, "Need at least two simulation replications (sim.reps) for statistics.")

# read deployments file
if haskey(sim.inputFiles, "deployments")
	(deployments, numStations2) = readDeploymentsFile(sim.inputFiles["deployments"])
	@assert(numStations == numStations2)
elseif haskey(sim.inputFiles, "deploymentPolicies") # compat
	(deployments, numStations2) = readDeploymentPoliciesFile(sim.inputFiles["deploymentPolicies"])
	@assert(numStations == numStations2)
else
	# generate unique random deployments
	deployments = makeRandDeployments(numAmbs, numStations, numDeployments; rng = deploymentRng)
	doPrint && println("Deployments input file not found in config, generated $numDeployments deployments.")
end
numDeployments = length(deployments)

writeDeploymentsFile(deploymentsOutputFilename, deployments, numStations)

# convert deployments, for ease of use
stationsNumAmbsList = [deploymentToStationsNumAmbs(deployment, numStations) for deployment in deployments]

# check which output files sim will write
simOutputFiles = filter(name -> haskey(sim.outputFiles, name), ["ambulancesStats", "callsStats", "hospitalsStats", "stationsStats", "statsDict"])
simOutputFilepaths = map(name -> sim.outputFiles[name].path, simOutputFiles)

# run simulation for each deployment, save stats
statsPeriodsLists = [] # statsPeriodsLists[i][j] is for deployment i, replication j
statsDicts = []
for (i, stationsNumAmbs) in enumerate(stationsNumAmbsList)
	doPrint && print("\rSimulating deployment $i (of $numDeployments).")
	function runRep!(rep::Simulation)
		resetRep!(rep)
		applyStationsNumAmbs!(rep, stationsNumAmbs)
		simulateRep!(rep)
		@assert(getStationsNumAmbs(rep) == stationsNumAmbs) # check that rep had the current deployment applied
	end
	if parallel == true
		runParallel!(runRep!, sim.reps...; numThreads = numThreads)
	else
		for rep in sim.reps runRep!(rep) end
	end

	writeStatsFiles(sim, sim.reps)

	# move output files into subfolder
	folder = joinpath(outputFolder, string(i))
	isdir(folder) || mkdir(folder)
	for filepath in simOutputFilepaths
		if isfile(filepath) mv(filepath, joinpath(folder, basename(filepath)); force = true) end
	end

	# get periods stats
	periods = getRepsPeriodStatsList(sim.reps)
	statsDict = statsDictFromPeriodStatsList(periods, conf = conf)
	push!(statsPeriodsLists, periods)
	push!(statsDicts, statsDict)

	# write deployment to file
	deploymentFilename = "deployment.csv"
	deployment = stationsNumAmbsToDeployment(stationsNumAmbs)
	writeDeploymentsFile(joinpath(folder, deploymentFilename), [deployment], numStations)
end
doPrint && println()
doPrint && println("Done.")

# write objective values (fraction of calls reached in time) to file
objVals = [statsDict["calls"]["all"]["fracResponsesInTime"] for statsDict in statsDicts]
objValsTable = Table("objVals", ["deploymentIndex", "mean", "halfWidth"];
	rows = [[string(i), objVals[i].mean, objVals[i].halfWidth] for i = 1:numDeployments])
writeTablesToFile(joinpath(outputFolder, "obj_vals.csv"), [objValsTable])


text:
"""
functions that enable modyfing of reactions to mouse using keyboard
for example by pressing f and s we can controll wheather we have fast or slow scroll
"""
module KeyboardMouseHelper
using  ..StructsManag, Logging, Setfield, ..PrepareWindow,   ..DataStructs , Rocket, GLFW, ..ForDisplayStructs, ..TextureManag,  ..OpenGLDisplayUtils,  ..Uniforms, Match, Parameters,DataTypesBasic

"""
For controlling the  fast scroll  on pressing normal f key and stopping it on n key
"""
function processKeysInfo(isTobeFast::Identity{Tuple{Bool,Bool}}
                        ,actor::SyncActor{Any, ActorWithOpenGlObjects}
                        ,keyInfo::KeyboardStruct
                        ,toBeSavedForBack::Bool = true) where T

                        isTobeFastVal = isTobeFast.value[1]
    #passing information to actor that we should do now fast scrolling
    # actor.actor.mainForDisplayObjects= setproperties .isFastScroll = isTobeFastVal
       actor.actor.mainForDisplayObjects=setproperties(actor.actor.mainForDisplayObjects, (isFastScroll=isTobeFastVal))
# for undoing action
if(toBeSavedForBack)
    addToforUndoVector(actor, ()-> processKeysInfo( Option((!isTobeFastVal,false)),actor, keyInfo,false ))
end

end#processKeysInfo

end#KeyboardMouseHelper


text:
using Plots
using DelimitedFiles

include("options.jl")

println()
print("Plotting results for $blas_num_threads BLAS thread")
blas_num_threads > 1 ? println("s") : println()

println()
println("Plotting results for $benchmarks")
println()

println("Plotting results for dimensions set to:")
display(maxdims)
println()

seperator = "#"^70

outputlevel = 1

for benchmark in benchmarks
  println("Plotting benchmark $benchmark for bond dimensions $(maxdims[benchmark])")
  println()

  #N = length(maxdims[benchmark])
  times_julia = Float64[]
  times_cpp = Float64[]

  for maxdim in maxdims[benchmark]
    if isnothing(which_version) || which_version == "julia"
      julia_dir = joinpath(@__DIR__, "bench_$benchmark", "julia")
      filename = "data_blas_num_threads_$(blas_num_threads)"
      filename *= "_maxdim_$(maxdim).txt"
      filepath = joinpath(julia_dir, "data", filename)
      time = readdlm(filepath)[]
      push!(times_julia, time)
    end
    if isnothing(which_version) || which_version == "c++"
      cpp_dir = joinpath(@__DIR__, "bench_$benchmark", "c++")
      filename = "data_blas_num_threads_$(blas_num_threads)"
      filename *= "_maxdim_$(maxdim).txt"
      filepath = joinpath(cpp_dir, "data", filename)
      time = readdlm(filepath)[]
      push!(times_cpp, time)
    end
  end
  p = plot(maxdims[benchmark], times_julia;
           title = descriptions[benchmark],
           legendfontsize = 14,
           xlabel = "Maximum bond dimension",
           ylabel = "Computation time (seconds)",
           xtickfontsize = 14,
           ytickfontsize = 14,
           xguidefontsize = 14,
           yguidefontsize = 14,
           line = (:solid, 4),
           marker = (:dot, 8),
           legend = :topleft,
           label = "Julia")
  plot!(p, maxdims[benchmark], times_cpp;
        line = (:solid, 4),
        marker = (:dot, 8),
        label = "C++")
  filename = "plot_blas_num_threads_$(blas_num_threads)"
  filename *= "_$benchmark.png"
  savefig(p, joinpath("plots", filename))
end


text:
# Plain arrays can be treated as images. Other types will have
# metadata associated, make yours a child of one of the following:
abstract AbstractImage{T,N} <: AbstractArray{T,N}         # image with metadata
abstract AbstractImageDirect{T,N} <: AbstractImage{T,N}   # each pixel has own value/color
abstract AbstractImageIndexed{T,N} <: AbstractImage{T,N}  # indexed images (i.e., lookup table)

# Direct image (e.g., grayscale, RGB)
type Image{T,N,A<:AbstractArray} <: AbstractImageDirect{T,N}
    data::A
    properties::Dict{ASCIIString,Any}
end
# Indexed image (colormap)
type ImageCmap{T<:Color,N,A<:AbstractArray} <: AbstractImageIndexed{T,N}
    data::A
    cmap::Vector{T}
    properties::Dict{ASCIIString,Any}
end
ImageCmap(data::AbstractArray, cmap::AbstractVector, props::Dict) = ImageCmap{eltype(cmap),ndims(data),typeof(data)}(data, cmap, props)
ImageCmap(data::AbstractArray, cmap::AbstractVector; kwargs...) = ImageCmap(data, cmap, kwargs2dict(kwargs))

Image(data::AbstractArray, props::Dict) = Image{eltype(data),ndims(data),typeof(data)}(data,props)
Image(data::AbstractArray; kwargs...) = Image(data, kwargs2dict(kwargs))

Base.eltype{T}(img::AbstractImage{T}) = T

Base.size(img::AbstractImage) = size(img.data)
Base.size(img::AbstractImage, i::Integer) = size(img.data, i)
Base.size(img::AbstractImage, dimname::String) = size(img.data, dimindex(img, dimname))

Base.ndims(img::AbstractImage) = ndims(img.data)

Base.strides(img::AbstractImage) = strides(img.data)

Base.copy(img::Image) = Image(copy(img.data), deepcopy(img.properties))
Base.copy(img::ImageCmap) = ImageCmap(copy(img.data), copy(img.cmap), deepcopy(img.properties))

const emptyset = Set()
function showim(io::IO, img::AbstractImageDirect)
    IT = typeof(img)
    print(io, colorspace(img), " ", IT.name, " with:\n  data: ", summary(img.data), "\n  properties:")
    showdictlines(io, img.properties, get(img, "suppress", emptyset))
end
function showim(io::IO, img::AbstractImageIndexed)
    IT = typeof(img)
    print(io, colorspace(img), " ", IT.name, " with:\n  data: ", summary(img.data), "\n  cmap: ", summary(img.cmap), "\n  properties:")
    showdictlines(io, img.properties, get(img, "suppress", emptyset))
end
Base.show(io::IO, img::AbstractImageDirect) = showim(io, img)
Base.writemime(io::IO, ::MIME"text/plain", img::AbstractImageDirect) = showim(io, img)
Base.show(io::IO, img::AbstractImageIndexed) = showim(io, img)
Base.writemime(io::IO, ::MIME"text/plain", img::AbstractImageIndexed) = showim(io, img)
function showdictlines(io::IO, dict::Dict, suppress::Set)
    for (k, v) in dict
        if k == "suppress"
            continue
        end
        if !in(k, suppress)
            print(io, "\n    ", k, ": ")
            printdictval(io, v)
        else
            print(io, "\n    ", k, ": <suppressed>")
        end
    end
end
macro get(img, k, default)
    quote
        img, k = $(esc(img)), $(esc(k))
        local val
        if !isa(img, AbstractImage)
            val = $(esc(default))
        else
            index = Base.ht_keyindex(img.properties, k)
            val = (index > 0) ? img.properties.vals[index] : $(esc(default))
        end
        val
    end
end

colorspace{C<:Color}(img::AbstractMatrix{C}) = string(C.name.name)
colorspace{C<:Color}(img::AbstractArray{C,3}) = string(C.name.name)
colorspace{C<:Color}(img::AbstractImage{C}) = string(C.name.name)
colorspace{C<:Color,T}(img::AbstractArray{AbstractAlphaColor{C,T},2}) = (S = string(C.name.name); S == "Gray" ? "GrayAlpha" : string(S, "A"))
colorspace{C<:Color,T}(img::AbstractImage{AbstractAlphaColor{C,T}}) = (S = string(C.name.name); S == "Gray" ? "GrayAlpha" : string(S, "A"))
colorspace(img::AbstractMatrix{Bool}) = "Binary"
colorspace(img::AbstractArray{Bool}) = "Binary"
colorspace(img::AbstractArray{Bool,3}) = "Binary"
colorspace(img::AbstractMatrix{Uint32}) = "RGB24"
colorspace(img::AbstractMatrix) = "Gray"
colorspace{T}(img::AbstractArray{T,3}) = (size(img, defaultarraycolordim) == 3) ? "RGB" : error("Cannot infer colorspace of Array, use an AbstractImage type")
colorspace(img::AbstractImage{Bool}) = "Binary"
colorspace{T,N,A<:AbstractArray}(img::ImageCmap{T,N,A}) = string(T.name.name)
colorspace(img::AbstractImageIndexed) = @get img "colorspace" csinfer(eltype(img.cmap))
colorspace{T}(img::AbstractImageIndexed{T,2}) = @get img "colorspace" csinfer(eltype(img.cmap))
csinfer{C<:Color}(::Type{C}) = string(C)
csinfer(C) = "Unknown"
colorspace(img::AbstractImage) = get(img.properties, "colorspace", "Unknown")

properties(A::AbstractArray) = @compat Dict(
    "colorspace" => colorspace(A),
    "colordim" => colordim(A),
    "timedim" => timedim(A),
    "pixelspacing" => pixelspacing(A),
    "spatialorder" => spatialorder(A))
properties{C<:Color}(A::AbstractArray{C}) = @compat Dict(
    "timedim" => timedim(A),
    "pixelspacing" => pixelspacing(A),
    "spatialorder" => spatialorder(A))
properties(img::AbstractImage) = img.properties

Base.haskey(a::AbstractArray, k::String) = false
Base.haskey(img::AbstractImage, k::String) = haskey(img.properties, k)

Base.get(img::AbstractArray, k::String, default) = default
Base.get(img::AbstractImage, k::String, default) = get(img.properties, k, default)

printdictval(io::IO, v) = print(io, v)
function printdictval(io::IO, v::Vector)
    for i = 1:length(v)
        print(io, " ", v[i])
    end
end


text:
using DataAssimilation
using Documenter

# NOTE: geostats.jl has good example documentation and structure in their make.jl

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

makedocs(;
    modules=[DataAssimilation],
    authors="John Waczak <john.louis.waczak@gmail.com>",
    repo="https://github.com/john-waczak/DataAssimilation.jl/blob/{commit}{path}#{line}",
    sitename="DataAssimilation.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://john-waczak.github.io/DataAssimilation.jl",
        assets=String[],
        mathengine = KaTeX(Dict(
            :macros => Dict(
                "\\cov" => "\\text{cov}",
                "\\R"   => "\\mathbb{R}",
                "\\E"   => "\\mathbb{E}",
            )
        ))
    ),
    pages=[
        "Home" => "index.md",
        "Data Assimilation" => [
            "Data Assimilation Oveview" => "dataAssimilation/dataassim.md",
            "Kalman Filter" => "dataAssimilation/kalman.md",
            "Extended Kalman Filter" => "dataAssimilation/ekf.md",
            "3d var" => "dataAssimilation/3dvar.md",
            "4d var" => "dataAssimilation/4dvar.md",
            "Ensemble Kalman Filter" => "dataAssimilation/enkf.md",
            "Sensativity Analysis" => "dataAssimilation/sensativity.md",
        ],
        "Function docs" => "fdocs.md",
    ],
)

deploydocs(;
    repo="github.com/john-waczak/DataAssimilation.jl",
    devbranch="main",
)


text:
using BlockDiagonalMatrices
using SparseArrays
using LinearAlgebra
using TimerOutputs

struct CollectiveDephasingModel <: Model
    params::ModelParameters
    Jx::SparseMatrixCSC
    Jy::SparseMatrixCSC
    Jz::SparseMatrixCSC
    Jx2::SparseMatrixCSC
    Jy2::SparseMatrixCSC
    Jz2::SparseMatrixCSC
    second_terms::Array{SparseMatrixCSC}
    M0::SparseMatrixCSC
    dM::SparseMatrixCSC
    measurement::Eigen

    function CollectiveDephasingModel(modelparams::ModelParameters, liouvillianfile::Union{String, Nothing}=nothing)
        Nj = modelparams.Nj
        dt = modelparams.dt
        kcoll = modelparams.kcoll
        Gamma = modelparams.Gamma
        ω = modelparams.omega
        η = modelparams.eta

        firstblock = block_sizes(Nj)[1]

        # Spin operators (only the first block)
        (Jx, Jy, Jz) = map(x -> sparse(x[1:firstblock, 1:firstblock]), jspin(Nj))

        # TODO: Find better name
        second_terms = [sqrt((1 - η) * dt * Gamma) * Jy, sqrt(dt * kcoll) * Jz]

        Jx2 = Jx^2
        Jy2 = Jy^2
        Jz2 = Jz^2

        H = ω * Jz
        dH = Jz

        # Kraus-like operator, trajectory-independent part
        M0 = I - 1im * H * dt - (Gamma/2) * Jy2 * dt - (kcoll/2) * Jz2 * dt

        # Derivative of the Kraus-like operator wrt to ω
        dM = -1im * dH * dt

        measurement = eigen(Matrix(Jy))

        new(modelparams, Jx, Jy, Jz, Jx2, Jy2, Jz2, second_terms, M0, dM, measurement)
    end
end

function measure_current(state::FixedjState, model::CollectiveDephasingModel)
    @inline dW() = sqrt(model.params.dt) * randn() # Define the Wiener increment
    # Homodyne current (Eq. 35)
    # dy = 2 sqrt(Gamma * eta) * tr(ρ * Jy) * dt + dW
    mul!(state._tmp1, model.Jy, state.ρ)
    return 2 * sqrt(model.params.Gamma * model.params.eta) * real(tr(state._tmp1)) * model.params.dt + dW()
end

function get_kraus(model::CollectiveDephasingModel, dy::Real)
    p = model.params
    kraus = similar(model.M0)

    # Kraus operator Eq. (36)
    kraus = (model.M0 + sqrt(p.eta * p.Gamma) * model.Jy * dy +
            p.eta * (p.Gamma / 2) * model.Jy2 * (dy^2 - p.dt))

    return kraus
end

function updatestate!(state::FixedjState, model::CollectiveDephasingModel, dy::Real)
    # Get Kraus operator for measured current dy
    @timeit_debug "get_kraus" M = get_kraus(model, dy)

    # Non-allocating code for
    # new_ρ = Mpre * Mpost * ρ + second_term * ρ

    @timeit_debug "rho" begin
        mul!(state._tmp1, state.ρ, M')
        mul!(state._new_ρ, M, state._tmp1)

        # Apply the second term of Eq. (34)
        for second_term in model.second_terms
            mul!(state._tmp1, state.ρ, second_term')
            mul!(state._new_ρ, second_term, state._tmp1, 1., 1.)
        end

        zchop!(state._new_ρ) # Round off elements smaller than 1e-14
        tr_ρ = tr(state._new_ρ)
        # Evolve the unnormalized derivative wrt ω
    end

    # Non-allocating code for
    # τ = (M * (τ * M' + ρ * dM') + dM * ρ * M' +
    #     ∑ second_term[i] * τ * second_term[i]') / tr_ρ;
    @timeit_debug "tau" begin

        fill!(state._tmp2, 0.0)
        # ∑ second_term[i] * τ * second_term[i]')
        for second_term in model.second_terms
            mul!(state._tmp1, second_term, state.τ)
            mul!(state._tmp2, state._tmp1, second_term', 1., 1.)
        end

        # (τ * M' + ρ * dM')
        mul!(state._tmp1, state.ρ, model.dM')
        mul!(state._tmp1, state.τ, M', 1., 1.)

        copy!(state.τ, state._tmp2)

        #  (M * (τ * M' + ρ * dM')
        mul!(state.τ, M, state._tmp1, 1., 1.)

        #  ... + dM * ρ * M'
        mul!(state._tmp1, state.ρ, M')
        mul!(state.τ, model.dM, state._tmp1, 1., 1.)

        state.τ ./=  tr_ρ
    end

    zchop!(state.τ) # Round off elements smaller than 1e-14

    tr_τ = tr(state.τ)
    # Now we can renormalize ρ and its derivative wrt ω
    # TODO: Use broadcasting when it is implemented
    @timeit_debug "norms" begin
        state.ρ  .= state._new_ρ ./ tr_ρ
        state.dρ .= state.τ .- tr_τ .* state.ρ
    end
    return tr_ρ, tr_τ
end

coherentspinstate(model::CollectiveDephasingModel) = fixedj_css(model.params.Nj)


text:
"""Returns the coriolis parameter on the q-grid for beta plane approximation."""
function beta_plane()

    if bc_x == "periodic"
        # points on the right edge needed too
        xx_q,yy_q = meshgrid(x_q_halo[3:end-1],y_q_halo[3:end-2])
    else
        xx_q,yy_q = meshgrid(x_q,y_q)
    end

    f₀ = 2*ω*sind(ϕ)      # Coriolis parameter [s^-1] at central latitude
    β = 2*ω/R*cosd(ϕ)     # Coriolis derivative wrt latitude [(ms)^-1]

    # for non-dimensional gradient operators f contains the grid spacing Δ
    f_q = Numtype.(Δ*(f₀ .+ β*yy_q))

    return f_q
end

max_stars_repo_name: eschnett/juls

max_issues_repo_name: eschnett/juls

max_forks_repo_name: eschnett/juls


text:
struct TemperatureProfile{AMP, FWHM, U, V}
    # amplitude and full-width half-max functions for generalized Gaussian temperature profile
    amplitude::AMP
    width::FWHM
    # parameters governing uncertainty calculation
    dT::U
    dx::V
end

# Constructors:
# defaults to parameters of Science Advances submission
function TemperatureProfile(dT::Real = 20, dx::Real = 10/1000)
    get_temperature_profile_Bi2O3_2021(dT, dx)
end

# template for t-profile constructor
function get_temperature_profile_template(dT::Real = 20, dx::Real = 10/1000)
    println("temperature_profile_template")
    amplitude(T_max, log10_τ) = 2
    function width(T, log10_τ)
        throw("not implemented - should be a functional expression from temperature and dwell time to the fwhm of the profile")
    end
    TemperatureProfile(amplitude, width, dT, dx)
end

# temperature profile for the Bi2O3 data in the Science Advances submission
function get_temperature_profile_Bi2O3_2021(dT::Real = 20, dx::Real = 10/1000)
    amplitude(T_max, log10_τ) = 2
    function width(T, log10_τ)
        fwhm_0 = 1.60742002e+03
        a = (-1.46900821e+00, -2.06060111e+02, 3.05703021e-01, 3.29193895e+02, -1.23361265e+13, 1.09165256e+01, -1.69547934e+04, -1.42207444e-04)
        fwhm_0 + a[1]*T + a[2]*log10_τ + a[3]*log10_τ*T + a[4]*log10_τ^4/T +
            a[5]/log10_τ^8/T^3 + a[6]*T/log10_τ^3 + a[7]/log10_τ^4 + a[8]*T^3/log10_τ^8
    end
    TemperatureProfile(amplitude, width, dT, dx)
end

# OLD temperature profile for the Bi2O3 data
function get_temperature_profile_Bi2O3_2020(dT::Real = 20, dx::Real = 10/1000)
    amplitude(T_max, log10_τ) = 2
    function width(T, log10_τ)
        fwhm_0 = 6.21599199e+01
        a = (-3.56990024e-01, 2.84613876e+02, -8.90507983e-03, -5.76367089e+02, -1.95522786e+13, 5.51350589e+00, -3.34684841e+03, -1.33972469e-04)
        fwhm_0 + a[1]*T + a[2]*log10_τ + a[3]*log10_τ*T + a[4]*log10_τ^4/T +
            a[5]/log10_τ^8/T^3 + a[6]*T/log10_τ^3 + a[7]/log10_τ^4 + a[8]*T^3/log10_τ^8
    end
    TemperatureProfile(amplitude, width, dT, dx)
end

# temperature profile for CHESS 2021
function get_temperature_profile_CHESS_2021(dT::Real = 20, dx::Real = 10/1000)
    amplitude(T_max, log10_τ) = 2
    # power = current in A
    # tau = dwell in us
    function _std(power, log10_τ)
        mm_per_pixel = 0.00153846153 # millimeters per pixel
        μm_per_pixel = mm_per_pixel * 1e3 # micrometers per pixel
        μm_per_pixel * (-6.0585 * power + 5.29285 * log10_τ + 1.432e-2 * power^2 - 1.59943e1 *log10_τ^2 + 7.6715e-1 * power*log10_τ + 726.0242)
    end
    _power(T_max::Int, log10_τ) = _power(float(T_max), float(log10_τ))
    _power(T_max, log10_τ) = (5000*sqrt(5)*sqrt(614000*T_max+1327591125*log10_τ^2 -8080505000*log10_τ+12142813538)-341554591*log10_τ+1065621710)/3070000
    _fwhm(power, log10_τ) = _std(power, log10_τ) / 2 # by the definitions of variables in the Gaussian-like function Max used
    width(T_peak, log10_τ) = _fwhm(_power(T_peak, log10_τ), log10_τ) # getting the power from the power profile
    TemperatureProfile(amplitude, width, dT, dx)
end

# Functor definitions:
# WARNING: temperature profile is a function of position in mm
function (P::TemperatureProfile)(T_max::Real, log10_τ::Real)
    a = P.amplitude(T_max, log10_τ)
    f = P.width(T_max, log10_τ)
    x->_T_helper(T_max, f, a, 1e3*x) # multiplying position by 1e3 to convert from mm to μm
end
function (P::TemperatureProfile)(T_max::Real, log10_τ::Real, x::Real)
    P(T_max, log10_τ)(x)
end
# WARNING: temperature profile is a function of position in μm
function _T_helper(T_max, fwhm, ampl, x::Real)
    return @. T_max * exp(-abs(2x / fwhm)^ampl)
end

profile(P::TemperatureProfile, T_max, log10_τ) = P(T_max, log10_τ)
# returns the positive side of the inverse mapping T to position in mm
function inverse_profile(P::TemperatureProfile, T_max::Real, log10_τ::Real)
    a = P.amplitude(T_max, log10_τ)
    f = P.width(T_max, log10_τ)
    function invprof(T)
        T > 0 || throw(DomainError("T is non-positive: $T"))
        T_max ≥ T || throw(DomainError("T exceeds T_max: T = $T > $T_max = T_max"))
        (T ≈ T_max) ? zero(eltype(T)) : log(T_max / T)^(1/a) * f/2 * (1e-3)
    end
end

# x is a position
# computes uncertainty of the temperature profile
# Two components:
# - multiplicative and
# - proportional to derivative
function temperature_uncertainty(P::TemperatureProfile, T_max, log10_τ, dT::Real = P.dT, dx::Real = P.dx)
    f = P(T_max, log10_τ)
    function temperature_variance(x::Real)
        T, dTdx = f(x), ForwardDiff.derivative(f, x)
        (dT * T / 1400)^2 + (dx * dTdx)^2 # error variance is additive
    end
end


text:
import Base:in, +, -, *, ==, convert, isempty, isequal, length, push!, sizeof, append!
import SeisIO: GeoLoc, PZResp
export CorrData

const corrfields = [:name, :loc, :comp, :rotated, :corr_type, :fs, :gain,
                    :freqmin, :freqmax, :cc_len, :cc_step, :whitened, :time_norm, :resp,
                    :maxlag, :dist, :azi, :baz]

"""
   CorrData

A structure for cross-correlations of ambient noise data.

## Fields: CorrData
| **Field** | **Description** |
|:------------|:-------------|
| :name       | Freeform channel name in NET1.STA1.LOC1.CHAN1.NET2.STA2.LOC2.CHAN2 form |
| :id         | Channel ids. use NET.STA.LOC.CHAN format when possible. |
| :loc        | Location (position) object for station 1. |
| :fs         | Sampling frequency in Hz. |
| :gain       | Scalar gain; divide data by the gain to convert to units  |
| :freqmin    | Minimum frequency for whitening.  |
| :freqmax    | Maximum frequency for whitening. |
| :cc_len     | Length of each correlation in seconds. |
| :cc_step    | Spacing between correlation windows in seconds. |
| :whitened   | Whitening applied.
| :time_norm  | Apply one-bit whitening with "one_bit". |
| :resp       | Instrument response object, format [zeros poles] |
| :misc       | Dictionary for non-critical information. |
| :notes      | Timestamped notes; includes automatically-logged acquisition and |
|             | processing information. |
| :dist       | Distance between station 1 and station 2 in Km. |
| :azi        | Azimuth from station 1 and station 2 in degrees. |
| :baz        | Back azimuth between station 1 and station 2 in degrees. |
| :maxlag     | Maximum lag time in seconds to keep from correlations. |
| :t          | Starttime of each correlation. |
| :corr       | Correlations stored in columns. |
"""
mutable struct CorrData
  name::String                                # name [Net1.Sta1.Loc1.Chan1.Net2.Sta2.Loc2.Chan2]
  id::String                                  # id [Y-mm-dd] this is date of corr
  loc::GeoLoc                                 # loc
  comp::String                                # 1st channel, 2nd channel [ZZ,RT,..]
  rotated::Bool                               # rotation applied?
  corr_type::String                           # "corr", "deconv", "coh"
  fs::Float64                                 # sampling rate [Hz]
  gain::Float64                               # gain
  freqmin::Float64                            # minumum frequency [Hz]
  freqmax::Float64                            # maximum frequency [Hz]
  cc_len::Int                                 # window_length [s]
  cc_step::Int                                # step between windows [s]
  whitened::Bool                              # whitening applied
  time_norm::Union{Bool,String}               # time normaliation
  resp::PZResp                                # response poles/zeros
  misc::Dict{String,Any}                      # misc
  notes::Array{String,1}                      # notes
  dist::Float64                               # distance between stations [km]
  azi::Float64                                # azimuth between stations [degrees]
  baz::Float64                                # back azimuth between stations [degrees]
  maxlag::Float64                             # maximum lag time [s]
  t::Array{Float64,1}                         # time
  corr::Array{<:Union{Float32,Float64},2}     # corr data

  function CorrData(
      name     ::String,
      id       ::String,
      loc      ::GeoLoc,
      comp     ::String,
      rotated  ::Bool,
      corr_type::String,
      fs       ::Float64,
      gain     ::Float64,
      freqmin  ::Float64,
      freqmax  ::Float64,
      cc_len   ::Int,
      cc_step  ::Int,
      whitened ::Bool,
      time_norm::Union{Bool,String},
      resp     ::PZResp,
      misc     ::Dict{String,Any},
      notes    ::Array{String,1},
      dist     ::Float64,
      azi      ::Float64,
      baz      ::Float64,
      maxlag   ::Float64,
      t        ::Array{Float64,1},
      corr     ::Array{<:Union{Float32,Float64},2}
      )

      return new(name, id, loc, comp, rotated, corr_type, fs, gain, freqmin,
                 freqmax, cc_len, cc_step, whitened, time_norm, resp, misc,
                 notes, dist, azi, baz, maxlag, t, corr)
    end
end

CorrData(;
          name     ::String                    = "",
          id       ::String                    = "",
          loc      ::GeoLoc                    = GeoLoc(),
          comp     ::String                    = "",
          rotated  ::Bool                      = false,
          corr_type::String                    = "",
          fs       ::Float64                   = zero(Float64),
          gain     ::Float64                   = one(Float64),
          freqmin  ::Float64                   = zero(Float64),
          freqmax  ::Float64                   = zero(Float64),
          cc_len   ::Int                       = zero(Int),
          cc_step  ::Int                       = zero(Int),
          whitened ::Bool                      = false,
          time_norm::Union{Bool,String}        = false,
          resp     ::PZResp                    = PZResp(),
          misc     ::Dict{String,Any}          = Dict{String,Any}(),
          notes    ::Array{String,1}           = Array{String,1}(undef, 0),
          dist     ::Float64                   = zero(Float64),
          azi      ::Float64                   = zero(Float64),
          baz      ::Float64                   = zero(Float64),
          maxlag   ::Float64                   = zero(Float64),
          t        ::Array{Float64,1}          = Array{Float64,1}(undef, 0),
          corr     ::Array{<:Union{Float32,Float64},2} = Array{Float32,2}(undef, 0, 2)
          ) = CorrData(name, id, loc, comp, rotated, corr_type, fs, gain,
                      freqmin, freqmax, cc_len, cc_step, whitened, time_norm,
                      resp, misc, notes, dist, azi, baz, maxlag, t, corr)

CorrData(F1::FFTData, F2::FFTData, comp::String, rotated::Bool, corr_type::String,
        maxlag::Float64, t::Array{Float64,1}, corr::Array{<:Union{Float32,Float64},2}
       ) = CorrData(F1.name * '.' * F2.name, F1.id, F1.loc, comp, rotated,
                  corr_type, F1.fs, F1.gain, F1.freqmin, F1.freqmax, F1.cc_len,
                  F1.cc_step, F1.whitened, F1.time_norm, F1.resp, F1.misc,
                  F1.notes, get_dist(F1,F2), get_azi(F1,F2), get_baz(F1,F2),
                  maxlag, t, corr)

in(s::String, C::CorrData) = C.id==s

isempty(C::CorrData) = min(C.name=="",C.id == "",isempty(C.loc), C.comp=="",
                          C.rotated==false,C.fs == zero(typeof(C.fs)),
                          C.corr_type=="",
                          C.gain == one(typeof(C.gain)),
                          C.freqmin == zero(typeof(C.freqmin)),
                          C.freqmax == zero(typeof(C.freqmax)),
                          C.cc_len == zero(typeof(C.cc_len)),
                          C.cc_step == zero(typeof(C.cc_step)),
                          C.whitened == false,
                          C.time_norm == false,
                          isempty(C.resp), isempty(C.misc),isempty(C.notes),
                          C.maxlag == zero(typeof(C.maxlag)), isempty(C.t),
                          isempty(C.corr))

isequal(C::CorrData, U::CorrData) = minimum([hash(getfield(C,i))==hash(getfield(U,i)) for i in corrfields]::Array{Bool,1})
==(C::CorrData, U::CorrData) = isequal(C,U)::Bool

append!(C::CorrData, U::CorrData)  = (if C == U;
  setfield!(C, :t, vcat(getfield(C,:t), getfield(U,:t)));
  setfield!(C, :corr, hcat(getfield(C,:corr), getfield(U,:corr)))
  elseif isempty(C)
  [setfield!(C, i, getfield(U,i)) for i in corrfields];
  [setfield!(C, i, getfield(U,i)) for i in [:t,:corr]];
  end;
  return C)
+(C::CorrData, U::CorrData) = (T = deepcopy(C); return append!(T, U))

function sizeof(C::CorrData)
s = sum([sizeof(getfield(C,f)) for f in corrfields])
if !isempty(C.notes)
  s += sum([sizeof(i) for i in C.notes])
end
if !isempty(C.misc)
  s += sum([sizeof(i) for i in values(C.misc)])
end
return s
end


text:
### LombScargle.jl ---  Perform Lomb–Scargle periodogram
#
# Copyright (C) 2016, 2017 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: periodogram, lomb scargle
#
# This file is a part of LombScargle.jl.
#
# License is BSD 3-clause "New" or "Revised".
#
### Code:

__precompile__()

module LombScargle

using FFTW, Measurements, LinearAlgebra, Statistics, SpecialFunctions

export lombscargle

# This is similar to Periodogram type of DSP.Periodograms module, but for
# unevenly spaced frequencies.
struct Periodogram{P<:AbstractFloat, F<:AbstractVector, T<:AbstractVector}
    power::Vector{P}
    freq::F
    # XXX: the `times' vector is only in the `M' function (see utils.jl), but
    # only maximum(times) and minimum(times) are used.  We could consider the
    # possibility to keep in this type only the extrema of t, instead of the
    # whole array.
    times::T
    norm::Symbol
end

abstract type PeriodogramPlan end

include("utils.jl")
include("bootstrap.jl")

include("townsend.jl")
include("gls.jl")
include("press-rybicki.jl")

include("planning.jl")

function normalize!(P::AbstractVector{<:Real},
                    signal::AbstractVector{<:Real},
                    psdfactor::Real,
                    N::Integer,
                    noise_level::Real,
                    normalization::Symbol)
    # CLH edit: I've used a complex Lomb-Scargle plan with :psd normalization as default, but if you're
    # normalizing you're likely trying to compute the periodogram. The following
    # reverses the :psd normalization
    if normalization == :psd
        return P
    else
      P .*= 2 ./ psdfactor
      if normalization == :standard
          return P
      elseif normalization == :model
          return P ./= (1 .- P)
      elseif normalization == :log
          return P .= -log.(1 .- P)
      elseif normalization == :Scargle
          return P ./= noise_level
      elseif normalization == :HorneBaliunas
          return P .*= (N .- 1) ./ 2
      elseif normalization == :Cumming
          M = maximum(P)
          return P .*= (N .- 3) ./ (1 .- M) ./ 2
      else
          error("normalization \"", string(normalization), "\" not supported")
      end
    end
end

normalize(P::AbstractVector{<:Complex}, p::PeriodogramPlan) =
    normalize!(Float64.(real.(abs2.(P))), p.signal, p.YY, length(p.signal), p.noise, p.norm)

function lombscargle(p::PeriodogramPlan)
    return Periodogram(normalize(_periodogram!(p), p), p.freq, p.times, p.norm)
end

lombscargle(args...; kwargs...) = lombscargle(plan(args...; kwargs...))

"""
    lombscargle(times::AbstractVector{<:Real}, signal::AbstractVector{<:Real},
                [errors::AbstractVector{<:Real}]; keywords...)

Compute the Lomb–Scargle periodogram of the `signal` vector, observed at
`times`.  You can also specify the uncertainties for each signal point with
`errors` argument.  All these vectors must have the same length.

All optional keywords are described in the docstring of
[`LombScargle.plan`](@ref).

If the signal has uncertainties, the `signal` vector can also be a vector of
`Measurement` objects (from
[`Measurements.jl`](https://github.com/giordano/Measurements.jl) package), in
which case you don’t need to pass a separate `errors` vector for the
uncertainties of the signal.
"""
lombscargle(::AbstractVector{<:Real}, rest...)

"""
    lombscargle(plan::PeriodogramPlan)

Compute the Lomb–Scargle periodogram for the given `plan`.  This method has no
other arguments.  See documentation of [`LombScargle.plan`](@ref) for how to
plan a Lomb–Scargle periodogram.
"""
lombscargle(::PeriodogramPlan)

end # module


text:
using Plots
using Omega
import Omega.Space
using Omega.Space: getdim

"""
Contour Plot of two dimensions of Ω

```
x = normal(0.0, 1.0)
y = normal(0.0, 1.0)
c1 = err(x ==ₛ y)
c2 = err(x >ₛ y)
ωcontour(c2)
```

"""
function ωcontour(xrv::RandVar;
                  ΩT = defΩ(),
                  xdim = 1,
                  ydim = 2,
                  xrng = 0:0.005:1,
                  yrng = 0:0.005:1,
                  plt = plot(),
                  fill = true,
                  kwargs...)
  ω = ΩT()
  xrv(ω)
  function f(x, y)
    ω = Space.update(ω, xdim, x)
    ω = Space.update(ω, ydim, y)
    xrv(ω)
  end
  contour!(plt, xrng, yrng, f; fill = fill, kwargs...)
end

# "x ∈ [0, 1]?"
# isunit(x) = 0.0 <= x <= 1.0

"Plot path in Omega Space from sequence of samples"
function ωtrace(ωs::Vector;
                   xdim = 1,
                   ydim = 2,
                   plt = plot())
  xs, ys = ntranspose(map(ω -> (getdim(ω, xdim), getdim(ω, ydim)), ωs))
  plot!(plt, xs, ys, arrow = :arrow, linealpha = 0.5, legend = false)
end

# function plottraces(qpdata, plt = plot())
#   for qp in qpdata
#     plottrace(qp, plt)
#   end
#   plt
# end

# function testcb2(;alg = HMC, n = 1000, kwargs...)
#   μ = normal(0.0, 1.0)
#   x = normal(μ, 1.0)
#   y = (x ==ₛ 0.0)
#   # y = (x == 0.0) | (μ < 0.0)
#   cb, cbdata = Omega.tracecb(Omega.QP, deepcopy)
#   # cb = [default_cbs(n); cb]
#   rand(μ, y, n; alg = alg, cb = cb, kwargs...)
#   qpdata = cbdata[2]
#   plt = plot()
#   ωcontour(y, x.id, μ.id, Omega.defΩ(HMC)(), plt = plt)
#   print(plottraces(qpdata, plt))
#   qpdata, plt
# end

# testcb(nsteps = 20, stepsize = 0.01, n = 10000)

# function testcb2(;kwargs...)
#   μ = normal(0.0, 1.0)
#   x = normal(μ, 1.0)
#   # y = (x == 0.0) | (μ < 0.0)
#   y = (x^2 + μ^2 == 1.0)
#   cb, cbdata = Omega.tracecb(Omega.QP)
#   cb = [default_cbs(n); cb]
#   rand(μ, y, HMC; n = n, cb = cb, kwargs...)
#   qpdata = cbdata[2]
#   plt = plot()
#   ωcontour(y, x.id, μ.id, Omega.defΩ(HMC)(), plt = plt)
#   plottraces(qpdata, plt)
#   qpdata, plt
# end
# # qpdata, plt = testcb2(nsteps = 20, stepsize = 0.01, n=5000)

# function to_output_space(q1, q2)
#   n1  = Distributions.Normal(0, 1)
#   μ = quantile(n1, q1)
#   n2 = Distributions.Normal(μ, 1)
#   μ, quantile(n2, q2)
# end

# # function plotoutput(data, i=0,plt = plot())
# #   transf(d) = to_output_space(d.q |> Omega.inv_transform)...)[2]
# #   ys = map(transf, data)
# #   k = 1.0/(data |> length)
# #   xs = i:k:(i+1-k)
# #   plot!(plt, xs, ys, arrow = :arrow, linealpha = 0.5, legend=false)
# # end

# function plotoutputs(data, plt=plot())
#   for (i, qp) in enumerate(data)
#     plotoutput(qp, i, plt)
#   end
#   plt
# end

# function plotoutputXY(data, plt = plot())
#   xs, ys = zip([g((d.q |> Omega.inv_transform)...)
#                 for d in data]...) |> collect
#   plot!(plt, xs |> collect, ys |> collect, arrow = :arrow, linealpha = 0.5, legend=false)
# end

# function plotoutputsXY(data, plt=plot())
#   foreach(data) do qp
#     plotoutputXY(qp, plt)
#   end
#   plt
# end

# function filtered_qp(qpdata)
#   indexes = filter(1:length(qpdata)-1) do i
#     qpdata[i][1].q != qpdata[i+1][1].q
#   end
#   qpdata[indexes]
# end


text:
module ConnectSDK

using HTTPClient.HTTPC
using JSON
using Compat

export ConnectClient
export push

# package code goes here
type ConnectClient
  projectid
  apikey
  baseurl::AbstractString

  function ConnectClient(projectid::AbstractString, apikey::AbstractString)
	baseurl = "https://api.getconnect.io/events/"
	projectid = ("X-Project-Id", projectid)
	apikey = ("X-Api-Key", apikey)

	new(projectid, apikey, baseurl)
  end

end #ConnectClient

function push(client::ConnectClient, collection::AbstractString, data::Dict{Any,Any})
  data = JSON.json(data)
  r = HTTPC.post(client.baseurl * collection, data, RequestOptions(headers=[client.projectid, client.apikey]))
  return r
end

function query(client::ConnectClient, collection::AbstractString, query::AbstractString)

end

# example push/post url: https://api.getconnect.io/events/purchases where purchases is the collection

# example query url: https://api.getconnect.io/events/purchases?query={select:{totalPrice:{sum:'price'}}}

end # module