text:
"""
    Storage

This struct is meant for holding storage parameters for a for a blosc2 container, allowing to specify,
for example, how to interpret the contents included in the schunk.
"""
struct Storage
    contiguous ::Bool
    urlpath ::Union{String, Nothing}
    cparams ::Union{CompressionParams, Nothing}
    dparams ::Union{DecompressionParams, Nothing}
end

 """
    Storage(;
        contiguous = false, urlpath::Union{AbstractString, Nothing} = nothing,
        cparams::Union{CompressionParams, Nothing} = nothing ,
        dparams::Union{DecompressionParams, Nothing} = nothing)

Create `Storage` struct

# Arguments

    - `contiguous` - Whether the chunks are contiguous or sparse.

    - `urlpath` - The path for persistent storage. If `nothing`, that means in-memory.

    - `cparams` - The compression params when creating a schunk.
    If `nothing`, sensible defaults are used depending on the context.

    - `dparams` - The compression params when creating a schunk.
    If `nothing`, sensible defaults are used depending on the context.
"""
function Storage(;
    contiguous = false, urlpath::Union{AbstractString, Nothing} = nothing,
    cparams::Union{CompressionParams, Nothing} = nothing ,
    dparams::Union{DecompressionParams, Nothing} = nothing)
    Storage(contiguous, urlpath, cparams, dparams)
end

function with_preserved_blosc2_storage(f, s::Storage)
    cparams = isnothing(s.cparams) ? nothing : make_cparams(s.cparams)
    dparams = isnothing(s.dparams) ? nothing : make_dparams(s.dparams)
    GC.@preserve s cparams dparams begin
        bstorage = Lib.blosc2_storage(
            s.contiguous,
            isnothing(s.urlpath) ? C_NULL : Base.unsafe_convert(Cstring, s.urlpath),
            isnothing(cparams) ? C_NULL : Base.unsafe_convert(Ptr{Lib.blosc2_cparams}, Ref(cparams)),
            isnothing(dparams) ? C_NULL : Base.unsafe_convert(Ptr{Lib.blosc2_dparams}, Ref(dparams)),
            C_NULL
        )
        f(bstorage)
    end
end


text:
using Ogre
using WaterData
using PyPlot
using Colors
using DataStructures: OrderedDict

function mr_diagrams()
    ms = linspace(0.5M_earth, 10M_earth, 3)

    h2o = WaterData.load_full_eos()["grid"]
    funcs = WaterData.load_piecewise_eoses()
    mgsio3 = funcs["mgsio3"]
    fe = funcs["fe"]

    eoses = [fe, mgsio3, h2o]
    mass_fractions = [0.15, 0.30, 0.55]
    total_points = 200
    R_bracket = [0, 10] * R_earth

    # Cₚ_const = WaterData.ConstantHeatCapacity(4200)  # constant heat capacity for liquid water
    # Cₚ_inf = WaterData.ConstantHeatCapacity(Inf)  # infinite heat capacity
    Cₚ = WaterData.load_heat_capacity()["heatcap_h2o"] # realistic h.c.

    isotherms = OrderedDict(# 300 => WaterData.slice(h2o, 300),
                            400 => WaterData.slice(h2o, 400),
                            500 => WaterData.slice(h2o, 500),
                            600 => WaterData.slice(h2o, 600),
                            700 => WaterData.slice(h2o, 700),
                            800 => WaterData.slice(h2o, 800),
                            # 900 => WaterData.slice(h2o, 900),
                            # 1000 => WaterData.slice(h2o, 1000)
                            )

    function R_adiabat(M, Tsurf, Psurf)
        Me = M / M_earth
        bvs = Ogre.ValueSet(M, R_earth, Psurf, Tsurf)
        grid = linspace(M, 0, total_points)
        planet = Ogre.PlanetSystem(M, h2o, Cₚ, bvs, grid, R_bracket)

        Ogre.find_radius!(planet)
    end

    function R_adiabat_cored(M, Tsurf, Psurf)
        Me = M / M_earth
        bvs = Ogre.ValueSet(M, R_earth, Psurf, Tsurf)
        grid = linspace(M, 0, total_points)
        coredeos = Ogre.MassPiecewiseEOS(eoses, M, mass_fractions)
        planet = Ogre.PlanetSystem(M, coredeos, Cₚ, bvs, grid, R_bracket)

        Ogre.find_radius!(planet)
    end

    function R_isotherm(M, Tsurf, Psurf)
        bvs = Ogre.ValueSet(M, R_earth, Psurf)
        grid = linspace(M, 0, total_points)
        planet = Ogre.PlanetSystem(M, isotherms[Tsurf], bvs, grid, R_bracket)

        Ogre.find_radius!(planet)
    end

    temps = keys(isotherms)
    dropfirst(x) = drop(x, 1)

    for f in [R_isotherm, R_adiabat, R_adiabat_cored], Psurf in [1e5, 1e6, 1e7]
        cmap = Colors.colormap("Reds", length(temps) + 1) |> dropfirst
        figure()
        for (T, color) in zip(temps, cmap)
            c = map(f -> f(color), [red, green, blue])
            @time rs = map(M -> f(M, T, Psurf), ms)
            plot(ms ./ M_earth, rs ./ R_earth, color=c, label="$T K")
        end

        method = string(f)[3:end]
        log10P = floor(Int, log10(Psurf))
        title("Psurf=1e$log10P Pa; $method")
        xlabel(L"Mass / M$_⊕$")
        ylabel(L"Radius / R$_⊕$")
        legend(loc="best")
        xlim((0, 10))
        ylim((1, 5))
        tight_layout()
        savefig("$(method)_$(log10P)Pa.png")
        close()
    end
end

function phaseplots(M, Psurf, Tsurf)
    M = M * M_earth
    R = cbrt(M)
    npoints = 200
    grid = linspace(M, 0, npoints)
    eos = WaterData.load_full_eos()["grid"]
    heatcap = WaterData.load_heat_capacity()["heatcap_h2o"]

    bcs = Ogre.ValueSet(M, R, Psurf, Tsurf)
    sys = Ogre.PlanetSystem(M, eos, heatcap, bcs, grid)

    soln = Ogre.find_structure!(sys)
    plot(soln, sys)
    # tight_layout()
    Ogre.phaseplot(soln)
end

close(:all)
phaseplots(1, 1e6, 400)
# mr_diagrams()

ext: jl


text:
using FHIRClientJSON
using Test

import BSON
import BrokenRecord
import HTTP
import JSON3

@testset "FHIRClientJSON.jl" begin
    include("unit-tests.jl")
    include("recorded-tests.jl")
end


text:
using Test
using SetProg

@testset "apply_matrix" begin
    SetProg.DynamicPolynomials.@polyvar x[1:2]
    SetProg.DynamicPolynomials.@polyvar y[1:2]
    SetProg.DynamicPolynomials.@polyvar z[1:1]
    Q = [1 2 4
         2 3 5
         4 5 6]
    q = SetProg.GramMatrix(Q, SetProg.MP.monomials(x, 2))
    p = SetProg.MP.polynomial(q)
    @testset "2x2 Float64" begin
        B = [2.0 3.0
             4.0 5.0]
        qB = SetProg.apply_matrix(q, B, y, 2)
        @test qB isa SetProg.GramMatrix{Float64}
        @test qB == p(x => B * y)
    end
    @testset "2x1 Int" begin
        A = reshape([2, 3], 2, 1)
        qA = SetProg.apply_matrix(q, A, z, 2)
        @test qA isa SetProg.GramMatrix{Int}
        @test qA == p(x => A * z)
    end
end


text:
#
# rules.jl -
#
# Implement rules for automatically simplifying expressions involving mappings.
#
#-------------------------------------------------------------------------------
#
# This file is part of LazyAlgebra (https://github.com/emmt/LazyAlgebra.jl)
# released under the MIT "Expat" license.
#
# Copyright (c) 2017-2021 Éric Thiébaut.
#

#------------------------------------------------------------------------------
# NEUTRAL ELEMENTS

# The neutral element ("zero") for the addition is zero times a mapping of the
# proper type.
zero(A::Mapping) = 0*A

iszero(A::Scaled) = iszero(multiplier(A))
iszero(::Mapping) = false

# The neutral element ("one") for the composition is the identity.
const Id = Identity()
one(::Mapping) = Id

isone(::Identity) = true
isone(::Mapping) = false

#------------------------------------------------------------------------------
# UNQUALIFIED OUTER CONSTRUCTORS

# Unqualified outer constructors are provided which call the corresponding
# inner constructors with all suitable parameters and rely on the inner
# constructors to check whether the call was allowed or not.  A constraint that
# must hold is that T(A), with T an unqualified type constructor, always yields
# an instance of T.
Direct(A::Mapping) = A # provided for completeness
Adjoint(A::T) where {T<:Mapping} = Adjoint{T}(A)
Inverse(A::T) where {T<:Mapping} = Inverse{T}(A)
InverseAdjoint(A::T) where {T<:Mapping} = InverseAdjoint{T}(A)
Gram(A::T) where {T<:Mapping} = Gram{T}(A)
Jacobian(A::M, x::T) where {M<:Mapping,T} = Jacobian{M,T}(A, x)
Scaled(α::S, A::T) where {S<:Number,T<:Mapping} = Scaled{T,S}(α, A)
Sum(ops::Mapping...) = Sum(ops)
Sum(ops::T) where {N,T<:NTuple{N,Mapping}} = Sum{N,T}(ops)
Composition(ops::Mapping...) = Composition(ops)
Composition(ops::T) where {N,T<:NTuple{N,Mapping}} = Composition{N,T}(ops)

# Qualified outer constructors to forbid decoration of mappings of specific
# types when, according to the simplification rules, another more simple
# construction should be built instead.  Everything not forbidden is allowed
# except that additional tests may be performed by the inner constructors
# (e.g., Adjoint check that its argument is linear).
#
# FIXME: Some restrictions may be bad ideas like adjoint of sums or
# compositions.

for (func, blacklist) in ((:Adjoint,        (:Identity,
                                             :Adjoint,
                                             :Inverse,
                                             :InverseAdjoint,
                                             :Scaled,
                                             :Sum,
                                             :Composition)),
                          (:Inverse,        (:Identity,
                                             :Adjoint,
                                             :Inverse,
                                             :InverseAdjoint,
                                             :Scaled,
                                             :Composition)),
                          (:InverseAdjoint, (:Identity,
                                             :Adjoint,
                                             :Inverse,
                                             :InverseAdjoint,
                                             :Scaled,
                                             :Composition)),
                          (:Gram,           (:Inverse,
                                             :InverseAdjoint,
                                             :Scaled,)),
                          (:Jacobian,       (:Scaled,)),
                          (:Scaled,         (:Scaled,)))
    for T in blacklist
        if func === :Scaled
            @eval $func{T,S}(α::S, A::T) where {S<:Number,T<:$T} =
                illegal_call_to($func, T)
        elseif func === :Jacobian
            @eval $func{M,T}(A::M, x::T) where {M<:$T,T} =
                illegal_call_to($func, M)
        else
            @eval $func{T}(A::T) where {T<:$T} = illegal_call_to($func, T)
        end
    end
end

@noinline illegal_call_to(::Type{Adjoint}, T::Type) =
    bad_argument("the `Adjoint` constructor cannot be applied to an instance of `",
                 brief(T), "`, use expressions like `A'` or `adjoint(A)`")

@noinline illegal_call_to(::Type{Inverse}, T::Type) =
    bad_argument("the `Inverse` constructor cannot be applied to an instance of `",
                 brief(T), "`, use expressions like `A\\B`, `A/B` or `inv(A)`")

@noinline illegal_call_to(::Type{InverseAdjoint}, T::Type) =
    bad_argument("the `InverseAdjoint` constructor cannot be applied to an instance of `",
                 brief(T), "`, use expressions like `A'\\B`, `A/(B')`, `inv(A')` or `inv(A)'`")

@noinline illegal_call_to(::Type{Gram}, T::Type) =
    bad_argument("the `Gram` constructor cannot be applied to an instance of `",
                 brief(T), "`, use expressions like `A'*A` or `gram(A)`")

@noinline illegal_call_to(::Type{Jacobian}, T::Type) =
    bad_argument("the `Jacobian` constructor cannot be applied to an instance of `",
                 brief(T), "`, use an expression like `∇(A,x)`")

@noinline illegal_call_to(::Type{Scaled}, T::Type) =
    bad_argument("the `Scaled` constructor cannot be applied to an instance of `",
                 brief(T), "`, use expressions like `α*A`")

brief(::Type{<:Adjoint}       ) = "Adjoint"
brief(::Type{<:Inverse}       ) = "Inverse"
brief(::Type{<:InverseAdjoint}) = "InverseAdjoint"
brief(::Type{<:Gram}          ) = "Gram"
brief(::Type{<:Jacobian}      ) = "Jacobian"
brief(::Type{<:Scaled}        ) = "Scaled"
brief(::Type{<:Sum}           ) = "Sum"
brief(::Type{<:Composition}   ) = "Composition"
brief(::Type{<:Identity}      ) = "Identity"
brief(T::Type) = repr(T)

#------------------------------------------------------------------------------
# SCALED TYPE

# Left-multiplication and left-division by a scalar.  The only way to
# right-multiply or right-divide a mapping by a scalar is to right multiply or
# divide it by the scaled identity.
*(α::Number, A::Mapping) = (α == 1 ? A : Scaled(α, A))
*(α::Number, A::Scaled) = (α*multiplier(A))*unscaled(A)
\(α::Number, A::Mapping) = inv(α)*A
\(α::Number, A::Scaled) = (multiplier(A)/α)*unscaled(A)
/(α::Number, A::Mapping) = α*inv(A)

