text:
export SatelliteSystemType

@enum SatelliteSystemType begin
    sstGPS
    sstGalileo
    sstGlonass
    sstGeosync
    sstLEO
    sstTransit
    sstBeiDou
    sstQZSS
    sstIRNSS
    sstMixed
    sstUserDefined
    sstUnknown
end

#  ObsID::map1to3sys\["(.)"\] = "(...)";
#  ObsID::map3to1sys\["(...)"\] = "(.)";
# $1=>$2,
const mapC2SS = Dict{Char,SatelliteSystemType}(
    'G' => sstGPS,
    'R' => sstGlonass,
    'E' => sstGalileo,
    'S' => sstGeosync,
    'C' => sstBeiDou,
    'J' => sstQZSS,
    'I' => sstIRNSS,
    'M' => sstMixed,
)

function char2SatelliteSystemType(c::Char,
    mapC2SS::Dict{Char,SatelliteSystemType}=mapC2SS,
)::SatelliteSystemType
    c = uppercase(c)
    if haskey(mapC2SS, c)
        mapC2SS[c]
    else
        print("$c>G")
        mapC2SS['G']
    end

end
const mapSS2C = Dict{SatelliteSystemType,Char}(
    sstGPS => 'G',
    sstGlonass => 'R',
    sstGalileo => 'E',
    sstGeosync => 'S',
    sstBeiDou => 'C',
    sstQZSS => 'J',
    sstIRNSS => 'I',
)
const mapSS2Str = Dict{SatelliteSystemType,String}(
    sstGPS => "GPS",
    sstGlonass => "GLO",
    sstGalileo => "GAL",
    sstGeosync => "GEO",
    sstBeiDou => "BDS",
    sstQZSS => "QZS",
    sstIRNSS => "IRN",
)

const validRinexFrequencies = "12356789";

# this defines the valid obs types
# NB these tc characters are ORDERED ~best to worst
# ObsID::validRinexTrackingCodes\[('.')\]\[('.')\] = (".+");
# $1 $2 $3
# $2 =>$3,
validRinexTrackingCodes = Dict{SatelliteSystemType,Dict{Char,String}}(
    sstGPS => Dict(
        '1' => "PYWLMIQSXCN* ",   # except no C1N
        '2' => "PYWLMIQSXCDN* ",  # except no C2N
        '5' => "IQX* ",),
    sstGlonass => Dict(
        '1' => "PC* ",
        '2' => "PC* ",
        '3' => "IQX* ",),

    sstGalileo => Dict(
        '1' => "ABCIQXZ* ",
        '5' => "IQX* ",
        '6' => "ABCIQXZ* ",
        '7' => "IQX* ",
        '8' => "IQX* ",),

    sstGeosync => Dict(
        '1' => "C* ",
        '5' => "IQX* ",),

    sstBeiDou => Dict(
        # NB 24Jun2013 MGEX data uses 2!  RINEX 3.03: 1 for 3.02, 2 for 3.0[013]
        '1' => "IQX* ",
        '2' => "IQX* ",
        '6' => "IQX* ",
        '7' => "IQX* ",),

    sstQZSS => Dict(
        '1' => "CSLXZ* ",
        '2' => "SLX* ",
        '5' => "IQX* ",
        '6' => "SLX* ",),

    sstIRNSS => Dict(
        '5' => "ABCX* ",
        '9' => "ABCX* ",),
)

export SatID
struct SatID
    str::String
    ss::SatelliteSystemType
    id::Int64
    function SatID(str::Union{String,SubString})
        length(str) == 3 || throw(error("not a valid satid: $str"))

        ss = char2SatelliteSystemType(str[1])
        if str[2] == ' '
            str = string(str[1], '0', str[3])
        end
        ids = str[2:3]
        id = parse(Int64, ids)
        new(str, ss, id)
    end
    function SatID(ss::SatelliteSystemType, id::Int64)
        str = string(mapSS2C[ss], @sprintf("%02d",id))
        new(str, ss, id)
    end
end
Base.show(io::IO, i::SatID) = print(io, i.str)


text:
using Test, DataFrameTools, DataFrames

include("io.jl")

ext: jl


text:
using Arrowhead, Base.Test

srand(0)

import Arrowhead: eigvals_FMM, synthesize_symarrow

@testset "Testing symmetric arrowhead methods" begin
    for n in (100, 300, 1000)

        A = SymArrow((1.0:n-1).*(2.0:n), n*randn(n-1), n*(n+1.0))

        @test issymmetric(A)

        B = SymArrow(ones(n-1), rand(n-1)./n^2, 1.0)

        @test issymmetric(B)

        R = chol(B)

        @test !issymmetric(R)

        @test norm(R'R-B)/norm(B) < n*eps()

        Ri = inv(R)

        @test !issymmetric(Ri)

        @test norm(Ri*R-I) < n*eps()

        SMA = Symmetric(Matrix(A))
        SMB = Symmetric(Matrix(B))

        λ = eigvals(A)
        Λ = eigvals(SMA)
        println("Testing eigenvalue error")
        println(norm(λ-Λ)/norm(Λ),"  ", 20*eps())
        @test norm(λ-Λ)/norm(Λ) < 20*eps()

        A1 = synthesize_symarrow(A, λ)
        Λ, Q = eig(A)

        println("Testing orthogonality")
        println(norm(Q'Q-I),"  ",20*eps())
        #@test norm(Q'Q-I) < 20*eps()

        println("Testing backward error")
        println(norm(A*Q-Q*Diagonal(Λ)),"  ",1.125*(norm(A.b-A1.b)+norm(A.c-A1.c)))
        @test norm(A*Q-Q*Diagonal(Λ)) < 1.125*(norm(A.b-A1.b)+norm(A.c-A1.c))

        λ = eigvals(A, B)
        Λ = eigvals(SMA, SMB)
        println("Testing generalized eigenvalue error")
        println(norm(λ-Λ)/norm(Λ),"  ", 20*eps())
        @test norm(λ-Λ)/norm(Λ) < 20*eps()
    end
end

@testset "Testing FMM acceleration" begin
    for n in (100, 300, 1000, 3000, 10000)
        println("n = ", n)
        A = SymArrow((1.0:n-1).*(2.0:n), n*ones(n-1), n*(n+1.0))
        if n ≤ 10_000
            @time λ = eigvals(A)
            @time λ_FMM = eigvals_FMM(A)

            println("Testing relative error")
            println(norm((λ-λ_FMM)./λ, Inf))
            @test norm((λ-λ_FMM)./λ, Inf) < 10*eps()
        else
            @time eigvals_FMM(A)
        end
    end
end


text:
module SHTOOLS

using SHTOOLS_jll

const Optional{T} = Union{Nothing,T}
# For optional arguments
optional(x::Optional, x0) = x !== nothing ? x : x0

include("legendre.jl")
include("transforms.jl")
include("storage.jl")
include("misc.jl")

end


text:
include("../src/g3.jl")

using PyPlot

@testset "g3_hn" begin

gamma = logspace(-10.0,0.0,1000)
gbig = big.(gamma)

for beta in [-1.5,-0.75,-0.1,-0.01,-0.001,-0.0001,-0.00001,-0.000001,0.000001,0.00001,0.0001,0.001,0.01,0.1,0.75,1.5]
  clf()
  betabig=big(beta)
  zbig = big(0.0)

  g3b = convert(Array{Float64,1},G3.(gbig,betabig;gc = zbig))
  h1b = convert(Array{Float64,1},H1.(gbig,betabig;gc = zbig))
  h2b = convert(Array{Float64,1},H2.(gbig,betabig;gc = zbig))
  h3b = convert(Array{Float64,1},H3.(gbig,betabig;gc = zbig))
  h5b = convert(Array{Float64,1},H5.(gbig,betabig;gc = zbig))
  h6b = convert(Array{Float64,1},H6.(gbig,betabig;gc = zbig))
  h7b = convert(Array{Float64,1},H7.(gbig,betabig;gc = zbig))
  h8b = convert(Array{Float64,1},H8.(gbig,betabig;gc = zbig))

  g3 = G3.(gamma,beta)
  h1 = H1.(gamma,beta)
  h2 = H2.(gamma,beta)
  h3 = H3.(gamma,beta)
  h5 = H5.(gamma,beta)
  h6 = H6.(gamma,beta)
  h7 = H7.(gamma,beta)
  h8 = H8.(gamma,beta)
  semilogy(gamma,abs.(g3./g3b-1),label="G3")
  plot(gamma,abs.(h1./h1b-1),label="H1")
  plot(gamma,abs.(h2./h2b-1),label="H2")
  plot(gamma,abs.(h3./h3b-1),label="H3")
  plot(gamma,abs.(h5./h5b-1),label="H5")
  plot(gamma,abs.(h6./h6b-1),label="H6")
  plot(gamma,abs.(h7./h7b-1),label="H7")
  plot(gamma,abs.(h8./h8b-1),label="H8")
  @test isapprox(g3b,g3;norm=maxabs)
  @test isapprox(h1b,h1;norm=maxabs)
  @test isapprox(h2b,h2;norm=maxabs)
  @test isapprox(h3b,h3;norm=maxabs)
  @test isapprox(h5b,h5;norm=maxabs)
  @test isapprox(h6b,h6;norm=maxabs)
  @test isapprox(h7b,h7;norm=maxabs)
  @test isapprox(h8b,h8;norm=maxabs)
  println("beta: ",beta)
  axis([0.0,1.0,1e-16,0.1])
  legend(loc="upper left")
  read(STDIN,Char)
end
end


text:
module HOFASM

using LinearAlgebra
using Combinatorics
using SparseArrays
using Statistics
using Random
using Arpack
using Distributed
using ImageFeatures  #used for demos

#
# Primary Experiment Drivers
#

export distributed_timing_experiments
export distributed_accuracy_experiments
export timing_experiments
export accuracy_experiments
export synthetic_HOM
export synthetic_HOFASM
export build_assignment

#
# Routines to help using our code in your own experiments
#

export align_photos
export align_embeddings

include("triangles.jl")
include("contraction.jl")
include("GraduatedAssignment.jl")
include("Experimental_code.jl")
include("HOFASM_impl.jl")
include("HOM_impl.jl")
include("experiments.jl")
#include("test/test_code.jl")

end # module


text:
"""
    unsafe_createcolortable(palette::GDALPaletteInterp)

Construct a new color table.
"""
unsafe_createcolortable(palette::GDALPaletteInterp)::ColorTable =
    ColorTable(GDAL.gdalcreatecolortable(palette))

"Destroys a color table."
function destroy(ct::ColorTable)::Nothing
    GDAL.gdaldestroycolortable(ct.ptr)
    ct.ptr = C_NULL
    return nothing
end

"""
    unsafe_clone(ct::ColorTable)

Make a copy of a color table.
"""
function unsafe_clone(ct::ColorTable)::ColorTable
    return if ct.ptr == C_NULL
        ColorTable(C_NULL)
    else
        ColorTable(GDAL.gdalclonecolortable(ct.ptr))
    end
end

"""
    paletteinterp(ct::ColorTable)

Fetch palette interpretation.

### Returns
palette interpretation enumeration value, usually `GPI_RGB`.
"""
paletteinterp(ct::ColorTable)::GDALPaletteInterp =
    GDAL.gdalgetpaletteinterpretation(ct.ptr)

"""
    getcolorentryasrgb(ct::ColorTable, i::Integer)

Fetch a table entry in RGB format.

In theory this method should support translation of color palettes in non-RGB
color spaces into RGB on the fly, but currently it only works on RGB color
tables.

### Parameters
* `i`   entry offset from zero to GetColorEntryCount()-1.

### Returns
`true` on success, or `false` if the conversion isn't supported.
"""
function getcolorentryasrgb(ct::ColorTable, i::Integer)::GDAL.GDALColorEntry
    colorentry = Ref{GDAL.GDALColorEntry}(GDAL.GDALColorEntry(0, 0, 0, 0))
    result = Bool(GDAL.gdalgetcolorentryasrgb(ct.ptr, i, colorentry))
    result || @warn("The conversion to RGB isn't supported.")
    return colorentry[]
end

"""
    setcolorentry!(ct::ColorTable, i::Integer, entry::GDAL.GDALColorEntry)

Set entry in color table.

Note that the passed in color entry is copied, and no internal reference to it
is maintained. Also, the passed in entry must match the color interpretation of
the table to which it is being assigned.

The table is grown as needed to hold the supplied offset.

### Parameters
* `i`     entry offset from `0` to `ncolorentry()-1`.
* `entry` value to assign to table.
"""
function setcolorentry!(ct::ColorTable, i::Integer, entry::GDAL.GDALColorEntry)
    GDAL.gdalsetcolorentry(ct.ptr, i, Ref{GDAL.GDALColorEntry}(entry))
    return ct
end


text:
# Linear ODE
linear = (u,p,t) -> (p*u)
linear_analytic = (u0,p,t) -> u0*exp(p*t)
"""
Linear ODE

```math
\\frac{du}{dt} = αu
```

with initial condition ``u0=1/2``, ``α=1.01``, and solution

```math
u(t) = u0e^{αt}
```

with Float64s. The parameter is ``α``
"""
prob_ode_linear = ODEProblem(
                  ODEFunction(linear,analytic=linear_analytic),
                  1/2,(0.0,1.0),1.01)

"""
Linear ODE

```math
\\frac{du}{dt} = αu
```

with initial condition ``u0=1/2``, ``α=1.01``, and solution

```math
u(t) = u0e^{αt}
```

with BigFloats
"""
prob_ode_bigfloatlinear = ODEProblem(
                          ODEFunction(linear,analytic=linear_analytic),
                          big(0.5),(0.0,1.0),big(1.01))

f_2dlinear = (du,u,p,t) -> (@. du = p*u)
f_2dlinear_analytic = (u0,p,t) -> @. u0*exp(p*t)
"""
4x2 version of the Linear ODE

```math
\\frac{du}{dt} = αu
```

with initial condition ``u0=1/2``, ``α=1.01``, and solution

```math
u(t) = u0e^{αt}
```

with Float64s
"""
prob_ode_2Dlinear = ODEProblem(
                    ODEFunction(f_2dlinear,analytic=f_2dlinear_analytic),
                    rand(4,2),(0.0,1.0),1.01)

"""
100x100 version of the Linear ODE

```math
\\frac{du}{dt} = αu
```

with initial condition ``u0=1/2``, ``α=1.01``, and solution

```math
u(t) = u0e^{αt}
```

with Float64s
"""
prob_ode_large2Dlinear = ODEProblem(
                    ODEFunction(f_2dlinear,analytic=f_2dlinear_analytic),
                    rand(100,100),(0.0,1.0),1.01)

"""
4x2 version of the Linear ODE

```math
\\frac{du}{dt} = αu
```

with initial condition ``u0=1/2``, ``α=1.01``, and solution

```math
u(t) = u0e^{αt}
```

with BigFloats
"""
prob_ode_bigfloat2Dlinear = ODEProblem(
                    ODEFunction(f_2dlinear,analytic=f_2dlinear_analytic),
                    BigFloat.(rand(4,2)).*ones(4,2)/2,(0.0,1.0),big(1.01))

f_2dlinear_notinplace = (u,p,t) -> p*u
"""
4x2 version of the Linear ODE

```math
\\frac{du}{dt} = αu
```

with initial condition ``u0=1/2``, ``α=1.01``, and solution

```math
u(t) = u0e^{αt}
```

on Float64. Purposefully not in-place as a test.
"""
prob_ode_2Dlinear_notinplace = ODEProblem(
                    ODEFunction(f_2dlinear_notinplace,analytic=f_2dlinear_analytic),
                    rand(4,2),(0.0,1.0),1.01)


text:
function hide(x) x end

for n = hide("foo") : 2
    println(n)
end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
# Tables and TableTraits interface

## Chains

function _check_columnnames(chn::Chains)
    for name in names(chn)
        symname = Symbol(name)
        if symname === :iteration || symname === :chain
            error("'$(name)' is a reserved column name. Please rename the parameter.")
        end
    end
end

Tables.istable(::Type{<:Chains}) = true

Tables.columnaccess(::Type{<:Chains}) = true

function Tables.columns(chn::Chains)
    _check_columnnames(chn)
    return chn
end

Tables.columnnames(chn::Chains) = (:iteration, :chain, names(chn)...)

function Tables.getcolumn(chn::Chains, i::Int)
    return Tables.getcolumn(chn, Tables.columnnames(chn)[i])
end
function Tables.getcolumn(chn::Chains, nm::Symbol)
    if nm == :iteration
        iterations = range(chn)
        nchains = size(chn, 3)
        return repeat(iterations, nchains)
    elseif nm == :chain
        chainids = chains(chn)
        niter = size(chn, 1)
        return repeat(chainids; inner = niter)
    else
        return vec(chn[nm])
    end
end

Tables.rowaccess(::Type{<:Chains}) = true

function Tables.rows(chn::Chains)
    _check_columnnames(chn)
    return chn
end

Tables.rowtable(chn::Chains) = Tables.rowtable(Tables.columntable(chn))

function Tables.namedtupleiterator(chn::Chains)
    return Tables.namedtupleiterator(Tables.columntable(chn))
end

function Tables.schema(chn::Chains)
    _check_columnnames(chn)
    nms = Tables.columnnames(chn)
    T = eltype(chn.value)
    types = (Int, Int, ntuple(_ -> T, size(chn, 2))...)
    return Tables.Schema(nms, types)
end

IteratorInterfaceExtensions.isiterable(::Chains) = true
function IteratorInterfaceExtensions.getiterator(chn::Chains)
    return Tables.datavaluerows(Tables.columntable(chn))
end

TableTraits.isiterabletable(::Chains) = true

## ChainDataFrame

Tables.istable(::Type{<:ChainDataFrame}) = true

Tables.columnaccess(::Type{<:ChainDataFrame}) = true

Tables.columns(cdf::ChainDataFrame) = cdf

Tables.columnnames(::ChainDataFrame{<:NamedTuple{names}}) where {names} = names

Tables.getcolumn(cdf::ChainDataFrame, i::Int) = cdf.nt[i]
Tables.getcolumn(cdf::ChainDataFrame, nm::Symbol) = cdf.nt[nm]

Tables.rowaccess(::Type{<:ChainDataFrame}) = true

Tables.rows(cdf::ChainDataFrame) = cdf

Tables.rowtable(cdf::ChainDataFrame) = Tables.rowtable(Tables.columntable(cdf))

function Tables.namedtupleiterator(cdf::ChainDataFrame)
    return Tables.namedtupleiterator(Tables.columntable(cdf))
end

function Tables.schema(::ChainDataFrame{NamedTuple{names,T}}) where {names,T}
    types = ntuple(i -> eltype(fieldtype(T, i)), fieldcount(T))
    return Tables.Schema(names, types)
end

IteratorInterfaceExtensions.isiterable(::ChainDataFrame) = true
function IteratorInterfaceExtensions.getiterator(cdf::ChainDataFrame)
    return Tables.datavaluerows(Tables.columntable(cdf))
end

TableTraits.isiterabletable(::ChainDataFrame) = true


text:
using Documenter, CuArrays

makedocs(
    modules = [CuArrays],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "CuArrays.jl",
    pages = [
        "Home" => "index.md",
    ],
    doctest = true
)


text:
###
### Gibbs samplers / compositional samplers.
###

"""
    Gibbs(algs...)

Compositional MCMC interface. Gibbs sampling combines one or more
sampling algorithms, each of which samples from a different set of
variables in a model.

Example:
```julia
@model gibbs_example(x) = begin
    v1 ~ Normal(0,1)
    v2 ~ Categorical(5)
end
```

# Use PG for a 'v2' variable, and use HMC for the 'v1' variable.
# Note that v2 is discrete, so the PG sampler is more appropriate
# than is HMC.
alg = Gibbs(1000, HMC(1, 0.2, 3, :v1), PG(20, 1, :v2))
```

Tips:
- `HMC` and `NUTS` are fast samplers, and can throw off particle-based
methods like Particle Gibbs. You can increase the effectiveness of particle sampling by including
more particles in the particle sampler.
"""
mutable struct Gibbs{A} <: InferenceAlgorithm
    algs      ::  A   # component sampling algorithms
    function Gibbs(algs...)
        return new{typeof(algs)}(algs)
    end
end

alg_str(::Sampler{<:Gibbs}) = "Gibbs"
transition_type(spl::Sampler{<:Gibbs}) = typeof(Transition(spl))

"""
    GibbsState{V<:VarInfo, S<:Tuple{Vararg{Sampler}}}

Stores a `VarInfo` for use in sampling, and a `Tuple` of `Samplers` that
the `Gibbs` sampler iterates through for each `step!`.
"""
mutable struct GibbsState{V<:VarInfo, S<:Tuple{Vararg{Sampler}}} <: AbstractSamplerState
    vi::V
    samplers::S
end

function GibbsState(model::Model, samplers::S) where S<:Tuple{Vararg{Sampler}}
    return GibbsState(VarInfo(model), samplers)
end

const GibbsComponent = Union{Hamiltonian,MH,PG}

function Sampler(alg::Gibbs, model::Model, s::Selector)
    info = Dict{Symbol, Any}()

    n_samplers = length(alg.algs)
    samplers = Array{Sampler}(undef, n_samplers)
    space = Set{Symbol}()

    for i in 1:n_samplers
        sub_alg = alg.algs[i]
        if isa(sub_alg, GibbsComponent)
            samplers[i] = Sampler(sub_alg, model, Selector(Symbol(typeof(sub_alg))))
        else
            @error("[Gibbs] Unsupported sampling algorithm $sub_alg")
        end
        space = (space..., getspace(sub_alg)...)
    end

    # Sanity check for space
    @assert issubset(get_pvars(model), space) "[Gibbs] symbols specified to samplers ($space) doesn't cover the model parameters ($(get_pvars(model)))"

    if !(issetequal(get_pvars(model), space))
        @warn("[Gibbs] extra parameters specified by samplers don't exist in model: $(setdiff(space, get_pvars(model)))")
    end

    # Create a state variable.
    state = GibbsState(model, tuple(samplers...))

    # Create the sampler.
    spl = Sampler(alg, info, s, state)

    # Add Gibbs to gids for all variables.
    for sym in keys(spl.state.vi.metadata)
        vns = getfield(spl.state.vi.metadata, sym).vns
        for vn in vns
            # Update the gid for the Gibbs sampler.
            Turing.RandomVariables.updategid!(spl.state.vi, vn, spl)

            # Try to store each subsampler's gid in the VarInfo.
            for local_spl in spl.state.samplers
                Turing.RandomVariables.updategid!(spl.state.vi, vn, local_spl)
            end
        end
    end

    return spl
end

# Initialize the Gibbs sampler.
function sample_init!(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:Gibbs},
    N::Integer;
    kwargs...
)
    # Initialize each local sampler.
    for local_spl in spl.state.samplers
        sample_init!(rng, model, local_spl, N; kwargs...)
    end
end

# Finalize the Gibbs sampler.
function sample_end!(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:Gibbs},
    N::Integer;
    kwargs...
)
    # Finalize each local sampler.
    for local_spl in spl.state.samplers
        sample_end!(rng, model, local_spl, N; kwargs...)
    end
end

# First step.
function step!(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:Gibbs},
    N::Integer;
    kwargs...
)
    Turing.DEBUG && @debug "Gibbs stepping..."

    time_elapsed = 0.0
    lp = nothing; ϵ = nothing; eval_num = nothing

    # Iterate through each of the samplers.
    for local_spl in spl.state.samplers
        Turing.DEBUG && @debug "$(typeof(local_spl)) stepping..."

        Turing.DEBUG && @debug "recording old θ..."

        # Update the sampler's VarInfo.
        local_spl.state.vi = spl.state.vi

        # Step through the local sampler.
        time_elapsed_thin =
            @elapsed step!(rng, model, local_spl, N; kwargs...)

        # After the step, update the master varinfo.
        spl.state.vi = local_spl.state.vi

        # Uncomment when developing thinning functionality.
        # Retrieve symbol to store this subsample.
        # symbol_id = Symbol(local_spl.selector.gid)

        # # Store the subsample.
        # spl.state.subsamples[symbol_id][] = trans

        # Record elapsed time.
        time_elapsed += time_elapsed_thin
    end

    return Transition(spl)
