text:
module GrassmannAverages

function grassmann_median(X, K)
    K >= 0 || throw(ArgumentError("Number of components must be positive"));

    N = size(X, 1);
    D = size(X, 2);
    eps = 1e-4;  "Too small of error is problematic"

    basis = Array{Real}(D, K);
    for k=1:K
        "q represents the k-th principal component"
        q = randn(D);
        q /= norm(q);

        while true
            dotsign = sign(X * q);

            "Compute the weighted average of the elements. Robustness is achieved by element-wise median."
            newQ = Array{Real}(D);
            for j=1:D
                newQ[j] = median(X[: , j] .* dotsign);
            end
            newQ /= norm(newQ);

            "If the L_1 difference between q and newQ is small, then terminate."
            max_diff = 0;
            for j=1:D
                max_diff = max(max_diff, abs(q[j] - newQ[j]));
            end
            if max_diff < eps
                break
            end

            q = newQ;
        end

        "Store the basis vector and subtract from X"
        basis[: , k] = q;
        X -= (X * q) * q';
    end

    return basis;
end

end # module

max_stars_repo_name: lvnguyen/GrassmannAverages

max_issues_repo_name: lvnguyen/GrassmannAverages

max_forks_repo_name: lvnguyen/GrassmannAverages


text:
function add_human_name_trait(T, d)
    ex = quote
        MMI.human_name(::Type{<:$T}) = $d
    end
    eval(ex)
    return
end

function meta(T; input=Unknown, target=Unknown, output=Unknown,
              weights::Bool=false, human_name::String="")
    ex = quote
        MMI.input_scitype(::Type{<:$T})    = $input
        MMI.output_scitype(::Type{<:$T})   = $output
        MMI.target_scitype(::Type{<:$T})   = $target
        MMI.supports_weights(::Type{<:$T}) = $weights
        isempty($human_name) || (MMI.human_name(::Type{<:$T}) = $human_name)
    end
    eval(ex)
    return
end

ext: jl


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

"""
    ValueAccessor{S<:AbstractValueShape}

A value accessor provides a means to access a value with a given shape
stored in a flat real-valued data vector with a given offset position.

Constructor:

    ValueAccessor{S}(shape::S, offset::Int)

The offset is relative to the first index of a flat data array, so if
the value is stored at the beginning of the array, the offset will be zero.

An `ValueAccessor` can be used to index into a given flat data array.

Example:

```julia
acc = ValueAccessor(ArrayShape{Real}(2,3), 2)
valshape(acc) == ArrayShape{Real,2}((2, 3))
data = [1, 2, 3, 4, 5, 6, 7, 8, 9]
data[acc] == [3 5 7; 4 6 8]
```
"""
struct ValueAccessor{S<:AbstractValueShape}
    shape::S
    offset::Int
    len::Int

    ValueAccessor{S}(shape::S, offset::Int) where {S<:AbstractValueShape} =
        new{S}(shape, offset, totalndof(shape))
end

export ValueAccessor

ValueAccessor(shape::S, offset::Int) where {S<:AbstractValueShape} = ValueAccessor{S}(shape, offset)

# Value accessors behave as scalars under broadcasting:
@inline Base.Broadcast.broadcastable(shape::ValueAccessor) = Ref(shape)

nonabstract_eltype(va::ValueAccessor) = nonabstract_eltype(va.shape)

Base.size(va::ValueAccessor) = size(va.shape)
Base.length(va::ValueAccessor) = length(va.shape)

valshape(va::ValueAccessor) = va.shape

# Would this be useful?
# AbstractValueShape(va::ValueAccessor) = valshape(va)
# Base.convert(::Type{AbstractValueShape}, va::ValueAccessor) = AbstractValueShape(va)

function view_range end

@inline function view_range(idxs::AbstractUnitRange{<:Integer}, va::ValueAccessor)
    from = first(idxs) + va.offset
    to = from + va.len - 1
    from:to
end

function view_idxs end

@inline view_idxs(idxs::AbstractUnitRange{<:Integer}, va::ValueAccessor) = view_range(idxs, va)

Base.@propagate_inbounds function Base.view(
    data::AbstractMatrix, va_row::ValueAccessor, va_col::ValueAccessor
)
    view(data, view_idxs(axes(data, 1), va_row), view_idxs(axes(data, 2), va_col))
end


text:
using Bridge
using StaticArrays
import Bridge: b, σ, B, β, a, constdiff
const ℝ = SVector{N,T} where {N,T}
import Base.display

"""
    FitzhughDiffusion <: ContinuousTimeProcess{ℝ{2}}

Struct defining FitzHugh-Nagumo model
"""
struct FitzhughDiffusion{T,TP} <: ContinuousTimeProcess{ℝ{2,T}}
    param::TP
    ϵ::T
    s::T
    γ::T
    β::T
    σ::T

    function FitzhughDiffusion(ϵ::T, s::T, γ::T, β::T, σ::T) where T
        TP = Val{:regular}
        new{T,TP}(TP(), ϵ, s, γ, β, σ)
    end

    function FitzhughDiffusion(::Val{S}, ϵ::T, s::T, γ::T, β::T, σ::T
                               ) where {T,S}
        TP = Val{S}
        new{T,TP}(TP(), ϵ, s, γ, β, σ)
    end

    function FitzhughDiffusion(sym::Symbol, ϵ::T, s::T, γ::T, β::T, σ::T
                               ) where T
        checkParamValid(sym)
        TP = Val{sym}
        new{T,TP}(TP(), ϵ, s, γ, β, σ)
    end

    function FitzhughDiffusion(param::String, ϵ::T, s::T, γ::T, β::T, σ::T
                               ) where T
        TP = stringToParam(param)
        new{T,TP}(TP(), ϵ, s, γ, β, σ)
    end
end

function checkParamValid(s::Symbol)
    sym = [:regular, :simpleAlter, :complexAlter, :simpleConjug, :complexConjug]
    (s in sym) || error("Invalid parametrisation of the FitzHugh-Nagumo model")
end

function stringToParam(s::String)
    s = lowercase(s)
    if s == "regular"
        param = Val{:regular}
    elseif s in ["simplealter", "simple alter", "simple alternative"]
        param = Val{:simpleAlter}
    elseif s in ["complexalter", "complex alter", "complex alternative"]
        param = Val{:complexAlter}
    elseif s in ["simpleconjug", "simple conjug", "simple conjugate"]
        param = Val{:simpleConjug}
    elseif s in ["complexconjug", "complex conjug", "complex conjugate"]
        param = Val{:complexConjug}
    else
        error("Invalid parametrisation of the FitzHugh-Nagumo model")
    end
    param
end

# REGULAR PARAMETRISATION
# -----------------------

function b(t, x, P::FitzhughDiffusion{T,Val{:regular}}) where T
    ℝ{2}((x[1]-x[2]-x[1]^3+P.s)/P.ϵ, P.γ*x[1]-x[2] + P.β)
end

function σ(t, x, P::FitzhughDiffusion{T,Val{:regular}}) where T
    ℝ{2}(0.0, P.σ)
end

# ALTERNATIVE PARAMETRISATION
# ---------------------------

function b(t, x, P::FitzhughDiffusion{T,TP}
           ) where {T,TP <: Union{Val{:simpleAlter},Val{:complexAlter}}}
    ℝ{2}(x[2], -( (P.γ-1.0)*x[1] + x[1]^3 + P.ϵ*x[2] - P.s + P.β
                  + (3.0*x[1]^2 - 1.0)*x[2])/P.ϵ )
end

function σ(t, x, P::FitzhughDiffusion{T,TP}
           ) where {T,TP <: Union{Val{:simpleAlter},Val{:complexAlter}}}
    ℝ{2}(0.0, P.σ/P.ϵ)
end

# CONJUGATE PARAMETRISATION
# -------------------------

function b(t, x, P::FitzhughDiffusion{T,TP}
           ) where {T,TP <: Union{Val{:simpleConjug},Val{:complexConjug}}}
    ℝ{2}(x[2], ((P.ϵ - P.γ)*x[1] - P.ϵ*(x[1]^3 + (3.0*x[1]^2 - 1.0)*x[2])
                + P.s - P.β - x[2]))
end

function σ(t, x, P::FitzhughDiffusion{T,TP}
           ) where {T,TP <: Union{Val{:simpleConjug},Val{:complexConjug}}}
    ℝ{2}(0.0, P.σ)
end

# APPLICABLE TO ALL PARAMETRISATIONS
# ----------------------------------

constdiff(::FitzhughDiffusion) = true
clone(P::FitzhughDiffusion, θ) = FitzhughDiffusion(P.param, θ...)
params(P::FitzhughDiffusion) = [P.ϵ, P.s, P.γ, P.β, P.σ]

"""
    struct FitzhughDiffusionAux <: ContinuousTimeProcess{ℝ{2}}

Struct defining proposal diffusion (proposal for sampling from FitzHugh-Nagumo
diffusion)
"""
struct FitzhughDiffusionAux{R,S1,S2,TP} <: ContinuousTimeProcess{ℝ{2,R}}
    param::TP
    ϵ::R
    s::R
    γ::R
    β::R
    σ::R
    t::Float64
    u::S1
    T::Float64
    v::S2

    function FitzhughDiffusionAux(ϵ::R, s::R, γ::R, β::R, σ::R, t::Float64,
                                  u::S1, T::Float64, v::S2) where {R,S1,S2}
        TP = Val{:regular}
        new{R,S1,S2,TP}(TP(), ϵ, s, γ, β, σ, t, u, T, v)
    end

    function FitzhughDiffusionAux(::Val{S}, ϵ::R, s::R, γ::R, β::R, σ::R,
                                  t::Float64, u::S1, T::Float64, v::S2
                                  ) where {R,S,S1,S2}
        TP = Val{S}
        new{R,S1,S2,TP}(TP(), ϵ, s, γ, β, σ, t, u, T, v)
    end

    function FitzhughDiffusionAux(sym::Symbol, ϵ::R, s::R, γ::R, β::R, σ::R,
                                  t::Float64, u::S1, T::Float64, v::S2
                                  ) where {R,S1,S2}
        checkParamValid(sym)
        TP = Val{sym}
        new{R,S1,S2,TP}(TP(), ϵ, s, γ, β, σ, t, u, T, v)
    end

    function FitzhughDiffusionAux(param::String, ϵ::R, s::R, γ::R, β::R, σ::R,
                                  t::Float64, u::S1, T::Float64, v::S2
                                  ) where {R,S1,S2}
        TP = stringToParam(param)
        new{R,S1,S2,TP}(TP(), ϵ, s, γ, β, σ, t, u, T, v)
    end
end

"""
    dependsOnParams(::FitzhughDiffusionAux)

Declare which parameters (1=>`ϵ`, 2=>`s`, 3=>`γ`, 4=>`β`, 5=>`σ`) the
auxiliary diffusion depends upon. Used for finding out which parameter
update requires also updating the values of the grid of `H`'s and `r`'s.
"""
function dependsOnParams end

# REGULAR PARAMETRISATION
# -----------------------

function B(t, P::FitzhughDiffusionAux{T,S1,S2,Val{:regular}}) where {T,S1,S2}
    @SMatrix [1/P.ϵ-3*P.v^2/P.ϵ  -1/P.ϵ; P.γ -1.0]
end

function β(t, P::FitzhughDiffusionAux{T,S1,S2,Val{:regular}}) where {T,S1,S2}
    ℝ{2}(P.s/P.ϵ+2*P.v^3/P.ϵ, P.β)
end

function σ(t, P::FitzhughDiffusionAux{T,S1,S2,Val{:regular}}) where {T,S1,S2}
    ℝ{2}(0.0, P.σ)
end

function dependsOnParams(::FitzhughDiffusionAux{T,S1,S2,Val{:regular}}) where {T,S1,S2}
    (1, 2, 3, 4, 5)
end

# SIMPLE ALTERNATIVE PARAMETRISATION
# ----------------------------------

function B(t, P::FitzhughDiffusionAux{T,S1,S2,Val{:simpleAlter}}) where {T,S1,S2}
    @SMatrix [0.0  1.0; 0.0 0.0]
end

function β(t, P::FitzhughDiffusionAux{T,S1,S2,Val{:simpleAlter}}) where {T,S1,S2}
    ℝ{2}(0.0, 0.0)
end

function σ(t, P::FitzhughDiffusionAux{T,S1,S2,Val{:simpleAlter}}) where {T,S1,S2}
    ℝ{2}(0.0, P.σ/P.ϵ)
end

function dependsOnParams(::FitzhughDiffusionAux{T,S1,S2,Val{:simpleAlter}}) where {T,S1,S2}
    (1, 5)
end

# COMPLEX ALTERNATIVE PARAMETRISATION
# -----------------------------------

function B(t, P::FitzhughDiffusionAux{T,S1,S2,Val{:complexAlter}}) where {T,S1,S2}
    @SMatrix [0.0  1.0; (1.0-P.γ-3.0*P.v[1]^2)/P.ϵ (1.0-P.ϵ-3.0*P.v[1]^2)/P.ϵ]
end

function β(t, P::FitzhughDiffusionAux{T,S1,S2,Val{:complexAlter}}) where {T,S1,S2}
    ℝ{2}(0.0, (2*P.v[1]^3+P.s-P.β)/P.ϵ)
end

function σ(t, P::FitzhughDiffusionAux{T,S1,S2,Val{:complexAlter}}) where {T,S1,S2}
    ℝ{2}(0.0, P.σ/P.ϵ)
end

function dependsOnParams(::FitzhughDiffusionAux{T,S1,S2,Val{:complexAlter}}) where {T,S1,S2}
    (1, 2, 3, 4, 5)
end

function B(t, P::FitzhughDiffusionAux{T,S1,SArray{Tuple{2},Float64,1,2},Val{:complexAlter}}) where {T,S1}
    @SMatrix [0.0  1.0;
              (1.0-P.γ-3.0*P.v[1]^2-6*P.v[1]*P.v[2])/P.ϵ (1.0-P.ϵ-3.0*P.v[1]^2)/P.ϵ]
end

function β(t, P::FitzhughDiffusionAux{T,S1,SArray{Tuple{2},Float64,1,2},Val{:complexAlter}}) where {T,S1}
    ℝ{2}(0.0, (2*P.v[1]^3+P.s-P.β+6*P.v[1]^2*P.v[2])/P.ϵ)#check later
end

# SIMPLE CONJUGATE PARAMETRISATION
# --------------------------------

function B(t, P::FitzhughDiffusionAux{T,S1,S2,Val{:simpleConjug}}) where {T,S1,S2}
    @SMatrix [0.0  1.0; 0.0 0.0]
end

function β(t, P::FitzhughDiffusionAux{T,S1,S2,Val{:simpleConjug}}) where {T,S1,S2}
    ℝ{2}(0.0, 0.0)
end