#------------------------------------------------------------------------------
# ADJOINT TYPE

# Adjoint for non-specific mappings.
adjoint(A::Mapping) = _adjoint(LinearType(A), A)

_adjoint(::Linear, A::Mapping) = _adjoint(Linear(), SelfAdjointType(A), A)
_adjoint(::Linear, ::SelfAdjoint, A::Mapping) = A
_adjoint(::Linear, ::NonSelfAdjoint, A::Mapping) = Adjoint(A)
_adjoint(::NonLinear, A::Mapping) =
    throw_forbidden_adjoint_of_non_linear_mapping()

# Adjoint for specific mapping types.
adjoint(A::Identity) = Id
adjoint(A::Scaled) = conj(multiplier(A))*adjoint(unscaled(A))
adjoint(A::Adjoint) = unveil(A)
adjoint(A::Inverse) = inv(adjoint(unveil(A)))
adjoint(A::InverseAdjoint) = inv(unveil(A))
adjoint(A::Jacobian) = Jacobian(A)
adjoint(A::Gram) = A
adjoint(A::Composition) =
    # It is assumed that the composition has already been simplified, so we
    # just apply the mathematical formula for the adjoint of a composition.
    Composition(reversemap(adjoint, terms(A)))

function adjoint(A::Sum{N}) where {N}
    # It is assumed that the sum has already been simplified, so we just apply
    # the mathematical formula for the adjoint of a sum and sort the resulting
    # terms.
    B = Vector{Mapping}(undef, N)
    @inbounds for i in 1:N
        B[i] = adjoint(A[i])
    end
    return Sum(to_tuple(sort!(B)))
end

@noinline throw_forbidden_adjoint_of_non_linear_mapping() =
    bad_argument("taking the adjoint of non-linear mappings is not allowed")

#------------------------------------------------------------------------------
# JACOBIAN

"""
    ∇(A, x)

yields a result corresponding to the Jacobian (first partial derivatives) of
the linear mapping `A` for the variables `x`.  If `A` is a linear mapping,
`A` is returned whatever `x`.

The call

    jacobian(A, x)

is an alias for `∇(A,x)`.

"""
∇(A::Mapping, x) = jacobian(A, x)
jacobian(A::Mapping, x) = _jacobian(LinearType(A), A, x)
_jacobian(::Linear, A::Mapping, x) = A
_jacobian(::NonLinear, A::Mapping, x) = Jacobian(A, x)
jacobian(A::Scaled, x) = multiplier(A)*jacobian(unscaled(A), x)

@doc @doc(∇) jacobian

#------------------------------------------------------------------------------
# INVERSE TYPE

# Inverse for non-specific mappings (a simple mapping or a sum or mappings).
inv(A::T) where {T<:Mapping} = Inverse{T}(A)

# Inverse for specific mapping types.
inv(A::Identity) = Id
inv(A::Scaled) = (is_linear(unscaled(A)) ?
                  inv(multiplier(A))*inv(unscaled(A)) :
                  inv(unscaled(A))*(inv(multiplier(A))*Id))
inv(A::Inverse) = unveil(A)
inv(A::InverseAdjoint) = adjoint(unveil(A))
inv(A::Adjoint) = InverseAdjoint(unveil(A))
inv(A::Composition) =
    # Even though the composition has already been simplified, taking the
    # inverse may trigger other simplifications, so we must rebuild the
    # composition term by term in reverse order (i.e. applying the mathematical
    # formula for the inverse of a composition).
    _merge_inv_mul(terms(A))

# `_merge_inv_mul([A,i,]B)` is recursively called to build the inverse of a
# composition.  Argument A is a mapping (initially not specified or the
# identity) of the resulting composition, argument `i` is the index of the next
# component to take (initially not specified or set to `N` the number of
# terms), argument `B` is a tuple (initially full) of the remaining terms.
_merge_inv_mul(B::NTuple{N,Mapping}) where {N} =
    # Initialize recursion.
    _merge_inv_mul(inv(last(B)), N - 1, B)

function _merge_inv_mul(A::Mapping, i::Int, B::NTuple{N,Mapping}) where {N}
    # Perform intermediate and last recursion step.
    C = A*inv(B[i])
    return (i > 1 ? _merge_inv_mul(C, i - 1, B) : C)
end

#------------------------------------------------------------------------------
# SUM OF MAPPINGS

# Unary minus and unary plus.
-(A::Mapping) = (-1)*A
-(A::Scaled) = (-multiplier(A))*unscaled(A)
+(A::Mapping) = A

# Subtraction.
-(A::Mapping, B::Mapping) = A + (-B)

# Simplify the sum of two mappings.
+(A::Mapping, B::Mapping) = add(A, B)

"""
    add(A, B)

performs the final stage of simplifying the sum `A + B` of mappings `A` and
`B`.  This method assumes that any other simplifications than those involving
sum of sums have been performed on `A + B` and that `A` and `B` have already
been simplified (individually).

Trivial simplifications of the composition `A + B` of mappings `A` and `B` must
be done by specializing the operator `+` and `sum(A,B)` may eventually be
called to simplify the sum of `A` and `B` when one of these may be a
sum.  This method just returns `Sum(A,B)` when none of its
arguments is a composition.

yields a simplified sum `A + B` of the mappings `A` and `B`.  This helper
function is intended to be called when at least one of `A` or `B` is a sum.

The ability to perform simplifications relies on implemented specializations of
`A + B` when neither `A` nor `B` are sums.  It is also assumed that `A` and `B`
have already been simplified if they are sums.

"""
add(A::Sum,     B::Mapping) = _add(A, B)
add(A::Mapping, B::Sum    ) = _add(A, B)
add(A::Sum,     B::Sum    ) = _add(A, B)
add(A::Mapping, B::Mapping) = begin
    # Neither `A` nor `B` is a sum.
    if identical(unscaled(A), unscaled(B))
        return (multiplier(A) + multiplier(B))*unscaled(A)
    elseif identifier(A) ≤ identifier(B)
        return Sum(A, B)
    else
        return Sum(B, A)
    end
end

_add(A::Mapping, B::Mapping) = begin
    V = add!(as_vector(+, A), B)
    length(V) == 1 ? V[1] : Sum(to_tuple(V))
end

# Add the terms of a sum one-by-one.  Since terms must be re-ordered, there are
# no obvious better ways to recombine.
function add!(A::Vector{Mapping}, B::Sum{N}) where {N}
    @inbounds for i in 1:N
        add!(A, B[i])
    end
    return A
end

function add!(A::Vector{Mapping}, B::Mapping)
    # Nothing to do if B is zero times anything.
    multiplier(B) == 0 && return A

    # If exact match found, update A in-place and return.
    n = length(A)
    @inbounds for i in 1:n
        if identical(unscaled(A[i]), unscaled(B))
            λ = multiplier(A[i]) + multiplier(B)
            if λ == 1
                A[i] = unscaled(B)
            elseif λ != 0
                A[i] = λ*unscaled(B)
            else
                # Multiplier is zero. Drop term if there are other terms
                # or keep the single term times zero.
                if n > 1
                    for j in i:n-1
                        A[j] = A[j+1]
                    end
                    resize!(A, n - 1)
                else
                    A[1] = 0*unscaled(A[1])
                end
            end
            return A
        end
    end

    # If no exact match found, insert B in A in order.
    id = identifier(B)
    i = 1
    while i ≤ n && @inbounds(identifier(A[i])) < id
        i += 1
    end
    resize!(A, n + 1)
    @inbounds for j in n:-1:i
        A[j+1] = A[j]
    end
    A[i] = B
    return A
end

#------------------------------------------------------------------------------
# COMPOSITION OF MAPPINGS

# Left and right divisions.
\(A::Mapping, B::Mapping) = inv(A)*B
/(A::Mapping, B::Mapping) = A*inv(B)

# Dot operator (\cdot + tab) involving a mapping acts as the multiply or
# compose operator.
⋅(A::Mapping, B::Mapping) = A*B
⋅(A::Mapping, B::Any    ) = A*B
⋅(A::Any,     B::Mapping) = A*B

# Compose operator (\circ + tab) beween mappings.
∘(A::Mapping, B::Mapping) = A*B

# Rules for the composition of 2 mappings.  Mappings that may behave
# specifically in a composition have type `Identity`, `Scaled` and
# `Composition`; all others have the same behavior.

# Composition with identity.
*(::Identity, ::Identity) = Id
for T in (Scaled, Composition, Sum, Mapping)
    @eval begin
        *(::Identity, A::$T) = A
        *(A::$T, ::Identity) = A
    end
end

# Simplify the composition of two mappings (including compositions).
*(A::Mapping, B::Mapping) = compose(A, B)

# Simplify compositions involving a scaled mapping.
*(A::Scaled, B::Mapping) = multiplier(A)*(unscaled(A)*B)
*(A::Mapping, B::Scaled) =
    if is_linear(A)
        multiplier(B)*(A*unscaled(B))
    else
        compose(A, B)
    end
*(A::Scaled, B::Scaled) =
    if is_linear(A)
        (multiplier(A)*multiplier(B))*(unscaled(A)*unscaled(B))
    else
        multiplier(A)*(unscaled(A)*B)
    end

# Simplify compositions involving an inverse mapping.
*(A::Inverse{T}, B::T) where {T<:Mapping} =
    identical(unveil(A), B) ? Id : compose(A, B)
*(A::T, B::Inverse{T}) where {T<:Mapping} =
    identical(A, unveil(B)) ? Id : compose(A, B)
*(A::Inverse, B::Inverse) = compose(A, B)
*(A::InverseAdjoint{T}, B::Adjoint{T}) where {T<:Mapping} =
    identical(unveil(A), unveil(B)) ? Id : compose(A, B)
*(A::Adjoint{T}, B::InverseAdjoint{T}) where {T<:Mapping} =
    identical(unveil(A), unveil(B)) ? Id : compose(A, B)
*(A::InverseAdjoint, B::InverseAdjoint) = compose(A, B)

# Automatically build Gram operators, Gram(A) ≡ A'*A.  The following automatic
# rules are implemented (for an "allowed" linear mapping A):
#
#     A'*A -> Gram(A)
#     A*A' -> Gram(A')
#     inv(A)*inv(A') -> inv(A'*A) -> inv(Gram(A))
#     inv(A')*inv(A) -> inv(A*A') -> inv(Gram(A'))
#
# other rules implemented elsewhere:
#
#     Gram(inv(A))  -> inv(Gram(A'))
#     Gram(inv(A')) -> inv(Gram(A))
#
# In principle, if forming the adjoint has been allowed, it is not needed to
# check whether operands are linear mappings.
*(A::Adjoint{T}, B::T) where {T<:Mapping} =
    identical(unveil(A), B) ? Gram(B) : compose(A, B)
*(A::T, B::Adjoint{T}) where {T<:Mapping} =
    identical(A, unveil(B)) ? Gram(B) : compose(A, B)
*(A::Inverse{T}, B::InverseAdjoint{T}) where {T<:Mapping} =
    identical(unveil(A), unveil(B)) ? Inverse(Gram(unveil(A))) :
    compose(A, B)
*(A::InverseAdjoint{T}, B::Inverse{T}) where {T<:Mapping} =
    identical(unveil(A), unveil(B)) ?
    Inverse(Gram(Adjoint(unveil(A)))) : compose(A, B)

"""
    compose(A,B)

performs the final stage of simplifying the composition `A*B` of mappings `A`
and `B`.  This method assumes that any other simplifications than those
involving composition of compositions have been performed on `A*B` and that `A`
and `B` have already been simplified (individually).

Trivial simplifications of the composition `A*B` of mappings `A` and `B` must
be done by specializing the operator `*` and `compose(A,B)` may eventually be
called to simplify the composition of `A` and `B` when one of these may be a
composition.  This method just returns `Composition(A,B)` when neither `A` nor
`B` is a composition.

""" compose

# Compose two mappings when at least one is a composition or when none is a
# composition.
compose(A::Composition, B::Composition) = _compose(A, B)
compose(A::Composition, B::Mapping    ) = _compose(A, B)
compose(A::Mapping,     B::Composition) = _compose(A, B)
compose(A::Mapping,     B::Mapping    ) = Composition(A, B)

_compose(A::Mapping, B::Mapping) = begin
    C = compose!(as_vector(*, A), B)
    n = length(C)
    return (n == 0 ? Id :
            n == 1 ? C[1] :
            Composition(to_tuple(C)))
end

"""
    compose!(A, B) -> A

overwrites `A` with a simplified composition of a left operand `A` and a right
operand `B`.  The left operand is a composition of (zero or more) mappings
whose terms are stored in the vector of mappings `A` (if `A` is empty, the left
operand is assumed to be the identity).  On return, the vector `A` is modified
to store the terms of the simplified composition of `A` and `B`.  The left
operand `B` may be itself a composition (as an instance of
`LazyAlgebra.Composition` or as a vector of mappings) or any other kind of
mapping.

""" compose!

function compose!(A::Vector{Mapping}, B::Composition{N}) where {N}
    @inbounds for i in 1:N
        # Build the simplified composition A*B[i].
        compose!(A, B[i])
        if identical(last(A), B[i])
            # The last term of the simplified composition A*B[i] is still B[i],
            # which indicates that composing A with B[i] did not yield any
            # simplifications.  It is sufficient to append all the other terms
            # of B to A as no further simplifications are expected.
            return append_terms!(A, B, (i+1):N)
        end
    end
    return A