end

# Steps 2:N
function step!(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:Gibbs},
    N::Integer,
    t::TransitionType;
    kwargs...
) where TransitionType<:AbstractTransition
    Turing.DEBUG && @debug "Gibbs stepping..."

    time_elapsed = 0.0
    lp = nothing
    ϵ = nothing
    eval_num = nothing

    # Iterate through each of the samplers.
    for local_spl in spl.state.samplers
        Turing.DEBUG && @debug "$(typeof(local_spl)) stepping..."

        Turing.DEBUG && @debug "recording old θ..."

        # Update the sampler's VarInfo.
        local_spl.state.vi = spl.state.vi

        # Step through the local sampler.
        time_elapsed_thin =
            @elapsed trans = step!(rng, model, local_spl, N, t; kwargs...)

        # After the step, update the master varinfo.
        spl.state.vi = local_spl.state.vi

        # Uncomment when developing thinning functionality.
        # Retrieve symbol to store this subsample.
        # symbol_id = Symbol(local_spl.selector.gid)
        #
        # # Store the subsample.
        # spl.state.subsamples[symbol_id][] = trans

        # Record elapsed time.
        time_elapsed += time_elapsed_thin
    end

    return Transition(spl)
end


text:
abstract type AbstractBlochMatrix{T<:Real} end

"""
    BlochMatrix(a11, a21, a31, a12, a22, a32, a13, a23, a33)
    BlochMatrix{T}()
    BlochMatrix()

Create a mutable `BlochMatrix` object representing a fixed-size 3×3 matrix.

# Properties
- `aij::Real`: Matrix entry (i,j) ∈ {1, 2, 3}²
"""
mutable struct BlochMatrix{T<:Real} <: AbstractBlochMatrix{T}
    a11::T
    a21::T
    a31::T
    a12::T
    a22::T
    a32::T
    a13::T
    a23::T
    a33::T
end

BlochMatrix(a...) = BlochMatrix(promote(a...)...)
BlochMatrix{T}() where {T} = BlochMatrix(zero(T), zero(T), zero(T), zero(T),
                                    zero(T), zero(T), zero(T), zero(T), zero(T))
BlochMatrix() = BlochMatrix{Float64}()

function Base.show(io::IO, ::MIME"text/plain", A::BlochMatrix{T}) where {T}

    print(io, "BlochMatrix{$T}:\n")
    haskey(io, :compact) || (io = IOContext(io, :compact => true))
    Base.print_array(io, Matrix(A))

end

function Base.fill!(A::BlochMatrix, v)

    A.a11 = v
    A.a21 = v
    A.a31 = v
    A.a12 = v
    A.a22 = v
    A.a32 = v
    A.a13 = v
    A.a23 = v
    A.a33 = v
    return nothing

end

function Base.copyto!(dst::BlochMatrix, src::BlochMatrix)

    dst.a11 = src.a11
    dst.a21 = src.a21
    dst.a31 = src.a31
    dst.a12 = src.a12
    dst.a22 = src.a22
    dst.a32 = src.a32
    dst.a13 = src.a13
    dst.a23 = src.a23
    dst.a33 = src.a33
    return nothing

end

function Base.Matrix(A::BlochMatrix{T}) where {T}

    mat = Matrix{T}(undef, 3, 3)
    mat[1,1] = A.a11
    mat[2,1] = A.a21
    mat[3,1] = A.a31
    mat[1,2] = A.a12
    mat[2,2] = A.a22
    mat[3,2] = A.a32
    mat[1,3] = A.a13
    mat[2,3] = A.a23
    mat[3,3] = A.a33

    return mat

end

"""
    BlochDynamicsMatrix(R1, R2, Δω)
    BlochDynamicsMatrix{T}()
    BlochDynamicsMatrix()

Create a mutable `BlochDynamicsMatrix` object.

# Properties
- `R1::Real`: Spin-lattice relaxation rate
- `R2::Real`: Spin-spin relaxation rate
- `Δω::Real`: Off-resonance frequency
"""
mutable struct BlochDynamicsMatrix{T<:Real} <: AbstractBlochMatrix{T}
    R1::T
    R2::T
    Δω::T
end

BlochDynamicsMatrix{T}() where {T} = BlochDynamicsMatrix(zero(T), zero(T), zero(T))
BlochDynamicsMatrix() = BlochDynamicsMatrix{Float64}()
BlochDynamicsMatrix(R1, R2, Δω) = BlochDynamicsMatrix(promote(R1, R2, Δω)...)

Base.convert(::Type{BlochDynamicsMatrix{T}}, A::BlochDynamicsMatrix) where {T} =
    BlochDynamicsMatrix(T(A.R1), T(A.R2), T(A.Δω))

function Base.Matrix(A::BlochDynamicsMatrix{T}) where {T}

    mat = Matrix{T}(undef, 3, 3)
    mat[1,1] = A.R2
    mat[2,1] = -A.Δω
    mat[3,1] = zero(T)
    mat[1,2] = A.Δω
    mat[2,2] = A.R2
    mat[3,2] = zero(T)
    mat[1,3] = zero(T)
    mat[2,3] = zero(T)
    mat[3,3] = A.R1

    return mat

end

"""
    FreePrecessionMatrix(E1, E2cosθ, E2sinθ)
    FreePrecessionMatrix{T}()
    FreePrecessionMatrix()

Create a mutable `FreePrecessionMatrix` object encoding the effects of
relaxation and off-resonance precession.

# Properties
- `E1::Real`: T1 relaxation
- `E2cosθ::Real`: T2 relaxation and off-resonance precession
- `E2sinθ::Real`: T2 relaxation and off-resonance precession

# Examples
```jldoctest
julia> T1 = 1000; T2 = 100; Δω = π/600; t = 100;

julia> F = FreePrecessionMatrix(exp(-t / T1), exp(-t / T2) * cos(Δω * t), exp(-t / T2) * sin(Δω * t))
FreePrecessionMatrix{Float64}:
 E1 = 0.9048374180359595
 E2cosθ = 0.31859294158449203
 E2sinθ = 0.18393972058572114

julia> Matrix(F)
3×3 Matrix{Float64}:
  0.318593  0.18394   0.0
 -0.18394   0.318593  0.0
  0.0       0.0       0.904837
```
"""
mutable struct FreePrecessionMatrix{T<:Real} <: AbstractBlochMatrix{T}
    E1::T
    E2cosθ::T
    E2sinθ::T
end

FreePrecessionMatrix{T}() where {T} = FreePrecessionMatrix(zero(T), zero(T), zero(T))
FreePrecessionMatrix() = FreePrecessionMatrix{Float64}()
FreePrecessionMatrix(E1, E2cosθ, E2sinθ) = FreePrecessionMatrix(promote(E1, E2cosθ, E2sinθ)...)

function Base.show(io::IO, ::MIME"text/plain", A::FreePrecessionMatrix{T}) where {T}

    print(io, "FreePrecessionMatrix{$T}:")
    print(io, "\n E1 = ", A.E1)
    print(io, "\n E2cosθ = ", A.E2cosθ)
    print(io, "\n E2sinθ = ", A.E2sinθ)

end

