text:
# SimpleMPSJulia
# Density matrix renormalization group (DMRG) in matrix product state (MPS)

# This file contains the definition of matrix product state
# For theoretical backgrounds, see the [reference]:
# Ulrich Schollwöck, The density-matrix renormalization group in the age of matrix product states,
# Annals of Physics, 326 (2011), 96-192

module SimpleMPS

export MatrixProductOperator, MatrixProductState, search_groundstate, expectation

using TensorOperations
using LinearAlgebra
using LinearMaps
using IterativeSolvers
using Statistics

"""
Base class for MPS and MPO
"""
abstract type MatrixProduct end

Base.length(mp::MatrixProduct) = length(mp.mp)
Base.iterate(mp::MatrixProduct, state=1) = iterate(mp.mp, state)
Base.show(io::IO, mp::MatrixProduct) = print(io, [size(m) for m in mp.mp])

struct MatrixProductOperator <: MatrixProduct
    mp::Array{Array}
    MatrixProductOperator(mp) = new(mp)
    function MatrixProductOperator(single_mo::Array, site_num::Int)
        # the first MPO site, only contains the last row
        mpo_1 = single_mo[end:end, :, :, :]
        # the last MPO site, only contains the first column
        mpo_n = single_mo[:, :, :, 1:1]
        new([[mpo_1]; fill(single_mo, site_num - 2); [mpo_n]])
    end
end

struct MatrixProductState <: MatrixProduct
    mp::Array{Array}
    MatrixProductState(mp) = new(mp)
    """
    Construct MPS with random matrices.
    The physical dimension is obtained from `mpo` and bond dimension is set to `m`
    """
    function MatrixProductState(mpo::MatrixProductOperator, m::Int)
        bond_dim = [[1]; fill(m, length(mpo) - 1); [1]]
        pdim = [size(m)[2] for m in mpo]
        # Firstly randomly fill in the matrices
        mp = [rand(bond_dim[i], pdim[i], bond_dim[i+1]) for i in 1:length(mpo)]
        mps = new(mp)
        # Then perform canonicalisation (and normalization)
        for i in 1:length(mpo)-1
            push_center!(mps, i)
        end
        return mps
    end
end

"""
Flip MPO and MPS. After the operation the left side of MPO (MPS) becomes the
right side and vice versa.
Note that <mps|mpo|mps> conserves after the operation.
"""
Base.reverse(mpo::MatrixProductOperator) = MatrixProductOperator(reverse([permutedims(m, (4, 2, 3, 1)) for m in mpo]))
Base.reverse(mps::MatrixProductState) = MatrixProductState(reverse([permutedims(m, (3, 2, 1)) for m in mps]))

"""
Perform left-canonicalisation on the `idx`th site of `mps` so that the canonical center is moved to `idx+1`.
The matrices are normalized by the way.
"""
function push_center!(mps::MatrixProductState, idx::Int)
    orig_shape = size(mps.mp[idx])
    ms = reshape(mps.mp[idx], (prod(orig_shape[1:2]), :))
    f = qr(ms)
    mps.mp[idx] = reshape(Array(f.Q), orig_shape[1], orig_shape[2], :)
    @tensor mps.mp[idx+1][a, c, d] := f.R[a, b] * mps.mp[idx+1][b, c, d]
    mps.mp[idx+1] /= norm(mps.mp[idx+1])
end

"""
Calculate environment on the next site based on previous environment `res`
`ms` and `mo` are local sites of the MPS and MPO.
"""
function update_env(res, ms, mo)
    """
    Graphical notation: (* for MPS and # for MPO)
    *--a--*--e
    |     |
    |     d
    |     |
    #--b--#--g
    |     |
    |     f
    |     |
    *--c--*--h
    """
    @tensor res[e, g, h] := res[a, b, c] * conj(ms)[a, d, e] * mo[b, d, f, g] * ms[c, f, h]
end

"""
Calculate an array of left-side enrivonment tensor based on `mps` and `mpo`.
The index `i` of the array represents the left environment of the `i`th site of the MPS.
"""
function get_env(mps::MatrixProductState, mpo::MatrixProductOperator)::Array{Array}
    sentinel = ones(1, 1, 1)
    env = [sentinel]
    for (ms, mo) in zip(mps, mpo)
        res = update_env(env[length(env)], ms, mo)
        push!(env, res)
    end
    push!(env, sentinel)
    return env
end

"""
Calculate the expectation value <mps|mpo|mps> assuming `mps` is normalized
"""
function expectation(mps::MatrixProductState, mpo::MatrixProductOperator)::AbstractFloat
    return get_env(mps, mpo)[end-1][1][1][1]
end

"""
Sweep `mps` from left to right to find the ground state of the Hamiltonian `mpo`.
The `mps` should be right-canonical
The energies calculatd during the sweep are appended to `energies`
and `env` is passed in as the environment from the right side
"""
function optimize_oneround!(mps::MatrixProductState, mpo::MatrixProductOperator, energies::Array, env::Array)
    sentinel = ones(1, 1, 1)
    # left-side environment and right-side environment. Corresponds to $L$ and $R$ in the reference literature.
    l_env = [sentinel]
    r_env = env

    for cur_idx in 1:length(mps)
        l, mo, r = l_env[cur_idx], mpo.mp[cur_idx], r_env[end-cur_idx-1]

        orig_shape = size(mps.mp[cur_idx])
        # don't optimize on the first site because it has probably been optmized during the previous sweep
        if cur_idx != 1
            # find the lowest eigenvalue of the local Hamitonian operator.
            # The computation can be accelerated by using a "sparse" representation of the operator
            # because some eigensolvers (such as `lobpcg` in this case) only need to know what would happen
            # if the operator(matrix) is applied on a certain state(vector).
            # See [power method](https://en.wikipedia.org/wiki/Power_iteration) as an example
            function h_map_f(ms)
                ms = reshape(ms, orig_shape)
                """
                Graphical notation: (* for MPS and # for MPO)
                *--a-- --h--*
                |     |     |
                |     f     |
                |     |     |
                #--b--#--g--#
                |     |     |
                |     d     |
                |     |     |
                *--c--*--e--*
                """
                @tensor res[a, f, h] := l[a, b, c] * ms[c, d, e] * mo[b, f, d, g] * r[h, g, e]
                return res
            end
            # some boilerplate to call the eigensolver
            new_dim = prod(orig_shape)
            h_map = LinearMap(h_map_f, new_dim, new_dim, issymmetric=true)
            # find the eigenvalue and eigenstate
            lobpcg_res = lobpcg(h_map, false, reshape(mps.mp[cur_idx], new_dim), 1)
            # update variables based on the result
            new_energy = lobpcg_res.λ[1]
            push!(energies, new_energy)
            mps.mp[cur_idx] = reshape(lobpcg_res.X, orig_shape)
        end

        # move the canonical center to the right side
        if cur_idx != length(mps)
            push_center!(mps, cur_idx)
        end
        # update the left-side environment
        push!(l_env, update_env(l, mps.mp[cur_idx], mpo.mp[cur_idx]))
    end

    push!(l_env, sentinel)

    return l_env
end

"""
Find the ground state of the Hamiltonian `mpo`
"""
function search_groundstate(mpo::MatrixProductOperator, bond_dimension::Int=20)
    # random initial guess
    mps = MatrixProductState(mpo, bond_dimension)
    l_env = get_env(mps, mpo)
    energies = [expectation(mps, mpo)]
    while true
        # do a right to left and then left to right sweep
        for i in 1:2
            # the mps is left canonical, so flip over to make it right canonical
            # so that we can achieve sweeping from right to left by sweeping from left to right
            mps, mpo = reverse.([mps, mpo])
            l_env = optimize_oneround!(mps, mpo, energies, l_env)
        end
        # consider it as converged if the lastest energies do not vary much
        if 10 < length(energies) && isapprox(std(energies[end-10:end]), 0, atol=1e-8)
            break
        end
    end
    return mps, energies
end

end


text:
export lookup

function lookup(w::Var, x::Var)
    ydata = lookup(w.data, x.data)
    Var(ydata, ∇lookup!, (w,x))
end
lookup(w::Node, x) = Node(lookup, (w,x))
lookup(w, x::Node) = Node(lookup, (w,x))

function lookup(w::Matrix{T}, x::Array{Int}) where T
    s = Base.setindex(size(x), size(x,1)*size(w,1), 1)
    y = fill!(similar(w,s...), 0)
    n = size(w, 1)
    for i = 1:length(x)
        x[i] == 0 && continue
        yi = (i-1) * n + 1
        wi = (x[i]-1) * n + 1
        copyto!(y, yi, w, wi, n)
    end
    y
end

@generated function lookup(w::CuMatrix{T}, x::CuArray{Cint}) where T
    Ct = cstring(T)
    k = Kernel("""
    __global__ void lookup($Ct *y, int sizeY, $Ct *w, int *x, int n) {
        int idxY = blockIdx.x * blockDim.x + threadIdx.x;
        if (idxY >= sizeY) return;

        int j = idxY / n;
        int i = idxY - n * j;
        y[idxY] = x[j] == 0 ? 0 : w[(x[j]-1) * n + i];
    }""")
    quote
        n = size(w, 1)
        ysize = Base.setindex(size(x), n*size(x,1), 1)
        y = similar(w, ysize)
        gdims, bdims = cudims(length(y))
        $k(gdims, bdims, pointer(y), length(y), pointer(w), pointer(x), n)
        y
    end
end

function ∇lookup!(y::Var, w::Var, x::Var)
    isnothing(w.grad) && return
    ∇lookup!(y.grad, w.grad, x.data)
end

function ∇lookup!(gy::Array, gw::Array, x::Array{Int})
    n = size(gw, 1)
    for i = 1:length(x)
        x[i] == 0 && continue
        yi = (i-1) * n + 1
        wi = (x[i]-1) * n + 1
        addto!(gw, wi, gy, yi, n)
    end
end

@generated function ∇lookup!(gy::CuArray{T}, gw::CuArray{T}, x::CuArray{Cint}) where T
    Ct = cstring(T)
    k = Kernel("""
    __global__ void lookup_grad($Ct *gy, int sizeY, $Ct *gw, int *x, int n) {
        int idxY = blockIdx.x * blockDim.x + threadIdx.x;
        if (idxY >= sizeY) return;

        int j = idxY / n;
        int i = idxY - n * j;
        if (x[j] > 0) {
            int idxW = (x[j]-1) * n + i;
            atomicAdd(&gw[idxW], gy[idxY]);
        }
    }""")
    quote
        n = size(gw, 1)
        gdims, bdims = cudims(length(gy))
        $k(gdims, bdims, pointer(gy), length(gy), pointer(gw), pointer(x), n)
    end
end

ext: jl


text:
"""
A `Panel` represents one of the panels in a Panel matrix.

It is used like a view into the matrix.
"""
struct Panel{T,D,P1,P2,C,F1,F2,L1,L2} <: AbstractMatrix{T}
    data::D # the actual data (typically a view)
    panel_size::Tuple{P1, P2} # (rows, columns) in memory
    panel_class::Val{C} # :CM = column major, :CN = row major
    pad_first::Tuple{F1, F2} # # (rows, columns) of padding before the first element of first tile
    pad_last::Tuple{L1, L2} # (rows, columns) of padding after the last element until tile ends

    function Panel{T}(data::D, panel_size::Tuple{P1, P2}, panel_class::Val{C},
            pad_first::Tuple{F1, F2}, pad_last::Tuple{L1, L2}) where {T, D<:Union{Ptr{T}, AbstractArray{T}},
            P1<:Integer, P2<:Integer, C, F1<:Integer, F2<:Integer, L1<:Integer, L2<:Integer}
        new{T,D,P1,P2,C,F1,F2,L1,L2}(data, panel_size, panel_class, pad_first, pad_last)
    end
end

Base.size(x::Panel) = x.panel_size .- x.pad_first .- x.pad_last
Base.eltype(::Type{Panel{T}}) where {T} = T
Base.eltype(::Panel{T}) where {T} = T

"""
    linear_index(panel_class, panel_size, i, j)

Compute a linear index into a panel without padding.
(This function is mainly for internal use. Currently, it does not check bounds,
but this may change.)
"""
@inline linear_index(::Val{:CM}, sz, i, j, pf=static.((0,0))) = i + pf[1] + (j + pf[2] - 1) * sz[1]
@inline linear_index(::Val{:RM}, sz, i, j, pf=static.((0,0))) = linear_index(Val(:CM), reverse(sz), j, i, reverse(pf))

@inline function Base.setindex!(x::Panel, v, i::Integer, j::Integer)
    @boundscheck checkbounds(x, i, j)
    @inbounds x.data[linear_index(x.panel_class, x.panel_size, i, j)] = v
end

@inline function Base.getindex(x::Panel, i::Integer, j::Integer)
    @boundscheck checkbounds(x, i, j)
    @inbounds x.data[linear_index(x.panel_class, x.panel_size, i, j)]
end

const panel_padding_error = ErrorException("Incorrect padding size.")
const panel_index_error = ErrorException("Panel index out of bounds.")

"""
   pie(i, b, p)

Pad-if-equal: returns p if i==b, otherwise 0.
"""
pie(i, b, p) = p === static(0) ? static(0) : i == b ? Int(p) : 0