end

function compose!(A::Vector{Mapping}, B::Mapping)
    # Compute the simplified composition of the last term of A with B.  The
    # result is either a simple mapping or a simplified composition.
    m = length(A); @certify m > 0
    C = A[m]*B

    # Replace the last term of A with C.
    if C isa Composition && identical(C[1], A[m])
        # Nothing has changed at the tail of the composition A.  No further
        # simplifications are expected.  Push all terms of C to A, but the
        # first term of C which is identical to the last term of A.
        append_terms!(A, C, 2:length(C))
    elseif m > 1
        # Drop the last term of A and compose the remaining terms with C.  This
        # may trigger further simplifications
        compose!(resize!(A, m - 1), C)
    elseif C isa Composition
        # Replace the only term of A by all the terms of the composition C.
        # This is the same as above but avoids calling `resize!` as a small
        # optimization.
        A[1] = C[1]
        append_terms!(A, C, 2:length(C))
    else
        # Replace the only term of A by the simple mapping C.
        A[1] = C
    end
    return A
end

#------------------------------------------------------------------------------
# UTILITIES FOR BUILDING SUMS AND COMPOSITIONS

"""
    as_vector(op, A)

yields a vector of mappings with the terms of the mapping `A`.  Argument `op`
is `+` or `*`.  If `op` is `+` (resp. `*`) and `A` is a sum (resp. a
composition) of mappings, the terms of `A` are extracted in the returned
vector; otherwise, the returned vector has just one element which is `A`.

"""
function as_vector(::Union{typeof(+),typeof(*)}, A::Mapping)
    V = Vector{Mapping}(undef, 1)
    V[1] = A
    return V
end

as_vector(::typeof(+), A::Sum        ) = collect_terms(A)
as_vector(::typeof(*), A::Composition) = collect_terms(A)

"""
    collect_terms(A)

collects the terms of the mapping `A` into a vector.  This is similar to
`collect(A)` except that the element type of the result is forced to be
`Mapping`.

"""
function collect_terms(A::Union{Sum{N},Composition{N}}) where {N}
    V = Vector{Mapping}(undef, N)
    @inbounds for i in 1:N
        V[i] = A[i]
    end
    return V
end

"""
    append_terms!(A, B, I=1:length(B)) -> A

pushes all terms `B[i]` for all `i ∈ I` to `A` and returns `A`.

"""
function append_terms!(A::Vector{Mapping},
                       B::Union{Vector{Mapping},Composition},
                       I::AbstractUnitRange{<:Integer} = Base.OneTo(length(B)))

    imin, imax = Int(first(I)), Int(last(I))
    (1 ≤ imin && imax ≤ length(B)) ||
        bad_argument("out of bounds indices in given range")
    if imin ≤ imax
        m = length(A)
        n = imax - imin + 1
        resize!(A, m + n)
        k = m + 1 - imin
        @inbounds for i in I
            A[k+i] = B[i]
        end
    end
    return A
end


text:
# ---
# title: 98. Validate Binary Search Tree
# id: problem98
# author: Indigo
# date: 2021-06-02
# difficulty: Medium
# categories: Tree, Depth-first Search, Recursion
# link: <https://leetcode.com/problems/validate-binary-search-tree/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, _determine if it is a valid binary search
# tree (BST)_.
#
# A **valid BST** is defined as follows:
#
#   * The left subtree of a node contains only nodes with keys **less than** the node's key.
#   * The right subtree of a node contains only nodes with keys **greater than** the node's key.
#   * Both the left and right subtrees must also be binary search trees.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/12/01/tree1.jpg)
#
#
#
#     Input: root = [2,1,3]
#     Output: true
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/12/01/tree2.jpg)
#
#
#
#     Input: root = [5,1,4,null,null,3,6]
#     Output: false
#     Explanation: The root node's value is 5 but its right child's value is 4.
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[1, 104]`.
#   * `-231 <= Node.val <= 231 - 1`
#
#
## @lc code=start
using LeetCode

function is_valid_BST(root)
    isnothing(root) && return true
    arr = Int[]
    function pre_tra(root::TreeNode{Int})
        pre_tra(root.left)
        push!(arr, root.val)
        pre_tra(root.right)
    end
    pre_tra(::Nothing) = nothing
    pre_tra(root)
    return issorted(arr)
end
## @lc code=end


text:
"""
    GaborKernel(; ell::Real=1.0, p::Real=1.0)

Gabor kernel with lengthscale `ell` and period `p`. Given by
```math
    κ(x,y) =  h(x-z), h(t) = exp(-sum(t.^2./(ell.^2)))*cos(pi*sum(t./p))
