text:
export SymbolicPermutation

"""
A probability estimator based on permutations.
"""
abstract type PermutationProbabilityEstimator <: SymbolicProbabilityEstimator end

"""
    SymbolicPermutation(; τ = 1, m = 3, lt = Entropies.isless_rand) <: ProbabilityEstimator
    SymbolicWeightedPermutation(; τ = 1, m = 3, lt = Entropies.isless_rand) <: ProbabilityEstimator
    SymbolicAmplitudeAwarePermutation(; τ = 1, m = 3, A = 0.5, lt = Entropies.isless_rand) <: ProbabilityEstimator

Symbolic, permutation-based probabilities/entropy estimators.
`m` is the permutation order (or the symbol size or the embedding dimension) and
`τ` is the delay time (or lag).

## Repeated values during symbolization

In the original implementation of permutation entropy [^BandtPompe2002], equal values are
ordered after their order of appearance, but this can lead to erroneous temporal
correlations, especially for data with low-amplitude resolution [^Zunino2017]. Here, we
resolve this issue by letting the user provide a custom "less-than" function. The keyword
`lt` accepts a function that decides which of two state vector elements are smaller. If two
elements are equal, the default behaviour is to randomly assign one of them as the largest
(`lt = Entropies.isless_rand`). For data with low amplitude resolution, computing
probabilities multiple times using the random approach may reduce these erroneous
effects.

To get the behaviour described in Bandt and Pompe (2002), use `lt = Base.isless`).

## Properties of original signal preserved

- **`SymbolicPermutation`**: Preserves ordinal patterns of state vectors (sorting information). This
    implementation is based on Bandt & Pompe et al. (2002)[^BandtPompe2002] and
    Berger et al. (2019) [^Berger2019].
- **`SymbolicWeightedPermutation`**: Like `SymbolicPermutation`, but also encodes amplitude
    information by tracking the variance of the state vectors. This implementation is based
    on Fadlallah et al. (2013)[^Fadlallah2013].
- **`SymbolicAmplitudeAwarePermutation`**: Like `SymbolicPermutation`, but also encodes
    amplitude information by considering a weighted combination of *absolute amplitudes*
    of state vectors, and *relative differences between elements* of state vectors. See
    description below for explanation of the weighting parameter `A`. This implementation
    is based on Azami & Escudero (2016) [^Azami2016].

## Probability estimation

### Univariate time series

To estimate probabilities or entropies from univariate time series, use the following methods:

- `probabilities(x::AbstractVector, est::SymbolicProbabilityEstimator)`. Constructs state vectors
    from `x` using embedding lag `τ` and embedding dimension `m`, symbolizes state vectors,
    and computes probabilities as (weighted) relative frequency of symbols.
- `genentropy(x::AbstractVector, est::SymbolicProbabilityEstimator; α=1, base = 2)` computes
    probabilities by calling `probabilities(x::AbstractVector, est)`,
    then computer the order-`α` generalized entropy to the given base.

#### Speeding up repeated computations

A pre-allocated integer symbol array `s` can be provided to save some memory
allocations if the probabilities are to be computed for multiple data sets.

*Note: it is not the array that will hold the final probabilities that is pre-allocated,
but the temporary integer array containing the symbolized data points. Thus, if
provided, it is required that `length(x) == length(s)` if `x` is a Dataset, or
`length(s) == length(x) - (m-1)τ` if `x` is a univariate signal that is to be embedded
first*.

Use the following signatures (only works for `SymbolicPermutation`).

```julia
probabilities!(s::Vector{Int}, x::AbstractVector, est::SymbolicPermutation) → ps::Probabilities
probabilities!(s::Vector{Int}, x::AbstractDataset, est::SymbolicPermutation) → ps::Probabilities
```

### Multivariate datasets

Although not dealt with in the original paper describing the estimators, numerically speaking,
permutation entropies can also be computed for multivariate datasets with dimension ≥ 2
(but see caveat below). Such datasets may be, for example, preembedded time series. Then,
just skip the delay reconstruction step, compute and symbols directly from the ``L``
existing state vectors ``\\{\\mathbf{x}_1, \\mathbf{x}_2, \\ldots, \\mathbf{x_L}\\}``.

- `probabilities(x::AbstractDataset, est::SymbolicProbabilityEstimator)`. Compute ordinal patterns of the
    state vectors of `x` directly (without doing any embedding), symbolize those patterns,
    and compute probabilities as (weighted) relative frequencies of symbols.
- `genentropy(x::AbstractDataset, est::SymbolicProbabilityEstimator)`. Computes probabilities from
    symbol frequencies using `probabilities(x::AbstractDataset, est::SymbolicProbabilityEstimator)`,
    then computes the order-`α` generalized (permutation) entropy to the given base.

*Caveat: A dynamical interpretation of the permutation entropy does not necessarily
hold if computing it on generic multivariate datasets. Method signatures for `Dataset`s are
provided for convenience, and should only be applied if you understand the relation
between your input data, the numerical value for the permutation entropy, and
its interpretation.*

## Description

All symbolic estimators use the same underlying approach to estimating probabilities.

### Embedding, ordinal patterns and symbolization

Consider the ``n``-element univariate time series ``\\{x(t) = x_1, x_2, \\ldots, x_n\\}``.
Let ``\\mathbf{x_i}^{m, \\tau} = \\{x_j, x_{j+\\tau}, \\ldots, x_{j+(m-1)\\tau}\\}``
for ``j = 1, 2, \\ldots n - (m-1)\\tau`` be the ``i``-th state vector in a delay
reconstruction with embedding dimension ``m`` and reconstruction lag ``\\tau``.
There are then ``N = n - (m-1)\\tau`` state vectors.

For an ``m``-dimensional vector, there are ``m!`` possible ways of sorting it in
ascending order of magnitude. Each such possible sorting ordering is called a
*motif*. Let ``\\pi_i^{m, \\tau}`` denote the motif associated with the
``m``-dimensional state vector ``\\mathbf{x_i}^{m, \\tau}``, and let ``R``
be the number of distinct motifs that can be constructed from the ``N`` state
vectors. Then there are at most ``R`` motifs; ``R = N`` precisely when all motifs
are unique, and ``R = 1`` when all motifs are the same.

Each unique motif ``\\pi_i^{m, \\tau}`` can be mapped to a unique integer
symbol ``0 \\leq s_i \\leq M!-1``. Let ``S(\\pi) : \\mathbb{R}^m \\to \\mathbb{N}_0`` be
the function that maps the motif ``\\pi`` to its symbol ``s``, and let ``\\Pi``
denote the set of symbols ``\\Pi = \\{ s_i \\}_{i\\in \\{ 1, \\ldots, R\\}}``.

### Probability computation

#### `SymbolicPermutation`

The probability of a given motif is its frequency of occurrence, normalized by the total
number of motifs (with notation from [^Fadlallah2013]),

```math
p(\\pi_i^{m, \\tau}) = \\dfrac{\\sum_{k=1}^N \\mathbf{1}_{u:S(u) = s_i} \\left(\\mathbf{x}_k^{m, \\tau} \\right) }{\\sum_{k=1}^N \\mathbf{1}_{u:S(u) \\in \\Pi} \\left(\\mathbf{x}_k^{m, \\tau} \\right)} = \\dfrac{\\sum_{k=1}^N \\mathbf{1}_{u:S(u) = s_i} \\left(\\mathbf{x}_k^{m, \\tau} \\right) }{N},
```

where the function ``\\mathbf{1}_A(u)`` is the indicator function of a set ``A``. That
    is, ``\\mathbf{1}_A(u) = 1`` if ``u \\in A``, and ``\\mathbf{1}_A(u) = 0`` otherwise.

#### `SymbolicAmplitudeAwarePermutation`

Amplitude-aware permutation entropy is computed analogously to regular permutation entropy
but probabilities are weighted by amplitude information as follows.

```math
p(\\pi_i^{m, \\tau}) = \\dfrac{\\sum_{k=1}^N \\mathbf{1}_{u:S(u) = s_i} \\left( \\mathbf{x}_k^{m, \\tau} \\right) \\, a_k}{\\sum_{k=1}^N \\mathbf{1}_{u:S(u) \\in \\Pi} \\left( \\mathbf{x}_k^{m, \\tau} \\right) \\,a_k} = \\dfrac{\\sum_{k=1}^N \\mathbf{1}_{u:S(u) = s_i} \\left( \\mathbf{x}_k^{m, \\tau} \\right) \\, a_k}{\\sum_{k=1}^N a_k}.
```

The weights encoding amplitude information about state vector ``\\mathbf{x}_i = (x_1^i, x_2^i, \\ldots, x_m^i)`` are

```math
a_i = \\dfrac{A}{m} \\sum_{k=1}^m |x_k^i | + \\dfrac{1-A}{d-1} \\sum_{k=2}^d |x_{k}^i - x_{k-1}^i|,
```

with ``0 \\leq A \\leq 1``. When ``A=0`` , only internal differences between the elements of
``\\mathbf{x}_i`` are weighted. Only mean amplitude of the state vector
elements are weighted when ``A=1``. With, ``0<A<1``, a combined weighting is used.

#### `SymbolicWeightedPermutation`

Weighted permutation entropy is also computed analogously to regular permutation entropy, but
adds weights that encode amplitude information too:

```math
p(\\pi_i^{m, \\tau}) = \\dfrac{\\sum_{k=1}^N \\mathbf{1}_{u:S(u) = s_i}
\\left( \\mathbf{x}_k^{m, \\tau} \\right)
\\, w_k}{\\sum_{k=1}^N \\mathbf{1}_{u:S(u) \\in \\Pi}
\\left( \\mathbf{x}_k^{m, \\tau} \\right) \\,w_k} = \\dfrac{\\sum_{k=1}^N
\\mathbf{1}_{u:S(u) = s_i}
\\left( \\mathbf{x}_k^{m, \\tau} \\right) \\, w_k}{\\sum_{k=1}^N w_k}.
```

The weighted permutation entropy is equivalent to regular permutation entropy when weights
are positive and identical (``w_j = \\beta \\,\\,\\, \\forall \\,\\,\\, j \\leq N`` and
``\\beta > 0)``. Weights are dictated by the variance of the state vectors.

Let the aritmetic mean of state vector ``\\mathbf{x}_i`` be denoted
by

```math
\\mathbf{\\hat{x}}_j^{m, \\tau} = \\frac{1}{m} \\sum_{k=1}^m x_{j + (k+1)\\tau}.
```

Weights are then computed as

```math
w_j = \\dfrac{1}{m}\\sum_{k=1}^m (x_{j+(k+1)\\tau} - \\mathbf{\\hat{x}}_j^{m, \\tau})^2.
```

*Note: in equation 7, section III, of the original paper, the authors write*

```math
w_j = \\dfrac{1}{m}\\sum_{k=1}^m (x_{j-(k-1)\\tau} - \\mathbf{\\hat{x}}_j^{m, \\tau})^2.
```
*But given the formula they give for the arithmetic mean, this is **not** the variance
of ``\\mathbf{x}_i``, because the indices are mixed: ``x_{j+(k-1)\\tau}`` in the weights
formula, vs. ``x_{j+(k+1)\\tau}`` in the arithmetic mean formula. This seems to imply
that amplitude information about previous delay vectors
are mixed with mean amplitude information about current vectors. The authors also mix the
terms "vector" and "neighboring vector" (but uses the same notation for both), making it
hard to interpret whether the sign switch is a typo or intended. Here, we use the notation
above, which actually computes the variance for ``\\mathbf{x}_i``*.

[^BandtPompe2002]: Bandt, Christoph, and Bernd Pompe. "Permutation entropy: a natural
    complexity measure for time series." Physical review letters 88.17 (2002): 174102.
[^Berger2019]: Berger, Sebastian, et al. "Teaching Ordinal Patterns to a Computer: Efficient Encoding Algorithms Based on the Lehmer Code." Entropy 21.10 (2019): 1023.
[^Fadlallah2013]: Fadlallah, Bilal, et al. "Weighted-permutation entropy: A complexity
    measure for time series incorporating amplitude information." Physical
    Review E 87.2 (2013): 022911.
[^Rényi1960]: A. Rényi, *Proceedings of the fourth Berkeley Symposium on Mathematics, Statistics and Probability*, pp 547 (1960)
[^Azami2016]: Azami, H., & Escudero, J. (2016). Amplitude-aware permutation entropy: Illustration in spike detection and signal segmentation. Computer methods and programs in biomedicine, 128, 40-51.
[^Fadlallah2013]: Fadlallah, Bilal, et al. "Weighted-permutation entropy: A complexity measure for time series incorporating amplitude information." Physical Review E 87.2 (2013): 022911.
[^Zunino2017]: Zunino, L., Olivares, F., Scholkmann, F., & Rosso, O. A. (2017). Permutation entropy based time series analysis: Equalities in the input signal can lead to false conclusions. Physics Letters A, 381(22), 1883-1892.
"""
struct SymbolicPermutation{F} <: PermutationProbabilityEstimator
    τ::Int
    m::Int
    lt::F
end
function SymbolicPermutation(; τ::Int = 1, m::Int = 3, lt::F=isless_rand) where {F <: Function}
    m >= 2 || error("Need m ≥ 2, otherwise no dynamical information is encoded in the symbols.")
    SymbolicPermutation{F}(τ, m, lt)
end