"""
    panel_view(x, i, j, pad_first, pad_last)

Get a `Panel` view of panel (i,j) in matrix `x`.

If `pad_first` and `pad_last` are not specified, then they are computed based
on the indices (i, j) and the padding of A.
"""
@inline function panel_view(x::PanelMatrix{T}, i::Integer, j::Integer,
        pad_first=(nothing, nothing), pad_last=(nothing, nothing)) where T
    pad_first=something.(pad_first, pie.((i, j), (1, 1), x.pad_first))
    pad_last=something.(pad_last, pie.((i, j), x.n_panels, x.pad_last))
    length = prod(x.panel_size)
    if all(isa.(StaticInteger, typeof.(x.panel_size)))
        length = static(length)
    end
    zeroth = ((i - 1) + (j - 1) * x.panel_stride) * length
    @boundscheck begin
        all(0 .<= pad_first) && all(0 .<= pad_last) && all(pad_first .+ pad_last .<= x.panel_size) || throw(panel_padding_error)
        all(1 .<= (i, j) .<= x.n_panels) || throw(panel_index_error)
    end
    Panel{T}(@inbounds(view(x.data, LengthUnitRange(zeroth, length))), x.panel_size, x.panel_class, pad_first, pad_last)
end

"""
    full_panel_view(x, i, j)

Get a full `Panel` view of panel (i,j) in matrix `x`. A view of an entire panel is
returned even if the matrix does not cover this entire panel.
"""
Base.@propagate_inbounds function full_panel_view(x::PanelMatrix{T}, i::Integer, j::Integer) where T
    panel_view(x, i, j, static.((0, 0)), static.((0, 0)))
end

"""
    get_panel(x, i, j)

Returns a copy of panel(i, j), whith the size depending on the location of the
panel.

    get_panel(x, i, j, pad_first, pad_last)

Returns a copy of panel (i,j) in matrix `x`. The amount of padding is specified
by input arguments, and may be different from the true size of this panel.
If the panel size and all of the padding is statically known, then an `SMatrix`
is returned. Otherwise a normal `Matrix` is returned.
"""
@inline function get_panel(x::PanelMatrix{T}, i::Integer, j::Integer,
        pad_first=(nothing, nothing), pad_last=(nothing, nothing)) where {T}
    v = panel_view(x, i, j, pad_first, pad_last)
    Matrix{T}(v)
end

@inline function get_panel(x::PanelMatrix{T,D,S1,S2,StaticInteger{P1},StaticInteger{P2}},
        i::Integer, j::Integer,
        pad_first::Tuple{StaticInteger{F1}, StaticInteger{F2}},
        pad_last::Tuple{StaticInteger{L1}, StaticInteger{L2}}) where {T,D,S1,S2,P1,P2,F1,F2,L1,L2}
    @boundscheck begin
        all(0 .<= pad_first) && all(0 .<= pad_last) && all(pad_first .+ pad_last .<= x.panel_size) || throw(panel_padding_error)
        all(1 .<= (i, j) .<= x.n_panels) || throw(panel_index_error)
    end
    v = @inbounds panel_view(x, i, j, pad_first, pad_last)
    SMatrix{P1-F1-L1, P2-F2-L2, T, (P1-F1-L1)*(P2-F2-L2)}(v)
end

"""
    get_full_panel(x, i, j)

Get a copy of panel (i, j) in matrix `x`. An entire panel is
returned even if the matrix does not cover this entire panel.
"""
Base.@propagate_inbounds get_full_panel(x, i, j) = get_panel(x, i, j, (static(0), static(0)), (static(0), static(0)))

"""
    set_panel!(x, y, i, j, pad_first, pad_last)

Set panel (i,j) in matrix `x` to `y`. The amount of padding is specified
by input arguments, and may be different from the true size of this panel.
"""
@inline function set_panel!(x::PanelMatrix, y::AbstractMatrix, i::Integer, j::Integer,
        pad_first=(nothing, nothing), pad_last=(nothing, nothing))
    v = panel_view(x, i, j, pad_first, pad_last)
    @boundscheck begin
        all(0 .<= pad_first) && all(0 .<= pad_last) && all(pad_first .+ pad_last .<= x.panel_size) || throw(panel_padding_error)
        checkbounds(y, 1:x.panel_size[1]-pad_first[1]-pad_last[1], 1:x.panel_size[2]-pad_first[2]-pad_last[2])
        all(1 .<= (i, j) .<= x.n_panels) || throw(panel_index_error)
    end
    @inbounds v .= y
end

"""
    set_full_panel!(x, y, i, j)

Set the full panel (i,j) in matrix `x` to `y`. An entire panel is
set even if the matrix does not cover this entire panel.
"""
Base.@propagate_inbounds set_full_panel!(x, y, i, j) = set_panel!(x, y, i, j, (static(0), static(0)), (static(0), static(0)))


text:
"""
    Params(AbstractType: TwiliteTimeline)

ApiParams is a abstruct type for any parameter subtype. So Any subtype of ApiParams should sataify all the requirements of the ApiParams.
"""
abstract type Params end

"""
    ParamsGetTweets(ConcreteType: TwiliteTimeline)

ParamGetTweets is a defined datatype to store the Resource Parameter mentioned in the twitter API document which we have to send when qurying the API in the request body.

To be able to use twitter API the fild value shoud satisfy the following criteria:
- When Providing valid `count` value the `include_rts` arguments must have to be 1; by default providing any valid `count` value of the `include_rts` will be set to 1 silintly
- When providing constructor argument value for the boolean field such as `trim_user`, `exclude_replies` and `include_rts` the value mush have to be 1 or 0. By the default the values are set to the default value according to Twitter API doc.

# Example
```
julia> resource_params = ParamsGetTweet("realDonaldTrump", count=100)
ParamsGetTweet("realDonaldTrump", nothing, nothing, 100, nothing, false, true)

julia> resource_params = ParamsGetTweet("realDonaldTrump", count=100, exclude_replies=true)
ParamsGetTweet("realDonaldTrump", nothing, nothing, 100, nothing, false, true)
```
"""
struct ParamsGetTweets <: Params
    screen_name::String
    user_id::Union{Integer, Nothing}
    since_id::Union{Integer, Nothing}
    count::Union{Integer, Nothing}
    max_id::Union{Integer, Nothing}
    trim_user::Union{Bool, Nothing}
    exclude_replies::Union{Bool, Nothing}

    # Custom Constructor for the Type Signature
    ParamsGetTweets(
        screen_name;
        user_id = nothing,
        science_id = nothing,
        count=1,
        max_id=nothing,
        trim_user=false,
        exclude_replies=true,
    ) = (
    count isa Integer && count < 1
    ? error("Count must have to be greter than or equal to 1")
    : new(screen_name, user_id, science_id, count, max_id, trim_user, exclude_replies)
    )
end

"""
    ParamsPostTweet(ConcreteType: TwiliteTimeline)

ParamsPostTweet is a defined datatype to store the Resource Parameter mentioned in the twitter API document which we have to send when qurying the API in the request body.

The status field is a mandatory field and rest of the field is set according to twitter api doc with default values. There is a lot of condition implied when using optional fields so it is suggested that use should first see the twitter version 1 api doc statuses/update.json.

# Example

```
julia> test = ParamsPostTweet("asdf asdf asdf a test tweet.")
ParamsPostTweet("asdf asdf asdf a test tweet.", nothing, false, nothing, nothing, nothing, false, nothing, nothing, nothing, nothing, false, false, true, nothing)

```
"""
struct ParamsPostTweet <: Params
    status::String
    in_reply_to_status_id::Union{Integer, Nothing}
    auto_populate_reply_metadata::Union{Bool, Nothing}
    # Comma Separeted List Not sure what?
    exclude_reply_user_ids::Union{String, Nothing}
    attachment_url::Union{String, Nothing}
    # Comma Separeted List Not sure what?
    media_ids::Union{String, Nothing}
    possibly_sensitive::Union{Bool, Nothing}
    lat::Union{Float32, Nothing}
    long::Union{Float32, Nothing}
    place_id::Union{String, Nothing}
    display_coordinates::Union{Bool, Nothing}
    trim_user::Union{Bool, Nothing}
    enable_dmcommands::Union{Bool, Nothing}
    fail_dmcommands::Union{Bool, Nothing}
    card_uri::Union{String, Nothing}

    ParamsPostTweet(
        status;
        in_reply_to_status_id = nothing,
        auto_populate_reply_metadata = false,
        exclude_reply_user_ids = nothing,
        attachment_url = nothing,
        media_ids = nothing,
        possibly_sensitive = false,
        lat = nothing,
        long = nothing,
        place_id = nothing,
        display_coordinates = nothing,
        trim_user = false,
        enable_dmcommands = false,
        fail_dmcommands = true,
        card_uri = nothing
        ) = new(
            status,
            in_reply_to_status_id,
            auto_populate_reply_metadata,
            exclude_reply_user_ids,
            attachment_url,
            media_ids,
            possibly_sensitive,
            lat,
            long,
            place_id,
            display_coordinates,
            trim_user,
            enable_dmcommands,
            fail_dmcommands,
            card_uri
        )
end

"""
    params_to_dict(p::T) where {T <: Params}

Is a internal function to convert the user provided Twitter API parameter into dictionary which `oauth_request_resource` from `OAuth` package expect. The method silently add `include_rts` paramater value to 1 when providing an valid `count paramter`.
"""
function params_to_dict(p::T) where {T <: Params}
    params = Dict{String, Union{String, Integer, Float32}}()
    for fn in fieldnames(T)
        fv = getfield(p, fn)
        if fv |> !isnothing
            fn_str = string(fn)
            if fv isa Bool
                params[fn_str] = string(fv)
            else
                params[fn_str] = fv
            end
        end
    end
    if haskey(params, "count")
        params["include_rts"] = 1
    end
    return params
end


text:
## Julia code to load-in chemostat data

using CSV, DataFrames, Dierckx, PyPlot

function deal_missing(x,t)
    # remove missing data (in csv file this is coded as NaN)
    n = length(x)
    x1 = zeros(n)
    t1 = zeros(n)

    j=0
    for i = 1:n
      if !isnan(x[i])
         j += 1
         x1[j] = x[i]
         t1[j] = t[i]
      end
    end
    x1[1:j], t1[1:j]
end

function get_data(exp_code, make_regular)
    # load-in data from csv-file
    # exp_code : string coding the experimental trial ("C1" to "C10")
    # make_regular : flag
    #                  true  --> interpolate to regular time spacing (needed for wavelets)
    #                  false --> keep the raw data

    forced = false
    if exp_code == "C8" || exp_code == "C9"  # these are the externally forced systems
        forced = true
    end

    # load-in the data
    df = CSV.read(exp_code * ".csv")  # read in data
    t  = df[:,1]  # time (days)
    m  = df[:,2]  # algae (10^6 cells / ml)
    b  = df[:,3]  # rotifers (animals / ml)
    er = df[:,4]  # egg ratio
    e  = df[:,5]  # eggs      (eggs / ml)
    d  = df[:,6]  # dead animals (dead animals / ml)
    s = df[:,7]   # nitrogen concentration in external medium (mumol N/l)

    if make_regular == false   # either return the raw data  ....
        return forced,t,m,b,er,e,d,s
    end

    # ... else, do some minimal pre-processing

    # eliminate missing values, this yields new time vectors for every state
    M,tM = deal_missing(m,t)
    B,tB = deal_missing(b,t)
    E,tE  = deal_missing(e,t)
    D,tD  = deal_missing(d,t)
    ER,tER  = deal_missing(er,t)
    S,tS  = deal_missing(s,t)

    # interpolation to a regular one-day sampling interval
    t_end = trunc(t[end])
    tt = 0:t_end
    mm = evaluate(Spline1D(tM,M),tt)
    bb = evaluate(Spline1D(tB,B),tt)
    ee = evaluate(Spline1D(tE,E),tt)
    er = evaluate(Spline1D(tER,ER),tt)
    dd = evaluate(Spline1D(tD,D),tt)

    # the splines artificially create some smoothing error
    # so we need to correct this
    # this is a simple hack, but works..
    if forced
        ss = evaluate(Spline1D(tS,S),tt)
        for i=1:length(ss)  # spline artifacts are largest for the external forcing signal
            if ss[i] < 100
                ss[i] = 0.0
            else
                ss[i] = 160.0
            end
        end
    else
        ss = S[1] * ones(length(tt))  # these are constants, so no splines needed
    end
    # if interpolation artificially created small negative values, we remove them
    mm = minimum(mm) < 0 ? mm .- minimum(mm) : mm
    bb = minimum(bb) < 0 ? bb .- minimum(bb) : bb
    ee = minimum(ee) < 0 ? ee .- minimum(ee) : ee
    er = minimum(er) < 0 ? er .- minimum(er) : er
    dd = minimum(dd) < 0 ? dd .- minimum(dd) : dd
    if forced
        ss = minimum(ss) < 0 ? ss .- minimum(ss) : ss
    end

    return forced,tt,mm,bb,er,ee,dd,ss
end

function test()
    # simple test code to load data and make a plot
    exp_code = "C10"
    make_regular = true
    forced, t, m, b, er, e, d, s = get_data(exp_code, make_regular)

    figure(exp_code)
    clf()
    subplot(3,2,1)
    plot(t,m)

    ylabel("algae")
    subplot(3,2,2)
    plot(t,b)
    ylabel("rotifers")
    subplot(3,2,3)
    plot(t,er)
    ylabel("egg-ratio")
    subplot(3,2,4)
    plot(t,e)
    ylabel("eggs")
    subplot(3,2,5)
    plot(t,d)
    ylabel("dead")
    xlabel("time")
    subplot(3,2,6)
    plot(t,s)
    ylabel("external nutrients")
    xlabel("time")

end

test()

max_stars_repo_name: berndblasius/WaveletAnalysis

max_issues_repo_name: berndblasius/WaveletAnalysis

max_forks_repo_name: berndblasius/WaveletAnalysis