function σ(t, P::FitzhughDiffusionAux{T,S1,S2,Val{:simpleConjug}}) where {T,S1,S2}
    ℝ{2}(0.0, P.σ)
end

function dependsOnParams(::FitzhughDiffusionAux{T,S1,S2,Val{:simpleConjug}}) where {T,S1,S2}
    (5,)
end

# COMPLEX CONJUGATE PARAMETRISATION
# ---------------------------------
function B(t, P::FitzhughDiffusionAux{T,S1,S2,Val{:complexConjug}}) where {T,S1,S2}
    @SMatrix [0.0  1.0; (P.ϵ-P.γ-3.0*P.ϵ*P.v[1]^2) (P.ϵ-1.0-3.0*P.ϵ*P.v[1]^2)]
end

function β(t, P::FitzhughDiffusionAux{T,S1,S2,Val{:complexConjug}}) where {T,S1,S2}
    ℝ{2}(0.0, 2*P.ϵ*P.v[1]^3+P.s-P.β) # check later
end

function σ(t, P::FitzhughDiffusionAux{T,S1,S2,Val{:complexConjug}}) where {T,S1,S2}
    ℝ{2}(0.0, P.σ)
end

function dependsOnParams(::FitzhughDiffusionAux{T,S1,S2,Val{:complexConjug}}) where {T,S1,S2}
    (1, 2, 3, 4, 5)
end

function B(t, P::FitzhughDiffusionAux{T,S1,SArray{Tuple{2},Float64,1,2},Val{:complexConjug}}) where {T,S1}
    @SMatrix [0.0  1.0;
              (P.ϵ-P.γ-3.0*P.ϵ*P.v[1]^2-6*P.ϵ*P.v[1]*P.v[2]) (P.ϵ-1.0-3.0*P.ϵ*P.v[1]^2)]
end

function β(t, P::FitzhughDiffusionAux{T,S1,SArray{Tuple{2},Float64,1,2},Val{:complexConjug}}) where {T,S1}
    ℝ{2}(0.0, 2*P.ϵ*P.v[1]^3+P.s-P.β+6*P.ϵ*P.v[1]^2*P.v[2]) # check later
end

# APPLICABLE TO ALL PARAMETRISATIONS
# ----------------------------------

constdiff(::FitzhughDiffusionAux) = true
b(t, x, P::FitzhughDiffusionAux) = B(t,P) * x + β(t,P)
a(t, P::FitzhughDiffusionAux) = σ(t,P) * σ(t, P)'

"""
    clone(P::FitzhughDiffusionAux, θ)

Clone the object `P`, but use a different vector of parameters `θ`.
"""
clone(P::FitzhughDiffusionAux, θ) = FitzhughDiffusionAux(P.param, θ..., P.t,
                                                         P.u, P.T, P.v)
# should copy starting point or sth, currently restricted by the same type of u and v
clone(P::FitzhughDiffusionAux, θ, v) = FitzhughDiffusionAux(P.param, θ..., P.t,
                                                            zero(v), P.T, v)
params(P::FitzhughDiffusionAux) = [P.ϵ, P.s, P.γ, P.β, P.σ]

"""
    regularToAlter(x, ϵ, offset=0)
Transform point from observation under :regular parametrisation to the one under
:alter(...) parametrisation
"""
function regularToAlter(x, ϵ, offset=0)
    ℝ{2}(x[1], (x[1] - x[1]^3 - x[2] + offset) / ϵ)
end

"""
    alterToRegular(x, ϵ, offset=0)
Transform point from observation under :alter(...) parametrisation to the one
under :regular parametrisation
"""
function alterToRegular(x, ϵ, offset=0)
    ℝ{2}(x[1], x[1] - x[1]^3 - x[2]*ϵ + offset)
end

"""
    regularToConjug(x, ϵ, offset=0)
Transform point from observation under :regular parametrisation to the one under
:conjug(...) parametrisation
"""
function regularToConjug(x, ϵ, offset=0)
    ℝ{2}(x[1], (x[1] - x[1]^3 - x[2] + offset) * ϵ)
end

"""
    conjugToRegular(x, ϵ, offset=0)
Transform point from observation under :conjug(...) parametrisation to the one
under :regular parametrisation
"""
function conjugToRegular(x, ϵ, offset=0)
    ℝ{2}(x[1], x[1] - x[1]^3 - x[2]/ϵ + offset)
end

_MSG = ["\n---------------------------------------------------\n",
        "The target diffusion model is set to FitzHugh-Nagumo diffusion, a solution ",
        " to an SDE:\n",
        "dYₜ = (Yₜ-Yₜ³-Xₜ+s)/ϵ dt,\ndXₜ = (γYₜ-Xₜ+β) dt + σ dWₜ,  t∈[0,T].\n",
        "dYₜ = Ẏₜ dt,\ndẎₜ = [(1-γ)Yₜ-Yₜ³-ϵẎₜ+s-β+(1-3Yₜ²)Ẏₜ] dt + σ/ϵ dWₜ,  t∈[0,T].\n",
        "dYₜ = Ẏₜ dt,\ndẎₜ = [(ϵ-γ)Yₜ-ϵYₜ³-Ẏₜ+s-β+ϵ(1-3Yₜ²)Ẏₜ] dt + σ dWₜ,  t∈[0,T].\n",
        "The auxiliary diffusion is set to ",
        ", a solution to an SDE:\n",
        "dỸₜ = [(1-3y²)Ỹₜ-X̃ₜ+s+2y³] dt,\ndX̃ₜ = (γỸₜ-X̃ₜ+β) dt + σ dWₜ,  t∈[0,T],\n",
        "dIₜ = Bₜ dt,\ndBₜ = σ/ϵ dWₜ,  t∈[0,T],\n",
        "dỸₜ = X̃ₜ dt,\ndX̃ₜ = [(1-γ-3y²-6yẏ)Ỹₜ+(1-ϵ-3y²)X̃ₜ+(2y³+s-β+6y²ẏ)]/ϵ dt + σ/ϵ dWₜ,  t∈[0,T],\n",
        "dIₜ = Bₜ dt,\ndBₜ = σ dWₜ,  t∈[0,T],\n",
        "dỸₜ = X̃ₜ dt,\ndX̃ₜ = {[ϵ(1-3y²-6yẏ)-γ]Ỹₜ+[ϵ(1-3y²)-1]X̃ₜ+[ϵ(2y³+6y²ẏ)+s-β]} dt + σ dWₜ,  t∈[0,T],\n",
        "where (y,ẏ) denotes an end-point of (Y,Ẏ).\n",
        "---------------------------------------------------\n",
        ]

function display(::FitzhughDiffusion{T,Val{:regular}}) where T
    print(_MSG[1], _MSG[2], "(Y,X)", _MSG[3], _MSG[4], _MSG[15])
end

function display(::FitzhughDiffusion{T,Val{:simpleAlter}}) where T
    print(_MSG[1], _MSG[2], "(Y,Ẏ)", _MSG[3], _MSG[5], _MSG[15])
end

function display(::FitzhughDiffusion{T,Val{:complexAlter}}) where T
    print(_MSG[1], _MSG[2], "(Y,Ẏ)", _MSG[3], _MSG[5], _MSG[15])
end

function display(::FitzhughDiffusion{T,Val{:simpleConjug}}) where T
    print(_MSG[1], _MSG[2], "(Y,Ẏ)", _MSG[3], _MSG[6], _MSG[15])
end

function display(::FitzhughDiffusion{T,Val{:complexConjug}}) where T
    print(_MSG[1], _MSG[2], "(Y,Ẏ)", _MSG[3], _MSG[6], _MSG[15])
end

function display(::FitzhughDiffusionAux{T,S1,S2,Val{:regular}}) where {T,S1,S2}
    print(_MSG[1], _MSG[7], "(Ỹ,X̃)", _MSG[8], _MSG[9], _MSG[14], _MSG[15])
end

function display(::FitzhughDiffusionAux{T,S1,S2,Val{:simpleAlter}}) where {T,S1,S2}
    print(_MSG[1], _MSG[7], "(I,B)", _MSG[8], _MSG[10], _MSG[15])
end

function display(::FitzhughDiffusionAux{T,S1,S2,Val{:complexAlter}}) where {T,S1,S2}
    print(_MSG[1], _MSG[7], "(Ỹ,X̃)", _MSG[8], _MSG[11], _MSG[14], _MSG[15])
end

function display(::FitzhughDiffusionAux{T,S1,S2,Val{:simpleConjug}}) where {T,S1,S2}
    print(_MSG[1], _MSG[7], "(I,B)", _MSG[8], _MSG[12], _MSG[15])
end

function display(::FitzhughDiffusionAux{T,S1,S2,Val{:complexConjug}}) where {T,S1,S2}
    print(_MSG[1], _MSG[7], "(Ỹ,X̃)", _MSG[8], _MSG[13], _MSG[14], _MSG[15])
end


text:
struct AssetModel
    X0
    sigma
end

# stochastic process interface

function stateSize(self::AssetModel)   # dimension of X(t)
    return 1
end

function factors(self::AssetModel)   # dimension of W(t)
    return 1
end

function initialValues(self::AssetModel)
    return [ 0.0 ]
end

# calculate volatility taking into account exogenous hybrid adjuster
function volatility(self::AssetModel, t, Y)
    return self.sigma + Y[3]
end

function evolve(self::AssetModel, t0, Y0, dt, dW, Y1)
    sigma = volatility(self,t0,Y0)
    Y1[1] = Y0[1] + (Y0[2] - 0.5*sigma*sigma)*dt + sigma*dW[1]*sqrt(dt)
end

function asset(self::AssetModel, t, Y, alias)
    return self.X0 * exp(Y[1])
end

max_stars_repo_name: sschlenkrich/HybridMonteCarlo

max_issues_repo_name: sschlenkrich/HybridMonteCarlo

max_forks_repo_name: sschlenkrich/HybridMonteCarlo


text:
import UncertainData: RandomSequences

x = rand(100)
y = rand(100)

test_ccm = ConvergentCrossMappingTest(timeseries_lengths = [45, 50], n_reps = 20)
test_cm = CrossMappingTest(n_reps = 10)
test_vf = VisitationFrequencyTest(binning = RectangularBinning(5), ηs = 1:5)
test_tog = TransferOperatorGridTest(binning = RectangularBinning(5), ηs = 1:5)
test_jdd = JointDistanceDistributionTest()
test_jddt = JointDistanceDistributionTTest()
predtest = VisitationFrequencyTest(binning = RectangularBinning(5), ηs = -5:5)
test_pa = PredictiveAsymmetryTest(predictive_test = predtest)
test_npa = NormalisedPredictiveAsymmetryTest(predtest, f = 1.0)
test_NNmi = NearestNeighbourMITest(ηs = -5:5)

Nseq = 80
nchunks = 3
rtest_cm = RandomSequencesTest(test_cm, RandomSequences(nchunks, Nseq))
rtest_ccm = RandomSequencesTest(test_ccm, RandomSequences(nchunks, Nseq))
rtest_vf = RandomSequencesTest(test_vf, RandomSequences(nchunks, Nseq))
rtest_tog = RandomSequencesTest(test_tog, RandomSequences(nchunks, Nseq))
rtest_jdd = RandomSequencesTest(test_jdd, RandomSequences(nchunks, Nseq))
rtest_jddt = RandomSequencesTest(test_jddt, RandomSequences(nchunks, Nseq))
rtest_pa = RandomSequencesTest(test_pa, RandomSequences(nchunks, Nseq))
rtest_npa = RandomSequencesTest(test_npa, RandomSequences(nchunks, Nseq))
rtest_NNmi = RandomSequencesTest(test_NNmi, RandomSequences(nchunks, Nseq))

@test rtest_cm isa RandomSequencesTest
@test rtest_ccm isa RandomSequencesTest
@test rtest_vf isa RandomSequencesTest
@test rtest_tog isa RandomSequencesTest
@test rtest_jdd isa RandomSequencesTest
@test rtest_jddt isa RandomSequencesTest
@test rtest_pa isa RandomSequencesTest
@test rtest_npa isa RandomSequencesTest
@test rtest_NNmi isa RandomSequencesTest

@test causality(rtest_cm, x, y) isa CausalAnalysis
@test causality(rtest_ccm, x, y) isa CausalAnalysis
@test causality(rtest_vf, x, y) isa CausalAnalysis
@test causality(rtest_tog, x, y) isa CausalAnalysis
@test causality(rtest_jdd, x, y) isa CausalAnalysis
@test causality(rtest_jddt, x, y) isa CausalAnalysis
@test causality(rtest_pa, x, y) isa CausalAnalysis
@test causality(rtest_npa, x, y) isa CausalAnalysis
@test causality(rtest_NNmi, x, y) isa CausalAnalysis


text:
using Cairo
c  = CairoRGBSurface(512, 128);
cr = CairoContext(c);
save(cr);
set_source_rgb(cr, 0.8, 0.8, 0.8);  	# light gray
rectangle(cr, 0.0, 0.0, 512.0, 128.0); 	# background
fill(cr);
restore(cr);
save(cr);

x0=51.2;  y0=64.0;
x1=204.8; y1=115.4;
x2=307.2; y2=12.8;
x3=460.8; y3=64.0;

move_to (cr, x0, y0);
curve_to (cr, x1, y1, x2, y2, x3, y3);
set_line_width (cr, 10.0);
stroke (cr);
restore(cr);

move_to(cr, 12.0, 12.0);
set_source_rgb (cr, 0, 0, 0);
show_text(cr,"Stroke");
write_to_png(c,"stroke.png");


text:
using wembed
using Base.Test

@test isfile("../res/GoogleNews-vectors-negative300.bin")

bin_path = "../res/GoogleNews-vectors-negative300.bin"
model_path = "../res/word2vec_model.hdf5"
@time @test wembed.bin2hdf5(bin_path, model_path)

@time @test wembed.init(model_path)
@time @test typeof(wembed.getvector("queen", false)) == Vector{Float32}
@time @test typeof(getsimilarity("king", "queen")) == Float32
@time @test typeof(wordnearest("king", 5)) == Array{Pair{AbstractString,Float32},1}

wordsAsString = "breakfast sushi dinner lunch"
@time @test typeof(odd1out(wordsAsString)) == String

wordset1 = ["king", "man"]
wordset2 = ["queen", "woman"]
@time @test typeof(nsimilarity(wordset1, wordset2)) == Float32

words = map(String, split("Julia is a developing language"))
getvector(words, fill(0.25,300))


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule xtrx_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("xtrx")
JLLWrappers.@generate_main_file("xtrx", UUID("f7ce82bc-bc1e-5125-bd1d-6a84dd39bfb4"))
end  # module xtrx_jll


text:
using StaticArrays: SVector

struct Cube{T}
    origin::SVector{3,T}
    shape::Tuple{Int,Int,Int}
    cell::PeriodicCell{T}
    density::Array{T,3}
end