function Base.Matrix(A::FreePrecessionMatrix{T}) where {T}

    mat = Matrix{T}(undef, 3, 3)
    mat[1,1] = A.E2cosθ
    mat[2,1] = -A.E2sinθ
    mat[3,1] = zero(T)
    mat[1,2] = A.E2sinθ
    mat[2,2] = A.E2cosθ
    mat[3,2] = zero(T)
    mat[1,3] = zero(T)
    mat[2,3] = zero(T)
    mat[3,3] = A.E1

    return mat

end

function Base.copyto!(dst::BlochMatrix{T}, src::FreePrecessionMatrix) where {T}

    dst.a11 = src.E2cosθ
    dst.a21 = -src.E2sinθ
    dst.a31 = zero(T)
    dst.a12 = src.E2sinθ
    dst.a22 = src.E2cosθ
    dst.a32 = zero(T)
    dst.a13 = zero(T)
    dst.a23 = zero(T)
    dst.a33 = src.E1
    return nothing

end

"""
    ExchangeDynamicsMatrix(r)
    ExchangeDynamicsMatrix{T}()
    ExchangeDynamicsMatrix()

Create a mutable `ExchangeDynamicsMatrix` object.

# Properties
- `r::Real`: Exchange rate from one compartment to another
"""
mutable struct ExchangeDynamicsMatrix{T<:Real} <: AbstractBlochMatrix{T}
    r::T
end

ExchangeDynamicsMatrix{T}() where {T} = ExchangeDynamicsMatrix(zero(T))
ExchangeDynamicsMatrix() = ExchangeDynamicsMatrix{Float64}()

Base.convert(::Type{ExchangeDynamicsMatrix{T}}, A::ExchangeDynamicsMatrix) where {T} =
    ExchangeDynamicsMatrix(T(A.r))

function Base.Matrix(A::ExchangeDynamicsMatrix{T}) where {T}

    mat = Matrix{T}(undef, 3, 3)
    mat[1,1] = A.r
    mat[2,1] = zero(T)
    mat[3,1] = zero(T)
    mat[1,2] = zero(T)
    mat[2,2] = A.r
    mat[3,2] = zero(T)
    mat[1,3] = zero(T)
    mat[2,3] = zero(T)
    mat[3,3] = A.r

    return mat

end

abstract type AbstractBlochMcConnellMatrix{T<:Real,N} end

"""
    BlochMcConnellDynamicsMatrix(A, E)
    BlochMcConnellDynamicsMatrix{T}(N)
    BlochMcConnellDynamicsMatrix(N)

Create a `BlochMcConnellDynamicsMatrix` object with `N` compartments.

# Properties
- `A::NTuple{N,BlochDynamicsMatrix{<:Real}}`: List of `BlochDynamicsMatrix`es
  that make up the main block diagonal of the `BlochMcConnellDynamicsMatrix`
- `E::NTuple{M,ExchangeDynamicsMatrix{<:Real}}`: List of
  `ExchangeDynamicsMatrix`es that describe exchange between the different
  compartments; these matrices make up the remaining `M = N^2 - N` blocks of the
  `BlochMcConnellDynamicsMatrix`, sorted by column-major ordering
"""
struct BlochMcConnellDynamicsMatrix{T<:Real,N,M} <: AbstractBlochMcConnellMatrix{T,N}
    A::NTuple{N,BlochDynamicsMatrix{T}}
    E::NTuple{M,ExchangeDynamicsMatrix{T}}

    function BlochMcConnellDynamicsMatrix(
        A::NTuple{N,BlochDynamicsMatrix{T1}},
        E::NTuple{M,ExchangeDynamicsMatrix{T2}}
    ) where {M,N,T1,T2}

        M == N * (N - 1) || error("exchange rates must be defined for each pair of compartments")
        Tnew = promote_type(T1, T2)
        A = convert.(BlochDynamicsMatrix{Tnew}, A)
        E = convert.(ExchangeDynamicsMatrix{Tnew}, E)
        new{Tnew,N,M}(A, E)

    end
end

function BlochMcConnellDynamicsMatrix{T}(N) where {T}

    A = ntuple(i -> BlochDynamicsMatrix{T}(), N)
    E = ntuple(i -> ExchangeDynamicsMatrix{T}(), N * (N - 1))
    BlochMcConnellDynamicsMatrix(A, E)

end

BlochMcConnellDynamicsMatrix(N) = BlochMcConnellDynamicsMatrix{Float64}(N)

# N = 3
#   [2,1]: (1-1)*3-1+2+(2<1) = 0-1+2+0 = 1 ✓
#   [3,1]: (1-1)*3-1+3+(3<1) = 0-1+3+0 = 2 ✓
#   [1,2]: (2-1)*3-2+1+(1<2) = 3-2+1+1 = 3 ✓
#   [3,2]: (2-1)*3-2+3+(3<2) = 3-2+3+0 = 4 ✓
#   [1,3]: (3-1)*3-3+1+(1<3) = 6-3+1+1 = 5 ✓
#   [2,3]: (3-1)*3-3+2+(2<3) = 6-3+2+1 = 6 ✓
function getblock(A::BlochMcConnellDynamicsMatrix{T,N,M}, i, j) where {T,N,M}

    return i == j ? A.A[i] : A.E[(j-1)*N-j+i+(i<j)]

end

function Base.Matrix(A::BlochMcConnellDynamicsMatrix{T,N,M}) where {T,N,M}

    mat = Matrix{T}(undef, 3N, 3N)
    index = 0
    for j = 1:N, i = 1:N

        if i == j

            mat[3i-2,3j-2] = A.A[i].R2
            mat[3i-1,3j-2] = -A.A[i].Δω
            mat[3i  ,3j-2] = zero(T)
            mat[3i-2,3j-1] = A.A[i].Δω
            mat[3i-1,3j-1] = A.A[i].R2
            mat[3i  ,3j-1] = zero(T)
            mat[3i-2,3j]   = zero(T)
            mat[3i-1,3j]   = zero(T)
            mat[3i  ,3j]   = A.A[i].R1

        else

            index += 1
            mat[3i-2,3j-2] = A.E[index].r
            mat[3i-1,3j-2] = zero(T)
            mat[3i  ,3j-2] = zero(T)
            mat[3i-2,3j-1] = zero(T)
            mat[3i-1,3j-1] = A.E[index].r
            mat[3i  ,3j-1] = zero(T)
            mat[3i-2,3j]   = zero(T)
            mat[3i-1,3j]   = zero(T)
            mat[3i  ,3j]   = A.E[index].r

        end

    end

    return mat

end

"""
    BlochMcConnellMatrix(A)
    BlochMcConnellMatrix{T}(N)
    BlochMcConnellMatrix(N)

Create a `BlochMcConnellMatrix` object with `N` compartments and representing a
fixed-size 3N×3N matrix.

# Properties
- `A::NTuple{N,NTuple{N,BlochMatrix{<:Real}}}`: List of 3×3 matrices that
  comprise the blocks of the `BlochMcConnellMatrix`; `A[i][j]` is the (i,j)th
  block

# Examples
```jldoctest
julia> B = BlochMcConnellMatrix(3)
BlochMcConnellMatrix{Float64,3}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

julia> fill!(B, 0.123456789)

julia> B
BlochMcConnellMatrix{Float64,3}:
 0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457
 0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457
 0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457
 0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457
 0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457
 0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457
 0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457
 0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457
 0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457  0.123457
```
"""
struct BlochMcConnellMatrix{T<:Real,N} <: AbstractBlochMcConnellMatrix{T,N}
    A::NTuple{N,NTuple{N,BlochMatrix{T}}}
end

function BlochMcConnellMatrix{T}(N) where {T}

    A = ntuple(i -> ntuple(i -> BlochMatrix{T}(), N), N)
    BlochMcConnellMatrix(A)

end

BlochMcConnellMatrix(N) = BlochMcConnellMatrix{Float64}(N)

function Base.show(io::IO, ::MIME"text/plain", A::BlochMcConnellMatrix{T,N}) where {T,N}

    print(io, "BlochMcConnellMatrix{$T,$N}:\n")
    haskey(io, :compact) || (io = IOContext(io, :compact => true))
    Base.print_array(io, Matrix(A))

end

getblock(A::BlochMcConnellMatrix, i, j) = A.A[i][j]

function Base.fill!(A::BlochMcConnellMatrix{T,N}, v) where {T,N}

    for i = 1:N, j = 1:N
        fill!(getblock(A, i, j), v)
    end

end

function Base.copyto!(dst::BlochMcConnellMatrix{T,N}, src::BlochMcConnellMatrix{S,N}) where {S,T,N}

    for j = 1:N, i = 1:N
        db = getblock(dst, i, j)
        sb = getblock(src, i, j)
        db.a11 = sb.a11
        db.a21 = sb.a21
        db.a31 = sb.a31
        db.a12 = sb.a12
        db.a22 = sb.a22
        db.a32 = sb.a32
        db.a13 = sb.a13
        db.a23 = sb.a23
        db.a33 = sb.a33
    end

end

# Used when computing matrix exponential
function Base.copyto!(dst::BlochMcConnellMatrix{T,N}, src::AbstractMatrix) where {T,N}

    for j = 1:N, i = 1:N
        b = getblock(dst, i, j)
        b.a11 = src[3i-2,3j-2]
        b.a21 = src[3i-1,3j-2]
        b.a31 = src[3i  ,3j-2]
        b.a12 = src[3i-2,3j-1]
        b.a22 = src[3i-1,3j-1]
        b.a32 = src[3i  ,3j-1]
        b.a13 = src[3i-2,3j]
        b.a23 = src[3i-1,3j]
        b.a33 = src[3i  ,3j]
    end

end

function Base.Matrix(A::BlochMcConnellMatrix{T,N}) where {T,N}

    mat = Matrix{T}(undef, 3N, 3N)
    for j = 1:N, i = 1:N

        b = getblock(A, i, j)
        mat[3i-2,3j-2] = b.a11
        mat[3i-1,3j-2] = b.a21
        mat[3i  ,3j-2] = b.a31
        mat[3i-2,3j-1] = b.a12
        mat[3i-1,3j-1] = b.a22
        mat[3i  ,3j-1] = b.a32
        mat[3i-2,3j]   = b.a13
        mat[3i-1,3j]   = b.a23
        mat[3i  ,3j]   = b.a33

    end

    return mat

end

function Base.eltype(
    ::Union{<:AbstractBlochMatrix{T},<:AbstractBlochMcConnellMatrix{T,N}}
) where {T,N}

    return T

end

function Base.:(==)(
    A::Union{<:AbstractBlochMatrix,<:AbstractBlochMcConnellMatrix},
    B::Union{<:AbstractBlochMatrix,<:AbstractBlochMcConnellMatrix}
)

    return Matrix(A) == Matrix(B)

end

function Base.isapprox(
    A::Union{<:AbstractBlochMatrix,<:AbstractBlochMcConnellMatrix},
    B::Union{<:AbstractBlochMatrix,<:AbstractBlochMcConnellMatrix};
    kwargs...
)

    return isapprox(Matrix(A), Matrix(B); kwargs...)

end

"""
    ExcitationMatrix(A)
    ExcitationMatrix{T}()
    ExcitationMatrix()

Create an `ExcitationMatrix` object. Multiplying by a `MagnetizationMC` object
has the effect of multiplying each component of the multi-compartment
magnetization by the `ExcitationMatrix`.

# Properties
- `A::BlochMatrix{<:Real}`: Matrix used to describe the excitation

# Examples
```jldoctest
julia> E = ExcitationMatrix(BlochMatrix(0, 1, 0, 1, 0, 0, 0, 0, 1))
ExcitationMatrix{Int64}:
 0  1  0
 1  0  0
 0  0  1

julia> M = MagnetizationMC((1, 2, 3), (4, 5, 6))
2-compartment Magnetization vector with eltype Int64:
 Compartment 1:
  Mx = 1
  My = 2
  Mz = 3
 Compartment 2:
  Mx = 4
  My = 5
  Mz = 6

julia> E * M
2-compartment Magnetization vector with eltype Int64:
 Compartment 1:
  Mx = 2
  My = 1
  Mz = 3
 Compartment 2:
  Mx = 5
  My = 4
  Mz = 6
```
"""
struct ExcitationMatrix{T<:Real}
    A::BlochMatrix{T}
end

ExcitationMatrix{T}() where {T} = ExcitationMatrix(BlochMatrix{T}())
ExcitationMatrix() = ExcitationMatrix(BlochMatrix())

function Base.show(io::IO, ::MIME"text/plain", A::ExcitationMatrix{T}) where {T}

    print(io, "ExcitationMatrix{$T}:\n")
    haskey(io, :compact) || (io = IOContext(io, :compact => true))
    Base.print_array(io, Matrix(A.A))

end

"""
    idealspoiling = IdealSpoilingMatrix()

Matrix representing ideal spoiling. Multiplying by a `Magnetization` or
`MagnetizationMC` has the effect of setting the x and y components to 0.

# Examples
```jldoctest
julia> idealspoiling * MagnetizationMC((1, 1, 1), (2, 2, 2))
2-compartment Magnetization vector with eltype Int64:
 Compartment 1:
  Mx = 0
  My = 0
  Mz = 1
 Compartment 2:
  Mx = 0
  My = 0
  Mz = 2
```
"""
struct IdealSpoilingMatrix end
const idealspoiling = IdealSpoilingMatrix()

# ==============================================================================
# Math operations
# ==============================================================================

# +
Base.:+(M1::Magnetization, M2::Magnetization) = Magnetization(M1.x + M2.x, M1.y + M2.y, M1.z + M2.z)
Base.:+(M1::MagnetizationMC{T1,N}, M2::MagnetizationMC{T2,N}) where {T1,T2,N} = MagnetizationMC(ntuple(i -> M1[i] + M2[i], N)...)

# -
Base.:-(M1::Magnetization, M2::Magnetization) = Magnetization(M1.x - M2.x, M1.y - M2.y, M1.z - M2.z)
Base.:-(M1::MagnetizationMC{T1,N}, M2::MagnetizationMC{T2,N}) where {T1,T2,N} = MagnetizationMC(ntuple(i -> M1[i] - M2[i], N)...)

function Base.:-(::UniformScaling, B::BlochMatrix{T}) where {T}

    return BlochMatrix(
        one(T) - B.a11,
        -B.a21,
        -B.a31,
        -B.a12,
        one(T) - B.a22,
        -B.a32,
        -B.a13,
        -B.a23,
        one(T) - B.a33
    )

end

# *
Base.:*(M::Magnetization, a) = Magnetization(M.x * a, M.y * a, M.z * a)
Base.:*(M::MagnetizationMC{T,N}, a) where {T,N} = MagnetizationMC(ntuple(i -> M[i] * a, N)...)

function Base.:*(A::BlochMatrix, M::Magnetization)

    return Magnetization(
        A.a11 * M.x + A.a12 * M.y + A.a13 * M.z,
        A.a21 * M.x + A.a22 * M.y + A.a23 * M.z,
        A.a31 * M.x + A.a32 * M.y + A.a33 * M.z
    )

end

function Base.:*(A::FreePrecessionMatrix, M::Magnetization)

    return Magnetization(
        A.E2cosθ * M.x + A.E2sinθ * M.y,
        A.E2cosθ * M.y - A.E2sinθ * M.x,
        A.E1 * M.z
    )

end

Base.:*(A::ExcitationMatrix, M::Magnetization) = A.A * M
Base.:*(::IdealSpoilingMatrix, M::Magnetization{T}) where {T} = Magnetization(zero(T), zero(T), M.z)