text:
########################################################################
# Uncertain values constructed from distributions with known parameters
########################################################################
o1 = UncertainValue(Normal, 0, 0.5)
o2 = UncertainValue(Normal, 2.0, 0.1)
o3 = UncertainValue(Uniform, 0, 4)
o4 = UncertainValue(Uniform, rand(100))
o5 = UncertainValue(Beta, 4, 5)
o6 = UncertainValue(Gamma, 4, 5)
o7 = UncertainValue(Frechet, 1, 2)
o8 = UncertainValue(BetaPrime, 1, 2)
o9 = UncertainValue(BetaBinomial, 10, 3, 2)
o10 = UncertainValue(Binomial, 10, 0.3)

@test resample(o1) isa Float64
@test resample(o2) isa Float64
@test resample(o3) isa Float64
@test resample(o4) isa Float64
@test resample(o5) isa Float64
@test resample(o6) isa Float64
@test resample(o7) isa Float64
@test resample(o8) isa Float64
@test resample(o9) isa Int
@test resample(o10) isa Int

########################################################################
# Resampling uncertain values with constraints
########################################################################

u = UncertainValue(Normal, 0, 1)

# Resampling many times
@test maximum(resample(u, TruncateLowerQuantile(0.2), 1000)) >= quantile(u.distribution, 0.2)
@test maximum(resample(u, TruncateUpperQuantile(0.7), 1000)) <= quantile(u.distribution, 0.7)
@test maximum(resample(u, TruncateMinimum(-0.1), 10)) >= -0.1
@test maximum(resample(u, TruncateMaximum(0.1), 10)) <= 0.1

# quantile truncation
c = resample(u, TruncateQuantiles(0.3, 0.9), 1000)
@test minimum(c) >= quantile(u.distribution, 0.3)
@test maximum(c) <= quantile(u.distribution, 0.9)

# standard deviation truncation
c = resample(u, TruncateStd(1))
@test maximum(c) <= 0 + std(u)
@test minimum(c) >= 0 - std(u)

# range truncation
c = resample(u, TruncateRange(-0.2, 0.2))
@test maximum(c) <= 0.2
@test minimum(c) >= -0.2

# range truncation
c = resample(u, TruncateRange(-0.2, 0.2), 100)
@test maximum(c) <= 0.2
@test minimum(c) >= -0.2

# Resampling once
using StatsBase
@test maximum(resample(u, TruncateLowerQuantile(0.2))) >= quantile(u.distribution, 0.2)
@test maximum(resample(u, TruncateUpperQuantile(0.7))) <= quantile(u.distribution, 0.7)
@test maximum(resample(u, TruncateMinimum(-0.1))) >= -0.1
@test maximum(resample(u, TruncateMaximum(0.1))) <= 0.1
c = resample(u, TruncateQuantiles(0.3, 0.9))
@test minimum(c) >= quantile(u.distribution, 0.3)
@test maximum(c) <= quantile(u.distribution, 0.9)


text:
module ChainRules

using Base.Broadcast: materialize, materialize!, broadcasted

include("markup.jl")
include("interface.jl")
include("rules.jl")

end # module


text:
module Convection
export flux,minf,maxf
function flux(x::Float64)
    x
end
function minf(a::Float64,b::Float64)
    if b<=0.0
        return flux(b)
    elseif a>=0.0
        return flux(a)
    else
        return 0.0
    end
end
function maxf(a::Float64,b::Float64)
    max(flux(a),flux(b))
end
end

max_stars_repo_name: danielcort/benchmarks-python-julia-c

max_issues_repo_name: danielcort/benchmarks-python-julia-c

max_forks_repo_name: danielcort/benchmarks-python-julia-c


text:
using InterpretableValidation
using Test
using Distributions
using Random

## Test construction of ConstrainedTimeseriesDistribution
x = collect(1:100.)
uniform = ConstrainedTimeseriesDistribution(IID(x, Uniform(2,3)))
v = rand(uniform)
@test all(v.>2) && all(v.<3)

normal = ConstrainedTimeseriesDistribution(IID(x, Normal(0.5, 1.)))
normal.lb = uniform.lb
v = rand(normal)
@test all(v.>2)

categorical = ConstrainedTimeseriesDistribution(IID(x, Categorical(5)))
v = rand(categorical)
@test all([v[i] in collect(1:5) for i=1:length(v)])

bernoulli = ConstrainedTimeseriesDistribution(IID(x, Bernoulli(0.1)))
v = rand(bernoulli)
@test sum(v) > 0
@test sum(v) < 20

gp = ConstrainedTimeseriesDistribution(GaussianProcess(m = (x) -> 0, k = squared_exp_kernel(l=4), x = x))
gp.lb = -1*ones(100)
v = rand(gp)
@test all(v .> -1)

@test N_pts(gp) == 100

## Test discrete vs. continuous
@test !isdiscrete(uniform)
@test isdiscrete(categorical)
@test !iscontinuous(categorical)
@test isdiscrete(bernoulli)
@test iscontinuous(gp)

## Test feasible check
@test isfeasible(uniform)
@test isfeasible(categorical)

old = categorical.feasible[10]
categorical.feasible[10] = Int64[]
@test !isfeasible(categorical)
categorical.feasible[10] = old

old = uniform.lb[1]
uniform.lb[1]=100
@test !isfeasible(uniform)
uniform.lb[1]=old

## Test construction of MvTimeseriesDistribution
mvts = MvTimeseriesDistribution(:x => uniform, :y =>normal, :z => categorical, :gp =>gp, :b => bernoulli)
@test N_pts(mvts) == 100

v = rand(mvts)
@test length(v) == 5

## Test the application of constraints
uniform.lb = -Inf*ones(100)
uniform.ub = Inf*ones(100)
greaterthan!(uniform, 3., [fill(true, 50)..., fill(false, 50)...])
@test all(uniform.lb[1:50] .== 3.0)
@test all(uniform.ub[51:100] .== 3.0)

greaterthan!(uniform, 2., [fill(true, 50)..., fill(false, 50)...])
@test all(uniform.lb[1:50] .== 3.0)
@test all(uniform.ub[51:100] .== 2.0)

greaterthan!(uniform, 4., [fill(true, 50)..., fill(false, 50)...])
@test all(uniform.lb[1:50] .== 4.0)
@test all(uniform.ub[51:100] .== 2.0)

greaterthan!(uniform, 5., [fill(true, 50)..., fill(:anybool, 50)...])
@test all(uniform.lb[1:50] .== 5.0)
@test all(uniform.ub[51:100] .== 2.0)

lessthan!(uniform, 4., [fill(true, 50)..., fill(false, 50)...])
@test all(uniform.ub[1:50] .== 4.0)
@test all(uniform.lb[51:100] .== 4.0)

uniform.lb = -Inf*ones(100)
uniform.ub = Inf*ones(100)
continuous_equality!(uniform, 3., [fill(true, 50)..., fill(:anybool, 50)...])
@test all(uniform.lb[1:50] .== uniform.ub[1:50])
@test all(rand(uniform)[1:50] .== 3.)

discrete_equality!(categorical, 3, [:anybool, fill(true, 50)..., fill(false, 49)...])
@test categorical.feasible[1] == [1,2,3,4,5]
@test all([categorical.feasible[i] == [3] for i=2:51])
@test all([!(3 in categorical.feasible[i]) for i in 52:100])

discrete_equality!(bernoulli, 1, [:anybool, fill(true, 50)..., fill(false, 49)...])
@test bernoulli.feasible[1] == [0,1]
@test all([bernoulli.feasible[i] == [1] for i=2:51])
@test all([bernoulli.feasible[i] == [0] for i=52:100])

## Test constraint on the time series
constrain_timeseries!(mvts, Meta.parse("x .<= 1."), Array{Any}(fill(true, 100)))
@test mvts[:x].ub == ones(100)

constraints = [(Meta.parse("y .<= 1."), Array{Any}(fill(true, 100)))]
constrain_timeseries!(mvts, constraints)
@test mvts[:y].ub == ones(100)

## Test the sampling based on expression
for (key, val) in mvts
    if iscontinuous(val)
        val.lb = -Inf*ones(100)
        val.ub = Inf*ones(100)
    elseif isdiscrete(val)
        val.feasible = fill([1,2,3,4,5], 100)
    end
end
v = Base.rand(Random.GLOBAL_RNG, Meta.parse("all(x .>= 1.) && all(y .>= 1.) && all(z .== 2)"), mvts)
@test all(v[:x] .> 1.)
@test all(v[:y] .> 1.)
@test all(v[:z] .== 2)

@test_throws InfeasibleConstraint Base.rand(Random.GLOBAL_RNG, Meta.parse("all(x .>= 1.) && all(x .<= 0.5)"), mvts)

mvts = MvTimeseriesDistribution(:x => IID(x, Uniform()), :y => IID(x, Normal()), :z => IID(x, Categorical(5)))

function target_expr(data::Dict{Symbol, Array})
    truth_vals = ((data[:x] .> 0.75) .& (data[:y] .< -0.75)) .| ((data[:x] .< 0.25) .& (data[:y] .> -0.25))
    truth_vals[1] = truth_vals[1] &&  (data[:z][1] == 5)
    return (length(truth_vals) - sum(truth_vals)) / length(truth_vals)
end

good_expr = Meta.parse("all(((x .>= 0.75) .& (y .<= -0.75)) .| ((x .<= 0.25) .& (y .>= -0.25)))")
good_vals = rand(Random.GLOBAL_RNG, good_expr, mvts)
@test target_expr(good_vals) <= 0.01

expr = Meta.parse("all_between((x .== 2.881512049363103) .| (x .<= 2.3141385919362447), 59, 14)")
constraints = sample_constraints(expr, 100, Random.GLOBAL_RNG)
mvts = MvTimeseriesDistribution(:x => IID(x, Uniform(2,3)))
good_vals = rand(Random.GLOBAL_RNG, expr, mvts)
x = good_vals[:x]
@test all_between((x .== 2.881512049363103) .| (x .<= 2.3141385919362447), 59, 14)


text:
immutable Lda
    c::Integer
    Astar
    data::Matrix
    kernel::Function
end

function project(lda::Lda, x::Vector, dim::Integer)
    n = size(lda.data, 2)
    Kt = Array(Float64, n)
    for i = 1:n
        Kt[i] = lda.kernel(lda.data[:, i], x)
    end
    return lda.Astar[:, 1:dim]' * Kt
end

function process(lda::Lda, data::Matrix, dim::Integer)
    n = size(data, 2)
    res = Array(Float64, dim, n)
    for i = 1:n
        res[:, i] = project(lda, data[:, i], dim)
    end
    return res
end

function lda(data::Matrix, labels::Vector, c::Integer = 2, kfun::Function = linear, karg = nothing)

    # Number of training samples
    n = size(data, 2)

    # Defines Kernel function
    if karg == nothing
        kernel(x, y) = kfun(x, y)
    else
        kernel(x, y) = kfun(x, y, karg)
    end

    # Separates c classes
    class = {}
    for i = 1:c
        push!(class, data[:, labels .== i])
    end

    # Compute Kj[m, n] = kernel(Xn, Xmj)
    K = {}
    for j = 1:c
        Kj = Array(Float64, n, size(class[j], 2))
        for i = 1:n
            for m = 1:size(class[j], 2)
                Kj[i, m] = kernel(data[:, i], class[j][:, m])
            end
        end
        push!(K, Kj)
    end

    # Compute Mi
    # Mij = 1/li * sum(k: 1 -> li) kernel(Xj, Xki)
    M = {}
    for i = 1:c
        Mi = Array(Float64, n)
        for j = 1:n
            for k = 1:size(class[i], 2)
                Mi[j] += kernel(data[:, j], class[i][:, k])
            end
        end
        Mi /= size(class[i], 2)
        push!(M, Mi)
    end

    # Compute N
    # sum(j = 1,2) Kj (I - 1lj) Kj'
    N = 0
    for j = 1:c
        ksize = size(K[j], 2)
        dsize = size(class[j], 2)
        N += K[j] * (eye(ksize) - ones(ksize, ksize) / dsize) * K[j]'
    end

    # Compute M*
    Mstar = Array(Float64, n)
    for j = 1:n
        for k = 1:n
            Mstar[k] = Mstar[k] + kernel(data[:, j], data[:, k])
        end
        Mstar[j] /= n
    end

    # M
    Mtot = 0
    for j = 1:c
        Mtot = Mtot + size(class[j], 2) * (M[j] - Mstar) * (M[j] - Mstar)'
    end

    println(Mstar)
    println(Mtot)
    Astar = eigvecs(inv(N + 0.1 *  eye(size(N, 1))) * Mtot)

    return Lda(c, Astar, data, kernel)
end


text:
using GHOST
using GHOST: @unpack
setup()
setup_parallel()
@unpack conn, schema = GHOST.PARALLELENABLER
data = execute(conn,
               String(read(joinpath(pkgdir(GHOST), "src", "assets", "sql", "branches_min_max.sql"))) |>
                   (obj -> replace(obj, "schema" => schema)) |>
                   (obj -> replace(obj, "min_lim" => 0)) |>
                   (obj -> replace(obj, "max_lim" => 100)),
               not_null = true) |>
    (obj -> getproperty.(obj, :branch))
time_start = now()

idx = (1:8:500)[6]
query_commits(view(data, idx:min(idx + 7, lastindex(data))), 100)

println(time_start)
# @sync @distributed for idx in 1:8:500
for idx in 1:8:500
    query_commits(view(data, idx:min(idx + 7, lastindex(data))), 100)
    sleep(0.25)
end
time_end = now()
canonicalize(CompoundPeriod(time_end - time_start))


text:
using DifferentialEquations
using Interpolations
using LinearAlgebra
using PyPlot
pygui(true)