function Cube(filename, ::Type{T}=Float64) where {T<:AbstractFloat}

    data = readdlm(filename; skipstart=2)
    natoms = convert(Int, data[1,1])
    origin = SVector{3,T}(data[1,2:4])
    shape = Tuple(data[2:4])
    vectors_rows = data[2:4,2:4] .* shape
    vectors = permutedims(vectors_rows, (2,1))
    cell = PeriodicCell{T}(vectors, [true, true, true])

    volumetric_data = permutedims(data[5+natoms:end,:], [2,1])
    indices = volumetric_data .!= ""
    density = convert(Vector{T}, volumetric_data[indices])
    density = permutedims(reshape(density, reverse(shape)), 3:-1:1)

    Cube(origin, shape, cell, density)
end

function (cube::Cube)(r::AbstractVector)
    r = r - cube.origin
    r = cube.cell.inverse * r

    indices = floor.(Int, r .* cube.shape) .+ 1

    if any(indices .< 1) || any(indices .> cube.shape)
        throw(DomainError(r, "cannot evaluate the density outside of the cell."))
    end

    return cube.density[indices...]
end


text:
push!(LOAD_PATH, "./src")
using DelimitedFiles
using Plots
using SetupDG
using UnPack
using UniformTriMesh
using PlotThemes
using ToeplitzMatrices
using Basis2DTri
using Basis1D
using LinearAlgebra

using WriteVTK
using Test
const FloatType = Float64
const vtk_filename_noext = "unstructured"

S_N(x) = @. sin(pi*x/h)/(2*pi/h)/tan(x/2)
function vandermonde_Sinc(h,r)
    N = convert(Int, 2*pi/h)
    V = zeros(length(r),N)
    for n = 1:N
        V[:,n] = S_N(r.-n*h)
    end
    V[1,1] = 1
    V[end,end] = 1
    return V
end

T = 1.0
N = 4
Nq_P = countlines("src/QuadratureData/quad_nodes_tri_N$(2*N).txt")
Nd = 5
K1D = 128
K = K1D*K1D*2
Np_F = 8
Nplot = 2
Nplot_f = 250
Mach = ".3"
delta = ".05"
extra_attribute = "varwidthsquared"

folder_name = "K$(K1D)NPF$(Np_F)N$(N)delta$(delta)Mach$(Mach)$(extra_attribute)"
file_name = "array_2"
Q = readdlm("./examples/ES_Fourier-DG/$folder_name/$file_name.txt",'\t',Float64,'\n')
xq = readdlm("./examples/ES_Fourier-DG/$folder_name/xq.txt",'\t',Float64,'\n')
yq = readdlm("./examples/ES_Fourier-DG/$folder_name/yq.txt",'\t',Float64,'\n')

xq = reshape(xq,Nq_P,Np_F,K)
yq = reshape(yq,Nq_P,Np_F,K)

h       = 2*pi/Np_F
column  = [0; .5*(-1).^(1:Np_F-1).*cot.((1:Np_F-1)*h/2)]
Dt      = Array{Float64,2}(Toeplitz(column,column[[1;Np_F:-1:2]]))

rd = init_reference_tri(N);
@unpack fv,Nfaces,r,s,VDM,V1,Vp,Dr,Ds,rf,sf,wf,nrJ,nsJ,rq,sq,wq,Vq,M,Pq,Vf,LIFT = rd
rp, sp = Basis2DTri.equi_nodes_2D(Nplot)
Vp = Basis2DTri.vandermonde_2D(N,rp,sp)/VDM
Np = length(rp)

"Mesh related variables"
# Initialize 2D triangular mesh
VX,VY,EToV = uniform_tri_mesh(K1D,K1D)
@. VX = VX
@. VY = VY
md   = init_mesh((VX,VY),EToV,rd)
# Intialize triangular prism
VX   = repeat(VX,2)
VY   = repeat(VY,2)
VZ   = [2/Np_F*ones((K1D+1)*(K1D+1),1); ones((K1D+1)*(K1D+1),1)]
@unpack rxJ,sxJ,ryJ,syJ = md

JP = 1/K1D^2
JF = 1/pi
J = JF*JP
Dr = JF*Dr
Ds = JF*Ds
Dt = JP*Dt
VpPq = Vp*Pq

Q = reshape(Q,Nq_P,Np_F,Nd,K)
u = reshape(Q[:,:,2,:],Nq_P,Np_F*K)
v = reshape(Q[:,:,3,:],Nq_P,Np_F*K)
w = reshape(Q[:,:,4,:],Nq_P,Np_F*K)
norm_v = @. u^2+v^2+w^2

dudr = Vq*Dr*Pq*u
dvdr = Vq*Dr*Pq*v
dwdr = Vq*Dr*Pq*w
duds = Vq*Ds*Pq*u
dvds = Vq*Ds*Pq*v
dwds = Vq*Ds*Pq*w

dudt = zeros(size(dudr))
dvdt = zeros(size(dvdr))
dwdt = zeros(size(dwdr))

for k = 1:K
   f_idx = 1+(k-1)*Np_F:k*Np_F
   dudt[:,f_idx] = u[:,f_idx]*Dt'
   dvdt[:,f_idx] = v[:,f_idx]*Dt'
   dwdt[:,f_idx] = w[:,f_idx]*Dt'
end

q_criterion = zeros(Nq_P,Np_F*K)
tz = 1/JF
for k = 1:K
   rx,ry,sx,sy = rxJ[1,k],ryJ[1,k],sxJ[1,k],syJ[1,k]
   for f = 1:Np_F
      for q = 1:Nq_P
         idx = f+(k-1)*Np_F
         v_tensor = [dudr[q,idx]*rx+duds[q,idx]*sx dudr[q,idx]*ry+duds[q,idx]*sy dudt[q,idx]*tz;
                     dvdr[q,idx]*rx+dvds[q,idx]*sx dvdr[q,idx]*ry+dvds[q,idx]*sy dvdt[q,idx]*tz;
                     dwdr[q,idx]*rx+dwds[q,idx]*sx dwdr[q,idx]*ry+dwds[q,idx]*sy dwdt[q,idx]*tz;]
         s_tensor = 1/2*(v_tensor+transpose(v_tensor))
         Ω_tensor = 1/2*(v_tensor-transpose(v_tensor))
         q_criterion[q,idx] = 1/2*(norm(Ω_tensor,2)^2-norm(s_tensor,2)^2)
      end
   end
end