function Base.:*(A::BlochMcConnellMatrix{T1,N}, M::MagnetizationMC{T2,N}) where {T1,T2,N}

    return MagnetizationMC(ntuple(N) do i
        Mc = getblock(A, i, 1) * M[1]
        for j = 2:N
            muladd!(Mc, getblock(A, i, j), M[j])
        end
        Mc
    end...)

end

Base.:*(A::ExcitationMatrix, M::MagnetizationMC) = MagnetizationMC((A * Mc for Mc in M)...)
Base.:*(A::IdealSpoilingMatrix, M::MagnetizationMC) = MagnetizationMC((A * Mc for Mc in M)...)

function Base.:*(A::BlochMatrix, B::BlochMatrix)

    return BlochMatrix(
        A.a11 * B.a11 + A.a12 * B.a21 + A.a13 * B.a31,
        A.a21 * B.a11 + A.a22 * B.a21 + A.a23 * B.a31,
        A.a31 * B.a11 + A.a32 * B.a21 + A.a33 * B.a31,
        A.a11 * B.a12 + A.a12 * B.a22 + A.a13 * B.a32,
        A.a21 * B.a12 + A.a22 * B.a22 + A.a23 * B.a32,
        A.a31 * B.a12 + A.a32 * B.a22 + A.a33 * B.a32,
        A.a11 * B.a13 + A.a12 * B.a23 + A.a13 * B.a33,
        A.a21 * B.a13 + A.a22 * B.a23 + A.a23 * B.a33,
        A.a31 * B.a13 + A.a32 * B.a23 + A.a33 * B.a33
    )

end

function Base.:*(A::BlochMatrix, B::FreePrecessionMatrix)

    return BlochMatrix(
        A.a11 * B.E2cosθ - A.a12 * B.E2sinθ,
        A.a21 * B.E2cosθ - A.a22 * B.E2sinθ,
        A.a31 * B.E2cosθ - A.a32 * B.E2sinθ,
        A.a11 * B.E2sinθ + A.a12 * B.E2cosθ,
        A.a21 * B.E2sinθ + A.a22 * B.E2cosθ,
        A.a31 * B.E2sinθ + A.a32 * B.E2cosθ,
        A.a13 * B.E1,
        A.a23 * B.E1,
        A.a33 * B.E1
    )

end

function Base.:*(A::BlochMatrix, ::IdealSpoilingMatrix)

    T = eltype(A)
    return BlochMatrix(zero(T), zero(T), zero(T), zero(T), zero(T), zero(T), A.a13, A.a23, A.a33)

end

Base.:*(A::ExcitationMatrix, ::IdealSpoilingMatrix) = A.A * idealspoiling

function Base.:*(A::BlochMcConnellMatrix{T1,N}, B::BlochMcConnellMatrix{T2,N}) where {T1,T2,N}

    C = BlochMcConnellMatrix{promote_type(T1, T2)}(N)
    mul!(C, A, B)
    return C

end

function Base.:*(A::ExcitationMatrix{T1}, B::BlochMcConnellMatrix{T2,N}) where {T1,T2,N}

    C = BlochMcConnellMatrix{promote_type(T1, T2)}(N)
    mul!(C, A, B)
    return C

end

function Base.:*(A::IdealSpoilingMatrix, B::BlochMcConnellMatrix{T,N}) where {T,N}

    C = BlochMcConnellMatrix{T}(N)
    mul!(C, A, B)
    return C

end

# /
Base.:/(M::Magnetization, a) = Magnetization(M.x / a, M.y / a, M.z / a)
Base.:/(M::MagnetizationMC{T,N}, a) where {T,N} = MagnetizationMC(ntuple(i -> M[i] / a, N)...)

# \
Base.:\(A::BlochMatrix, M::Magnetization) = Magnetization((Matrix(A) \ Vector(M))...)

function Base.:\(A::BlochMcConnellMatrix{T1,N}, M::MagnetizationMC{T2,N}) where {T1,T2,N}

    result = Matrix(A) \ Vector(M)
    return MagnetizationMC(ntuple(N) do i
        @inbounds r1 = result[3i-2]
        @inbounds r2 = result[3i-1]
        @inbounds r3 = result[3i]
        Magnetization(r1, r2, r3)
    end...)

end

"""
    add!(A, B)

Compute `A + B`, storing the result in `A`.
"""
function add!(M1::Magnetization, M2::Magnetization)

    M1.x += M2.x
    M1.y += M2.y
    M1.z += M2.z
    return nothing

end

function add!(M1::MagnetizationMC{T1,N}, M2::MagnetizationMC{T2,N}) where {T1,T2,N}

    for i = 1:N
        add!(M1[i], M2[i])
    end

end

"""
    add!(C, A, B)

Compute `A + B`, storing the result in `C`.
"""
function add!(C::AbstractVector, M1::Magnetization, M2::Magnetization)

    C[1] = M1.x + M2.x
    C[2] = M1.y + M2.y
    C[3] = M1.z + M2.z
    return nothing

end

function add!(C::AbstractVector, M1::MagnetizationMC{T1,N}, M2::MagnetizationMC{T2,N}) where {T1,T2,N}

    for i = 1:N
        add!(view(C, 3i-2:3i), M1[i], M2[i])
    end

end

function add!(C::AbstractMatrix, A::BlochMcConnellMatrix{T1,N}, B::BlochMcConnellMatrix{T2,N}) where {T1,T2,N}

    for j = 1:N, i = 1:N
        Ab = getblock(A, i, j)
        Bb = getblock(B, i, j)
        C[3i-2,3j-2] = Ab.a11 + Bb.a11
        C[3i-1,3j-2] = Ab.a21 + Bb.a21
        C[3i  ,3j-2] = Ab.a31 + Bb.a31
        C[3i-2,3j-1] = Ab.a12 + Bb.a12
        C[3i-1,3j-1] = Ab.a22 + Bb.a22
        C[3i  ,3j-1] = Ab.a32 + Bb.a32
        C[3i-2,3j]   = Ab.a13 + Bb.a13
        C[3i-1,3j]   = Ab.a23 + Bb.a23
        C[3i  ,3j]   = Ab.a33 + Bb.a33
    end

end

"""
    subtract!(A, B)

Compute `A - B`, storing the result in `A`.
"""
function subtract!(M1::Magnetization, M2::Magnetization)

    M1.x -= M2.x
    M1.y -= M2.y
    M1.z -= M2.z
    return nothing

end

function subtract!(M1::MagnetizationMC{T1,N}, M2::MagnetizationMC{T2,N}) where {T1,T2,N}

    for i = 1:N
        subtract!(M1[i], M2[i])
    end

end

"""
    subtract!(C, A, B)

Compute `A - B`, storing the result in `C`.
"""
function subtract!(C::AbstractMatrix{T}, ::UniformScaling, B::BlochMatrix) where {T}

    C[1,1] = one(T) - B.a11
    C[2,1] = -B.a21
    C[3,1] = -B.a31
    C[1,2] = -B.a12
    C[2,2] = one(T) - B.a22
    C[3,2] = -B.a32
    C[1,3] = -B.a13
    C[2,3] = -B.a23
    C[3,3] = one(T) - B.a33
    return nothing

end

function subtract!(C::AbstractMatrix, A::BlochMcConnellMatrix{T1,N}, B::BlochMcConnellMatrix{T2,N}) where {T1,T2,N}

    for j = 1:N, i = 1:N
        Ab = getblock(A, i, j)
        Bb = getblock(B, i, j)
        C[3i-2,3j-2] = Ab.a11 - Bb.a11
        C[3i-1,3j-2] = Ab.a21 - Bb.a21
        C[3i  ,3j-2] = Ab.a31 - Bb.a31
        C[3i-2,3j-1] = Ab.a12 - Bb.a12
        C[3i-1,3j-1] = Ab.a22 - Bb.a22
        C[3i  ,3j-1] = Ab.a32 - Bb.a32
        C[3i-2,3j]   = Ab.a13 - Bb.a13
        C[3i-1,3j]   = Ab.a23 - Bb.a23
        C[3i  ,3j]   = Ab.a33 - Bb.a33
    end

end

function subtract!(C::AbstractMatrix{T1}, A::UniformScaling, B::BlochMcConnellMatrix{T2,N}) where {T1,T2,N}

    for j = 1:N, i = 1:N
        b = getblock(B, i, j)
        if i == j
            C[3i-2,3j-2] = one(T1) - b.a11
            C[3i-1,3j-2] = -b.a21
            C[3i  ,3j-2] = -b.a31
            C[3i-2,3j-1] = -b.a12
            C[3i-1,3j-1] = one(T1) - b.a22
            C[3i  ,3j-1] = -b.a32
            C[3i-2,3j]   = -b.a13
            C[3i-1,3j]   = -b.a23
            C[3i  ,3j]   = one(T1) - b.a33
        else
            C[3i-2,3j-2] = -b.a11
            C[3i-1,3j-2] = -b.a21
            C[3i  ,3j-2] = -b.a31
            C[3i-2,3j-1] = -b.a12
            C[3i-1,3j-1] = -b.a22
            C[3i  ,3j-1] = -b.a32
            C[3i-2,3j]   = -b.a13
            C[3i-1,3j]   = -b.a23
            C[3i  ,3j]   = -b.a33
        end
    end

end

"""
    mul!(A, a)

Compute `A * a`, storing the result in `A`.
"""
LinearAlgebra.mul!(M::Magnetization, a) = (M.x *= a; M.y *= a; M.z *= a; nothing)
LinearAlgebra.mul!(M::MagnetizationMC{T,N}, a) where {T,N} = foreach(M -> mul!(M, a), M)

function LinearAlgebra.mul!(A::BlochDynamicsMatrix, t::Real)

    A.R1 *= t
    A.R2 *= t
    A.Δω *= t
    return nothing

end

function LinearAlgebra.mul!(E::ExchangeDynamicsMatrix, t::Real)

    E.r *= t
    return nothing

end

function LinearAlgebra.mul!(A::BlochMcConnellDynamicsMatrix, t::Real)

    for A in A.A
        mul!(A, t)
    end
    for E in A.E
        mul!(E, t)
    end

end

"""
    mul!(C, A, B)

Compute `A * B`, storing the result in `C`.
"""
function LinearAlgebra.mul!(M2::Magnetization, A::BlochMatrix, M1::Magnetization)

    M2.x = A.a11 * M1.x + A.a12 * M1.y + A.a13 * M1.z
    M2.y = A.a21 * M1.x + A.a22 * M1.y + A.a23 * M1.z
    M2.z = A.a31 * M1.x + A.a32 * M1.y + A.a33 * M1.z
    return nothing

end

function LinearAlgebra.mul!(M2::Magnetization, A::FreePrecessionMatrix, M1::Magnetization)

    M2.x = A.E2cosθ * M1.x + A.E2sinθ * M1.y
    M2.y = A.E2cosθ * M1.y - A.E2sinθ * M1.x
    M2.z = A.E1 * M1.z
    return nothing

end

function LinearAlgebra.mul!(M2::Magnetization, A::ExcitationMatrix, M1::Magnetization)

    mul!(M2, A.A, M1)

end

function LinearAlgebra.mul!(M2::Magnetization{T}, ::IdealSpoilingMatrix, M1::Magnetization) where {T}

    M2.x = zero(T)
    M2.y = zero(T)
    M2.z = M1.z
    return nothing

end

function LinearAlgebra.mul!(M2::Vector, A::BlochMatrix, M1::Magnetization)

    M2[1] = A.a11 * M1.x + A.a12 * M1.y + A.a13 * M1.z
    M2[2] = A.a21 * M1.x + A.a22 * M1.y + A.a23 * M1.z
    M2[3] = A.a31 * M1.x + A.a32 * M1.y + A.a33 * M1.z
    return nothing

end

function LinearAlgebra.mul!(M2::Vector, A::ExcitationMatrix, M1::Magnetization)

    mul!(M2, A.A, M1)

end

function LinearAlgebra.mul!(M2::MagnetizationMC{T1,N}, A::BlochMcConnellMatrix{T2,N}, M1::MagnetizationMC{T3,N}) where {T1,T2,T3,N}

    for i = 1:N
        M = M2[i]
        mul!(M, getblock(A, i, 1), M1[1])
        for j = 2:N
            muladd!(M, getblock(A, i, j), M1[j])
        end
    end

end

function LinearAlgebra.mul!(M2::MagnetizationMC{T1,N}, A::ExcitationMatrix, M1::MagnetizationMC{T2,N}) where {T1,T2,N}

    for i = 1:N
        mul!(M2[i], A, M1[i])
    end

end

function LinearAlgebra.mul!(M2::MagnetizationMC{T1,N}, ::IdealSpoilingMatrix, M1::MagnetizationMC{T2,N}) where {T1,T2,N}

    for i = 1:N
        mul!(M2[i], idealspoiling, M1[i])
    end

end

function LinearAlgebra.mul!(C::BlochMatrix, A::BlochMatrix, t::Real)

    C.a11 = A.a11 * t
    C.a21 = A.a21 * t
    C.a31 = A.a31 * t
    C.a12 = A.a12 * t
    C.a22 = A.a22 * t
    C.a32 = A.a32 * t
    C.a13 = A.a13 * t
    C.a23 = A.a23 * t
    C.a33 = A.a33 * t
    return nothing

end

function LinearAlgebra.mul!(C::BlochMatrix, A::BlochMatrix, B::BlochMatrix)

    C.a11 = A.a11 * B.a11 + A.a12 * B.a21 + A.a13 * B.a31
    C.a21 = A.a21 * B.a11 + A.a22 * B.a21 + A.a23 * B.a31
    C.a31 = A.a31 * B.a11 + A.a32 * B.a21 + A.a33 * B.a31
    C.a12 = A.a11 * B.a12 + A.a12 * B.a22 + A.a13 * B.a32
    C.a22 = A.a21 * B.a12 + A.a22 * B.a22 + A.a23 * B.a32
    C.a32 = A.a31 * B.a12 + A.a32 * B.a22 + A.a33 * B.a32
    C.a13 = A.a11 * B.a13 + A.a12 * B.a23 + A.a13 * B.a33
    C.a23 = A.a21 * B.a13 + A.a22 * B.a23 + A.a23 * B.a33
    C.a33 = A.a31 * B.a13 + A.a32 * B.a23 + A.a33 * B.a33
    return nothing

end

function LinearAlgebra.mul!(C::BlochMatrix, A::BlochMatrix, B::FreePrecessionMatrix)

    C.a11 = A.a11 * B.E2cosθ - A.a12 * B.E2sinθ
    C.a21 = A.a21 * B.E2cosθ - A.a22 * B.E2sinθ
    C.a31 = A.a31 * B.E2cosθ - A.a32 * B.E2sinθ
    C.a12 = A.a11 * B.E2sinθ + A.a12 * B.E2cosθ
    C.a22 = A.a21 * B.E2sinθ + A.a22 * B.E2cosθ
    C.a32 = A.a31 * B.E2sinθ + A.a32 * B.E2cosθ
    C.a13 = A.a13 * B.E1
    C.a23 = A.a23 * B.E1
    C.a33 = A.a33 * B.E1
    return nothing

end

function LinearAlgebra.mul!(C::BlochMatrix, A::FreePrecessionMatrix, B::BlochMatrix)

    C.a11 = A.E2cosθ * B.a11 + A.E2sinθ * B.a21
    C.a21 = A.E2cosθ * B.a21 - A.E2sinθ * B.a11
    C.a31 = A.E1 * B.a31
    C.a12 = A.E2cosθ * B.a12 + A.E2sinθ * B.a22
    C.a22 = A.E2cosθ * B.a22 - A.E2sinθ * B.a12
    C.a32 = A.E1 * B.a32
    C.a13 = A.E2cosθ * B.a13 + A.E2sinθ * B.a23
    C.a23 = A.E2cosθ * B.a23 - A.E2sinθ * B.a13
    C.a33 = A.E1 * B.a33
    return nothing