function historymaker(si,τ,N)
    #itph = interpolate(si, BSpline(Constant()))#Constant, Linear, Quadratic, and Cubic corresponding to B-splines of degree 0, 1, 2, and 3 respectively.
    itph = interpolate(si, BSpline(Linear()))#Constant, Linear, Quadratic, and Cubic corresponding to B-splines of degree 0, 1, 2, and 3 respectively.
    #itph = interpolate(si, BSpline(Quadratic()))#Constant, Linear, Quadratic, and Cubic corresponding to B-splines of degree 0, 1, 2, and 3 respectively.
    #itph = interpolate(si, BSpline(Cubic(Line(OnGrid()))))#Constant, Linear, Quadratic, and Cubic corresponding to B-splines of degree 0, 1, 2, and 3 respectively.
# sitp = scale(itp, A_x)
  h(p, t; idxs) =itph((1.0-(-t/τ))*(N-1)+1.0)
  # h([1],-0.1,1)
  # t=-6:0.01:0
  # (tt->h(1,tt)).(t)
end

function delay_mathieu_model(du,u,h,p,t)
  δ, ϵ, b, κ, τ = p
  # println(h(p, t-τ))
  # println(u[:,1])
  # println(u[:,2])
  # println(du[:,1])
  # println(du[:,2])
  # du[:,1] .= - κ .* u[:,1] - (δ + ϵ * cos(t)) .* u[:,2] - b .* h(p, t-τ)
  # du[:,2] .= u[:,1]
  du[1] = - κ * u[1] - (δ + ϵ * cos(t)) * u[2] + b * h(p, t-τ, idxs=2)
  du[2] = u[1]
end

function mappinggeneration(Si,p)
    # figure(3)
    # clf()
    Vi = zero(Si)
    τ=p[5]
    N=size(Si,1)
    for i in 1:size(Si,2)

        hist = historymaker(Si[:,i],τ,N)

        # p = [1.0,2.0,1,0.1]

        T = 2 * pi
        lags = [τ]
        tspan = (0.0, T)
        Δt = 1e-6
        u0 = [(hist(p, 0, idxs = 2) - hist(p, -Δt, idxs = 2)) / Δt, hist(p, 0, idxs = 2)]
        prob = DDEProblem(delay_mathieu_model, u0, hist, tspan, p; constant_lags = lags)

        alg = MethodOfSteps(Tsit5())
        sol = solve(prob, alg,reltol=1e-3,abstol=1e-3)#,reltol=1e-8,abstol=1e-8)

        tplot = LinRange(T - τ, T, N)
        solplot = [sol(ttt)[2] for ttt in tplot]
        Vi[:,i] .= solplot

        # figure(3)
        # # clf()
        # tplot = LinRange(-τ, 0, N)
        # solplot = [hist(p, ttt, idxs = 1) for ttt in tplot]
        # plot(tplot, real.(solplot), marker="o",color="green", markerfacecolor="blue", markersize=4)
        #
        # # tplot = LinRange(0, T - τ, N)
        # # solplot = [sol(ttt)[1] for ttt in tplot]
        # # plot(tplot, real.(solplot), marker="o",color="green", markerfacecolor="blue", markersize=4)
        #
        # tplot = LinRange(T - τ, T, N)
        # solplot = [sol(ttt)[2] for ttt in tplot]
        # plot(tplot, real.(solplot), marker="o",color="green", markerfacecolor="blue", markersize=4)

    end

    return Vi

end

function delayoscillmapp(δ, b)

# initialization of the Eigen vector
N = 50
eigvecnum = 6
# δ, ϵ, b, κ, τ = [2.0,1.0,0.02,0.0,2*pi] # [1.0,2.0,1,0.1]
# p=[δ, ϵ, b, κ, τ]
p=[δ,1.0,b,0.0,2*pi]
#
Si = rand(N, eigvecnum)
eigsol_inter=Array{Array{ComplexF64}}(undef,0)
#@time begin
for k in 1:6
    #global(Si)
    Vi=mappinggeneration(Si,p)
    Hi=Si\Vi
    # figure(9)
    # clf()
    # surf(abs.(Hi))
    Eigs=eigen(Hi)
    Si=Vi*Eigs.vectors
    for iloc in 1:eigvecnum
        sinorminv=1/norm(Si[:,iloc])
        Si[:,iloc] *= sinorminv
    end
    push!(eigsol_inter,Eigs.values)

    # figure(14)
    # plot(k*ones(size(Eigs.values)),abs.(Eigs.values), linestyle="none", marker="o",color="green", markerfacecolor="blue", markersize=4)

    #figure(4)
    # clf()
    #plot(real(Eigs.values),imag(Eigs.values), linestyle="none", marker="o",color="green", markerfacecolor="blue", markersize=4)
    # fi=0:0.01:2*pi
    # plot(sin.(fi),cos.(fi))
    #pause(0.0001)
end

#end

# figure(16)
# clf()
# for keig in 1:min(10,eigvecnum)#1:2
#     scatter(1:length(eigsol_inter),log.(abs.([abs(e[keig])-abs(eigsol_inter[end][keig]) for e in eigsol_inter])) ./ log(10))
#     #scatter(1:length(eigsol_inter),log.(abs.([e[keig] for e in eigsol_inter])))
# end
#
#
# figure(4)
# clf()
# plot(real(eigsol_inter[end]),imag(eigsol_inter[end]), linestyle="none", marker="o",color="green", markerfacecolor="blue", markersize=4)
# fi=0:0.01:2*pi
# plot(sin.(fi),cos.(fi))
# pause(0.1)

return abs(eigsol_inter[end][1])-1

end
delayoscillmapp(2.0, 0.02)

using MDBM

using PyPlot;
pygui(true);

mymdbm=MDBM_Problem(delayoscillmapp,[-1:0.25:9,-1.5:0.25:1.5])

for kinter in 1:4
MDBM.solve!(mymdbm,1)

x_eval,y_eval=getevaluatedpoints(mymdbm)

x_sol,y_sol=getinterpolatedsolution(mymdbm)

fig = figure(112);clf()
scatter(x_eval,y_eval,s=5)
scatter(x_sol,y_sol,s=5)
pause(0.01)
end


text:
# User defined variable

files = ["src","test"]
interval = 0.25
quotes = readdlm(Pkg.dir("JulieTest/res/quotes.txt"), '\n')
reporter = ""

# end

# do not polute global space with function definition
function julieTest()
  global files
  global quotes
  global interval
  global reporter

  # we need to call runTests method from this module
  require(joinpath(Pkg.dir(),"JulieTest/src/JulieTest.jl"))

  lastRun = time()
  running = false
  watcher = {}
  PWD = pwd()
  currentFile = "testConfig.jl"

  currentFile in files || push!(files, currentFile)

  function line()
    "\n" * "=" ^ Base.tty_size()[2] * "\n"
  end

  function surround(txt)
    line() * txt * line()
  end

  function isTestFile(str::String)
    ismatch(r"^([^/]+/)*[^_]+_test\.jl$",str)
  end

  function getQuotes()
    quotes[floor(rand() * length(quotes))]
  end

  function runAllTest()
    for file in readdir(joinpath(PWD,"test"))
      isTestFile(file) || continue
      reload(joinpath(PWD, "test", file))
    end
    JulieTest.runTests()
  end

  function runSingleTest(filepath::String)
    reload(filepath)
    JulieTest.runTests()
  end

  function runTestFor(filepath::String)
    # Try to find and run the test file of the given file
    if ismatch(r"src/.*\.jl", filepath)
      testpath = joinpath("test", replace(basename(filepath), r"\.jl$", "_test.jl"))
      isfile(testpath) && return runSingleTest(testpath)
    elseif isTestFile(filepath)
      return runSingleTest(filepath)
    end
    # If no test file found for this file just run all test file
    runAllTest()
  end

  function watch(dir::String, isDir::Bool)
    # watch recursively
    if isDir
      for file in readdir(dir)
        file = joinpath(dir,file)
        isdir(file) && push!(watcher, watch(file,true))
      end
    end

    watch_file(dir) do filename,event,status
      (running || time() - lastRun < interval) && return
      running = true
      if isDir
        filename =  dir == "." ? filename : joinpath(dir,filename)
      else
        filename = dir
      end
      print_with_color(:green, filename, " "); println("has modified")
      try
        filename == currentFile ? restart() : runTestFor(filename)
      catch e
        print("\033[91m"); dump(e); print("\033[0m")
      end
      running, lastRun = false, time()
    end
  end

  function watch(dir::String)
    watch(dir, isdir(dir))
  end

  function startWatch()
    watcher = [watch(file) for file in files]
  end

  function restart()
    println("\033[91m---Reloading configuration---\033[0m")
    [close(file) for file in watcher]
    reload(currentFile)
  end

  print(surround(getQuotes()))

  for arg in ARGS
    if ismatch(r"--reporter=.+",arg)
      reporter = match(r"(?<=--reporter=).+", arg).match
    end
  end

  JulieTest.reporter(reporter)

  "--single" in ARGS  && return Base.exit(runAllTest())
  "--skip-init" in ARGS || runAllTest()

  startWatch()
  while true; sleep(typemax(Int32)); end

end


text:
function select(db::Database, elnames::Vector{<:AbstractString}, phasids::Vector{Int})
    elems = Element[]
    funcs = GFunction[]
    phass = Phase[]

    for elem in db.elems
        if elem.name in vcat(elnames, ["/-", "Va"])
            push!(elems, elem)
        end
    end

    K = length(db.phass)
    for k in 1:K
        !in(k, phasids) && continue
        phas_selected = select(db.phass[k], elnames)
        isnothing(phas_selected) && continue
        push!(phass, phas_selected)
    end

    # Functions used in parameters
    for func in db.funcs
        if any(phas -> isused(func, phas), phass)
            push!(funcs, func)
        end
    end

    # Functions used in functions
    function push_func!(func)
        func in funcs && return false
        if any(_func -> isused(func, _func), funcs)
            push!(funcs, func)
            return true
        else
            return false
        end
    end
    found = true
    while found
        found = any(push_func!, db.funcs)
    end
    return Database(elems, funcs, phass, db.types)
end

function select(db::Database, elnames::AbstractString = "*", args::Union{Int,Function}...)
    if elnames == "*"
        elnames = filter(e -> !in(e, ["/-", "Va"]), map(e -> e.name, db.elems))
    else
        elnames = split(elnames)
    end
    phasids = Int[]
    exclude = Int[]
    selectall = false
    for arg in args
        if arg == *
            selectall = true
            continue
        elseif typeof(arg) == Int
            arg > 0 ? push!(phasids, arg) : push!(exclude, -arg)
        else
            @error "Unsupported argument: $arg"
        end
    end
    if selectall
        K = length(db.phass)
        phasids = vcat(phasids, [ k for k in 1:K if !in(k, exclude) ])
    end
    select(db, elnames, phasids)
end

function select(phas::Phase, elnames::Vector{<:AbstractString})
    cons, nonempty = select(phas.cons, elnames)
    [] in cons && return nothing

    sites = [ phas.sites[s] for s in nonempty ]

    params = Parameter[]
    for param in phas.params
        if iscomposedof(param.comb, vcat(elnames, "Va"))
            comb = [ param.comb[s] for s in nonempty ]
            push!(params, Parameter(param.name,
                                    param.symbol,
                                    comb,
                                    param.order,
                                    param.temp,
                                    param.temps,
                                    param.exprs,
                                    param.funcname,
                                    param.funcstr))
        end
    end

    return Phase(phas.name, phas.state, phas.type, sites, cons, params)
end

function select(cons::Constitution, elnames::Vector{<:AbstractString})
    latts = Sublattice[]
    nonempty = Int[]
    S = length(cons)
    for s in 1:S
        latt = filter(spec -> iscomposedof(spec, vcat(elnames, "Va")), cons[s])
        latt == [] || latt == ["Va"] && continue
        push!(latts, latt)
        push!(nonempty, s)
    end
    return latts, nonempty
end

function Base.merge!(db_master::Database, db_source::Database,
                     args::Union{Pair{Int,Int},Int,Function}...)
    # Elements
    elemnames_master = map(elem -> getfield(elem, :name), db_master.elems)
    for elem in db_source.elems
        elem.name in elemnames_master && continue
        push!(db_master.elems, elem)
    end

    # Phases
    merged = Int[]
    for arg in args
        if typeof(arg) == Pair{Int,Int}
            ks = arg.first
            km = arg.second
            merge!(db_master.phass[km], db_source.phass[ks])
            push!(merged, ks)
        elseif typeof(arg) == Int
            if arg < 0
                push!(merged, -arg)
                continue
            end
            arg in merged && continue
            push!(db_master, db_source.phass[arg])
            push!(merged, arg)
        elseif arg == *
            K = length(db_source.phass)
            for k in 1:K
                k in merged && continue
                push!(db_master, db_source.phass[k])
                push!(merged, k)
            end
        else
            @error "unsupported type of argument: $arg"
        end
    end

    # Functions used in parameters
    for func in db_source.funcs
        if any(phas -> isused(func, phas), db_master.phass)
            push!(db_master.funcs, func)
        end
    end

    # Functions used in functions
    function push_func!(func)
        func in db_master.funcs && return false
        if any(_func -> isused(func, _func), db_master.funcs)
            push!(db_master.funcs, func)
            return true
        else
            return false
        end
    end
    found = true
    while found
        found = any(push_func!, db_source.funcs)
    end

    return nothing
end

function Base.merge!(phas_master::Phase, phas_source::Phase)
    if phas_master.sites ≠ phas_source.sites
        @error "Stoichiometry of both phases must match ($(phas_master.name))"
    end

    if length(phas_master.cons) ≠ length(phas_source.cons)
        @error "Constitution of both phases must match ($(phas_master.name))"
    end

    # Constitution
    S = length(phas_master.cons)
    for s in 1:S
        latt = phas_source.cons[s]
        for spec in latt
            if !(spec in phas_master.cons[s])
                push!(phas_master.cons[s], spec)
            end
        end
    end

    # Parameters
    paramnames_master = map(localname, phas_master.params)
    for param in phas_source.params
        name = localname(param)
        if name in paramnames_master
            @info "$(param.symbol)($name) already exists in $(phas_master.name)"
            continue
        end
        push!(phas_master.params, param)
    end