```

"""
struct GaborKernel{K<:Kernel} <: Kernel
    kernel::K
    function GaborKernel(;ell=nothing, p=nothing)
        k = _gabor(ell=ell, p=p)
        return new{typeof(k)}(k)
    end
end

@functor GaborKernel

(κ::GaborKernel)(x, y) = κ.kernel(x ,y)

function _gabor(; ell = nothing, p = nothing)
    if ell === nothing
        if p === nothing
            return SqExponentialKernel() * CosineKernel()
        else
            return SqExponentialKernel() * transform(CosineKernel(), 1 ./ p)
        end
    elseif p === nothing
        return transform(SqExponentialKernel(), 1 ./ ell) * CosineKernel()
    else
        return transform(SqExponentialKernel(), 1 ./ ell) * transform(CosineKernel(), 1 ./ p)
    end
end

function Base.getproperty(k::GaborKernel, v::Symbol)
    if v == :kernel
        return getfield(k, v)
    elseif v == :ell
        kernel1 = k.kernel.kernels[1]
    	if kernel1 isa TransformedKernel
    	    return 1 ./ kernel1.transform.s[1]
    	else
    	    return 1.0
    	end
    elseif v == :p
        kernel2 = k.kernel.kernels[2]
        if kernel2 isa TransformedKernel
            return 1 ./ kernel2.transform.s[1]
        else
            return 1.0
        end
    else
        error("Invalid Property")
    end
end

Base.show(io::IO, κ::GaborKernel) = print(io, "Gabor Kernel (ell = ", κ.ell, ", p = ", κ.p, ")")

kernelmatrix(κ::GaborKernel, x::AbstractVector) = kernelmatrix(κ.kernel, x)

function kernelmatrix(κ::GaborKernel, x::AbstractVector, y::AbstractVector)
    return kernelmatrix(κ.kernel, x, y)
end

kerneldiagmatrix(κ::GaborKernel, x::AbstractVector) = kerneldiagmatrix(κ.kernel, x)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Fontconfig"
version = v"2.13.1"

# Collection of sources required to build Fontconfig
sources = [
    "https://www.freedesktop.org/software/fontconfig/release/fontconfig-2.13.1.tar.gz" =>
    "9f0d852b39d75fc655f9f53850eb32555394f36104a044bb2b2fc9e66dbbfa7f",

]

# Bash recipe for building across all platforms
script = raw"""
	if [ $target = "aarch64-linux-gnu" ]; then
	  ln -s /opt/aarch64-linux-gnu/aarch64-linux-gnu/sys-root/lib/ld-linux-aarch64.so.1 /lib/ld-linux-aarch64.so.1
	  export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/x86_64-linux-gnu/x86_64-linux-gnu/sys-root/lib64:/opt/aarch64-linux-gnu/aarch64-linux-gnu/sys-root/lib64
	fi

	if [ $target = "powerpc64le-linux-gnu" ]; then
	  ln -s /opt/powerpc64le-linux-gnu/powerpc64le-linux-gnu/sys-root/lib64/ld64.so.2 /lib64/ld64.so.2
	  export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/x86_64-linux-gnu/x86_64-linux-gnu/sys-root/lib64:/opt/powerpc64le-linux-gnu/powerpc64le-linux-gnu/sys-root/lib64
	fi

	if [ $target = "arm-linux-gnueabihf" ]; then
	  ln -s /opt/arm-linux-gnueabihf/arm-linux-gnueabihf/sys-root/lib/ld-linux-armhf.so.3 /lib/ld-linux-armhf.so.3
	  export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/x86_64-linux-gnu/x86_64-linux-gnu/sys-root/lib64:/opt/arm-linux-gnueabihf/arm-linux-gnueabihf/sys-root/lib
	fi

	cd $WORKSPACE/srcdir/fontconfig-2.13.1/
	export CFLAGS=-I/workspace/destdir/include
	./configure --prefix=$prefix --host=$target
	make
	make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:aarch64, libc=:glibc),
    Linux(:powerpc64le, libc=:glibc),
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libfontconfig", :libfontconfig),
    ExecutableProduct(prefix, "fc-query", :fc_query),
    ExecutableProduct(prefix, "fc-list", :fc_list),
    ExecutableProduct(prefix, "fc-validate", :fc_validate),
    ExecutableProduct(prefix, "fc-cat", :fc_cat),
    ExecutableProduct(prefix, "fc-conflist", :fc_conflist),
    ExecutableProduct(prefix, "fc-match", :fc_match),
    ExecutableProduct(prefix, "fc-scan", :fc_scan),
    ExecutableProduct(prefix, "fc-pattern", :fc_pattern),
    ExecutableProduct(prefix, "fc-cache", :fc_cache)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/bjarthur/ExpatBuilder/releases/download/v2.2.6/build_Expat.v2.2.6.jl",
    "https://github.com/bjarthur/UtilLinuxBuilder/releases/download/v2.32/build_UtilLinux.v2.32.0.jl",
    "https://github.com/JuliaGraphics/FreeTypeBuilder/releases/download/v2.9.1-1/build_FreeType2.v2.9.1.jl",
    "https://github.com/bjarthur/GperfBuilder/releases/download/v3.1.0/build_GperfBuilder.v3.1.0.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: bjarthur/FontconfigBuilder

max_issues_repo_name: bjarthur/FontconfigBuilder

max_forks_repo_name: JuliaGraphics/FontconfigBuilder


text:
using JuMP, Compat.Test, ReSHOP

## Solve test problem with non-linear binary variables
 #
 #  min   100 * (x2 - (0.5 + x1) ^ 2) ^ 2 + (1 - x1) ^ 2
 #  s.t.  x1, x2 binary
 #
 #  The solution is (0, 0).
 ##

if !isdef(:solver_sbb); solver_sbb = ReSHOP.ReSHOPSolver("sbb"); end

@testset "example: jump_nonlinearbinary" begin
    m = Model(solver=solver_sbb)
    @variable(m, x[1:2], Bin)

    # Set some non-binary bounds on x1 and x2. These should be ignored.
    # The optimal solution if x is Int is (1, 2) which is allowed by these bounds
    setupperbound(x[1], 2)
    setupperbound(x[2], 2)

    @NLobjective(m, Min, 100*(x[2] - (0.5 + x[1])^2)^2 + (1 - x[1])^2)

    @test solve(m) == :Optimal
    @test getvalue(x)[:] == [0.0, 0.0]
    @test getobjectivevalue(m) == 7.25
end


text:
include("multexpe.jl")
using MultRoot
using MultExperimental

import PyPlot; using PyPlot

p = ptest53(0)


text:
"""
    generate_imports(src_name)

We import the appropriate `Artifacts` module based on Julia version; Julia 1.3+ contains
an `Artifacts` module in `Pkg`, while Julia 1.6+ has it as a builtin standard library.
"""
function generate_imports(src_name)
    # We lie a bit in the registry that JLL packages are usable on Julia 1.0-1.2.
    # This is to allow packages that might want to support Julia 1.0 to get the
    # benefits of a JLL package on 1.3 (requiring them to declare a dependence on
    # the JLL package in their Project.toml) but engage in heroic hacks to do
    # something other than actually use a JLL package on 1.0-1.2.  By allowing
    # this package to be installed (but not loaded) on 1.0-1.2, we enable users
    # to avoid splitting their package versions into pre-1.3 and post-1.3 branches
    # if they are willing to engage in the kinds of hoop-jumping they might need
    # to in order to install binaries in a JLL-compatible way on 1.0-1.2. One
    # example of this hoop-jumping being to express a dependency on this JLL
    # package, then import it within a `VERSION >= v"1.3"` conditional, and use
    # the deprecated `build.jl` mechanism to download the binaries through e.g.
    # `BinaryProvider.jl`.  This should work well for the simplest packages, and
    # require greater and greater heroics for more and more complex packages.
    if VERSION < v"1.3.0-rc4"
        return quote
            error("Unable to use $(src_name)_jll on Julia versions older than 1.3!")
        end
    elseif VERSION < v"1.6.0-DEV"
        # Use slow Pkg-based Artifacts
        return quote
            using Libdl, Pkg, Pkg.BinaryPlatforms, Pkg.Artifacts
            using Pkg.Artifacts: load_artifacts_toml, unpack_platform
            using Pkg.BinaryPlatforms: triplet, select_platform
            HostPlatform() = platform_key_abi()
        end
    else
        # Use fast stdlib-based Artifacts
        return quote
            using Libdl, Artifacts, Base.BinaryPlatforms
            using Artifacts: load_artifacts_toml, unpack_platform
            using Base.BinaryPlatforms: triplet, select_platform
        end
    end
end

"""
    generate_compiler_options(src_name)

Because JLL packages do not contain code that benefits much from compiler optimizations,
we disable them for a sizable boost in first load time.
"""
function generate_compiler_options(src_name)
    # Newer Julias have `@compiler_options` that can enable interpreted mode
    if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@compiler_options"))
        return quote
            Core.eval($(Symbol("$(src_name)_jll")), :(Base.Experimental.@compiler_options compile=min optimize=0 infer=false))
        end
    end

    # Older Julias only have `@optlevel`
    if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
        return quote
            @eval Base.Experimental.@optlevel 0
        end
    end

    # If none of these are available, TOUGH BEANS.
    return nothing
end

"""
    generate_toplevel_definitions(src_name)

This method generates the toplevel definitions common to all JLL packages, such as
`is_available()` and `find_artifact_dir()`, etc...
"""
function generate_toplevel_definitions(src_name, __source__)
    pkg_dir = dirname(String(__source__.file))
    return quote
        """
            is_available()

        Return whether the artifact is available for the current platform.
        """
        function is_available end

        PATH = ""
        LIBPATH = ""
        # We put these inter-JLL-package API values here so that they are always defined, even if there
        # is no underlying wrapper held within this JLL package.
        const PATH_list = String[]
        const LIBPATH_list = String[]

        function find_artifact_dir()
            # We determine at compile-time whether our JLL package has been dev'ed and overridden
            @static if isdir(joinpath(dirname($(pkg_dir)), "override"))
                return joinpath(dirname($(pkg_dir)), "override")
            else
                # We explicitly use `macrocall` here so that we can manually pass the `__source__`
                # argument, to avoid `@artifact_str` trying to lookup `Artifacts.toml` here.
                return $(Expr(:macrocall, Symbol("@artifact_str"), __source__, src_name))
            end
        end

        # We sub off to JLLWrappers' dev_jll, but avoid backedges
        dev_jll() = Base.invokelatest(JLLWrappers.dev_jll, $(src_name))
    end
end

"""
    generate_wrapper_load(src_name, pkg_uuid)

Because each platform could have completely disjoint products, we embed the information
for each in its own `.jl` file, named by triplet, and stored within the `src/wrappers`
directory.  This method generates the code to load `Artifacts.toml`, parse it for
artifacts, find the one that matches the host platform, then load the matching wrapper.
"""
function generate_wrapper_load(src_name, pkg_uuid, __source__)
    jll_name = "$(src_name)_jll"
    pkg_dir = dirname(String(__source__.file))

    function platform_parse_compat()
        if VERSION < v"1.6.0-DEV"
            return :(parse_wrapper_platform(x) = platform_key_abi(x))
        else
            return :(parse_wrapper_platform(x) = parse(Platform, x))
        end
    end

    return quote
        global best_wrapper
        # Load Artifacts.toml file and select best platform at compile-time, since this is
        # running at toplevel, and therefore will be run completely at compile-time.  We use
        # a `let` block here to avoid storing unnecessary data in our `.ji` files
        best_wrapper = let
            artifacts_toml = joinpath($(pkg_dir), "..", "Artifacts.toml")
            valid_wrappers = Dict{Platform,String}()
            artifacts = load_artifacts_toml(artifacts_toml; pkg_uuid=$(pkg_uuid))[$(src_name)]

            # Helper function to parse triplets for us
            $(platform_parse_compat())
            function make_wrapper_dict(dir, x)
                return Dict(
                    parse_wrapper_platform(basename(f)[1:end-3]) =>
                        joinpath(dir, "wrappers", f) for f in x
                )
            end

            # If it's a Dict, that means this is an AnyPlatform artifact, act accordingly:
            if isa(artifacts, Dict)
                joinpath($(pkg_dir), "wrappers", "any.jl")
            else
                # Otherwise, it's a Vector, and we must select the best platform
                # First, find all wrappers on-disk, parse their platforms, and match:
                wrapper_files = filter(x -> endswith(x, ".jl"), readdir(joinpath($(pkg_dir), "wrappers")))
                wrappers = make_wrapper_dict($(pkg_dir), wrapper_files)
                for e in artifacts
                    platform = unpack_platform(e, $(jll_name), artifacts_toml)

                    # Filter platforms based on what wrappers we've generated on-disk.
                    # Because the wrapper file naming strategy relies upon BB's assumptions of triplet
                    # parsing, it can be somewhat fragile.  So we have loaded all the wrapper filenames
                    # in, parsed them, and are now using `select_platform()` to match them, which is
                    # much more robust.  This step avoids a disconnect between what is recorded in the
                    # `Artifacts.toml` file and what wrappers are available on-disk.
                    wrapper_file = select_platform(wrappers, platform)
                    if wrapper_file !== nothing
                        valid_wrappers[platform] = wrapper_file
                    end
                end

                # From the available options, choose the best wrapper script
                select_platform(valid_wrappers)
            end
        end

        # Load in the wrapper, if it's not `nothing`!
        if best_wrapper === nothing
            @debug(string("Unable to load ", $(src_name), "; unsupported platform ", triplet(HostPlatform())))
            is_available() = false
        else
            Base.include($(Symbol("$(src_name)_jll")), best_wrapper)
            is_available() = true
        end
    end
end

macro generate_main_file_header(src_name)
    return excat(
        # Declare this module as interpreted
        generate_compiler_options(src_name),
        # import Artifacts module
        generate_imports(src_name),
    )
end

"""
    @generate_main_file(src_name, pkg_uuid)

Generate the main JLL file to Import `Artifacts`, set compiler options, perform forward
declarations for JLL-internal APIs, load the wrappers, etc...
"""
macro generate_main_file(src_name, pkg_uuid)
    return excat(
        # `is_available()` forward declaration, `PATH_list` and `LIBPATH_list` forward definitions,
        # `find_artifact_dir()` definition.
        generate_toplevel_definitions(src_name, __source__),
        # Select and load the best wrapper file
        generate_wrapper_load(src_name, pkg_uuid, __source__),
    )
end


text:
@testset "setup object" begin
    param = :complexConjug
    θ₀ = [10.0, -5.0, 5.0, 0.0, 3.0]
    P˟ = BSI.FitzhughDiffusion(param, θ₀...)
    obs = [[0.0, 0.0], [-0.3, -0.2], [0.4, 0.2]]
    tt = [0.0, 1.0, 1.5]
    P̃ = [BSI.FitzhughDiffusionAux(param, θ₀..., tt[1], obs[1], tt[2], obs[2]),
         BSI.FitzhughDiffusionAux(param, θ₀..., tt[2], obs[2], tt[3], obs[3])]
    setup = BSI.MCMCSetup(P˟, P̃, PartObs())

    @testset "initialisation" begin
        @test setup.P˟ == P˟
        @test setup.P̃ == P̃
        @test setup.blocking == BSI.NoBlocking()
        @test setup.blocking_params == ([], 0.1, BSI.NoChangePt())
        @test !any(values(setup.setup_completion))
        @test @suppress !BSI.check_if_complete(setup, [:obs])
        @test @suppress !BSI.check_if_complete(setup, [:imput])
        @test @suppress !BSI.check_if_complete(setup, [:tkern])
        @test @suppress !BSI.check_if_complete(setup, [:prior])
        @test @suppress !BSI.check_if_complete(setup, [:mcmc])
        @test @suppress !BSI.check_if_complete(setup, [:solv])
    end

    L = [1. 0.; 0. 1.]
    Σ = [0.5 0.0; 0.0 1.0]
    BSI.set_observations!(setup, [L, L], [Σ, 2*Σ], obs, tt)

    @testset "setting observations" begin
        @test setup.Ls == [L, L]
        @test setup.Σs == [Σ, 2*Σ]
        @test setup.obs == obs
        @test setup.obs_times == tt
        @test setup.fpt == [nothing, nothing]
        @test @suppress !BSI.check_if_complete(setup, [:imput])
        @test @suppress !BSI.check_if_complete(setup, [:tkern])
        @test @suppress !BSI.check_if_complete(setup, [:prior])
        @test @suppress !BSI.check_if_complete(setup, [:mcmc])
        @test @suppress !BSI.check_if_complete(setup, [:solv])
        @test @suppress BSI.check_if_complete(setup, [:obs])
    end

    dt = 0.01
    τ(t₀,T) = (x) ->  t₀ + (x-t₀) * (2-(x-t₀)/(T-t₀))
    BSI.set_imputation_grid!(setup, dt)

    @testset "setting imputation grid" begin
        @test setup.dt == dt
        @test setup.τ(tt[1], tt[2])(0.5*(tt[1]+tt[2])) == τ(tt[1], tt[2])(0.5*(tt[1]+tt[2]))
        @test @suppress !BSI.check_if_complete(setup, [:tkern])
        @test @suppress !BSI.check_if_complete(setup, [:prior])
        @test @suppress !BSI.check_if_complete(setup, [:mcmc])
        @test @suppress !BSI.check_if_complete(setup, [:solv])
        @test @suppress BSI.check_if_complete(setup, [:obs, :imput])
    end

    t_kernels = [BSI.RandomWalk([0.002, 0.1], [true, true]),
                 BSI.RandomWalk([0.2, 1.0], 2)]
    ρ = 0.5
    param_updt = true
    updt_coord_true = (Val((true,true,false,false,false)),
                       Val((false,true,true,false,false)))
    updt_coord = ((1,2),(2,3))
    updt_type=(BSI.MetropolisHastingsUpdt(),
               BSI.ConjugateUpdt())
    set_transition_kernels!(setup, t_kernels, ρ, param_updt, updt_coord,
                            updt_type)
    @testset "setting transition kernels" begin
        @test setup.t_kernel == t_kernels
        @test setup.ρ == ρ
        @test setup.param_updt == param_updt
        @test setup.updt_coord == updt_coord_true
        @test setup.updt_type == updt_type
        @test !BSI.check_if_adapt(setup.adaptive_prop)
        @test @suppress !BSI.check_if_complete(setup, [:prior])
        @test @suppress !BSI.check_if_complete(setup, [:mcmc])
        @test @suppress !BSI.check_if_complete(setup, [:solv])
        @test @suppress BSI.check_if_complete(setup, [:obs, :imput, :tkern])
    end

    @testset "reformat updt_coord object" begin
        @test BSI.reformat_updt_coord((1,3),[1,2,3]) == (Val((true,false,true)),)
        @test BSI.reformat_updt_coord(nothing,[1,2,3]) == (Val((true,)),)
        @test BSI.reformat_updt_coord([(1,2),(2,3)],[1,2,3]) == (Val{(true, true, false)}(), Val{(false, true, true)}())
        @test_throws AssertionError BSI.reformat_updt_coord([4],[1,2,3])
    end

    priors = BSI.Priors((BSI.ImproperPrior(), BSI.ImproperPrior()))
    x0_prior = BSI.KnownStartingPt(obs[1])
    set_priors!(setup, priors, x0_prior)
    @testset "setting priors" begin
        @test setup.priors == priors
        @test setup.x0_prior == x0_prior
        @test @suppress !BSI.check_if_complete(setup, [:mcmc])
        @test @suppress !BSI.check_if_complete(setup, [:solv])
        @test @suppress BSI.check_if_complete(setup, [:obs, :imput, :tkern, :prior])
    end

    num_mcmc_steps = 100
    BSI.set_mcmc_params!(setup, num_mcmc_steps)
    @testset "setting  mcmc parameters" begin
        @test setup.num_mcmc_steps == num_mcmc_steps
        @test isnan(setup.save_iter)
        @test isnan(setup.verb_iter)
        @test setup.skip_for_save == 1
        @test setup.warm_up == 0
        @test @suppress !BSI.check_if_complete(setup, [:solv])
        @test @suppress BSI.check_if_complete(setup, [:obs, :imput, :tkern, :prior,
                                                  :mcmc])
    end

    BSI.set_solver!(setup)
    @testset "setting solver" begin
        @test setup.solver == BSI.Ralston3()
        @test setup.change_pt == BSI.NoChangePt()
        @test @suppress BSI.check_if_complete(setup)
    end

    @testset "determining data type" begin
        @test BSI.determine_data_type(setup) == (SArray{Tuple{2},Float64,1,2},
                                             Float64)
    end

    BSI.prepare_containers!(setup)
    @testset "setting internal containers" begin
        @test setup.Wnr == Wiener{Float64}()
        @test eltype(setup.WW) == SamplePath{Float64}
        @test length(setup.WW) == 2
        @test eltype(setup.XX) == SamplePath{SArray{Tuple{2},Float64,1,2}}
        @test length(setup.XX) == 2
        @test setup.Ls == [SMatrix{2,2}(L), SMatrix{2,2}(L)]
        @test setup.Σs == [SMatrix{2,2}(Σ), SMatrix{2,2}(2*Σ)]
        @test setup.obs == map(x->SVector{2}(x), obs)
    end

    BSI.initialise!(Float64, setup)
    @testset "initialisation of proposal law" begin
        @test length(setup.P) == 2
        @test typeof(setup.P[1].Pt) == typeof(setup.P̃[1])
        @test setup.P[1].Target == setup.P˟
    end
end


text:
struct DiscontinuousGalerkin1DGrid{N, C}
                node_coordinates :: N
    element_to_node_connectivity :: C
end

"""
Generate simple equidistant grid with N elements.
"""
function DiscontinuousGalerkin1DGrid(domain, N)
    xₗ, xᵣ = domain

    node_coordinates = [xₗ + (xᵣ - xₗ) * (i - 1) / N for i in 1:N+1]

    # read element to node connectivity
    element_to_node_connectivity = zeros(N, 2)
    for n in 1:N
        element_to_node_connectivity[n, 1] = n
        element_to_node_connectivity[n, 2] = n+1
    end

    return DiscontinuousGalerkin1DGrid{typeof(node_coordinates), typeof(element_to_node_connectivity)}(
        node_coordinates, element_to_node_connectivity)
end


text:
# This file was generated, do not modify it. # hide
colnames = names(weather)

filter_by_name(name, cols) =
    filter(cn -> occursin(name, String(cn)), cols)

wind   = weather[:, filter_by_name("windspeed", colnames)]
temp   = weather[:, filter_by_name("temperature", colnames)]
raddir = weather[:, filter_by_name("radiation_direct", colnames)]
raddif = weather[:, filter_by_name("radiation_diffuse", colnames)];


text:
using FFTW
using LinearAlgebra
using DSP
using Plots
using Statistics
using ProgressMeter
include("pcfm.jl")

function autocorr(x)
  return conv(x, conj.(x[end:-1:1]))
end

function ∇J(B, x, u, l)
  """
  ∇J(B,Bb,x,m,u,a)
  Compute the frequency template error (FTE) and its associated gradient

  # Arguments
  - `B::Array`: M X M Orthogonal basis
  - `x::Vector`: M X 1 phase code vector
  - `u::Vector`: 2M-1 X 1 Frequency template
  - `l::Integer`: Norm to use for error calculation
  """
  # Oversampled phase code length
  m = size(B, 1)
  # Zero-pad basis matrix to length 2M-1 (for FFT)
  Bb = vcat(B, zeros(m - 1, trunc(Int, size(B, 2))))
  # PCFM representation of the input phase code vector
  s = exp.(im .* B * x)
  # Pad the waveform to length 2M-1
  sb = vcat(s, zeros(m - 1, size(s, 2)))
  # Compute the (normalized) PSD of the PCFM waveform
  sbf = fftshift(fft(sb))
  sbf = sbf ./ maximum(abs.(sbf))
  # FTE calculation
  J = norm(abs.(sbf) .^ 2 .- u, l)
  # Return the error and gradient
  return (J,
    2 / (J) .* transpose(Bb) * imag.(conj.(sb) .* ifft(ifftshift((abs.(sbf) .^ 2 .- u) .* sbf))))
end

function ∇logJ(B, x, u, a, l)
  """
  ∇logJ(B,Bb,x,m,u,a)
  Compute the log frequency template error (log-FTE) and its associated gradient

  # Arguments
  - `B::Array`: M X M Orthogonal basis
  - `x::Vector`: M X 1 phase code vector
  - `u::Vector`: 2M-1 X 1 Frequency template
  - `a::Integer`: Log base for error computation
  - `l::Integer`: Norm to use for error calculation

  """
  # Oversampled phase code length
  m = size(B, 1)
  # Zero-pad basis matrix to length 2M-1 (for FFT)
  Bb = vcat(B, zeros(m - 1, size(B, 2)))
  # PCFM representation of the input phase code vector
  s = exp.(im .* B * x)
  # Pad the waveform to length 2M-1
  sb = vcat(s, zeros(m - 1, size(s, 2)))
  # Compute the (normalized) PSD of the PCFM waveform
  sbf = fftshift(fft(sb))
  sbf = sbf ./ maximum(abs.(sbf), dims = 1)
  # log-FTE calculation
  J = norm(log.(a, abs.(sbf) .^ 2) .- log.(a, u), l)
  # Return the error and gradient
  return (J,
    (2 / (log(a) * J)) .* transpose(Bb) * imag.(conj.(sb) .* ifft(ifftshift((log.(a, abs.(sbf) .^ 2) .- log.(a, u)) .* sbf))))
end

function profm(u, iter)
  """
  profm(u,iter)

  Iteratively optimize the PSD using alternating projections as described in the
  PRO-FM paper.
  """
  pk = exp.(im .* angle.(ifft(ifftshift(u))))
  for ii = 1:iter
    rk = ifft(ifftshift(abs.(u) .* exp.(im .* angle.(fftshift(fft(pk))))))
    pk = exp.(im .* angle.(rk))
  end
  return pk
end

function optimize(u, nWaveforms, k; a = 10, tol = 1e-5, maxIter = 1000, savePlots = false, filename = "anim_fps60.gif")
  """
  optimize(u,a,tol,maxIter)

  # Arguments
  - `u::Matrix`: 2M-1 X 1 Frequency template
  - `a::Int`: Log base for log-FTE (ignored if regular FTE is used)
  - `tol::Float64`: Tolerance for early stopping
  - `maxIter::Integer`: Maximum number of iterations
  """
  # TODO: Allow the user to decide which optimization method to use

  #Calculate m from u.
  m = trunc(Int, (length(u) + 1) / 2)
  # Get a randomly initialized phase change vector and phase shaping basis
  # functions from PCFM generator
  (_, x, B) = pcfm(m, k, nWaveforms)
  # Gradient descent Parameters
  μ = 0.5
  β = 0.9
  # Store the error at each iteration
  Jvec = ones(maxIter - 1, 1)
  pkOld = 0
  sbf = zeros(length(u), nWaveforms)
  anim = @animate for ii = 1:maxIter
    # Heavy-ball gradient descent
    (J, ∇) = ∇J(B, x, u, 2)
    # (J, ∇) = ∇logJ(B, x, u, a, 2)
    Jvec[ii:end] .= J
    # Update the search direction
    if ii == 1
      pk = ∇
    else
      pk = ∇ .+ β .* pkOld
    end
    # Update x
    x -= μ .* pk
    # Stopping condition
    if all(abs.(pk .- pkOld) .< tol)
      break
    end
    pkOld = pk
    # Plots
    # Compute and plot the PSD
    s = exp.(im .* B * x)
    sb = vcat(s, zeros(m - 1, nWaveforms))
    sbf = fftshift(fft(sb))
    sbf = mean(sbf, dims = 2)
    sbf = sbf ./ maximum(abs.(sbf), dims = 1)
    p1 = plot(10 * log10.(abs.(sbf) .^ 2), ylim = (-50, 0),
      xlabel = "Sample Index", label = "Actual PSD")
    plot!(10 * log10.(u), ylim = (-50, 0), label = "PSD Template")
    # Compute and plot the autocorrelation
    corr = abs.(autocorr(s)) ./ maximum(abs.(autocorr(s)))
    p2 = plot(10 * log10.(corr), ylim = (-30, 0), xlabel = "Sample Index",
      ylabel = "Magnitude (dB)", label = "Autocorrelation")
    # Compute and plot the current error
    p3 = plot(Jvec, ylim = (0, maximum(Jvec)), xlabel = "Iteration Number", ylabel = "Error", label="Error")
    plot(p1, p2, p3, layout = (3, 1), linewidth = 2)
  end
  if savePlots
    gif(anim, filename, fps = 60)
  end
  s = exp.(im .* B * x)
  return (x, s)

end


text:
using MRIReco, FFTW, LinearAlgebra, BenchmarkTools

FFTW.set_num_threads(1);BLAS.set_num_threads(1)

N = 256
numCoils = 8
coilsens = birdcageSensitivity(N, numCoils, 1.5)
cmap = 1im*quadraticFieldmap(N,N,125*2pi)

f = ISMRMRDFile(joinpath(@__DIR__(),"simparallel.h5"))
raw = RawAcquisitionData(f)
acqData = AcquisitionData(raw)

# reco parameters
params = Dict{Symbol, Any}()
params[:reco] = "multiCoil"
params[:reconSize] = (N,N)
params[:regularization] = "L2"
params[:λ] = 1.e-3
params[:iterations] = 40
params[:solver] = "cgnr"
params[:senseMaps] = coilsens
#params[:correctionMap] = cmap

t = @belapsed global img = reconstruction(acqData, params)
filename = joinpath(@__DIR__(),"reco.png")
exportImage(filename, abs.(img[:,:,1,1,1]) )

println(t)


text:
module CachesimOut

export init_cachesimout, take_available_addr_range, add_cachesim_array, add_cachesim_tmp_array, remove_all_tmp_arrays,
        cachesim_load_range, cachesim_store_range, cachesim_load, cachesim_store, build_cache_level, build_cache,
        build_cache_auto, analyze_cache_with_grid
export CachesimArray

use_lib = true;
if use_lib
    include("Jpycachesim.jl");
    using .Jpycachesim

else

end

addr_offset = 0;
line_size = 0;
geofac_set = false;
arrays = [];
tmparrays = [];
tmpflags = [];
N=0; # number of global nodes
Np=0;# number of elemental nodes
Nn=0;# number of global dofs
Nd=0;# number of elemental dofs
Nv=0;# number of global vertices
Nel = 0;# number of elements

active_cache = nothing; # keep track of the active cache

# A special nonexistant array that also has an identifier, fake address range, size
mutable struct CachesimArray
    id::Int             # identifying unique number
    addr::Array{Int,1}  # fake address range
    size::Array{Int,1}  # size of the array(in elements)
    elsize::Int         # size of the elements in the array(in bytes)
    bytesize::Int       # total number of bytes in array

    CachesimArray(i, sz, elsz) = new(
        i,
        take_available_addr_range(sz, elsz),
        size_tuple_to_array(sz),
        elsz
    )
end

function init_cachesimout(n, refel, nel, dof, vars)
    global addr_offset = 0;
    global geofac_set = false;
    global arrays = [];
    global tmparrays = [];
    global tmpflags = [];

    global N = n; # number of global nodes
    global Nn = N*dof; # number of global dof
    global Nd = refel.Np*dof; # elemental dof
    global Np = refel.Np; # number of elemental nodes
    global Nel = nel; # number of elements
    #=
    The arrays are indexed as such:
    1 A
    2 b
    3 u
    4 Ak
    5 bk
    6 Q
    7 Q1
    8 Q2
    9 Q3
    10 D1
    11 D2
    12 D3
    13 Jx
    14 Jy
    15 Jz
    16 detJ
    17 w
    18+ variable values, then everything else
    =#
    push!(arrays, CachesimArray(1, [Nel*Nd*Nd, 3], 8));
    push!(arrays, CachesimArray(2, [Nn], 8));
    push!(arrays, CachesimArray(3, [Nn], 8));
    push!(arrays, CachesimArray(4, [Nd, Nd], 8));
    push!(arrays, CachesimArray(5, [Nd], 8));
    push!(arrays, CachesimArray(6, size(refel.Q), 8));
    push!(arrays, CachesimArray(7, size(refel.Q), 8));
    push!(arrays, CachesimArray(8, size(refel.Q), 8));
    push!(arrays, CachesimArray(9, size(refel.Q), 8));
    push!(arrays, CachesimArray(10, size(refel.Q), 8));
    push!(arrays, CachesimArray(11, size(refel.Q), 8));
    push!(arrays, CachesimArray(12, size(refel.Q), 8));
    push!(arrays, CachesimArray(13, size(refel.Q), 8));
    push!(arrays, CachesimArray(14, size(refel.Q), 8));
    push!(arrays, CachesimArray(15, size(refel.Q), 8));
    push!(arrays, CachesimArray(16, [Nd], 8));
    push!(arrays, CachesimArray(17, [Nd], 8));

    for i=1:length(vars)
        push!(arrays, CachesimArray(17+i, size(vars[i].values), 8));
    end

    if use_lib
        global active_cache =  pcs_get_cachesim_from_file("src/cachesim/cachedef"); # Just to test things. Should be set later by user.
    else
        println("Not using the cachesim library. Change use_lib in CachesimOut.jl to use.");
        println("Instead, writing access sequence to cachesim_output.out");
        println("Analyze with cachesim_read.py");
        global output = open("cachesim/cachesim_output.out","w");
    end

end

function build_cache_level(level, sets, ways, cl_size, policy)
    if use_lib
        label = "L"*string(level);
        return Cache(label, sets, ways, cl_size, policy);
    else
        println("Can't build cache without library.");
    end
end

function build_cache(levels)
    if use_lib
        global active_cache = pcs_build_cache(levels);
        global line_size = Int(levels[1].cl_size);
        return active_cache;
    else
        println("Can't build cache without library.");
    end
end

# Build a single-level cache with enough capacity to hold el_num*np values
function build_cache_auto(el_num, grid; bytes_per_val=8, bytes_per_line=64, assoc=1)
    Np = size(grid.loc2glb,1);
    nel = size(grid.loc2glb,2);
    dim = size(grid.allnodes,1);

    lines = Int(ceil((el_num*Np*bytes_per_val)/bytes_per_line));
    lines = max(4, lines);

    # assoc=0 -> fully associative
    # assoc=1 -> direct mapped
    if assoc == 0
        assoc = lines;
    end
    assoc = min(assoc, lines);

    l1 = build_cache_level(1, lines-assoc+1, assoc, bytes_per_line, "LRU");
    return build_cache([l1]);
end

function set_active_cache(c)
    global active_cache = c;
end

####################################################################################

function add_cachesim_array(sz, elsz)
    id = length(arrays)+1;
    push!(arrays, CachesimArray(id, sz, elsz));
    return id;
end

function add_cachesim_array_aligned(sz, elsz)
    # find the offset to align with the next cache line
    offset = mod(addr_offset, line_size);
    if offset > 0
        add_cachesim_array(line_size - offset,1);
    end
    return add_cachesim_array(sz, elsz);
end

function add_cachesim_tmp_array()
    # search the flags for an open tmp array
    id = -1;
    for i=1:length(tmpflags)
        if tmpflags[i] == false
            stupidjulia = id;
            id = tmparrays[i];
            tmpflags[i] = true;
            break;
        end
    end
    # get a new array if needed
    if id < 0
        id = add_cachesim_array([Np],8);
        push!(tmparrays, id);
        push!(tmpflags, true);
    end
    return id;
end

function remove_tmp_array(id)
    for i=1:length(tmparrays)
        if tmparrays[i] == id
            tmpflags[i] = false;
            break;
        end
    end
end

function remove_all_tmp_arrays()
    for i=1:length(tmparrays)
        tmpflags[i] = false;
    end
end

function take_available_addr_range(sz, elsz)
    range = [0,0];
    totalsz = sz[1];
    for i=2:length(sz)
        totalsz = totalsz*sz[i];
    end
    bytesz = totalsz*elsz;

    global addr_offset;
    range[1] = addr_offset;
    range[2] = range[1] + bytesz - 1;
    addr_offset += bytesz;

    return range;
end

function size_tuple_to_array(sz)
    if typeof(sz) <: Array
        return sz;
    elseif typeof(sz) <: Tuple
        ar = [];
        for i=1:length(sz)
            push!(ar,sz[i]);
        end
        return ar;
    elseif typeof(sz) <: Number
        return [sz];
    end
end

function ind_to_addr(id, row, col=0)
    start = arrays[id].addr[1];
    elsize = arrays[id].elsize;
    if col == 0
        return start + (row-1)*elsize;
    else
        colsize = arrays[id].size[1];
        return start + (col-1)*colsize + (row-1)*elsize;
    end
end

function finalize()
    if use_lib
        pcs_print_stats();
    else
        close(output);
    end
end

##########################################################################################
# Interact with cache - load/store
#
function cachesim_load_range(id, rows = [], cols = [])
    if length(rows) == 0
        rows = 1:arrays[id].size[1];
    end
    if length(cols) == 0 && length(arrays[id].size) > 1
        cols = 1:arrays[id].size[2];
    end
    if length(cols) > 0
        for j=1:length(cols)
            for i=1:length(rows)
                cachesim_load(id, rows[i], cols[j]);
            end
        end
    else
        for i=1:length(rows)
            cachesim_load(id, rows[i]);
        end
    end
end
function cachesim_store_range(id, rows = [], cols = [])
    if length(rows) == 0
        rows = 1:arrays[id].size[1];
    end
    if length(cols) == 0 && length(arrays[id].size) > 1
        cols = 1:arrays[id].size[2];
    end
    if length(cols) > 0
        for j=1:length(cols)
            for i=1:length(rows)
                cachesim_store(id, rows[i], cols[j]);
            end
        end
    else
        for i=1:length(rows)
            cachesim_store(id, rows[i]);
        end
    end
end

function cachesim_load(id, row, col=0)
    addr = ind_to_addr(id, row, col);
    if use_lib
        pcs_load(addr, 8, active_cache);
    else
        write(output,Int8(0));
        write(output,addr);
    end
end

function cachesim_store(id, row, col=0)
    addr = ind_to_addr(id, row, col);
    if use_lib
        pcs_store(addr, 8, active_cache);
    else
        write(output,Int8(1));
        write(output,addr);
    end
end

function cachesim_load_sparse(id, range=-1, rowoffset=0)
    a = arrays[id];
    alen = a.addr[2] - a.addr[1] + 1;
    if range == -1
        range = 1:alen;
    end

    for i=1:length(range)
        ad = a.addr[1] + rowoffset + (range[i]-1)*a.elsize;
        if use_lib
            pcs_load(ad, 8, active_cache);
        else
            write(output,Int8(0));
            write(output,ad);
        end
    end

end

function cachesim_store_sparse(id, range=-1, rowoffset=0)
    a = arrays[id];
    alen = a.addr[2] - a.addr[1] + 1;
    if range == -1
        range = 1:alen;
    end

    for i=1:length(range)
        ad = a.addr[1] + rowoffset + (range[i]-1)*a.elsize;
        if use_lib
            pcs_store(ad, 8, active_cache);
        else
            write(output,Int8(1));
            write(output,ad);
        end
    end
end
#########################################################################################

#########################################

function analyze_cache_with_grid(cache, grid, el_order, dofs)
    set_active_cache(cache);
    cshobj = unsafe_load(cache);

    bytes_per_val = 8;
    bytes_per_line = Int(cshobj.cl_size);
    bytes_per_node = bytes_per_val * dofs;

    nodes_per_line = bytes_per_line / bytes_per_node;

    np = size(grid.loc2glb,1);
    dim = size(grid.allnodes,1);
    nnodes = size(grid.allnodes,2);
    nel = length(el_order);

    bigarray = add_cachesim_array_aligned([nnodes*dofs],bytes_per_val);

    coveredlines = zeros(Int,np*2);
    lineBins = zeros(Int,np*2);
    ave_lines = 0;

    # Loop ever elements
    for ei=el_order
        # Find number of lines touched by each element and average
        lines = 0;
        # Loop over nodes
        for ni=1:np
            tmp = lines;
            already1 = false;
            already2 = false;
            gind = grid.loc2glb[ni,ei];
            startaddr = (gind-1)*bytes_per_node;
            endaddr = startaddr + bytes_per_node-1;
            lin1 = floor(startaddr/bytes_per_line);
            lin2 = floor(endaddr/bytes_per_line);
            coveredlines[ni*2-1] = lin1;
            coveredlines[ni*2] = lin2;
            # Check if line is already covered
            for li=1:((ni-1)*2)
                if coveredlines[li] == lin1
                    tmp = already1;
                    already1 = true;
                end
                if coveredlines[li] == lin2
                    tmp = already2;
                    already2 = true;
                end
                if already1 && already2
                    break;
                end
            end
            if !already1
                lines = lines+1;
            end
            if !already2
                lines = lines+1;
            end

            # do a simple load of this node
            cachesim_load(bigarray, gind);
        end
        # add to totals
        lineBins[lines] += 1;

    end

    # average
    for li=1:(np*2)
        ave_lines = ave_lines + lineBins[li]*li;
    end
    ave_lines /= nel;

    #println("cache # = "*string(cache));

    println("[lines] - number of elements");
    for i=1:length(lineBins)
        if lineBins[i] > 0
            println("["*string(i)*"] - "*string(lineBins[i]));
        end
    end

    println("ave lines/element = "*string(ave_lines));
    # print some stats
    (misses, hits) = pcs_get_stats(cache);

    println("simplified stats:");
    pcs_print_stats_simple(cache);
end

end# module


text:
@testset "Current step" begin
    N = 10
    lc = LinearColorant(1, N, [crayon"light_red", crayon"yellow", crayon"light_green"])
    c = CurrentStep(N, lc=lc, header="Current step")
    for (i,color) in enumerate([crayon"light_red",
                                crayon"light_red",
                                crayon"light_red",
                                crayon"yellow",
                                crayon"yellow",
                                crayon"yellow",
                                crayon"yellow",
                                crayon"light_green",
                                crayon"light_green",
                                crayon"light_green"])
        si = i == 10 ? "10" : " $i"
        @test format(c.fmt, c(i)...) == "[$(color)$(si)/$(N)$(crayon"reset")]     "
    end
end


text:
using JuMP
using CPLEX
include("parser.jl")

function AutoBenders(MyFileName::String, nn = 0, bendersbool = false)
  n, f, c, w, d = read_instance(MyFileName)
  if nn != 0
    n = nn
    d = n÷2
    f = Array{Int64,1}(zeros(n))
    c = Array{Int64,1}(zeros(n))
    f[1], c[1] = 7, 8
    for i in 2:n
      f[i] = (f[i-1]*f[1])%159
      c[i] = (c[i-1]*c[1])%61
    end
  end
  if bendersbool
    w = 1
  end

  # Create the model
  m = Model(CPLEX.Optimizer)

  ## Variables
  @variable(m, y[1:n-1], Bin)
  @variable(m, x[1:n-1] >= 0)

  ## Constraints
  @constraint(m, sum(x[i] for i in 1:n-1) == d)
  @constraint(m, [i in 1:n-1], x[i] <= y[i])
  @constraint(m, y[1] >= y[2])
  @constraint(m, y[1] >= y[3])

  ## Objective
  @objective(m, Min, sum(f[i]*y[i] + c[i]*x[i] for i in 1:n-1))

  #resolution
  set_optimizer_attribute(m, "CPXPARAM_Benders_Strategy",3)
  optimize!(m)
end

max_stars_repo_name: boumpteryx/MLA

max_issues_repo_name: boumpteryx/MLA

max_forks_repo_name: boumpteryx/MLA


text:
#Load colormap
using Colors
if (false)
    include("cubeYF.jl")
else
    include("viridis.jl")
end

#Load specific backend
using LaTeXStrings
include("backends/plotsGRBackend.jl")
include("backends/plotsMakieBackend.jl")
include("backends/plotlyBackend.jl")

function isaPlot3d()
    runPlotlyBackend()
end

function isaPlot3d(z::AMFMmodel, t::Vector{Float64}; backend="PlotsGR", FreqUnits="rad/s")
    if backend=="Makie"
        isaPlot3d_Makie(z,t )
    elseif backend=="Plotly"
        dict = Dict()
        components = []
        for i in 1:length(z.comps)
            componentDict = Dict()
            y = -z.comps[i].ω.(t)
            if FreqUnits == "Hz" y = y .* (1/(2π)) end
            colorMap = cmap[max.(min.(round.(Int, abs.(z.comps[i].a.(t)) .* 256 / 1), 256), 1)]
            colorMap = map(color -> "rgb($(red(color)),$(green(color)),$(blue(color)))", colorMap)
            componentDict["x"] = -t .* -1
            componentDict["y"] = y .* -1
            componentDict["z"] = real.(z.comps[i](t))
            componentDict["colorMap"] = colorMap
            push!(components, componentDict)
        end
        dict["components"] = components
        dict["freqUnits"] = FreqUnits
        runPlotlyBackend(dict)
    else
        isaPlot3d_PlotsGR(z,t, FreqUnits=FreqUnits)
    end
end

function isaPlot3d(z::AMFMmodel, t::StepRangeLen; backend="PlotsGR", FreqUnits="rad/s")
    if backend=="Makie"
        isaPlot3d_Makie(z,collect(t)  )
    elseif backend=="Plotly"
        runPlotlyBackend() # use runPlotlyBackend(dict) to specify data
    else
        isaPlot3d_PlotsGR(z,collect(t), FreqUnits=FreqUnits)
    end
end

#Construction of an isaPlot3d from a component set
function isaPlot3d(S::Array{Tuple{Function,Function,Real},1}, t::Vector{Float64}; backend="PlotsGR", FreqUnits="rad/s")
    if backend=="Makie"
        isaPlot3d_Makie(AMFMmodel(S),t  )
    elseif backend=="Plotly"
        runPlotlyBackend() # use runPlotlyBackend(dict) to specify data
    else
        isaPlot3d_PlotsGR(AMFMmodel(S),t, FreqUnits=FreqUnits)
    end
end
function isaPlot3d(S::Array{Tuple{Function,Function,Real},1}, t::StepRangeLen; backend="PlotsGR", FreqUnits="rad/s")
    if backend=="Makie"
        isaPlot3d_Makie(AMFMmodel(S),collect(t)  )
    elseif backend=="Plotly"
        runPlotlyBackend() # use runPlotlyBackend(dict) to specify data
    else
        isaPlot3d_PlotsGR(AMFMmodel(S),collect(t), FreqUnits=FreqUnits)
    end
end

#Construction of an isaPlot3d from an AMFMcomp
function isaPlot3d(ψ::AMFMcomp, t::Vector{Float64}; backend="PlotsGR", FreqUnits="rad/s")
    if backend=="Makie"
        isaPlot3d_Makie(AMFMmodel([ψ]), t )
    elseif backend=="Plotly"
        runPlotlyBackend() # use runPlotlyBackend(dict) to specify data
    else
        isaPlot3d_PlotsGR(AMFMmodel([ψ]), t, FreqUnits=FreqUnits)
    end
end
function isaPlot3d(ψ::AMFMcomp, t::StepRangeLen; backend="PlotsGR", FreqUnits="rad/s")
    if backend=="Makie"
        isaPlot3d_Makie(AMFMmodel([ψ]),collect(t)  )
    elseif backend=="Plotly"
        runPlotlyBackend() # use runPlotlyBackend(dict) to specify data )
    else
        isaPlot3d_PlotsGR(AMFMmodel([ψ]),collect(t), FreqUnits=FreqUnits)
    end
end

#Construction of an isaPlot3d from an canonical triplet
function isaPlot3d(C::Tuple{Function,Function,Real}, t::Vector{Float64}; backend="PlotsGR", FreqUnits="rad/s")
    if backend=="Makie"
        isaPlot3d_Makie(AMFMmodel([AMFMcomp(C)]), t  )
    elseif backend=="Plotly"
        runPlotlyBackend() # use runPlotlyBackend(dict) to specify data
    else
        isaPlot3d_PlotsGR(AMFMmodel([AMFMcomp(C)]), t, FreqUnits=FreqUnits)
    end
end
function isaPlot3d(C::Tuple{Function,Function,Real}, t::StepRangeLen; backend="PlotsGR", FreqUnits="rad/s")
    if backend=="Makie"
        isaPlot3d_Makie(AMFMmodel([AMFMcomp(C)]),collect(t)  )
    elseif backend=="Plotly"
        runPlotlyBackend() # use runPlotlyBackend(dict) to specify data
    else
        isaPlot3d_PlotsGR(AMFMmodel([AMFMcomp(C)]),collect(t), FreqUnits=FreqUnits)
    end
end

#Construction of an isaPlot3d from a numerical component
function isaPlot3d(𝚿ₖ::Array{AMFMcompN,1}; backend="PlotsGR")
    if backend=="Makie"
        isaPlot3d_Makie(𝚿ₖ)
    else
        isaPlot3d_PlotsGR(𝚿ₖ)
    end
end
function isaPlot3d(𝚿::AMFMcompN; backend="PlotsGR", FreqUnits="rad/s")
    if backend=="Makie"
        isaPlot3d_Makie([𝚿])
    else
        isaPlot3d_PlotsGR([𝚿])
    end
end
function isaPlot3d(Ψ::Vector{ComplexF64}, t::Vector{Float64}; backend="PlotsGR", FreqUnits="rad/s")
    if backend=="Makie"
        isaPlot3d_Makie([AMFMdemod(Ψ,t)])
    else
        isaPlot3d_PlotsGR([AMFMdemod(Ψ,t)])
    end
end
function isaPlot3d(Ψ::Vector{ComplexF64}, t::StepRangeLen; backend="PlotsGR", FreqUnits="rad/s")
    if backend=="Makie"
        isaPlot3d_Makie([AMFMdemod(Ψ,collect(t))])
    else
        isaPlot3d_PlotsGR([AMFMdemod(Ψ,collect(t))],)
    end
end

max_stars_repo_name: kricheso/ISA

max_issues_repo_name: kricheso/ISA

max_forks_repo_name: kricheso/ISA


text:
@testset "OffsetArrays" begin
    @testset "centered" begin
        if isdefined(OffsetArrays, :centered) # OffsetArrays v1.9
            check_range(r, f, l) = (@test first(r) == f; @test last(r) == l)
            check_range_axes(r, f, l) = check_range(Base.axes(r)[1], f, l)

            check_range(Base.axes(OffsetArrays.centered(1:3))[1], -1, 1)
            a = AxisArray(rand(3, 3), Axis{:y}(0.1:0.1:0.3), Axis{:x}(1:3))

            ca = OffsetArrays.centered(a)
            axs = Base.axes(ca)
            check_range(axs[1], -1, 1)
            check_range(axs[2], -1, 1)
            @test ca[OffsetArrays.center(ca)...] == a[OffsetArrays.center(a)...]

            axs = AxisArrays.axes(ca)
            check_range(axs[1].val, 0.1, 0.3)
            check_range(axs[2].val, 1, 3)
            check_range_axes(axs[1].val, -1, 1)
            check_range_axes(axs[1].val, -1, 1)
        end
    end
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

dt = 0.01
t_start = 0.0
t_stop = 8.0

pathToFMU = joinpath(dirname(@__FILE__), "../model/SpringFrictionPendulum1D.fmu")

myFMU = fmiLoad(pathToFMU)

@test fmiInstantiate!(myFMU; loggingOn=false) != 0

@test fmiSetupExperiment(myFMU, 0.0) == 0
@test fmiEnterInitializationMode(myFMU) == 0
@test fmiExitInitializationMode(myFMU) == 0

solution = fmiSimulateME(myFMU, dt, t_start, t_stop)
@test length(solution.u) > 0
@test length(solution.t) > 0
# ToDo: Time series comparision

fmiUnload(myFMU)


text:
# imports
include("../iron_line_profile.jl")
using Measures
gr()

"""
Automatically determines the labels for various lines to go into the legend
"""
function label(value, name, units)
    scale = floor(log(10, value)) - 1
    if scale > 4
        scale = Int32(scale)
        value = "\$10^{$scale}\$"
    else
        value = "\$$value\$"
    end
    label = "$name = $value $(units) "
end

# values to change
# vars = 5:40:85 # obs angles
# vars=[5,15]
vars = [0.0,0.25,0.50,0.75,0.998] # spins
# vars = [10, 30, 60, 75]

# legend labels, colours and styles
labels = permutedims(label.(vars, "Spin",""))
styles = [:solid :dash :dot :dashdot :dashdotdot]
colors = [:black :blue :red :green :purple]

# values to be plotted
bins_vals = []
x_vals_vals = []

# generating data for different spins
for var in vars
    fvar = Float64(var)
    x_vals, bins = iron_line_profile(
                                    tolerance=1e-12,
                                    size_multiplier=3,
                                    fov=12,
                                    dtmax=5,
                                    spin=fvar,
                                    output="data",
                                    normalised=false,
                                    obs_angle=40
                                    )
    push!(x_vals_vals, x_vals)
    push!(bins_vals, bins)
end

# adding a line to represent the lab frame 6.4 keV emission line
verticle_line_x = fill(6.4, 100)
verticle_line_y = LinRange(0, 1.37e55, 100)

# plotting
plt = plot(
        x_vals_vals[1],
        bins_vals,
        framestyle=:box,
        grid=false,
        xlabel="Energy (keV)",
        ylabel="Flux (Arbitrary Units)",
        label=labels,
        linestyles=styles,
        linecolor=colors,
        legend=:topleft
        )
# vertical plot
# plot!(
#     verticle_line_x,
#     verticle_line_y,
#     linestyle=:dash,
#     linecolor=:black,
#     label=false,
#     ylims=(0,1.37e55),
#     lw=2
#     )

# saving image
png(plt, "iron_line_plot_spin.png")
display(plt)


text:
using Oceananigans.Advection:
    _advective_momentum_flux_Uu,
    _advective_momentum_flux_Uv,
    _advective_momentum_flux_Vu,
    _advective_momentum_flux_Vv

using Oceananigans.Grids: AbstractGrid
using Oceananigans.Operators: Ax_qᶠᶜᶜ, Ay_qᶜᶠᶜ

#####
##### Momentum flux operators
#####

@inline momentum_flux_huu(i, j, k, grid, advection, solution) =
    @inbounds _advective_momentum_flux_Uu(i, j, k, grid, advection, solution.uh, solution.uh) / solution.h[i, j, k]

@inline momentum_flux_hvu(i, j, k, grid, advection, solution) =
    @inbounds _advective_momentum_flux_Vu(i, j, k, grid, advection, solution.vh, solution.uh) / ℑxyᶠᶠᵃ(i, j, k, grid, solution.h)

@inline momentum_flux_huv(i, j, k, grid, advection, solution) =
    @inbounds _advective_momentum_flux_Uv(i, j, k, grid, advection, solution.uh, solution.vh) / ℑxyᶠᶠᵃ(i, j, k, grid, solution.h)

@inline momentum_flux_hvv(i, j, k, grid, advection, solution) =
    @inbounds _advective_momentum_flux_Vv(i, j, k, grid, advection, solution.vh, solution.vh) / solution.h[i, j, k]

#####
##### Momentum flux divergence operators
#####

@inline div_hUu(i, j, k, grid, advection, solution) =
    1 / Vᶠᶜᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, momentum_flux_huu, advection, solution) +
                               δyᵃᶜᵃ(i, j, k, grid, momentum_flux_hvu, advection, solution))

@inline div_hUv(i, j, k, grid, advection, solution) =
    1 / Vᶜᶠᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, momentum_flux_huv, advection, solution) +
                               δyᵃᶠᵃ(i, j, k, grid, momentum_flux_hvv, advection, solution))

# Support for no advection
@inline div_hUu(i, j, k, grid::AbstractGrid{FT}, ::Nothing, solution) where FT = zero(FT)
@inline div_hUv(i, j, k, grid::AbstractGrid{FT}, ::Nothing, solution) where FT = zero(FT)

#####
##### Mass transport divergence operator
#####

"""
    div_Uh(i, j, k, grid, solution)

Calculates the divergence of the mass flux into a cell,

    1/V * [δxᶜᵃᵃ(Ax * uh) + δyᵃᶜᵃ(Ay * vh)]

which will end up at the location `ccc`.
"""
@inline function div_Uh(i, j, k, grid, solution)
    1/Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_qᶠᶜᶜ, solution.uh) +
                             δyᵃᶜᵃ(i, j, k, grid, Ay_qᶜᶠᶜ, solution.vh))
end

#####
##### Tracer advection operator
#####

@inline transport_tracer_flux_x(i, j, k, grid, advection, uh, h, c) =
    @inbounds advective_tracer_flux_x(i, j, k, grid, advection, uh, c) / h[i, j, k]

@inline transport_tracer_flux_y(i, j, k, grid, advection, vh, h, c) =
    @inbounds advective_tracer_flux_y(i, j, k, grid, advection, vh, c) / h[i, j, k]

"""
    div_Uc(i, j, k, grid, advection, U, c)

Calculates the divergence of the flux of a tracer quantity c being advected by
a velocity field U = (u, v), ∇·(Uc),

    1/V * [δxᶜᵃᵃ(Ax * uh * ℑxᶠᵃᵃ(c) / h) + δyᵃᶜᵃ(Ay * vh * ℑyᵃᶠᵃ(c) / h)]

which will end up at the location `ccc`.
"""
@inline function div_Uc(i, j, k, grid, advection, solution, c)
    1/Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, transport_tracer_flux_x, advection, solution.uh, solution.h, c) +
                             δyᵃᶜᵃ(i, j, k, grid, transport_tracer_flux_y, advection, solution.vh, solution.h, c))
end

# Support for no advection
@inline div_Uc(i, j, k, grid::AbstractGrid{FT}, ::Nothing, solution, c) where FT = zero(FT)

@inline u(i, j, k, grid, solution) = @inbounds solution.uh[i, j, k] / solution.h[i, j, k]
@inline v(i, j, k, grid, solution) = @inbounds solution.vh[i, j, k] / solution.h[i, j, k]

"""
    c_div_U(i, j, k, grid, advection, U)

Calculates the product of the tracer concentration c with
the horizontal divergence of the velocity field U = (u, v), c ∇·(U),

    1/V * [δxᶜᵃᵃ(Ax * uh / h) + δyᵃᶜᵃ(Ay * vh / h]

which will end up at the location `ccc`.
"""
@inline function c_div_U(i, j, k, grid, solution, c)
    @inbounds c[i, j, k] * 1/Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, u, solution) + δyᵃᶜᵃ(i, j, k, grid, v, solution))
end

# Support for no advection
@inline c_div_Uc(i, j, k, grid::AbstractGrid{FT}, ::Nothing, solution, c) where FT = zero(FT)


text:
"""
Delegates loading of code with functions and parameters specific to different instruments.
Subdirectories of src/instruments include provide functions specialized for each instrument,
typically the file I/O and pre-processing, so data ends up in a common format.
src/instruments/common.jl provides routines that can be shared by instruments.
"""

include("common.jl")   # Mostly trait functions to be specialized by instruments

include("io.jl")
export read_manifest, read_header, read_metadata_from_fits

include("linelists.jl")
export read_linelist_espresso, read_linelist_vald

include("masks.jl")
export ChunkWidthFixedΔlnλ
export read_mask_espresso, read_mask_vald

#export predict_intrinsic_stellar_line_width
export λ_vac_to_air, λ_air_to_vac
export find_overlapping_chunks
export merge_chunks

include("param.jl")

#=
For each instrument/data file type, users need to create a sub-type of either AbstractInstrument2D or AbstractInstrument1D and
to implement the functions in neid/traits.jl for each instrument trait.
=#
function filter_line_list end
function find_worst_telluric_in_each_chunk end
function get_inst_module end

include("theory/theory.jl")
import .TheoreticalInstrument: TheoreticalInstrument1D, TheoreticalInstrument2D, AnyTheoreticalInstrument
import .TheoreticalInstrument: get_inst_module, filter_line_list, find_worst_telluric_in_each_chunk
export TheoreticalInstrument, TheoreticalInstrument1D, TheoreticalInstrument2D, AnyTheoreticalInstrument

include("neid/neid.jl")
import .NEID: NEID1D, NEID2D, AnyNEID
import .NEID: get_inst_module, filter_line_list, find_worst_telluric_in_each_chunk
export NEID, NEID1D, NEID2D, AnyNEID

include("expres/expres.jl")
import .EXPRES: EXPRES1D, EXPRES2D, AnyEXPRES, get_inst_module
import .EXPRES: get_inst_module, filter_line_list, find_worst_telluric_in_each_chunk
export EXPRES, EXPRES1D, EXPRES2D, AnyEXPRES

# TODO: Add more instruments:  HARPS-N, HPF, etc.
include("harps-n/harps-n.jl")
export HARPSN, HARPSN1D, HARPSN2D, AnyHARPSN
import .HARPSN: HARPSN1D, HARPSN2D, AnyHARPSN, get_inst_module
#import .HARPSN: get_inst_module, filter_line_list, find_worst_telluric_in_each_chunk

# using .TheoreticalInstrument, .NEID, .EXPRES, .HARPSN

include("tellurics.jl")
export find_worst_telluric_in_each_chunk, make_clean_line_list_from_tellurics_expres
export filter_line_list

export min_order, max_order, min_pixel_in_order, max_pixel_in_order, min_col_default, max_col_default
export orders_all, pixels_all, max_pixels_in_spectra       # generic implementations avaliable
export metadata_symbols_default, metadata_strings_default  # need to specialize
export default_ccf_mask_v_width
export get_inst_module

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
include("CG_Convergence_Utilities.jl")

# ==============================================================================
# Exact Arithematic error bound.
# ==============================================================================

function PerformCGFor(A, b, immitate_exact=true)
    cgm = ConjGradModified(A, b)
    if !immitate_exact
        TurnOffReorthgonalize(cgm)
    end
    ResidualNorm = GetResidualNorm(cgm)
    ResidualNorms = Vector{Float64}()
    Guesses = Vector{Vector{Float64}}()
    push!(ResidualNorms, ResidualNorm)
    while ResidualNorm > 1e-10
        ResidualNorm = cgm()
        push!(ResidualNorms, ResidualNorm)
        push!(Guesses, cgm.x)
    end
    return ResidualNorms, Guesses
end

function Run1()

return end

N = 512
A = GetUniformPSDMatrix(N)
b = rand(N)
ResidualNorms, Guesses = PerformCGFor(A, b)
ResidualEnergy = ResRelEnergyNorm(A, b,Guesses)
k = length(ResidualNorms)
ErrorsBound = [TheoreticalErrorBound(A, idx) for idx in 1: k]

fig1 = plot(log10.(ResidualEnergy), label="Relative Energy (exact)", legend=:bottomleft)
plot!(fig1, log10.(ErrorsBound), label="Theoretical Error Bound")
display(fig1)
SaveFigToCurrentScriptDir(fig1, "fig1.png")

# ==============================================================================
# Floating Points Error Bound.
# ==============================================================================

A = GetUniformPSDMatrix(N)
b = rand(N)
ResidualNorms, Guesses = PerformCGFor(A, b, false)
ResidualEnergy = ResRelEnergyNorm(A, b,Guesses)
k = length(ResidualNorms)
ErrorsBound = [TheoreticalErrorBound(A, idx) for idx in 1: k]

fig2 = plot(log10.(ResidualEnergy), label="Relative Energy (floats)", legend=:bottomleft)
plot!(fig2, log10.(ErrorsBound), label="Theoretical Error Bound(exact)")
display(fig2)
SaveFigToCurrentScriptDir(fig2, "fig2.png")


text:
# Copyright (c) 2020 California Institute of Technology (“Caltech”). U.S.
# Government sponsorship acknowledged.
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#   • Redistributions of source code must retain the above copyright notice,
#     this list of conditions and the following disclaimer.
#   • Redistributions in binary form must reproduce the above copyright notice,
#     this list of conditions and the following disclaimer in the documentation
#     and/or other materials provided with the distribution.
#   • Neither the name of Caltech nor its operating division, the Jet
#   Propulsion Laboratory, nor the names of its contributors may be used to
#   endorse or promote products derived from this software without specific
#   prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
# THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
# PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
# OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
# OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
# ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

module Covariance

export calc_covariance_modecounting,
       calc_covariance_efstathiou,
       calc_covariance_exact_chain

using ..Modes
using ..Windows
using ..SeparableArrays
using ..Theory
using ..WindowChains

function calc_covariance_modecounting(CBlnn, bcmodes, fvol, Δℓ, Δn)
    lnnsize = getlnnsize(bcmodes)
    V_mc = fill(0.0, lnnsize, lnnsize)
    @time for j=1:lnnsize, i=1:lnnsize
        L, K, K′ = getlkk(bcmodes, j)
        l, k, k′ = getlkk(bcmodes, i)
        C1 = 0.0
        C2 = 0.0
        C3 = 0.0
        C4 = 0.0
        # Need to ensure that what we calculate exists in 'CBlnn'
        idx = getidxapprox(bcmodes, l, k, K)
        if typeof(idx) <: Integer
            C1 = CBlnn[idx]
        end
        idx = getidxapprox(bcmodes, L, k′, K′)
        if typeof(idx) <: Integer
            C2 = CBlnn[idx]
        end
        idx = getidxapprox(bcmodes, l, k, K′)
        if typeof(idx) <: Integer
            C3 = CBlnn[idx]
        end
        idx = getidxapprox(bcmodes, L, k′, K)
        if typeof(idx) <: Integer
            C4 = CBlnn[idx]
        end
        #(i == j) && @show i,j,(l,k,k′),(L,K,K′),C1,C2,C3,C4
        if l == L
            Nmodes = fvol * (2*L+1) * Δℓ * Δn
            V_mc[i,j] = (C1*C2 + C3*C4) / Nmodes
        end
    end
    return V_mc
end

function calc_covariance_efstathiou(CBlnn, win, wmodes, w̃mat, vmat, bcmodes)
    # calculate covariance matrix
    win² = exponentiate(win, 2)
    @time bcmix⁽²⁾ᵃ = power_win_mix(win², w̃mat, vmat, wmodes, bcmodes, div2Lp1=true)
    @time bcmix⁽²⁾ᵇ = power_win_mix(win², w̃mat, vmat, wmodes, bcmodes, div2Lp1=true, interchange_NN′=true)
    #plot_cmix(bcmix⁽²⁾ᵃ, "\$M^{(2)a}\$", id*"_bcmix2a")
    #plot_cmix(bcmix⁽²⁾ᵇ, "\$M^{(2)b}\$", id*"_bcmix2b")

    lnnsize = getlnnsize(bcmodes)
    CBlnN_CBlnN_bcmix⁽²⁾ = fill(0.0, lnnsize, lnnsize)
    @time for j=1:lnnsize, i=1:lnnsize
        L, K, K′ = getlkk(bcmodes, j)
        l, k, k′ = getlkk(bcmodes, i)
        C1 = 0.0
        C2 = 0.0
        C3 = 0.0
        C4 = 0.0
        # Need to ensure that what we calculate exists in 'CBlnn'
        idx = getidxapprox(bcmodes, l, k, K)
        if typeof(idx) <: Integer
            C1 = CBlnn[idx]
        end
        idx = getidxapprox(bcmodes, L, k′, K′)
        if typeof(idx) <: Integer
            C2 = CBlnn[idx]
        end
        idx = getidxapprox(bcmodes, l, k, K′)
        if typeof(idx) <: Integer
            C3 = CBlnn[idx]
        end
        idx = getidxapprox(bcmodes, L, k′, K)
        if typeof(idx) <: Integer
            C4 = CBlnn[idx]
        end
        if i == j
            @show i,j,(l,k,k′),(L,K,K′),C1,C2,C3,C4
        end
        CBlnN_CBlnN_bcmix⁽²⁾[i,j] = C1*C2*bcmix⁽²⁾ᵃ[i,j] + C3*C4*bcmix⁽²⁾ᵇ[i,j]
    end
    return CBlnN_CBlnN_bcmix⁽²⁾
end

################ exact calculation ################

function calc_covariance_exact_chain(Clnn, nbar, win, wmodes, cmodes; Δℓ=1, Δn=1)
    if any(@. Clnn != 0) && nbar != 0
        @time VW = calc_covariance_exact_A1(Clnn, win, wmodes, cmodes, Δℓ=Δℓ, Δn=Δn)
        @time VW .+= calc_covariance_exact_A2(Clnn, nbar, win, wmodes, cmodes, Δℓ=Δℓ, Δn=Δn)
        @time VW .+= calc_covariance_exact_A3(nbar, win, wmodes, cmodes, Δℓ=Δℓ, Δn=Δn)
    else
        @time VW = calc_covariance_exact_A3(nbar, win, wmodes, cmodes, Δℓ=Δℓ, Δn=Δn)
    end
    return VW
end

function calc_covariance_exact_A1(Clnn, win, wmodes, cmodes; Δℓ=1, Δn=1)
    T = Float64
    amodes = cmodes.amodes
    wccache = WindowChainsCache(win, wmodes, cmodes.amodes)
    lnnsize = getlnnsize(cmodes)
    symmetries = [1=>1, 2=>2, 3=>1]
    A1 = fill(T(0), lnnsize, lnnsize)
    for j=1:lnnsize, i=j:lnnsize
        l, n, n′ = getlnn(cmodes, i)
        L, N, N′ = getlnn(cmodes, j)
        if abs(L-l) > Δℓ || abs(n-N) > Δn || abs(n′-N′) > Δn
            continue
        end
        @show i,(l,n,n′),j,(L,N,N′)
        ell = [0, L, 0, l]
        enn = [0, N, 0, n′]
        enn′ = [0, N′, 0, n]
        A = T(0)
        for k=1:lnnsize
            l1, n1, n2 = getlnn(cmodes, k)
            ell[1] = l1
            enn[1] = n1
            enn′[1] = n2
            W4 = fill(NaN, lnnsize)
            for m=1:lnnsize
                l3, n4, n3 = getlnn(cmodes, m)
                ell[3] = l3
                enn[3] = n4
                enn′[3] = n3
                W4[m] = window_chain(ell, enn, enn′, wccache, symmetries)
            end
            A += Clnn[k] * (Clnn' * W4)
        end
        A /= (2*l+1) * (2*L+1)
        A1[i,j] = A1[j,i] = A
    end
    return A1
end

function calc_covariance_exact_A2(Clnn, nbar, win, wmodes, cmodes; Δℓ=1, Δn=1)
    T = Float64
    amodes = cmodes.amodes
    wccache = WindowChainsCache(win, wmodes, cmodes.amodes)
    lnnsize = getlnnsize(cmodes)
    symmetries = [1=>1, 2=>2, 3=>2]
    A2 = fill(T(0), lnnsize, lnnsize)
    for j=1:lnnsize, i=j:lnnsize
        l, n, n′ = getlnn(cmodes, i)
        L, N, N′ = getlnn(cmodes, j)
        if abs(L-l) > Δℓ || abs(n-N) > Δn || abs(n′-N′) > Δn
            continue
        end
        @show i,(l,n,n′),j,(L,N,N′)
        ell = [0, L, l]
        enn = [0, N, n′]
        enn′ = [0, N′, n]
        W3 = fill(NaN, lnnsize)
        for k=1:lnnsize
            l1, n1, n2 = getlnn(cmodes, k)
            ell[1] = l1
            enn[1] = n1
            enn′[1] = n2
            W3[k] = window_chain(ell, enn, enn′, wccache, symmetries)
        end
        A = (Clnn' * W3) / ((2*l+1) * (2*L+1) * nbar)
        A2[i,j] = A2[j,i] = A
    end
    return A2
end

function calc_covariance_exact_A3(nbar, win, wmodes, cmodes; Δℓ=1, Δn=1)
    T = Float64
    amodes = cmodes.amodes
    wccache = WindowChainsCache(win, wmodes, cmodes.amodes)
    lnnsize = getlnnsize(cmodes)
    A3 = fill(T(0), lnnsize, lnnsize)
    for j=1:lnnsize, i=j:lnnsize
        l, n, n′ = getlnn(cmodes, i)
        L, N, N′ = getlnn(cmodes, j)
        if abs(L-l) > Δℓ || abs(n-N) > Δn || abs(n′-N′) > Δn
            continue
        end
        @show i,(l,n,n′),j,(L,N,N′)
        ell = [L, l]
        enn = [N, n′]
        enn′ = [N′, n]
        W2 = window_chain(ell, enn, enn′, wccache)
        if N == N′
            W2 *= 2
        else
            enn[1], enn′[1] = N′, N
            W2 += window_chain(ell, enn, enn′, wccache)
        end
        A = W2 / ((2*l+1) * (2*L+1) * nbar^2)
        A3[i,j] = A3[j,i] = A
    end
    return A3
end

end

# vim: set sw=4 et sts=4 :


text:
factor_ab(x) = logpdf(MvNormal(PDMat([1.0 0.5; 0.5 1.0])), x)
factor_bc(x) = logpdf(MvNormal([2.0, 2.0], PDMat([2.0 -1.0; -1.0 2.0])), x)

struct FactorProposal
    addresses :: Vector{String}
    last :: Vector{String}
    std :: Float64
end
FactorProposal(addresses) = FactorProposal(addresses, [addresses[1]], 0.5)
function (fp::FactorProposal)(x::Dict)
    new = deepcopy(x)
    address = rand(fp.addresses)
    new[address] = randn() * fp.std + new[address]
    fp.last[1] = address
    new
end
function (fp::FactorProposal)(x_prime::Dict, x::Dict)
    address = fp.last[1]
    logpdf(Normal(x[address], fp.std), x_prime[address]) - log(length(fp.addresses))
end

#=
@testset "make field and sample" begin
    factors = Dict(["a", "b"] => factor_ab, ["b", "c"] => factor_bc)
    evidence = Dict("b" => 3.0)
    #evidence = Dict()
    field = RandomField(factors, evidence)
    addresses = ["a", "c"]
    proposal = FactorProposal(addresses)
    initial_values = Dict("a" => 0.0, "b" => 3.0, "c" => 0.0)
    samples = mh(field, [proposal], initial_values; burn=1000, thin=100, num_iterations=11000)
    for address in ["a", "b", "c"]
        plot_marginal(samples, address, "out", "factor-marginal-$address.png")
    end
end
=#

@testset "using generative field" begin
    factors = Dict(["a", "b"] => factor_ab, ["b", "c"] => factor_bc)
    field = RandomField(factors)
    addresses = ["a", "b", "c"]
    proposal = FactorProposal(addresses)
    init = Dict("a" => 0.0, "b" => 0.0, "c" => 0.0)
    gf = GenerativeField(field, proposal, init)

    # generative fields can be used in generic stochastic functions
    t = trace()
    sample(t, "field", gf)
    @info t

    # note that we can pass upstream samped values as evidence
    t = trace()
    upstream_value = sample(t, "a", Normal())
    field = RandomField(factors, Dict("a" => upstream_value))
    proposal = FactorProposal(["b", "c"])
    init = Dict("a" => 0.0, "b" => 0.0, "c" => 0.0)
    gf = GenerativeField(field, proposal, init)
    sample(t, "field", gf)
    @info t
end