end

function LinearAlgebra.mul!(C::BlochMatrix, A::BlochDynamicsMatrix, B::BlochMatrix)

    C.a11 = A.R2 * B.a11 + A.Δω * B.a21
    C.a21 = A.R2 * B.a21 - A.Δω * B.a11
    C.a31 = A.R1 * B.a31
    C.a12 = A.R2 * B.a12 + A.Δω * B.a22
    C.a22 = A.R2 * B.a22 - A.Δω * B.a12
    C.a32 = A.R1 * B.a32
    C.a13 = A.R2 * B.a13 + A.Δω * B.a23
    C.a23 = A.R2 * B.a23 - A.Δω * B.a13
    C.a33 = A.R1 * B.a33
    return nothing

end

function LinearAlgebra.mul!(C::BlochMatrix, A::ExchangeDynamicsMatrix, B::BlochMatrix)

    C.a11 = A.r * B.a11
    C.a21 = A.r * B.a21
    C.a31 = A.r * B.a31
    C.a12 = A.r * B.a12
    C.a22 = A.r * B.a22
    C.a32 = A.r * B.a32
    C.a13 = A.r * B.a13
    C.a23 = A.r * B.a23
    C.a33 = A.r * B.a33
    return nothing

end

function LinearAlgebra.mul!(C::BlochMatrix, A::ExcitationMatrix, B::BlochMatrix)

    mul!(C, A.A, B)

end

function LinearAlgebra.mul!(C::BlochMatrix{T}, ::IdealSpoilingMatrix, B::BlochMatrix) where {T}

    C.a11 = zero(T)
    C.a21 = zero(T)
    C.a31 = B.a31
    C.a12 = zero(T)
    C.a22 = zero(T)
    C.a32 = B.a32
    C.a13 = zero(T)
    C.a23 = zero(T)
    C.a33 = B.a33
    return nothing

end

function LinearAlgebra.mul!(C::BlochMatrix{T}, A::FreePrecessionMatrix, B::FreePrecessionMatrix) where {T}

    C.a11 = A.E2cosθ * B.E2cosθ - A.E2sinθ * B.E2sinθ
    C.a12 = A.E2cosθ * B.E2sinθ + A.E2sinθ * B.E2cosθ
    C.a13 = zero(T)
    C.a21 = -C.a12
    C.a22 = C.a11
    C.a23 = zero(T)
    C.a31 = zero(T)
    C.a32 = zero(T)
    C.a33 = A.E1 * B.E1
    return nothing

end

function LinearAlgebra.mul!(C::BlochMatrix, A::FreePrecessionMatrix, B::ExcitationMatrix)

    mul!(C, A, B.A)

end

function LinearAlgebra.mul!(C::BlochMatrix, A::ExcitationMatrix, B::FreePrecessionMatrix)

    mul!(C, A.A, B)

end

function LinearAlgebra.mul!(C::BlochMatrix{T}, ::IdealSpoilingMatrix, B::FreePrecessionMatrix) where {T}

    C.a11 = zero(T)
    C.a21 = zero(T)
    C.a31 = zero(T)
    C.a12 = zero(T)
    C.a22 = zero(T)
    C.a32 = zero(T)
    C.a13 = zero(T)
    C.a23 = zero(T)
    C.a33 = B.E1
    return nothing

end

function LinearAlgebra.mul!(C::BlochMatrix{T}, A::BlochDynamicsMatrix, B::BlochDynamicsMatrix) where {T}

    C.a11 = A.R2 * B.R2 - A.Δω * B.Δω
    C.a21 = -A.Δω * B.R2 - A.R2 * B.Δω
    C.a31 = zero(T)
    C.a12 = -C.a21
    C.a22 = C.a11
    C.a32 = zero(T)
    C.a13 = zero(T)
    C.a23 = zero(T)
    C.a33 = A.R1 * B.R1
    return nothing

end

function LinearAlgebra.mul!(C::BlochMatrix{T}, A::BlochDynamicsMatrix, B::ExchangeDynamicsMatrix) where {T}

    C.a11 = A.R2 * B.r
    C.a21 = -A.Δω * B.r
    C.a31 = zero(T)
    C.a12 = A.Δω * B.r
    C.a22 = A.R2 * B.r
    C.a32 = zero(T)
    C.a13 = zero(T)
    C.a23 = zero(T)
    C.a33 = A.R1 * B.r
    return nothing

end

function LinearAlgebra.mul!(C::BlochMatrix, A::ExchangeDynamicsMatrix, B::BlochDynamicsMatrix)

    mul!(C, B, A)

end

function LinearAlgebra.mul!(C::BlochMatrix{T}, A::ExchangeDynamicsMatrix, B::ExchangeDynamicsMatrix) where {T}

    C.a11 = A.r * B.r
    C.a21 = zero(T)
    C.a31 = zero(T)
    C.a12 = zero(T)
    C.a22 = A.r * B.r
    C.a32 = zero(T)
    C.a13 = zero(T)
    C.a23 = zero(T)
    C.a33 = A.r * B.r
    return nothing

end

function LinearAlgebra.mul!(C::BlochMcConnellMatrix{T1,N}, A::BlochMcConnellMatrix{T2,N}, t::Real) where {T1,T2,N}

    for i = 1:N, j = 1:N
        mul!(getblock(C, i, j), getblock(A, i, j), t)
    end

end

function LinearAlgebra.mul!(
    C::BlochMcConnellMatrix{T1,N},
    A::ExcitationMatrix,
    B::BlochMcConnellMatrix{T2,N}
) where {T1,T2,N}

    for j = 1:N, i = 1:N
        mul!(getblock(C, i, j), A.A, getblock(B, i, j))
    end

end

function LinearAlgebra.mul!(
    C::BlochMcConnellMatrix{T1,N},
    ::IdealSpoilingMatrix,
    B::BlochMcConnellMatrix{T2,N}
) where {T1,T2,N}

    for j = 1:N, i = 1:N
        mul!(getblock(C, i, j), idealspoiling, getblock(B, i, j))
    end

end

function LinearAlgebra.mul!(
    C::AbstractBlochMcConnellMatrix{T1,N},
    A::AbstractBlochMcConnellMatrix{T2,N},
    B::AbstractBlochMcConnellMatrix{T3,N}
) where {T1,T2,T3,N}

    for j = 1:N, i = 1:N
        mul!(getblock(C, i, j), getblock(A, i, 1), getblock(B, 1, j))
        for k = 2:N
            muladd!(getblock(C, i, j), getblock(A, i, k), getblock(B, k, j))
        end
    end

end

"""
    div!(A, a)

Compute `A / a`, storing the result in `A`.
"""
div!(M::Magnetization, a) = (M.x /= a; M.y /= a; M.z /= a; nothing)
div!(M::MagnetizationMC{T,N}, a) where {T,N} = foreach(M -> div!(M, a), M)

"""
    muladd!(C, A, B)

Compute `A * B + C`, storing the result in `C`.
"""
function muladd!(M2::Magnetization, A::BlochMatrix, M1::Magnetization)

    M2.x += A.a11 * M1.x + A.a12 * M1.y + A.a13 * M1.z
    M2.y += A.a21 * M1.x + A.a22 * M1.y + A.a23 * M1.z
    M2.z += A.a31 * M1.x + A.a32 * M1.y + A.a33 * M1.z
    return nothing

end

function muladd!(M2::Magnetization, A::FreePrecessionMatrix, M1::Magnetization)

    M2.x += A.E2cosθ * M1.x + A.E2sinθ * M1.y
    M2.y += A.E2cosθ * M1.y - A.E2sinθ * M1.x
    M2.z += A.E1 * M1.z
    return nothing

end

function muladd!(M2::MagnetizationMC{T1,N}, A::BlochMcConnellMatrix{T2,N}, M1::MagnetizationMC{T3,N}) where {T1,T2,T3,N}

    for i = 1:N
        M = M2[i]
        for j = 1:N
            muladd!(M, getblock(A, i, j), M1[j])
        end
    end

end

function muladd!(C::BlochMatrix, A::BlochMatrix, t::Real)

    C.a11 += A.a11 * t
    C.a21 += A.a21 * t
    C.a31 += A.a31 * t
    C.a12 += A.a12 * t
    C.a22 += A.a22 * t
    C.a32 += A.a32 * t
    C.a13 += A.a13 * t
    C.a23 += A.a23 * t
    C.a33 += A.a33 * t
    return nothing

end

function muladd!(C::BlochMatrix, ::UniformScaling, t::Real)

    C.a11 += t
    C.a22 += t
    C.a33 += t
    return nothing

end

function muladd!(C::BlochMatrix, A::BlochMatrix, B::BlochMatrix)

    C.a11 += A.a11 * B.a11 + A.a12 * B.a21 + A.a13 * B.a31
    C.a21 += A.a21 * B.a11 + A.a22 * B.a21 + A.a23 * B.a31
    C.a31 += A.a31 * B.a11 + A.a32 * B.a21 + A.a33 * B.a31
    C.a12 += A.a11 * B.a12 + A.a12 * B.a22 + A.a13 * B.a32
    C.a22 += A.a21 * B.a12 + A.a22 * B.a22 + A.a23 * B.a32
    C.a32 += A.a31 * B.a12 + A.a32 * B.a22 + A.a33 * B.a32
    C.a13 += A.a11 * B.a13 + A.a12 * B.a23 + A.a13 * B.a33
    C.a23 += A.a21 * B.a13 + A.a22 * B.a23 + A.a23 * B.a33
    C.a33 += A.a31 * B.a13 + A.a32 * B.a23 + A.a33 * B.a33
    return nothing

end

function muladd!(C::BlochMatrix, A::BlochDynamicsMatrix, B::BlochMatrix)

    C.a11 += A.R2 * B.a11 + A.Δω * B.a21
    C.a21 += A.R2 * B.a21 - A.Δω * B.a11
    C.a31 += A.R1 * B.a31
    C.a12 += A.R2 * B.a12 + A.Δω * B.a22
    C.a22 += A.R2 * B.a22 - A.Δω * B.a12
    C.a32 += A.R1 * B.a32
    C.a13 += A.R2 * B.a13 + A.Δω * B.a23
    C.a23 += A.R2 * B.a23 - A.Δω * B.a13
    C.a33 += A.R1 * B.a33
    return nothing

end

function muladd!(C::BlochMatrix, A::ExchangeDynamicsMatrix, B::BlochMatrix)

    C.a11 += A.r * B.a11
    C.a21 += A.r * B.a21
    C.a31 += A.r * B.a31
    C.a12 += A.r * B.a12
    C.a22 += A.r * B.a22
    C.a32 += A.r * B.a32
    C.a13 += A.r * B.a13
    C.a23 += A.r * B.a23
    C.a33 += A.r * B.a33
    return nothing

end

function muladd!(C::BlochMatrix, A::BlochDynamicsMatrix, B::BlochDynamicsMatrix)

    C.a11 += A.R2 * B.R2 - A.Δω * B.Δω
    C.a21 -= A.Δω * B.R2 + A.R2 * B.Δω
    C.a12 += A.R2 * B.Δω + A.Δω * B.R2
    C.a22 += A.R2 * B.R2 - A.Δω * B.Δω
    C.a33 += A.R1 * B.R1
    return nothing

end

function muladd!(C::BlochMatrix, A::BlochDynamicsMatrix, B::ExchangeDynamicsMatrix)

    C.a11 += A.R2 * B.r
    C.a21 -= A.Δω * B.r
    C.a12 += A.Δω * B.r
    C.a22 += A.R2 * B.r
    C.a33 += A.R1 * B.r
    return nothing

end

function muladd!(C::BlochMatrix, A::ExchangeDynamicsMatrix, B::BlochDynamicsMatrix)

    muladd!(C, B, A)

end

function muladd!(C::BlochMatrix, A::ExchangeDynamicsMatrix, B::ExchangeDynamicsMatrix)

    C.a11 += A.r * B.r
    C.a22 += A.r * B.r
    C.a33 += A.r * B.r
    return nothing

end

function muladd!(C::BlochMcConnellMatrix{T1,N}, A::BlochMcConnellMatrix{T2,N}, t::Real) where {T1,T2,N}

    for i = 1:N, j = 1:N
        muladd!(getblock(C, i, j), getblock(A, i, j), t)
    end

end

function muladd!(C::BlochMcConnellMatrix{T,N}, I::UniformScaling, t::Real) where {T,N}

    for i = 1:N
        muladd!(getblock(C, i, i), I, t)
    end

end

"""
    subtractmul!(C, X, A, B)

Compute `(X - A) * B`, storing the result in `C`.
"""
function subtractmul!(
    M2::Magnetization,
    ::Nothing,
    A::BlochMatrix,
    M1::Magnetization
)

    M2.x = -A.a11 * M1.x - A.a12 * M1.y - A.a13 * M1.z
    M2.y = -A.a21 * M1.x - A.a22 * M1.y - A.a23 * M1.z
    M2.z = -A.a31 * M1.x - A.a32 * M1.y - A.a33 * M1.z
    return nothing

end

function subtractmul!(
    M2::Magnetization,
    ::UniformScaling,
    A::BlochMatrix,
    M1::Magnetization
)

    M2.x =  (1 - A.a11) * M1.x - A.a12 * M1.y - A.a13 * M1.z
    M2.y = -A.a21 * M1.x + (1 - A.a22) * M1.y - A.a23 * M1.z
    M2.z = -A.a31 * M1.x - A.a32 * M1.y + (1 - A.a33) * M1.z
    return nothing

end

function subtractmul!(
    M2::MagnetizationMC{T1,N},
    I::UniformScaling,
    A::BlochMcConnellMatrix{T2,N},
    M1::MagnetizationMC{T3,N}
) where {T1,T2,T3,N}

    for i = 1:N
        M = M2[i]
        if i == 1
            subtractmul!(M, I, getblock(A, i, 1), M1[1])
        else
            subtractmul!(M, nothing, getblock(A, i, 1), M1[1])
        end
        for j = 2:N
            if j == i
                subtractmuladd!(M, I, getblock(A, i, j), M1[j])
            else
                subtractmuladd!(M, nothing, getblock(A, i, j), M1[j])
            end
        end
    end

end

"""
    subtractmuladd!(C, X, A, B)

Compute `(X - A) * B + C`, storing the result in `C`.
"""
function subtractmuladd!(
    M2::Magnetization,
    ::Nothing,
    A::BlochMatrix,
    M1::Magnetization
)

    M2.x -= A.a11 * M1.x + A.a12 * M1.y + A.a13 * M1.z
    M2.y -= A.a21 * M1.x + A.a22 * M1.y + A.a23 * M1.z
    M2.z -= A.a31 * M1.x + A.a32 * M1.y + A.a33 * M1.z
    return nothing

end

function subtractmuladd!(
    M2::Magnetization,
    ::UniformScaling,
    A::BlochMatrix,
    M1::Magnetization
)

    M2.x +=  (1 - A.a11) * M1.x - A.a12 * M1.y - A.a13 * M1.z
    M2.y += -A.a21 * M1.x + (1 - A.a22) * M1.y - A.a23 * M1.z
    M2.z += -A.a31 * M1.x - A.a32 * M1.y + (1 - A.a33) * M1.z
    return nothing

end

"""
    absolutesum(A)

Compute the sum of the absolute values of the elements of `A`.
"""
function absolutesum(A::BlochDynamicsMatrix)

    return 2 * abs(A.R2) + 2 * abs(A.Δω) + abs(A.R1)

end

function absolutesum(E::ExchangeDynamicsMatrix)

    return 3 * abs(E.r)

end