end

function Base.push!(db::Database, phas::Phase)
    push!(db.phass, phas)
end

function Base.merge(db_master::Database, db_source::Database,
                    args::Union{Pair{Int,Int},Int,Function}...)
    db = deepcopy(db_master)
    merge!(db, db_source, args...)
    return db
end


text:
function solve()
    n, s, k = [parse(Int, x) for x in split(readline())]
    s = n-s
    d = gcd(n, s, k)
    n ÷= d
    s ÷= d
    k ÷= d
    if gcd(k, n) ≠ 1
        -1
    else
        i = invmod(k, n)
        mod(i*s, n)
    end
end

function main()
    t = parse(Int, readline())
    for _=1:t
        println(solve())
    end
end

main()

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
export PS_IDLE,
    PS_DND,
    PS_ONLINE,
    PS_OFFLINE

const PRESENCE_STATUSES = ["idle", "dnd", "online", "offline"]

"""
A [`User`](@ref)'s status sent in a [`Presence`](@ref). Available values are `PS_IDLE`,
`PS_DND`, `PS_ONLINE`, and `PS_OFFLINE`.
More details [here](https://discordapp.com/developers/docs/resources/channel#message-object-message-application-structure).
"""
@enum PresenceStatus PS_IDLE PS_DND PS_ONLINE PS_OFFLINE PS_UNKNOWN

function PresenceStatus(x::AbstractString)
    i = findfirst(s -> s == x, PRESENCE_STATUSES)
    return i === nothing ? PS_UNKNOWN : PresenceStatus(i - 1)
end

Base.string(x::PresenceStatus) = x === PS_UNKNOWN ? "" : PRESENCE_STATUSES[Int(x) + 1]
JSON.lower(x::PresenceStatus) = string(x)

"""
A [`User`](@ref)'s presence.
More details [here](https://discordapp.com/developers/docs/topics/gateway#presence-update).
"""
struct Presence
    user::User
    roles::Union{Vector{Snowflake}, Missing}
    game::Union{Activity, Nothing}
    guild_id::Union{Snowflake, Missing}
    status::PresenceStatus
    activities::Vector{Activity}
end
@boilerplate Presence :dict :docs :lower :merge


text:
using Genie.Configuration, Logging

const config = Settings(
    server_port=8000,
    server_host="127.0.0.1",
    log_level=Logging.Info,
    log_to_file=false,
    server_handle_static_files=true
)

ENV["JULIA_REVISE"] = "auto"

max_stars_repo_name: MLH-Fellowship/ZoomCord

max_issues_repo_name: MLH-Fellowship/ZoomCord

max_forks_repo_name: MLH-Fellowship/ZoomCord


text:
#!/usr/bin/env julia

# Calculate max count for some huffman tree depth.

function fibonacci(n)
	f = zeros(Int, n)
	f[1] = 1
	f[2] = 1
	for i in 3:n
		f[i] = f[i-1] + f[i-2]
	end
	f
end

if length(ARGS) ≠ 1
	println("""
Usage:
	huffman_tree_depth.jl <block_len_log2>

""")
	exit(1)
end

block_len_log2 = parse(Int, ARGS[1])
@show block_len_log2

block_len = 1 << block_len_log2

max_depth = 0
for d in 1:32
	f = fibonacci(d+1)
	max_freq = cumsum(f)[end]
	if max_freq >= block_len
		max_depth = d
		break
	end
end

@show max_depth


text:
using .NCDatasets

export NCDarray, NCDstack, NCDdimMetadata, NCDarrayMetadata, NCDstackMetadata

"""
    NCDdimMetadata(val::Dict)

[`Metadata`](@ref) wrapper for [`NCDarray`](@ref) dimensions.
"""
struct NCDdimMetadata{K,V} <: DimMetadata{K,V}
    val::Dict{K,V}
end

"""
    NCDarrayMetadata(val::Dict)

[`Metadata`](@ref) wrapper for [`NCDarray`](@ref) metadata.
"""
struct NCDarrayMetadata{K,V} <: ArrayMetadata{K,V}
    val::Dict{K,V}
end

"""
    NCDstackMetadata(val::Dict)

[`Metadata`](@ref) wrapper for [`NCDarray`](@ref) metadata.
"""
struct NCDstackMetadata{K,V} <: StackMetadata{K,V}
    val::Dict{K,V}
end

const UNNAMED_NCD_KEY = "unnamed"

# Utils ########################################################################

ncread(f, path::String) = NCDatasets.Dataset(f, path)

nondimkeys(dataset) = begin
    dimkeys = keys(dataset.dim)
    removekeys = if "bnds" in dimkeys
        dimkeys = setdiff(dimkeys, ("bnds",))
        boundskeys = map(k -> dataset[k].attrib["bounds"], dimkeys)
        union(dimkeys, boundskeys)
    else
        dimkeys
    end
    setdiff(keys(dataset), removekeys)
end

# Add a var array to a dataset before writing it.
ncwritevar!(dataset, A::AbstractGeoArray{T,N}) where {T,N} = begin
    A = reorder(A, ForwardIndex()) |> a -> reorder(a, ForwardRelation())
    if ismissing(missingval(A))
        # TODO default _FillValue for Int?
        fillvalue = get(metadata(A), "_FillValue", NaN)
        A = replace_missing(A, convert(T, fillvalue))
    end
    # Define required dim vars
    for dim in dims(A)
        key = lowercase(string(name(dim)))
        haskey(dataset.dim, key) && continue

        if dim isa Lat || dim isa Lon
            dim = convertmode(Converted, dim)
        end
        index = ncshiftindex(dim)
        md = metadata(dim)
        attribvec = [] #md isa Nothing ? [] : [val(md)...]
        defDim(dataset, key, length(index))
        println("        key: \"", key, "\" of type: ", eltype(index))
        defVar(dataset, key, Vector(index), (key,); attrib=attribvec)
    end
    # TODO actually convert the metadata type
    attrib = if metadata isa NCDarrayMetadata
        deepcopy(val(metadata(A)))
    else
        Dict()
    end
    # Remove stack metdata if it is attached
    pop!(attrib, "dataset", nothing)
    # Set missing value
    if !ismissing(missingval(A))
        try
            fv = convert(T, missingval(A))
            attrib["_FillValue"] = fv
        catch
            @warn "`missingval` $(missingval(A)) was invalid for data of type $T."
        end
    end
    key = if string(name(A)) == ""
        UNNAMED_NCD_KEY
    else
        string(name(A))
    end
    println("        key: \"", key, "\" of type: ", T)
    dimnames = lowercase.(string.(name.(dims(A))))
    attribvec = [attrib...]
    var = defVar(dataset, key, eltype(A), dimnames; attrib=attribvec)

    var[:] = data(A)
end

ncshiftindex(dim::Dimension) = ncshiftindex(mode(dim), dim)
ncshiftindex(mode::AbstractSampled, dim::Dimension) = val(dim)
# As with plotting, we really should shift the index
# to `Center`, but reprojection introduces errors
    # if span(mode) isa Regular
    #     # that are tricky to handle currently.
    #     if dim isa TimeDim
    #         if eltype(dim) isa Dates.AbstractDateTime
    #             val(dim)
    #         else
    #             shiftindexloci(Center(), dim)
    #         end
    #     else
    #         shiftindexloci(Center(), dim)
    #     end
    # else
    #     dim
    # end |> val

ncshiftindex(mode::IndexMode, dim::Dimension) = val(dim)

# CF standards don't enforce dimension names.
# But these are common, and should take care of most dims.
const dimmap = Dict("lat" => Lat,
                    "latitude" => Lat,
                    "lon" => Lon,
                    "long" => Lon,
                    "longitude" => Lon,
                    "time" => Ti,
                    "lev" => Vert,
                    "mlev" => Vert,
                    "level" => Vert,
                    "vertical" => Vert,
                    "x" => X,
                    "y" => Y,
                    "z" => Z,
                   )

# Array ########################################################################
"""
    NCDarray(filename::AbstractString; name=nothing, refdims=(),
             dims=nothing, metadata=nothing, crs=EPSG(4326), dimcrs=EPSG(4326))

A [`DiskGeoArray`](@ref) that loads that loads NetCDF files lazily from disk.

The first non-dimension layer of the file will be used as the array. Dims are usually
detected as [`Lat`](@ref), [`Lon`](@ref), [`Ti`]($DDtidocs), and [`Vert`] or
possibly `X`, `Y`, `Z` when detected. Undetected dims will use the generic `Dim{:name}`.

This is an incomplete implementation of the NetCDF standard. It will currently
handle simple files in lattitude/longitude projections, or projected formats
if you manually specify `crs` and `dimcrs`. How this is done may also change in
future, including detecting and converting the native NetCDF projection format.

## Arguments

- `filename`: `String` pointing to a netcdf file.

## Keyword arguments

- `crs`: defaults to lat/lon `EPSG(4326)` but may be any `GeoFormat` like `WellKnownText`.
  If the underlying data is in a different projection `crs` will need to be set to allow
  `write` to a different file format. In future this may be detected automatically.
- `dimcrs`: The crs projection actually present in the Dimension index `Vector`, which
  may be different to the underlying projection. Defaults to lat/lon `EPSG(4326)` but
  may be any crs `GeoFormat`.
- `name`: `Symbol` name for the array. Will use array key if not supplied.
- `dims`: `Tuple` of `Dimension`s for the array. Detected automatically, but can be passed in.
- `refdims`: `Tuple of` position `Dimension`s the array was sliced from.
- `missingval`: Value reprsenting missing values. Detected automatically when possible, but
  can be passed it.
- `metadata`: [`Metadata`](@ref) object for the array. Detected automatically as
  [`NCDarrayMetadata`](@ref), but can be passed in.

## Example

```julia
A = NCDarray("folder/file.ncd")
# Select Australia from the default lat/lon coords:
A[Lat(Between(-10, -43), Lon(Between(113, 153)))
```
"""
struct NCDarray{T,N,A,D<:Tuple,R<:Tuple,Na<:Symbol,Me,Mi,S,K
               } <: DiskGeoArray{T,N,D,LazyArray{T,N}}
    filename::A
    dims::D
    refdims::R
    name::Na
    metadata::Me
    missingval::Mi
    size::S
    key::K
end
NCDarray(filename::AbstractString, key...; kwargs...) = begin
    isfile(filename) || error("File not found: $filename")
    ncread(dataset -> NCDarray(dataset, filename, key...; kwargs...), filename)
end
NCDarray(dataset::NCDatasets.Dataset, filename, key=nothing;
         crs=EPSG(4326),
         dimcrs=EPSG(4326),
         name=nothing,
         dims=nothing,
         refdims=(),
         metadata=nothing,
         missingval=missing,
        ) = begin
    keys_ = nondimkeys(dataset)
    key = (key isa Nothing || !(string(key) in keys_)) ? first(keys_) : string(key) |> Symbol
    var = dataset[string(key)]
    dims = dims isa Nothing ? GeoData.dims(dataset, key, crs, dimcrs) : dims
    name = Symbol(name isa Nothing ? key : name)
    metadata_ = metadata isa Nothing ? GeoData.metadata(var, GeoData.metadata(dataset)) : metadata
    size_ = map(length, dims)
    T = eltype(var)
    N = length(dims)

    NCDarray{T,N,typeof.((filename,dims,refdims,name,metadata_,missingval,size_,key))...
       }(filename, dims, refdims, name, metadata_, missingval, size_, key)
end

key(A::NCDarray) = A.key

# AbstractGeoArray methods

withsourcedata(f, A::NCDarray) =
    ncread(dataset -> f(dataset[string(key(A))]), filename(A))

# Base methods

Base.size(A::NCDarray) = A.size

"""
    Base.write(filename::AbstractString, ::Type{NCDarray}, s::AbstractGeoArray)

Write an NCDarray to a NetCDF file using NCDatasets.jl

Returns `filename`.
"""
Base.write(filename::AbstractString, ::Type, A::AbstractGeoArray) = begin
    meta = metadata(A)
    # if meta isa Nothing
    # else
        # Remove the dataset metadata
        # stackmd = pop!(deepcopy(val(meta)), "dataset", Dict())
    #    dataset = NCDatasets.Dataset(filename, "c"; attrib=stackmd)
    # end
    dataset = NCDatasets.Dataset(filename, "c")
    try
        println("    Writing netcdf...")
        ncwritevar!(dataset, A)
    finally
        close(dataset)
    end
    return filename
end

# Stack ########################################################################

"""
    NCDstack(filenames; keys, kwargs...)
    NCDstack(filenames...; keys, kwargs...)
    NCDstack(files::NamedTuple; refdims=(), window=(), metadata=nothing, childkwargs=())
    NCDstack(filename::String; refdims=(), window=(), metadata=nothing, childkwargs=())

A lazy [`AbstractGeoStack`](@ref) that uses NCDatasets.jl to load NetCDF files.
Can load a single multi-layer netcdf file, or multiple single-layer netcdf
files. In multi-file mode it returns a regular `GeoStack` with a `childtype`
of [`NCDarray`](@ref).

Indexing into `NCDstack` with layer keys (`Symbol`s) returns a [`GeoArray`](@ref).
Dimensions are usually detected as [`Lat`](@ref), [`Lon`](@ref), [`Ti`]($DDtidocs),
and [`Vert`] or `X`, `Y`, `Z` when detected. Undetected dims use the generic `Dim{:name}`.

# Arguments

- `filename`: `Tuple` or `Vector` or splatted arguments of `String`,
  or single `String` path, to NetCDF files.

# Keyword arguments

- `keys`: Used as stack keys when a `Tuple`, `Vector` or splat of filenames are passed in.
  These default to the first non-dimension data key in each NetCDF file.
- `refdims`: Add dimension position array was sliced from. Mostly used programatically.
- `window`: A `Tuple` of `Dimension`/`Selector`/indices that will be applied to the
  contained arrays when they are accessed.
- `metadata`: A [`StackMetadata`](@ref) object.
- `childkwargs`: A `NamedTuple` of keyword arguments to pass to the
  [`NCDarray`](@ref) constructor.

# Examples

```julia
stack = NCDstack(filename; window=(Lat(Between(20, 40),))
# Or
stack = NCDstack([fn1, fn1, fn3, fn4])
# And index with a layer key
stack[:soiltemp]
```
"""
struct NCDstack{T,R,W,M,K} <: DiskGeoStack{T}
    filename::T
    refdims::R
    window::W
    metadata::M
    childkwargs::K
