text:
# # Literated Julia Script

# This page was generated from a **julia script**, automatically converted to markdown with [Literate.jl](https://github.com/fredrikekre/Literate.jl). The conversion happens automatically when the menu is constructed.

# ## Julia Dots

# PNG with path relative to site, i.e. `/assets/images/julia-logo-dots-small.png`, but this only works in the Franklin-rendered site:

# ![Julia dots](/assets/images/julia-logo-dots-small.png)

# ## Julia Speeder

# GIF with path relative to file and within `_assets/`, i.e. `../../_assets/images/juliaspeeder32x32.gif`, and postprocessed accordingly, so it works both in Franklin and in Literate:

# ![Julia speeder](../../_assets/images/juliaspeeder32x32.gif)

# ## Math

# Inline, as in $e^{i\pi} = -1$, or in a block:
# $$
# \exp(i\pi) + 1 = 0
# $$

# ## Code chunks

# Here is a julia chunk:

x = 1

# and another:

println("hello world!")

# ## Load file

# Keep in mind that Literate renders this from the root dir, which is different from the way Weave handles it (*c.f.* Section [weavejl.jl](/pages/weaved/weavejl)):

pwd()

# With that settled, we can use the proper path (but for some weird reason, I first need to check whether the file exists before attempting to open it, otherwise Literate complains the file does not exist 🤷):

filename = "_assets/data/wow.txt"

if isfile(filename)
    open(filename) do io
        read(io, String)
    end
end

# ## Plot

using Plots

x = 0.0:0.01:2π
y = sin.(2x) + sin.(5x)

plot(x, y, title = "A nice wave", titlefont=12, label = false)

max_stars_repo_name: rmsrosa/booksjl-franklin-template

max_issues_repo_name: rmsrosa/booksjl-franklin-template

max_forks_repo_name: rmsrosa/booksjl-franklin-template


text:
"""
`prior(m::AbstractModel{T})`

Calculates log joint prior density of m.parameters.
"""
function prior{T<:AbstractFloat}(m::AbstractModel{T})
    x = zero(T)
    for θ in m.parameters
        if !θ.fixed
            x += logpdf(θ)
        end
    end
    return x
end

"""
```
posterior{T<:AbstractFloat}(m::AbstractModel{T}, data::Matrix{T};
                             mh::Bool = false, catch_errors::Bool = false)
```

Calculates and returns the log of the posterior distribution for m.parameters:
```
log posterior = log likelihood + log prior
log Pr(Θ|data)  = log Pr(data|Θ)   + log Pr(Θ)
```

### Arguments
-`m`: the model object
-`data`: matrix of data for observables

### Optional Arguments
-`mh`: Whether metropolis_hastings is the caller. If `mh=true`, the log likelihood and the
  transition matrices for the zero-lower-bound period are also returned.
-`catch_errors`: Whether or not to catch errors of type `GensysError` or `ParamBoundsError`
"""
function posterior{T<:AbstractFloat}(m::AbstractModel{T},
                                     data::Matrix{T};
                                     mh::Bool = false,
                                     catch_errors::Bool = false)
    catch_errors = catch_errors | mh
    like, out = likelihood(m, data; mh=mh, catch_errors=catch_errors)
    post = like + prior(m)
    if mh
        return Posterior(post, like, out)
    else
        return Posterior(post, like)
    end
end

"""
```
posterior!{T<:AbstractFloat}(m::AbstractModel{T}, parameters::Vector{T}, data::Matrix{T};
                              mh::Bool = false, catch_errors::Bool = false)
```

Evaluates the log posterior density at `parameters`.

### Arguments
-`m`: The model object
-`parameters`: New values for the model parameters
- `data`: Matrix of input data for observables

### Optional Arguments
-`mh`: Whether metropolis_hastings is the caller. If `mh=true`, the log likelihood and the
  transition matrices for the zero-lower-bound period are also returned.
-`catch_errors`: Whether or not to catch errors of type `GensysError` or `ParamBoundsError`.
  If `mh = true`, both should always be caught.
"""
function posterior!{T<:AbstractFloat}(m::AbstractModel{T},
                                      parameters::Vector{T},
                                      data::Matrix{T};
                                      mh::Bool = false,
                                      catch_errors::Bool = false)
    catch_errors = catch_errors | mh
    if mh
        try
            update!(m, parameters)
        catch err
            return Posterior()
        end
    else
        update!(m, parameters)
    end
    return posterior(m, data; mh=mh, catch_errors=catch_errors)

end

# Empty outputs from `likelihood`.
const LIKE_NULL_DICT   = Dict{Symbol, Matrix{AbstractFloat}}()
const LIKE_NULL_OUTPUT = (-Inf, LIKE_NULL_DICT)

"""
```
likelihood{T<:AbstractFloat}(m::AbstractModel, data::Matrix{T};
                              mh::Bool = false, catch_errors::Bool = false)
```

Evaluate the DSGE likelihood function. Can handle "two part" estimation where the observed
sample contains both a normal stretch of time (in which interest rates are positive) and
a stretch of time in which interest rates reach the zero lower bound. If there is a
zero-lower-bound period, then we filter over the 2 periods separately.  Otherwise, we
filter over the main sample all at once.

### Arguments
-`m`: The model object
-`data`: matrix of data for observables

### Optional Arguments
-`mh`: Whether metropolis_hastings is the caller. If `mh=true`, the transition matrices for
  the zero-lower-bound period are returned in a dictionary.
-`catch_errors`: If `mh = true`, `GensysErrors` should always be caught.
"""
function likelihood{T<:AbstractFloat}(m::AbstractModel,
                                      data::Matrix{T};
                                      mh::Bool = false,
                                      catch_errors::Bool = false)
    catch_errors = catch_errors | mh

    # During Metropolis-Hastings, return -∞ if any parameters are not within their bounds
    if mh
        for θ in m.parameters
            (left, right) = θ.valuebounds
            if !θ.fixed && !(left <= θ.value <= right)
                return LIKE_NULL_OUTPUT
            end
        end
    end

    # Partition sample into three regimes, and store associated matrices:
    # - R1: presample
    # - R2: normal
    # - R3: zero lower bound and beyond
    R1 = Dict{Symbol, Matrix{T}}()
    R2 = Dict{Symbol, Matrix{T}}()
    R3 = Dict{Symbol, Matrix{T}}()
    regime_mats = [R1, R2, R3]
    regime_likes = zeros(T, 3)

    n_T0         = n_presample_periods(m)
    n_ant        = n_anticipated_shocks(m)
    t_zlb_start  = zlb_start_index(m)
    n_obs_no_ant = n_observables(m) - n_anticipated_shocks(m)
    n_obs        = n_observables(m)
    n_exo        = n_shocks_exogenous(m)

    n_states_no_ant = n_states_augmented(m) - n_anticipated_shocks(m)
    n_states_aug    = n_states_augmented(m)
    n_states        = DSGE.n_states(m)
    regime_states   = [n_states_no_ant, n_states_no_ant, n_states_aug]

    R1[:data] = data[1:n_T0, 1:n_obs_no_ant]
    R2[:data] = data[(n_T0+1):t_zlb_start-1, 1:n_obs_no_ant]
    R3[:data] = data[t_zlb_start:end, :]

    # Step 1: solution to DSGE model - delivers transition equation for the state variables
    # transition equation: S_t = TC+TTT S_{t-1} +RRR eps_t, where var(eps_t) = QQ
    # If we are in MH, then any errors coming out of gensys should be caught and a -Inf
    # posterior should be returned.
    try
        R3[:TTT], R3[:RRR], R3[:CCC] = solve(m)
    catch err
        if catch_errors && isa(err, GensysError)
            info(err.msg)
            return LIKE_NULL_OUTPUT
        else
            rethrow(err)
        end
    end

    # Get normal, no ZLB matrices
    state_inds = [1:(n_states-n_ant); (n_states+1):n_states_aug]
    shock_inds = 1:(n_exo-n_ant)

    R2[:TTT] = R3[:TTT][state_inds, state_inds]
    R2[:RRR] = R3[:RRR][state_inds, shock_inds]
    R2[:CCC] = R3[:CCC][state_inds, :]

    ## step 2: define the measurement equation: X_t = ZZ*S_t + D + u_t
    ## where u_t = eta_t+MM* eps_t with var(eta_t) = EE
    ## where var(u_t) = HH = EE+MM QQ MM', cov(eps_t,u_t) = VV = QQ*MM'

    # Get measurement equation matrices set up for normal and zlb periods
    measurement_R2 = measurement(m, R2[:TTT], R2[:RRR], R2[:CCC]; shocks=false)
    measurement_R3 = measurement(m, R3[:TTT], R3[:RRR], R3[:CCC]; shocks=true)
    for d in (:ZZ, :DD, :QQ, :VVall)
        R2[d] = measurement_R2[d]
        R3[d] = measurement_R3[d]
    end

    # Presample measurement & transition equation matrices are same as normal period
    for d in (:TTT, :RRR, :QQ, :ZZ, :DD, :VVall)
        R1[d] = R2[d]
    end

    ## step 3: compute log-likelihood using Kalman filter
    ##         note that kalman_filter function assumes a transition equation written as:
    ##         S_t = TTT S_{t-1} +eps2_t, where eps2_t = RRReps_t
    ##         therefore redefine QQ2 = var(eps2_t) = RRR*QQ*RRR'
    ##         and  VV2 = cov(eps2_t,u_u) = RRR*VV
    ##         define VVall as the joint variance of the two shocks VVall = var([eps2_tu_t])

    # Run Kalman filter on presample
    R1[:A0]         = zeros(T, n_states_no_ant, 1)
    R1[:P0]         = solve_discrete_lyapunov(R1[:TTT], R1[:RRR]*R1[:QQ]*R1[:RRR]')
    out             = kalman_filter(R1[:data]', 1, zeros(T, n_states_no_ant, 1), R1[:TTT], R1[:DD], R1[:ZZ], R1[:VVall], R1[:A0], R1[:P0])
    regime_likes[1] = out[:L]
    R1[:zend]       = out[:zend]
    R1[:Pend]       = out[:Pend]

    # Run Kalman filter on normal period
    zprev           = R1[:zend]
    Pprev           = R1[:Pend]
    out             = kalman_filter(R2[:data]', 1, zeros(regime_states[2], 1), R2[:TTT], R2[:DD], R2[:ZZ], R2[:VVall], zprev, Pprev)
    regime_likes[2] = out[:L]
    R2[:zend]       = out[:zend]
    R2[:Pend]       = out[:Pend]

    # Run Kalman filter on ZLB period
    # This section expands the number of states to accomodate extra states for the
    # anticipated policy shocks. It does so by taking the zend and Pend for the
    # state space without anticipated policy shocks, then shoves in nant
    # zeros in the middle of zend and Pend in the location of
    # the anticipated shock entries.
    before_shocks    = 1:(n_states-n_ant)
    after_shocks_old = (n_states-n_ant+1):(n_states_aug-n_ant)
    after_shocks_new = (n_states+1):n_states_aug

    zprev = [R2[:zend][before_shocks, :];
             zeros(T, n_ant, 1);
             R2[:zend][after_shocks_old, :]]

    Pprev                                     = zeros(T, n_states_aug, n_states_aug)
    Pprev[before_shocks, before_shocks]       = R2[:Pend][before_shocks, before_shocks]
    Pprev[before_shocks, after_shocks_new]    = R2[:Pend][before_shocks, after_shocks_old]
    Pprev[after_shocks_new, before_shocks]    = R2[:Pend][after_shocks_old, before_shocks]
    Pprev[after_shocks_new, after_shocks_new] = R2[:Pend][after_shocks_old, after_shocks_old]

    out             = kalman_filter(R3[:data]', 1, zeros(regime_states[3], 1), R3[:TTT], R3[:DD], R3[:ZZ], R3[:VVall], zprev, Pprev)
    regime_likes[3] = out[:L]
    R3[:zend]       = out[:zend]
    R3[:Pend]       = out[:Pend]

    # Return total log-likelihood, excluding the presample
    like = regime_likes[2] + regime_likes[3]
    if mh
        return like, R3
    else
        return like, LIKE_NULL_DICT
    end
end

# Type returned by posterior
immutable Posterior{T<:AbstractFloat}
    post::T
    like::T
    mats::Dict{Symbol, Matrix{T}}
end
function Posterior{T<:AbstractFloat}(post::T = -Inf,
                                     like::T = -Inf,
                                     mats::Dict{Symbol,Matrix{T}}    = Dict{Symbol,Matrix{T}}())
    return Posterior{T}(post, like, mats)
end
function Base.getindex(P::Posterior, d::Symbol)
    if d in (:post, :like, :mats)
        return getfield(P, d)
    else
        throw(KeyError(d))
    end
end

ext: jl


text:
### part of package SupremumInfimum

import Base:  Callable, reduce_empty, mapreduce_empty, union, intersect

for (sup, supremum, inf, infimum, max, typemin, union, ubound, lbound) in
((:sup,:supremum,:inf,:infimum,:max,:typemin,:union,     "least upper bound",    "lower bound"),
 (:inf,:infimum,:sup,:supremum,:min,:typemax,:intersect, "greatest lower bound", "upper bound"))

  global supstring = string(supremum)
  global infstring = string(infimum)
  global ubound = ubound
  global lbound = lbound

@eval begin
"""
    $supstring(a::S, b::T)

Calculate $ubound of `a` and `b`. If the type `promote_type(S,T)`
has a canonical partial order (implements `<` but does not implement `isless`),
the method must be overwritten for this type. If `isless` is defined, `$max` is called.
"""
($sup)(a::T, b::T) where {T<:Any} = ($max)(a, b)
($sup)(a::T, b::T) where {T<:Real} = ($max)(a, b)
($sup)(a::S, b::T) where {S<:Real,T<:Real} = ($max)(a, b)
($sup)(a::T, b::T) where {T<:AbstractString} = ($max)(a, b)
($sup)(a::S, b::T) where {S<:AbstractSet,T<:AbstractSet} = ($union)(a, b)
($sup)(a::S, b::T) where {S<:AbstractArray,T<:AbstractArray} = ($union)(a, b)

$sup(x::Real) = x

# This schould go to operators.jl:419
    ($sup)(a, b, c, xs...) = Base.afoldl(($sup), ($sup)(($sup)(a,b), c), xs...)

"""
    $supstring(itr)

Calculate $ubound of the elements of collection `itr`.
If itr is empty, return a $lbound of the appropriate type.
"""
($supremum)(itr) = mapreduce(identity, ($sup), itr)

"""
    $supstring(f::Callable, itr)

Calculate $ubound of the elements of `f.(itr)`.
If itr is empty, return a $lbound of the return type of `f`.
"""
($supremum)(f::Callable, itr) = mapreduce(f, ($sup), itr)

######
# this is inspired by reduce.jl:243 ff
# note that we do not enforce the return type of ($infimum)
reduce_empty(::typeof(($sup)), T) = ($infimum)(T)

# derive inferred return type
function mapreduce_empty(f, op::typeof(($supremum)), T)
  S = Core.Inference.return_type(f, (T,))
  ($infimum)(ifelse(S<:Union{}, Any, S))
end

"""
    $infstring(::Type{T})

Calculate $lbound for all elements of type `T`.
"""
($infimum)(::Type{T}) where {T<:Any} = ($typemin)(leaf_type(T)) # default behaviour

end # eval begin
end # for

infimum(::Type{T}) where {T<:AbstractSet} = leaf_type(T)()
infimum(::Type{String}) = ""
infimum(::Type{BigInt}) = IntInf

supremum(::Type{T}) where {T<:AbstractSet} = error("no upper bound defined (yet) for Set")
supremum(::Type{String}) = error("no upper bound defined for String")
supremum(::Type{BigInt}) = IntSup

"""
  leaf_type(::Type) -> ::Type

Return a default concrete subtype of argument.
"""

leaf_type(::Type{T}) where {T<:AbstractSet{U}} where {U<:Any} = isconcretetype(T) ? T : Set{U}
leaf_type(::Type{T}) where {T<:AbstractArray{U,N}} where {U,N} = isconcretetype(T) ? T : Array{U,N}
leaf_type(::Type{BigInt}) = BigInt
leaf_type(::Type{T}) where {T<:Signed} = isconcretetype(T) ? T : Int
leaf_type(::Type{T}) where {T<:Unsigned} = isconcretetype(T) ? T : UInt
leaf_type(::Type{T}) where {T<:Integer} = isconcretetype(T) ? T : Int
leaf_type(::Type{T}) where {T<:Real} = isconcretetype(T) ? T : Float64


text:
struct SimulationState
    current_time::Base.RefValue{Dates.DateTime}
    last_decision_model::Base.RefValue{Symbol}
    decision_states::DatasetContainer{DataFrameDataset}
    system_states::DatasetContainer{DataFrameDataset}
end

function SimulationState()
    return SimulationState(
        Ref(UNSET_INI_TIME),
        Ref(:None),
        DatasetContainer{DataFrameDataset}(),
        DatasetContainer{DataFrameDataset}(),
    )
end

get_current_time(s::SimulationState) = s.current_time[]
get_last_decision_model(s::SimulationState) = s.last_decision_model[]
get_decision_states(s::SimulationState) = s.decision_states
get_system_states(s::SimulationState) = s.system_states

# Not to be used in hot loops
function get_system_states_resolution(s::SimulationState)
    system_state = get_system_states(s)
    # All the system states have the same resolution
    return get_data_resolution(first(values(system_state.variables)))
end

function set_current_time!(s::SimulationState, val::Dates.DateTime)
    s.current_time[] = val
    return
end

function set_last_decision_model!(s::SimulationState, val::Symbol)
    s.last_decision_model[] = val
    return
end

const STATE_TIME_PARAMS = NamedTuple{(:horizon, :resolution), NTuple{2, Dates.Millisecond}}

function _get_state_params(models::SimulationModels, simulation_step::Dates.Millisecond)
    params = OrderedDict{OptimizationContainerKey, STATE_TIME_PARAMS}()
    for model in get_decision_models(models)
        container = get_optimization_container(model)
        model_resolution = get_resolution(model)
        model_interval = get_interval(model)
        horizon_length = get_horizon(model) * model_resolution
        # This is the portion of the Horizon that "overflows" into the next step
        time_residual = horizon_length - model_interval
        @assert_op time_residual >= zero(Dates.Millisecond)
        num_runs = simulation_step / model_interval
        total_time = (num_runs - 1) * model_interval + horizon_length
        for type in fieldnames(DatasetContainer)
            field_containers = getfield(container, type)
            for key in keys(field_containers)
                !should_write_resulting_value(key) && continue
                if !haskey(params, key)
                    params[key] = (
                        horizon=max(simulation_step + time_residual, total_time),
                        resolution=model_resolution,
                    )
                else
                    params[key] = (
                        horizon=max(params[key].horizon, total_time),
                        resolution=min(params[key].resolution, model_resolution),
                    )
                end
                @debug get_name(model) key params[key]
            end
        end
    end
    return params
end

function _initialize_model_states!(
    sim_state::SimulationState,
    model::OperationModel,
    simulation_initial_time::Dates.DateTime,
    simulation_step::Dates.Millisecond,
    params::OrderedDict{OptimizationContainerKey, STATE_TIME_PARAMS},
)
    states = get_decision_states(sim_state)
    container = get_optimization_container(model)
    for field in fieldnames(DatasetContainer)
        field_containers = getfield(container, field)
        field_states = getfield(states, field)
        for (key, value) in field_containers
            !should_write_resulting_value(key) && continue
            value_counts = params[key].horizon ÷ params[key].resolution
            column_names = get_column_names(key, value)
            if !haskey(field_states, key) || length(field_states[key]) < value_counts
                field_states[key] = DataFrameDataset(
                    DataFrames.DataFrame(
                        fill(NaN, value_counts, length(column_names)),
                        column_names,
                    ),
                    collect(
                        range(
                            simulation_initial_time,
                            step=params[key].resolution,
                            length=value_counts,
                        ),
                    ),
                    params[key].resolution,
                    Int(simulation_step / params[key].resolution),
                )
            end
        end
    end
    return
end

function _initialize_system_states!(
    sim_state::SimulationState,
    ::Nothing,
    simulation_initial_time::Dates.DateTime,
    params::OrderedDict{OptimizationContainerKey, STATE_TIME_PARAMS},
)
    decision_states = get_decision_states(sim_state)
    emulator_states = get_system_states(sim_state)
    min_res = minimum([v.resolution for v in values(params)])
    for key in get_dataset_keys(decision_states)
        cols = get_column_names(key, get_dataset(decision_states, key))
        set_dataset!(
            emulator_states,
            key,
            make_system_state(
                DataFrames.DataFrame(cols .=> NaN),
                simulation_initial_time,
                min_res,
            ),
        )
    end
    return
end

function _initialize_system_states!(
    sim_state::SimulationState,
    emulation_model::EmulationModel,
    simulation_initial_time::Dates.DateTime,
    params::OrderedDict{OptimizationContainerKey, STATE_TIME_PARAMS},
)
    decision_states = get_decision_states(sim_state)
    emulator_states = get_system_states(sim_state)
    emulation_container = get_optimization_container(emulation_model)

    for field in fieldnames(DatasetContainer)
        field_containers = getfield(emulation_container, field)
        for (key, value) in field_containers
            !should_write_resulting_value(key) && continue
            column_names = get_column_names(key, value)
            set_dataset!(
                emulator_states,
                key,
                make_system_state(
                    DataFrames.DataFrame(column_names .=> NaN),
                    simulation_initial_time,
                    get_resolution(emulation_model),
                ),
            )
        end
    end

    for key in get_dataset_keys(decision_states)
        if has_dataset(emulator_states, key)
            dm_cols = get_column_names(key, get_dataset(decision_states, key))
            em_cols = get_column_names(key, get_dataset(emulator_states, key))
            @assert_op dm_cols == em_cols
            continue
        end
        cols = get_column_names(key, get_dataset(decision_states, key))
        set_dataset!(
            emulator_states,
            key,
            make_system_state(
                DataFrames.DataFrame(cols .=> NaN),
                simulation_initial_time,
                get_resolution(emulation_model),
            ),
        )
    end
    return
end

function initialize_simulation_state!(
    sim_state::SimulationState,
    models::SimulationModels,
    simulation_step::Dates.Millisecond,
    simulation_initial_time::Dates.DateTime,
)
    params = _get_state_params(models, simulation_step)
    for model in get_decision_models(models)
        _initialize_model_states!(
            sim_state,
            model,
            simulation_initial_time,
            simulation_step,
            params,
        )
    end
    set_last_decision_model!(sim_state, get_name(last(get_decision_models(models))))
    em = get_emulation_model(models)
    _initialize_system_states!(sim_state, em, simulation_initial_time, params)
    return
end

function update_decision_state!(
    state::SimulationState,
    key::OptimizationContainerKey,
    store_data::DataFrames.DataFrame,
    simulation_time::Dates.DateTime,
    model_params::ModelStoreParams,
)
    state_data = get_decision_state_data(state, key)
    model_resolution = get_resolution(model_params)
    state_resolution = get_data_resolution(state_data)
    resolution_ratio = model_resolution ÷ state_resolution
    state_timestamps = state_data.timestamps
    @assert_op resolution_ratio >= 1

    if simulation_time > get_end_of_step_timestamp(state_data)
        state_data_index = 1
        state_data.timestamps[:] .=
            range(simulation_time, step=state_resolution, length=length(state_data))
    else
        state_data_index = find_timestamp_index(state_timestamps, simulation_time)
    end

    offset = resolution_ratio - 1
    result_time_index = axes(store_data)[1]
    set_update_timestamp!(state_data, simulation_time)
    for t in result_time_index
        state_range = state_data_index:(state_data_index + offset)
        for name in DataFrames.names(store_data), i in state_range
            # TODO: We could also interpolate here
            state_data.values[i, name] = store_data[t, name]
        end
        set_last_recorded_row!(state_data, state_range[end])
        state_data_index += resolution_ratio
    end
    return
end

function update_decision_state!(
    state::SimulationState,
    key::AuxVarKey{EnergyOutput, T},
    store_data::DataFrames.DataFrame,
    simulation_time::Dates.DateTime,
    model_params::ModelStoreParams,
) where {T <: PSY.Component}
    state_data = get_decision_state_data(state, key)
    model_resolution = get_resolution(model_params)
    state_resolution = get_data_resolution(state_data)
    resolution_ratio = model_resolution ÷ state_resolution
    state_timestamps = state_data.timestamps
    @assert_op resolution_ratio >= 1

    if simulation_time > get_end_of_step_timestamp(state_data)
        state_data_index = 1
        state_data.timestamps[:] .=
            range(simulation_time, step=state_resolution, length=length(state_data))
    else
        state_data_index = find_timestamp_index(state_timestamps, simulation_time)
    end

    offset = resolution_ratio - 1
    result_time_index = axes(store_data)[1]
    set_update_timestamp!(state_data, simulation_time)
    for t in result_time_index
        state_range = state_data_index:(state_data_index + offset)
        for name in DataFrames.names(store_data), i in state_range
            state_data.values[i, name] = store_data[t, name] / resolution_ratio
        end
        set_last_recorded_row!(state_data, state_range[end])
        state_data_index += resolution_ratio
    end

    return
end

function update_decision_state!(
    state::SimulationState,
    key::AuxVarKey{S, T},
    store_data::DataFrames.DataFrame,
    simulation_time::Dates.DateTime,
    model_params::ModelStoreParams,
) where {T <: PSY.Component, S <: Union{TimeDurationOff, TimeDurationOn}}
    state_data = get_decision_state_data(state, key)
    model_resolution = get_resolution(model_params)
    state_resolution = get_data_resolution(state_data)
    resolution_ratio = model_resolution ÷ state_resolution
    @assert_op resolution_ratio >= 1

    if simulation_time > get_end_of_step_timestamp(state_data)
        state_data_index = 1
        state_data.timestamps[:] .=
            range(simulation_time, step=state_resolution, length=length(state_data))
    else
        state_data_index = find_timestamp_index(state_data.timestamps, simulation_time)
    end

    offset = resolution_ratio - 1
    result_time_index = axes(store_data)[1]
    set_update_timestamp!(state_data, simulation_time)

    if resolution_ratio == 1.0
        increment_per_period = 1.0
    elseif state_resolution < Dates.Hour(1) && state_resolution > Dates.Minute(1)
        increment_per_period = Dates.value(Dates.Minute(state_resolution))
    else
        error("Incorrect Problem Resolution specification")
    end

    for t in result_time_index
        state_range = state_data_index:(state_data_index + offset)
        @assert_op state_range[end] <= length(state_data)
        for name in DataFrames.names(store_data), i in state_range
            if t == 1 && i == 1
                state_data.values[i, name] = store_data[t, name] * resolution_ratio
            else
                state_data.values[i, name] =
                    store_data[t, name] > 0 ?
                    state_data.values[i - 1, name] + increment_per_period : 0
            end
        end
        set_last_recorded_row!(state_data, state_range[end])
        state_data_index += resolution_ratio
    end

    return
end

function get_decision_state_data(state::SimulationState, key::OptimizationContainerKey)
    return get_dataset(get_decision_states(state), key)
end

function get_decision_state_value(state::SimulationState, key::OptimizationContainerKey)
    return get_dataset_values(get_decision_states(state), key)
end

function get_decision_state_value(
    state::SimulationState,
    key::OptimizationContainerKey,
    date::Dates.DateTime,
)
    return get_dataset_values(get_decision_states(state), key, date)
end

function get_system_state_data(state::SimulationState, key::OptimizationContainerKey)
    return get_dataset(get_system_states(state), key)
end

function get_system_state_value(state::SimulationState, key::OptimizationContainerKey)
    return get_dataset_values(get_system_states(state), key)[1, :]
end

function update_system_state!(
    state::DatasetContainer{DataFrameDataset},
    key::OptimizationContainerKey,
    store::SimulationStore,
    model_name::Symbol,
    simulation_time::Dates.DateTime,
)
    em_data = get_em_data(store)
    ix = get_last_recorded_row(em_data, key)
    res = read_result(DataFrames.DataFrame, store, model_name, key, ix)
    dataset = get_dataset(state, key)
    set_update_timestamp!(dataset, simulation_time)
    set_dataset_values!(state, key, 1, res)
    set_last_recorded_row!(dataset, 1)
    return
end

function update_system_state!(
    state::DatasetContainer{DataFrameDataset},
    key::OptimizationContainerKey,
    decision_state::DatasetContainer{DataFrameDataset},
    simulation_time::Dates.DateTime,
)
    decision_dataset = get_dataset(decision_state, key)
    # Gets the timestamp of the value used for the update, which might not match exactly the
    # simulation time since the value might have not been updated yet
    ts = get_value_timestamp(decision_dataset, simulation_time)
    system_dataset = get_dataset(state, key)
    get_update_timestamp(system_dataset)
    if ts == get_update_timestamp(system_dataset)
        # Uncomment for debugging
        #@warn "Skipped overwriting data with the same timestamp \\
        #       key: $(encode_key_as_string(key)), $(simulation_time), $ts"
        return
    end

    if get_update_timestamp(system_dataset) > ts
        error("Trying to update with past data a future state timestamp \\
            key: $(encode_key_as_string(key)), $(simulation_time), $ts")
    end

    # Writes the timestamp of the value used for the update
    set_update_timestamp!(system_dataset, ts)
    # Keep coordination between fields. System state is an array of size 1
    system_dataset.timestamps[1] = ts
    set_dataset_values!(state, key, 1, get_dataset_value(decision_dataset, simulation_time))
    # This value shouldn't be other than one and after one execution is no-op.
    set_last_recorded_row!(system_dataset, 1)
    return
end

function update_system_state!(
    state::DatasetContainer{DataFrameDataset},
    key::AuxVarKey{T, PSY.ThermalStandard},
    decision_state::DatasetContainer{DataFrameDataset},
    simulation_time::Dates.DateTime,
) where {T <: Union{TimeDurationOn, TimeDurationOff}}
    decision_dataset = get_dataset(decision_state, key)
    # Gets the timestamp of the value used for the update, which might not match exactly the
    # simulation time since the value might have not been updated yet

    ts = get_value_timestamp(decision_dataset, simulation_time)
    system_dataset = get_dataset(state, key)
    system_state_resolution = get_data_resolution(system_dataset)
    decision_state_resolution = get_data_resolution(decision_dataset)

    decision_state_value = get_dataset_value(decision_dataset, simulation_time)

    if ts == get_update_timestamp(system_dataset)
        # Uncomment for debugging
        #@warn "Skipped overwriting data with the same timestamp \\
        #       key: $(encode_key_as_string(key)), $(simulation_time), $ts"
        return
    end

    if get_update_timestamp(system_dataset) > ts
        error("Trying to update with past data a future state timestamp \\
            key: $(encode_key_as_string(key)), $(simulation_time), $ts")
    end

    # Writes the timestamp of the value used for the update
    set_update_timestamp!(system_dataset, ts)
    # Keep coordination between fields. System state is an array of size 1
    system_dataset.timestamps[1] = ts
    time_ratio = (decision_state_resolution / system_state_resolution)
    # Don't use set_dataset_values!(state, key, 1, decision_state_value).
    # For the time variables we need to grab the values to avoid mutation of the
    # dataframe row
    set_value!(system_dataset, values(decision_state_value) .* time_ratio, 1)
    # This value shouldn't be other than one and after one execution is no-op.
    set_last_recorded_row!(system_dataset, 1)
    return
end

function get_system_state_value(
    state::SimulationState,
    ::T,
    ::Type{U},
) where {T <: VariableType, U <: Union{PSY.Component, PSY.System}}
    return get_system_state_value(state, VariableKey(T, U))
end

function get_system_state_value(
    state::SimulationState,
    ::T,
    ::Type{U},
) where {T <: AuxVariableType, U <: Union{PSY.Component, PSY.System}}
    return get_system_state_value(state, AuxVarKey(T, U))
end

function get_system_state_value(
    state::SimulationState,
    ::T,
    ::Type{U},
) where {T <: ConstraintType, U <: Union{PSY.Component, PSY.System}}
    return get_system_state_value(state, ConstraintKey(T, U))
end

function get_system_state_data(
    state::SimulationState,
    ::T,
    ::Type{U},
) where {T <: VariableType, U <: Union{PSY.Component, PSY.System}}
    return get_system_state_data(state, VariableKey(T, U))
end

function get_system_state_data(
    state::SimulationState,
    ::T,
    ::Type{U},
) where {T <: AuxVariableType, U <: Union{PSY.Component, PSY.System}}
    return get_system_state_data(state, AuxVarKey(T, U))
end

function get_system_state_data(
    state::SimulationState,
    ::T,
    ::Type{U},
) where {T <: ConstraintType, U <: Union{PSY.Component, PSY.System}}
    return get_system_state_data(state, ConstraintKey(T, U))
end


text:
"""
A module which provides visualizations
for `Trajectory` and `Manifold`
instances, as well as other
common astrodynamical visualizations.

# Extended Help

**Exports**

$(EXPORTS)

**Imports**

$(IMPORTS)
"""
module Visualizations

export zerovelocityplot, zerovelocityplot!

using Plots, RecipesBase
using DifferentialEquations
using ..Calculations
using ..CoordinateFrames
using ..States, ..Propagation
using DocStringExtensions

@template (FUNCTIONS, METHODS, MACROS) =
    """
    $(SIGNATURES)
    $(DOCSTRING)
    """

@template (TYPES, CONSTANTS) =
    """
    $(TYPEDEF)
    $(DOCSTRING)
    """

"""
Rather than supply indices, users can use `Symbol`
instances to select indices for plotting.

# Example

traj = propagate(orbit, period)
plot(traj; vars=:tx)
plot(traj; vars=xyz)
plot(traj; vars=xẋ)
"""
function process_vars(vars::Symbol)
    @assert 2 ≤ length(string(vars)) ≤ 3 "Keyword argument `vars` has an inproper length $(length(vars)). Choose a length in [2, 3]."
    @assert sum(c -> count(c, lowercase(string(vars))), ("t", "x", "y", "z", "ẋ", "ẏ", "ż")) == length(string(vars)) "Invalid character provided in `vars`."

    indices = Dict(
        "t" => 0,
        "x" => 1,
        "y" => 2,
        "z" => 3,
        "ẋ" => 4,
        "ẏ" => 5,
        "ż" => 6
    )

    return [indices[string(char)] for char ∈ lowercase(string(vars))]
end

include(joinpath("Trajectories", "Trajectories.jl"))
include(joinpath("Manifolds", "Manifolds.jl"))
include(joinpath("Energy", "Energy.jl"))

end # module


text:
using Lasso
using Random, SparseArrays
DATADIR = joinpath(dirname(@__FILE__), "data")

function diffmatslow(order, n)
    D = sparse(SparseArrays.spdiagm_internal(0 => fill(-1., n-1), 1=>fill(1., n-1))...)
    for i = 1:order
        D = sparse(SparseArrays.spdiagm_internal(0 => fill(-1., n-i-1), 1=>fill(1., n-i-1))...)*D
    end
    D
end

# Test that DifferenceMatrix behaves as expected
@testset "DifferenceMatrix" begin
    @testset "order = $(order)" for order in (1, 2, 3)
        D1 = diffmatslow(order, 100)
        D2 = Lasso.TrendFiltering.DifferenceMatrix{Float64}(order, 100)
        x = randn(testrng, 100)
        @test D2'D2 == D1'D1
        @test D2*x ≈ D1*x
        @test D2'*x[1:size(D1, 1)] ≈ D1'*x[1:size(D1, 1)]
    end
end

# Test against results from glmgen
lakehuron = readcsvmat(joinpath(DATADIR, "LakeHuron.csv"); header=true)[:, 3]
@testset "TrendFiltering" begin
    @testset "order = $(order)" for order in (1, 2, 3)
        @testset "lambda = $(lambda)" for lambda in (1., 10., 100.)
            @test coef(fit(TrendFilter, lakehuron, order, lambda; tol=1e-9)) ≈
                  vec(readcsvmat(joinpath(DATADIR, "LakeHuron_order_$(order)_lambda_$(lambda).csv")))
        end
    end
end


text:
if basename(pwd()) == "aoc"
    cd("2017/2")
end

function part1(filename::AbstractString)
    sum(-(reverse(extrema(parse.(Int, split(line))))...) for line in eachline(filename))
end

@assert part1("example.txt") == 18
part1("input.txt")

findround(str::AbstractString) = findround(parse.(Int, split(str)))

function findround(xs::Vector{Int})
    for i in firstindex(xs):lastindex(xs)-1
        for j in i+1:lastindex(xs)
            lo, hi = minmax(xs[i], xs[j])
            if isinteger(hi / lo)
                return hi ÷ lo
            end
        end
    end
end

part2(filename::AbstractString) = sum(findround.(eachline(filename)))

@assert part2("example2.txt") == 9
part2("input.txt")

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
using AdventOfCode2018
using AdventOfCode2018.Day16

inp = get_aoc_input(16)

examples, test_program = parse_input(inp)

println("Solution part1: ", part_1(examples))
println("Solution part2: ", part_2(examples, test_program))


text:
export  HyperHtml,
    html_, head_, title_, meta_,
    # sections https://html.spec.whatwg.org/#sections
    body_, article_, section_, nav_, aside_,
    h1_, h2_, h3_, h4_, h5_, h6_,
    hgroup_, header_, footer_, address_,
    # grouping https://html.spec.whatwg.org/#grouping-content
    p_, hr_, pre_, blockquote_, ol_, ul_,
    menu_, li_, dl_, dt_, dd_, figure_,
    figcaption_, main_, div_,
    # table https://html.spec.whatwg.org/#tables
    table_, caption_, colgroup_, col_, tbody_,
    thead_, tfoot_, tr_, td_,  th_,
    # text level https://html.spec.whatwg.org/#text-level-semantics
    a_, em_, strong_, small_, s_, cite_,
    q_, dfn_, abbr_, ruby_, rt_, rp_,
    data_, code_, var_, samp_, kbd_,
    sub_, sup_, i_, b_, u_, mark_,
    bdi_, bdo_, span_, br_, wbr_,
    # input element types
    button_, checkbox_, color_, date_,
    datetimeLocal_, email_, file_, hidden_,
    image_, month_, number_, password_,
    radio_, range_, reset_, search_, submit_,
    tel_, text_, time_, url_, week_,
    # other
    area_

# hyperhtml = Dict{Symbol, Hyperscript.Node{Hyperscript.HTMLSVG}}([])
HyperHtmlDict = Dict{Symbol, Function}([])

for (T,F) in (("html", :html_), ("head", :head_), ("title", :title_), ("meta", :meta_),
     # sections https://html.spec.whatwg.org/#sections
    ("body", :body_),  ("article", :article_), ("section", :section_), ("nav", :nav_), ("aside", :aside_),
    ("h1", :h1_), ("h2", :h2_), ("h3", :h3_), ("h4", :h4_), ("h5", :h5_), ("h6", :h6_),
    ("hgroup", :hgroup_), ("header", :header_), ("footer", :footer_), ("address", :address_),
    # grouping https://html.spec.whatwg.org/#grouping-content
    ("p", :p_), ("hr", :hr_), ("pre", :pre_), ("blockquote", :blockquote_), ("ol", :ol_), ("ul", :ul_),
    ("menu", :menu_), ("li", :li_), ("dl", :dl_), ("dt", :dt_), ("dd", :dd_), ("figure", :figure_),
    ("figcaption", :figcaption_), ("main", :main_), ("div", :div_),
    # table https://html.spec.whatwg.org/#tables
    ("table", :table_), ("caption", :caption_), ("colgroup", :colgroup_), ("col", :col_), ("tbody", :tbody_),
    ("thead", :thead_), ("tfoot", :tfoot_), ("tr", :tr_), ("td", :td_),  ("th", :th_),
    # text level https://html.spec.whatwg.org/#text-level-semantics
    ("a", :a_), ("em", :em_), ("strong", :strong_), ("small", :small_), ("s", :s_), ("cite", :cite_),
    ("q", :q_), ("dfn", :dfn_), ("abbr", :abbr_), ("ruby", :ruby_), ("rt", :rt_), ("rp", :rp_),
    ("data", :data_), ("code", :code_), ("var", :var_), ("samp", :samp_), ("kbd", :kbd_),
    ("sub", :sub_), ("sup", :sup_), ("i", :i_), ("b", :b_), ("u", :u_), ("mark", :mark_),
    ("bdi", :bdi_), ("bdo", :bdo_), ("span", :span_), ("br", :br_), ("wbr", :wbr_),
    # input element types
    ("button", :button_), ("checkbox", :checkbox_), ("color", :color_), ("date", :date_),
    ("datetimeLocal", :datetimeLocal_), ("email", :email_), ("file", :file_), ("hidden", :hidden_),
    ("image", :image_), ("month", :month_), ("number", :number_), ("password", :password_),
    ("radio", :radio_), ("range", :range_), ("reset", :reset_), ("search", :search_), ("submit", :submit_),
    ("tel", :tel_), ("text", :text_), ("time", :time_), ("url", :url_), ("week", :week_),
    # other
    ("area", :area_),
    )
  @eval begin
     $F() = m($T)
     HyperHtmlDict[Symbol($T)] = $F
     $F(class::String) = m($T, class=class)
     function $F(class::String, xs::Vararg{Pair,N}) where {N}
         dict = Dict(xs)
         dict[:class] = class
         return Hyperscript.Node(Hyperscript.DEFAULT_HTMLSVG_CONTEXT, $T, [], dict)
     end
     function $F(xs::Vararg{Pair,N}) where {N}
         dict = Dict(xs)
         return Hyperscript.Node(Hyperscript.DEFAULT_HTMLSVG_CONTEXT, $T, [], dict)
     end
  end
end

const HyperHtml = copy(HyperHtmlDict)
HyperHtmlDict = nothing

# link types
const a_alternate_ = a_(:rel=>"alternate")
const a_author_ = a_(:rel=>"author")
const a_bookmark_ = a_(:rel=>"bookmark")
const a_canonical_ = a_(:rel=>"canonical")
const a_dnsprefetch_ = a_(:rel=>"dns-prefetch")
const a_external_ = a_(:rel=>"external")
const a_help_ = a_(:rel=>"help")
const a_icon_ = a_(:rel=>"icon")
const a_license_ = a_(:rel=>"license")
const a_modulepreupload_ = a_(:rel=>"modulepreupload")
const a_nofollow_ = a_(:rel=>"nofollow")
const a_noopener_ = a_(:rel=>"noopener")
const a_noreferrer_ = a_(:rel=>"noreferrer")
const a_opener_ = a_(:rel=>"opener")
const a_pingback_ = a_(:rel=>"pingback")
const a_preconnect_ = a_(:rel=>"preconnect")
const a_prefetch_ = a_(:rel=>"prefetch")
const a_preload_ = a_(:rel=>"preload")
const a_prerender_ = a_(:rel=>"prerender")
const a_search_ = a_(:rel=>"search")
const a_stylesheet_ = a_(:rel=>"stylesheet")
const a_tag_ = a_(:rel=>"tag")
const a_next_ = a_(:rel=>"next")
const a_prev_ = a_(:rel=>"prev")

# area https://html.spec.whatwg.org/#the-area-element
const area_circle_ = area_(:shape=>"circle")
const area_rectangle_ = area_(:shape=>"rect")
const area_polygon_ = area_(:shape=>"poly")

### html elements

# html elements
const HTML_global_elements = Set([:html, :head, :body])
const HTML_head_elements = Set([:title, :meta])
const HTML_section_elements = Set([
    :body, :article, :section, :nav, :aside,
    :h1, :h2, :h3, :h4, :h5, :h6,
    :hgroup, :header, :footer, :address])
const HTML_grouping_elements = Set([
    :p, :hr, :pre, :blockquote, :ol, :ul,
    :menu, :li, :dl, :dt, :dd, :figure,
    :figcaption, :main, :div])
const HTML_table_elements = Set([
    :table, :caption, :colgroup, :col, :tbody,
    :thead, :tfoot, :tr, :td, :th])
const HTML_text_elements = Set([
    :a, :em, :strong, :small, :s, :cite,
    :q, :dfn, :abbr, :ruby, :rt, :rp,
    :data, :code, :var, :samp, :kbd,
    :sub, :sup, :i, :b, :u, :mark,
    :bdi, :bdo, :span, :br, :wbr])
const HTML_input_elements = Set([
    :button, :checkbox, :color, :date,
    :datetimeLocal, :email, :file, :hidden,
    :image, :month, :number, :password,
    :radio, :range, :reset, :search, :submit,
    :tel, :text, :time, :url, :week])
const HTML_other_elements = Set([:area])

const HTML_elements = union(HTML_global_elements, HTML_head_elements, HTML_section_elements, HTML_grouping_elements,
    HTML_table_elements, HTML_text_elements, HTML_input_elements, HTML_other_elements)

### html attributes

# global attributes

const HTML_global_attributes = Set([
    :accesskey
    :autocapitalize
    :autofocus
    :contenteditable
    :dir
    :draggable
    :enterkeyhint
    :hidden
    :inputmode
    :is
    :itemid
    :itemprop
    :itemref
    :itemscope
    :itemtype
    :lang
    :nonce
    :spellcheck
    :style
    :tabindex
    :title
    :translate
  ])

const HTML_eventhandler_attributes = Set([
    :onabort
    :onauxclick
    :onblur
    :oncancel
    :oncanplay
    :oncanplaythrough
    :onchange
    :onclick
    :onclose
    :oncontextmenu
    :oncopy
    :oncuechange
    :oncut
    :ondblclick
    :ondrag
    :ondragend
    :ondragenter
    :ondragexit
    :ondragleave
    :ondragover
    :ondragstart
    :ondrop
    :ondurationchange
    :onemptied
    :onended
    :onerror
    :onfocus
    :onformdata
    :oninput
    :oninvalid
    :onkeydown
    :onkeypress
    :onkeyup
    :onload
    :onloadeddata
    :onloadedmetadata
    :onloadstart
    :onmousedown
    :onmouseenter
    :onmouseleave
    :onmousemove
    :onmouseout
    :onmouseover
    :onmouseup
    :onpaste
    :onpause
    :onplay
    :onplaying
    :onprogress
    :onratechange
    :onreset
    :onresize
    :onscroll
    :onsecuritypolicyviolation
    :onseeked
    :onseeking
    :onselect
    :onslotchange
    :onstalled
    :onsubmit
    :onsuspend
    :ontimeupdate
    :ontoggle
    :onvolumechange
    :onwaiting
    :onwheel
  ])

const HTML_image_attributes = Set([
    :alt
    :crossorigin
    :height
    :ismap
    :longdesc
    :sizes
    :src
    :srcset
    :usemap
    :width
  ])

#= kinds of HTML Input Attributes =#

const HTML_mostgeneral_input_attributes = Set([
    :autofocus
    :disabled
    :form
    :name
  ])

const HTML_broadlygeneral_input_attributes = Set([
    :autocomplete
    :placeholder
    :readonly
    :required
    :type
    :value
  ])

const HTML_moregeneral_input_attributes = union(HTML_mostgeneral_input_attributes, HTML_broadlygeneral_input_attributes)

const HTML_lessgeneral_input_attributes = Set([
    :list
    :multiple
    :pattern
    :size
    :step
  ])

const HTML_checkbox_input_attributes = Set([
    :checked
  ])

const HTML_file_input_attributes = Set([
    :accept
  ])

const HTML_form_input_attributes = Set([
    :formaction
    :formenctype
    :formmethod
    :formnovalidate
    :formtarget
  ])

 const HTML_image_input_attributes = Set([
    :alt
    :formaction
    :formenctype
    :formmethod
    :formtarget
    :height
    :src
    :width
  ])

const HTML_submit_input_attributes = Set([
    :formenctype
    :formmethod
    :formtarget
  ])

const HTML_input_attributes =
    union(HTML_moregeneral_input_attributes, HTML_lessgeneral_input_attributes,
          HTML_checkbox_input_attributes, HTML_file_input_attributes, HTML_form_input_attributes,
          HTML_image_input_attributes, HTML_submit_input_attributes)

# select

const HTML_select_attributes =
    Set([:autofocus, :disabled, :form, :multiple, :name, :required, :size])

const HTML_option_attributes =
    Set([:disabled, :label, :selected, :value])

const HTML_optgroup_attributes =
    Set([:disabled, :label])

# list

const HTML_ol_attributes =
    Set([:reversed, :start, :type])

const HTML_ul_attributes = Set([])

# style
const HTML_style_attributes =
    Set([:media, :type])

# elements with global_attributes only

const HTML_nav_attributes = Set([])
# table
const HTML_table_attributes = Set([])
const HTML_tbody_attributes = Set([])
const HTML_td_attributes = Set([])
const HTML_tfoot_attributes = Set([])
const HTML_th_attributes = Set([])
const HTML_thead_attributes = Set([])
const HTML_tr_attributes = Set([])


text:
using DFUtils
using Test
using DataFrames

@testset "DFUtils.jl" begin
    df = DataFrame(day = [1,2,1], group = ["a", "a", "b"], value = [1, 2, 3])
    @test isequal(complete(df, :day, :group).value, [1,2,3,missing])
    @test complete(df, :day, :group; replace_missing = 4).value == [1,2,3,4]

    @test toReal([1, 1., 0]) isa Array{Int}
    @test toReal(["1", "1.0", "1."]) isa Array{Int}
    @test toReal([1, "1", "1.0"]) isa Array{Int}
    @test toReal(["1", 1, "", 1.0]) isa Array{Union{Missing, Int}}
    @test toReal(["1.1", 1, 1]) isa Array{Float64}
    @test toReal(["1", 1, "", 1.1]) isa Array{Union{Missing, Float64}}
    @test toReal([missing, 1,1]) isa Array{Union{Missing, Int}}
    @test toReal(["abc", 1,1]) isa Array{Union{Missing, Int}}
    @test toReal([string(typemax(Int) + BigInt(1)), 1, 1., 1.]) isa Array{BigInt}
    @test toReal([string(typemax(Int) + BigInt(1)), 1, 1.1, 1.]) isa Array{BigFloat}
end


text:
__precompile__(true)

"""
A graph and network analysis package for julia.
"""
module Erdos

using Random
using SparseArrays
using LinearAlgebra
using Arpack
import GZip # I/O
import StatsFuns # randgraphs
import Printf
using EzXML # I/O  graphml, gexf  #TODO import instead of using
import Clustering: kmeans # community detection
using IterTools: distinct # all_neighbors
using IterTools: nth # EdgeIter
using Base.Iterators: flatten
import Dates

import DataStructures: MutableBinaryHeap, update!, compare,  # push_relabel
                        PriorityQueue, dequeue!, peek, heappush!, heappop!,
                        enqueue!, compare, top

import SparseArrays: sparse, blockdiag
import LinearAlgebra: issymmetric

import Base: write, ==, <, *, ≈, isless, union, intersect,
            reverse, reverse!, getindex, setindex!, show, print, copy, in,
            sum, size, eltype, length, ndims,
            join, iterate, eltype, get,
            sort, push!, pop!, IteratorSize, values, valtype,
            SizeUnknown, IsInfinite, #iterators
            HasLength, HasShape,     #iterators
            haskey, Matrix #edgemap

#interface
export AGraph, ADiGraph, AGraphOrDiGraph,
ANetwork, ADiNetwork, ANetOrDiNet,
graphtype, digraphtype, edgetype, vertextype,

# edge
AEdge, Edge, reverse,
AIndexedEdge, IndexedEdge, idx,

# core
vertices, edges, src, dst,
graph, digraph,
adjacency_list,
neighbors, in_neighbors, out_neighbors, all_neighbors,
in_edges, out_edges, all_edges,
has_vertex, has_edge, is_directed,
nv, ne, add_edge!, rem_edge!, add_vertex!, add_vertices!,
in_degree, out_degree, degree, has_self_loops, num_self_loops,
rem_vertex!, edge, clean_vertex!,
unsafe_add_edge!, rebuild!,
rem_vertices!, swap_vertices!, pop_vertex!,

# graph types (factory)
reverse!, Graph, DiGraph, DiNetwork, Network,

# distance
eccentricities, eccentricity, diameter, periphery, radius, center,

# distance between graphs
spectral_distance, edit_distance,

# edit path cost functions
MinkowskiCost, BoundedMinkowskiCost,

# operators
complement, union, intersect,
difference, symmetric_difference,
join, tensor_product, cartesian_product, crosspath,
subgraph, egonet, complete, complete!,
subnetwork, contract!,

# graph visit
discover_vertex!, open_vertex!, close_vertex!,
examine_neighbor!, visited_vertices, traverse_graph!, traverse_graph_withlog,

# bfs
gdistances, gdistances!, bfs_tree, is_bipartite, bipartite_map,
has_path,

# dfs
has_cycles, topological_sort_by_dfs, dfs_tree,
#is_tree, #TODO add back

# random
randomwalk, self_avoiding_randomwalk, nonbacktracking_randomwalk,

# connectivity
connected_components, strongly_connected_components, weakly_connected_components,
is_connected, is_strongly_connected, is_weakly_connected, period,
condensation, attracting_components, neighborhood, is_graphical, density,

# maximum_adjacency_visit
MaximumAdjacency, maximum_adjacency_visit,

# a-star, dijkstra, bellman-ford, floyd-warshall
a_star, dijkstra_shortest_paths,
bellman_ford_shortest_paths, has_negative_edge_cycle, enumerate_paths,
floyd_warshall_shortest_paths,
shortest_paths,

# centrality
betweenness_centrality, closeness_centrality, degree_centrality,
in_degree_centrality, out_degree_centrality, katz_centrality, pagerank,
cores, kcore,

# spectral
adjacency_matrix,laplacian_matrix,
CombinatorialAdjacency, nonbacktracking_matrix, incidence_matrix,
nonbacktrack_embedding,

# astar
a_star,

# persistence
readgraph, writegraph, readnetwork, writenetwork,

# flow
maximum_flow, EdmondsKarpAlgorithm, DinicAlgorithm, BoykovKolmogorovAlgorithm, PushRelabelAlgorithm,
multiroute_flow, KishimotoAlgorithm, ExtendedMultirouteFlowAlgorithm,
minimum_cut,

# randgraphs
erdos_renyi, watts_strogatz, random_regular_graph, random_regular_digraph, random_configuration_model,
StochasticBlockModel, make_edgestream, nearbipartiteSBM, blockcounts, blockfractions,
stochastic_block_model, barabasi_albert, barabasi_albert!, static_fitness_model, static_scale_free,

#community
modularity, core_periphery_deg,
local_clustering,local_clustering_coefficient, global_clustering_coefficient, triangles,
label_propagation, maximal_cliques,
community_detection_nback, community_detection_bethe,

#generators
CompleteGraph, StarGraph, PathGraph, WheelGraph, CycleGraph,
CompleteBipartiteGraph, CompleteDiGraph, StarDiGraph, PathDiGraph, Grid,
WheelDiGraph, CycleDiGraph, BinaryTree, DoubleBinaryTree, RoachGraph, CliqueGraph,

#smallgraphs
# graph, digraph

# Euclidean graphs
euclidean_graph,

# Spanning Trees
minimum_spanning_tree, count_spanning_trees,

# dismantling
dismantle_ci, dismantle_ci_init, dismantle_ci_oneiter!,

# maps
AVertexMap, ConstVertexMap, VertexMap,
AEdgeMap, ConstEdgeMap, EdgeMap,
edgemap2adjlist, weights,

# properties
PropertyStore,

graph_property, set_graph_property!, rem_graph_property!, has_graph_property,
vertex_property, add_vertex_property!, rem_vertex_property!, has_vertex_property,
edge_property, add_edge_property!, rem_edge_property!, has_edge_property,

#short forms
gprop, gprop!, rem_gprop!, has_gprop, gprop_names,
vprop, vprop!, rem_vprop!, has_vprop, vprop_names,
eprop, eprop!, rem_eprop!, has_eprop, eprop_names

include("utils.jl")
include("core/interface_graph.jl")
    include("core/interface_network.jl")
    include("core/core.jl")
    include("core/edge.jl")
    include("core/edgeiter.jl")
    include("core/misc.jl")
include("maps/vertexmap.jl")
    include("maps/edgemap.jl")
    include("maps/property_store.jl")
include("operators/operators.jl")
include("traversals/bipartition.jl")
    include("traversals/bfs.jl")
    include("traversals/dfs.jl")
    include("traversals/maxadjvisit.jl")
    include("traversals/randomwalks.jl")
    include("traversals/connectivity.jl")
include("distances/distance.jl")
    include("distances/edit_distance.jl")
include("shortestpaths/shortestpaths.jl")
include("linalg/nonbacktracking.jl")
    include("linalg/spectral.jl")
include("generators/staticgraphs.jl")
    include("generators/randgraphs.jl")
    include("generators/euclideangraphs.jl")
    include("generators/smallgraphs.jl")
include("centrality/betweenness.jl")
    include("centrality/closeness.jl")
    include("centrality/degree.jl")
    include("centrality/katz.jl")
    include("centrality/pagerank.jl")
    include("centrality/kcore.jl")
include("community/modularity.jl")
    include("community/core-periphery.jl")
    include("community/clustering.jl")
    include("community/cliques.jl")
    include("community/detection.jl")
include("flow/maximum_flow.jl")
    include("flow/edmonds_karp.jl")
    include("flow/dinic.jl")
    include("flow/boykov_kolmogorov.jl")
    include("flow/push_relabel.jl")
    include("flow/multiroute_flow.jl")
    include("flow/kishimoto.jl")
    include("flow/ext_multiroute_flow.jl")
include("spanningtrees/spanningtrees.jl")
include("dismantling/ci.jl")
include("factory/graph.jl")
    include("factory/net.jl")
include("persistence/common.jl")
    include("persistence/dot.jl")
    include("persistence/gexf.jl")
    include("persistence/gml.jl")
    include("persistence/graphml.jl")
    include("persistence/net.jl")
    include("persistence/gt.jl")
include("deprecate.jl")
end # module


text:
"""
    gridscore(Xtrain, Ytrain, X, Y; score, fun, pars, verbose = FALSE)
Model validation over a grid of parameters.
* `Xtrain` : Training X-data.
* `Ytrain` : Training Y-data.
* `X` : Validation X-data.
* `Y` : Validation Y-data.
* `score` : Function computing the prediction score (= error rate; e.g. msep).
* `fun` : Function computing the prediction model.
* `pars` : tuple of named vectors (arguments of fun) of same length
    involved in the calculation of the score.
* `verbose` : If true, fitting information are printed.

Compute a prediction score for a given model over a grid of parameters.

The score is computed over X and Y for each combination of the
grid defined in `pars`.

The vectors in `pars` must have same length.
"""
function gridscore(Xtrain, Ytrain, X, Y; score, fun,
        pars, verbose = false)
    q = size(Ytrain, 2)
    nco = length(pars[1]) # nb. combinations in pars
    verbose ? println("-- Nb. combinations = ", nco) : nothing
    res = map(values(pars)...) do v...
        verbose ? println(Pair.(keys(pars), v)...) : nothing
        fm = fun(Xtrain, Ytrain ; Pair.(keys(pars), v)...)
        pred = predict(fm, X).pred
        score(pred, Y)
    end
    verbose ? println("-- End.") : nothing
    nco == 1 ? res = res[1] : res = reduce(vcat, res)
    dat = DataFrame(pars)
    namy = map(string, repeat(["y"], q), 1:q)
    res = DataFrame(res, Symbol.(namy))
    hcat(dat, res)
end

"""
    gridscorelv(Xtrain, Ytrain, X, Y; score, fun, nlv, pars, verbose = FALSE)
* `nlv` : Nb., or collection of nb., of latent variables (LVs).

Same as [`gridscore`](@ref) but specific to (and much faster for) models
using latent variables (e.g. PLSR).

Argument `pars` must not contain `nlv`.
"""
function gridscorelv(Xtrain, Ytrain, X, Y; score, fun, nlv,
        pars = nothing, verbose = false)
    q = size(Ytrain, 2)
    nlv = max(minimum(nlv), 0):maximum(nlv)
    le_nlv = length(nlv)
    if isnothing(pars)
        verbose ? println("-- Nb. combinations = 0.") : nothing
        fm = fun(Xtrain, Ytrain, nlv = maximum(nlv))
        pred = predict(fm, X; nlv = nlv).pred
        le_nlv == 1 ? pred = [pred] : nothing
        res = zeros(le_nlv, q)
        @inbounds for i = 1:le_nlv
            res[i, :] = score(pred[i], Y)
        end
        dat = DataFrame(nlv = nlv)
    else
        nco = length(pars[1])  # nb. combinations in pars
        verbose ? println("-- Nb. combinations = ", nco) : nothing
        res = map(values(pars)...) do v...
            verbose ? println(Pair.(keys(pars), v)...) : nothing
            fm = fun(Xtrain, Ytrain ; nlv = maximum(nlv), Pair.(keys(pars), v)...)
            pred = predict(fm, X; nlv = nlv).pred
            le_nlv == 1 ? pred = [pred] : nothing
            zres = zeros(le_nlv, q)
            for i = 1:le_nlv
                zres[i, :] = score(pred[i], Y)
            end
            zres
        end

        nco == 1 ? res = res[1] : res = reduce(vcat, res)
        ## Make dat
        if le_nlv == 1
            dat = DataFrame(pars)
        else
            zdat = DataFrame(pars)
            dat = list(nco)
            @inbounds for i = 1:nco
                dat[i] = reduce(vcat, fill(zdat[i:i, :], le_nlv))
            end
            dat = reduce(vcat, dat)
        end
        znlv = repeat(nlv, nco)
        dat = hcat(dat, DataFrame(nlv = znlv))
        ## End
    end
    verbose ? println("-- End.") : nothing
    namy = map(string, repeat(["y"], q), 1:q)
    res = DataFrame(res, Symbol.(namy))
    hcat(dat, res)
end

"""
    gridscorelb(Xtrain, Ytrain, X, Y; score, fun, lb, pars, verbose = FALSE)
* `lb` : Value, or collection of values, of the regularization parameter "lambda".

Same as [`gridscore`](@ref) but specific to (and much faster for) models
using ridge regularization (e.g. RR).

Argument `pars` must not contain `lb`.
"""
function gridscorelb(Xtrain, Ytrain, X, Y; score, fun, lb,
        pars = nothing, verbose = false)
    q = size(Ytrain, 2)
    lb = sort(unique(lb))
    le_lb = length(lb)
    if isnothing(pars)
        verbose ? println("-- Nb. combinations = 0.") : nothing
        fm = fun(Xtrain, Ytrain, lb = maximum(lb))
        pred = predict(fm, X; lb = lb).pred
        le_lb == 1 ? pred = [pred] : nothing
        res = zeros(le_lb, q)
        @inbounds for i = 1:le_lb
            res[i, :] = score(pred[i], Y)
        end
        dat = DataFrame(lb = lb)
    else
        nco = length(pars[1])  # nb. combinations in pars
        verbose ? println("-- Nb. combinations = ", nco) : nothing
        res = map(values(pars)...) do v...
            verbose ? println(Pair.(keys(pars), v)...) : nothing
            fm = fun(Xtrain, Ytrain ; lb = maximum(lb), Pair.(keys(pars), v)...)
            pred = predict(fm, X; lb = lb).pred
            le_lb == 1 ? pred = [pred] : nothing
            zres = zeros(le_lb, q)
            for i = 1:le_lb
                zres[i, :] = score(pred[i], Y)
            end
            zres
        end
        nco == 1 ? res = res[1] : res = reduce(vcat, res)
        ## Make dat
        if le_lb == 1
            dat = DataFrame(pars)
        else
            zdat = DataFrame(pars)
            dat = list(nco)
            @inbounds for i = 1:nco
                dat[i] = reduce(vcat, fill(zdat[i:i, :], le_lb))
            end
            dat = reduce(vcat, dat)
        end
        zlb = repeat(lb, nco)
        dat = hcat(dat, DataFrame(lb = zlb))
        ## End
    end
    verbose ? println("-- End.") : nothing
    namy = map(string, repeat(["y"], q), 1:q)
    res = DataFrame(res, Symbol.(namy))
    hcat(dat, res)
end


text:
if basename(pwd()) == "aoc"
    cd("2016/9")
end

function decomplength(str::AbstractString)
    total = 0
    i = 1
    @views while i <= lastindex(str)
        c = str[i]
        if c == '('
            m = match(r"(\d+)x(\d+)\)", str[i+1:end])
            if !isnothing(m)
                x, y = parse.(Int, m)
                i += length(m.match) + 1
                total += y * count(!isspace, str[i:i+x-1])
                i += x - 1
            end
        elseif !isspace(c)
            total += 1
        end
        i += 1
    end
    total
end

@assert decomplength.(eachline("example.txt")) == [6, 7, 9, 11, 6, 18]
decomplength(readline("input.txt"))

function deepdecomplength(str::AbstractString)
    total = 0
    i = 1
    @views while i <= lastindex(str)
        c = str[i]
        if c == '('
            m = match(r"(\d+)x(\d+)\)", str[i+1:end])
            if !isnothing(m)
                x, y = parse.(Int, m)
                i += length(m.match) + 1
                total += y * deepdecomplength(str[i:i+x-1])
                i += x - 1
            end
        elseif !isspace(c)
            total += 1
        end
        i += 1
    end
    total
end

@assert deepdecomplength.(eachline("example2.txt")) == [9, 20, 241920, 445]
deepdecomplength(readline("input.txt"))


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

module CircuitOptimization

import Zygote

using QuantumCircuits.QML.Optimization
using QuantumCircuits.QCircuits.QBase
using QuantumCircuits.QCircuits.Circuit
using QuantumCircuits.QCircuits.Math
using QuantumCircuits.QCircuits.Gates: ParamT

export findparam

const ϵ = 1e-7
const accept_ϵ = 1e-3

function multipleGradientDescent(qc, loss, addGlobalPhase; uϵ=ϵ, N=100, debug=false)
    maxItr=20000
    process = true

    # Get parameters
    params = getparameters(qc)
    if addGlobalPhase
        pushfirst!(params, rand() * 2π)
    end

    # best solution
    bestVal = 1e+5
    bestParams = []

    debug && println("Calculate the derivate.")
    dloss(params) = real(Zygote.gradient(loss, params)[1])
    debug && dloss(params)

    i = 1
    while process
        debug && println("Calculate gradientDescent for iteration $i.")
        val, xparams, itr = gradientDescent(loss, dloss, params, α=1.0, maxItr=maxItr,
                                      argsArePeriodic=true, isExpectedZero=true, ϵ=uϵ)
        debug && println("Solution with error $val found in $itr iteration.")

        # save best solution
        if val < bestVal
            @assert abs(val - loss(xparams)) <= 1e-4 "Unexpected error in multipleGradientDescent: val=$val, loss=$(loss(xparams)), err=$(abs(val - loss(xparams)))"
            bestVal = val
            bestParams = xparams
        end

        if val < uϵ || i > N
            process = false
        else
            if i % 2 == 0
                params = 2π .- params
            else
                params = getRandParameters(qc)
                if addGlobalPhase
                    pushfirst!(params, rand() * 2π)
                end
            end
        end

        # Next iteration
        i += 1
    end

    if addGlobalPhase
        popfirst!(bestParams)
    end

    return bestVal, bestParams
end

function findparam(expmat, qc; N=100, addGlobalPhase=false, debug=false, trystandard=true)
    #println("getparameters")
    #params = getparameters(qc)

    lossPhase(params) = real(unitary_error(expmat, cis(params[1]) * tomatrix(qc, @view(params[2:end]))))
    lossNoPhase(params) = real(unitary_error(expmat, tomatrix(qc, params)))

    if addGlobalPhase
        debug && println("Add global phase to unitary.")
        loss = lossPhase
    else
        loss = lossNoPhase
    end

    val, params = multipleGradientDescent(qc, loss, addGlobalPhase, N=N, debug=debug)
    sge = standardGateError(qc, params)
    params2 = nothing
    val2 = nothing

    if trystandard && val < ϵ && sge > ϵ
        debug && println("====================================================")
        debug && println("The optimal solution uses not standard gate.")

        #params2 = getparameters(qc)
        #debug && println("loss2")
        loss2(params) = real(unitary_error(expmat, tomatrix(qc, params))) + standardGateError(qc, params)
        #debug && loss2(params)

        val2, params2 = multipleGradientDescent(qc, loss2, false, N=N, debug=debug)
        val2std = loss(params2)
        sge2 = standardGateError(qc, params2)

        if (val2std - val) <= accept_ϵ && sge2 < sge
            params = params2
        else
            #@show val2std, val, ϵ, sge2, sge
            println("ERROR: The solution with constraint is worse than without them.")
        end
    end

    debug && println("setparameters")
    setparameters!(qc, params)
    debug && println("end")

    return params, params2, val, val2
end

end  # module CircuitOptimization


text:
struct DefaultWalker <: Walker; end

function trav(w::Walker, x::CSTParser.LeafNode, s, S::State)
    if x isa CSTParser.IDENTIFIER && !S.isquotenode
        binding = find_ref(CSTParser.str_value(x), S)
        push!(S.refs, Reference(x, binding, Location(S.loc.path, S.loc.offset + x.span)))
    end
    S.loc.offset += x.fullspan
end

function trav(w::Walker, x::CSTParser.EXPR{T}, s, S::State) where T <: Union{CSTParser.MacroName,CSTParser.Quote}
    S.loc.offset += x.fullspan
end

function trav(w::Walker, x, s, S::State)
    x isa CSTParser.EXPR{CSTParser.Quotenode} && (S.isquotenode = true)
    for a in x
        get_external_binding(a, s, S)
        if enter_scope(x, s, S)
            create_scope(a, s, S)
            lint_call(w, a, s, S)
            trav(w, a, S.current_scope, S)
            S.current_scope = s
        else
            S.loc.offset += a.fullspan
        end
    end
    x isa CSTParser.EXPR{CSTParser.Quotenode} && (S.isquotenode = false)
    s
end

trav(w::Walker, x::OverlayNode, s, S::State) = trav(w, x.expr, s, S)

function trav(w::Walker, x)
    S = State()
    trav(w, x, S.current_scope, S)
    find_bad_refs(S)
    return S
end

function trav(w::Walker, path::String)
    S = State{FileSystem}(Scope(), Location(path, 0), "", [], [], 0:0, false, Dict(path => File(path, nothing, [])), FileSystem())
    x = CSTParser.parse(readstring(path), true)
    trav(w, x, S.current_scope, S)
    find_bad_refs(S)
    return S
end

function enter_scope(x, s, S)
    isempty(S.target_file) && return true
    ns =  CSTParser.defines_function(x) ||
    # CSTParser.defines_module(x) ||
    CSTParser.defines_macro(x) ||
    CSTParser.defines_datatype(x) ||
    x isa CSTParser.EXPR{CSTParser.Let} ||
    x isa CSTParser.EXPR{CSTParser.Do} ||
    x isa CSTParser.EXPR{CSTParser.While} ||
    x isa CSTParser.EXPR{CSTParser.For} ||
    x isa CSTParser.EXPR{CSTParser.Try} ||
    x isa CSTParser.WhereOpCall ||
    x isa CSTParser.EXPR{CSTParser.Generator} ||
    x isa CSTParser.EXPR{CSTParser.Quote} ||
    CSTParser.defines_anon_function(x) ||
    CSTParser.is_assignment(x) && x.arg1 isa CSTParser.EXPR{CSTParser.Curly}

    if ns
        if S.loc.path == S.target_file
            return true
        else
            return false
        end
    else
        return true
    end

end


text:
export SquareCSop, RombergOP, forward_model

import Base.*
import Base.-
import Base.adjoint
using Flux: cpu, gpu

mutable struct SquareCSop
    A::AbstractArray{Float32,2}
    nx::Int64
    ny::Int64
    nc::Int64
    s::Float32
end

function SquareCSop(nc::Int; s::Float32=3f0)

    A = randn(Float32, nc, nc) + s*I
    A = A/opnorm(A)

    return SquareCSop(A, 1, 1, nc, s)
end

function SquareCSop(A::AbstractArray{Float32,2}, nc::Int; s::Float32=3f0)
    return SquareCSop(A, 1, 1, nc, s)
end

function adjoint(A::SquareCSop)
    return SquareCSop(adjoint(A.A), A.nc; s=A.s)
end

function -(A::SquareCSop)
    return SquareCSop(-A.A, A.nc; s=A.s)
end

struct RombergOP
    A::joLinearOperator{Float32,Float32}
    nx::Int64
    ny::Int64
    nc::Int64
    sr::Float32
end

function RombergOP(nx::Int, ny::Int, nc::Int, sr::Float32)

    M = joRomberg(prod([nx, ny]); DDT=Float32, RDT=Float32)
    R = joRestriction(
        prod([nx, ny]), randperm(prod([nx, ny]))[1:Int(cld(prod([nx, ny]), sr))];
        DDT=Float32, RDT=Float32
    )

    A_flat = R*M
    A = A_flat'*A_flat

    return RombergOP(A, nx, ny, nc, sr)
end

function RombergOP(nc::Int, sr::Float32)

    M = joRomberg(nc; DDT=Float32, RDT=Float32)
    R = joRestriction(
        nc, randperm(nc)[1:Int(cld(nc, sr))];
        DDT=Float32, RDT=Float32
    )

    A_flat = R*M
    A = A_flat'*A_flat

    return RombergOP(A, 1, 1, nc, sr)
end

function adjoint(A::RombergOP)
    return RombergOP(A.A', A.ny, A.nx, A.nc, A.sr)
end

function -(A::RombergOP)
    return RombergOP(-A.A, A.nx, A.ny, A.nc, A.sr)
end

function *(CS::SquareCSop , X::AbstractArray{Float32,4})
    nb = size(X, 4)
    return reshape(CS.A*reshape(X, :, nb), CS.nx, CS.ny, CS.nc, :)
end

function *(CS::RombergOP , X::Array{Float32,4})
    nb = size(X, 4)
    return reshape(CS.A*reshape(X, :, nb), CS.nx, CS.ny, CS.nc, :)
end

function forward_model(X::Array{Float32,4})

    nx, ny, nc, nb = size(X)

    M = joRomberg(prod([nx, ny]); DDT=Float32, RDT=Float32)
    R = joRestriction(
        prod([nx, ny]), randperm(prod([nx, ny]))[1:cld(prod([nx, ny]), 2)];
        DDT=Float32, RDT=Float32
    )

    A_flat = R*M
    A = A_flat'*A_flat

    return reshape(A*reshape(X, :, nb), nx, ny, nc, nb)
end


text:
using Base.Test

using OPCSPs.MVNTools

rho = 0.1*rand()
mu1 = rand()
mu2 = rand()
mu3 = rand()
x3 = rand()

sigma = eye(3)
sigma[1,2] = sigma[2,1] = rho
sigma[1,3] = sigma[3,1] = rho^2

@assert isposdef(sigma)

d = MVN([mu1,mu2,mu3], sigma)
apply_measurement!(d, 3, x3)

@test_approx_eq d.mean[1] mu1+rho^2*(x3-mu3)
@test_approx_eq d.mean[2] mu2
@test_approx_eq d.covariance[1,1] 1-rho^4
@test_approx_eq d.covariance[1,2] rho
@test_approx_eq d.covariance[2,1] rho
@test_approx_eq d.covariance[2,2] 1.0


text:
using Documenter, Muon

makedocs(sitename="Muon Documentation")

deploydocs(
    repo = "github.com/scverse/Muon.jl.git",
    devbranch = "main",
)


text:
using GeometricIntegrators
using GeometricIntegrators.Solutions
using SafeTestsets
using Test

include("dataseries_tests.jl")
include("timeseries_tests.jl")

@safetestset "Deterministic Solutions                                                         " begin include("deterministic_solutions_tests.jl") end
@safetestset "Stochastic Solutions                                                            " begin include("stochastic_solutions_tests.jl") end


text:
using MAT

## Colors
C = [1 0 0;
0 .75 0;
0 0 1;
0 .5 .5;
.5 0 .5;
1 .5 0;
0 0 0;
.75 .75 0;
0.1 0.5 0.6;
0 .75 .75;
.4 .5 .3]

using Plots
# Must be kept in this order, unless you are careful what you do
labels = [1; 2; 3; 12; 18; 17; 25; 15; 24]
names = ["Amazon Fashion", "All Beauty", "Appliances",
"Gift Cards", "Magazine Subsciptions", "Luxury Beauty", "Software",
"Industrial and Scientific","Prime Pantry" ]
sizes = [31; 85; 48; 148; 157; 1581; 802; 5334; 4970]
lnum = length(labels)

# See outer parameters
seednum = 10   # 5
grownum = 300  #[200; 300; 500]
epsis = [10.0 1.0 0.1 0.01 0.001]

plot()
for lab = [1;3;4;5]                 # For each cluster
label = labels[lab]
outputmat = "Output_VaryEps/Label_$label"*"_$seednum"*"_$grownum.mat"

mat = matread(outputmat)
hl_time = mean(mat["hl_time"],dims = 2)
hl = round.(mean(mat["hl_f1"],dims = 2),digits = 2)
b1 = round.(mat["b1_f1"],digits = 2)
b2 = round.(mat["b2_f1"],digits = 2)
r = round.(mat["r_f1"],digits = 2)
re = round.(mat["r_re"],digits = 2)
pr = round.(mat["r_pr"],digits = 2)
nS = mat["newS"]

lw = 3
plot!(epsis',hl,linewidth = lw, grid = false, markershape =:circle,
    legend = :bottomleft, color = RGBA(C[lab,1],C[lab,2],C[lab,3],1),
    xaxis=:log10, label = "label $(label)")

end
savefig("Plots/Smallest_VaryEps_seednum_$seednum"*"_grownum_$grownum.pdf")

max_stars_repo_name: nveldt/HypergraphFlowClustering

max_issues_repo_name: nveldt/HypergraphFlowClustering

max_forks_repo_name: nveldt/HypergraphFlowClustering


text:
"""
    Gradient()

Construct a gradient tensor.
"""
struct Gradient <: Tensor{1,'∇'} end

"""
    system(::Type{Gradient})

The gradient only acts on the coordinates ``r``, ``\\theta``, and
``\\phi``.
"""
system(::Type{Gradient}) = SpatialSubSystem()

"""
    ReducedGradient()

Construct a gradient tensor acting on reduced wavefunctions.
"""
struct ReducedGradient <: Tensor{1,'∂'} end

"""
    system(::Type{ReducedGradient})

The reduced gradient only acts on the coordinates ``r``, ``\\theta``,
and ``\\phi``.
"""
system(::Type{ReducedGradient}) = SpatialSubSystem()

@doc raw"""
    RadialGradientMatrixElement(k)

This represents the matrix element of the radial component of the
gradient operator:

```math
\begin{equation}
\tag{V13.2.23}
\int_0^\infty\diff{r}r^2
\conj{\Psi}_{n'\ell'}(r)
\left(\partial_r+\frac{k}{r}\right)
\Psi_{n\ell}(r)
\end{equation}
```
"""
struct RadialGradientMatrixElement <: OneBodyOperator
    k::Int
end

function Base.show(io::IO, me::RadialGradientMatrixElement)
    if iszero(me.k)
        write(io, "∂ᵣ")
    else
        write(io, "(∂ᵣ ")
        write(io, me.k > 0 ? "+" : "-")
        write(io, " $(abs(me.k))/r)")
    end
end

@tensor(Gradient) do
    begin
        n′ ~ n # The gradient couples orbitals of different n, but
               # there is no selection rule.
        ℓ′ == ℓ ± 1
    end

    raw"""
    rme((n′,ℓ′), ::Gradient, (n,ℓ))

Computes the reduced matrix element of `∇` in terms of
[`RadialGradientMatrixElement`](@ref).
"""
    if ℓ′ == ℓ+1
        √(ℓ+1)*RadialGradientMatrixElement(-ℓ)
    elseif ℓ′==ℓ-1
        - √(ℓ)*RadialGradientMatrixElement(ℓ+1)
    end
end

@tensor(ReducedGradient) do
    begin
        n′ ~ n # The gradient couples orbitals of different n, but
               # there is no selection rule.
        ℓ′ == ℓ ± 1
    end

    raw"""
    rme((n′,ℓ′), ::ReducedGradient, (n,ℓ))

Computes the reduced matrix element of `∂` in terms of
[`RadialGradientMatrixElement`](@ref).
"""
    if ℓ′ == ℓ+1
        √(ℓ+1)*RadialGradientMatrixElement(-(ℓ+1))
    elseif ℓ′==ℓ-1
        - √(ℓ)*RadialGradientMatrixElement(ℓ)
    end
end

module LinearMomenta
import ..cartesian_tensor_component, ..Gradient, ..ReducedGradient

@doc raw"""
    𝐩

The linear momentum operator ``\vec{p}=-\im\nabla``; the elements
correspond to `[px,py,pz]`, i.e. the Cartesian tensor components. Can
be entered as `\bfp`.
"""
const 𝐩 = [-im*cartesian_tensor_component(Gradient(), c)
           for c in [:x, :y, :z]]

@doc raw"""
    𝐩̃

The linear momentum operator ``\vec{p}=-\im\nabla``, but _evaluated in
the basis of reduced wavefunctions_; the elements correspond to
`[px,py,pz]`, i.e. the Cartesian tensor components. Can be entered as
`\bfp\tilde`.
"""
const 𝐩̃ = [-im*cartesian_tensor_component(ReducedGradient(), c)
           for c in [:x, :y, :z]]
export 𝐩, 𝐩̃
end

export Gradient, ReducedGradient


text:
#!/usr/bin/env julia

#load path to qjulia home directory
push!(LOAD_PATH, joinpath(@__DIR__, "..", "core"))
push!(LOAD_PATH, joinpath(@__DIR__, "..", "libs/scidac-routines"))
push!(LOAD_PATH, joinpath(@__DIR__, "..", "libs/quda-routines"))

import QJuliaUtils
import QJuliaEnums
import QJuliaInterface
import QJuliaGaugeUtils
import QJuliaComms
import QUDARoutines
import SCIDACRoutines

using Random
using LinearAlgebra
using MPI

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

load_config_from_file = "/home/astrel/Configs/wl_5p5_x2p38_um0p4125_cfg_1000.lime"

[QJuliaUtils.gridsize_from_cmdline[i] = 1 for i = 1:length(QJuliaUtils.gridsize_from_cmdline)]
QJuliaUtils.get_rank_order("col")

#initialize MPI
MPI.Init()

QUDARoutines.initCommsGridQuda_qj(length(QJuliaUtils.gridsize_from_cmdline), QJuliaUtils.gridsize_from_cmdline, QJuliaUtils.lex_rank_from_coords_t_c, C_NULL)

QUDARoutines.initQuda_qj(0)

Random.seed!(2019)

solve_unit_source = true

const lx = 16
const ly = 16
const lz = 16
const lt = 64
const ls = 1

const dim = 4
const vol = lx*ly*lz*lt*ls

#field latt point sizes
const ssize = 12
const gsize = 9

const splen = vol*ssize
const gflen = vol*gsize

sp_in = Vector{Complex{Float64}}(undef, splen)
sp_ou = Vector{Complex{Float64}}(undef, splen)
gauge = Matrix{Complex{Float64}}(undef, gflen, 4)

if solve_unit_source == false
  QJuliaUtils.gen_random_spinor!(sp_in)
else
  QJuliaUtils.gen_unit_spinor!(sp_ou)
end

gauge_param = QJuliaInterface.QJuliaGaugeParam_qj()

gauge_param.X = (lx, ly, lz, lt)
gauge_param.cpu_prec   = QJuliaEnums.QJULIA_DOUBLE_PRECISION
gauge_param.t_boundary = QJuliaEnums.QJULIA_ANTI_PERIODIC_T
gauge_param.gtype      = QJuliaEnums.QJULIA_WILSON_LINKS
gauge_param.anisotropy = 2.38

gauge_param.cuda_prec                     = QJuliaEnums.QJULIA_DOUBLE_PRECISION
gauge_param.reconstruct                   = QJuliaEnums.QJULIA_RECONSTRUCT_12
gauge_param.cuda_prec_sloppy              = QJuliaEnums.QJULIA_SINGLE_PRECISION
gauge_param.reconstruct_sloppy            = QJuliaEnums.QJULIA_RECONSTRUCT_12
gauge_param.cuda_prec_precondition        = QJuliaEnums.QJULIA_HALF_PRECISION
gauge_param.reconstruct_precondition      = QJuliaEnums.QJULIA_RECONSTRUCT_12
gauge_param.reconstruct_refinement_sloppy = QJuliaEnums.QJULIA_RECONSTRUCT_12
gauge_param.cuda_prec_refinement_sloppy   = QJuliaEnums.QJULIA_HALF_PRECISION

println("======= Gauge parameters =======")
QUDARoutines.printQudaGaugeParam_qj(gauge_param)

inv_param = QJuliaInterface.QJuliaInvertParam_qj()
inv_param.residual_type = QJuliaEnums.QJULIA_L2_RELATIVE_RESIDUAL

println("======= Invert parameters =======")
QUDARoutines.printQudaInvertParam_qj(inv_param)

gauge_load_type = 1
if load_config_from_file != ""
  Xdims = Vector{Cint}(undef, 4)
  for i in 1:length(Xdims); Xdims[i] = gauge_param.X[i] ; end
  qio_prec = Cint(8) #gauge_param.cuda_prec

  SCIDACRoutines.QMPInitComms_qj(0, C_NULL, QJuliaUtils.gridsize_from_cmdline)
  SCIDACRoutines.read_gauge_field_qj(load_config_from_file, gauge, qio_prec, Xdims, 0, C_NULL)
  gauge_load_type = 2
end
QJuliaGaugeUtils.construct_gauge_field!(gauge, gauge_load_type, gauge_param)

gauge_param.gtype       = QJuliaEnums.QJULIA_SU3_LINKS 		#currently cannot set QJULIA_WILSON_LINKS (=QJULIA_SU3_LINKS)  for QUDA

x_face_size = gauge_param.X[2]*gauge_param.X[3]*Int(gauge_param.X[4]/2);
y_face_size = gauge_param.X[1]*gauge_param.X[3]*Int(gauge_param.X[4]/2);
z_face_size = gauge_param.X[1]*gauge_param.X[2]*Int(gauge_param.X[4]/2);
t_face_size = gauge_param.X[1]*gauge_param.X[2]*Int(gauge_param.X[3]/2);

gauge_param.ga_pad = max(x_face_size, y_face_size, z_face_size, t_face_size)

@time QUDARoutines.loadGaugeQuda_qj(gauge, gauge_param)

#Check plaquette
plaq = Array{Cdouble, 1}(undef, 3)
@time QUDARoutines.plaqQuda_qj(plaq)
println("Computed plaquette is ", plaq[1], ", (spatial = ",  plaq[2], ", temporal = ", plaq[3], ")")

mass = -0.4125
#mass = -0.95

#inv_param.inv_type = QJuliaEnums.QJULIA_CG_INVERTER
inv_param.inv_type = QJuliaEnums.QJULIA_PIPEPCG_INVERTER

inv_param.mass = mass
inv_param.kappa = 1.0 / (2.0 * (1.0 + 3.0/gauge_param.anisotropy + mass))
inv_param.maxiter = 2000
inv_param.Nsteps = 2
inv_param.tol = 1e-9
inv_param.reliable_delta = ( inv_param.inv_type == QJuliaEnums.QJULIA_CG_INVERTER ) ? 1e-2 : 1e+2;

inv_param.cuda_prec = QJuliaEnums.QJULIA_DOUBLE_PRECISION
inv_param.cuda_prec_sloppy = QJuliaEnums.QJULIA_SINGLE_PRECISION
#inv_param.cuda_prec_sloppy = QJuliaEnums.QJULIA_DOUBLE_PRECISION
#inv_param.solve_type = QJuliaEnums.QJULIA_DIRECT_PC_SOLVE
# this are default params but we set it explicitly:
inv_param.solution_type = QJuliaEnums.QJULIA_MAT_SOLUTION
inv_param.solve_type    = QJuliaEnums.QJULIA_NORMOP_PC_SOLVE

# set up the preconditioner:
inv_param.inv_type_precondition  = QJuliaEnums.QJULIA_MR_INVERTER
#inv_param.inv_type_precondition  = QJuliaEnums.QJULIA_INVALID_INVERTER
inv_param.schwarz_type           = QJuliaEnums.QJULIA_ADDITIVE_SCHWARZ
inv_param.precondition_cycle     = 1
inv_param.tol_precondition       = 1e-1
inv_param.maxiter_precondition   = 8
inv_param.cuda_prec_precondition = QJuliaEnums.QJULIA_HALF_PRECISION
inv_param.omega                  = 1.0

inv_param.Nsteps     = 2
inv_param.gcrNkrylov = 8

inv_param.ca_basis      = QJuliaEnums.QJULIA_INVALID_BASIS
inv_param.ca_lambda_min = 0.0
inv_param.ca_lambda_max = 1.0

println("Kappa = ",  inv_param.kappa)

mat(out, inp) = QUDARoutines.MatQuda_qj(out, inp, inv_param)

tmpl_src_norm = norm(sp_ou)

if solve_unit_source == true
  mat(sp_in, sp_ou)
  sp_ou .=@. 0.0
end

init_src_norm = norm(sp_in)
println("Initial source norm:: ", init_src_norm, " , template src norm is: ", tmpl_src_norm)

@time QUDARoutines.invertQuda_qj(sp_ou, sp_in, inv_param)

res_vector = Vector{Complex{Float64}}(undef, splen)
mat(res_vector, sp_ou)
res_vector .=@. sp_in - res_vector
res_norm = norm(res_vector)
println("Residual norm:: ", res_norm)

QUDARoutines.endQuda_qj()

MPI.Finalize()

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
# uses the Fast Iterative Soft-Thresholding Algorithm (FISTA) to
# minimize f(x) + g(x) = ½(|Ax - y|² + β|x|²) + ½α|Ψx|₁

fista(A, y, α, β, iters, tol, reg::Regularizer) = fista(A, y, α, β, iters, tol, reg::Regularizer, zeros(size(A)[2]))

# TODO: structure similarly to IterativeSolvers
function fista(A, y, α, β, iters, tol, reg::Regularizer, xstart)
    n, p = size(A)

    maxeig = abs(powm(A' * A, maxiter=100)[1]) # TODO: handle number of eigiters
    L = maxeig + β # Lipschitz constant of f (strongly convex term)
    η = 1 / (L * lipschitz_scale(reg))

    x = xstart[:]
    z = x[:]
    xold = similar(x)
    res = similar(y)
    grad = similar(x)
    t = 1

    iters_done = iters
    xupdates = Float64[]
    convdists = Float64[]
    evals = Float64[]

    Ψ = transform_op(reg)
    Fold = Inf

    for i = 1:iters

        xold .= x

        res .= mul!(res, A, z) .- y # TODO: maybe use five-arg mul! instead. (but IterativeSolvers sticks to three-arg)
        grad .= mul!(grad, A', res) .+ β .* z

        x .= z .- η .* grad
        proximal!(x, 1/2 * η * α, reg)

        restart = dot(z .- x, x .- xold)
        restart > 0 && (t = 1)

        told = t
        t = 1/2 * (1 + √(1 + 4t^2))

        z .= x .+ (told - 1)/t .* (x .- xold)

        xupdate = norm(z .- xold) / norm(xold)
        append!(xupdates, xupdate)

        if xupdate < tol
            iters_done = i
            break
        end
    end

    x, (;iters=iters_done, final_tol=norm(x .- xold) / norm(x), xupdates)
end


text:
using Plots ; pyplot()

occupancyAnalytic(n,r) =  sum([(-1)^k*binomial(n,k)*(1 - k/n)^r for k in 0:n])

function occupancyMC(n,r,N)
    fullCount = 0
    for _ in 1:N
        envelopes = zeros(Int,n)
        for k in 1:r
            target = rand(1:n)
            envelopes[target] += 1
        end
        numFilled = sum(envelopes .> 0)
        if numFilled == n
            fullCount += 1
        end
    end
    return fullCount/N
end

max_n, N, Kvals = 100, 10^3, [2,3,4]

analytic = [[occupancyAnalytic(big(n),big(k*n)) for n in 1:max_n] for k in Kvals]
monteCarlo = [[occupancyMC(n,k*n,N) for n in 1:max_n] for k in Kvals]

plot(1:max_n, analytic, c=[:blue :red :green],
	label=["K=2" "K=3" "K=4"])
scatter!(1:max_n, monteCarlo, mc=:black, shape=:+,
	label="", xlims=(0,max_n),ylims=(0,1),
	xlabel="n Envelopes", ylabel="Probability", legend=:topright)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
#Datasets used for testing the functions

x=[1.672064, 0.7876588, 0.317322, 0.9721646, 0.4004206, 1.665123, 3.059971, 0.09459603, 1.27424, 3.522148,
   0.8211308, 1.328767, 2.825956, 0.1102891, 0.06314285, 2.59152, 8.624108, 0.6516885, 5.770285, 0.5154299]

srand(1)
x2=randn(21)
srand(10)
x3=randn(100, 1, 2);

srand(1)
m = randn(20);
srand(2)
y = randn(20) + 2.0;

srand(2)
xnew = rand(6, 4)
xold = reshape(x, 5, 4)

srand(1)
c = cnorm(2000)

srand(3)
y2 = randn(20)+3;


text:
# AUTO GENERATED FILE - DO NOT EDIT

export row

"""
    row(;kwargs...)
    row(children::Any;kwargs...)
    row(children_maker::Function;kwargs...)

A Row component.
Grid Row
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional): The children of the element
- `style` (Dict; optional): The inline styles
- `id` (String; optional): The ID of this component, used to identify dash components
in callbacks. The ID needs to be unique across all of the
components in an app.
- `className` (String; optional): The class of the element
"""
function row(; kwargs...)
        available_props = Symbol[:children, :style, :id, :className]
        wild_props = Symbol[]
        return Component("row", "Row", "dash_carbon_components", available_props, wild_props; kwargs...)
end

row(children::Any; kwargs...) = row(;kwargs..., children = children)
row(children_maker::Function; kwargs...) = row(children_maker(); kwargs...)

max_stars_repo_name: Matheus-Rangel/dash-carbon-components

max_issues_repo_name: Matheus-Rangel/dash-carbon-components

max_forks_repo_name: Matheus-Rangel/dash-carbon-components


text:
@test   Options() == Options( 1e9, :ScatteringParameters, :MagnitudeAngle, 50.0 )
@test   TouchstoneData( DataPoint[], NoiseDataPoint[], Options( 1e9, :ScatteringParameters, :MagnitudeAngle, 50.0 ) ) == TouchstoneData( DataPoint[] )

@test DataPoint( 1.0, [ 1 2; 3 4 ] ) ==  DataPoint( 1.0, [ 1 2; 3 4 ] )
@test DataPoint( 2.0, [ 1 2; 3 4 ] ) !=  DataPoint( 1.0, [ 1 2; 3 4 ] )
@test DataPoint( 1.0, [ 2 2; 3 4 ] ) !=  DataPoint( 1.0, [ 1 2; 3 4 ] )

@test TouchstoneData( [ DataPoint( 1.0, [ 1 2; 3 4] ) ] ) == TouchstoneData( [ DataPoint( 1.0, [ 1 2; 3 4]) ] )
@test TouchstoneData( [ DataPoint( 2.0, [ 1 2; 3 4] ) ] ) != TouchstoneData( [ DataPoint( 1.0, [ 1 2; 3 4]) ] )
@test TouchstoneData( [ DataPoint( 1.0, [ 2 2; 3 4] ) ] ) != TouchstoneData( [ DataPoint( 1.0, [ 1 2; 3 4]) ] )

@test TouchstoneData( [
    DataPoint( 1.0, [ 1 2; 3 4]),
    DataPoint( 2.0, [ 1 2; 3 4])
  ] ) == TouchstoneData( [
    DataPoint( 1.0, [ 1 2; 3 4]),
    DataPoint( 2.0, [ 1 2; 3 4])
  ] )

@test TouchstoneData( [
    DataPoint( 1.0, [ 1 2; 3 4]),
    DataPoint( 2.0, [ 1 2; 3 4])
  ] ) != TouchstoneData( [
    DataPoint( 1.0, [ 1 2; 3 4]),
    DataPoint( 3.0, [ 1 2; 3 4])
  ] )

@test_throws ErrorException(
    "Frequencies of data vector not in ascending order."
  ) TouchstoneData( [
    DataPoint( 2.0, [ 1 2; 3 4]),
    DataPoint( 1.0, [ 1 2; 3 4])
  ] )

@test NoiseDataPoint( 1.0, 2.0, 3.0 + 4im, 5.0 ) ==  NoiseDataPoint( 1.0, 2.0, 3.0 + 4im, 5.0 )

@test NoiseDataPoint( 1.0, 2.0, 3.0 + 4im, 5.0 ) !=  NoiseDataPoint( 2.0, 2.0, 3.0 + 4im, 5.0 )
@test NoiseDataPoint( 1.0, 2.0, 3.0 + 4im, 5.0 ) !=  NoiseDataPoint( 2.0, 2.0, 3.0 + 4im, 5.0 )
@test NoiseDataPoint( 1.0, 2.0, 3.0 + 4im, 5.0 ) !=  NoiseDataPoint( 2.0, 1.0, 3.0 + 4im, 5.0 )
@test NoiseDataPoint( 1.0, 2.0, 3.0 + 4im, 5.0 ) !=  NoiseDataPoint( 2.0, 2.0, 5.0 + 4im, 5.0 )
@test NoiseDataPoint( 1.0, 2.0, 3.0 + 4im, 5.0 ) !=  NoiseDataPoint( 2.0, 2.0, 3.0 + 4im, 4.0 )

@test TouchstoneData( [
    DataPoint( 1.0, [ 1 2; 3 4])
  ], [
    NoiseDataPoint( 1.0, 2.0, 3.0 + 4im, 5.0 ),
    NoiseDataPoint( 2.0, 2.0, 3.0 + 4im, 5.0 )
  ] ) == TouchstoneData( [
    DataPoint( 1.0, [ 1 2; 3 4])
  ], [
    NoiseDataPoint( 1.0, 2.0, 3.0 + 4im, 5.0 ),
    NoiseDataPoint( 2.0, 2.0, 3.0 + 4im, 5.0 )
  ] )

@test TouchstoneData( [
    DataPoint( 2.0, [ 1 2; 3 4])
  ], [
    NoiseDataPoint( 1.0, 2.0, 3.0 + 4im, 5.0 )
  ] ) != TouchstoneData( [
    DataPoint( 2.0, [ 1 2; 3 4])
  ], [
    NoiseDataPoint( 2.0, 2.0, 3.0 + 4im, 5.0 )
  ] )

@test_throws ErrorException(
    "Frequencies of noise data vector not in ascending order."
  ) TouchstoneData( [
    DataPoint( 2.0, [ 1 2; 3 4])
  ], [
    NoiseDataPoint( 2.0, 2.0, 3.0 + 4im, 5.0 ),
    NoiseDataPoint( 1.0, 2.0, 3.0 + 4im, 5.0 )
  ] )

@test_throws ErrorException(
    "Last frequency of data vector is smaller than first frequency of noise data vector."
  ) TouchstoneData( [
    DataPoint( 1.0, [ 1 2; 3 4])
  ], [
    NoiseDataPoint( 2.0, 2.0, 3.0 + 4im, 5.0 )
  ] )