function absolutesum(A::BlochMcConnellDynamicsMatrix{T,N,M}) where {T,N,M}

    result = absolutesum(A.A[1])
    for i = 2:N
        result += absolutesum(A.A[i])
    end
    for E in A.E
        result += absolutesum(E)
    end
    return result

end


text:
using Colors

include("helpers/page.jl")
include("helpers/doc.jl")

intro = md"""
$(title(2, "Signal"))
$(vskip(1em))

# Functions
$(vskip(1em))

$(
   showdocs([subscribe, intent, sampler,
             watch!, trigger!, bubble])
)

$(vskip(4em))

"""

function main(window)
    docpage(intro)
end


text:
@with_kw struct IF2Parameter <: AbstractIFParameter
    Ee::SNNFloat = 0mV
    Ei::SNNFloat = 0mV
end

@with_kw mutable struct IF2 <: AbstractIF
    param::IF2Parameter = IF2Parameter()
    N::SNNInt = 100
end

function integrate!(p::IF2, param::IF2Parameter, dt::SNNFloat)
    @unpack N = p
    @unpack Ee, Ei = param
    @inbounds for i = 1:N
        v[i] += dt * (ge[i] * (Ee - v[i]) + gi[i] * (Ei - v[i]) - (v[i] - El)) / τm
        ge[i] += dt * -ge[i] / τe
        gi[i] += dt * -gi[i] / τi
    end
    @inbounds for i = 1:N
        fire[i] = v[i] > Vt
        v[i] = ifelse(fire[i], Vr, v[i])
    end
end


text:
"""
Canonical polyadic N-mode tensor decomposition (CANDECOMP/PARAFAC).
"""
struct CANDECOMP{T<:Number, N} <: TensorDecomposition{T, N}
  factors::NTuple{N, Matrix{T}} # rank X i-th tensor mode factor matrices
  lambdas::Vector{T}
  props::Dict{Symbol, Any}      # extra properties

  CANDECOMP{T, N}(factors::NTuple{N, Matrix{T}}, lambdas::Vector{T}) where {T<:Number, N} =
    new(factors, lambdas, Dict{Symbol, Any}())
end

CANDECOMP(factors::NTuple{N, Matrix{T}}, lambdas::Vector{T}) where {T<:Number, N} =
    CANDECOMP{T, N}(factors, lambdas)

"""
Returns CANDECOMP rank.
"""
LinearAlgebra.rank(decomp::CANDECOMP) = length(decomp.lambdas)

"""
Re-composes the tensor from CANDECOMP decomposition.
"""
@generated function compose!(dest::Array{T,N}, factors::NTuple{N, Matrix{T}}, lambdas::Vector{T}) where {T <: Number, N}
  quote
    csize = @ntuple($N, k -> size(factors[k], 2))
    @nall($N, k -> csize[k] == length(lambdas)) ||
        throw(DimensionMismatch("length of lambdas ($(length(lambdas))) doesn't match the factors columns ($csize)"))
    rsize = @ntuple($N, k -> size(factors[k], 1))
    size(dest) == rsize ||
        throw(DimensionMismatch("dest dimensions $(size(dest)) do not match the factors rows ($rsize)"))
    @nloops $N i dest begin
        elm = zero(T)
        # FIXME traversal of factors is not very efficient, would be better to have them transposed
        @inbounds for j in eachindex(lambdas)
            elm += lambdas[j] * (*(@ntuple($N, k -> factors[k][i_k, j])...))
        end
        @inbounds @nref($N, dest, i) = elm
    end
    return dest
  end
end

compose(factors::NTuple{N, Matrix{T}}, lambdas::Vector{T}) where {T <: Number, N} =
    compose!(Array{T, N}(undef, ntuple(i -> size(factors[i], 1), N)), factors, lambdas)

compose(decomp::CANDECOMP) = compose(decomp.factors, decomp.lambdas)

compose!(dest::Array{T,N}, decomp::CANDECOMP{T,N}) where {T <: Number, N} =
    compose!(dest, decomp.factors, decomp.lambdas)

"""
Calculates canonical polyadic tensor decomposition (CANDECOMP/PARAFAC).

Returns:
  `CANDECOMP` object
"""
function candecomp(tnsr::AbstractArray{T,N},
                   r::Integer,
                   initial_guess::NTuple{N, Matrix{T}};
                   method::Symbol=:ALS,
                   tol::Float64=1e-5,
                   maxiter::Integer=100,
                   compute_error::Bool=false,
                   verbose::Bool=true) where {T <: Number, N}

    _check_tensor(tnsr, r)
    verbose && @info("initializing factor matrices...")
    all([(size(tnsr, i), r) == size(initial_guess[i]) for i in 1:N]) || throw(ArgumentError("dimension of initial guess does not match input tensor."))
    factors = collect(Matrix{T}, initial_guess)
    verbose && @info("applying candecomp $method method...")
    res = _candecomp(Val{method}, tnsr, r, factors, tol, maxiter, verbose)
    compute_error && _set_rel_residue(res, tnsr)
    return res
end

"Stub for non-implemented CANDECOMP algorithms."
_candecomp(
  method::Val{S},
  tnsr::AbstractArray{T,N},
  r::Integer,
  factors::Vector{Matrix{T}},
  tol::Float64,
  maxiter::Integer,
  verbose::Bool) where {T,N,S<:Symbol} =
    throw(ArgumentError("Unknown CANDECOMP method: $method"))