end
NCDstack(filename::AbstractString;
         refdims=(),
         window=(),
         metadata=ncread(metadata, filename),
         childkwargs=()) =
    NCDstack(filename, refdims, window, metadata, childkwargs)
# These actually return a GeoStack
NCDstack(filenames...; kwargs...) = NCDstack(filenames; kwargs...)
NCDstack(filenames::Union{Tuple{AbstractString,Vararg},Vector{AbstractString}},
         keys=ncfilenamekeys(filenames); kwargs...) =
    NCDstack(NamedTuple{keys}(filenames); kwargs...)
NCDstack(files::NamedTuple;
         refdims=(),
         window=(),
         metadata=nothing,
         childkwargs=()) =
    GeoStack(NamedTuple{keys}(filenames), refdims, window, metadata,
             childtype=NCDarray, childkwargs)

ncfilenamekeys(filenames) =
    (ncread(ds -> first(nondimkeys(ds)), fn) for fn in filenames) |> cleankeys

childtype(::NCDstack) = NCDarray
childkwargs(stack::NCDstack) = stack.childkwargs
crs(stack::NCDarray) = get(childkwargs(stack), :crs, EPSG(4326))
dimcrs(stack::NCDarray) = get(childkwargs(stack), :dimcrs, nothing)

# AbstractGeoStack methods
withsource(f, ::Type{NCDarray}, path::AbstractString, key=nothing) = ncread(f, path)
withsourcedata(f, ::Type{NCDarray}, path::AbstractString, key) =
    ncread(d -> f(d[string(key)]), path)

# Override the default to get the dims of the specific key,
# and pass the crs and dimscrs from the stack
dims(stack::NCDstack, dataset, key::Key) =
    dims(dataset, key, crs(stack), dimcrs(stack))

missingval(stack::NCDstack) = missing

# Base methods

Base.keys(stack::NCDstack{<:AbstractString}) =
    cleankeys(ncread(nondimkeys, getsource(stack)))

"""
    Base.write(filename::AbstractString, ::Type{NCDstack}, s::AbstractGeoStack)

Write an NCDstack to a single netcdf file, using NCDatasets.jl.

Currently [`DimMetadata`](@ref) is not handled, and [`ArrayMetadata`](@ref)
from other [`AbstractGeoArray`](@ref) @types is ignored.
"""
Base.write(filename::AbstractString, ::Type{<:NCDstack}, s::AbstractGeoStack) = begin
    dataset = NCDatasets.Dataset(filename, "c"; attrib=val(metadata(s)))
    try
        map(key -> ncwritevar!(dataset, s[key]), keys(s))
    finally
        close(dataset)
    end
end

# DimensionalData methods for NCDatasets types ###############################

dims(dataset::NCDatasets.Dataset, key::Key, crs=nothing, dimcrs=nothing) = begin
    v = dataset[string(key)]
    dims = []
    for (i, dimname) in enumerate(NCDatasets.dimnames(v))
        if haskey(dataset, dimname)
            dvar = dataset[dimname]
            # Find the matching dimension constructor. If its an unknown name use
            # the generic Dim with the dim name as type parameter
            dimtype = get(dimmap, dimname, Dim{Symbol(dimname)})
            index = dvar[:]
            meta = NCDdimMetadata(Dict{String,Any}(dvar.attrib))
            mode = _ncdmode(index, dimtype, crs, dimcrs, meta)

            # Add the dim containing the dimension var array
            push!(dims, dimtype(index, mode, meta))
        else
            # The var doesn't exist. Maybe its `complex` or some other marker,
            # so make it a custom `Dim` with `NoIndex`
            push!(dims, Dim{Symbol(dimname)}(1:size(v, i), NoIndex(), nothing))
        end
    end
    (dims...,)
end

_ncdmode(index::AbstractArray{<:Number}, dimtype, crs, dimcrs, metadata) = begin
    # Assume the locus is at the center of the cell if boundaries aren't provided.
    # http://cfconventions.org/cf-conventions/cf-conventions.html#cell-boundaries
    # Unless its a time dimension.
    order = _ncdorder(index)
    span = _ncdspan(index, order)
    sampling = Intervals(Center())
    if dimtype in (Lat, Lon)
        Converted(order, span, sampling, crs, dimcrs)
    else
        Sampled(order, span, sampling)
    end
end
_ncdmode(index::AbstractArray{<:Dates.AbstractTime}, dimtype, crs, dimcrs, metadata) = begin
    order = _ncdorder(index)
    span, sampling  = _get_period(index, metadata)
    Sampled(order, span, sampling)
end
_ncdmode(index, dimtype, crs, dimcrs, mode) = Categorical()

_ncdorder(index) = index[end] > index[1] ? Ordered(ForwardIndex(), ForwardArray(), ForwardRelation()) :
                                           Ordered(ReverseIndex(), ReverseArray(), ForwardRelation())

_ncdspan(index, order) = begin
    # Handle a length 1 index
    length(index) == 1 && return Regular(zero(eltype(index)))

    step = index[2] - index[1]
    for i in 2:length(index) -1
        # If any step sizes don't match, its Irregular
        if !(index[i+1] - index[i] ≈ step)
            bounds = if length(index) > 1
                beginhalfcell = abs((index[2] - index[1]) * 0.5)
                endhalfcell = abs((index[end] - index[end-1]) * 0.5)
                if isrev(indexorder(order))
                    index[end] - endhalfcell, index[1] + beginhalfcell
                else
                    index[1] - beginhalfcell, index[end] + endhalfcell
                end
            else
                index[1], index[1]
            end
            return Irregular(bounds)
        end
    end
    # Otherwise regular
    return Regular(step)
end

function _get_period(index, mode)
    if haskey(mode, "delta_t")
        period = _parse_period(mode["delta_t"])
        period isa Nothing || return Regular(period), Points()
    elseif haskey(mode, "avg_period")
        period = _nc_parse_period(mode["avg_period"])
        period isa Nothing || return Regular(period), Intervals(Center())
    end
    return sampling = Irregular(), Points()
end

function _parse_period(period_str::String)
    regex = r"(\d\d\d\d)-(\d\d)-(\d\d) (\d\d):(\d\d):(\d\d)"
    mtch = match(regex, period_str)
    if mtch isa Nothing
        nothing
    else
        vals = Tuple(parse.(Int, mtch.captures))
        periods = (Year, Month, Day, Hour, Minute, Second)
        if length(vals) == length(periods)
            sum(map((p, v) -> p(v), periods, vals))
        end
    end
end

metadata(dataset::NCDatasets.Dataset) = NCDstackMetadata(Dict{String,Any}(dataset.attrib))
metadata(dataset::NCDatasets.Dataset, key::Key) = metadata(dataset[string(key)])
metadata(var::NCDatasets.CFVariable) = NCDarrayMetadata(Dict{String,Any}(var.attrib))
metadata(var::NCDatasets.CFVariable, stackmetadata::NCDstackMetadata) = begin
    md = metadata(var)
    md["dataset"] = stackmetadata
    md
end

missingval(var::NCDatasets.CFVariable) = missing

# Direct loading: better memory handling?
# readwindowed(A::NCDatasets.CFVariable) = readwindowed(A, axes(A)...)
# readwindowed(A::NCDatasets.CFVariable, i, I...) = begin
#     var = A.var
#     indices = to_indices(var, (i, I...))
#     shape = Base.index_shape(indices...)
#     dest = Array{eltype(var),length(shape)}(undef, map(length, shape)...)
#     NCDatasets.load!(var, dest, indices...)
#     dest
# end


text:
function get_realization(n::Integer, prep_fn = "prep.h5")
    habitat = load_realization(prep_fn, "habitat", n)
    movement = load_movement(prep_fn, "movement", n)
    init_pop = load_spatial_eq(prep_fn, "spatial_eq", n)
    comm_catchability = load_realization(prep_fn, "catchability_devs", n)
    log_qdevs = load_realization(prep_fn, "log_catchability_devs", n)

    habitat, movement, init_pop, comm_catchability, log_qdevs
end

struct SpatQSimPrep{M, P, Q, F}
    realization::Int
    movement::M
    init_pop::P
    log_qdevs::Q
    prep_file::F

    function SpatQSimPrep(realization::Int,
                          movement::M,
                          init_pop::P,
                          log_qdevs::Q,
                          prep_file::F = nothing) where {M<:MovementModel,
                                                         P<:PopState,
                                                         Q<:Matrix,
                                                         F<:Union{String, Nothing}}
        new{M, P, Q, F}(realization, movement, init_pop, log_qdevs, prep_file)
    end
end

function SpatQSimPrep(realization::Int, prep_file::String)
    _, movement, init_pop, _, log_qdevs =
        get_realization(realization, prep_file)
    SpatQSimPrep(realization, movement, init_pop, log_qdevs, prep_file)
end

movement(prep::SpatQSimPrep) = prep.movement
init_pop(prep::SpatQSimPrep) = prep.init_pop
log_qdevs(prep::SpatQSimPrep) = prep.log_qdevs
realization(prep::SpatQSimPrep) = prep.realization
prep_file(prep::SpatQSimPrep) = prep.prep_file


text:
module Utils
export
    # Traverse tree
    getchildindex,
    getparentindex,
    getleaf,
    getdepth,
    gettreelength,
    # Extraction
    getbasiscoef,
    getbasiscoefall,
    getrowrange,
    getcolrange,
    nodelength,
    coarsestscalingrange,
    finestdetailrange,
    # Metrics
    relativenorm,
    psnr,
    snr,
    ssim,
    # Datasets
    duplicatesignals,
    generatesignals,
    ClassData,
    generateclassdata

using
    Wavelets,
    LinearAlgebra,
    Distributions,
    Random,
    ImageQualityIndexes

# `maxtransformlevels` of a specific dimension
"""
    maxtransformlevels(n)

    maxtransformlevels(x[, dims])

Extension function from Wavelets.jl. Finds the max number of transform levels for an array
`x`. If `dims` is provided for a multidimensional array, then it finds the max transform
level for that corresponding dimension.

# Arguments
- `x::AbstractArray`: Input array.
- `dims::Integer`: Dimension used for wavelet transform.
- `n::Integer`: Length of input vector.

# Returns
`::Integer`: Max number of transform levels.

# Examples
```@repl
using Wavelets, WaveletsExt

# Define random signal
x = randn(64, 128)

# Max transform levels and their corresponding return values
maxtransformlevels(128)     # 7
maxtransformlevels(x)       # 6
maxtransformlevels(x, 2)    # 7
```
"""
function Wavelets.Util.maxtransformlevels(x::AbstractArray, dims::Integer)
    # Sanity check
    @assert 1 ≤ dims ≤ ndims(x)
    # Compute max transform levels of relevant dimension
    return maxtransformlevels(size(x)[dims])
end

"""
    getbasiscoef(Xw, tree)

Get the basis coefficients for the decomposed signal `Xw` with respect to the tree `tree`.

# Arguments
- `Xw::AbstractArray{T,2} where T<:Number`: Decomposed 1D-signal.
- `tree::BitVector`: The corresponding basis tree.

# Returns
- `::Array{T,1}`: Basis coefficients of signal.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate signal and wavelet
x = generatesignals(:heavysine)
wt = wavelet(WT.db4)

# Decompose signal
Xw = iwpd(x, wt)
tree = maketree(128, 6, :dwt)

# Get basis coefficients
xw = getbasiscoef(Xw, tree)
```
"""
function getbasiscoef(Xw::AbstractArray{T}, tree::BitVector) where T<:Number
    # Setup and Sanity Check
    N = ndims(Xw)
    @assert 2 ≤ N ≤ 3
    sz = size(Xw)[1:end-1]
    k = size(Xw)[end]
    x = N==2 ? view(Xw,:,1) : view(Xw,:,:,1)
    L = maxtransformlevels(x)
    @assert isvalidtree(x, tree)
    @assert k-1 ≤ L
    leaf = N==2 ? getleaf(tree, :binary) : getleaf(tree, :quad)
    leaf_len = N==2 ? gettreelength(1<<(L+1)) : gettreelength(1<<(L+1),1<<(L+1))
    @assert leaf_len == length(leaf)
    xw = Array{T}(undef, sz)

    # Extract basis coefficients
    for (i, isleaf) in enumerate(leaf)
        if isleaf
            d = N==2 ? getdepth(i,:binary) : getdepth(i, :quad) # Depth of node (0 for root node)
            d < k || throw(ArgumentError("Not enough decomposition levels in Xw."))
            if N==2
                nn = i-1<<d                     # Node number (0 for leftmost node)
                n₀ = nodelength(sz[1], d)           # Length of node
                rng = (nn*n₀+1):((nn+1)*n₀)
                @inbounds xw[rng] = Xw[rng, d+1]
            elseif N==3
                rng₁ = getrowrange(sz[1], i)
                rng₂ = getcolrange(sz[2], i)
                @inbounds xw[rng₁,rng₂] = Xw[rng₁,rng₂,d+1]
            end
        end
    end
    return xw
end

"""
    getbasiscoefall(Xw, tree)

Get the basis coefficients for all decomposed signals in `Xw` with respect to the tree(s)
`tree`.

# Arguments
- `Xw::AbstractArray{T,3} where T<:Number`: A set of decomposed 1D-signals.
- `tree::BitVector` or `tree::BitArray{2}`: The corresponding basis tree(s). If input is a
  `BitMatrix`, each column corresponds to a signal in `Xw`, and therefore the number of
  columns must be equal to the number of signals.

# Returns
- `::Array{T,2}`: Basis coefficients of signals.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate signals and wavelet
c = ClassData(:cbf, 10, 10, 10)
X = generateclassdata(c)
wt = wavelet(WT.db4)

# Decompose signals
Xw = iwpdall(X, wt)
tree = maketree(128, 6, :dwt)

# Get basis coefficients
getbasiscoefall(Xw, tree)
```
"""
function getbasiscoefall(Xw::AbstractArray{T}, tree::BitVector) where T<:Number
    # Setup and Sanity Check
    N = ndims(Xw)
    sz = size(Xw)[1:end-2]
    k = size(Xw)[end-1]
    m = size(Xw)[end]
    x = N==3 ? view(Xw,:,1,1) : view(Xw,:,:,1,1)
    @assert 3 ≤ N ≤ 4
    @assert k-1 ≤ maxtransformlevels(x)
    @assert isvalidtree(x, tree)
    nₜ = length(tree)
    @assert nₜ == gettreelength(sz...)
    leaf = N==3 ? getleaf(tree, :binary) : getleaf(tree, :quad)
    xw = Array{T}(undef, (sz...,m))

    # Extract basis coefficients
    for (i, isleaf) in enumerate(leaf)
        if isleaf
            d = N==3 ? getdepth(i,:binary) : getdepth(i, :quad) # Depth of node (0 for root node)
            d < k || throw(ArgumentError("Not enough decomposition levels in Xw."))
            if N==3
                nn = i-1<<d                     # Node number (0 for leftmost node)
                n₀ = nodelength(sz[1], d)           # Length of node
                rng = (nn*n₀+1):((nn+1)*n₀)
                @inbounds xw[rng,:] = @view Xw[rng, d+1,:]
            elseif N==4
                rng₁ = getrowrange(sz[1], i)
                rng₂ = getcolrange(sz[2], i)
                @inbounds xw[rng₁,rng₂,:] = @view Xw[rng₁,rng₂,d+1,:]
            end
        end
    end
    return xw
end

function getbasiscoefall(Xw::AbstractArray{T}, tree::BitArray{2}) where T<:Number
    # Setup and Sanity Check
    N = ndims(Xw)
    sz = size(Xw)[1:end-2]
    k = size(Xw)[end-1]
    m = size(Xw)[end]
    x = N==3 ? view(Xw,:,1,1) : view(Xw,:,:,1,1)
    @assert 3 ≤ N ≤ 4
    nₜ, mₜ = size(tree)
    @assert k-1 ≤ maxtransformlevels(x)
    @assert all( mapslices(treeᵢ -> isvalidtree(x,treeᵢ), tree, dims=1) )
    @assert m == mₜ
    @assert nₜ == gettreelength(sz...)
    xw = Array{T}(undef, (sz...,m))

    # Extract basis coefficients
    for (i, (Xwᵢ, xwᵢ)) in enumerate(zip(eachslice(Xw, dims=N), eachslice(xw, dims=N-1)))
        @inbounds treeᵢ = tree[:,i]
        @inbounds xwᵢ[:] = getbasiscoef(Xwᵢ, treeᵢ)
    end
    return xw
end

# Length of node at level L
"""
    nodelength(N, L)

Returns the node length at level L of a signal of length N when performaing wavelet packet
decomposition. Level L == 0 corresponds to the original input signal.

# Arguments
- `N::Integer`: Length of signal.
- `L::Integer`: Level of signal.

# Returns
`::Integer`: Length of nodes at level `L`.
"""
nodelength(N::Integer, L::Integer) = N >> L

# Packet table Indexing
"""
    packet(d, b, n)

Packet table indexing.

# Arguments
- `d::Integer`: Depth of splitting in packet decomposition. *Note: Depth of root node is 0.*
- `b::Integer`: Block index among 2ᵈ possibilities at depth `d`. *Note: Block indexing
  starts from 0.*
- `n::Integer`: Length of signal.

# Returns
`::UnitRange{Int64}`: Index range of block `b` of signal of length `n` at level `d`.

# Examples
```julia
using WaveletsExt

Utils.packet(0, 0, 8)       # 1:8
Utils.packet(2, 1, 8)       # 3:4
```

Translated from Wavelab850 by Zeng Fung Liew.
"""
function packet(d::Integer, b::Integer, n::Integer)
    npack = 1 << d                                  # Number of blocks in current level
    p = (b * (n÷npack) + 1) : ((b+1) * (n÷npack))   # Range of specified block
    return p
end

# Shift decomposed into BitVector
"""
    main2depthshift(sm, L)

Given the overall shift `sm`, compute the cumulative shift at each depth. Useful for
computing the shift based inverse redundant wavelet transforms.

# Arguments
- `sm::Integer`: Overall shift.
- `L::Integer`: The total number of depth.

# Returns
`::Vector{Int}`: Vector where each entry `i` describes the shift at depth `i-1`.

# Examples
```
using WaveletsExt

Utils.main2depthshift(10, 4)      # [0, 0, 2, 2, 10]
Utils.main2depthshift(5, 5)       # [0, 1, 1, 5, 5, 5]
```
"""
function main2depthshift(sm::Integer, L::Integer)
    @assert sm < 1<<L
    sb = Array{Bool,1}(undef, L)
    digits!(sb, sm, base=2)         # Compute if a shift is necessary at each depth
    d = collect(0:(L-1))            # Collect all depths
    sd = sb .<< d |> cumsum         # Compute overall shift at each depth
    pushfirst!(sd, 0)               # At depth 0, there will be no shift
    return sd
end

# Index range of coarsest scaling coefficients
"""
    coarsestscalingrange(x, tree[, redundant])
    coarsestscalingrange(n, tree[, redundant])

Given a binary tree, returns the index range of the coarsest scaling coefficients.

# Arguments
- `x::AbstractArray{T} where T<:Number`: Decomposed 1D-signal.
- `n::Integer`: Length of signal of interest.
- `tree::BitVector`: Binary tree.
- `redundant::Bool`: (Default: `false`) Whether the wavelet decomposition is redundant.
  Examples of redundant wavelet transforms are the Autocorrelation wavelet transform (ACWT),
  Stationary wavelet transform (SWT), and the Maximal Overlap wavelet transform (MOWT).

# Returns
`UnitRange{Integer}` or `::Tuple{UnitRange{Integer}, Integer}`: The index range of the
coarsest scaling subspace based on the input binary tree.

# Examples
```@repl
using Wavelets, WaveletsExt

x = randn(8)
wt = wavelet(WT.haar)
tree = maketree(x)

# Non-redundant wavelet transform
xw = wpd(x, wt)
coarsestscalingrange(xw, tree)          # 1:1

# Redundant wavelet transform
xw = swpd(x, wt)
coarsestscalingrange(xw, tree, true)    # (1:8, 8)
```

*See also:* [`finestdetailrange`](@ref)
"""
function coarsestscalingrange(x::AbstractArray{T},
                              tree::BitVector,
                              redundant::Bool=false) where T<:Number
    return coarsestscalingrange(size(x,1), tree, redundant)
end

function coarsestscalingrange(n::Integer, tree::BitVector, redundant::Bool=false)
    L = getdepth(length(tree),:binary)
    @assert L+1 == maxtransformlevels(n)
    if !redundant          # regular wt
        i = 1
        j = 0
        while i<length(tree) && tree[i]       # has children
            i = getchildindex(i,:left)
            j += 1
        end
        rng = 1:(n>>j)
    else                   # redundant wt
        i = 1
        while i<length(tree) && tree[i]       # has children
            i = getchildindex(i,:left)
        end
        rng = (1:n, i)
    end
    return rng
end

# Index range of finest detail coefficients
"""
    finestdetailrange(x, tree[, redundant])
    finestdetailrange(n, tree[, redundant])

Given a binary tree, returns the index range of the finest detail coefficients.

# Arguments
- `x::AbstractArray{T} where T<:Number`: Decomposed 1D-signal.
- `n::Integer`: Length of signal of interest.
- `tree::BitVector`: Binary tree.
- `redundant::Bool`: (Default: `false`) Whether the wavelet decomposition is redundant.
  Examples of redundant wavelet transforms are the Autocorrelation wavelet transform (ACWT),
  Stationary wavelet transform (SWT), and the Maximal Overlap wavelet transform (MOWT).

# Returns
`UnitRange{Integer}` or `::Tuple{UnitRange{Integer}, Integer}`: The index range of the
finest detail subspace based on the input binary tree.

# Examples
```@repl
using Wavelets, WaveletsExt

x = randn(8)
wt = wavelet(WT.haar)
tree = maketree(x)

# Non-redundant wavelet transform
xw = wpd(x, wt)
finestdetailrange(xw, tree)          # 8:8

# Redundant wavelet transform
xw = swpd(x, wt)
finestdetailrange(xw, tree, true)    # (1:8, 15)
```

*See also:* [`coarsestscalingrange`](@ref)
"""
function finestdetailrange(x::AbstractArray{T}, tree::BitVector,
        redundant::Bool=false) where T<:Number

    return finestdetailrange(size(x,1), tree, redundant)
end

function finestdetailrange(n::Integer, tree::BitVector, redundant::Bool=false)
    L = getdepth(length(tree),:binary)
    @assert L+1 == maxtransformlevels(n)
    if !redundant      # regular wt
        i = 1
        j = 0
        while i≤length(tree) && tree[i]
            i = getchildindex(i,:right)
            j += 1
        end
        n₀ = nodelength(n, j)
        rng = (n-n₀+1):n
    else               # redundant wt
        i = 1
        while i≤length(tree) && tree[i]
            i = getchildindex(i,:right)
        end
        rng = (1:n, i)
    end
    return rng
end

# TODO: Build a generalizable function `getslicerange` that does the same thing as
# TODO: `getrowrange` and `getcolrange` but condensed in to 1 function and generalizes into
# TODO: any arbitrary number of dimensions.
# Get range of particular row from a given signal length and quadtree index
"""
    getrowrange(n, idx)

Get the row range from a matrix with `n` rows that corresponds to `idx` from a quadtree.

# Arguments
- `n::Integer`: Number of rows in matrix.
- `idx::T where T<:Integer`: Index from a quadtree corresponding to the matrix.

# Returns
`::UnitRange{Int64}`: Row range in matrix that corresponds to `idx`.

# Examples
```@repl
using WaveletsExt

x = randn(8,8)
tree = maketree(x, 3, :full)
getrowrange(8,3)            # 1:4
```

**See also:** [`Wavelets.Util.maketree`](@ref), [`getcolrange`](@ref)
"""
function getrowrange(n::T, idx::T) where T<:Integer
    # Sanity check
    L₀ = maxtransformlevels(n)
    k = (1<<(2*L₀+2)-1)÷3       # (= sum(4^(0:L₀)))
    @assert 0 < idx ≤ k

    # Root node => range is of entire array
    if idx == 1
        return 1:n
    # Slice from parent node
    else
        # Get parent node's row range and midpoint of the range.
        parent_idx = floor(T, (idx+2)/4)
        parent_rng = getrowrange(n, parent_idx)
        midpoint = (parent_rng[begin]+parent_rng[end]) ÷ 2
        # Children nodes of subspaces LL and HL have indices 4i-2 and 4i-1 for any parent of
        # index i, these nodes take the upper half of the parent's range.
        if idx < 4*parent_idx
            return parent_rng[begin]:midpoint
        # Children nodes of subspaces LH and HH have indices 4i and 4i+1 for any parent of
        # index i, these nodes take the lower half of the parent's range.
        else
            return (midpoint+1):parent_rng[end]
        end
    end
end

# Get range of particular column from a given signal length and quadtree index
"""
    getcolrange(n, idx)

Get the column range from a matrix with `n` columns that corresponds to `idx` from a
quadtree.

# Arguments
- `n::Integer`: Number of columns in matrix.
- `idx::T where T<:Integer`: Index from a quadtree corresponding to the matrix.

# Returns
`::UnitRange{Int64}`: Column range in matrix that corresponds to `idx`.

# Examples
```@repl
using WaveletsExt

x = randn(8,8)
tree = maketree(x, 3, :full)
getcolrange(8,3)            # 5:8
```

**See also:** [`Wavelets.Util.maketree`](@ref), [`getrowrange`](@ref)
"""
function getcolrange(n::Integer, idx::T) where T<:Integer
    # Sanity check
    L₀ = maxtransformlevels(n)
    k = (1<<(2*L₀+2)-1)÷3       # (= sum(4^(0:L₀)))
    @assert 0 < idx ≤ k

    # Root node => range is of entire array
    if idx == 1
        return 1:n
    # Slice from parent node
    else
        # Get parent node's column range and midpoint of the range.
        parent_idx = floor(T, (idx+2)/4)
        parent_rng = getcolrange(n, parent_idx)
        midpoint = (parent_rng[begin]+parent_rng[end]) ÷ 2
        # Children nodes of subspaces LL and LH have indices 4i-2 and 4i for any parent of
        # index i, these nodes take the left half of the parent's range.
        if iseven(idx)
            return parent_rng[begin]:midpoint
        # Children nodes of subspaces HL and HH have indices 4i-1 and 4i+1 for any parent of
        # index i, these nodes take the right half of the parent's range.
        else
            return (midpoint+1):parent_rng[end]
        end
    end
end

include("utils_tree.jl")
include("utils_metrics.jl")
include("utils_dataset.jl")

end # end module


text:
"""
    test_grad_updates(experiment; num_games=5000)

Run a gradient update phase using dummy data.

This is useful to ensure that the chosen hyperparameters do not
lead to _Out of Memory_ errors during training.
"""
function test_grad_updates(exp::Experiment; num_games=5000)
  session = Session(exp, autosave=false, dir="sessions/grad-updates-$(exp.name)")
  UI.Log.section(session.logger, 1, "Generating $num_games some random traces")
  for i in 1:num_games
    trace = play_game(exp.gspec, RandomPlayer())
    AlphaZero.push_trace!(session.env.memory, trace, 1.0)
  end
  UI.Log.section(session.logger, 1, "Starting gradient updates")
  AlphaZero.learning_step!(session.env, session)
end


text:
function expand(x::Tuple, n::Int)
    ntuple(i -> i > length(x) ? deepcopy(x[end]) : x[i], Val(n))
end

function expand(x::Tx, n::Int) where {Tx<:AbstractVector}
    Tx[i > length(x) ? deepcopy(x[end]) : x[i] for i in 1:n]
end

function expand(x, n::Int)
    [deepcopy(x) for _ in 1:n]
end


text:
const timestamp_raw, buses_raw = readlines("day13.txt")

const timestamp = parse(Int, timestamp_raw)

const buses_p1 = parse.(Int, [b for b in split(buses_raw, ',') if b != "x"])

function part1(buses::Array{Int}, timestamp::Int)
    remainders = mod.(-timestamp, buses)
    bus_min = argmin(remainders)
    return buses[bus_min] * remainders[bus_min]
end

println("Part 1: $(part1(buses_p1, timestamp))")

buses_p2 = [[parse(Int, b), -i + 1] for (i, b) in enumerate(split(buses_raw, ",")) if b != "x"]

function chineseremainder(Ns::Vector{BigInt}, remainders::Vector{BigInt})
    n = BigInt(prod(Ns))
    Ni = n .÷ Ns
    invs = invmod.(Ni, Ns) .* Ni
    s = sum(remainders .* invs) % n
    return s >= 0 ? s : s + n
end

function part2(buses)
    Ns = @. BigInt(first(buses))
    remainders = @. BigInt(last(buses))
    return chineseremainder(Ns, remainders)
end

println("Part 2: $(part2(buses_p2))")


text:
using RiskAdjustedLinearizations, FastGaussQuadrature, Test

@testset "Gauss-Hermite Quadrature for Expectations of Functions of Independent Normally Distributed Random Variables/Vectors" begin
    f(x) = x  # calculate the expected value
    g(x) = 1. # calculate the probability

    ϵᵢ, wᵢ = RiskAdjustedLinearizations.standard_normal_gausshermite(3)
    true_eps, true_wts = gausshermite(3)
    @test ϵᵢ == true_eps .* sqrt(2.)
    @test wᵢ == true_wts ./ sqrt(π)

    mean51 = gausshermite_expectation(f, 5., 1., 5)
    mean01 = gausshermite_expectation(f, 0., 1., 5)
    mean53 = gausshermite_expectation(f, 5., 3., 5)
    mean03 = gausshermite_expectation(f, 0., 3., 5)
    prob   = gausshermite_expectation(g, -5., .1)

    @test mean51 ≈ 5.
    @test mean53 ≈ 5.
    @test isapprox(mean01, 0., atol = 1e-14)
    @test isapprox(mean03, 0., atol = 1e-14)
    @test prob ≈ 1.

    h1(x) = x[1]
    h2(x) = x[2]

    prob1 = gausshermite_expectation(g, [.5, 5.], [1., 1.], 5)
    mean11 = gausshermite_expectation(h1, [.5, 5.], [1., 1.], 5)
    mean21 = gausshermite_expectation(h2, [.5, 5.], [1., 1.], 5)
    prob2 = gausshermite_expectation(g, [5., -1.], [1., 1.], (5, 5))
    mean12 = gausshermite_expectation(h1, [.5, 5.], [1., 1.], (5, 5))
    mean22 = gausshermite_expectation(h2, [.5, 5.], [1., 1.], (5, 5))
    prob3 = gausshermite_expectation(g, [5., -1., 2.], [1., 1., 2.], [5, 5, 3])
    mean13 = gausshermite_expectation(h1, [.5, 5., 2.], [1., 1., 1.], [5, 5, 3])
    mean23 = gausshermite_expectation(h2, [.5, 5., 2.], [1., 1., 1.], [5, 5, 3])

    @test prob1 ≈ prob2 ≈ prob3 ≈ 1.
    @test mean11 ≈ mean12 ≈ mean13 ≈ .5
    @test mean21 ≈ mean22 ≈ mean23 ≈ 5.
end

nothing


text:
# Autogenerated wrapper script for fpzip_jll for x86_64-apple-darwin
export libfpzip

JLLWrappers.@generate_wrapper_header("fpzip")
JLLWrappers.@declare_library_product(libfpzip, "@rpath/libfpzip.1.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libfpzip,
        "lib/libfpzip.1.3.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#---------------------------------------------------------
# # [FBP Overview](@id 1-fbp)
#---------------------------------------------------------

#=
This example illustrates how to perform filtered back-projection (FBP)
image reconstruction in CT
using the Julia language.
=#

#=
This entire page was generated using a single Julia file:
[1-fbp.jl](@__REPO_ROOT_URL__/ct/1-fbp.jl).
=#
#md # In any such Julia documentation,
#md # you can access the source code
#md # using the "Edit on GitHub" link in the top right.

#md # The corresponding notebook can be viewed in
#md # [nbviewer](http://nbviewer.jupyter.org/) here:
#md # [`1-fbp.ipynb`](@__NBVIEWER_ROOT_URL__/ct/1-fbp.ipynb),
#md # and opened in [binder](https://mybinder.org/) here:
#md # [`1-fbp.ipynb`](@__BINDER_ROOT_URL__/ct/1-fbp.ipynb),

# This is under construction.

#=
First we add the Julia packages that are need for these examples.
Change `false` to `true` in the following code block
if you are using any of the following packages for the first time.
=#

if false
    import Pkg
    Pkg.add([
        "ImagePhantoms"
        "Plots"
        "Unitful"
        "UnitfulRecipes"
        "LaTeXStrings"
        "MIRTjim"
        "MIRT"
        "Sinograms"
        "InteractiveUtils"
    ])
end

# Now tell this Julia session to use the following packages for this example.
# Run `Pkg.add()` in the preceding code block first, if needed.

#=
using ImagePhantoms: shepp_logan, SheppLogan, radon, phantom
using Plots; default(label="", markerstrokecolor=:auto)
using Unitful: mm
using UnitfulRecipes
using LaTeXStrings
using MIRT: diffl_map, ncg
using Sinograms: todo
=#
using MIRTjim: jim, prompt
using InteractiveUtils: versioninfo

# The following line is helpful when running this jl-file as a script;
# this way it will prompt user to hit a key after each image is displayed.

isinteractive() && jim(:prompt, true);

#
isinteractive() && prompt();

# Get the ellipse parameters for a CT-suitable version of the Shepp-Logan phantom
# and calculate (analytically) its sinogram.

#=
if false
object = shepp_logan(SheppLogan(); fovs=(FOV,FOV))
sino = radon(object).(r,ϕ')
data = data / oneunit(eltype(data)) # abandon units at this point
jim(kr, kϕ, abs.(data), title="k-space data magnitude",
    xlabel=L"k_r",
    ylabel=L"k_{\phi}",
    xticks = (-1:1) .* maximum(abs, kr),
    yticks = (0,π),
    ylims = (0,π),
    aspect_ratio = :none,
)
end
=#

# ## Reproducibility

# This page was generated with the following version of Julia:

io = IOBuffer()
versioninfo(io)
split(String(take!(io)), '\n')

# And with the following package versions

import Pkg; Pkg.status()

max_stars_repo_name: JuliaImageRecon/Examples

max_issues_repo_name: JuliaImageRecon/Examples

max_forks_repo_name: JuliaImageRecon/Examples


text:
# Import the module. Using "using" keyword to make all export functions accessible.
using MORFEInvariantManifold
using MATLAB

# Name of the mesh file. The one of this example is a COMSOL mesh format.
# Meshes associated to the blade are used to benchmark the scalability of the code.
mesh_file = "beam_13.mphtxt"

### DOMAINS INFO
# domain_list is a vector that stores vectors of integers.
# Each subvector is a domain. Each integer is a COMSOL volume.
domains_list = [
[1,2,3,4,5,6,7,8]
]

### MATERIAL
# define material properties
material_name = "polysilicon";
density = 2.33e-3;
young_modulus = 160e3;
poisson_ratio = 0.22;
# add material to the library.
# Once added is saved and can be used without redefining it
MORFE_add_material(material_name,density,young_modulus,poisson_ratio)
# assign materials to the domain:
# materials is an array of strings.
# materials[i] embeds the material associated to the domain defined in domains_list[j]
materials = [
"polysilicon"
]

### BOUNDARIES INFO
# boundaries_list is a vector that stores vectors of integers.
# Each subvector is a boundary. Each integer is a COMSOL face
boundaries_list = [[1, 6, 48, 49]]

# constrained degrees of freedom of the surfaces defined above.
# It is an array n_{surface}*dim, with n_{surface} number of boundaries
# defined above and dim dimension of the problem (3D = 3). The triplet
# identifies which direction is constrained in the model.
constrained_dof = [
[1, 1, 1]
]

# It prescribes the Dirichlet boundary conditions of the mechanical
# problem. Currently it handles only homogeneous boundary conditions
# so leave it to zero.
bc_vals = [
[0.0, 0.0, 0.0]
]

# Rayleigh damping coefficients.
# The damping matrix is defined as C = α M + β K
α = 1.000890603046891/3000
β = 0.0

# FORCING INFO
# list of excitation frequencies used to excite the system.
# A parametrisation is performed for each frequency so
# avoid using more thatn 1 or 2 frequencies
Ω_list = [1.0]
# For the moment we have modal-proportional forcing.
# k_modes is a vector with length equal to the number of
# excitation frequencies. For each frequency, a list of modes
# used to define the excitation forcing is reported.
# For instance, for the first excitation
# frequency Ω_list[i] the system is excited with a forcing vector
# defined by the combination of κ_modes[i] modes, wighted by the
# load multipliers defined in κ_list[i], i.e.:
# F = sum_{j = 1}^{n_frequencies} M κ_modes[i][j]*κ_list[i][j]
# the phase of the input force is given by κ_phase. 'c' stands for cosine,
# while 's' stands for sine.
κ_modes = [[1]]
κ_list = [[0.002]]
κ_phase = [['c']]

# PARAMETRISATION INFO
# Φₗᵢₛₜ: vector of integer. Each integer correspond to which mode is included in the reduced model.
Φₗᵢₛₜ = [1,4]
# style: parametrisation style. 'g' graph, 'r' real normal form, 'c' complex normal form
style = 'g'
# max_order_a: maximum order of the asymptotic expansion of the autonomous problem
max_order_a = 3
# max_order_na: maximum order of the asymptotic expansion of the nonautonomous problem
max_order_na = 0

@time odir, Cp, Cp_na, rdyn = MORFE_mech_nonautonomous(mesh_file,domains_list,materials,
                                       boundaries_list,constrained_dof,bc_vals,
                                       Ω_list,κ_modes,κ_list,κ_phase,
                                       α,β,
                                       Φₗᵢₛₜ,style,max_order_a,
                                       max_order_na
                                       );

# post proc
ω₀ = imag(Cp[2].f[1,1]) # extract eigenfrequency
T₀ = 2*pi /ω₀  # extract period
#
L = 1.02  # characteristic length of the vibration amplitude
N = size(Cp[2].W)[1]; # number of physical DOFs
ϕ = maximum(abs.(Cp[2].W[Int(N/2)+1:N,1])) # maximum value of the mode
#
# param_init:  [artificial damping parameter = 0,  κ multiplier = 1, initial Ω = 0.8*ω₀]
param_init = [0.0,1.0,0.95*ω₀]
cont_param = 3.0 # continuation parameter is Ω
#
zero_amplitude = [0.0,0.0,0.0,0.0]
harmonics_init = [1.0,0.0]
time_integration_length = 1000*T₀
forward = true
MaxNumPoints = 450.0
minstep = 1e-8
maxstep = 0.75
initstep = 0.1
ncol = 4.0
ntst = 40.0

ms = MORFE_integrate_rdyn_frc(odir,zero_amplitude,harmonics_init,param_init,cont_param,
                                       time_integration_length,forward,MaxNumPoints,
                                       minstep,maxstep,initstep,ncol,ntst)
#
frf = MORFE_compute_frc_modal(odir,Ω_list)

x = frf[:,1]./ω₀
y = frf[:,2].*(ϕ/L)
# plot results
put_variable(ms,:x,x)
put_variable(ms,:y,y)
put_variable(ms,:max_order_a,max_order_a)

show_msession(ms) # do not close the pop up matlab windows until done with the analyses
eval_string(ms,"
figure(1);hold on
plot(x,y,'DisplayName',strcat(\"Order \",num2str(max_order_a)))
xlabel('\$\\Omega/\\omega_1\$','Interpreter','latex');
ylabel('max[\$u_1 \\phi_1\$]/\$L\$','Interpreter','latex');
xlim([0.99,1.03])
ylim([0.,1.5])
legend();
")
close(ms)