"""
    symbolize(x::AbstractVector{T}, est::SymbolicPermutation) where {T} → Vector{Int}
    symbolize!(s, x::AbstractVector{T}, est::SymbolicPermutation) where {T} → Vector{Int}

If `x` is a univariate time series, first `x` create a delay reconstruction of `x`
using embedding lag `est.τ` and embedding dimension `est.m`, then symbolizing the resulting
state vectors with [`encode_motif`](@ref).

Optionally, the in-place `symbolize!` can be used to put symbols in a pre-allocated
integer vector `s`, where `length(s) == length(x)-(est.m-1)*est.τ`.

    symbolize(x::AbstractDataset{m, T}, est::SymbolicPermutation) where {m, T} → Vector{Int}
    symbolize!(s, x::AbstractDataset{m, T}, est::SymbolicPermutation) where {m, T} → Vector{Int}

If `x` is an `m`-dimensional dataset, then motif lengths are determined by the dimension of
the input data, and `x` is symbolized by converting each `m`-dimensional
state vector as a unique integer in the range ``1, 2, \\ldots, m-1``, using
[`encode_motif`](@ref).

Optionally, the in-place `symbolize!` can be used to put symbols in a pre-allocated
integer vector `s`, where `length(s) == length(s)`.

## Examples

Symbolize a 7-dimensional dataset. Motif lengths (or order of the permutations) are
inferred to be 7.

```julia
using DelayEmbeddings, Entropies
D = Dataset([rand(7) for i = 1:1000])
s = symbolize(D, SymbolicPermutation())
```

Symbolize a univariate time series by first embedding it in dimension 5 with embedding lag 2.
Motif lengths (or order of the permutations) are therefore 5.

```julia
using DelayEmbeddings, Entropies
n = 5000
x = rand(n)
s = symbolize(x, SymbolicPermutation(m = 5, τ = 2))
```

The integer vector `s` now has length `n-(m-1)*τ = 4992`, and each `s[i]` contains
the integer symbol for the ordinal pattern of state vector `x[i]`.

[^Berger2019]: Berger, Sebastian, et al. "Teaching Ordinal Patterns to a Computer: Efficient Encoding Algorithms Based on the Lehmer Code." Entropy 21.10 (2019): 1023.
"""
function symbolize(x::AbstractDataset{m, T}, est::PermutationProbabilityEstimator) where {m, T}
    m >= 2 || error("Data must be at least 2-dimensional to symbolize. If data is a univariate time series, embed it using `genembed` first.")
    s = zeros(Int, length(x))
    symbolize!(s, x, est)
    return s
end

function symbolize(x::AbstractVector{T}, est::PermutationProbabilityEstimator) where {T}
    τs = tuple([est.τ*i for i = 0:est.m-1]...)
    x_emb = genembed(x, τs)

    s = zeros(Int, length(x_emb))
    symbolize!(s, x_emb, est)
    return s
end

function fill_symbolvector!(s, x, sp, m::Int; lt::Function = isless_rand)
    @inbounds for i = 1:length(x)
        sortperm!(sp, x[i], lt = lt)
        s[i] = encode_motif(sp, m)
    end
end

function symbolize!(s::AbstractVector{Int}, x::AbstractDataset{m, T}, est::SymbolicPermutation) where {m, T}
    @assert length(s) == length(x)
    #=
    Loop over embedding vectors `E[i]`, find the indices `p_i` that sort each `E[i]`,
    then get the corresponding integers `k_i` that generated the
    permutations `p_i`. Those integers are the symbols for the embedding vectors
    `E[i]`.
    =#
    sp = zeros(Int, m) # pre-allocate a single symbol vector that can be overwritten.
    fill_symbolvector!(s, x, sp, m, lt = est.lt)

    return s
end

function symbolize!(s::AbstractVector{Int}, x::AbstractVector{T}, est::SymbolicPermutation) where T
    τs = tuple([est.τ*i for i = 0:est.m-1]...)
    x_emb = genembed(x, τs)
    symbolize!(s, x_emb, est)
end

function probabilities!(s::AbstractVector{Int}, x::AbstractDataset{m, T}, est::SymbolicPermutation) where {m, T}
    length(s) == length(x) || throw(ArgumentError("Need length(s) == length(x), got `length(s)=$(length(s))` and `length(x)==$(length(x))`."))
    m >= 2 || error("Data must be at least 2-dimensional to compute the permutation entropy. If data is a univariate time series embed it using `genembed` first.")

    @inbounds for i = 1:length(x)
        s[i] = encode_motif(x[i], m)
    end
    probabilities(s)
end

function probabilities!(s::AbstractVector{Int}, x::AbstractVector{T}, est::SymbolicPermutation) where {T<:Real}
    L = length(x)
    N = L - (est.m-1)*est.τ
    length(s) == N || error("Pre-allocated symbol vector `s`needs to have length `length(x) - (m-1)*τ` to match the number of state vectors after `x` has been embedded. Got length(s)=$(length(s)) and length(x)=$(L).")

    τs = tuple([est.τ*i for i = 0:est.m-1]...)
    x_emb = genembed(x, τs)

    probabilities!(s, x_emb, est)
end

function probabilities(x::AbstractDataset{m, T}, est::SymbolicPermutation) where {m, T}
    s = zeros(Int, length(x))
    probabilities!(s, x, est)
end

function probabilities(x::AbstractVector{T}, est::SymbolicPermutation) where {T<:Real}
    τs = tuple([est.τ*i for i = 0:est.m-1]...)
    x_emb = genembed(x, τs)

    s = zeros(Int, length(x_emb))
    probabilities!(s, x_emb, est)
end

function genentropy!(
    s::AbstractVector{Int}, x::AbstractDataset{m, T}, est::SymbolicPermutation;
    q = 1.0, α = nothing, base::Real = MathConstants.e
    ) where {m, T}

    if α ≠ nothing
        @warn "Keyword `α` is deprecated in favor of `q`."
        q = α
    end
    length(s) == length(x) || error("Pre-allocated symbol vector s need the same number of elements as x. Got length(s)=$(length(s)) and length(x)=$(L).")
    ps = probabilities!(s, x, est)

    genentropy(ps, α = α, base = base)
end

function genentropy!(
        s::AbstractVector{Int}, x::AbstractVector{T}, est::SymbolicPermutation;
        q::Real = 1.0, α = nothing, base::Real = MathConstants.e
    ) where {T<:Real}
    if α ≠ nothing
        @warn "Keyword `α` is deprecated in favor of `q`."
        q = α
    end
    L = length(x)
    N = L - (est.m-1)*est.τ
    length(s) == N || error("Pre-allocated symbol vector `s` needs to have length `length(x) - (m-1)*τ` to match the number of state vectors after `x` has been embedded. Got length(s)=$(length(s)) and length(x)=$(L).")

    ps = probabilities!(s, x, est)
    genentropy(ps, α = α, base = base)
end


text:
using DataFrames

export EVI, log_likelihood_per_instance, log_likelihood, log_likelihood_avg

"""
Compute the likelihood of the PC given each individual instance in the data
"""
function log_likelihood_per_instance(spc::SharedProbCircuit, data, weights::Union{AbstractArray, Nothing} = nothing; use_gpu::Bool = isgpu(data), component_idx = 0)::AbstractVector
    if component_idx != 0
        ll, _, _ = EVI(ParamBitCircuit(spc, data; component_idx = component_idx), data)
        return ll
    end
    ll_per_instance = zeros(Float64, num_examples(data))
    v, f = nothing, nothing
    if isnothing(weights)
        for component_idx = 1 : num_components(spc)
            bc = ParamBitCircuit(spc, data; component_idx = component_idx)
            ll, v, f = log_likelihood_per_instance(bc, data, v, f; use_gpu)
            ll = isgpu(ll) ? to_cpu(ll) : ll
            @inbounds @views ll_per_instance[:] .= ll_per_instance[:] .+ ll[:]
        end
        @inbounds @views ll_per_instance[:] .= ll_per_instance[:] ./ num_components(spc)
        return ll_per_instance
    end
    @assert length(weights) == num_components(spc) "Weights must have same dimension as number of components."
    if weights isa AbstractVector weights = reshape(weights, (1, length(weights))) end
    n = length(weights)
    lls = Array{Float64, 2}(undef, num_examples(data), n)
    for i in 1:n
        pbc = ParamBitCircuit(spc, data; component_idx = i)
        ll, v, f = log_likelihood_per_instance(pbc, data, v, f; use_gpu)
        lls[:,i] .= isgpu(ll) ? to_cpu(ll) : ll
    end
    return logsumexp(lls .+ log.(weights), 2)
end

function log_likelihood_per_instance(pc::ProbCircuit, data; use_gpu::Bool = isgpu(data))
    if isweighted(data)
        data, weights = split_sample_weights(data)
    end
    bc = ParamBitCircuit(pc, data)

    if isbatched(data)
        v, f = nothing, nothing
        lls = map(data) do d
            ll, v, f = log_likelihood_per_instance(bc, d; use_gpu)

            ll
        end

        ll = cat(lls...; dims = 1)
    else
        ll, v, f = log_likelihood_per_instance(bc, data; use_gpu)
    end

    if use_gpu
        CUDA.unsafe_free!(v) # save the GC some effort
        CUDA.unsafe_free!(f) # save the GC some effort
    end

    ll
end

function log_likelihood_per_instance(bc::ParamBitCircuit, data, reuse_v = nothing, reuse_f = nothing;
                                     use_gpu::Bool = isgpu(data))
    @assert isbinarydata(data) || isfpdata(data) "Probabilistic circuit likelihoods are for binary/float data only"
    if isgpu(data)
        use_gpu = true
    end
    if use_gpu
        ll, v, f = log_likelihood_per_instance_gpu(to_gpu(bc), to_gpu(data), reuse_v, reuse_f)
    else
        ll, v, f = log_likelihood_per_instance_cpu(bc, data, reuse_v, reuse_f)
    end

    ll, v, f
end

function log_likelihood_per_instance_cpu(bc, data, reuse_v = nothing, reuse_f = nothing)
    ll::Vector{Float64} = zeros(Float64, num_examples(data))
    ll_lock::Threads.ReentrantLock = Threads.ReentrantLock()

    @inline function on_edge(flows, values, prime, sub, element, grandpa, single_child, weights)
        if !single_child
            lock(ll_lock) do # TODO: move lock to inner loop? change to atomic float?
                for i = 1:size(flows,1)
                    @inbounds edge_flow = values[i, prime] & values[i, sub] & flows[i, grandpa]
                    first_true_bit = trailing_zeros(edge_flow)+1
                    last_true_bit = 64-leading_zeros(edge_flow)
                    @simd for j = first_true_bit:last_true_bit
                        ex_id = ((i-1) << 6) + j
                        if get_bit(edge_flow, j)
                            @inbounds ll[ex_id] += bc.params[element]
                        end
                    end
                end
            end
        end
        nothing
    end
    @inline function on_edge(flows, values::Matrix{<:AbstractFloat}, prime, sub, element, grandpa, single_child, weights)
        if !single_child
            lock(ll_lock) do # TODO: move lock to inner loop? change to atomic float?
                for i = 1:size(flows,1)
                    @inbounds edge_flow = values[i, prime] * values[i, sub] / values[i, grandpa] * flows[i, grandpa]
                    @inbounds ll[i] += bc.params[element] * edge_flow
                end
            end
        end
        nothing
    end

    v, f = satisfies_flows(bc.bitcircuit, data, reuse_v, reuse_f; on_edge)

    if isbinarydata(data)
        # when the example is outside of the support, give 0 likelihood
        in_support = AbstractBitVector(v[:,end], num_examples(data))
        ll[.! in_support] .= -Inf
    else
        @assert isfpdata(data)

        @inbounds @views in_support = v[:,end]
        ll[in_support .< 1e-8] .= -Inf
    end

    return ll, v, f
end

function log_likelihood_per_instance_gpu(bc, data, reuse_v = nothing, reuse_f = nothing)
    params_device = CUDA.cudaconvert(bc.params)
    ll::CuVector{Float64} = CUDA.zeros(Float64, num_examples(data))
    ll_device = CUDA.cudaconvert(ll)

    @inline function on_edge(flows, values, prime, sub, element, grandpa, chunk_id, edge_flow, single_child, weights)
        if !single_child
            first_true_bit = 1+trailing_zeros(edge_flow)
            last_true_bit = 64-leading_zeros(edge_flow)
            for j = first_true_bit:last_true_bit
                ex_id = ((chunk_id-1) << 6) + j
                if get_bit(edge_flow, j)
                    CUDA.@atomic ll_device[ex_id] += params_device[element]
                end
            end
        end
        nothing
    end
    @inline function on_edge(flows, values, prime, sub, element, grandpa, sample_idx, edge_flow::AbstractFloat,
                             single_child, weights)
        if !single_child
            c::Float64 = edge_flow * params_device[element]
            CUDA.@atomic ll_device[sample_idx] += c
        end
        nothing
    end

    v, f = satisfies_flows(bc.bitcircuit, data, reuse_v, reuse_f; on_edge)

    if isbinarydata(data)
        # when the example is outside of the support, give 0 likelihood
        # lazy programmer: do the conversion to a Vector{Bool} on CPU
        #   so that CUDA.jl can build a quick kernel
        # TODO: write a function to do this conversion on GPU directly
        in_support = AbstractBitVector(to_cpu(v[:,end]), num_examples(data))
        in_support = to_gpu(convert(Vector{Bool}, in_support))
    else
        @assert isfpdata(data)
        @inbounds @views in_support = v[:,end] .> 1e-8
    end
    ll2 = map((x,s) -> s ? x : -Inf, ll, in_support)

    CUDA.unsafe_free!(ll) # save the GC some effort

    return ll2, v, f
end

"""
    EVI(pc, data)

Computes the log likelihood data given full evidence.
Outputs `` \\log{p(x)} `` for each datapoint.
"""
const EVI = log_likelihood_per_instance

"""
    log_likelihood(pc, data)

Compute the likelihood of the PC given the data
"""
log_likelihood(pc, data; use_gpu::Bool = isgpu(data)) = begin
    if isweighted(data)
        # `data' is weighted according to its `weight' column
        data, weights = split_sample_weights(data)

        log_likelihood(pc, data, weights; use_gpu)
    else
        likelihoods = log_likelihood_per_instance(pc, data; use_gpu)

        isgpu(likelihoods) ? sum(to_cpu(likelihoods)) : sum(likelihoods)
    end
end
log_likelihood(pc, data, weights::DataFrame; use_gpu::Bool = isgpu(data)) = log_likelihood(pc, data, weights[:, 1]; use_gpu)
log_likelihood(pc, data, weights::AbstractArray; use_gpu::Bool = isgpu(data)) = begin
    if isgpu(weights)
        weights = to_cpu(weights)
    end
    likelihoods = log_likelihood_per_instance(pc, data; use_gpu)
    likelihoods = isgpu(likelihoods) ? to_cpu(likelihoods) : likelihoods
    mapreduce(*, +, likelihoods, weights)
end
log_likelihood(pc, data::Vector{DataFrame}; use_gpu::Bool = isgpu(data)) = begin
    if pc isa SharedProbCircuit
        total_ll = 0.0
        for component_idx = 1 : num_components(pc)
            total_ll += log_likelihood_batched(pc, data; use_gpu, component_idx)
        end
        total_ll / num_components(pc)
    else
        log_likelihood_batched(pc, data; use_gpu)
    end
end
log_likelihood_batched(pc, data::Vector{DataFrame}; use_gpu::Bool = isgpu(data), component_idx::Integer = 0) = begin
    # mapreduce(d -> log_likelihood(pc, d; use_gpu), +, data)
    if pc isa SharedProbCircuit
        pbc = ParamBitCircuit(pc, data; component_idx)
    else
        pbc = ParamBitCircuit(pc, data)
    end

    total_ll::Float64 = 0.0
    if isweighted(data)
        data, weights = split_sample_weights(data)

        v, f = nothing, nothing
        for idx = 1 : length(data)
            likelihoods, v, f = log_likelihood_per_instance(pbc, data[idx], v, f; use_gpu)
            if isgpu(likelihoods)
                likelihoods = to_cpu(likelihoods)
            end
            w = weights[idx]
            if isgpu(w)
                w = to_cpu(w)
            end
            total_ll += mapreduce(*, +, likelihoods, w)
        end
    else
        v, f = nothing, nothing
        for idx = 1 : length(data)
            likelihoods, v, f = log_likelihood_per_instance(pbc, data[idx], v, f; use_gpu)
            total_ll += sum(isgpu(likelihoods) ? to_cpu(likelihoods) : likelihoods)
        end
    end

    if use_gpu
        CUDA.unsafe_free!(v) # save the GC some effort
        CUDA.unsafe_free!(f) # save the GC some effort
    end

    total_ll
end

"""
    log_likelihood_avg(pc, data)

Compute the likelihood of the PC given the data, averaged over all instances in the data
"""
log_likelihood_avg(pc, data; use_gpu::Bool = isgpu(data)) = begin
    if isweighted(data)
        # `data' is weighted according to its `weight' column
        data, weights = split_sample_weights(data)

        log_likelihood_avg(pc, data, weights; use_gpu)
    else
        log_likelihood(pc, data; use_gpu) / num_examples(data)
    end
end
log_likelihood_avg(pc, data, weights::DataFrame; use_gpu::Bool = isgpu(data)) = log_likelihood_avg(pc, data, weights[:, 1]; use_gpu)
log_likelihood_avg(pc, data, weights; use_gpu::Bool = isgpu(data)) = begin
    if isgpu(weights)
        weights = to_cpu(weights)
    end
    log_likelihood(pc, data, weights; use_gpu) / sum(weights)
end
log_likelihood_avg(pc, data::Vector{DataFrame}; use_gpu::Bool = isgpu(data)) = begin
    if isweighted(data)
        weights = get_weights(data)
        if isgpu(weights)
            weights = to_cpu(weights)
        end
        log_likelihood(pc, data; use_gpu) / mapreduce(sum, +, weights)
    else
        log_likelihood(pc, data; use_gpu) / num_examples(data)
    end
end

ext: jl


text:
module PandorasBoxAirBubbles

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/airBubbles" AirBubbles(x::Integer, y::Integer, oneUse::Bool=false)

const placements = Ahorn.PlacementDict(
    "Air Bubbles (Pandora's Box)" => Ahorn.EntityPlacement(
        AirBubbles
    )
)

sprite = "objects/pandorasBox/airBubbles/idle00"

function Ahorn.selection(entity::AirBubbles)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::AirBubbles, room::Maple.Room)
    Ahorn.drawSprite(ctx, sprite, 0, 0)
end

end

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    TauModel()

Probabilistic rank aggregation with the Tau model.

## References

* Journel 2002. Combining Knowledge From Diverse Sources:
  An Alternative to Traditional Data Independence Hypotheses.
"""
struct TauModel <: AggregationMethod end

function rank_impl(objects, scores::NTuple{N,Symbol},
                   method::TauModel) where {N}
  # score columns
  cols = []
  for col in propertynames(objects)
    if col ∈ scores
      push!(cols, getproperty(objects, col))
    end
  end

  # conditional probabilities
  S = reduce(hcat, cols)
  P = S ./ sum(S, dims=1)
  n = size(P, 1)

  # uniform prior
  xₒ = (1 - 1/n) / (1/n)

  # odds with no redundancy
  X = (1 .- P) ./ P
  x = xₒ * prod(X/xₒ, dims=2)

  # posterior probabilities
  p = 1 ./ (1 .+ x)

  # final rank
  sortperm(vec(p), rev=true)
end


text:
import Pkg
Pkg.activate(joinpath(@__DIR__, "../../"))

using BSON
using CSV
using Distributions
using DelimitedFiles
using DataFrames
using Random
using RCall

@rimport base
@rimport FlowSOM
@rimport flowCore

include(joinpath(@__DIR__, "../PlotUtils/PlotUtils.jl"))

import PyPlot
const plt = PyPlot.plt
PyPlot.matplotlib.use("Agg")

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
using WriteVTK

function writeVtu(fileName, fields...)
    println("Writing output in vtu format")
    x = [mesh.nodes[i].centroid[1] for i in 1:length(mesh.nodes)]
    y = [mesh.nodes[i].centroid[2] for i in 1:length(mesh.nodes)]
    z = [mesh.nodes[i].centroid[3] for i in 1:length(mesh.nodes)]
    nodes = transpose(hcat(x, y, z))

    cells = VTKPolyhedron[]
    for iCell in 1:mesh.nCells
        cell = mesh.cells[iCell]
        faces = []
        for iFace in cell.iFaces
            push!(faces, mesh.faces[iFace].iNodes)
        end
        polyCell = VTKPolyhedron(cell.iNodes, faces...)
        push!(cells, polyCell)
    end

    for iFace in 1:mesh.nBoundaryFaces
        face = mesh.faces[mesh.nInternalFaces + iFace]
        lines = []
        for i in 1:length(face.iNodes)
            iStart = face.iNodes[i]
            iEnd = i < length(face.iNodes) ? face.iNodes[i + 1] : face.iNodes[1]
            push!(lines, [iEnd, iStart, iStart, iEnd])
        end
        polyFace = VTKPolyhedron(face.iNodes, lines...)
        push!(cells, polyFace)
    end

    vtu = vtk_grid(fileName, nodes, cells)
    for field in fields
        if field.type == "scalar"
            vtu["$(field.name)"] = field.cellValues
        elseif field.type == "vector"
            vtu["$(field.name)"] = transpose(mapreduce(permutedims, vcat, field.cellValues))
        end
    end
    vtk_save(vtu)
end


text:
# # Wavefield integration
#
# The process of mode conversion between two different segments of
# [`HomogeneousWaveguide`](@ref)'s
# requires the computation of electromagnetic wavefields from ground to a great
# height in the ionosphere.
#
# [Pitteway (1965)](https://doi.org/10.1098/rsta.1965.0004) presents one well-known method of
# integrating the wavefields that breaks the solution for the differential wave
# equations in the anisotropic ionosphere into solutions for a _penetrating_ and
# _non-penetrating_ mode. Because of the amplitude of the two waves varies greatly
# during the integration, Pitteway scales and orthogonalizes the solutions to maintain
# their independence.
#
# In this example we'll reproduce the figures in Pitteway's paper.
#
# First we'll import the necessary packages.

using CSV, Interpolations
using Plots
using Plots.Measures
using OrdinaryDiffEq

using TimerOutputs

using LongwaveModePropagator
using LongwaveModePropagator: QE, ME
const LMP = LongwaveModePropagator
nothing  #hide

# ## The ionosphere

# Pitteway uses an ionosphere presented in
# [Piggot et. al., 1965](https://doi.org/10.1098/rsta.1965.0005).
#
# ```@raw html
# <img src="../../images/Piggott_ionosphere.png"/>
# ```
#
# He begins with the midday profile with a 16 kHz radio wave at an angle of incidence of
# 40° from normal. We'll also assume the magnetic field has a strength of 50,000 nT
# at a dip angle of 68° and azimuth of 111°. Piggott was concerned with a path from Rugby
# to Cambridge, UK, which, according to the
# [Westinghouse VLF effective-conductivity map](https://apps.dtic.mil/sti/citations/AD0675771),
# has a ground conductivity index of 8.
# This can be specified as `GROUND[8]`.

ea = EigenAngle(deg2rad(40))
frequency = Frequency(16e3)
bfield = BField(50e-6, deg2rad(68), deg2rad(111))
ground = GROUND[8]

# To define the electron density and collision frequency profile, we will read in a
# [digitized table](@__REPO_ROOT_URL__/examples/piggott1965_data.csv)
# of the curves from Piggott and linearly interpolate them.
# While we're at it, we'll also prepare the nighttime ionosphere.

root_dir = dirname(dirname(pathof(LongwaveModePropagator)))
examples_dir = joinpath(root_dir, "examples")
data = CSV.File(joinpath(examples_dir, "piggott1965_data.csv"))

## interpolation object
day_itp = interpolate((data.day_ne_z,), data.day_ne, Gridded(Linear()))
day_etp = extrapolate(day_itp, Line())
dayfcn(z) = (v = day_etp(z); v > 0 ? v : 0.001)

## night has some rows of `missing` data
filtnight_z = collect(skipmissing(data.night_ne_z))
filtnight = collect(skipmissing(data.night_ne))
night_itp = interpolate((filtnight_z,), filtnight, Gridded(Linear()))
night_etp = extrapolate(night_itp, Line())
nightfcn(z) = (v = night_etp(z); v > 0 ? v : 0.001)

## so does the collision frequency
filtnu_z = collect(skipmissing(data.nu_z))
filtnu = collect(skipmissing(data.nu))
nu_itp = interpolate((filtnu_z,), filtnu, Gridded(Linear()))
nu_etp = extrapolate(nu_itp, Line())
nufcn(z) = (v = nu_etp(z); v > 0 ? v : 0.001)

day = Species(QE, ME, dayfcn, nufcn)
night = Species(QE, ME, nightfcn, nufcn)

nothing  #hide

# ## Scaled, integrated wavefields
#
# We use the coordinate frame where ``z`` is directed upward into the ionosphere, ``x`` is
# along the propagation direction, and ``y`` is perpendicular to complete the
# right-handed system.
# Where the ionosphere only varies in the ``z`` direction, the ``E_z``
# and ``H_z`` fields can be eliminated so that we only need to study the remaining four.
# The differential equations for the wave in the ionosphere can be written in matrix
# form for our coordinate system as
# ```math
# \frac{\mathrm{d}\bm{e}}{\mathrm{d}z} = -ik\bm{T}\bm{e}
# ```
# where
# ```math
# \bm{e} = \begin{pmatrix}
# E_x \\ -E_y \\ H_x \\ H_y
# \end{pmatrix}
# ```
# and ``\bm{T}`` is the ``4 \times 4`` matrix presented in
# [Clemmow and Heading (1954)](https://doi.org/10.1017/S030500410002939X).
# ``\bm{T}`` consists of elements of the susceptibility matrix for the ionosphere
# (and can be calculated with [`LongwaveModePropagator.tmatrix`](@ref)).
#
# To compute the wavefields, we calculate them at some height ``z``, use
# ``\mathrm{d}\bm{e}/\mathrm{d}z`` to step to a new height, and repeat from a point high
# in the ionosphere down to the ground.
# The initial solution comes from the Booker Quartic, which is a solution for the
# wavefields in a homogeneous ionosphere.
# At a great height in the ionosphere, we are interested in the two quartic solutions
# corresponding to upward-going waves.
# Therefore, we are integrating two sets of 4 complex variables simultaneously.
# Throughout the integration, the two solutions lose their independence because of
# numerical accuracy limitations over a wide range of field magnitudes.
# To maintain accuracy, the wavefields are orthonomalized repeatedly during the downward
# integration, and the scaling values are stored so that the fields can be "recorrected"
# after the integration is complete.
#
# Here are what the real component of the ``E_{x,1}`` and ``H_{x,2}`` wavefields
# looks like with and without the "recorrection".

zs = 110e3:-50:0
zskm = zs/1000

e = LMP.integratewavefields(zs, ea, frequency, bfield, day; unscale=true)
e_unscaled = LMP.integratewavefields(zs, ea, frequency, bfield, day; unscale=false)

ex1 = getindex.(e, 1)
ex1_unscaled = getindex.(e_unscaled, 1)
hx2 = getindex.(e, 7)
hx2_unscaled = getindex.(e_unscaled, 7)

p1 = plot(real(ex1), zskm; title="\$E_{x,1}\$",
          ylims=(0, 90), xlims=(-1.2, 1.2), linewidth=1.5, ylabel="altitude (km)",
          label="corrected", legend=:topleft);
plot!(p1, real(ex1_unscaled),
      zskm; linewidth=1.5, label="scaled only");

p2 = plot(real(hx2), zskm; title="\$H_{x,2}\$",
          ylims=(0, 90), xlims=(-1.2, 1.2), linewidth=1.5, legend=false);
plot!(p2, real(hx2_unscaled),
      zskm; linewidth=1.5);

plot(p1, p2; layout=(1,2))
#md savefig("wavefields_scaling.png"); nothing #hide
#md # ![](wavefields_scaling.png)
#

# ## Differential equations solver
#
# Pitteway (1965) used a basic form of a Runge-Kutta integrator with fixed step size.
# Julia has a fantastic [DifferentialEquations suite](https://diffeq.sciml.ai/stable/)
# for integrating many different forms of differential equations.
# Although wavefields are only computed two times for every transition in a
# [`SegmentedWaveguide`](@ref), we would still like to choose an efficient
# solver that requires
# relatively few function calls while still maintaining good accuracy.
#
# Let's try a few different solvers and compare their runtime for the
# day and night ionospheres.
#
# We can pass [`IntegrationParams`](@ref) through the [`LMPParams`](@ref) struct.
# Where necessary, we set the `lazy` interpolant option  of the solvers to `false`
# because we discontinuously scale the fields in the middle of the integration and the
# `lazy` option is not aware of this.

TO = TimerOutput()

zs = 110e3:-50:0
solvers = [Tsit5(), BS5(lazy=false), OwrenZen5(),
           Vern6(lazy=false), Vern7(lazy=false), Vern8(lazy=false), Vern9(lazy=false)]

solverstrings = replace.(string.(solvers), "OrdinaryDiffEq."=>"")

day_es = []
night_es = []
for s in eachindex(solvers)
    ip = IntegrationParams(solver=solvers[s], tolerance=1e-6)
    params = LMPParams(wavefieldintegrationparams=ip)

    ## make sure method is compiled
    LMP.integratewavefields(zs, ea, frequency, bfield, day; params=params);
    LMP.integratewavefields(zs, ea, frequency, bfield, night; params=params);

    solverstring = solverstrings[s]
    let day_e, night_e
        ## repeat 200 times to average calls
        for i = 1:200
            ## day ionosphere
            @timeit TO solverstring begin
                day_e = LMP.integratewavefields(zs, ea, frequency, bfield, day; params=params)
            end
            ## night ionosphere
            @timeit TO solverstring begin
                night_e = LMP.integratewavefields(zs, ea, frequency, bfield, night; params=params)
            end
        end
        push!(day_es, day_e)
        push!(night_es, night_e)
    end
end

# A quick plot to ensure none of the methods have problems.

day_e1s = [getindex.(e, 1) for e in day_es]

plot(real(day_e1s), zs/1000;
     label=permutedims(solverstrings), legend=:topleft)
#md savefig("wavefields_day.png"); nothing #hide
#md # ![](wavefields_day.png)

#
# And at night...

night_e1s = [getindex.(e, 1) for e in night_es]

plot(real(night_e1s), zs/1000;
     label=permutedims(solverstrings), legend=:topright)
#md savefig("wavefields_night.png"); nothing #hide
#md # ![](wavefields_night.png)

#
# The times to run each...

TO

# The `Tsit5`, `BS5`, and 6th, 7th, and 8th order Vern methods all have similar performance.
# In fact, rerunning these same tests multiple times can result in different solvers
# being "fastest".
#
# We use `Tsit5()`as the default in LongwaveModePropagator, which is similar to MATLAB's
# `ode45`, but usually more efficient.
#
# Why not a strict `RK4` or a more basic method? I tried them and they produce some
# discontinuities around the reflection height. It is admittedly difficult to tell
# if this is a true failure of the methods at this height or a problem related
# to the scaling and saving callbacks that occur during the integration.
# In any case, none of the methods tested above exhibit that issue.

# ## Pitteway figure 2
#
# Let's reproduce the wavefields in figure 2 of Pitteway (1965).

zs = 110e3:-500:50e3
zskm = zs/1000

e = LMP.integratewavefields(zs, ea, frequency, bfield, day; unscale=true)

ex1 = getindex.(e, 1)
ey1 = getindex.(e, 2)
ex2 = getindex.(e, 5)
hx2 = getindex.(e, 7)

function plotfield(field; kwargs...)
    p = plot(real(field), zskm; color="black", linewidth=1.5, legend=false,
             xlims=(-0.8, 0.8), label="real",
             framestyle=:grid; kwargs...)
    plot!(p, imag(field), zskm; color="black",
          linewidth=1.5, linestyle=:dash, label="imag")
    plot!(p, abs.(field), zskm; color="black", linewidth=3, label="abs")
    plot!(p, -abs.(field), zskm; color="black", linewidth=3, label="")
    return p
end

fs = 10
ex1p = plotfield(ex1; ylims=(49, 81), yaxis=false, yformatter=_->"");
annotate!(ex1p, 0.2, 79, text("\$E_{x,1}\$", fs));
ey1p = plotfield(ey1; ylims=(49, 81));
annotate!(ey1p, 0.2, 79, text("\$E_{y,1}\$", fs));
annotate!(ey1p, -0.98, 83, text("height (km)", fs));
ex2p = plotfield(ex2; ylims=(49, 86), yaxis=false, yformatter=_->"");
annotate!(ex2p, 0.3, 84, text("\$E_{x,2}\$", fs));
hx2p = plotfield(hx2; ylims=(49, 86));
annotate!(hx2p, 0.35, 84, text("\$H_{x,2}\$", fs, :center));
plot(ex1p, ey1p, ex2p, hx2p; layout=(2,2), size=(400,600), top_margin=5mm)
#md savefig("wavefields_fig2.png"); nothing #hide
#md # ![](wavefields_fig2.png)
#

# ```@raw html
# <img src="../../images/Pitteway1965_fig2.png"/>
# ```

# The envelopes of the two are very similar.
# The precise position of the real and imaginary wave components are not important
# because they each represent an instant in time and change based on the starting
# height of the integration.
#
# ## Pitteway figure 3
#
# Figure 3 of Pitteway (1965) uses a different scenario.
# A wave of frequency 202 kHz is at normal incidence through the nighttime ionosphere
# presented in Piggott (1965).
#
# First let's set up the new scenario.

ea = EigenAngle(0)
frequency = Frequency(202e3)
bfield = BField(50e-6, deg2rad(68), deg2rad(111))
ground = GROUND[8]
nothing  #hide

#  Now integrating the wavefields.

zs = 110e3:-50:70e3
zskm = zs/1000

e = LMP.integratewavefields(zs, ea, frequency, bfield, night)

ey1 = getindex.(e, 2)
hx2 = getindex.(e, 7)

ey1p = plotfield(ey1; ylims=(75, 102), title="\$E_{y,1}\$");
hx2p = plotfield(hx2; ylims=(75, 102), title="\$H_{x,2}\$");
plot(ey1p, hx2p; layout=(1,2), size=(400,500))
#md savefig("wavefields_fig3.png"); nothing #hide
#md # ![](wavefields_fig3.png)
#

# ```@raw html
# <img src="../../images/Pitteway1965_fig3.png"/>
# ```


text:
# mapping symbols to Human-readable label
activities_h = (
                (:work,         "Work"),
                (:sleep,        "Sleep"),
                (:tv,           "TV/Radio/Computer"),
                (:relax,        "Socalizing, Relaxing, Leisure"),
                (:household,    "Household chores"),
                (:sport,        "Sports, Exercise, Recreation"),
                (:edu,          "Education"),
                (:religion,     "Religious & Spiritual activities"),
                (:buy,          "Consumer Purchases"),
                (:travel,       "Travelling"),
                (:eat,          "Eating and drinking"),
                (:personalcare, "Personal Care"),
                (:carehh,       "Caring for household members"),
                (:carenhh,      "Caring for non HH members"),
                (:volunteer,    "Volunteering"),
                (:telephone,    "Telephone"));

# Formatter for y-axis
function timeformatter(x)
  if abs(x)<=60
    minutes=round(Int,x)
    return string(minutes)*"m"
  else
    hours, minutes = fldmod(round(Int,x),60)
    if minutes == 0
      return string(hours)*"h:00m"
    else
      return string(hours)*"h:"*string(minutes)*"m"
    end
  end
end

# Tick selection helper function
import PlotUtils
function optimize_minutes_ticks(mmin,mmax,k_min=2,k_max=4)
  msec_per_min=60000
  return PlotUtils.optimize_datetime_ticks(mmin*msec_per_min-1, mmax*msec_per_min, k_min=k_min,k_max=k_max)[1]/msec_per_min
end

# Annotation overlay
# usage:
# overlay = add_annotation_overlay!(p)
# annotate!(overlay, [(0.95, 0.15,"annotation")])
function add_annotation_overlay!(p)
  overlay=plot!(inset = (bbox(0.0,0.0,1,1, :bottom,:right)), ticks=nothing, subplot=length(p)+1, bg_inside=nothing,grid=false,border=nothing,leg=false)
  return p[end]
end

# Apply loess smoothing to (x,y)
import Loess
function loess_filter(x,y)
  model=Loess.loess(Float64.(x),Float64.(y))
  return Loess.predict(model,Float64.(x))
end

function create_plot_workvsfree(work_df, nwork_df)

    # a list of available fonts: PyPlot.matplotlib[:font_manager][:findSystemFonts]()
    # or in python:
    # import matplotlib.font_manager
    # set([f.name for f in matplotlib.font_manager.fontManager.ttflist])
    # I picked Arial as the closest Helvetica alternative available
    tickfont = Plots.font("Arial", 18)
    titlefont = Plots.font("Arial", 20)
    annotation_font = Plots.font("Arial", 16, :left)
    default(xtickfont=tickfont, ytickfont=tickfont, titlefont=titlefont)
    #gr()
    #4x4 plot layout plus an extended plot to create space for annotations
    l = @layout  [grid(4,4);b{0.01h} ]
    p=plot(layout=l,axis=true,size=(2400,2000),left_margin=50px, bottom_margin=50px,top_margin=50px, right_margin=50px)
    overlay=add_annotation_overlay!(p)

    #add actual plots
    for i in 1:size(activities_h,1)
      sym=activities_h[i][1]
      title_txt=activities_h[i][2]

      #non workdays
      x=convert(Array,nwork_df[:age])
      y=convert(Array,nwork_df[sym])
      y_smooth = loess_filter(x,y)
      scatter!(p[i], x, y, leg=false, label="", title=title_txt, color=:orange, marker=:rect, markersize=3,markerstrokewidth=0)
      plot!(p[i], x, y_smooth, leg=false, label="Non-workday", title=title_txt, lw=4, color=:orange)

      #workdays
      x=convert(Array,work_df[:age])
      y=convert(Array,work_df[sym])
      y_smooth = loess_filter(x,y)
      scatter!(p[i], x, y, leg=false, label="", title=title_txt, color=:green, marker=:circle, markersize=3,markerstrokewidth=0)
      plot!(p[i], x, y_smooth, leg=false, label="Workday", title=title_txt,ls=:dash, color=:green, lw=4)

      #adjust y-limits and optimize tick marks for minutes/hours
      ylims!(p[i], 0, Inf)
      expand_factor=1.2
      upper_y=ylims(p[i])[2]*expand_factor
      ylims!(p[i],0,upper_y)
      ticks=optimize_minutes_ticks(0,upper_y)
      yticks!(p[i], ticks, formatter=timeformatter)
      xticks!(p[i], 20:20:80 ,string.(20:20:80).*"yrs",grid=false)
      # this also applies to xaxis in PyPlot, probably a bug
      yaxis!(p[i],foreground_color_border=:white)
    end

    # put legend in a selected subplot
    plot!(p[4], leg=true,legendfont=tickfont)

    # add data source information in plot "17"
    plot!(p[17],axis=false)
    annotate!(overlay, [(0.75, 0.03, text("\nData source: American Time Use Survey (ATUS) 2003-2015\nBureau of Labor Statistics\nhttps://www.bls.gov/tus/datafiles_0315.htm\nWorkdays are days with more than 1 hour of work.",annotation_font))])

    # save file using Plots interface
    savefig("../figures/workfvsfree.pdf")

    # PyPlot fine-tuning by working on Plots.o
    p=current()

    # remove top and right ticksmarks, add x-spine, change titlefont weight to "bold"
    for i in 1:16
      ax=p.o[:axes][i]
      ax[:xaxis][:set_ticks_position]("bottom")
      ax[:yaxis][:set_ticks_position]("left")
      ax[:spines]["bottom"][:set_color]("grey")
      ax[:spines]["bottom"][:set_linewidth](1.5)
      title_text=ax[:get_title]()
      ax[:set_title](title_text, fontweight="bold", fontsize=titlefont.pointsize)
    end

    # save figure -- use PyPlot.savefig because we have modified the PyPlot object
    p.o[:set_size_inches](610/25.4, 508/25.4)
    PyPlot.savefig("../figures/workfvsfree.pdf")
end

function create_plot_diff(work_df, nwork_df)

  # a list of available fonts: PyPlot.matplotlib[:font_manager][:findSystemFonts]()
  # or in python:
  # import matplotlib.font_manager
  # set([f.name for f in matplotlib.font_manager.fontManager.ttflist])
  # I picked Arial as the closest Helvetica alternative available
  tickfont = Plots.font("Arial", 18)
  titlefont = Plots.font("Arial", 20)
  annotation_font = Plots.font("Arial", 16, :left)
  default(xtickfont=tickfont, ytickfont=tickfont, titlefont=titlefont)
  #gr()
  #4x4 plot layout plus an extended plot to create space for annotations
  l = @layout  [grid(4,4);b{0.01h} ]
  p=plot(layout=l,axis=true,size=(2400,2000),left_margin=50px, bottom_margin=50px,top_margin=50px, right_margin=50px)
  overlay=add_annotation_overlay!(p)

  #add actual plots
  for i in 1:size(activities_h,1)
    sym=activities_h[i][1]
    title_txt=activities_h[i][2]

    # don't plot work
    if sym==:work
      plot!(p[i],ticks=nothing, bg_inside=nothing,grid=false,border=nothing,leg=false)
      continue
    end;

    #non workdays
    x=convert(Array,nwork_df[:age])
    y1=convert(Array,nwork_df[sym])
    y1_smooth = loess_filter(x,y1)
    y2=convert(Array,work_df[sym])
    y2_smooth = loess_filter(x,y2)
    y=y1-y2
    y_smooth=y1_smooth-y2_smooth
    scatter!(p[i], x, y, leg=false, label="", title=title_txt, color=:orange, marker=:rect, markersize=3,markerstrokewidth=0)
    plot!(p[i], x, y_smooth, leg=false, label="additional time spent\non days off", title=title_txt, lw=4, color=:orange)

    #adjust y-limits and optimize tick marks for minutes/hours
    #ylims!(p[i], 0, Inf)
    expand_factor=1.2
    upper_y=ylims(p[i])[2]*expand_factor
    lower_y=ylims(p[i])[1]*expand_factor
    #ylims!(p[i],0,upper_y)
    ticks=optimize_minutes_ticks(lower_y,upper_y)
    yticks!(p[i], ticks, formatter=timeformatter)
    xticks!(p[i], 20:20:80 ,string.(20:20:80).*"yrs",grid=false)
    # this also applies to xaxis in PyPlot, probably a bug
    yaxis!(p[i],foreground_color_border=:white)
  end

  # put legend in a selected subplot
  plot!(p[4], leg=true,legendfont=tickfont)

  # add data source information in plot "17"
  plot!(p[17],axis=false)
  annotate!(overlay, [(0.75, 0.03, text("\nData source: American Time Use Survey (ATUS) 2003-2015\nBureau of Labor Statistics\nhttps://www.bls.gov/tus/datafiles_0315.htm\nWorkdays are days with more than 1 hour of work.",annotation_font))])

  # save file using Plots interface
  savefig("../figures/difference.pdf")

  # PyPlot fine-tuning by working on Plots.o
  p=current()

  # remove top and right ticksmarks, add x-spine, change titlefont weight to "bold"
  for i in 2:16
    ax=p.o[:axes][i]
    ax[:xaxis][:set_ticks_position]("bottom")
    ax[:yaxis][:set_ticks_position]("left")
    ax[:spines]["bottom"][:set_color]("grey")
    ax[:spines]["bottom"][:set_linewidth](1.5)
    title_text=ax[:get_title]()
    ax[:set_title](title_text, fontweight="bold", fontsize=titlefont.pointsize)
  end

  # save figure -- use PyPlot.savefig because we have modified the PyPlot object
  p.o[:set_size_inches](610/25.4, 508/25.4)
  PyPlot.savefig("../figures/difference.pdf")
end


text:
module FCFP
        export FCFPInstance

        struct FCFPInstance
                network_objective::Real
                all_edges::Array{Tuple}
                production_nodes::Array{}
                building_nodes::Array{}
                junctions_nodes::Array{}
                fixed_costs::Dict{Tuple, Real}
                heat_demand::Dict{}
                production::Dict{}
                capacities::Dict{Tuple, Real}
                old_buildings::Array{}
        end

        """
        Initialize the FCFP instance struct object with given arguments. All edges
        are specified within the `costs` dictionnary.
        """
        function initialize_data(network_objective, V, costs, demands, productions, capacities, old_buildings)
                production_nodes = [p for p=keys(productions)]
                building_nodes = [b for b=keys(demands)]
                junctions_nodes = setdiff(V, union(production_nodes, building_nodes))
                all_edges = collect(keys(costs))
                old_buildings = collect(old_buildings)
                inst = FCFPInstance(network_objective, all_edges, production_nodes, building_nodes, junctions_nodes, costs, demands, productions, capacities, old_buildings)
                return inst
        end

        FCFPInstance(network_objective, V, costs, demands, productions, capacities, old_buildings) = initialize_data(network_objective, V, costs, demands, productions, capacities, old_buildings)
end

max_stars_repo_name: Planheat/Planheat-Tool

max_issues_repo_name: Planheat/Planheat-Tool

max_forks_repo_name: Planheat/Planheat-Tool


text:
using Seismic

download("http://certmapper.cr.usgs.gov/nersl/NPRA/seismic/1979/616_79/PROCESSED/616_79_PR.SGY","616_79_PR.SGY")
SegyToSeis("616_79_PR.SGY","616_79_PR.seis")
d,h,e = SeisRead("616_79_PR.seis")
imx = Seismic.ExtractHeader(h,"imx")
println("size(d)=",size(d))
println("imx=",imx)
SeisWrite("616_79_PR_b",d,h,e)

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
mutable struct Rake
    stop_words::Vector{String}
    min_char_length::Int
    max_words_length::Int
    min_keyword_frequency::Int
    min_words_length_adj::Int
    max_words_length_adj::Int
    min_phrase_freq_adj::Int
end

struct Score
    rank::Float64
    frequency::Int
end

function Rake(stop_words; min_char_length=3, max_words_length=3, min_keyword_frequency=10,
              min_words_length_adj=1, max_words_length_adj=1, min_phrase_freq_adj=2)
    Rake(stop_words, min_char_length, max_words_length, min_keyword_frequency,
         min_words_length_adj, max_words_length_adj, min_phrase_freq_adj)
end

import Base.run

"Take vector of strings and split each string into all possible keyphrases"
function find_all_possible_keyphrases(phrases, max_phrase_length::Int)
    all_possible_keyphrases = Array(String,0)
    for sentence in phrases
        split_sentence = split.(sentence, ' ')
        for j = 0: max_phrase_length-1
            for i = 1:(length(split_sentence)-j)
                push!(all_possible_keyphrases, join(split_sentence[i:i+j], " "))
            end
        end
    end
    return all_possible_keyphrases
end

function run(text, stop_words, max_phrase_length)

    punct = ['.', '?', '!', '\n', ';', ',', ':', '\u2019', '\u2013']
    whitespace = false
    buffer = ""
    current_word = ""

    phrases = SubString{String}[]
    for i in text
        if i in punct # punctuation
            append!(phrases, find_all_possible_keyphrases([buffer], max_phrase_length))
            buffer = ""
            current_word = ""
            whitespace = false

        elseif (i == ' ') && (!whitespace) # hit whitespace, last character wasnt whitespace
            if current_word in stop_words # stopword
                append!(phrases, find_all_possible_keyphrases([buffer], max_phrase_length))
                buffer = ""
                current_word = ""
                whitespace = false
            else
                if buffer != ""
                    buffer = buffer*" "*string(current_word)
                else
                    buffer = buffer*string(current_word)
                end
                current_word = ""
                whitespace = true
            end

        elseif (i == ' ') && (whitespace) # hit whitspace, last character was whitespace
            continue

        else
            current_word = current_word*string(lowercase(i))
            whitespace = false
        end
    end

    if buffer != ""
        buffer = buffer*" "*string(current_word)
    else
        buffer = buffer*string(current_word)
    end
    print(buffer)
    append!(phrases, find_all_possible_keyphrases([buffer], max_phrase_length))

    keyphrases = unique(phrases)
    d = Dict([(i, count(x->x==i,phrases)) for i in phrases])
    return d
end


text:
using LaTeXStrings

function draw_axes()
    axvline(x=0, color="0.7", lw=2, zorder=1)
    axhline(y=0, color="0.7", lw=2, zorder=1)
end

"""
    viz_response(t, y,
                 plot_title="", plot_xlabel="time, t",
                 plot_ylabel="output, y(t)")

plot `y` vs. `t` to visualize the response of a system to an input. typically `t` and `y` are outputs of [`simulate`](@ref).

# Arguments
* `t::Array{Float64}`: array of times
* `y::Array{Float64}`: array of values of response variables at the corresponding times in `t`

# Example
```
julia> g = 4 / (4 * s ^ 2 + 0.8 * s + 1)
julia> u = 1 / s
julia> t, y = simulate(g * u, (0.0, 50.0))
julia> viz_response(t, y)
```
"""
function viz_response(t::Array{Float64}, y::Array{Float64};
        plot_title::Union{String, LaTeXString}="",
        plot_xlabel::Union{String, LaTeXString}=L"time, $t$",
        plot_ylabel::Union{String, LaTeXString}=L"output, $y(t)$",
    )

    figure()
    plot(t, y, zorder=100)
    xlabel(plot_xlabel)
    ylabel(plot_ylabel)
    title(plot_title)
    draw_axes()
    tight_layout()
    return nothing
end

"""
    viz_poles_and_zeros(tf)

plot the zeros and poles of the transfer function `tf` in the complex plane.
"""
function viz_poles_and_zeros(tf::TransferFunction)
    z, p, k = zeros_poles_gain(tf)

    figure()
    scatter(real.(z), imag.(z), marker="o", label="zeros", color="C1", zorder=100, s=50)
    scatter(real.(p), imag.(p), marker="x", label="poles", color="C2", s=50, zorder=100)
    legend()
    draw_axes()
    xlabel("Re")
    ylabel("Im")
    title("poles and zeros")
    tight_layout()
    return nothing
end

"""
    nyquist_diagram(tf)

plot the Nyquist diagram for a transfer function `tf` to visualize its frequency response.
"""
function nyquist_diagram(tf::TransferFunction; nb_pts::Int=300)
    ω = range(-10.0, 10.0, length=nb_pts)

    g_iω = [evaluate(tf, ω_i * im) for ω_i in ω]

    figure()
    plot(real(g_iω), imag(g_iω), zorder=100)
    draw_axes()
    xlabel("Re[G(iω)]")
    ylabel("Im[G(iω)]")
    title("Nyquist diagram")
    tight_layout()
    return nothing
end

"""
    root_locus(g_ol)

visualize the root locus plot of an open-loop transfer function `g_ol`.
"""
function root_locus(g_ol::TransferFunction)
    z, p, k = zeros_poles_k(g_ol)

    # keep gain of G_ol(s) +ve, scale with gain K
    Kcs = k * 10.0 .^ range(-6, 2, length=500)
    pushfirst!(Kcs, 0.0)

    rloc = zeros(Complex, length(Kcs), length(p))

    for (i_k, Kc) in enumerate(Kcs)
        c_poly = characteristic_polynomial(Kc * g_ol)
        roots_c_poly = roots(c_poly)
        for i_p = 1:length(roots_c_poly)
            rloc[i_k, i_p] = roots_c_poly[i_p]
        end
    end

    figure()
    for i = 1:length(p)
        plot(real.(rloc[:, i]), imag.(rloc[:, i]), zorder=100, color="C$(i-1)")
        scatter(real.([p[i]]), imag.([p[i]]), marker="x", label="poles", color="C$(i-1)", s=50, zorder=100)
    end
    xlabel("Re")
    ylabel("Im")
    draw_axes()
    title("root locus")
    tight_layout()
    return nothing
end

"""
    bode_plot(tf, log10_ω_min=-4.0, log10_ω_max=4.0)

draw the Bode plot of a transfer function `tf` to visualize its frequency response.
"""
function bode_plot(g::TransferFunction; log10_ω_min::Float64=-4.0, log10_ω_max::Float64=4.0)
    ω = 10.0 .^ range(log10_ω_min, log10_ω_max, length=300)
    g_iω = [evaluate(g, im * ω_i) for ω_i in ω]

    fig, (ax1, ax2) = subplots(2, 1, sharex=true, figsize=(8, 7))
    ax1.plot(ω, abs.(g_iω))
    ax1.set_ylabel(L"$|g(i\omega)|$")
    ax1.set_title("Bode plot")
    ax1.set_xscale("log")
    ax1.set_yscale("log")
    ax2.set_xscale("log")
    ax2.plot(ω, angle.(g_iω), color="C1")
    ax2.set_ylabel(L"$\angle g(i\omega)$")
    xlabel(L"frequency, $\omega$")
    tight_layout()
    return nothing
end


text:
using Test
using BenchmarkTools
using Plots
using SparseArrays

opts = TrajectoryOptimization.SolverOptions()
opts.verbose = false
opts.cost_tolerance = 1e-5
opts.cost_tolerance_intermediate = 1e-5
opts.constraint_tolerance = 1e-5
opts.square_root = false
opts.active_constraint_tolerance = 0.0
opts.outer_loop_update_type = :default
opts.penalty_max = 1e8
opts.live_plotting = false

# Parallel Park
model, = TrajectoryOptimization.Dynamics.dubinscar_parallelpark
# initial and goal states
x0 = [0.;0.;0.]
xf = [0.;1.;0.]

# costs
Q = (1e-2)*Diagonal(I,model.n)
Qf = 100.0*Diagonal(I,model.n)
R = (1e-2)*Diagonal(I,model.m)

# simulation
tf = 5.0
dt = 0.01

obj_uncon = LQRObjective(Q, R, Qf, tf, x0, xf)
obj = LQRObjective(Q, R, Qf, tf, x0, xf)

x_min = [-0.25; -0.001; -Inf]
x_max = [0.25; 1.001; Inf]

obj_con = TrajectoryOptimization.ConstrainedObjective(obj,x_min=x_min,x_max=x_max)

solver = TrajectoryOptimization.Solver(model,obj_con,integration=:rk4,N=31,opts=opts)
U = rand(solver.model.m, solver.N)

results, stats = TrajectoryOptimization.solve(solver,U)
λ_update_default!(results,solver)
update_constraints!(results,solver)
@assert max_violation(results) < opts.constraint_tolerance

J_prev = cost(solver,results)
c_max_prev = stats["c_max"][end]

p,pI,pE = get_num_constraints(solver)
p_N,pI_N,pE_N = get_num_terminal_constraints(solver)
## Newton 2 ###############
results_new = copy(results)

# newton_solve!(results_new,solver)

newton_results = NewtonResults(solver)
newton_active_set!(newton_results,results_new,solver)
# sum(newton_results.active_set)
# sum(vcat(results_new.active_set...))
# sum(newton_results.active_set_ineq)
# newton_results.s[findall(x->x != 0.0, newton_results.s)]
# findall(x->x < 0.0, vcat(results.C...)[newton_results.active_set_ineq])
update_newton_results!(newton_results,results_new,solver)
newton_step!(results_new,newton_results,solver,1.0)
max_violation(results_new)

a = 1


text:
# ---
# title: 433. Minimum Genetic Mutation
# id: problem433
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories:
# link: <https://leetcode.com/problems/minimum-genetic-mutation/description/>
# hidden: true
# ---
#
# A gene string can be represented by an 8-character long string, with choices
# from `"A"`, `"C"`, `"G"`, `"T"`.
#
# Suppose we need to investigate about a mutation (mutation from "start" to
# "end"), where ONE mutation is defined as ONE single character changed in the
# gene string.
#
# For example, `"AACCGGTT"` -> `"AACCGGTA"` is 1 mutation.
#
# Also, there is a given gene "bank", which records all the valid gene
# mutations. A gene must be in the bank to make it a valid gene string.
#
# Now, given 3 things - start, end, bank, your task is to determine what is the
# minimum number of mutations needed to mutate from "start" to "end". If there
# is no such a mutation, return -1.
#
# **Note:**
#
#   1. Starting point is assumed to be valid, so it might not be included in the bank.
#   2. If multiple mutations are needed, all mutations during in the sequence must be valid.
#   3. You may assume start and end string is not the same.
#
#
#
# **Example 1:**
#
#
#
#     start: "AACCGGTT"
#     end:   "AACCGGTA"
#     bank: ["AACCGGTA"]
#
#     return: 1
#
#
#
#
# **Example 2:**
#
#
#
#     start: "AACCGGTT"
#     end:   "AAACGGTA"
#     bank: ["AACCGGTA", "AACCGCTA", "AAACGGTA"]
#
#     return: 2
#
#
#
#
# **Example 3:**
#
#
#
#     start: "AAAAACCC"
#     end:   "AACCCCCC"
#     bank: ["AAAACCCC", "AAACCCCC", "AACCCCCC"]
#
#     return: 3
#
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
import Blosc: compress, decompress

function parallel_density_grid{T}(C::Matrix{T}, width::T, d::Int, σ::T, ϵ::T = 1e-4, num_workers::Int = nworkers())
    num_points = size(C, 2)

    num_points_per_worker = fill(fld(num_points, num_workers), num_workers)
    num_points_per_worker[1:mod(num_points, num_workers)] += 1

    start_indices = unshift!(cumsum(num_points_per_worker), 0) .+ 1
    pop!(start_indices)
    params = map(zip(start_indices, num_points_per_worker)) do r
        C[:, range(r[1], r[2])], width, d, σ, ϵ
    end

    #ρs_compressed = pmap(params, pids = collect(take(workers(), num_workers))) do p
    ρs_compressed = pmap(params) do p
        compress(density_grid(p...), level = 9, shuffle = true, itemsize = sizeof(T))
    end

    ρ = mapreduce(.+, zeros(T, d^3), ρs_compressed) do ρ_compressed
        decompress(T, ρ_compressed)
    end

    return reshape(ρ, d, d, d)
end

function density_grid{T}(C::Matrix{T}, width::T, d::Int, σ::T, ϵ::T= 1e-4)
    hw = width / 2
    ρ = zeros(T, d, d, d)

    num_points = size(C, 2)

    r_max_sq = max(0, -σ^2 * log(ϵ))
    r_max = sqrt(r_max_sq)

    b = floor(Int, d * r_max / width) + 1

    g = collect(linspace(-hw, hw, d))

    @inbounds for i in 1:num_points
        ix = floor(Int, d * (C[1, i] + hw) / width) + 1
        iy = floor(Int, d * (C[2, i] + hw) / width) + 1
        iz = floor(Int, d * (C[3, i] + hw) / width) + 1

        x_min = max(1, ix - b)
        y_min = max(1, iy - b)
        z_min = max(1, iz - b)

        x_max = min(d, ix + b)
        y_max = min(d, iy + b)
        z_max = min(d, iz + b)

        for jx in x_min:x_max, jy in y_min:y_max, jz in z_min:z_max
            r_sq = (C[1, i] - g[jx])^2 + (C[2, i] - g[jy])^2 + (C[3, i] - g[jz])^2

            if r_sq < r_max_sq
                ρ[jx, jy, jz] += exp(-r_sq / σ^2)
            end
        end
    end

    return ρ
end

function density_grid{T}(x::Array{T, 3}, y::Array{T, 3}, z::Array{T, 3}, C::Matrix{T}, σ::T, ϵ::T = 1e-4)
    D = size(x)
    ρ = zeros(T, D)

    num_points = size(C, 2)

    R = CartesianRange(size(ρ))
    I_min, I_max = first(R), last(R)

    x_spread = x[I_max] - x[I_min]
    y_spread = y[I_max] - y[I_min]
    z_spread = z[I_max] - z[I_min]

    r_max_sq = max(0, -σ^2 * log(ϵ))
    r_max = sqrt(r_max_sq)

    dx = floor(Int, D[1] * r_max / x_spread) + 1
    dy = floor(Int, D[2] * r_max / y_spread) + 1
    dz = floor(Int, D[3] * r_max / z_spread) + 1

    B = CartesianIndex((dx, dy, dz))

    @inbounds for i in 1:num_points
        ix = floor(Int, D[1] * (C[1, i] - x[I_min]) / x_spread) + 1
        iy = floor(Int, D[2] * (C[2, i] - y[I_min]) / y_spread) + 1
        iz = floor(Int, D[3] * (C[3, i] - z[I_min]) / z_spread) + 1
        I = CartesianIndex((ix, iy, iz))

        for J in CartesianRange(max(I_min, I - B), min(I_max, I + B))
            r_sq = (C[1, i] - x[J])^2 + (C[2, i] - y[J])^2 + (C[3, i] - z[J])^2
            if r_sq < r_max_sq
                ρ[J] += exp(-r_sq / σ^2)
            end
        end
    end

    return ρ
end

function density_grid_stats{T}(width::T, d::Int, σ::T, ϵ::T)
    r_max = sqrt(max(0, -σ^2 * log(ϵ)))
    b = floor(Int, r_max / (width / d)) + 1
    num_blocks = (2b + 1)^3

    println("Max radius: $(r_max)")
    println("Blocks per dimension: $b")
    println("Total blocks: $(num_blocks)")
    flush(STDOUT)
end


text:
# ## Stem plot

#md # ![FILE_NAME.png](images/FILE_NAME.png)

x = sin.(0:0.3:2pi)
plot(x, line=:stem, marker=:star, markersize=20)


text:
# This file was generated, do not modify it. # hide
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)
mcr = misclassification_rate(ŷ, y[test])
accuracy = 1 - mcr

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
using DataFrames,CSV,Pipe
using DocStringExtensions

#include("structs.jl")

"""
$(FUNCTIONNAME)(xdataset::Dataset,ydataset::Dataset,comps::Int64,incsamples::Array{Int64,1} = collect(1:size(xdataset.X)[1]))

    Calculates a PLS model
"""
function calcPLS(xdataset::Dataset,ydataset::Dataset,comps::Int64,incsamples::Array{Int64,1} = collect(1:size(xdataset.X)[1]))

    X::Array{Union{Missing, Float64},2} = xdataset.X[incsamples,:]
    Y::Array{Union{Missing, Float64},2} = ydataset.X[incsamples,:]

    xmask::BitArray{2} = xdataset.xmask[incsamples,:]
    ymask::BitArray{2} = ydataset.xmask[incsamples,:]

    hasmv = sum(.~xdataset.xmask) > 0

    n,k = size(X)
    T = Array{Float64,1}[]
    P = Array{Float64,1}[]
    C = Array{Float64,1}[]
    W = Array{Float64,1}[]
    U = Array{Float64,1}[]

    # initialize by setting zeroes in missing values
    X[.~xmask] .= 0
    Y[.~ymask] .= 0

    dCrit = 1e-12
    maxIter = 1000

    ssx_orig = sum(X.^2)

    for i in range(1,stop=comps)

        iter = 0
        diff = 1

        t = []
        p = []
        c = []
        w = []
        u = convert(Array{Float64,1},X[:,1])

        while diff > dCrit && iter < maxIter

            t,w,u,c,diff = nipals(X, Y, xdataset, u, xmask,hasmv)

            iter+=1
        end

        #println("Iter:$(iter) ConvValue:$(diff)")

        p = X't

        # if X has missing values
        if hasmv
            pcorr = ((t.^2)'*xmask)'
            p ./= pcorr

            p[pcorr.==0].=0
        else
            p /= t't
        end

        Xcomp = t*p'
        Xcomp[.~xmask] .= 0
        X -= Xcomp
        X[.~xmask] .= 0

        Ycomp = t*c'
        Ycomp[.~ymask] .= 0
        Y -= Ycomp
        Y[.~ymask] .= 0

        matrixType = Array{Float64,1}

        push!(T,convert(Array{Float64,1},t))
        push!(P,convert(Array{Float64,1},p))
        push!(C,convert(Array{Float64,1},c))
        push!(W,convert(Array{Float64,1},w))
        push!(U,convert(Array{Float64,1},u))

    end

    tocompindices = string.(range(1,stop=comps))

    PLS(
        DataFrame(hcat(T...),"t".*tocompindices),
        DataFrame(hcat(P...),"p".*tocompindices),
        DataFrame(hcat(C...),"c".*tocompindices),
        DataFrame(hcat(W...),"w".*tocompindices),
        DataFrame(hcat(U...),"u".*tocompindices)
    )
end

function numcomps(model::MultivariateModel)
    size(model.P,2)
end

function predictY(xdataset::Dataset, model::PLS, incsamples::Array{Int64,1} = collect(1:size(xdataset.X,1)); comps::Int64 = numcomps(model))

    P = convert(Array{Float64,2},model.W)
    W = convert(Array{Float64,2},model.P)
    C = convert(Array{Float64,2},model.C)

    predictY(xdataset,W,P,C,incsamples,comps=comps)
end

function predictY(xdataset::Dataset, W::Array{Float64,2}, P::Array{Float64,2}, C::Array{Float64,2}, incsamples::Array{Int64,1} = collect(1:size(xdataset.X,1));comps::Int64=numcomps(model))

    X::Array{Union{Missing, Float64},2} = xdataset.X[incsamples,:]
    xmask::BitArray{2} = xdataset.xmask[incsamples,:]

    Ycvs::Array{Array{Union{Missing, Float64},2}} = []

    for a in 1:comps
        tpred = X*W[:,a]

        # if X has missing values
        tpred ./= xmask*W[:,a].^2

        Ycv = tpred*C[:,a]'

        push!(Ycvs,Ycv)

        Xpred = tpred * P[:,a]'

        X -= Xpred

        X[.~xmask] .= 0
    end

    YpredVar = cat(Ycvs...,dims=3)
    Xres = X

    Xres[.~xmask] .= missing

    return (YpredVar=YpredVar,Xres=Xres)
end

function setmissing(model::T, dataset::Dataset) where T <: MultivariateModel
    model.Xres[.~dataset.xmask] .= missing
end

function crossvalidate(xdataset::Dataset, ydataset::Dataset, comps::Int64, cvgroups::Int64=7)

    cvgroups = createCVgroups(xdataset.X,cvgroups);

    Ycv = Array{Union{Missing, Float64}}(missing, (size(ydataset.X)...,comps))

    for (train,test) in cvgroups
        cvmodel = calcPLS(xdataset,ydataset,comps, train)

        Ycv[test,:,:] = predictY(xdataset,cvmodel,test,comps=comps).YpredVar
    end

    return calcQ2(Ycv,ydataset.X)
end

function calcQ2(Ycvs,y)
    Ycv_cum = cumsum(Ycvs,dims=3)

    R2cv_var = sum((y.-Ycv_cum).^2,dims=1)
    R2y_var = sum(y.^2,dims=1)

    Q2y_var = (R2y_var .- R2cv_var) ./ R2y_var

    return Q2y_var
end

function createCVgroups(df, numgroups::Int64)

    n,k = size(df)

    testsets = [collect(range(step,step=numgroups,stop=n)) for step in 1:numgroups]
    cvgroups = []

    for testset in testsets
        trainset = trues(n)
        trainset[testset] .= false
        push!(cvgroups,(findall(trainset),testset))
    end

    cvgroups
end

function nipals(X::Array{Union{Missing, Float64},2}, Y::Array{Union{Missing, Float64},2}, dataset::Dataset, u::Array{Float64,1}, xmask::BitArray{2},mv::Bool)

    w=X'*u

    # if X has missing values
    # valueCols=vec(sum(xmask,dims=1).>0)
    # novalueCols=.~valueCols
    if mv
        wcorr = ((u.^2)'*xmask)'
        w ./= wcorr

        w[wcorr.==0].=0
    end

    w/=norm(w)

    t =X*w

    # if X has missing values
    if mv
        t ./= xmask*w.^2
    end

    c = Y't/(t'*t)

    unew = Y*c / (c'*c)

    diff =  rootsq(u .- unew)/rootsq(u)

    t,w,unew,c,diff

end

function rootsq(array)
    array .^2 |> sum |>sqrt
end


text:
Base.@kwdef mutable struct GlobalOptimizerOptions{T<:AbstractFloat}

    # Subsolver Options
    "Subsolver used to solve mixed-integer linear programs"
    mip_optimizer::MOI.OptimizerWithAttributes = MOI.OptimizerWithAttributes(GLPK.Optimizer)
    "Subsolver used to solve second-order conic programs"
    socp_optimizer::MOI.OptimizerWithAttributes = MOI.OptimizerWithAttributes(Hypatia.Optimizer)
    "Semidefinite Programming Optimizer"
    semidefinite_optimizer::MOI.OptimizerWithAttributes = MOI.OptimizerWithAttributes(Hypatia.Optimizer)
    "Subsolver used to locally solve nonlinear programs"
    nlp_optimizer::MOI.OptimizerWithAttributes = MOI.OptimizerWithAttributes(Ipopt.Optimizer)
    "Subsolver used to locally solve mixed-integer nonlinear programs "
    minlp_optimizer::MOI.OptimizerWithAttributes = MOI.OptimizerWithAttributes(Pavito.Optimizer)

    # Extensions (set as user-specified option)
    enable_optimize_hook::Bool = false
    ext::Dict{Symbol, Any} = Dict{Symbol,Any}()
    ext_type::ExtensionType = DefaultExt()

    # Default optimizer parameters
    # Presolving options
    "Should EAGO attempt to remove type assert issues for user-defined functions (default = false)"
    presolve_scrubber_flag::Bool = false
    "Create and use DAG representations of user-defined function (default = false)."
    presolve_to_JuMP_flag::Bool = false
    "[FUTURE FEATURE, NOT CURRENTLY IMPLEMENTED] Apply the epigraph reformulation
    to the problem (default = false)."
    presolve_epigraph_flag::Bool = false
    "Rerranges the DAG using registered transformations (default = false)"
    presolve_flatten_flag::Bool = false

    # Conic reformulations
    "Attempt to bridge convex constraint to second order cone"
    conic_convert_quadratic::Bool = false

    # Iteration logging options
    "Turns logging on records global bounds, node count and run time. Additional
     options are available for recording information specific to subproblems (default = false)."
    log_on::Bool = false
    "Turns on logging of times and feasibility of subproblems (default = false)"
    log_subproblem_info::Bool = false
    "Log data every `log_interval` iterations (default = 1)."
    log_interval::Int = 1

    # Optimizer display options
    "The amount of information that should be printed to console while solving
    values range from 0 - 4: 0 is silent, 1 shows iteration summary statistics
    only, 2-4 show varying degrees of details about calculations within each
    iteration (default = 1)."
    verbosity::Int = 1
    "Display summary of iteration to console every `output_iterations` (default = 10)"
    output_iterations::Int = 1000
    "Display header for summary to console every `output_iterations` (default = 100)"
    header_iterations::Int = 10000

    # Node branching options
    "Convex coefficient used to select branch point. Branch point is given by
    `branch_cvx_factor*xmid + (1-branch_cvx_factor)*xsol` (default = 0.25)"
    branch_cvx_factor::T = 0.25
    "Minimum distance from bound to have branch point normalized by width of
    dimension to branch on (default = 0.15)"
    branch_offset::T = 0.15
    "Variables to branch on (default is all nonlinear)."
    branch_variable::Vector{Bool} = Bool[]
    "[FUTURE FEATURE, NOT CURRENTLY IMPLEMENTED] Number of times repeat node
    processing priorto branching (default = 4)."
    branch_max_repetitions::Int = 4
    "[FUTURE FEATURE, NOT CURRENTLY IMPLEMENTED] Volume ratio tolerance required
    to repeat processing the current node (default = 0.9)"
    branch_repetition_tol::T = 0.9

    # Termination limits
    "Maximum number of nodes (default = 1E-7)"
    node_limit::Int = 1*10^7
    "Maximum CPU time in seconds (default = 1000)"
    time_limit::T = 1000.0
    "Maximum number of iterations (default 3E6)"
    iteration_limit::Int = 3*10^6
    "Absolute tolerance for termination (default = 1E-3)"
    absolute_tolerance::T = 1E-3
    "Relative tolerance for termination (default = 1E-3)"
    relative_tolerance::T = 1E-3
    "Absolute constraint feasibility tolerance"
    absolute_constraint_feas_tolerance::T = 1E-6
    "Perform only a local solve of the problem (default = false)."
    local_solve_only::Bool = false
    "[TO BE REMOVED] Flag stops B&B loop if feasible point found."
    feasible_local_continue::Bool = false

    # Options for constraint propagation
    "Depth in B&B tree above which constraint propagation should be disabled (default = 1000)"
    cp_depth::Int = 20
    "Number of times to repeat forward-reverse pass routine (default = 3)"
    cp_repetitions::Int = 3
    "Disable constraint propagation if the ratio of new node volume to beginning node volume exceeds
    this number (default = 0.99)"
    cp_tolerance::T = 0.99
    "Use only valid interval bounds during constraint propagation (default = false)"
    cp_interval_only::Bool = false

    # obbt options
    "Depth in B&B tree above which OBBT should be disabled (default = 6)"
    obbt_depth::Int = 4
    "Number of repetitions of OBBT to perform in preprocessing (default = 3)"
    obbt_repetitions::Int = 20
    "Turn aggresive OBBT on (default = false)"
    obbt_aggressive_on::Bool = true
    "Maximum iteration to perform aggresive OBBT (default = 2)"
    obbt_aggressive_max_iteration::Int = 2
    "Minimum dimension to perform aggresive OBBT (default = 2)"
    obbt_aggressive_min_dimension::Int = 2
    "Tolerance to consider bounds equal (default = 1E-9)"
    obbt_tolerance::T = 1E-9
    obbt_variable_values::Vector{Bool} = Bool[]

    # Options for linear bound tightening
    "Depth in B&B tree above which linear FBBT should be disabled (default = 1000)"
    fbbt_lp_depth::Int  = 1000
    "Number of repetitions of linear FBBT to perform in preprocessing (default = 3)"
    fbbt_lp_repetitions::Int  = 3

    # Options for quadratic bound tightening
    "[FUTURE FEATURE, NOT CURRENTLY IMPLEMENTED] Depth in B&B tree above which univariate quadratic FBBT should be disabled (default = -1)"
    quad_uni_depth::Int = -1
    "[FUTURE FEATURE, NOT CURRENTLY IMPLEMENTED] Number of repetitions of univariate quadratic FBBT to perform in preprocessing (default = 2)"
    quad_uni_repetitions::Int = 2
    "[FUTURE FEATURE, NOT CURRENTLY IMPLEMENTED] Depth in B&B tree above which bivariate
    quadratic FBBT should be disabled (default = -1)"
    quad_bi_depth::Int = -1
    "[FUTURE FEATURE, NOT CURRENTLY IMPLEMENTED] Number of repetitions of bivariate quadratic FBBT to perform in preprocessing (default = 2)."
    quad_bi_repetitions::Int = 2

    # Duality-based bound tightening (DBBT) options
    "Depth in B&B tree above which duality-based bound tightening should be disabled (default = 1E10)"
    dbbt_depth::Int = 10^10
    "New bound is considered equal to the prior bound if within dbbt_tolerance (default = 1E-9)."
    dbbt_tolerance::T = 1E-8

    # Subgradient tightening flag
    "Relax Tag used to specify type of McCormick operator"
    relax_tag::RelaxTag = NS()
    "Perform tightening of interval bounds using subgradients at each factor in
    each nonlinear tape during a forward pass (default = true)."
    subgrad_tighten::Bool = true
    "Perform tightening of interval bounds using subgradients at each factor in
    each nonlinear tape during a reverse pass (default = false)."
    reverse_subgrad_tighten::Bool = false
    "Outer round computed subgradient bounds by this amount"
    subgrad_tol::T = 1E-10

    # Tolerance to add cuts and max number of cuts
    "Minimum number of cuts at each node to attempt (unsafe cuts not necessarily added)"
    cut_min_iterations::Int = 1
    "Maximum number of cuts at each node to attempt"
    cut_max_iterations::Int = 3
    "Convex coefficient used to select point for new added cuts. Branch point is
    given by `(1-cut_cvx)*xmid + cut_cvx*xsol` (default = 0.9)."
    cut_cvx::T = 0.9
    "Add cut if the L1 distance from the prior cutting point to the new cutting
    point normalized by the box volume is greater than the tolerance (default = 0.05)."
    cut_tolerance::T = 0.05

    "Use tolerances to determine safe cuts in a Khajavirad 2018 manner"
    cut_safe_on::Bool = true
    "Lower tolerance for safe-lp cut, Khajavirad 2018"
    cut_safe_l::T = 1E-8
    "Upper tolerance for safe-lp cut, Khajavirad 2018"
    cut_safe_u::T = 1E8
    "Constant tolerance for safe-lp cut, Khajavirad 2018"
    cut_safe_b::T = 1E9

    "Solve upper problem for every node with depth less than `upper_bounding_depth`
    and with a probabilityof (1/2)^(depth-upper_bounding_depth) otherwise (default = 6)"
    upper_bounding_depth::Int = 6

    # handling for domain violations
    "Amount about a domain violation to ignore when propagating bounds."
    domain_violation_guard_on::Bool = false
    "Amount about a domain violation to ignore when propagating bounds."
    domain_violation_ϵ::T = 1E-9
end


text:
using OptStoic
using Test

@testset "OptStoic.jl test suite" begin
    include("thermodynamics.jl")
    include("analysis.jl")
    include("build_model.jl")
end

max_stars_repo_name: VictorVeraFrazao/OptStoic

max_issues_repo_name: VictorVeraFrazao/OptStoic

max_forks_repo_name: VictorVeraFrazao/OptStoic


text:
const FlatField{P,T,M} = Union{FlatS0{P,T,M},FlatS2{P,T,M},FlatS02{P,T,M}}

const FlatFieldMap{P,T,M} = Union{FlatMap{P,T,M},FlatS2Map{P,T,M},FlatS02Map{P,T,M}}
const FlatFieldFourier{P,T,M} = Union{FlatFourier{P,T,M},FlatS2{P,T,M},FlatS02Fourier{P,T,M}}

### pretty printing
@show_datatype show_datatype(io::IO, t::Type{F}) where {N,θ,∂mode,D,T,M,F<:FlatField{Flat{N,θ,∂mode,D},T,M}} =
print(io, "$(pretty_type_name(F)){$(join(N .* (1,1), "×"))$(D==1 ? "" : "(×$D)") map, $(θ)′ pixels, $(∂mode.name.name), $M}")
for F in (:FlatMap, :FlatFourier,
          :FlatQUMap, :FlatQUFourier, :FlatEBMap, :FlatEBFourier,
          :FlatIQUMap, :FlatIQUFourier, :FlatIEBMap, :FlatIEBFourier)
    @eval pretty_type_name(::Type{<:$F}) = $(string(F))
end
function Base.summary(io::IO, f::FlatField{<:Flat{N,<:Any,<:Any,D}}) where {N,D}
    print(io, "$(length(f))", (D==1 ? "" : "(×$D)"), "-element ")
    showarg(io, f, true)
end

### field info
@memoize fieldinfo(::Type{P},::Type{T}=Float32,::Type{M}=Matrix) where {Nside,θpix,∂mode,D,P<:Flat{Nside,θpix,∂mode,D},T,M} =
    (;FlatInfo(T,basetype(M),θpix,Nside,D)..., ∂mode=∂mode)
@memoize fieldinfo(::Type{F}) where {P<:Flat,T,M,F<:FlatField{P,T,M}} =
    (;fieldinfo(P,T,M)..., @namedtuple(P,M,B=basis(F),S=spin(F))...)
fieldinfo(::F) where {F<:FlatField} = fieldinfo(F)

### promotion & conversion
# note: we don't need to promote the eltype T here since that will be
# automatically handled in broadcasting
function promote(f1::F1, f2::F2) where {T1,θ1,N1,∂mode1,F1<:FlatField{<:Flat{N1,θ1,∂mode1},T1},T2,θ2,N2,∂mode2,F2<:FlatField{<:Flat{θ2,N2,∂mode2},T2}}
    B     = promote_type(basis(F1),basis(F2))
    ∂mode = promote_type(∂mode1,∂mode2)
    B(∂mode(f1)), B(∂mode(f2))
end
(::Type{∂mode})(f::F) where {∂mode<:∂modes,N,θ,D,F<:FlatS0{<:Flat{N,θ,<:Any,D}}} = basetype(F){Flat{N,θ,∂mode,D}}(fieldvalues(f)...)
(::Type{∂mode})(f::FieldTuple{B}) where {∂mode<:∂modes,B} = FieldTuple{B}(map(∂mode,f.fs))

### basis-like definitions
LenseBasis(::Type{<:FlatS0})  =    Map
LenseBasis(::Type{<:FlatS2})  =  QUMap
LenseBasis(::Type{<:FlatS02}) = IQUMap
DerivBasis(::Type{<:FlatS0{<:Flat{<:Any,<:Any,fourier∂}}})  =    Fourier
DerivBasis(::Type{<:FlatS2{<:Flat{<:Any,<:Any,fourier∂}}})  =  QUFourier
DerivBasis(::Type{<:FlatS02{<:Flat{<:Any,<:Any,fourier∂}}}) = IQUFourier
DerivBasis(::Type{<:FlatS0{<:Flat{<:Any,<:Any,map∂}}})      =    Map
DerivBasis(::Type{<:FlatS2{<:Flat{<:Any,<:Any,map∂}}})      =  QUMap
DerivBasis(::Type{<:FlatS02{<:Flat{<:Any,<:Any,map∂}}})     = IQUMap

### derivatives

## Fourier-space
# the use of @generated here is for memoization so that we don't have do the
# `prefactor * im * ...` each time. actually, since for the first two cases
# these are only 1-d arrays, that's a completely unnecessary optimization, but
# without the @generated it triggers an order-dependenant compilation bug which
# e.g. slows down LenseFlow by a factor of ~4 so we gotta keep it for now ¯\_(ツ)_/¯
@generated function broadcastable(::Type{F}, ::∇diag{1,<:Any,prefactor}) where {P,T,M,prefactor,F<:FlatFourier{P,T,M}}
    @unpack kx = fieldinfo(F)
    basetype(M)(@. prefactor * im * kx')
end
@generated function broadcastable(::Type{F}, ::∇diag{2,<:Any,prefactor}) where {P,T,M,prefactor,F<:FlatFourier{P,T,M}}
    @unpack ky, Ny = fieldinfo(F)
    basetype(M)(@. prefactor * im * ky[1:Ny÷2+1])
end
@generated function broadcastable(::Type{F}, ::∇²diag) where {P,T,M,F<:FlatFourier{P,T,M}}
    @unpack kx, ky, Ny = fieldinfo(F)
    basetype(M)(@. kx'^2 + ky[1:Ny÷2+1]^2)
end

## Map-space
function copyto!(f′::F, bc::Broadcasted{<:Any,<:Any,typeof(*),Tuple{∇diag{coord,covariant,prefactor},F}}) where {coord,covariant,prefactor,T,θ,N,D,F<:FlatMap{Flat{N,θ,map∂,D},T}}
    D!=1 && error("Gradients of batched map∂ flat maps not implemented yet.")
    f = bc.args[2]
    n,m = size(f.Ix)
    α = 2 * prefactor * fieldinfo(f).Δx
    if coord==1
        @inbounds for j=2:m-1
            @simd for i=1:n
                f′[i,j] = (f[i,j+1] - f[i,j-1])/α
            end
        end
        @inbounds for i=1:n
            f′[i,1] = (f[i,2]-f[i,end])/α
            f′[i,end] = (f[i,1]-f[i,end-1])/α
        end
    elseif coord==2
        @inbounds for j=1:n
            @simd for i=2:m-1
                f′[i,j] = (f[i+1,j] - f[i-1,j])/α
            end
            f′[1,j] = (f[2,j]-f[end,j])/α
            f′[end,j] = (f[1,j]-f[end-1,j])/α
        end
    end
    f′
end

### bandpass
HarmonicBasis(::Type{<:FlatS0}) = Fourier
HarmonicBasis(::Type{<:FlatQU}) = QUFourier
HarmonicBasis(::Type{<:FlatEB}) = EBFourier
broadcastable(::Type{F}, bp::BandPass) where {P,T,F<:FlatFourier{P,T}} = Cℓ_to_2D(P,T,bp.Wℓ)

### logdets
logdet(L::Diagonal{<:Complex,<:FlatFourier}) = batch(real(sum_kbn(nan2zero.(log.(L.diag[:Il,full_plane=true])),dims=(1,2))))
logdet(L::Diagonal{<:Real,   <:FlatMap})     = batch(real(sum_kbn(nan2zero.(log.(complex.(L.diag.Ix))),dims=(1,2))))
### traces
tr(L::Diagonal{<:Complex,<:FlatFourier}) = batch(real(sum_kbn(L.diag[:Il,full_plane=true],dims=(1,2))))
tr(L::Diagonal{<:Real,   <:FlatMap})     = batch(real(sum_kbn(complex.(L.diag.Ix),dims=(1,2))))

### misc
Cℓ_to_Cov(f::FlatField{P,T,M}, args...; kwargs...) where {P,T,M} = adapt(M, Cℓ_to_Cov(P,T,spin(f),args...; kwargs...))

function pixwin(f::FlatField)
    @unpack θpix,P,T,k = fieldinfo(f)
    Diagonal(FlatFourier{P,T}((pixwin.(θpix,k) .* pixwin.(θpix,k'))[1:end÷2+1,:]))
end

global_rng_for(::Type{<:FlatField{<:Any,<:Any,M}}) where {M} = global_rng_for(M)

"""
    fixed_white_noise(rng, F)

Like white noise but the amplitudes are fixed to unity, only the phases are
random. Currently only implemented when F is a Fourier basis. Note that unlike
[`white_noise`](@ref), fixed white-noise generated in EB and QU Fourier bases
are not statistically the same.
"""
fixed_white_noise(rng, F::Type{<:FlatFieldFourier}) =
     exp.(im .* angle.(basis(F)(white_noise(rng,F)))) .* fieldinfo(F).Nside


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#    This function computes the sunlight, penumbra, and umbra times of the
#    satellite for one orbit every day in a year.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Longo, C. R. O., Rickman, S. L (1995). Method for the Calculation of
#       Spacecraft Umbra and Penumbra Shadow Terminator Points. NASA Technical
#       Paper 3547.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export satellite_eclipse_time

"""
    function satellite_eclipse_time(JD0::Number, a::Number, e::Number, i::Number, w::Number, RAAN::Number, numDays::Integer, relative::Bool = false)

Compute the eclipse time of a satellite.

# Args

* `JD0`: Initial instant for the analysis [Julian day].
* `a`: Semi-major axis of the orbit [m].
* `e`: Orbit eccentricity.
* `i`: Orbit inclination [rad].
* `w`: Argument of perigee [rad].
* `RAAN`: Right ascension of the ascending node at `JD0` [rad].
* `numDays`: Number of days of the analysis.
* `relative`: Compute the eclipse time relative to the nodal period.

# Returns

The following table:

        day | Sunlight Time | Penumbra Time | Umbra Time
       -----+---------------+---------------+------------

"""
function satellite_eclipse_time(JD0::Number,
                                a::Number,
                                e::Number,
                                i::Number,
                                w::Number,
                                RAAN::Number,
                                numDays::Integer,
                                relative::Bool = false)
    # Constants
    deg2rad = pi/180.0
    rad2deg = 180.0/pi
    day2sec = 24.0*60.0*60.0

    # Step of the orbit propagation (mean anomaly) [rad].
    step = 0.1*deg2rad

    # Initialization of variables.
    theta = 0.0                   # Sun angle relative to the inertial
                                  # coordinate frame.

    days = collect(0:1:numDays-1) # Vector of the days in which the beta angle
                                  # will be computed.

    # Mean anomaly.
    M = collect(0:step:2*pi)

    # Penumbra time.
    p_time = zeros(numDays)

    # Umbra time.
    u_time = zeros(numDays)

    # Sunlight time.
    s_time = zeros(numDays)

    # Semi-lactum rectum.
    p = a*(1.0-e^2)

    # Angular velocity of the orbit [rad/s].
    n = angvel(a, e, i, :J2)

    # Step in time
    tstep = step/n

    # Perturbations.
    #
    # RAAN rotation rate [rad/s].
    dOmega = dRAAN(a, e, i, :J2)

    # Perturbation of the argument of perigee [rad/s].
    dw = dArgPer(a, e, i, :J2)

    # Loop.
    for d in days
        # Get the sun position represented in the Inertial coordinate frame.
        s_i = sun_position_i(JD0+d)

        # Compute the new orbit parameters due to perturbations.
        w_d = w + dw*(d*day2sec)
        RAAN_d = RAAN + dOmega*(d*day2sec)

        for m_k in M
            # Get the satellite position vector represented in the Inertial coordinate
            # frame.
            f = M_to_f(e, m_k)

            (r_i, rt_i) = satellite_position_i(a, e, i, RAAN_d, w_d, f)

            # Check the lighting conditions.
            lighting = satellite_lighting_condition(r_i, s_i)

            if (lighting == SAT_LIGHTING_SUNLIGHT)
                s_time[d+1] += tstep
            elseif (lighting == SAT_LIGHTING_UMBRA)
                u_time[d+1] += tstep
            else
                p_time[d+1] += tstep
            end
        end
    end

    if (!relative)
        [days s_time p_time u_time]
    else
        total_time = s_time + p_time + u_time
        [days s_time./total_time p_time./total_time u_time./total_time]
    end
end


text:
# Command: Back
"""
	back!(session::Session)

This command causes the browser to traverse one step backward in the joint session history of the current top-level browsing context.
This is equivalent to pressing the back button in the browser chrome or invoking window.history.back.
"""
function back!(session::Session)
    @unpack addr, id = session
    response = HTTP.post("$addr/session/$id/back", [("Content-Type" => "application/json")])
    @assert response.status == 200
    nothing
end


text:
# Temporarily add this repo's Maple to the module path.
# Do not use this if you installed Maple using Pkg.clone.
push!(LOAD_PATH, "$(pwd())/src")

using Maple

# '0' and ' ' will be confidered "nothing"
lvl_1_fg = """
333333333333333333333333333333333333333333333333
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
333333333333333333333333333333333333333333333333
"""

fgTiles = FgTiles(lvl_1_fg)

println("Map type")
@time map = Map(
    "Testing",
    Room[
        Room(
            name = "lvl_1",

            fgTiles = fgTiles,

            position = (0, 0),
            size = size(fgTiles),

            entities = Entity[
                Player(16, 144),
                ChapterCheckpoint(16, 144),

                CustomMemorial(160, 160, true, "memorial_test", "decals/playerR"),
            ],
        )
    ]
)

println("Dict of map $(map.package)")
@time dMap = Dict(map)

println("Total encode time")
@time encodeMap(dMap, "Testing.bin")

max_stars_repo_name: CelestialCartographers/Maple

max_issues_repo_name: CelestialCartographers/Maple

max_forks_repo_name: CelestialCartographers/Maple


text:
#!/usr/bin/env julia
using LazySets, Reachability, MathematicalSystems,
      LinearAlgebra, SparseArrays

# fix namespace conflicts with MathematicalSystems
using LazySets: LinearMap

# common aliases
const CLCCS = ConstrainedLinearControlContinuousSystem
const CLCDS = ConstrainedLinearControlDiscreteSystem
const LCS = LinearContinuousSystem
const LDS = LinearDiscreteSystem

# compatibility between Julia versions
include("../src/compat.jl")
using Compat.Test

# Note: CDDLib is only required for v0.6; for v0.7 or later the default Polyhedra
# library is used and CDDLib can be removed from REQUIRE

include("Systems/alltests.jl")
include("ReachSets/alltests.jl")
include("Reachability/alltests.jl")
include("Properties/alltests.jl")


text:
# Get displacement at a given iteration
function get_fields( fname::String , itt::Int)

    #Get Variables
    fid         = h5open( fname , "r")

    max_iter = read_attribute(fid, "max_iter")

    xmin = read_attribute(fid, "xmin")
    xmax = read_attribute(fid, "xmax")
    Nx   = read_attribute(fid, "xnodes")

    ymin = read_attribute(fid, "ymin")
    ymax = read_attribute(fid, "ymax")
    Ny   = read_attribute(fid, "ynodes")

    #Discretize space
    dx  = (xmax - xmin)/(Nx-1)
    x   = [xmin + (i-1)*dx for i in 1:Nx]

    dy  = (ymax - ymin)/(Ny-1)
    y   = [ymin + (i-1)*dy for i in 1:Ny]

    #Check it time iteration is OK
    if( itt > max_iter)
        println("ERROR: out of bounds (max iteration is" , max_iter,")")
        return
    end

    #Get Displacement
    g1      = open_group(fid, "ψ")

    dname   = string("iteration_" , itt)
    dset    = open_dataset(g1, dname )
    ttt     = read_attribute(dset,"time")

    ψ = reshape(Array(dset) , Nx,Ny )

    #Get Derivative
    g1      = open_group(fid, "dψ")

    dname   = string("iteration_" , itt)
    dset    = open_dataset(g1, dname )

    dψ = reshape(Array(dset),Nx,Ny)

    close(fid)

    #Return all data
    return ttt , x , y , ψ , dψ

end