"""
Computes CANDECOMP by ALS (Alternating Least Squares) method.
"""
function _candecomp(
    method::Type{Val{:ALS}},
    tnsr::AbstractArray{T,N},
    r::Integer,
    factors::Vector{Matrix{T}},
    tol::Float64,
    maxiter::Integer,
    verbose::Bool) where {T <: Number, N}

    gram = [F'F for F in factors]
    tnsr_norm = norm(tnsr)
    tnsr_flat = _col_unfold(tnsr, N)
    tnsr_size = size(tnsr)
    niters = 0
    converged = false
    resid = tnsr_norm
    lbds = Matrix{T}(undef, 1, r)
    V = Matrix{T}(undef, length(tnsr) ÷ minimum(tnsr_size), r)
    pb = Progress(maxiter, "ALS iterations ")
    while !converged && niters < maxiter
        update!(pb, niters)
        nVi = 0
        @inbounds for i in 1:N
            idx = [N:-1:i + 1; i - 1:-1:1]
            nVi = length(tnsr) ÷ tnsr_size[i]
            Vi = view(V, 1:nVi, :)
            khatrirao!(Vi, factors[idx])
            factors[i] .= _row_unfold(tnsr, i) * Vi / reduce((x, y) -> x .* y, gram[idx])
            sum!(abs, lbds, factors[i])
            factors[i] ./= lbds
            mul!(gram[i], factors[i]', factors[i])
        end
        resid_old = resid
        resid = norm(view(V, 1:nVi, :) * permutedims((factors[N] .* lbds)) .- tnsr_flat)
        converged = abs(resid - resid_old) < tol * resid_old
        niters += 1
    end
    finish!(pb)

    verbose && _iter_status(converged, niters, maxiter)

    return CANDECOMP(tuple(factors...), dropdims(lbds, dims=1))
end

"""
Computes CANDECOMP by SGSD (Simultaneous Generalized Schur Decomposition) method.
"""
function _candecomp(
    method::Type{Val{:SGSD}},
    tnsr::AbstractArray{T, N},
    r::Int,
    factors::Vector{<:StridedMatrix},
    tol::Float64,
    maxiter::Integer,
    verbose::Bool) where {T <: Number, N}

    ndims(tnsr) == 3 || throw(ArgumentError("This algorithm only applies to 3-mode tensors."))
    length(factors) == 3 || throw(ArgumentError("3 factor matrices expected."))
    (n1, n2, n3) = size(tnsr)
    IB1 = min(n1 - 1, r)
    IB2 = (n2 == r) ? 2 : 1
    Q = qr(factors[1]).Q * Matrix(I, n1, n1)
    Z = qr(reverse(factors[2], dims=2)).Q * Matrix(I, n2, n2)
    q = Matrix{T}(undef, n1, n1)
    z = Matrix{T}(undef, n2, n2)

    @tensor R[4,5,3] := tnsr[1,2,3] * Q[1,4] * Z[2,5]
    R2 = similar(R)

    res = norm(tnsr)
    converged = false
    niters = 0
    pb = Progress(maxiter, "SGSD iterations ")
    @inbounds while !converged && niters < maxiter
        update!(pb, niters)
        # reset to diagonal
        @inbounds for i in 1:n1, j in 1:n1
            q[i, j] = ifelse(i == j, 1.0, 0.0)
        end
        @inbounds for i in 1:n2, j in 1:n2
            z[i, j] = ifelse(i == j, 1.0, 0.0)
        end

        for i in 1:IB1
            q[:, i:n1] *= svd(view(q, :, i:n1)' * view(R, :, n2 - r + i, :)).U
        end

        @tensor R2[4,2,3] = R[1,2,3] * q[1,4]
        for i in r:-1:IB2
            z[:, 1:n2 - r + i] *= reverse(svd(view(R2, i, :, :)' * view(z, :, 1:n2 - r + i)).V, dims=2)
        end
        Q *= q
        Z *= z

        @tensor R[4,5,3] = tnsr[1,2,3] * Q[1,4] * Z[2,5]

        res_old = res
        res = norm(tril(dropdims(sum(abs2, R, dims=3), dims=3), n2 - r - 1))
        converged = abs(res - res_old) < tol * res_old
        niters += 1
    end
    finish!(pb)
    verbose && _iter_status(converged, niters, maxiter)

    R = R[1:r, n2-r+1:n2, :]
    M = Array{T,3}(undef, r, r, 2)
    @inbounds for k in 1:2, j in 1:r, i in 1:r
        M[i, j, k] = i == j ? 1.0 : 0.0
    end
    @inbounds for i in (r - 1):-1:1, j = (i + 1):r
        if i + 1 < j
            d = (i + 1):(j - 1)
            M[i, j, :] = hcat(view(R, j, j, :), view(R, i, i, :)) \
                (view(R, i, j, :) .- mapslices(R3 -> sum(view(M, i, d, 1) .* diag(R3)[d] .* view(M, d, j, 2)), R, dims=(1, 2)))
        end
    end

    factors[1] = view(Q, :, 1:r) * view(M, :, :, 1)
    factors[2] = view(Z, :, n2 - r + 1:n2) * view(M, :, :, 2)'
    factors[3] = _row_unfold(tnsr, 3) * khatrirao(factors[2], factors[1]) /
        ((factors[2]'factors[2]) .* (factors[1]'factors[1]))

    # normalize
    lbds = ones(1, r)
    for i in 1:3
        lbd = mapslices(norm, factors[i], dims=1)
        factors[i] ./= lbd
        lbds .*= lbd
    end

    return CANDECOMP(tuple(factors...), dropdims(lbds, dims=1))
end

"""
Codes of implemented CANDECOMP methods.
"""
const CANDECOMP_methods = Symbol[t.parameters[1] for t in filter(t -> isa(t, DataType),
                                                                 [Base.unwrap_unionall(m.sig).parameters[2].parameters[1] for m in methods(_candecomp)])] |> Set{Symbol}


text:
using BenchmarkTools, Test

N = 2^20
x = fill(1.0f0, N)  # a vector filled with 1.0 (Float32)
y = fill(2.0f0, N)  # a vector filled with 2.0

y .+= x             # increment each element of y with the corresponding element of x

@test all(y .== 3.0f0)

function sequential_add!(y, x)
    for i in eachindex(y, x)
        @inbounds y[i] += x[i]
    end
    return nothing
end

fill!(y, 2)
sequential_add!(y, x)
@test all(y .== 3.0f0)

function parallel_add!(y, x)
    Threads.@threads for i in eachindex(y, x)
        @inbounds y[i] += x[i]
    end
    return nothing
end

fill!(y, 2)
parallel_add!(y, x)
@test all(y .== 3.0f0)

@btime sequential_add!($y, $x)
@btime parallel_add!($y, $x)
@btime y .+= x             # increment each element of y with the corresponding element of x

using CuArrays

x_d = cufill(1.0f0, N)  # a vector stored on the GPU filled with 1.0 (Float32)
y_d = cufill(2.0f0, N)  # a vector stored on the GPU filled with 2.0

y_d .+= x_d
@test all(Array(y_d) .== 3.0f0)
@btime y_d .+= x_d

function add_broadcast!(y, x)
    CuArrays.@sync y .+= x
    return
end
@btime add_broadcast!(y_d, x_d)


text:
module HospitalizationData

using CSV
using Dates
using DataFrames

basepath = normpath(@__DIR__, "../../")

export hospitalizations

function hospitalizations(locations::Array, start_date::Date, end_date::Date;
		level::Symbol=:state, source::Symbol=:cdc,
		forecast_type::Symbol=:active, patient_type::Symbol=:all, bound_type::Symbol=:mean,
		hospitalized_days::Int=-1, active_from_admitted::Bool=false,
)
	if source == :cdc || (source == :gt && forecast_type == :active)
		@assert level == :state
		@assert forecast_type == :active
		@assert patient_type  == :all
		gt = load_cdc_gt()
		gt = cdc_gt_filter!(gt, locations, start_date, end_date)
		patients = cdc_gt_history(gt, locations, start_date, end_date, bound_type=bound_type)
	elseif source == :covidtracking || (source == :gt && forecast_type == :admitted)
		@assert level == :state
		@assert patient_type  == :all
		@assert bound_type    == :mean
		forecast = load_covidtracking()
		forecast = covidtracking_filter!(forecast, locations, start_date, end_date)
		patients = covidtracking_history(forecast, forecast_type=forecast_type)
	end

	return patients
end

########################
#### COVIDTracking #####
########################

function load_covidtracking()
	gt = DataFrame(CSV.File(joinpath(basepath, "data/hospitalizations/covidtracking/2020-06-30/states_daily.csv"), silencewarnings=true), copycols=true)
	gt.date = map(d -> Date(string(d), "yyyymmdd"), gt.date)
	return gt
end

function covidtracking_filter!(gt::DataFrame, locations::Array{String,1}, start_date::Date, end_date::Date)
	@assert locations == sort(locations)
	@assert Date(2020, 1, 22) <= start_date <= end_date <= Date(2020, 6, 6)

	filter!(row -> start_date <= row.date <= end_date, gt)
	filter!(row -> row.state in locations, gt)
	sort!(gt, [:state, :date])
	return gt
end

function covidtracking_history(gt::DataFrame; forecast_type::Symbol=:active, keepmissings::Bool=false)
	groups = groupby(gt, :state, sort=true)
	if forecast_type == :active
		patients = vcat([g.hospitalizedCurrently' for g in groups]...)
	elseif forecast_type == :admitted
		patients = vcat([g.hospitalizedIncrease' for g in groups]...)
	else
		error("Invalid forecast type for covidtracking")
	end
	if !keepmissings
		patients[ismissing.(patients)] .= 0
	end
	return Float32.(patients)
end

########################
######## CDC GT ########
########################

function load_cdc_gt()
	goodcols = Dict(
		:state => :state,
		:collectionDate  => :date,
		:ICUBeds_Occ_AnyPat_Est  => :icu_active,
		:ICUBeds_Occ_AnyPat_LoCI => :icu_active_lb,
		:ICUBeds_Occ_AnyPat_UpCI => :icu_active_ub,
		:InpatBeds_Occ_AnyPat_Est  => :all_active,
		:InpatBeds_Occ_AnyPat_LoCI => :all_active_lb,
		:InpatBeds_Occ_AnyPat_UpCI => :all_active_ub,
		:InpatBeds_Occ_COVID_Est  => :all_covid_active,
		:InpatBeds_Occ_COVID_LoCI => :all_covid_active_lb,
		:InpatBeds_Occ_COVID_UpCI => :all_covid_active_ub,
	)
	goodcolnames = sort(collect(keys(goodcols)),rev=true)

	fn = joinpath(basepath, "data/hospitalizations/cdc/2020-06-16.csv")
	gt = DataFrame(CSV.File(fn, skipto=3, dateformat="dduuuyyyy", select=goodcolnames))
	rename!(gt, goodcols)
	return gt
end

function cdc_gt_filter!(gt::DataFrame, locations::Array{String,1}, start_date::Date, end_date::Date)
	@assert locations == sort(locations)
	@assert Date(2020, 4, 1) <= start_date <= end_date <= Date(2020, 6, 2)

	filter!(row -> row.state in locations, gt)
	filter!(row -> start_date <= row.date <= end_date, gt)
	sort!(gt, [:state, :date])
	return gt
end

function cdc_gt_history(
		gt::DataFrame,
		locations::Array{String,1}, start_date::Date, end_date::Date;
		patient_type::Symbol=:all, bound_type::Symbol=:mean,
)
	colselect = Dict(
		(:all, :mean) => :all_covid_active,
		(:all, :lb)   => :all_covid_active_lb,
		(:all, :ub)   => :all_covid_active_ub,
	)
	col = colselect[(patient_type, bound_type)]

	gt_by_state = groupby(gt, :state, sort=true)
	gt_by_state_date = Dict(s.state[1] => Dict(zip(s.date, s[!,col])) for s in gt_by_state)
	patients = [haskey(gt_by_state_date[s], d) ? gt_by_state_date[s][d] : missing for s in locations, d in start_date:Day(1):end_date]
	@assert !any(ismissing.(patients))

	return Float32.(patients)
end

end;

max_stars_repo_name: flixpar/covid-resource-allocation

max_issues_repo_name: flixpar/covid-resource-allocation

max_forks_repo_name: flixpar/covid-resource-allocation


text:
#
# Functions related to Halo orbits
#


text:
@testset "Range Interface" begin
    for r in (OneToMRange(10),
              OneToSRange(10),
              UnitMRange(1, 10),
              UnitSRange(1, 10),
              StepMRange(1, 1, 10),
              StepSRange(1, 1, 10),
              LinMRange(1, 10, 10),
              LinSRange(1, 10, 10),
              StepMRangeLen(1, 1, 10),
              StepSRangeLen(1, 1, 10)
             )
        @testset "$(typeof(r))" begin
            @testset "first" begin
                @test @inferred(first(r)) == 1
            end
            @testset "firstindex" begin
                @test @inferred(firstindex(r)) == 1
            end
            @testset "step" begin
                @test @inferred(step(r)) == 1
            end
            @testset "last" begin
                @test @inferred(last(r)) == 10
            end
            @testset "lastindex" begin
                @test @inferred(lastindex(r)) == 10
            end
            @testset "length" begin
                @test @inferred(length(r)) == 10
            end

            if r isa StaticRanges.StaticUnitRange
                io = IOBuffer()
                show(io, r)
                str = String(take!(io))
                @test str == string(typeof(r).name, "(", first(r), ":", last(r), ")")
            end

            if r isa StaticRanges.AbstractStepRangeLen
                io = IOBuffer()
                show(io, r)
                str = String(take!(io))
                @test str == string(typeof(r).name, "(", first(r), ":", step(r), ":", last(r), ")")
            end

            if r isa StaticRanges.AbstractStepRange
                @testset "start-property" begin
                    @test @inferred(getproperty(r, :start)) == 1
                end
                @testset "step-property" begin
                    @test @inferred(getproperty(r, :step)) == 1
                end
                @testset "stop-property" begin
                    @test @inferred(getproperty(r, :stop)) == 10
                end
                io = IOBuffer()
                show(io, r)
                str = String(take!(io))
                @test str == string(typeof(r).name, "(", first(r), ":", step(r), ":", last(r), ")")
            end

            if r isa StaticRanges.AbstractLinRange
                @testset "start-property" begin
                    @test getproperty(r, :start) == 1
                end
                @testset "stop-property" begin
                    @test getproperty(r, :stop) == 10
                end
                @testset "lendiv-property" begin
                    @test getproperty(r, :lendiv) == 9
                    @test StaticRanges.lendiv(r) == 9
                end
                @testset "len-property" begin
                    @test getproperty(r, :len) == 10
                end
            end
        end
    end
end


text:
macro cache(expr)
    name = expr.args[2].args[1].args[1]
    fields = expr.args[3].args[2:2:end]
    cache_vars = Expr[]
    rand_vars = Expr[]
    jac_vars = Pair{Symbol,Expr}[]
    ratenoise_vars = Expr[]
    for x in fields
        if x.args[2] == :uType || x.args[2] == :rateType ||
         x.args[2] == :kType || x.args[2] == :uNoUnitsType # || x.args[2] == :possibleRateType
            push!(cache_vars, :(c.$(x.args[1])))
        elseif x.args[2] == :JCType
        push!(cache_vars, :(c.$(x.args[1]).duals...))
      elseif x.args[2] == :GCType
        push!(cache_vars, :(c.$(x.args[1]).duals))
      elseif x.args[2] == :DiffCacheType
        push!(cache_vars, :(c.$(x.args[1]).du))
        push!(cache_vars, :(c.$(x.args[1]).dual_du))
      elseif x.args[2] == :JType || x.args[2] == :WType
        push!(jac_vars, x.args[1] => :(c.$(x.args[1])))
      elseif x.args[2] == :randType
        push!(rand_vars, :(c.$(x.args[1])))
      elseif x.args[2] == :rateNoiseType || x.args[2] == :rateNoiseCollectionType
        # Should be a pair for handling non-diagonal
        push!(ratenoise_vars, :(c.$(x.args[1])))
        end
    end
    quote
        $expr
      $(esc(:full_cache))(c::$name) = tuple($(cache_vars...))
      $(esc(:jac_iter))($(esc(:c))::$name) = tuple($(jac_vars...))
      $(esc(:rand_cache))($(esc(:c))::$name) = tuple($(rand_vars...))
      $(esc(:ratenoise_cache))($(esc(:c))::$name) = tuple($(ratenoise_vars...))
    end
end

max_stars_repo_name: pgunnink/LLGTools

max_issues_repo_name: pgunnink/LLGTools

max_forks_repo_name: pgunnink/LLGTools


text:
"""
    simulate0()

Obtain the percentile ranks that would be expected for the set of molecules if none were bound together. This state is
simulated by randomly displacing all molecules within a given radius. This is done to simulate freedom of motion but
to roughly maintain local density. This radius should be the same as that used for montecarloaffinity().
"""
function simulate0(molecules1, molecules2, localradius, mc_iterations)
    displacements1 = randomcoordinates2d(length(molecules1), localradius)
    displacements2 = randomcoordinates2d(length(molecules2), localradius)

    newcoordinates1 = extractcoordinates(molecules1) .+ displacements1
    newcoordinates2 = extractcoordinates(molecules2) .+ displacements2

    newmolecules1 = [Molecule(Localization(i, "", newcoordinates1[1,i], newcoordinates1[2,i], newcoordinates1[3,i], 0, 1, 1)) for i ∈ 1:length(molecules1)]
    newmolecules2 = [Molecule(Localization(i, "", newcoordinates2[1,i], newcoordinates2[2,i], newcoordinates2[3,i], 0, 1, 1)) for i ∈ 1:length(molecules2)]

    neighbors1, neighbors2, distances = exclusivenearestneighbors(newmolecules1, newmolecules2)

    percentileranks = montecarloaffinity(newmolecules1, newmolecules2, neighbors1, neighbors2, distances, localradius, mc_iterations)
    percentileranks, distances
end

"""
    simulate100()

Obtain the percentile ranks that would be expected for the set of molecules of all were bound together. This state is
simulated by removing one random member of each pair and regenerating them within a given radius
of their original partner. This radius should be your best guess at how far apart fluorophores of a bound pair of your
molecule pairs should be, in nanometers. 80 nm would be a reasonable guess for average proteins with an antibody stack.
"""
function simulate100(molecules1, molecules2, neighbors1, neighbors2, bindingradius, localradius, iterations)
    replace1 = rand([true false], length(neighbors1))
    replace2 = .!replace1

    neighborindexes1_replace = [n.index for n ∈ neighbors1[replace1]]
    neighborindexes1_keep = [n.index for n ∈ neighbors1[replace2]]
    neighborindexes2_replace = [n.index for n ∈ neighbors2[replace2]]
    neighborindexes2_keep = [n.index for n ∈ neighbors2[replace1]]

    newmolecules1 = filter(m -> m.index ∉ neighborindexes1_replace, molecules1)
    newmolecules2 = filter(m -> m.index ∉ neighborindexes2_replace, molecules2)

    displacements1 = randomcoordinates2d(length(neighborindexes1_replace), bindingradius)
    displacements2 = randomcoordinates2d(length(neighborindexes2_replace), bindingradius)

    newcoordinates1 = extractcoordinates(filter(m -> m.index ∈ neighborindexes2_keep, newmolecules2)) .+ displacements1
    newcoordinates2 = extractcoordinates(filter(m -> m.index ∈ neighborindexes1_keep, newmolecules1)) .+ displacements2

    append!(newmolecules1, [Molecule(Localization(neighborindexes1_replace[i], "", newcoordinates1[1,i], newcoordinates1[2,i], newcoordinates1[3,i], 0, 1, 1)) for i ∈ 1:length(neighborindexes1_replace)])
    append!(newmolecules2, [Molecule(Localization(neighborindexes2_replace[i], "", newcoordinates2[1,i], newcoordinates2[2,i], newcoordinates2[3,i], 0, 1, 1)) for i ∈ 1:length(neighborindexes2_replace)])

    newneighbors1, newneighbors2, distances = exclusivenearestneighbors(newmolecules1, newmolecules2)

    percentileranks = montecarloaffinity(newmolecules1, newmolecules2, newneighbors1, newneighbors2, distances, localradius, iterations)
    percentileranks, distances
end


text:
using Test

@testset "Backscatter - Elemental" begin
    @test isapprox(η(LoveScott1978η, n"C", 10.0e3), 0.08432227, atol = 1.0e-5)
    @test isapprox(η(LoveScott1978η, n"Fe", 30.0e3), 0.2804824, atol = 1.0e-5)
    @test isapprox(η(LoveScott1978η, n"Ce", 10.0e3), 0.44044972, atol = 1.0e-5)
    @test isapprox(η(LoveScott1978η, n"U", 12.0e3), 0.5104889, atol = 1.0e-5)

    @test η(Tomlin1963, n"H", 10.0e3) > 0.0
    @test isapprox(η(Tomlin1963, n"C", 10.0e3), 0.04862657, atol = 1.0e-5)
    @test isapprox(η(Tomlin1963, n"Fe", 30.0e3), 0.2930160, atol = 1.0e-5)
    @test isapprox(η(Tomlin1963, n"Ce", 10.0e3), 0.42674050, atol = 1.0e-5)
    @test isapprox(η(Tomlin1963, n"U", 12.0e3), 0.5036314, atol = 1.0e-5)

    @test all(map(z -> η(August1989η, elements[z], 15.0e3) == η(elements[z], 15.0e3), 1:99))
end

@testset "Backscatter - Material" begin

    @test isapprox(
        η(LoveScott1978η, n"C", 10.0e3),
        η(LoveScott1978η, mat"0.9999*C+0.0001*H", 10.0e3),
        atol = 1.0e-5,
    )
    @test isapprox(
        η(LoveScott1978η, n"U", 30.0e3),
        η(LoveScott1978η, mat"0.99999*U+0.00001*H", 30.0e3),
        atol = 1.0e-5,
    )
    @test isapprox(
        η(LoveScott1978η, n"C", 10.0e3),
        η(LoveScott1978η, mat"0.9999*C+0.0001*H", 10.0e3),
        atol = 1.0e-5,
    )

    @test isapprox(
        NeXLCore.electronfraction(n"Au", mat"0.756*Au+0.244*Cu"),
        0.731,
        atol = 0.001,
    )
    @test isapprox(NeXLCore.electronfraction(n"Th", mat"ThSiO4"), 0.6618, atol = 0.0001)
    @test isapprox(NeXLCore.electronfraction(n"Si", mat"ThSiO4"), 0.1029, atol = 0.0001)

    @test isapprox(η(LoveScott1978η, mat"0.224*Au+0.775*Ag", 20.0e3), 0.4271, atol = 1.0e-3)
    @test isapprox(η(LoveScott1978η, mat"0.80*Au+0.199*Ag", 20.0e3), 0.4761, atol = 1.0e-3)

    @test isapprox(η(LoveScott1978η, mat"0.201*Au+0.798*Cu", 20.0e3), 0.3511, atol = 1.0e-3)
    @test isapprox(η(LoveScott1978η, mat"0.801*Au+0.198*Cu", 20.0e3), 0.4610, atol = 1.0e-3)
end


text:
using Test, ModelingToolkit

@parameters t

@connector function Foo(;name)
    @variables x(t)
    ODESystem(Equation[], t, [x], [], defaults=Dict(x=>1.0))
end

@connector function Goo(;name)
    @variables x(t)
    @parameters p
    ODESystem(Equation[], t, [x], [p], defaults=Dict(x=>1.0, p=>1.0))
end

ModelingToolkit.connect(::Type{<:Foo}, sys1, sys2) = [sys1.x ~ sys2.x]
@named f1 = Foo()
@named f2 = Foo()
@named g = Goo()

@test isequal(connect(f1, f2), [f1.x ~ f2.x])
@test_throws ArgumentError connect(f1, g)

# Note that since there're overloadings, these tests are not re-runable.
ModelingToolkit.promote_connect_rule(::Type{<:Foo}, ::Type{<:Goo}) = Foo
@test isequal(connect(f1, g), [f1.x ~ g.x])
ModelingToolkit.promote_connect_rule(::Type{<:Goo}, ::Type{<:Foo}) = Foo
@test isequal(connect(f1, g), [f1.x ~ g.x])
# test conflict
ModelingToolkit.promote_connect_rule(::Type{<:Goo}, ::Type{<:Foo}) = Goo
@test_throws ArgumentError connect(f1, g)


text:
function vecOuterProd!(C, a, b)
    @turbo for i ∈ eachindex(b)
        for j ∈ eachindex(a)
            C[j,i] = a[j]*b[i]
        end
    end
end

max_stars_repo_name: GrantHecht/OptimalEstimationProject


text:
"""
    FeedbackNets

Implements deep networks with feedback operations from higher to lower layers.
Uses Flux as a backend.
"""
module FeedbackNets
using Reexport

include("Splitters.jl")
include("AbstractMergers.jl")
include("Mergers.jl")
include("AbstractFeedbackNets.jl")
include("FeedbackChains.jl")
include("FeedbackTrees.jl")
include("modelfactory/ModelFactory.jl")
@reexport using .Splitters
@reexport using .AbstractMergers
@reexport using .Mergers
@reexport using .AbstractFeedbackNets
@reexport using .FeedbackChains
@reexport using .FeedbackTrees
@reexport using .ModelFactory

# add definition for dictionaries to Flux._truncate in order to make truncate!
# work states of Recur-ed FeedbackChains and FeedbackTrees
import Flux: _truncate
_truncate(d::Dict) = Dict(key => _truncate(val) for (key, val) in pairs(d))

end # module FeedbackNets


text:
struct DummySolver <: SemialgebraicSets.AbstractAlgebraicSolver end

@testset "Basic semialgebraic set" begin
    Mod.@polyvar x y
    @test isa(FullSpace(), FullSpace)
    V = @set x * y == 1
    @test V isa AlgebraicSet
    @test_throws ArgumentError addinequality!(V, x*y)
    @testset "Basic" begin
        S = @set x - y == 0 && x^2*y >= 1
        addequality!(S, x^2 - y)
        addinequality!(S, x + y - 1)
        # Algebraic set forces `Float64`
        @test S isa BasicSemialgebraicSet{Float64}
        @test sprint(show, S) == "{ (x, y) | x - y = 0, x^2 - y = 0, x^2*y - 1.0 ≥ 0, x + y - 1.0 ≥ 0 }"
        @test sprint(show, MIME"text/plain"(), S) == "Basic semialgebraic Set defined by 2 equalities\n x - y = 0\n x^2 - y = 0\n2 inequalities\n x^2*y - 1.0 ≥ 0\n x + y - 1.0 ≥ 0\n"
        @test S.V isa AlgebraicSet{Float64}
        @test sprint(show, S.V) == "{ (x, y) | x - y = 0, x^2 - y = 0 }"
        @test sprint(show, MIME"text/plain"(), S.V) == "Algebraic Set defined by 2 equalities\n x - y = 0\n x^2 - y = 0\n"
        @test S === MultivariatePolynomials.changecoefficienttype(S, Float64)
        @test S.V === MultivariatePolynomials.changecoefficienttype(S.V, Float64)
        @test S.V.I === convert(typeof(S.V.I), S.V.I)
        @test BasicSemialgebraicSet{Int, polynomialtype(x, Int)}() isa BasicSemialgebraicSet{Float64, polynomialtype(x, Float64)}
        @test Int32(2)*x^2*y isa MultivariatePolynomials.AbstractTerm{Int32}
        Sf = MultivariatePolynomials.changecoefficienttype(S, Float32)
        @test Sf isa BasicSemialgebraicSet{Float32}
        @test Sf.V isa AlgebraicSet{Float32}

        @testset "Mixed types" begin
            S = (@set Int32(2)*x^2*y == 0 && 1.0x^2*y >= 0 && (6//3)*x^2*y == -y && 1.5x+y >= 0)
            S2 = S ∩ V
            S3 = V ∩ S
            @test inequalities(S2) == inequalities(S3) == S.p
            @test equalities(S2) == S3.V.I.p
        end

        T = (@set x*y^2 == -1 && x^2 + y^2 <= 1)
        V2 = @set T.V && V && x + y == 2.0
        @test V2 isa AlgebraicSet
        @test V2.I.p == [equalities(T); equalities(V); x + y - 2.0]
        S4 = @set S && T
        @test S4.p == [S.p; inequalities(T)]
        @test equalities(S4) == [S.V.I.p; T.V.I.p]

        @testset "Different variables" begin
            T = (@set x == x^2 && y <= y^2)
            @test sprint(show, T) == "{ (x, y) | -x^2 + x = 0, y^2 - y ≥ 0 }"
            @test sprint(show, MIME"text/plain"(), T) == "Basic semialgebraic Set defined by 1 equalitty\n -x^2 + x = 0\n1 inequalitty\n y^2 - y ≥ 0\n"
        end
    end
    @testset "Basic with no equality" begin
        S = @set x + y ≥ 1 && x ≤ y
        @test sprint(show, S) == "{ (x, y) | x + y - 1 ≥ 0, -x + y ≥ 0 }"
        @test sprint(show, MIME"text/plain"(), S) == """
Basic semialgebraic Set defined by no equality
2 inequalities
 x + y - 1 ≥ 0
 -x + y ≥ 0
"""
        @test sprint(show, S.V) == "R^n"
        @test sprint(show, MIME"text/plain"(), S.V) == """
Algebraic Set defined by no equality
"""
        @test S isa BasicSemialgebraicSet{Int}
        @test S === MultivariatePolynomials.changecoefficienttype(S, Int)
        @test S.V isa FullSpace
        Sf = MultivariatePolynomials.changecoefficienttype(S, Float64)
        @test Sf isa BasicSemialgebraicSet{Float64}
        @test Sf.V isa FullSpace

        @test S ∩ FullSpace() === S
        @test S.V ∩ FullSpace() === S.V
        @test FullSpace() ∩ S  === S
        @test FullSpace() ∩ S.V  === S.V
        @test S.V === MultivariatePolynomials.changecoefficienttype(S.V, Float64)
    end
    @testset "Basic with fixed variables" begin
        S = @set x == 1 && x ≤ x^2
        @test sprint(show, S) == "{ (x) | x - 1 = 0, x^2 - x ≥ 0 }"
        @test sprint(show, MIME"text/plain"(), S) == "Basic semialgebraic Set defined by 1 equalitty\n x - 1 = 0\n1 inequalitty\n x^2 - x ≥ 0\n"
        @test sprint(show, S.V) == "{ (x) | x - 1 = 0 }"
        @test sprint(show, MIME"text/plain"(), S.V) == "Algebraic Set defined by 1 equalitty\n x - 1 = 0\n"

        S = @set x == 1 && x ≤ y && 2 == y
        @test S isa BasicSemialgebraicSet{Int}
        @test S.V isa FixedVariablesSet{<:AbstractVariable, Int}
        @test rem(x + y, ideal(S.V)) == 3
        @test S === MultivariatePolynomials.changecoefficienttype(S, Int)
        @test S.V === MultivariatePolynomials.changecoefficienttype(S.V, Int)
        @test S.V.ideal === convert(typeof(S.V.ideal), S.V.ideal)
        Sf = MultivariatePolynomials.changecoefficienttype(S, Float64)
        @test Sf isa BasicSemialgebraicSet{Float64}
        @test Sf.V isa FixedVariablesSet{<:AbstractVariable, Float64}
        @test rem(x + y, ideal(Sf.V)) == 3

        S = @set x == 1 && x ≤ y && 2 == y && 1 == x
        @test S isa BasicSemialgebraicSet{Int}
        @test S.V isa FixedVariablesSet{<:AbstractVariable, Int}
        @test rem(x + y, ideal(S.V)) == 3
        Sf = MultivariatePolynomials.changecoefficienttype(S, Float64)
        @test Sf isa BasicSemialgebraicSet{Float64}
        @test Sf.V isa FixedVariablesSet{<:AbstractVariable, Float64}
        @test rem(x + y, ideal(Sf.V)) == 3

        @testset "Empty" begin
            for S in (@set(x == 1 && x ≤ y && 2 == y && 2 == x),
                      @set(x == 1 && x ≤ y && 2 == x),
                      @set(x == 1 && x ≤ y && 2 == x && y == 3),
                      @set(x ≤ y && y == 3 && x == 1 && y == 1),
                      @set(x ≤ y && y == 3 && x == 1 && y == 3 && x == 2))
                @test S isa BasicSemialgebraicSet{Int}
                @test S.V isa FixedVariablesSet{<:AbstractVariable, Int}
                @test isempty(S.V)
                @test iszero(length(S.V))
                @test isempty(collect(S.V))
                @test iszero(rem(x + y, ideal(S.V)))
            end
        end
    end
    @testset "Basic mixed fixed variables and equalities" begin
        for S in [@set(x == 1 && x ≤ y && 2 + y == x),
                  @set(x == 1 && 2 + y == x && x ≤ y),
                  @set(x ≤ y && 2 + y == x && x == 1),
                  @set(x ≤ y && x == 1 && 2 + y == x),
                  @set(2 + y == x && x ≤ y && x == 1),
                  @set(2 + y == x && x ≤ y && x == 1)]
            @test S isa BasicSemialgebraicSet{Float64}
            @test S.V isa AlgebraicSet{Float64}
        end

        solver = DummySolver()
        for S in [@set(x == 1 && x ≤ y && 2 + y == x, solver),
                  @set(x == 1 && 2 + y == x && x ≤ y, solver),
                  @set(x ≤ y && 2 + y == x && x == 1, solver),
                  @set(x ≤ y && x == 1 && 2 + y == x, solver),
                  @set(2 + y == x && x ≤ y && x == 1, solver),
                  @set(2 + y == x && x ≤ y && x == 1, solver)]
            @test S isa BasicSemialgebraicSet{Float64}
            @test S.V isa AlgebraicSet{Float64}
            @test S.V.solver === solver
        end
    end
end


text:
#!/usr/local/bin/julia

using PyPlot, JLD
# Plot the learning curve for the MNIST tutorial to verify
# the solver is progressing toward convergence

stats = load("snapshots/statistics.jld")
# println(typeof(stats))

tables = stats["statistics"]
ov = tables["obj_val"]
xy = sort(collect(ov))
datapoints = size(xy[:])
println("XY contains $datapoints datapoints")
x = [i for (i,j) in xy]
y = [j for (i,j) in xy]
x = convert(Array{Int64}, x)
y = convert(Array{Float64}, y)

raw = plot(x, y, linewidth=1, label="Raw")
xlabel("Iterations")
ylabel("Objective Value")
title("MNIST Learning Curve")
grid("on")

function low_pass{T <: Real}(x::Vector{T}, window::Int)
  len = length(x)
  y = Vector{Float64}(len)
  for i in 1:len
      # I want the mean of the first i terms up to width of window
      # Putting some numbers to this with window 4
      # i win lo  hi
      # 1  4  1   1
      # 2  4  1   2
      # 3  4  1   3
      # 4  4  1   4
      # 5  4  1   5
      # 6  4  2   6  => window starts to slide
      lo = max(1, i - window)
      hi = i
      y[i] = mean(x[lo:hi])
  end
  return y
end

window = Int64(round(length(xy)/4.0))
y_avg = low_pass(y, window)
avg = plot(x, y_avg, linewidth=2, label="Low Pass")
legend(handles=[raw; avg])
show()  #Shows the final plot


text: include("test.jl")


text:
include("functions.jl")
include("meixner_ortho_poly.jl")
include("hypergeometric.jl")
include("types.jl")


text:
#=
    The file for keeping track of intermediate results.
    It mirrors some f4 functions to collect statistics
    while computing the basis.
=#

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


text:
@inline function C_se(CellData,s,z,Def)
    """
    Concentration Solid-Electrolyte Transfer Function
    # Add Ins and Outs
        # Cell Data
        # Frequency Vector
        # Discretisation Locations
        # Electrode Definition
    """

 if Def == "Pos"
    Electrode = CellData.Pos #Electrode Length
 else
    Electrode = CellData.Neg #Electrode Length
 end

κ_eff = CellData.Const.κ*Electrode.ϵ_e^Electrode.κ_brug #Effective Electrolyte Conductivity
σ_eff = Electrode.σ*Electrode.ϵ_s^Electrode.σ_brug #Effective Electrode Conductivity

#Defining SOC
θ = CellData.Const.SOC * (Electrode.θ_100-Electrode.θ_0) + Electrode.θ_0

#Beta's
β = @. Electrode.Rs*sqrt(s/Electrode.Ds)

#Prepare for j0
cs0 = Electrode.cs_max * θ

#Current Flux Density
if CellData.Const.CellTyp == "Doyle_94"
    κ = Electrode.k_norm/Electrode.cs_max/CellData.Const.ce0^(1-Electrode.α)
    j0 = κ*(CellData.Const.ce0*(Electrode.cs_max-cs0))^(1-Electrode.α)*cs0^Electrode.α
else
    j0 = Electrode.k_norm*(CellData.Const.ce0*(cs0/Electrode.cs_max*(1-cs0/Electrode.cs_max)))^(1-Electrode.α)
end

#Resistances
Rtot = R*CellData.Const.T/(j0*F^2) + Electrode.RFilm

#∂Uocp_Def
∂Uocp_elc = CellData.Const.∂Uocp(Def,θ)/Electrode.cs_max

cse_res = -3/(Electrode.as*F*Electrode.L*CellData.Const.CC_A*Electrode.Rs) #Residual Variable for Pole Removal - eq. 4.45
ν = @. Electrode.L*sqrt((Electrode.as/σ_eff+Electrode.as/κ_eff)/(Rtot+∂Uocp_elc*(Electrode.Rs/(F*Electrode.Ds))*(tanh(β)/(tanh(β)-β)))) #Condensing Variable - eq. 4.13

cse_tf = @. ν*Electrode.Rs*(σ_eff*cosh(ν*z)+κ_eff*cosh(ν*(z-1)))*tanh(β)/(Electrode.as*F*Electrode.L*CellData.Const.CC_A*Electrode.Ds*(κ_eff+σ_eff)*sinh(ν)*(tanh(β)-β)) #Transfer Function - eq. 4.17
cse_tf = @. cse_tf-cse_res/s  #Pole removal - eq. 4.43
zero_tf = @. (5*Electrode.as*Electrode.Ds*F*Electrode.L^2*(κ_eff*(2-6*z+3*z^2)+(3*z^2-1)*σ_eff)-6*∂Uocp_elc*Electrode.Rs*κ_eff*σ_eff)/(30*CellData.Const.CC_A*Electrode.as*Electrode.Ds*∂Uocp_elc*F*Electrode.L*κ_eff*σ_eff) #For s = 0 / Wolfram Alpha
cse_tf[:,findall(s.==0)] .= zero_tf[:,findall(s.==0)]
D = zeros(length(z))
D_term = "zeros(length($z))"

if Def == "Pos"
    cse_tf = -cse_tf
    cse_res = -cse_res

    if Debug == 1
        println("D:C_se:Pos",D)
        println("z:C_se:Pos",z)
        println("zero_tf:C_se:Pos",zero_tf[:,1])
        println("κ_eff:C_se:Pos",κ_eff)
        println("σ_eff:C_se:Pos",σ_eff)
        #println("ν_∞:C_se:Pos",ν_∞)
        println("j0:C_se:Pos",j0)
        println("as:C_se:Pos",Electrode.as)
        println("Rtot:C_se:Pos",Rtot)
        println("σ_eff:C_se:Pos",σ_eff)
        println("∂Uocp_elc:C_se:Pos",∂Uocp_elc)
        println("Electrode.L:C_se:Pos",Electrode.L)
    end
end
cse_res = cse_res*ones(length(z))
return cse_tf, D, cse_res, D_term

end

max_stars_repo_name: BradyPlanden/RealisationAlgorithms

max_issues_repo_name: BradyPlanden/RealisationAlgorithms

max_forks_repo_name: BradyPlanden/RealisationAlgorithms


text:
module DataVoyager

using Electron, DataValues

import IteratorInterfaceExtensions, TableTraits, IterableTables, JSON, VegaLite

export Voyager

app = nothing

mutable struct Voyager
    w::Window

    function Voyager()
        main_html_uri = string("file:///", replace(joinpath(@__DIR__, "htmlui", "main.html"), '\\' => '/'))

        global app

        if app==nothing
            app = Application()
        end

        w = Window(app, URI(main_html_uri), options=Dict("title"=>"Data Voyager"))

        new(w)
    end
end

function (v::Voyager)(source)
    TableTraits.isiterabletable(source)===false && error("'source' is not a table.")

    it = IteratorInterfaceExtensions.getiterator(source)

    data_dict = Dict{String,Any}()

    data_dict["values"] = [Dict{Symbol,Any}(c[1]=>isa(c[2], DataValue) ? (isna(c[2]) ? nothing : get(c[2])) : c[2] for c in zip(keys(r), values(r))) for r in it]

    data = JSON.json(data_dict)

    code = "voyagerInstance.updateData($data)"

    run(v.w, code)

    return v
end

function Base.getindex(v::Voyager)
    code = "voyagerInstance.getSpec(true)"

    content = run(v.w, code)

    return VegaLite.VLSpec{:plot}(content)
end

# function Base.getindex(v::Voyager, index::Int)
#     code = "voyagerInstance.getBookmarkedSpecs()"

#     content = run(v.w, code)

#     info(content)

#     return VegaLite.VLSpec{:plot}(JSON.parse(content[index]))
# end

function Voyager(source)
    v = Voyager()

    v(source)

    return v
end

end # module