q_criterion = q_criterion/maximum(abs.(q_criterion))
@show maximum(abs.(q_criterion))
q_criterion = reshape(q_criterion,Nq_P,Np_F,K)
zq = reshape(repeat(collect(-1+2/Np_F:(2/Np_F):1),inner=(1,Nq_P),outer=(K,1))',Nq_P,Np_F,K)
zp = LinRange(-1+2/Np_F,1,Nplot_f)
tp = LinRange(h,2*pi,Nplot_f)
VDM_F = vandermonde_Sinc(h,tp)
V2 = vandermonde_1D(1,tp)/vandermonde_1D(1,LinRange(h,2*pi,Np_F))
p'p
x_plot = zeros(size(Vp,1),Nplot_f,K)
y_plot = zeros(size(Vp,1),Nplot_f,K)
z_plot = zeros(size(Vp,1),Nplot_f,K)
q_criterion_plot = zeros(size(Vp,1),Nplot_f,K)
norm_v = reshape(norm_v,Nq_P,Np_F,K)
norm_v_plot = zeros(size(Vp,1),Nplot_f,K)
for k = 1:K
   x_plot[:,:,k] = Vp*Pq*xq[:,:,k]*V2'
   y_plot[:,:,k] = Vp*Pq*yq[:,:,k]*V2'
   z_plot[:,:,k] = Vp*Pq*zq[:,:,k]*V2'
   q_criterion_plot[:,:,k] = Vp*Pq*q_criterion[:,:,k]*VDM_F'
   norm_v_plot[:,:,k] = Vp*Pq*norm_v[:,:,k]*VDM_F'
end

@show maximum(abs.(q_criterion_plot))

# plotting for 2 linear-quadratic wedge elements
Np = 12
Nplot_f = Nplot_f - 1
pts = Array{FloatType}(undef,3,Np*Nplot_f*K)
order_v = [1;6;3]
order_m = [4;5;2]
h_f = 1.75/Nplot_f

qcriterion_data = Array{FloatType}(undef,Np*Nplot_f*K)
norm_v_data = Array{FloatType}(undef,Np*Nplot_f*K)
for k = 1:K
    for f = 1:Nplot_f
        trivx = x_plot[order_v,f,k]
        trivy = y_plot[order_v,f,k]
        trimx = x_plot[order_m,f,k]
        trimy = y_plot[order_m,f,k]
        triz = z_plot[order_v,f,k]
        qv = q_criterion_plot[order_v,f,k]
        qm = q_criterion_plot[order_m,f,k]
        qv_next = q_criterion_plot[order_v,f+1,k]
        qm_next = q_criterion_plot[order_m,f+1,k]

        normv = norm_v_plot[order_v,f,k]
        normm = norm_v_plot[order_m,f,k]
        normv_next = norm_v_plot[order_v,f+1,k]
        normm_next = norm_v_plot[order_m,f+1,k]
        pts[1,(k-1)*(Np*Nplot_f)+(f-1)*Np+1:(k-1)*(Np*Nplot_f)+f*Np] = [trivx trivx trimx trimx]
        pts[2,(k-1)*(Np*Nplot_f)+(f-1)*Np+1:(k-1)*(Np*Nplot_f)+f*Np] = [trivy trivy trimy trimy]
        pts[3,(k-1)*(Np*Nplot_f)+(f-1)*Np+1:(k-1)*(Np*Nplot_f)+f*Np] = [triz triz.+h_f triz triz.+h_f]
        qcriterion_data[(k-1)*(Np*Nplot_f)+(f-1)*Np+1:(k-1)*(Np*Nplot_f)+f*Np] = [qv; qv_next; qm; qm_next]
        norm_v_data[(k-1)*(Np*Nplot_f)+(f-1)*Np+1:(k-1)*(Np*Nplot_f)+f*Np] = [normv; normv_next; normm; normm_next]
    end
end

# Create cells
celltype = VTKCellTypes.VTK_QUADRATIC_LINEAR_WEDGE
cells = MeshCell[]
@test VTKCellType(celltype.vtk_id) === celltype
@test_throws BoundsError   VTKCellType(-3)   # ids start at 0
@test_throws ArgumentError VTKCellType(200)  # 200 is an unknown cell type
inds = 1:Np
for i = 1:Nplot_f*K
    c = MeshCell(celltype,inds.+(i-1)*Np)
    push!(cells, c)
end

# Initialise new vtu file (unstructured grid).
fname = "$(vtk_filename_noext)_3D_$(Mach)_T=$(T)"
outfile = vtk_grid(fname, pts, cells, compress=3) do vtk
    vtk["q-criterion"] = qcriterion_data
    vtk["norm-velocity"] = norm_v_data
end

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
using StatsBase
using Test, Random

Random.seed!(1234)

n = 100000

# test that rng specification is working correctly
# a) if the same rng is passed to a sample function twice,
#    the results should be the same (repeatability)
# b) not specifying a rng should be the same as specifying Random.GLOBAL_RNG
function test_rng_use(func, non_rng_args...)
    # some sampling methods mutate a passed array and return it
    # so that the tests don't pass trivially, we need to copy those
    # pre-allocated storage arrays

    # repeatability
    @test func(MersenneTwister(1), deepcopy(non_rng_args)...) ==
          func(MersenneTwister(1), deepcopy(non_rng_args)...)
    # default RNG is Random.GLOBAL_RNG
    Random.seed!(47)
    x = func(deepcopy(non_rng_args)...)
    Random.seed!(47)
    y = func(Random.GLOBAL_RNG, deepcopy(non_rng_args)...)
    @test x == y
end

#### sample with replacement

function check_sample_wrep(a::AbstractArray, vrgn, ptol::Real; ordered::Bool=false)
    vmin, vmax = vrgn
    (amin, amax) = extrema(a)
    @test vmin <= amin <= amax <= vmax
    n = vmax - vmin + 1
    p0 = fill(1/n, n)
    if ordered
        @test issorted(a)
        if ptol > 0
            @test isapprox(proportions(a, vmin:vmax), p0, atol=ptol)
        end
    else
        @test !issorted(a)
        ncols = size(a,2)
        if ncols == 1
            @test isapprox(proportions(a, vmin:vmax), p0, atol=ptol)
        else
            for j = 1:ncols
                aj = view(a, :, j)
                @test isapprox(proportions(aj, vmin:vmax), p0, atol=ptol)
            end
        end
    end
end

import StatsBase: direct_sample!

a = direct_sample!(1:10, zeros(Int, n, 3))
check_sample_wrep(a, (1, 10), 5.0e-3; ordered=false)

a = direct_sample!(3:12, zeros(Int, n, 3))
check_sample_wrep(a, (3, 12), 5.0e-3; ordered=false)

a = direct_sample!([11:20;], zeros(Int, n, 3))
check_sample_wrep(a, (11, 20), 5.0e-3; ordered=false)

test_rng_use(direct_sample!, 1:10, zeros(Int, 6))

a = sample(3:12, n)
check_sample_wrep(a, (3, 12), 5.0e-3; ordered=false)

a = sample(3:12, n; ordered=true)
check_sample_wrep(a, (3, 12), 5.0e-3; ordered=true)

a = sample(3:12, 10; ordered=true)
check_sample_wrep(a, (3, 12), 0; ordered=true)

test_rng_use(sample, 1:10, 10)

@testset "sampling pairs" begin

    rng = Random.MersenneTwister(1)

    @test samplepair(rng, 2)  ===  (1, 2)
    @test samplepair(rng, 10) === (8, 2)

    @test samplepair(rng, [3, 4, 2, 6, 8]) === (2, 6)
    @test samplepair(rng, [1, 2])          === (1, 2)
end

test_rng_use(samplepair, 1000)

#### sample without replacement

function check_sample_norep(a::AbstractArray, vrgn, ptol::Real; ordered::Bool=false)
    # each column of a for one run

    vmin, vmax = vrgn
    (amin, amax) = extrema(a)
    @test vmin <= amin <= amax <= vmax
    n = vmax - vmin + 1

    for j = 1:size(a,2)
        aj = view(a,:,j)
        @assert allunique(aj)
        if ordered
            @assert issorted(aj)
        end
    end

    if ptol > 0
        p0 = fill(1/n, n)
        if ordered
            @test isapprox(proportions(a, vmin:vmax), p0, atol=ptol)
        else
            b = transpose(a)
            for j = 1:size(b,2)
                bj = view(b,:,j)
                @test isapprox(proportions(bj, vmin:vmax), p0, atol=ptol)
            end
        end
    end
end

import StatsBase: knuths_sample!, fisher_yates_sample!, self_avoid_sample!
import StatsBase: seqsample_a!, seqsample_c!, seqsample_d!

a = zeros(Int, 5, n)
for j = 1:size(a,2)
    knuths_sample!(3:12, view(a,:,j))
end
check_sample_norep(a, (3, 12), 5.0e-3; ordered=false)

test_rng_use(knuths_sample!, 1:10, zeros(Int, 6))

a = zeros(Int, 5, n)
for j = 1:size(a,2)
    fisher_yates_sample!(3:12, view(a,:,j))
end
check_sample_norep(a, (3, 12), 5.0e-3; ordered=false)

test_rng_use(fisher_yates_sample!, 1:10, zeros(Int, 6))

a = zeros(Int, 5, n)
for j = 1:size(a,2)
    self_avoid_sample!(3:12, view(a,:,j))
end
check_sample_norep(a, (3, 12), 5.0e-3; ordered=false)

test_rng_use(self_avoid_sample!, 1:10, zeros(Int, 6))

a = zeros(Int, 5, n)
for j = 1:size(a,2)
    seqsample_a!(3:12, view(a,:,j))
end
check_sample_norep(a, (3, 12), 5.0e-3; ordered=true)

test_rng_use(seqsample_a!, 1:10, zeros(Int, 6))

a = zeros(Int, 5, n)
for j = 1:size(a,2)
    seqsample_c!(3:12, view(a,:,j))
end
check_sample_norep(a, (3, 12), 5.0e-3; ordered=true)

test_rng_use(seqsample_c!, 1:10, zeros(Int, 6))

a = zeros(Int, 5, n)
for j = 1:size(a,2)
    seqsample_d!(3:12, view(a,:,j))
end
check_sample_norep(a, (3, 12), 5.0e-3; ordered=true)

test_rng_use(seqsample_d!, 1:10, zeros(Int, 6))

a = sample(3:12, 5; replace=false)
check_sample_norep(a, (3, 12), 0; ordered=false)

a = sample(3:12, 5; replace=false, ordered=true)
check_sample_norep(a, (3, 12), 0; ordered=true)

# tests of multidimensional sampling

a = sample(3:12, (2, 2); replace=false)
check_sample_norep(a, (3, 12), 0; ordered=false)

@test sample(1:1, (2, 2); replace=true) == ones(Int, 2, 2)

# test of weighted sampling without replacement
a = [1:10;]
wv = Weights([zeros(6); 1:4])
x = vcat([sample(a, wv, 1, replace=false) for j in 1:100000]...)
@test minimum(x) == 7
@test maximum(x) == 10
@test maximum(abs, proportions(x) .- (1:4)/10) < 0.01

x = vcat([sample(a, wv, 2, replace=false) for j in 1:50000]...)
exact2 = [0.117261905, 0.220634921, 0.304166667, 0.357936508]
@test minimum(x) == 7
@test maximum(x) == 10
@test maximum(abs, proportions(x) .- exact2) < 0.01

x = vcat([sample(a, wv, 4, replace=false) for j in 1:10000]...)
@test minimum(x) == 7
@test maximum(x) == 10
@test maximum(abs, proportions(x) .- 0.25) == 0

@test_throws DimensionMismatch sample(a, wv, 5, replace=false)

wv = Weights([zeros(5); 1:4; -1])
@test_throws ErrorException sample(a, wv, 1, replace=false)

#### weighted sampling with dimension
Random.seed!(1);

@test sample([1, 2], Weights([1, 1]), (2,2)) == ones(2,2)
@test sample([1, 2], Weights([0, 1]), (2,2)) == [2 2 ; 2 2]
@test sample(collect(1:4), Weights(1:4), (2,2), replace=false) == [4 1; 3 2]


text:
abstract type AbstractPenalization end
abstract type AbstractPenalizer end
# Penalization API #
# ------------------------------------------------------------
penalty(ph::AbstractProgressiveHedgingSolver) = penalty(ph, ph.penalization)
initialize_penalty!(ph::AbstractProgressiveHedgingSolver) = initialize_penalty!(ph, ph.penalization)
update_penalty!(ph::AbstractProgressiveHedgingSolver) = update_penalty!(ph, ph.penalization)
# ------------------------------------------------------------
include("fixed.jl")
include("adaptive.jl")


text:
# ---
# title: 1523. Count Odd Numbers in an Interval Range
# id: problem1523
# author: zhwang
# date: 2022-03-18
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/count-odd-numbers-in-an-interval-range/description/>
# hidden: true
# ---
#
# Given two non-negative integers `low` and `high`. Return the _count of odd
# numbers between_`low` _and_` high` _  (inclusive)_.
#
#
#
# **Example 1:**
#
#
#
#     Input: low = 3, high = 7
#     Output: 3
#     Explanation: The odd numbers between 3 and 7 are [3,5,7].
#
# **Example 2:**
#
#
#
#     Input: low = 8, high = 10
#     Output: 1
#     Explanation: The odd numbers between 8 and 10 are [9].
#
#
#
# **Constraints:**
#
#   * `0 <= low <= high <= 10^9`
#
#
## @lc code=start
using LeetCode

count_odd_numbers(low::Int, high::Int) = (high - low + 1 + (low & 1)) >> 1

## @lc code=end

## @lc test=start
@testset "1523.count-odd-numbers-in-an-interval-range.jl" begin
    @test count_odd_numbers(3, 7) == count_odd_numbers(3, 8) == 3
    @test count_odd_numbers(2, 7) == 3
    @test count_odd_numbers(0, 10 ^ 9) == 5 * 10 ^ 8
end
## @lc test=end

## @lc info=start
# link: [solution to 1523]()
## @lc info=end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using LinearAlgebra
using ValueShapes, Distributions, ArraysOfArrays
using ForwardDiff, Zygote, DistributionsAD

@testset "test_distribution_transform" begin
    function test_back_and_forth(trg_d, src_d)
        @testset "transform $(typeof(trg_d).name) <-> $(typeof(src_d).name)" begin
            src_v = rand(src_d)
            prev_ladj = 7.9

            @test @inferred(BAT.apply_dist_trafo(trg_d, src_d, src_v, prev_ladj)) isa NamedTuple{(:v,:ladj)}
            trg_v, trg_ladj = BAT.apply_dist_trafo(trg_d, src_d, src_v, prev_ladj)

            @test BAT.apply_dist_trafo(src_d, trg_d, trg_v, trg_ladj) isa NamedTuple{(:v,:ladj)}
            src_v_reco, prev_ladj_reco = BAT.apply_dist_trafo(src_d, trg_d, trg_v, trg_ladj)

            @test src_v ≈ src_v_reco
            @test prev_ladj ≈ prev_ladj_reco
            @test trg_ladj ≈ logabsdet(ForwardDiff.jacobian(x -> unshaped(BAT.apply_dist_trafo(trg_d, src_d, x, prev_ladj).v, varshape(trg_d)), src_v))[1] + prev_ladj
        end
    end

    function get_trgxs(trg_d, src_d, X)
        return (x -> BAT.apply_dist_trafo(trg_d, src_d, x, 0.0).v).(nestedview(X))
    end

    function get_trgxs(trg_d, src_d::Distribution{Univariate}, X)
        return (x -> BAT.apply_dist_trafo(trg_d, src_d, x, 0.0).v).(X)
    end

    function test_dist_trafo_moments(trg_d, src_d)
        @testset "check moments of trafo $(typeof(trg_d).name) <- $(typeof(src_d).name)" begin
            let trg_d = trg_d, src_d = src_d
                X = flatview(rand(src_d, 10^5))
                trgxs = get_trgxs(trg_d, src_d, X)
                unshaped_trgxs = map(unshaped, trgxs)
                @test isapprox(mean(unshaped_trgxs), mean(unshaped(trg_d)), atol = 0.1)
                @test isapprox(cov(unshaped_trgxs), cov(unshaped(trg_d)), rtol = 0.1)
            end
        end
    end

    uniform1 = Uniform(-5.0, -0.01)
    uniform2 = Uniform(0.01, 5.0)

    normal1 = Normal(-10, 1)
    normal2 = Normal(10, 5)

    uvnorm = BAT.StandardUvNormal()

    standnorm1 = BAT.StandardMvNormal(1)
    standnorm2 = BAT.StandardMvNormal(2)

    standuni2 = BAT.StandardMvUniform(2)

    standnorm2_reshaped = ReshapedDist(standnorm2, varshape(standnorm2))

    mvnorm = MvNormal([0.3, -2.9], [1.7 0.5; 0.5 2.3])
    beta = Beta(3,1)
    gamma = Gamma(0.1,0.7)
    dirich = Dirichlet([0.1,4])

    ntdist = NamedTupleDist(
        a = uniform1,
        b = mvnorm,
        c = [4.2, 3.7],
        x = beta,
        y = gamma
    )

    test_back_and_forth(beta, standnorm1)
    test_back_and_forth(gamma, standnorm1)

    test_back_and_forth(mvnorm, mvnorm)
    test_back_and_forth(mvnorm, standnorm2)
    test_back_and_forth(mvnorm, standuni2)

    test_dist_trafo_moments(normal2, normal1)
    test_dist_trafo_moments(uniform2, uniform1)

    test_dist_trafo_moments(beta, gamma)
    test_dist_trafo_moments(uvnorm, standnorm1)

    test_dist_trafo_moments(beta, standnorm1)
    test_dist_trafo_moments(gamma, standnorm1)

    test_dist_trafo_moments(mvnorm, standnorm2)
    test_dist_trafo_moments(dirich, standnorm1)

    test_dist_trafo_moments(mvnorm, standuni2)
    test_dist_trafo_moments(standuni2, mvnorm)

    test_dist_trafo_moments(standnorm2, standuni2)

    test_dist_trafo_moments(mvnorm, standnorm2_reshaped)
    test_dist_trafo_moments(standnorm2_reshaped, mvnorm)
    test_dist_trafo_moments(standnorm2, standnorm2_reshaped)
    test_dist_trafo_moments(standnorm2_reshaped, standnorm2_reshaped)

    test_back_and_forth(ntdist, BAT.StandardMvNormal(5))
    test_back_and_forth(ntdist, BAT.StandardMvUniform(5))

    let
        mvuni = product_distribution([Uniform(), Uniform()])

        x = rand()
        @test_throws ArgumentError BAT.apply_dist_trafo(uvnorm, mvnorm, x, 0.0)
        @test_throws ArgumentError BAT.apply_dist_trafo(uvnorm, standnorm1, x, 0.0)
        @test_throws ArgumentError BAT.apply_dist_trafo(uvnorm, standnorm2, x, 0.0)

        x = rand(2)
        @test_throws ArgumentError BAT.apply_dist_trafo(mvuni, mvnorm, x, 0.0)
        @test_throws ArgumentError BAT.apply_dist_trafo(mvnorm, mvuni, x, 0.0)
        @test_throws ArgumentError BAT.apply_dist_trafo(uvnorm, mvnorm, x, 0.0)
        @test_throws ArgumentError BAT.apply_dist_trafo(uvnorm, standnorm1, x, 0.0)
        @test_throws ArgumentError BAT.apply_dist_trafo(uvnorm, standnorm2, x, 0.0)
    end

    let
        primary_dist = NamedTupleDist(x = Normal(2), c = 5)
        f = x -> NamedTupleDist(y = Normal(x.x, 3), z = MvNormal([1.3 0.5; 0.5 2.2]))
        trg_d = @inferred(HierarchicalDistribution(f, primary_dist))
        src_d = BAT.StandardMvNormal(totalndof(varshape(trg_d)))
        test_back_and_forth(trg_d, src_d)
        test_dist_trafo_moments(trg_d, src_d)
    end

    #=
    using Cuba
    function integrate_over_unit(density::AbstractDensity)
        vs = varshape(density)
        f_cuba(source_x, y) = y[1] = exp(logdensityof(density)(vs(source_x)))
        Cuba.vegas(f_cuba, 1, 1).integral[1]
    end
    =#

    @testset "Custom cdf and quantile for dual numbers" begin
        Dual = ForwardDiff.Dual

        @test BAT._trafo_cdf(Normal(Dual(0, 1, 0, 0), Dual(1, 0, 1, 0)), Dual(0.5, 0, 0, 1)) == cdf(Normal(Dual(0, 1, 0, 0), Dual(1, 0, 1, 0)), Dual(0.5, 0, 0, 1))
        @test BAT._trafo_cdf(Normal(0, 1), Dual(0.5, 1)) == cdf(Normal(0, 1), Dual(0.5, 1))

        @test BAT._trafo_quantile(Normal(0, 1), Dual(0.5, 1)) == quantile(Normal(0, 1), Dual(0.5, 1))
        @test BAT._trafo_quantile(Normal(Dual(0, 1, 0, 0), Dual(1, 0, 1, 0)), Dual(0.5, 0, 0, 1)) == quantile(Normal(Dual(0, 1, 0, 0), Dual(1, 0, 1, 0)), Dual(0.5, 0, 0, 1))
    end

    @testset "trafo broadcasting" begin
        dist = NamedTupleDist(a = Weibull(), b = Exponential())
        smpls = bat_sample(dist, IIDSampling(nsamples = 100)).result
        trafo = BAT.DistributionTransform(Normal, dist)
        @inferred(broadcast(trafo, smpls)) isa DensitySampleVector
        smpls_tr = trafo.(smpls)
        smpls_tr_cmp = [trafo(s) for s in smpls]
        @test smpls_tr == smpls_tr_cmp
	    @test @inferred(varshape(trafo)) == @inferred(varshape(dist)) == trafo._varshape
    end

    @testset "trafo composition" begin
        dist1 = @inferred(NamedTupleDist(a = Normal(), b = Uniform(), c = Cauchy()))
        dist2 = @inferred(NamedTupleDist(a = Exponential(), b = Weibull(), c = Beta()))
        normal1 = Normal()
        normal2 = Normal(2)

        trafo = @inferred(BAT.DistributionTransform(dist1, dist2))
        inv_trafo = @inferred(inv(trafo))

        composed_trafo = @inferred(∘(trafo, inv_trafo))
        @test composed_trafo.source_dist == composed_trafo.target_dist == dist1
        @test composed_trafo ∘ trafo == trafo
        @test_throws ArgumentError  trafo ∘ composed_trafo

        trafo = @inferred(BAT.DistributionTransform(normal1, normal2))
        @test_throws ArgumentError trafo ∘ trafo
    end

    @testset "full density transform" begin
        likelihood = @inferred(NamedTupleDist(a = Normal(), b = Exponential()))
        prior = product_distribution([Normal(), Gamma()])
        posterior_density = PosteriorDensity(likelihood, prior)

        posterior_density_trafod = @inferred(bat_transform(PriorToUniform(), posterior_density, FullDensityTransform()))

        @test posterior_density_trafod.result.orig.likelihood.dist == likelihood
        @test posterior_density_trafod.result.orig.prior.dist == prior

        @test posterior_density_trafod.result.trafo.target_dist isa BAT.StandardMvUniform

        lower_bounds = Float32.([-10, -10, -10])
        upper_bounds = Float32.([10, 10, 10])
        rect_bounds = @inferred(BAT.HyperRectBounds(lower_bounds, upper_bounds))
            mvn = @inferred(product_distribution(Normal.(randn(3))))
        dist_density = @inferred(BAT.DistributionDensity(mvn, rect_bounds))

        dist_density_trafod = @inferred(bat_transform(PriorToUniform(), dist_density, FullDensityTransform()))

        @test dist_density_trafod.trafo.target_dist isa BAT.StandardMvUniform
        @test dist_density_trafod.result.orig == dist_density
        @test dist_density_trafod.trafo.source_dist == dist_density_trafod.result.trafo.source_dist == mvn

        @test dist_density_trafod.trafo._valshape == @inferred(varshape(dist_density))

        dist_density_trafod = @inferred(bat_transform(PriorToGaussian(), dist_density, FullDensityTransform()))

        @test dist_density_trafod.trafo.target_dist isa BAT.StandardMvNormal
        @test dist_density_trafod.result.orig == dist_density
        @test dist_density_trafod.trafo.source_dist == dist_density_trafod.result.trafo.source_dist == mvn

        @test dist_density_trafod.trafo._valshape == @inferred(varshape(dist_density))
    end

    @testset "trafo autodiff pullbacks" begin
        # ToDo: Test for type stability and fix where necessary.

        xs = rand(5)
        @test Zygote.jacobian(BAT._pushfront, xs, 42)[1] ≈ ForwardDiff.jacobian(xs -> BAT._pushfront(xs, 1), xs)
        @test Zygote.jacobian(BAT._pushfront, xs, 42)[2] ≈ vec(ForwardDiff.jacobian(x -> BAT._pushfront(xs, x[1]), [42]))
        @test Zygote.jacobian(BAT._pushback, xs, 42)[1] ≈ ForwardDiff.jacobian(xs -> BAT._pushback(xs, 1), xs)
        @test Zygote.jacobian(BAT._pushback, xs, 42)[2] ≈ vec(ForwardDiff.jacobian(x -> BAT._pushback(xs, x[1]), [42]))
        @test Zygote.jacobian(BAT._rev_cumsum, xs)[1] ≈ ForwardDiff.jacobian(BAT._rev_cumsum, xs)
        @test Zygote.jacobian(BAT._exp_cumsum_log, xs)[1] ≈ ForwardDiff.jacobian(BAT._exp_cumsum_log, xs) ≈ ForwardDiff.jacobian(cumprod, xs)

        src_v = [0.6, 0.7, 0.8, 0.9]
        f = inv(BAT.DistributionTransform(Uniform, DistributionsAD.TuringDirichlet([3.0, 4.0, 5.0, 6.0, 7.0])))
        @test isapprox(ForwardDiff.jacobian(f, src_v), Zygote.jacobian(f, src_v)[1], rtol = 10^-4)
        f = inv(BAT.DistributionTransform(Uniform, Dirichlet([3.0, 4.0, 5.0, 6.0, 7.0])))
        @test isapprox(ForwardDiff.jacobian(f, src_v), Zygote.jacobian(f, src_v)[1], rtol = 10^-4)
        f = inv(BAT.DistributionTransform(Normal, Dirichlet([3.0, 4.0, 5.0, 6.0, 7.0])))
        @test isapprox(ForwardDiff.jacobian(f, src_v), Zygote.jacobian(f, src_v)[1], rtol = 10^-4)
    end
end

@testset "bat_transform_defaults" begin
    mvn = @inferred(product_distribution([Normal(-1), Normal(), Normal(1)]))
    uniform_prior = @inferred(product_distribution([Uniform(-3, 1), Uniform(-2, 2), Uniform(-1, 3)]))

    posterior_uniform_prior = @inferred(PosteriorDensity(mvn, uniform_prior))
    posterior_gaussian_prior = @inferred(PosteriorDensity(mvn, mvn))

    @test @inferred(bat_transform(PriorToGaussian(), posterior_uniform_prior)).result.prior.dist == @inferred(BAT.StandardMvNormal(3))
    @test @inferred(bat_transform(PriorToUniform(), posterior_gaussian_prior)).result.prior.dist == @inferred(BAT.StandardMvUniform(3))
    @test @inferred(bat_transform(NoDensityTransform(), posterior_uniform_prior)).result.prior.dist == uniform_prior
    pd = @inferred(product_distribution([Uniform() for i in 1:3]))
    density = @inferred(BAT.DistributionDensity(pd))
    @test @inferred(bat_transform(NoDensityTransform(), density)).result.dist == density.dist

    # ToDo: Improve comparison for bounds so `.dist` is not required here:
    @inferred(bat_transform(PriorToUniform(), convert(AbstractDensity, BAT.StandardUvUniform()))).result.dist == convert(AbstractDensity, BAT.StandardUvUniform()).dist
    @inferred(bat_transform(PriorToUniform(), convert(AbstractDensity, BAT.StandardMvUniform(4)))).result.dist == convert(AbstractDensity, BAT.StandardMvUniform(4)).dist
    @inferred(bat_transform(PriorToGaussian(), convert(AbstractDensity, BAT.StandardUvNormal()))).result.dist == convert(AbstractDensity, BAT.StandardUvNormal()).dist
    @inferred(bat_transform(PriorToGaussian(), convert(AbstractDensity, BAT.StandardMvNormal(4)))).result.dist == convert(AbstractDensity, BAT.StandardMvNormal(4)).dist
end


text:
################################################# FILE DESCRIPTION #########################################################

# This file contains methods to export data to LightGraphs.jl

################################################# IMPORT/EXPORT ############################################################

export export_adjacency, export_vertex_property, export_edge_property

################################################# HELPERS ##################################################################

process_array(x::Vector) = x

process_array(x::DataArrays.DataArray) = x.data
###
# TODO: Figure out what to do with nonstandard arrays like NullableArrays or DataArrays.
###

"""
This method provides compatibilty with LightGraphs.jl, by returning the graph's adjacency
matrix.

```julia
import LightGraphs
g = Graph(10^5, 10^6)

# Construct a LightGraphs DiGraph from the exported adjacency matrix
LightGraphs.DiGraph(export_adjacency(g))
```
"""
export_adjacency(g::Graph) = copy(indxs(g))

"""
This method provides compatibilty with LightGraphs.jl, by returning an array containing
all values for a vertex property
"""
function export_vertex_property(g::Graph, vprop::Symbol)
   process_array(getvprop(g, :, vprop))
end

"""
This method provides compatibilty with LightGraphs.jl, by returning a SparseMatrixCSC
containing all values for an edge property
"""
function export_edge_property(g::Graph, eprop::Symbol)
   eit = edges(g)
   arr = process_array(geteprop(g, :, eprop))
   sparse(eit.vs, eit.us, arr, nv(g), nv(g))
end


text:
init_bounds(bounds, pomdp, sol) = bounds
init_bounds(t::Tuple, pomdp, sol) = (init_bound(first(t), pomdp, sol), init_bound(last(t), pomdp, sol))
bounds(f::Function, pomdp::POMDP, b::ScenarioBelief) = f(pomdp, b)
bounds(t::Tuple, pomdp::POMDP, b::ScenarioBelief) = (lbound(t[1], pomdp, b), ubound(t[2], pomdp, b))

function bounds_sanity_check(pomdp::POMDP, sb::ScenarioBelief, L_0, U_0)
    if L_0 > U_0
        @warn("L_0 ($L_0) > U_0 ($U_0)   |ϕ| = $(length(sb.scenarios))")
        @info("Try e.g. `IndependentBounds(l, u, consistency_fix_thresh=1e-5)`.", maxlog=1)
    end
    if all(isterminal(pomdp, s) for s in particles(sb))
        if L_0 != 0.0 || U_0 != 0.0
            error(@sprintf("If all states are terminal, lower and upper bounds should be zero (L_0=%-10.2g, U_0=%-10.2g). (try IndependentBounds(l, u, check_terminal=true))", L_0, U_0))
        end
    end
    if isinf(L_0) || isnan(L_0)
        @warn("L_0 = $L_0. Infinite bounds are not supported.")
    end
    if isinf(U_0) || isnan(U_0)
        @warn("U_0 = $U_0. Infinite bounds are not supported.")
    end
end

"""
    IndependentBounds(lower, upper, check_terminal=false, consistency_fix_thresh=0.0)

Specify lower and upper bounds that are independent of each other (the most common case).

# Keyword Arguments
- `check_terminal::Bool=false`: if true, then if all the states in the belief are terminal, the upper and lower bounds will be overridden and set to 0.
- `consistency_fix_thresh::Float64=0.0`: if `upper < lower` and `upper >= lower-consistency_fix_thresh`, then `upper` will be bumped up to `lower`.
"""
struct IndependentBounds{L, U}
    lower::L
    upper::U
    check_terminal::Bool
    consistency_fix_thresh::Float64
end

function IndependentBounds(l, u;
                           check_terminal=false,
                           consistency_fix_thresh=0.0)
    return IndependentBounds(l, u, check_terminal, consistency_fix_thresh)
end

function bounds(bounds::IndependentBounds, pomdp::POMDP, b::ScenarioBelief)
    if bounds.check_terminal && all(isterminal(pomdp, s) for s in particles(b))
        return (0.0, 0.0)
    end
    l = lbound(bounds.lower, pomdp, b)
    u = ubound(bounds.upper, pomdp, b)
    if u < l && u >= l-bounds.consistency_fix_thresh
        u = l
    end
    return (l,u)
end

function init_bounds(bounds::IndependentBounds, pomdp::POMDP, sol::DESPOTSolver)
    return IndependentBounds(init_bound(bounds.lower, pomdp, sol),
                             init_bound(bounds.upper, pomdp, sol),
                             bounds.check_terminal,
                             bounds.consistency_fix_thresh
                            )
end
init_bound(bound, pomdp, sol) = bound

ubound(n::Number, pomdp, b) = convert(Float64, n)
lbound(n::Number, pomdp, b) = convert(Float64, n)

ubound(f::Function, pomdp, b) = f(pomdp, b)
lbound(f::Function, pomdp, b) = f(pomdp, b)

# Value if Fully Observed Under a Policy
struct FullyObservableValueUB{P<:Union{Solver, Policy}}
    p::P
end

ubound(ub::FullyObservableValueUB, pomdp::POMDP, b::ScenarioBelief) = mean(value(ub.p, s) for s in particles(b)) # assumes that all are weighted equally

function init_bound(ub::FullyObservableValueUB{S}, pomdp::POMDP, sol::DESPOTSolver) where S <: Solver
    return FullyObservableValueUB(solve(ub.p, pomdp))
end

# Default Policy Lower Bound

"""
    DefaultPolicyLB(policy; max_depth=nothing, final_value=(m,x)->0.0)
    DefaultPolicyLB(solver; max_depth=nothing, final_value=(m,x)->0.0)

A lower bound calculated by running a default policy on the scenarios in a belief.

# Keyword Arguments
- `max_depth::Union{Nothing,Int}=nothing`: max depth to run the simulation. The depth of the belief will be automatically subtracted so simulations for the bound will be run for `max_depth-b.depth` steps. If `nothing`, the solver's max depth will be used.
- `final_value=(m,x)->0.0`: a function (or callable object) that specifies an additional value to be added at the end of the simulation when `max_depth` is reached. This function will be called with two arguments, a `POMDP`, and a `ScenarioBelief`. It will not be called when the states in the belief are terminal.
"""
struct DefaultPolicyLB{P<:Union{Solver, Policy}, D<:Union{Nothing,Int}, T}
    policy::P
    max_depth::D
    final_value::T
end

function DefaultPolicyLB(policy_or_solver::T;
                         max_depth=nothing,
                         final_value=(m,x)->0.0) where T <: Union{Solver, Policy}
    return DefaultPolicyLB(policy_or_solver, max_depth, final_value)
end

function lbound(lb::DefaultPolicyLB, pomdp::POMDP, b::ScenarioBelief)
    rsum = branching_sim(pomdp, lb.policy, b, lb.max_depth-b.depth, lb.final_value)
    return rsum/length(b.scenarios)
end

function init_bound(lb::DefaultPolicyLB{S}, pomdp::POMDP, sol::DESPOTSolver) where S <: Solver
    policy = solve(lb.policy, pomdp)
    return init_bound(DefaultPolicyLB(policy, lb.max_depth, lb.final_value), pomdp, sol)
end

function init_bound(lb::DefaultPolicyLB{P}, pomdp::POMDP, sol::DESPOTSolver) where P <: Policy
    max_depth = something(lb.max_depth, sol.D)
    return DefaultPolicyLB(lb.policy, max_depth, lb.final_value)
end


text:
export AbstractGrid

abstract type AbstractGrid end

include("fourier_grid.jl")
include("transform.jl")

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: upiterbarg/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
using HDF5
using Statistics

# ΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛ
# Data Defs
# ΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛ

"""
lattice info stored in this object:
number of spins, energy of interaction,
temperature, spin configuration and energy.
"""
mutable struct Lattice{M <: AbstractArray}
    nspins::Int64
    eps::Float64
    T::Float64
    pbc::Bool
    spins::M
    E::Float64
    M::Float64

    # object constructor, throw random spins
    function Lattice(nspins, eps, temp, pbc)
        probs = rand( nspins )
        spins = Vector( [ prob < 0.5 ? -1 : 1 for prob in probs ] )

        new{typeof(spins)}(nspins, eps, temp, pbc, spins, 0.0, 0.0)
    end
end

# ΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛ
# Physics Defs
# ΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛ

"""
Calculate the energy (Hamiltonian) of the lattice.
"""
function energy(lat::Lattice)
    nspins = lat.nspins
    spins  = lat.spins
    enrg   = [ spins[i] * spins[i+1] for i in 1:(nspins-1) ]
    if (lat.pbc)
        append!( enrg, spins[1] * spins[end] )
    end
    lat.E = -lat.eps * sum(enrg)
    return
end

function magnetization(lat::Lattice)
    return sum(lat.spins)
end

"""
The all-knowing Bolztmann factor
"""
boltz( E, T ) = ( k = 1.38064852e-23; exp( -E/(k*T) ) )

function corr(lat::Lattice)
    if (lat.pbc)
        nspins = lat.nspins
        spins  = lat.spins
        savg   = mean(spins)

        c = ones( (nspins-1)÷2 )

        for r in 1:(nspins-1)÷2
            for i in 1:nspins
                ir    = (i + r) == nspins ? nspins : mod( (i + r), nspins )
                il    = (i - r) == 0      ? nspins : mod( (i - r), nspins )
                c[r] += spins[i]*( spins[ir] + spins[il] )
            end
            c[r] *= 1.0 / ( nspins * 2.0 )
            c[r] -= savg * savg
        end
    end
    return c
end

# ΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛ
# Algorithm Defs
# ΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛ

function MC_flip(lat::Lattice, idx)
    nspins = lat.nspins
    spins  = lat.spins

    # These are the calculations to obtain ΔE
    # Ef = Ei + eps * ( spins[idx]*spins[idx-1] + spins[idx]*spins[idx+1] ) -
    #           ( (-1)*spins[idx]*spins[idx-1] + (-1)*spins[idx]*spins[idx+1] )
    # Ef = Ei + 2 * ( spins[idx]*spins[idx-1] + spins[idx]*spins[idx+1] )

    if lat.pbc
        if idx == nspins
            ΔE = 2 * spins[idx] * ( spins[idx-1] + spins[1] )
        elseif idx == 1
            ΔE = 2 * spins[idx] * ( spins[end]   + spins[idx+1] )
        else
            ΔE = 2 * spins[idx] * ( spins[idx-1] + spins[idx+1] )
        end
    else
        if idx == nspins
            ΔE = 2 * spins[idx] * spins[idx-1]
        elseif idx == 1
            ΔE = 2 * spins[idx] * spins[idx+1]
        else
            ΔE = 2 * spins[idx] * ( spins[idx-1] + spins[idx+1] )
        end
    end
    ΔM = -2 * spins[idx]
    return ΔE, ΔM
end

function MC_step!(lat::Lattice)
    nspins = lat.nspins
    idx    = ceil(Int64, rand() * nspins )

    ΔE, ΔM = MC_flip(lat, idx)

    if ΔE <= 0.0
        lat.spins[idx] *= -1
        lat.E += ΔE
        lat.M += ΔM
    else
        if ( rand() < boltz(ΔE, lat.T) )
            lat.spins[idx] *= -1
            lat.E += ΔE
            lat.M += ΔM
        end
    end
    return
end

function MC_run!(iter, lat::Lattice, data)

    energy(lat)
    magnetization(lat)

    data["mag"][1]     = lat.M
    data["energy"][1]  = lat.E
    data["spins"][:,1] = lat.spins

    Eavg = zeros(iter[2])

    i = 1
    for step in 1:iter[1]
        MC_step!(lat)
        Eavg[i] = lat.E
        i += 1

        if step%iter[2] == 0
            k = step ÷ iter[2]
            println("Step $((k-1)*iter[2])")
            data["energy"][k]  = mean(Eavg)
            data["mag"][k]     = lat.M
            data["spins"][:,k] = lat.spins
            i = 1
        end
    end

    data["corr"][:] = corr(lat)
    return
end

# ΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛ
# Interface
# ΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛΛ

function main(args)
    nspins = parse(Float64, args[1]) |> Int
    eps    = parse(Float64, args[2])
    temp   = parse(Float64, args[3])
    pbc    = true
    iters  = parse(Float64, args[4]) |> Int
    wrt_it = parse(Float64, args[5]) |> Int

    lattice_pbc = Lattice(nspins, eps, temp, pbc)

    data = h5open("./Ising1D.hdf5", "w")
    data_dims = iters ÷ wrt_it
    create_dataset( data, "energy", datatype(1.0), dataspace( (data_dims,) ) )
    create_dataset( data, "mag",    datatype(1.0), dataspace( (data_dims,) ) )
    create_dataset( data, "corr",   datatype(1.0), dataspace( ((nspins-1) ÷ 2,) ) )
    create_dataset( data, "spins",  datatype(1.0), dataspace( (nspins, data_dims) ) )

    create_dataset( data, "params", datatype(1.0), dataspace( (5,) ) )

    data["params"][1] = nspins
    data["params"][2] = eps
    data["params"][3] = temp
    data["params"][4] = iters
    data["params"][5] = wrt_it

    MC_run!( (iters, wrt_it), lattice_pbc, data )

    io = open("julia.out", "w")
    println(io, "Initial Energy = $(data["energy"][1])")
    println(io, "Final Energy = $(data["energy"][end])")
    close(io)

    "Initial Energy = $(data["energy"][1])" |> display
    "Final Energy = $(data["energy"][end])" |> display

    close(data)
end

main(ARGS)

max_stars_repo_name: fysixs/IsingModel

max_issues_repo_name: fysixs/IsingModel

max_forks_repo_name: fysixs/IsingModel


text:
#=

  From BLOG examples/mixture-of-gaussian-infinite.blog
  """
  Gaussian mixture model with infinite number of components

  @author leili
  @date 2014-07-10
  """

  Note: The BLOG model remove numComponent == 0 via a condition.
        Here we remove it using Truncated(Poisson(2),0,Inf)

  Summary Statistics
    parameters      mean       std   naive_se       mcse         ess       rhat   ess_per_sec
        Symbol   Float64   Float64    Float64   Float64?    Float64?   Float64?      Float64?

  numComponent    2.2097    1.2555     0.0126     0.0179   4752.1622     0.9999     1866.5209
          z[1]    1.6089    0.9319     0.0093     0.0121   4800.4385     0.9999     1885.4825
          z[2]    1.6034    0.9092     0.0091     0.0128   5013.1564     0.9999     1969.0324
          z[3]    1.6042    0.9161     0.0092     0.0134   5163.1008     0.9999     2027.9265
          z[4]    1.5902    0.8956     0.0090     0.0136   4864.7405     1.0001     1910.7386
         mu[1]    0.2789    0.4861     0.0049     0.0064   5152.2290     1.0000     2023.6563
         mu[2]    0.2106    0.5170     0.0064     0.0090     missing    missing       missing
         mu[3]    0.1611    0.5371     0.0093     0.0134     missing    missing       missing
         mu[4]    0.1342    0.5243     0.0138     0.0170     missing    missing       missing
     x_post[1]    0.3065    1.1066     0.0111     0.0121   8080.3308     1.0000     3173.7356
     x_post[2]    0.3533    1.1055     0.0111     0.0116   8837.3254     1.0001     3471.0626
     x_post[3]    0.3147    1.0961     0.0110     0.0105   9327.8597     0.9999     3663.7312
     x_post[4]    0.3229    1.0999     0.0110     0.0120   8585.7202     0.9999     3372.2389
         mu[5]    0.0978    0.5476     0.0237     0.0348     missing    missing       missing
         mu[6]    0.0579    0.5407     0.0401    missing     missing    missing       missing
         mu[7]    0.0565    0.4821     0.0568    missing     missing    missing       missing
         mu[8]    0.3198    0.4958     0.0891    missing     missing    missing       missing
         mu[9]    0.1723    0.4504     0.2014    missing     missing    missing       missing

  Distributions of variable numComponent (num:0)
  1.00000 =>    3517  (0.351700)
  2.00000 =>    3139  (0.313900)
  3.00000 =>    1900  (0.190000)
  4.00000 =>     908  (0.090800)
  5.00000 =>     354  (0.035400)
  6.00000 =>     110  (0.011000)
  7.00000 =>      41  (0.004100)
  8.00000 =>      26  (0.002600)
  9.00000 =>       5  (0.000500)

  Cf ~/webppl/mixture_of_gaussian_infinite.wppl

=#

using Turing, StatsPlots
include("jl_utils.jl")

@model function mixture_of_gaussian_infinite(x=[0.2,1.0,0.5,0.6])
    n = length(x)

    # Note: We have to set 0..Inf (not 1..Inf), otherwise
    # the value 1 will not be drawn.
    numComponent ~ truncated(Poisson(2),0,Inf)
    z ~ filldist(DiscreteUniform(1,numComponent),n)
    mu ~ filldist(Uniform(-1,1), numComponent)

    for i in 1:n
        x[i] ~ Normal(mu[z[i]],1.0)
    end

    x_post = tzeros(n) # Vector{Real}(undef, n)
    for i in 1:n
        x_post[i] ~ Normal(mu[z[i]],1.0)
    end

end

x=[0.2,1.0,0.5,0.6]
model = mixture_of_gaussian_infinite(x)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(15), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, NUTS(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns, :numComponent)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
#!/usr/bin/julia

# Change working directory to test dir
cd(joinpath(@__DIR__, "..", "test"))

tests = map(realpath, readlines(`find . -name "*.jl"`))
julia = joinpath(Sys.BINDIR, Base.julia_exename())

import Base.Threads: @threads

@threads for test in tests
    runlines = String[]
    content = String[]
    open(test, "r") do io
        last_empty = false
        while !eof(io)
            line = readline(io)
            # Skip multiple empty lines
            if isempty(line)
                if last_empty
                    continue
                else
                    last_empty = true
                end
            else
                last_empty = false
            end

            if startswith(line, r"#\s*+RUN:\s*")
                push!(runlines, line)
            elseif startswith(line, r"#\s*CHECK")
                # ignore
            else
                push!(content, line)
            end
        end
    end
    @assert length(runlines) == 1
    runline = first(runlines)
    m = match(r"#\s*RUN:\s*(.*?)\s?\|.*", runline)
    @assert m !== nothing
    cmd = m.captures[1]
    cmd = replace(cmd, "%s"    => test)
    cmd = replace(cmd, "julia" => julia)
    cmd = replace(cmd, "2>&1" => "")

    in = IOBuffer()
    process = run(pipeline(Cmd(Base.shell_split(cmd)), stdout=in, stderr=in))
    seekstart(in)

    if !success(process)
        @error "Updating test failed" test
        write(stderr, in)
        continue
    end

    fname, io = mktemp()

    write(io, runline, "\n")
    for line in content
        write(io, line, "\n")
    end

    first_check = true
    while !eof(in)
        line = readline(in)
        isempty(line) && continue

        if startswith(line, "module {")
            first_check = true
        elseif startswith(line, "after" ) || contains(line, r"loc\(")
            first_check = true
            continue
        end

        if first_check
            println(io)
            write(io, "# CHECK: ")
            first_check = false
        else
            write(io, "# CHECK-NEXT: ")
        end

        # These are runtime constants
        line = replace(line, r"llvm\.mlir\.constant\(\d+ : i64\)" => "llvm.mlir.constant({{[0-9]+}} : i64)")

        write(io, line, "\n")
    end
    close(io)

    mv(fname, test, force=true)
end

max_stars_repo_name: femtomc/brutus

max_issues_repo_name: femtomc/brutus

max_forks_repo_name: femtomc/brutus


text:
using NCDatasets, DataStructures

fname = "filename_fv.nc"

if isfile(fname)
    rm(fname)
end

ds = NCDataset(fname,"c")

# Dimensions

sz = (1000,500,100)

ds.dim["longitude"] = sz[1]
ds.dim["latitude"] = sz[2]
ds.dim["time"] = sz[3]

# Declare variables

ncv1 = defVar(ds,"v1", UInt8, ("longitude", "latitude", "time"), fillvalue = UInt8(255), attrib = [
    "add_offset"                => -1.0,
    "scale_factor"              => 0.5,
])

# Define variables

for n = 1:sz[3]
    @show n
    ncv1[:,:,n] = rand(1:100,sz[1],sz[2])
    ncv1[:,1,n] = missing
end

close(ds)


text:
#module Spike

include("PSTH.jl") # export PSTH, PSTHPlot

#end # Spike


text:
# GPU implementation of generator update of ADMM.
function generator_kernel(
    baseMVA::Float64, ngen::Int, gen_start::Int,
    u::CuDeviceArray{Float64,1}, v::CuDeviceArray{Float64,1},
    l::CuDeviceArray{Float64,1}, rho::CuDeviceArray{Float64,1},
    pgmin::CuDeviceArray{Float64,1}, pgmax::CuDeviceArray{Float64,1},
    qgmin::CuDeviceArray{Float64,1}, qgmax::CuDeviceArray{Float64,1},
    c2::CuDeviceArray{Float64,1}, c1::CuDeviceArray{Float64,1}
)
    I = threadIdx().x + (blockDim().x * (blockIdx().x - 1))
    if (I <= ngen)
        pg_idx = gen_start + 2*(I-1)
        qg_idx = gen_start + 2*(I-1) + 1
        u[pg_idx] = max(pgmin[I],
                        min(pgmax[I],
                            (-(c1[I]*baseMVA + l[pg_idx] -
                            rho[pg_idx]*v[pg_idx])) / (2*c2[I]*(baseMVA^2) + rho[pg_idx])))
        u[qg_idx] = max(qgmin[I],
                        min(qgmax[I],
                            (-(l[qg_idx] - rho[qg_idx]*v[qg_idx])) / rho[qg_idx]))
    end

    return
end

function generator_kernel(
    baseMVA::Float64, ngen::Int, gen_start::Int,
    u, v, l, rho,
    pgmin::Array{Float64,1}, pgmax::Array{Float64,1},
    qgmin::Array{Float64,1}, qgmax::Array{Float64,1},
    c2::Array{Float64,1}, c1::Array{Float64,1}
)
    for I=1:ngen
        pg_idx = gen_start + 2*(I-1)
        qg_idx = gen_start + 2*(I-1) + 1
        u[pg_idx] = max(pgmin[I],
                        min(pgmax[I],
                            (-(c1[I]*baseMVA + l[pg_idx] -
                               rho[pg_idx]*v[pg_idx])) / (2*c2[I]*(baseMVA^2) + rho[pg_idx])))
        u[qg_idx] = max(qgmin[I],
                        min(qgmax[I],
                            (-(l[qg_idx] - rho[qg_idx]*v[qg_idx])) / rho[qg_idx]))
    end

    return
end

function generator_kernel(
    gen_mod::GeneratorModel{CuArray{Float64,1}},
    baseMVA::Float64, u, v, l, rho
)
    nblk = div(gen_mod.ngen, 32, RoundUp)
    tgpu = CUDA.@timed @cuda threads=32 blocks=nblk generator_kernel(baseMVA, gen_mod.ngen, gen_mod.gen_start,
                u, v, l, rho, gen_mod.pgmin, gen_mod.pgmax, gen_mod.qgmin, gen_mod.qgmax, gen_mod.c2, gen_mod.c1)
    return tgpu
end

function generator_kernel(
    gen_mod::GeneratorModel{Array{Float64,1}},
    baseMVA::Float64, u, v, l, rho
)
    tcpu = @timed generator_kernel(baseMVA, gen_mod.ngen, gen_mod.gen_start,
                u, v, l, rho, gen_mod.pgmin, gen_mod.pgmax, gen_mod.qgmin, gen_mod.qgmax, gen_mod.c2, gen_mod.c1)
    return tcpu
end

function generator_kernel_two_level(
    baseMVA::Float64, ngen::Int, gen_start::Int,
    u::CuDeviceArray{Float64,1}, x::CuDeviceArray{Float64,1}, z::CuDeviceArray{Float64,1},
    l::CuDeviceArray{Float64,1}, rho::CuDeviceArray{Float64,1},
    pgmin::CuDeviceArray{Float64,1}, pgmax::CuDeviceArray{Float64,1},
    qgmin::CuDeviceArray{Float64,1}, qgmax::CuDeviceArray{Float64,1},
    c2::CuDeviceArray{Float64,1}, c1::CuDeviceArray{Float64,1}
)
    I = threadIdx().x + (blockDim().x * (blockIdx().x - 1))
    if (I <= ngen)
        pg_idx = gen_start + 2*(I-1)
        qg_idx = gen_start + 2*(I-1) + 1
        u[pg_idx] = max(pgmin[I],
                        min(pgmax[I],
                            (-(c1[I]*baseMVA + l[pg_idx] + rho[pg_idx]*(-x[pg_idx] + z[pg_idx]))) / (2*c2[I]*(baseMVA^2) + rho[pg_idx])))
        u[qg_idx] = max(qgmin[I],
                        min(qgmax[I],
                            (-(l[qg_idx] + rho[qg_idx]*(-x[qg_idx] + z[qg_idx]))) / rho[qg_idx]))
    end

    return
end

function generator_kernel_two_level(
    baseMVA::Float64, ngen::Int, gen_start::Int,
    u, x, z, l, rho,
    pgmin::Array{Float64,1}, pgmax::Array{Float64,1},
    qgmin::Array{Float64,1}, qgmax::Array{Float64,1},
    c2::Array{Float64,1}, c1::Array{Float64,1}
)
    for I=1:ngen
        pg_idx = gen_start + 2*(I-1)
        qg_idx = gen_start + 2*(I-1) + 1
        u[pg_idx] = max(pgmin[I],
                        min(pgmax[I],
                            (-(c1[I]*baseMVA + l[pg_idx] + rho[pg_idx]*(-x[pg_idx] + z[pg_idx]))) / (2*c2[I]*(baseMVA^2) + rho[pg_idx])))
        u[qg_idx] = max(qgmin[I],
                        min(qgmax[I],
                            (-(l[qg_idx] + rho[qg_idx]*(-x[qg_idx] + z[qg_idx]))) / rho[qg_idx]))
    end

    return
end

function generator_kernel_two_level(
    gen_mod::GeneratorModel{CuArray{Float64,1}},
    baseMVA::Float64, u::CuArray{Float64,1}, xbar::CuArray{Float64,1},
    zu::CuArray{Float64,1}, lu::CuArray{Float64,1}, rho_u::CuArray{Float64,1}
)
    nblk = div(gen_mod.ngen, 32, RoundUp)
    tgpu = CUDA.@timed @cuda threads=32 blocks=nblk generator_kernel_two_level(baseMVA, gen_mod.ngen, gen_mod.gen_start,
                u, xbar, zu, lu, rho_u, gen_mod.pgmin, gen_mod.pgmax, gen_mod.qgmin, gen_mod.qgmax, gen_mod.c2, gen_mod.c1)
    return tgpu
end

function generator_kernel_two_level(
    gen_mod::GeneratorModel{Array{Float64,1}},
    baseMVA::Float64, u, xbar, zu, lu, rho_u
)
    tcpu = @timed generator_kernel_two_level(baseMVA, gen_mod.ngen, gen_mod.gen_start,
                u, xbar, zu, lu, rho_u, gen_mod.pgmin, gen_mod.pgmax, gen_mod.qgmin, gen_mod.qgmax, gen_mod.c2, gen_mod.c1)
    return tcpu
end


text:
# # Tutorials
# A suite of concrete examples are provided here as a guidance for constructing experiments.
# ## Balance Law
# An introduction on components within a balance law is provided.
# ## Atmos
# Showcase drivers for atmospheric modelling in GCM, single stack, and LES simulations are provided.
# - Dry Idealzed GCM:  The Held-Suarez configuration is used as a guidance to create a driver that runs a simple GCM simulation.
# - Single Element Stack:  The Burgers Equations with a passive tracer is used as a guidance to run the simulation on a single element stack.
# - LES Experiment:  The dry rising bubble case is used as a quigance in creating an LES driver.
# - Topography:  Experiments of dry flow over prescirbe topography (Agnesi mountain) are provided for:
#     * Linear Hydrostatic Mountain
#     * Linear Non-Hydrostatic Mountain
# ## Ocean
# A showcase for Ocean model is still under construction.
# ## Land
# Examples are provided in constructing balance law and solving for fundemental equations in land modelling.
# - Heat:  A tutorial shows how to create a HeatModel to solve the heat equation and visualize the outputs.
# - Soil:  Examples of solving fundemental equations in the soil model are provided.
#     * Hydraulic Functions:  a tutorial to specify the hydraulic function in the Richard's equation.
#     * Soil Heat Equations:  a tutorial for solving the heat equation in the soil.
#     * Coupled Water and Heat:  a tutorial for solving interactive heat and wateri in the soil model.
# ## Numerics (need to be moved to How-to-Guide)
# - System Solvers:  Two numerical methods to solve the linear system Ax=b are provided.
#     * Conjugate Gradient
#     * Batched Generalized Minimal Residual
# - DG Methods
#     * Filters
# ## Diagnostics
# A diagnostic tool that can
# - generate statistics for MPIStateArrays
# - validate with reference values
# for debugging purposes.


text:
##################
# time series data
##################

# The time series object and interface are intended to model possibly
# high-dimensional arrays where the first dimesion represents and remaining
# dimension represent a scalar, vector, etc. that evolves in time.

# As such, this package seeks to treat time-series data as representing arrays
# with ambient and often common time axes. Alternatively, time series could be
# thought of as modeling functions that map time to the space of the array
# cross-sections. In this sense, this package seeks to lift operations on arrays
# to operate on time series.

# There are other packages such as JuliaStats/TimeSeries or JuliaData/DataFrames
# that could be applicable to time-series data. However, none of these seem to
# be quite applicable due to inflexibility in representation of time and
# inflexibility with data dimension respectively

struct TimeSeries{D, N, AT <: AbstractVector, AD <: AbstractArray} <: AbstractArray{D,N}

  time::AT
  data::AD

  # Here we enforce the constraint that AT <: AbstractArray{D,N}
  function TimeSeries(time::AT, data::AD) where {AT <: AbstractVector, AD <: AbstractArray}
    if size(time, 1) != size(data, 1)
      error("first dimension of time and data do not match")
    end

    assert_sorted(time)

    new{eltype(AD), ndims(AD), AT, AD}(time, data)
  end
end

# Map style constructor
TimeSeries(f, time::AbstractArray) = TimeSeries(time, map(f, time))

TimeSeries(data) = TimeSeries(Base.OneTo(size(data,1)), data)

# Copy constructor
TimeSeries(x::TimeSeries) = TimeSeries(get_time(x), get_data(x))

#######################################
# Internal type queries for convenience
#######################################
time_type(::Type{TimeSeries{<:Any,<:Any,AT}}) where {AT} = AT
time_type(::T) where {T <: TimeSeries} = time_type(T)
data_type(::Type{TimeSeries{D,N,AT,AD}}) where {D,N,AT,AD} = AD
data_type(::T) where {T <: TimeSeries} = data_type(T)

time_eltype(::Type{TimeSeries{<:Any,<:Any,AT}}) where {AT} = eltype(AT)
time_eltype(::T) where {T <: TimeSeries} = time_eltype(T)

######################
# TimeSeries interface
######################

# Note: "_data" methods can be deprecated given the array interface although
# some may still be usefull or necessary for odd jobs.

# A good example of when get_data would be useful is when it makes sense to
# discard the time array for the underlying data array. This can happen when the
# "data" encodes a range or an abstract notion of time such as an iteration
# number. While "iterations" may have been published in real time and been used
# to construct a time array, they may also be used as the time component of new
# TimeSeries objects. In such cases, it generally makes sense to use an
# appropriate Array for the time data rather than defining the TimeSeries in
# terms of another TimeSeries.

get_time(x::TimeSeries) = x.time
get_data(x::TimeSeries) = x.data
get_data(x::AbstractArray) = x

mutate_time(f, x::TimeSeries) = TimeSeries(f(x.time), x.data)
mutate_data(f, x::TimeSeries) = TimeSeries(x.time, f(x.data))
map_time(f, x::TimeSeries) = mutate_time(x->map(f, x), x)
map_data(f, x::TimeSeries) = mutate_data(x->map(f, x), x)

# Do progressively more broad equality checks for comparisons of indices. It is
# possible that the first two can be combined.
indices_match(a, b) = a === b || a == b || all(a .== b)
indices_match(a::TimeSeries, b::TimeSeries) = indices_match(get_time(a),
                                                            get_time(b))
#################
# Array interface
#################

# See documentation at:
# https://docs.julialang.org/en/latest/manual/interfaces/#man-interface-array-1

import Base.size, Base.IndexStyle, Base.getindex, Base.setindex!,
Base.similar

# This deviates somewhat from the documentation but also fits the actual
# definitions in the code so these seem to be the right way to define everything
size(x::TimeSeries) = size(get_data(x))
IndexStyle(::Type{T}) where {T<:TimeSeries} = IndexStyle(data_type(T))
getindex(x::TimeSeries, i::Int) = getindex(get_data(x), i)
getindex(x::TimeSeries{<:Any,N}, I::Vararg{Int,N}) where {N} = getindex(get_data(x),I...)
setindex!(x::TimeSeries, v, I...) = setindex!(get_data(x), v, I...)

# "similar" will produce a new TimeSeries rather than an Array so that we can
# keep the time data around. The data is in turn defined recursively according
# to its type with another call to similar.
function similar(x::TimeSeries, ::Type{S}, dims::Dims{N}) where {S, N}
  if size(x, 1) != dims[1]
    error("similar not defined for varying first dimension of a TimeSeries")
  end

  TimeSeries(get_time(x), similar(get_data(x), S, dims))
end

##########
# Printing
##########

import Base.string, Base.print, Base.show
string(x::TimeSeries) = string("Time:\n", get_time(x), "\nData:\n", get_data(x))
print(io::IO, x::TimeSeries) = print(io, string(x))
show(io::IO, x::TimeSeries) = print(io, string(x))

#######
# Other
#######

function assert_sorted(x::AbstractArray, prefix = "")
  if !issorted(x)
    error(prefix * "sample times should be sorted")
  end
end

max_stars_repo_name: mcorah/RosDataProcess

max_issues_repo_name: mcorah/RosDataProcess

max_forks_repo_name: mcorah/RosDataProcess


text:
include("regression.jl")
include("classification.jl")
include("multiclass.jl")
include("event.jl")
# include("generic_likelihood.jl")
function (l::AbstractLikelihood)(::Real, ::Real)
    return error("pdf not implemented for likelihood $(typeof(l))")
end

Distributions.loglikelihood(l::AbstractLikelihood, y::Real, f) = log(l(y, f))

## Default function for getting gradient ##
function ∇loglikehood(l::AbstractLikelihood, y::Real, f::Real)
    return first(ForwardDiff.gradient(x -> loglikelihood(l, y, x[1]), [f]))
end

function ∇loglikehood(l::AbstractLikelihood, y::Real, f::AbstractVector)
    return ForwardDiff.gradient(x -> loglikelihood(l, y, x), f)
end

function hessloglikehood(l::AbstractLikelihood, y::Real, f::Real)
    return first(ForwardDiff.hessian(x -> loglikelihood(l, y, x[1]), [f]))
end

function hessloglikelihood(l::AbstractLikelihood, y::Real, f::AbstractVector)
    return ForwardDiff.hessian(x -> loglikelihood(l, y, x), f)
end

implemented(::AbstractLikelihood, ::AbstractInference) = false

isaugmentable(::AbstractLikelihood) = false

Base.length(::AbstractLikelihood) = 1

num_latent(::AbstractLikelihood) = 1

Base.iterate(l::AbstractLikelihood) = (l, nothing)
Base.iterate(::AbstractLikelihood, ::Any) = nothing


text:
"""
MPResults

This module has the code you need for the mixed precision example.
"""
module MPResults

export
# Functions
heqtest,
PlotHist,
Fname,
data_harvest,
data_populate,
readmpdata,
writempdata,
MakeTable,
PlotData

using PyPlot
using LinearAlgebra
using Printf
using SIAMFANLEquations
using SIAMFANLEquations.TestProblems
using SIAMFANLEquations.Examples

include("heqtest.jl")
include("data_populate.jl")
include("Data_IO.jl")
include("PlotData.jl")
include("TableData.jl")

end #module

max_stars_repo_name: ctkelley/MPResults


text:
module TopOpt

const PENALTY_BEFORE_INTERPOLATION = true
using Requires, Reexport

macro cuda_only(mod, code)
    return esc(
        quote
            @init @require CUDASupport = "97986420-7ec3-11e9-24cd-4f0e301eb539" @eval $mod begin
                $code
            end
        end,
    )
end

abstract type AbstractDevice end
struct CPU <: AbstractDevice end
struct GPU <: AbstractDevice end
whichdevice(::Any) = CPU()

# GPU utilities
module GPUUtils end

@reexport using Nonconvex, NonconvexMMA, NonconvexSemidefinite, NonconvexPercival

# Utilities
include(joinpath("Utilities", "Utilities.jl"))
using .Utilities

# Topopology optimization problem definitions
include(joinpath("TopOptProblems", "TopOptProblems.jl"))

using LinearAlgebra, Statistics
using Reexport, Parameters, Setfield
@reexport using .TopOptProblems

# Truss Topopology optimization problem definitions
include(joinpath("TrussTopOptProblems", "TrussTopOptProblems.jl"))
@reexport using .TrussTopOptProblems

using Ferrite, StaticArrays

using ForwardDiff, IterativeSolvers#, Preconditioners
@reexport using VTKDataTypes

const DEBUG = Base.RefValue(false)

# FEA solvers
include(joinpath("FEA", "FEA.jl"))
using .FEA

# Chequeurboard filter
include(joinpath("CheqFilters", "CheqFilters.jl"))
using .CheqFilters

# Objective and constraint functions
@reexport using Flux
include(joinpath("Functions", "Functions.jl"))
@reexport using .Functions

# Various topology optimization algorithms
include(joinpath("Algorithms", "Algorithms.jl"))
using .Algorithms

macro init_cuda()
    return esc(
        quote
            const CuArrays = CUDASupport.CuArrays
            const CUDAdrv = CUDASupport.CUDAdrv
            const CUDAnative = CUDASupport.CUDAnative
            const GPUArrays = CUDASupport.GPUArrays
            CuArrays.allowscalar(false)
            const dev = CUDAdrv.device()
            const ctx = CUDAdrv.CuContext(dev)
            using .CuArrays, .CUDAnative
            using .GPUArrays: GPUVector, GPUArray
        end,
    )
end

@cuda_only GPUUtils include("GPUUtils/GPUUtils.jl")
@cuda_only Utilities include("Utilities/gpu_utilities.jl")
@cuda_only TopOptProblems include("TopOptProblems/gpu_support.jl")
@cuda_only FEA include("FEA/gpu_solvers.jl")
@cuda_only CheqFilters include("CheqFilters/gpu_cheqfilter.jl")
@cuda_only Functions include("Functions/gpu_support.jl")
@cuda_only Algorithms include("Algorithms/SIMP/gpu_simp.jl")

export TopOpt,
    simulate,
    TopOptTrace,
    DirectDisplacementSolver,
    PCGDisplacementSolver,
    StaticMatrixFreeDisplacementSolver,
    SensFilter,
    DensityFilter,
    Displacement,
    Compliance,
    CG,
    Direct,
    Assembly,
    MatrixFree,
    FEASolver,
    Optimizer,
    SIMP,
    ContinuationSIMP,
    BESO,
    GESO,
    PowerContinuation,
    ExponentialContinuation,
    LogarithmicContinuation,
    CubicSplineContinuation,
    SigmoidContinuation,
    Continuation,
    save_mesh,
    CPU,
    GPU,
    DefaultCriteria,
    EnergyCriteria,
    PowerPenalty,
    RationalPenalty,
    SinhPenalty,
    MMA87,
    MMA02,
    HeavisideProjection,
    ProjectedPenalty,
    PowerPenalty
end


text:
# TODO: Need to finish
#       Need to test
function fit_fs_imcmc_pt!(cfs::ConstantsFS,
                          dfs::DataFS;
                          nmcmc::Int, nburn::Int,
                          # Args for PT:
                          tempers::Vector{Float64},
                          inits=nothing,
                          save_all_states=false,
                          randpair=0.0,
                          # End of PT Args.
                          # Args are for iMCMC:
                          batchprop::Float64=0.1,
                          batchsizes=nothing,
                          prior_thin::Int=2,
                          imcmc_burn_prop=0.6,
                          swap_freq::Float64=1.0,
                          # End of iMCMC args.
                          tfs::Union{Nothing, Vector{TunersFS}}=nothing,
                          monitors=[monitor1, monitor2],
                          fix::Vector{Symbol}=Symbol[],
                          thins::Vector{Int}=[2, nsamps_to_thin(10, nmcmc)],
                          ndden_samps::Int=200,
                          printFreq::Int=0,
                          checkpoint=0,
                          computeDIC::Bool=true, computeLPML::Bool=true,
                          computedden::Bool=true,
                          sb_ibp::Bool=false,
                          use_repulsive::Bool=true,
                          Z_marg_lamgam::Float64=1.0,
                          Z_marg_lamgam_decay_rate::Float64=100.0,
                          Z_marg_lamgam_min::Float64=0.05,
                          verbose::Int=1,
                          time_updates=false,
                          seed::Int=-1)

  printMsg(iter::Int, msg::String) = if printFreq > 0 && iter % printFreq == 0
    print(msg)
  end

  # How frequently to thin dden
  thin_dden = nsamps_to_thin(ndden_samps, nmcmc)

  # Number of temperatures
  num_tempers = length(tempers)

  @assert mod(num_tempers, 2) == 0
  if tfs != nothing
    @assert num_tempers == length(tfs)
  end

  # Assert swap frequency is in unit interval
  @assert 0 < swap_freq <= 1.0

  # Cache for swap counts
  swapcounts = zeros(Int, num_tempers, num_tempers)

  # Cache for pair counts
  paircounts = zeros(Int, num_tempers, num_tempers)

  # Set random seed if needed
  if seed >= 0
    Random.seed!(seed)
  end

  if verbose >= 1
    fixed_vars_str = join(fix, ", ")
    if fixed_vars_str == ""
      fixed_vars_str = "nothing"
    end
    println("fixing: $fixed_vars_str")
    println("Use stick-breaking IBP: $(sb_ibp)")
    println("Z_marg_lamgam: $(Z_marg_lamgam)")
    println("Z_marg_lamgam_decay_rate: $(Z_marg_lamgam_decay_rate)")
    println("Z_marg_lamgam_min: $(Z_marg_lamgam_min)")
    println("use_repulsive: $(use_repulsive)")
    println("batchprop: $(batchprop)")
    println("batchsizes: $(batchsizes)")
    flush(stdout)
  end

  @assert printFreq >= -1
  if printFreq == 0
    numPrints = 10
    printFreq = Int(ceil((nburn + nmcmc) / numPrints))
  end

  # Instantiate (but not initialize) CPO stream
  if computeLPML
    cpoStream = MCMC.CPOstream{Float64}()
  end

  # DIC
  if computeDIC
    dicStream = initDicStream(dfs.data)
    loglikeDIC(param::DICparam) = computeLoglikeDIC(dfs.data, param)

    convertStateToDicParam(s::State)::DICparam = let
      _convertStateToDicParam(s, cfs.constants, dfs.data)
    end
  end

  # Cache for data density
  dden = Matrix{Vector{Float64}}[]

  # Update function
  function update!(states, args, iter::Int)
    # Whether or not to marginalize over lambda and gamma.
    # We want to do this more often at the beginning, and less at the end.
    zmarg = ((Z_marg_lamgam - Z_marg_lamgam_min) *
             exp(-iter/Z_marg_lamgam_decay_rate) +
             Z_marg_lamgam_min) > rand()

    s_arg_vec = pmap(states, args) do s, arg
      tu = (arg[:c].constants.temper == 1) && time_updates

      # For reproducibility.
      if seed > - 1
        temper_idx = findfirst(tau -> tau == arg[:c].constants.temper, tempers)
        Random.seed!(iter + temper_idx + seed)
      end

      # Update state using trained prior
      imcmc_all_params = let
        iter < imcmc_burn_prop * nburn
      end

      update_via_trained_prior!(s, dfs, arg[:c], arg[:t],
                                batchprop, prior_thin,
                                batchsizes=batchsizes,
                                fix=fix, use_repulsive=use_repulsive,
                                Z_marg_lamgam=zmarg, sb_ibp=sb_ibp,
                                time_updates=time_updates, temper=1.0,
                                minibatch_update_all_params=imcmc_all_params,
                                verbose=verbose-2)
      (s, arg)
    end

    states = [sa[1] for sa in s_arg_vec]
    args = [sa[2] for sa in s_arg_vec]

    # Swap Chains
    if verbose > 0
      println()
    end

    llf(s, tuner) = compute_marg_loglike(s, cfs.constants, dfs.data, tuner)

    if swap_freq > rand()
      swapchains!(states, llf, tempers,
                  paircounts=paircounts, swapcounts=swapcounts,
                  randpair=randpair, verbose=verbose)
    end
    if iter == nburn
      println("swapcounts / paircounts:")
      println(swapcounts ./ (paircounts .+ 1e-6))
      println("Resetting swapcounts, paircounts ...")
      swapcounts .= 0.0
      paircounts .= 0.0
    end

    # Pull out inner componenets for convenience
    s = states[1]
    c = args[1][:c]
    d = dfs
    ll = args[1][:ll]

    # # Append loglike
    append!(ll, compute_marg_loglike(s.theta, c.constants, d.data, 1.0))

    if computedden && iter > nburn && (iter - nburn) % thin_dden == 0
      # NOTE: `datadensity(s, c, d)` returns an (I x J) matrix of vectors of
      # length g.
      append!(dden, [datadensity(s.theta, c.constants, d.data)])
    end

    if computeLPML && iter > nburn
      # Inverse likelihood for each data point
      like = [[compute_like(i, n, s.theta, c.constants, d.data)
               for n in 1:d.data.N[i]] for i in 1:d.data.I]

      # Update (or initialize) CPO
      MCMC.updateCPO(cpoStream, vcat(like...))

      # Add to printMsg
      printMsg(iter, " -- LPML: $(MCMC.computeLPML(cpoStream))")
    end

    if computeDIC && iter > nburn
      # Update DIC
      MCMC.updateDIC(dicStream, s.theta, updateParams,
                     loglikeDIC, convertStateToDicParam)

      # Add to printMsg
      printMsg(iter, " -- DIC: $(MCMC.computeDIC(dicStream, loglikeDIC,
                                                 paramMeanCompute))")

      DICg = MCMC.DIC_gelman(ll[(nburn+1):end])
      printMsg(iter, " -- DIC_gelman: $(DICg)")
    end

    printMsg(iter, "\n")
    flush(stdout)

    return states, args
  end

  # Create vectors of states
  states = if inits == nothing
    println("Using random inits!")
    [let
       s = genInitialState(cfs.constants, dfs.data)
       s.eps .= 0.0
       StateFS{Float64}(s, dfs)
     end for _ in tempers]
  else
    @assert length(inits) == num_tempers
    for _init in inits
      _init.theta.eps .= 0.0
    end
    inits
  end

  # Create Args
  args = [let
            ll = Float64[]  # FIXME
            c = deepcopy(cfs)
            c.constants.temper = tempers[i]
            t = if tfs == nothing
              TunersFS(Tuners(dfs.data.y, cfs.constants.K),
                       states[1].theta, dfs.X)
            else
              tfs[i]
            end
            Dict(:ll => ll, :c => c, :t => t)
          end
          for i in 1:num_tempers]

  println("Running Gibbs sampler ...")
  samples, states, args = let
    MCMC.gibbs_pt(states, args, update!, monitors=monitors,
                  thins=thins, nmcmc=nmcmc, nburn=nburn,
                  printFreq=printFreq,
                  save_all_states=save_all_states,
                  printlnAfterMsg=false)
  end

  out = Dict(:samples => samples,
             :lastState => states[1],
             :inits => inits,
             :c => cfs,
             :d => dfs,
             :lls => [arg[:ll] for arg in args],
             :save_all_states => save_all_states,
             :tempers => tempers,
             :swap_freq => swap_freq,
             :randpair => randpair,
             :paircounts => paircounts,
             :swapcounts => swapcounts)

  if computeDIC || computeLPML
    LPML = computeLPML ? MCMC.computeLPML(cpoStream) : NaN
    Dmean, pD = computeDIC ? MCMC.computeDIC(dicStream, loglikeDIC,
                                             paramMeanCompute,
                                             return_Dmean_pD=true) : (NaN, NaN)

    ll1 = args[1][:ll]
    DICg = MCMC.DIC_gelman(ll1[(nburn+1):end])

    metrics = Dict(:LPML => LPML,
                   :DIC => Dmean + pD,
                   :Dmean => Dmean,
                   :pD => pD,
                   :DICg => DICg)
    println()
    println("metrics:")
    for (k, v) in metrics
      out[k] = v
      println("$k => $v")
    end
    flush(stdout)
  end

  if computedden
    out[:dden] = dden
  end

  out[:nburn] = nburn
  out[:nmcmc] = nmcmc
  out[:batchprop] = batchprop
  out[:batchsizes] = batchsizes
  out[:prior_thin] = prior_thin

  # Return all states if requested
  if save_all_states
    out[:all_last_states] = deepcopy(states)
  end

  return out
end


text:
using ControlBenchmarks
using ControlSystems

benchmarkProb = controlbenchmark( DrumBoiler() )

@test ControlSystems.nstates( benchmarkProb.sys ) == 9
@test ControlSystems.ninputs( benchmarkProb.sys ) == 3


text:
"""
    solve(problem::Karger, gradient)

Solve the finite pulse Karger model (FPK) using precomputed effective diffusion tensors
`difftensors`.
"""
function solve(problem::Karger, gradient::ScalarGradient)
    (; model, difftensors, odesolver, timestep) = problem
    (; mesh, T₂, ρ, κ, γ) = model

    dt = timestep
    f = gradient.profile
    d = gradient.dir
    g = gradient.amplitude
    TE = echotime(f)

    # Deduce sizes
    ncompartment, nboundary = size(mesh.facets)

    # Volumes
    volumes = get_cmpt_volumes(mesh)

    # Compute surface areas
    surface_areas = spzeros(ncompartment, nboundary)
    for icmpt = 1:ncompartment
        # Finite elements
        points = mesh.points[icmpt]
        facets = mesh.facets[icmpt, :]

        # Surface area
        for iboundary = 1:nboundary
            if !isempty(facets[iboundary])
                surface_areas[icmpt, iboundary], =
                    get_mesh_surface(points, facets[iboundary])
            end
        end
    end

    # Associate surface areas to the compartment interfaces they measure
    tau_inv = spzeros(ncompartment, ncompartment)
    for iboundary = 1:nboundary
        inds = surface_areas[:, iboundary].nzind
        if length(inds) == 2
            i₁ = inds[1]
            i₂ = inds[2]
            tmp = κ[iboundary] * surface_areas[i₁, iboundary]
            tau_inv[i₁, i₂] = tmp
            tau_inv[i₂, i₁] = tmp
        end
    end
    tau_inv = tau_inv ./ volumes'
    A = tau_inv - spdiagm((tau_inv * volumes) ./ volumes)

    # Relaxation tensor
    R = spdiagm(1 ./ T₂)

    # Initial signal
    S₀ = volumes .* ρ

    # Update function for linear ODE operator
    function update_func(J, u, p, t)
        (; A, ADC_diag, R, f, γ, g) = p
        J .= A .- (integral(f, t)^2 * (γ * g)^2) .* ADC_diag .- R
    end

    ADC_diag = spdiagm([d' * D * d for D ∈ difftensors])

    p = (; A, ADC_diag, R, f, γ, g)
    J_prototype = A - ADC_diag - R
    J = DiffEqArrayOperator(J_prototype; update_func)

    prob = ODEProblem(J, S₀, (0.0, TE), p)
    sol = OrdinaryDiffEq.solve(prob, odesolver, dt = dt)

    sol.u[end]
end


text:
mcmcresultsname(::Val{:ols}) = "ans_$(PARAM[:crspdatalabel])" *
  "_$(PARAM[:mcmcresultssuffix])" *
  "_$(PARAM[:mcmcmodel])" *
  "_$(PARAM[:mcmcprior])" *
  "_$(Dates.format(now(),"yyyymmdd_HHMMSS"))"

#optresultsname() = optresultsname(Val(PARAM[:itermodel]))
#primary function for estimating the measure using optimizations
#note that this still uses the same volume parts construct as flux
function formmcmcregression(panel::DataFrame, ms::MeasureSpec,
  MCMCModel::Val = Val(PARAM[:mcmcmodel]))

  #issorted(panel, [:permno, :date]) || error("panel must be sorted")
  sort!(panel, :date)

  #in contrast to other methods, drop any column where we don't have
  #the lag weight (improves performance I think)

  #conditionforturing!(panel::SubDataFrame, ms::MeasureSpec)
  zs::ZSpec = condition!(panel, ms)

  #testing
  if PARAM[:testmeasure]
    testmcmc(panel, zs)
  end

  if PARAM[:refreshmeasure]
    results::DataFrame = mcmcmodel(panel, zs, MCMCModel)
    #destandardize!(results, zs, intercept=false)

    #this will be the name of the results file
    #resultname::String = mcmcresultsname(λ⁺)

    #results |> CSV.write("$(PARAM[:outputpath])\\$(resultname).csv")

  end

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA