text:
"""
```
measurement(m::Model1002{T}, TTT::Matrix{T}, RRR::Matrix{T},
    CCC::Vector{T}; reg::Int = 1) where {T<:AbstractFloat}
```

Assign measurement equation

```
y_t = ZZ*s_t + DD + u_t
```

where

```
Var(ϵ_t) = QQ
Var(u_t) = EE
Cov(ϵ_t, u_t) = 0
```
"""
function measurement(m::Model1002{T},
                     TTT::Matrix{T},
                     RRR::Matrix{T},
                     CCC::Vector{T}; reg::Int = 1,
                     TTTs::Vector{<: AbstractMatrix{T}} = Matrix{T}[],
                     CCCs::Vector{<: AbstractVector{T}} = Vector{T}[],
                     information_set::UnitRange = reg:reg,
                     memo::Union{ForwardMultipleExpectationsMemo, Nothing} = nothing) where {T <: AbstractFloat}

    endo     = m.endogenous_states
    endo_new = m.endogenous_states_augmented
    exo      = m.exogenous_shocks
    obs      = m.observables

    _n_observables = n_observables(m)
    _n_states = n_states_augmented(m)
    _n_shocks_exogenous = n_shocks_exogenous(m)

    ZZ = zeros(_n_observables, _n_states)
    DD = zeros(_n_observables)
    EE = zeros(_n_observables, _n_observables)
    QQ = zeros(_n_shocks_exogenous, _n_shocks_exogenous)

    for para in m.parameters
        if !isempty(para.regimes)
            if haskey(get_settings(m), :model2para_regime) && haskey(get_setting(m, :model2para_regime), para.key)
                ModelConstructors.toggle_regime!(para, reg, get_setting(m, :model2para_regime)[para.key])
            else
                ModelConstructors.toggle_regime!(para, reg)
            end
        end
    end

    # Set up for calculating k-periods ahead expectations and expected sums
    permanent_t = length(information_set[findfirst(x -> x == reg, information_set):end]) - 1 + reg
    if information_set[1] == information_set[end] # this condition is probably wrong
        # In this case, we do not need to pass the TTTs, CCCs in, so we redefine them as empty.
        # This step is also necessary to ensure the memo is properly used.
        TTTs = Matrix{T}[]
        CCCs = Vector{T}[]

        # for now, assume that in this case, we do not want to use the memo.
        # The complicating problem is that, in the historical regimes, we don't want to use the memo
        # (at least, not the one that is currently being automatically constructed).
        # But in the forecast regimes, we do want to use the memo if we ever reach a case where
        # the information set becomes myopic.
        memo = nothing
        # TODO: maybe instead of emptying TTTs, CCCs, we add
        # a step to recompute the memo since we will still likely be recalculating
        # products/powers of TTT multiple times that could be pre-computed.
    end
    use_fwd_exp_sum = haskey(get_settings(m), :use_forward_expected_sum_memo) && get_setting(m, :use_forward_expected_sum_memo)
    use_fwd_exp     = haskey(get_settings(m), :use_forward_expectations_memo) && get_setting(m, :use_forward_expectations_memo)

    flex_ait_2020Q3 = haskey(get_settings(m), :flexible_ait_policy_change) &&
        get_setting(m, :flexible_ait_policy_change)

    # With flexible_ait_policy_change: there is a regime-break in 2020:Q3, so before 2020:Q2,
    # the final regime should be considered to be the regime corresponding to 2020:Q2, namely for regimes
    # before 2020:Q3, the measurement equation should reflect the belief that the "permanent" regime is the
    # regime before 2020:Q3 starts.
    if flex_ait_2020Q3
        if get_setting(m, :regime_dates)[4] != get_setting(m, :flexible_ait_policy_change_date)
            reg_2020Q3 = 0
            for i in 1:get_setting(m, :n_regimes)
                if get_setting(m, :regime_dates)[i] == get_setting(m, :flexible_ait_policy_change_date)
                    reg_2020Q2 = i
                    break
                end
            end
            @assert reg_2020Q2 != 0 "The setting :regime_dates does not contain the date 2020:Q3, which is required if the setting :flexible_ait_policy_change is true."
        else
            reg_2020Q3 = 4
        end

        permanent_t = reg < reg_2020Q3 ? reg_2020Q3 - 1 : permanent_t
    end

    # Remove integrated states (e.g. states w/unit roots)
    no_integ_inds = inds_states_no_integ_series(m)
    if ((haskey(get_settings(m), :add_altpolicy_pgap) && haskey(get_settings(m), :pgap_type)) ?
        (get_setting(m, :add_altpolicy_pgap) && get_setting(m, :pgap_type) == :ngdp) : false) ||
        (haskey(get_settings(m), :ait_Thalf) ? (exp(log(0.5) / get_setting(m, :ait_Thalf)) ≈ 0.) : false)
        no_integ_inds = setdiff(no_integ_inds, [m.endogenous_states[:pgap_t]])
    end
    if (haskey(get_settings(m), :gdp_Thalf) ? (exp(log(0.5) / get_setting(m, :gdp_Thalf)) ≈ 0.) : false)
        no_integ_inds = setdiff(no_integ_inds, [m.endogenous_states[:ygap_t]])
    end
    if haskey(get_settings(m), :ρ_rw) ? (get_setting(m, :ρ_rw) ≈ 1.) : false
        no_integ_inds = setdiff(no_integ_inds, [m.endogenous_states[:rw_t]])
    end
    if haskey(get_settings(m), :rw_ρ_smooth) ? (get_setting(m, :rw_ρ_smooth) ≈ 1.) : false
        no_integ_inds = setdiff(no_integ_inds, [m.endogenous_states[:Rref_t]])
    end
    integ_series = length(no_integ_inds) != n_states_augmented(m) # Are the series integrated?

    ## GDP growth - Quarterly!
    ZZ[obs[:obs_gdp], endo[:y_t]]          = 1.0
    ZZ[obs[:obs_gdp], endo_new[:y_t1]]     = -1.0
    ZZ[obs[:obs_gdp], endo[:z_t]]          = 1.0
    ZZ[obs[:obs_gdp], endo_new[:e_gdp_t]]  = 1.0
    ZZ[obs[:obs_gdp], endo_new[:e_gdp_t1]] = -m[:me_level]
    DD[obs[:obs_gdp]]                      = 100*(exp(m[:z_star])-1)

    if haskey(get_settings(m), :add_iid_cond_obs_gdp_meas_err) ?
        get_setting(m, :add_iid_cond_obs_gdp_meas_err) : false
        ZZ[obs[:obs_gdp], endo_new[:e_condgdp_t]] = 1.
    end

    ## GDI growth- Quarterly!
    ZZ[obs[:obs_gdi], endo[:y_t]]          = m[:γ_gdi]
    ZZ[obs[:obs_gdi], endo_new[:y_t1]]     = -m[:γ_gdi]
    ZZ[obs[:obs_gdi], endo[:z_t]]          = m[:γ_gdi]
    ZZ[obs[:obs_gdi], endo_new[:e_gdi_t]]  = 1.0
    ZZ[obs[:obs_gdi], endo_new[:e_gdi_t1]] = -m[:me_level]
    DD[obs[:obs_gdi]]                      = 100*(exp(m[:z_star])-1) + m[:δ_gdi]

    ## Hours growth
    ZZ[obs[:obs_hours], endo[:L_t]] = 1.0
    DD[obs[:obs_hours]]             = m[:Lmean]

    ## Labor Share/real wage growth
    if subspec(m) in ["ss16", "ss17"]
        # log(labor_share) = log(wage) + log(hours) - log(GDP)
        ZZ[obs[:obs_laborshare], endo[:w_t]] = 1.
        ZZ[obs[:obs_laborshare], endo[:L_t]] = 1.
        DD[obs[:obs_laborshare]] = 100. * log(m[:wstar] * m[:Lstar] / m[:ystar])
        ZZ[obs[:obs_laborshare], endo[:y_t]] = -1.
    else
        ZZ[obs[:obs_wages], endo[:w_t]]      = 1.0
        ZZ[obs[:obs_wages], endo_new[:w_t1]] = -1.0
        ZZ[obs[:obs_wages], endo[:z_t]]      = 1.0
        DD[obs[:obs_wages]]                  = 100*(exp(m[:z_star])-1)
    end

    ## Inflation (GDP Deflator)
    ZZ[obs[:obs_gdpdeflator], endo[:π_t]]            = m[:Γ_gdpdef]
    ZZ[obs[:obs_gdpdeflator], endo_new[:e_gdpdef_t]] = 1.0
    DD[obs[:obs_gdpdeflator]]                        = 100*(m[:π_star]-1) + m[:δ_gdpdef]

    ## Inflation (Core PCE)
    ZZ[obs[:obs_corepce], endo[:π_t]]             = 1.0
    ZZ[obs[:obs_corepce], endo_new[:e_corepce_t]] = 1.0
    DD[obs[:obs_corepce]]                         = 100*(m[:π_star]-1)

    if haskey(get_settings(m), :add_iid_cond_obs_corepce_meas_err) ?
        get_setting(m, :add_iid_cond_obs_corepce_meas_err) : false
        ZZ[obs[:obs_corepce], endo_new[:e_condcorepce_t]] = 1.
    end

    ## Nominal interest rate
    ZZ[obs[:obs_nominalrate], endo[:R_t]] = 1.0
    DD[obs[:obs_nominalrate]]             = m[:Rstarn]

    ## Consumption Growth
    ZZ[obs[:obs_consumption], endo[:c_t]]      = 1.0
    ZZ[obs[:obs_consumption], endo_new[:c_t1]] = -1.0
    ZZ[obs[:obs_consumption], endo[:z_t]]      = 1.0
    DD[obs[:obs_consumption]]                  = 100*(exp(m[:z_star])-1)

    ## Investment Growth
    ZZ[obs[:obs_investment], endo[:i_t]]      = 1.0
    ZZ[obs[:obs_investment], endo_new[:i_t1]] = -1.0
    ZZ[obs[:obs_investment], endo[:z_t]]      = 1.0
    DD[obs[:obs_investment]]                  = 100*(exp(m[:z_star])-1)

    ## Spreads
    ZZ[obs[:obs_spread], endo[:ERktil_t]] = 1.0
    ZZ[obs[:obs_spread], endo[:R_t]]       = -1.0
    DD[obs[:obs_spread]]                   = 100*log(m[:spr])

    ## 10 yrs infl exp
    TTT10, CCC10 = k_periods_ahead_expected_sums(TTT, CCC, TTTs, CCCs, reg, 40, permanent_t;
                                                 integ_series = integ_series,
                                                 memo = use_fwd_exp_sum ? memo : nothing)
    TTT10        = TTT10 ./ 40. # divide by 40 to average across 10 years
    CCC10        = CCC10 ./ 40.
    ZZ[obs[:obs_longinflation], :] = view(TTT10, endo[:π_t], :)
    # ZZ[obs[:obs_longinflation], :] = TTT10[endo[:π_t], :]
    DD[obs[:obs_longinflation]]    = 100*(m[:π_star]-1) + CCC10[endo[:π_t]]

    ## Long Rate
    # ZZ_long_rate = ZZ[obs[:obs_nominalrate], :]'
    # ZZ[obs[:obs_longrate], :]                 = ZZ_long_rate * TTT10 # TODO: this is slightly inefficient, do what we do with long inflation
    ZZ[obs[:obs_longrate], :]                 =  view(TTT10, endo[:R_t], :) # TODO: this is slightly inefficient, do what we do with long inflation
    ZZ[obs[:obs_longrate], endo_new[:e_lr_t]] = 1.0
    DD[obs[:obs_longrate]]                    = m[:Rstarn] + CCC10[endo[:R_t]]
    # DD[obs[:obs_longrate]]                    = m[:Rstarn] + ZZ_long_rate * CCC10

    ## TFP
    ZZ[obs[:obs_tfp], endo[:z_t]] = (1-m[:α])*m[:Iendoα] + 1*(1-m[:Iendoα])
    if subspec(m) in ["ss14", "ss15", "ss16", "ss18", "ss19"]
        ZZ[obs[:obs_tfp], endo_new[:e_tfp_t]]  = 1.0
        ZZ[obs[:obs_tfp], endo_new[:e_tfp_t1]] = -m[:me_level]
    else
        ZZ[obs[:obs_tfp], endo_new[:e_tfp_t]] = 1.0
    end
    if !(subspec(m) in ["ss15", "ss16"])
        ZZ[obs[:obs_tfp], endo[:u_t]]       = m[:α]/( (1-m[:α])*(1-m[:Iendoα]) + 1*m[:Iendoα] )
        ZZ[obs[:obs_tfp], endo_new[:u_t1]]  = -(m[:α]/( (1-m[:α])*(1-m[:Iendoα]) + 1*m[:Iendoα]) )
    end

    # ygap and pgap for Flexible AIT rule
    if (haskey(get_settings(m), :add_initialize_pgap_ygap_pseudoobs) ? get_setting(m, :add_initialize_pgap_ygap_pseudoobs) : false)
        ZZ[obs[:obs_pgap], endo[:pgap_t]] = 1.
        ZZ[obs[:obs_ygap], endo[:ygap_t]] = 1.
    end

    ## Set up structural shocks covariance matrix
    QQ[exo[:g_sh], exo[:g_sh]]             = m[:σ_g]^2
    QQ[exo[:b_sh], exo[:b_sh]]             = m[:σ_b]^2
    QQ[exo[:μ_sh], exo[:μ_sh]]             = m[:σ_μ]^2
    QQ[exo[:ztil_sh], exo[:ztil_sh]]       = m[:σ_ztil]^2
    QQ[exo[:λ_f_sh], exo[:λ_f_sh]]         = m[:σ_λ_f]^2
    QQ[exo[:λ_w_sh], exo[:λ_w_sh]]         = m[:σ_λ_w]^2
    QQ[exo[:rm_sh], exo[:rm_sh]]           = m[:σ_r_m]^2
    QQ[exo[:σ_ω_sh], exo[:σ_ω_sh]]         = m[:σ_σ_ω]^2
    QQ[exo[:μ_e_sh], exo[:μ_e_sh]]         = m[:σ_μ_e]^2
    QQ[exo[:γ_sh], exo[:γ_sh]]             = m[:σ_γ]^2
    QQ[exo[:π_star_sh], exo[:π_star_sh]]   = m[:σ_π_star]^2
    QQ[exo[:lr_sh], exo[:lr_sh]]           = m[:σ_lr]^2
    QQ[exo[:zp_sh], exo[:zp_sh]]           = m[:σ_z_p]^2
    QQ[exo[:tfp_sh], exo[:tfp_sh]]         = m[:σ_tfp]^2
    QQ[exo[:gdpdef_sh], exo[:gdpdef_sh]]   = m[:σ_gdpdef]^2
    QQ[exo[:corepce_sh], exo[:corepce_sh]] = m[:σ_corepce]^2
    QQ[exo[:gdp_sh], exo[:gdp_sh]]         = m[:σ_gdp]^2
    QQ[exo[:gdi_sh], exo[:gdi_sh]]         = m[:σ_gdi]^2

    if subspec(m) in ["ss59", "ss60", "ss61", "ss62", "ss63", "ss64", "ss65", "ss66", "ss67", "ss68", "ss69", "ss70", "ss71", "ss72", "ss73", "ss74", "ss75", "ss76", "ss77", "ss78", "ss79", "ss80", "ss81", "ss82", "ss83", "ss84", "ss85", "ss86"]
        QQ[exo[:ziid_sh], exo[:ziid_sh]]   = m[:σ_ziid]^2
        QQ[exo[:biidc_sh], exo[:biidc_sh]] = m[:σ_biidc]^2
        QQ[exo[:φ_sh], exo[:φ_sh]]         = m[:σ_φ]^2
    end
    if subspec(m) in ["ss86"]
        QQ[exo[:λ_f_iid_sh], exo[:λ_f_iid_sh]] = m[:σ_λ_f_iid]^2
    end

    if subspec(m) in ["ss67", "ss68", "ss69", "ss70", "ss71", "ss72", "ss73", "ss74", "ss75", "ss76", "ss77", "ss78", "ss80", "ss82", "ss83"]
        QQ[exo[:g_covid_sh], exo[:g_covid_sh]]   = m[:σ_g_covid]^2
        QQ[exo[:μ_covid_sh], exo[:μ_covid_sh]]   = m[:σ_μ_covid]^2
        QQ[exo[:λ_f_covid_sh], exo[:λ_f_covid_sh]]   = m[:σ_λ_f_covid]^2
        QQ[exo[:σ_ω_covid_sh], exo[:σ_ω_covid_sh]]   = m[:σ_σ_ω_covid]^2
        QQ[exo[:lr_covid_sh], exo[:lr_covid_sh]]   = m[:σ_lr_covid]^2
        QQ[exo[:tfp_covid_sh], exo[:tfp_covid_sh]]   = m[:σ_tfp_covid]^2
        QQ[exo[:gdp_covid_sh], exo[:gdp_covid_sh]]   = m[:σ_gdp_covid]^2
        QQ[exo[:gdi_covid_sh], exo[:gdi_covid_sh]]   = m[:σ_gdi_covid]^2
    end
    if subspec(m) in ["ss69", "ss70", "ss73", "ss74", "ss77", "ss78"]
        QQ[exo[:zp_covid_sh], exo[:zp_covid_sh]]   = m[:σ_z_p_covid]^2
    end

    if subspec(m) in ["ss62"] # TODO: add ss63
        ## Set up structural shocks covariance matrix
        QQ[exo[:g_sh], exo[:g_sh]]             *= m[:damp_standard_shocks]^2
        QQ[exo[:b_sh], exo[:b_sh]]             *= m[:damp_standard_shocks]^2
        QQ[exo[:μ_sh], exo[:μ_sh]]             *= m[:damp_standard_shocks]^2
        QQ[exo[:ztil_sh], exo[:ztil_sh]]       *= m[:damp_standard_shocks]^2
        QQ[exo[:λ_f_sh], exo[:λ_f_sh]]         *= m[:damp_standard_shocks]^2
        QQ[exo[:λ_w_sh], exo[:λ_w_sh]]         *= m[:damp_standard_shocks]^2
        QQ[exo[:rm_sh], exo[:rm_sh]]           *= m[:amplify_σ_r_m]^2
        QQ[exo[:σ_ω_sh], exo[:σ_ω_sh]]         *= m[:damp_standard_shocks]^2
        QQ[exo[:μ_e_sh], exo[:μ_e_sh]]         *= m[:damp_standard_shocks]^2
        QQ[exo[:γ_sh], exo[:γ_sh]]             *= m[:damp_standard_shocks]^2
        QQ[exo[:π_star_sh], exo[:π_star_sh]]   *= m[:damp_standard_shocks]^2
        QQ[exo[:lr_sh], exo[:lr_sh]]           *= m[:damp_standard_shocks]^2
        # QQ[exo[:zp_sh], exo[:zp_sh]]           = m[:σ_z_p]^2 # σ_z_p is handled on its own via regime-switching
        QQ[exo[:tfp_sh], exo[:tfp_sh]]         *= m[:damp_standard_shocks]^2
        QQ[exo[:gdpdef_sh], exo[:gdpdef_sh]]   *= m[:amplify_inflation_me]^2
        QQ[exo[:corepce_sh], exo[:corepce_sh]] *= m[:amplify_inflation_me]^2
        QQ[exo[:gdp_sh], exo[:gdp_sh]]         *= m[:damp_standard_shocks]^2
        QQ[exo[:gdi_sh], exo[:gdi_sh]]         *= m[:damp_standard_shocks]^2
    end

    # ygap and pgap shocks
    if (haskey(get_settings(m), :add_initialize_pgap_ygap_pseudoobs) ? get_setting(m, :add_initialize_pgap_ygap_pseudoobs) : false)
        QQ[exo[:pgap_sh], exo[:pgap_sh]] = m[:σ_pgap]^2
        QQ[exo[:ygap_sh], exo[:ygap_sh]] = m[:σ_ygap]^2
    end

    # measurement errors for "conditional" observations of GDP
    if haskey(get_settings(m), :add_iid_cond_obs_gdp_meas_err) ?
        get_setting(m, :add_iid_cond_obs_gdp_meas_err) : false
        QQ[exo[:condgdp_sh], exo[:condgdp_sh]] = m[:σ_condgdp] ^ 2
    end

    if haskey(get_settings(m), :add_iid_anticipated_obs_gdp_meas_err) ?
        get_setting(m, :add_iid_anticipated_obs_gdp_meas_err) : false
        QQ[exo[:gdpexp_sh], exo[:gdpexp_sh]] = m[:σ_gdpexp] ^ 2
    end

    # measurement errors for "conditional" observations of Core PCE
    if haskey(get_settings(m), :add_iid_cond_obs_corepce_meas_err) ?
        get_setting(m, :add_iid_cond_obs_corepce_meas_err) : false
        QQ[exo[:condcorepce_sh], exo[:condcorepce_sh]] = m[:σ_condcorepce] ^ 2
    end

    # Automated addition of anticipated shocks to QQ
    for (k, v) in get_setting(m, :antshocks)
        for i = 1:v
            QQ[exo[Symbol(k, "_shl$i")], exo[Symbol(k, "_shl$i")]] = m[Symbol("σ_", k, "$i")]^2
        end
    end

    ## Anticipated observables
    use_current_regime = haskey(get_settings(m), :measurement_use_current_regime_matrices) ?
        get_setting(m, :measurement_use_current_regime_matrices) : true

    # Anticipated monetary policy shocks
    # ZZ_obs_nomrate = ZZ[obs[:obs_nominalrate], :]'
    for i = 1:n_mon_anticipated_shocks(m)
        TTT_accum, CCC_accum = k_periods_ahead_expectations(TTT, CCC, TTTs, CCCs, reg, i, permanent_t;
                                                            integ_series = integ_series,
                                                            memo = (isnothing(memo) || !use_fwd_exp) ? nothing :
                                                            ForwardExpectationsMemo(memo.time_varying_memo[min(reg + i, permanent_t, 17)],
                                                                                memo.permanent_memo))

        # ZZ[obs[Symbol("obs_nominalrate$i")], :] = ZZ_obs_nomrate * TTT_accum
        ZZ[obs[Symbol("obs_nominalrate$i")], :] = view(TTT_accum, endo[:R_t], :)
        # DD[obs[Symbol("obs_nominalrate$i")]]    = m[:Rstarn] + ZZ_obs_nomrate * CCC_accum
        DD[obs[Symbol("obs_nominalrate$i")]]    = m[:Rstarn] + CCC_accum[endo[:R_t]]
        if subspec(m) == "ss11"
            QQ[exo[Symbol("rm_shl$i")], exo[Symbol("rm_shl$i")]] = m[:σ_r_m]^2 / n_mon_anticipated_shocks(m)
        else
            QQ[exo[Symbol("rm_shl$i")], exo[Symbol("rm_shl$i")]] = m[Symbol("σ_r_m$i")]^2
        end
    end

    # Anticipated GDP growth
    if haskey(get_settings(m), :add_anticipated_obs_gdp)
        if get_setting(m, :add_anticipated_obs_gdp)
            ZZ_obs_gdp = ZZ[obs[:obs_gdp], :]'
            meas_err = haskey(get_settings(m), :meas_err_anticipated_obs_gdp) ?
                get_setting(m, :meas_err_anticipated_obs_gdp) : 0. # Ignore measurement error for anticipated GDP growth
            ZZ_obs_gdp[endo_new[:e_gdp_t]]  = meas_err
            ZZ_obs_gdp[endo_new[:e_gdp_t1]] = -meas_err * m[:me_level]

            for i = 1:get_setting(m, :n_anticipated_obs_gdp)
                TTT_accum, CCC_accum = k_periods_ahead_expectations(TTT, CCC, TTTs, CCCs, reg, i, permanent_t;
                                                                    integ_series = integ_series,
                                                                    memo = (isnothing(memo) || !use_fwd_exp) ? nothing :
                                                                    ForwardExpectationsMemo(memo.time_varying_memo[min(reg + i, permanent_t, 17)],
                                                                                memo.permanent_memo))
                ZZ[obs[Symbol("obs_gdp$i")], :] = ZZ_obs_gdp * TTT_accum
                DD[obs[Symbol("obs_gdp$i")]]    = 100. * (exp(m[:z_star]) - 1.) + ZZ_obs_gdp * CCC_accum
                if haskey(get_settings(m), :add_iid_anticipated_obs_gdp_meas_err) ?
                    get_setting(m, :add_iid_anticipated_obs_gdp_meas_err) : false
                    ZZ[obs[Symbol("obs_gdp$i")], endo_new[:e_gdpexp_t]] = 1.
                end
            end
        end
    end

    # Adjustment to DD because measurement equation assumes CCC is the zero vector
    if any(CCC .!= 0.)
        # Are we using the zero rate alternative policy or not?
        is_zero_rate_rule = (haskey(get_settings(m), :regime_eqcond_info) ?
                             (!isempty(get_setting(m, :regime_eqcond_info)) && get_setting(m, :replace_eqcond)) : false) ||
        alternative_policy(m).eqcond == zero_rate_eqcond

        # If we are using the zero rate rule, then we don't adjust DD. If we run the code block below,
        # then the zero rate rule will cause CCC to be nonzero in multiple places, leading to the
        # wrong steady state measurement vector DD.
        #
        # The current coding of is_zero_rate_rule assumes that, when using a zero_rate rule (perhaps with gensys2 and temp alt policies)
        # any regimes in which the zero rate rule does not apply also will not have to update DD to reflect nonzero CCC.
        # This is typically the case for m1002.
        #
        # Currently it looks like we modify CCC in augment_states to track expected inflation,
        # but it doesn't look like we use this variable anywhere anyway. So ignoring this code block is fine on that front.
        # However, all the observables and pseudo-observables were not thoroughly checked, so it is possible
        # that the measurement equation is flawed in some way given how this has been coded.
        if !is_zero_rate_rule
            DD += ZZ * ((UniformScaling(1) - TTT) \ CCC)
        end
    end

    for para in m.parameters
        if !isempty(para.regimes)
            ModelConstructors.toggle_regime!(para, 1)
        end
    end

    return Measurement(ZZ, DD, QQ, EE)
end


text:
module PerformanceRNN

using Flux, MIDI, NoteSequences
using NoteSequences.PerformanceRepr

include("perfrnn.jl")
include("pretrain.jl")

@init register_configs(configs)

end

ext: jl

max_stars_repo_name: VasanthManiVasi/PerformanceRNN

max_issues_repo_name: VasanthManiVasi/PerformanceRNN

max_forks_repo_name: VasanthManiVasi/PerformanceRNN


text:
# # Quadrotor
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/Quadrotor.ipynb)
#
#md # !!! note "Overview"
#md #     System type: polynomial continuous system\
#md #     State dimension: 2\
#md #     Application domain: Chemical kinetics
#
# ## Model description
#
# We study the dynamics of a quadrotor as derived in [xxx]. Let us first introduce
# the variables required to describe the model: the inertial (north) position
# ``x_1``, the inertial (east) position ``x_2``, the altitude ``x_3``, the
# longitudinal velocity ``x_4``, the lateral velocity ``x_5``, the vertical
# velocity ``x_6``, the roll angle ``x_7``, the pitch angle ``x_8``, the yaw
# angle ``x_9``, the roll rate ``x_{10}``, the pitch rate ``x_{11}``, and the
# yaw rate ``x_{12}``. We further  require the following parameters: gravity
# constant ``g = 9.81`` [m/s``^2``], radius of center mass ``R = 0.1`` [m],
# distance of motors to center mass ``l = 0.5`` [m], motor mass
# ``M_{rotor} = 0.1`` [kg], center mass ``M = 1`` [kg], and total mass
# ``m = M + 4M_{rotor}``.
#
# From the above parameters we can compute the moments of inertia as
#
# ```math
# \begin{array}{lcll}
# J_x & = & \frac{2}{5}, M, R^2 + 2, l^2, M_{rotor}, \\
# J_y & = & J_x, \\
# J_z & = & \frac{2}{5}, M, R^2 + 4, l^2, M_{rotor}.
# \end{array}
# ```
#
# Finally, we can write the set of ordinary differential equations for the
# quadrotor according to [xxx]:
# ```math
# \left\{
# \begin{array}{lcl}
# \dot{x}_1 & = & \cos(x_8)\cos(x_9)x_4 + \Big(\sin(x_7)\sin(x_8)\cos(x_9) - \cos(x_7)\sin(x_9)\Big)x_5 \\
# & & + \Big(\cos(x_7)\sin(x_8)\cos(x_9) + \sin(x_7)\sin(x_9)\Big)x_6 \\
# \dot{x}_2 & = & \cos(x_8)\sin(x_9)x_4 + \Big(\sin(x_7)\sin(x_8)\sin(x_9) + \cos(x_7)\cos(x_9)\Big)x_5 \\
# & & + \Big(\cos(x_7)\sin(x_8)\sin(x_9) - \sin(x_7)\cos(x_9)\Big)x_6 \\
# \dot{x}_3 & = & \sin(x_8)x_4 - \sin(x_7)\cos(x_8)x_5 - \cos(x_7)\cos(x_8)x_6 \\
# \dot{x}_4 & = & x_{12}x_5 - x_{11}x_6 - g\sin(x_8) \\
# \dot{x}_5 & = & x_{10}x_6 - x_{12}x_4 + g\cos(x_8)\sin(x_7) \\
# \dot{x}_6 & = & x_{11}x_4 - x_{10}x_5 + g\cos(x_8)\cos(x_7) - \frac{F}{m} \\
# \dot{x}_7 & = & x_{10} + \sin(x_7)\tan(x_8)x_{11} + \cos(x_7)\tan(x_8)x_{12} \\
# \dot{x}_8 & = & \cos(x_7)x_{11} - \sin(x_7)x_{12} \\
# \dot{x}_9 & = & \frac{\sin(x_7)}{\cos(x_8)}x_{11} + \frac{\cos(x_7)}{\cos(x_8)}x_{12} \\
# \dot{x}_{10} & = & \frac{J_y - J_z}{J_x}x_{11}x_{12} + \frac{1}{J_x}\tau_\phi \\
# \dot{x}_{11} & = & \frac{J_z - J_x}{J_y}x_{10}x_{12} + \frac{1}{J_y}\tau_\theta \\
# \dot{x}_{12} & = & \frac{J_x - J_y}{J_z}x_{10}x_{11} + \frac{1}{J_z}\tau_\psi
# \end{array}
# \right.
# ```
#
# To check interesting control specifications, we stabilize the quadrotor using
# simple PD controllers for height, roll, and pitch. The inputs to the
# controller are the desired values for height, roll, and pitch ``u_1``, ``u_2``,
# and ``u_3``, respectively. The equations of the controllers are:
# ```math
# \begin{array}{lcll}
# F & = & m \, g - 10(x_3 - u_1) + 3x_6 \; & (\text{height control}), \\
# \tau_\phi & = & -(x_7 - u_2) - x_{10} & (\text{roll control}), \\
# \tau_\theta & = & -(x_8 - u_3) - x_{11} & (\text{pitch control}).
# \end{array}
# ```

# We leave the heading uncontrolled so that we set ``\tau_\psi = 0``.

using ReachabilityAnalysis
using ReachabilityAnalysis: is_intersection_empty

# ## Reachability settings
#
#The task is to change the height from ``0``~[m] to ``1``~[m] within ``5``~[s].
# A goal region ``[0.98,1.02]`` of the height ``x_3`` has to be reached within
# ``5``~[s] and the height has to stay below ``1.4`` for all times. After
# ``1``~[s] the height should stay above ``0.9``~[m]. The initial value for the
# position and velocities (i.e., from ``x_1`` to ``x_6``) is uncertain and given
# by ``[-\Delta,\Delta]``~[m], with ``\Delta=0.4``. All other variables are
# initialized to ``0``. This preliminary analysis must be followed by a
# corresponding evolution for ``\Delta = 0.1`` and ``\Delta= 0.8`` while keeping
# all the settings the same. No goals are specified for these cases: the
# objective instead is to understand the scalability of each tool with fixed
# settings.

const Tspan = (0.0, 5.0)

@taylorize function drone_10d!(dx, x, params, t)
    # dx[1] = 1.4*x[3] - 0.9*x[1]
    # dx[2] = 2.5*x[5] - 1.5*x[2]
    # dx[3] = 0.6*x[7] - 0.8*(x[2]*x[3])
    # dx[4] = 2 - 1.3*(x[3]*x[4])
    # dx[5] = 0.7*x[1] - (x[4]*x[5])
    # dx[6] = 0.3*x[1] - 3.1*x[6]
    # dx[7] = 1.8*x[6] - 1.6*(x[2]*x[7])
    # dx[0+1]=x[1+1];
    # dx[1+1]=9.81*tan(x[2+1]);
    # dx[2+1]=-8*x[2+1] + x[3+1];
    # dx[3+1]=-31.6227766*x[0+1] - 45.1663196*x[1+1] - 52.5011267*x[2+1] - 13.6015533*x[3+1] + 8.81028908e-17*x[4+1] - 4.23825517e-16*x[5+1] - 6.54099623e-16*x[6+1] + 1.59750882e-16*x[7+1] - 5.67884552e-16*x[8+1] - 1.48006779e-15*x[9+1] + 158.113883;
    # dx[4+1]=x[5+1];
    # dx[5+1]=9.81*tan(x[6+1]);
    # dx[6+1]=-8*x[6+1] + x[7+1];
    # dx[7+1]=2.05880071e-15*x[0+1] + 2.64242752e-15*x[1+1] + 2.62071672*x[2+1] + 1.59750882e-16*x[3+1] - 10.0*x[4+1] - 18.322492*x[5+1] - 24.6496768*x[6+1] - 11.3709874*x[7+1] + 1.9321227e-15*x[8+1] + 2.67619048e-15*x[9+1] - 1.029400355e-14;
    # dx[8+1]=x[9+1];
    # dx[9+1]=-5.0652413175e-17*x[0+1] - 6.7108633774e-17*x[1+1] - 1.194713975e-16*x[2+1] - 1.2256441379e-17*x[3+1] + 9.306933363e-17*x[4+1] + 1.2985133525e-16*x[5+1] + 1.0396072232e-16*x[6+1] + 2.2161533394e-17*x[7+1] - 0.91*x[8+1] - 1.6272983727*x[9+1];
    dx[0+1]=x[1+1];
    dx[1+1]=9.81*tan(x[2+1]);
    dx[2+1]=-8*x[2+1] + x[3+1];
    dx[3+1]=-45.1663196*x[1+1] - 52.5011267*x[2+1] - 13.6015533*x[3+1] + 8.81028908e-17*x[4+1]/sqrt((5 - x[0+1])*(5 - x[0+1]) + (x[4+1])*(x[4+1]) + (x[8+1])*(x[8+1])) - 4.23825517e-16*x[5+1] - 6.54099623e-16*x[6+1] + 1.59750882e-16*x[7+1] - 5.67884552e-16*x[8+1]/sqrt((5 - x[0+1])*(5 - x[0+1]) + (x[4+1])*(x[4+1]) + (x[8+1])*(x[8+1])) - 1.48006779e-15*x[9+1] + 31.6227766*(5 - x[0+1])/sqrt((5 - x[0+1])*(5 - x[0+1]) + (x[4+1])*(x[4+1]) + (x[8+1])*(x[8+1]));
    dx[4+1]=x[5+1];
    dx[5+1]=9.81*tan(x[6+1]);
    dx[6+1]=-8*x[6+1] + x[7+1];
    dx[7+1]=2.64242752e-15*x[1+1] + 2.62071672*x[2+1] + 1.59750882e-16*x[3+1] - 10.0*x[4+1]/sqrt((5 - x[0+1])*(5 - x[0+1]) + (x[4+1])*(x[4+1]) + (x[8+1])*(x[8+1])) - 18.322492*x[5+1] - 24.6496768*x[6+1] - 11.3709874*x[7+1] + 1.9321227e-15*x[8+1]/sqrt((5 - x[0+1])*(5 - x[0+1]) + (x[4+1])*(x[4+1]) + (x[8+1])*(x[8+1])) + 2.67619048e-15*x[9+1] - 2.05880071e-15*(5 - x[0+1])/sqrt((5 - x[0+1])*(5 - x[0+1]) + (x[4+1])*(x[4+1]) + (x[8+1])*(x[8+1]));
    dx[8+1]=x[9+1];
    dx[9+1]=-6.7108633774e-17*x[1+1] - 1.194713975e-16*x[2+1] - 1.2256441379e-17*x[3+1] + 9.306933363e-17*x[4+1]/sqrt((5 - x[0+1])*(5 - x[0+1]) + (x[4+1])*(x[4+1]) + (x[8+1])*(x[8+1])) + 1.2985133525e-16*x[5+1] + 1.0396072232e-16*x[6+1] + 2.2161533394e-17*x[7+1] - 0.91*x[8+1]/sqrt((5 - x[0+1])*(5 - x[0+1]) + (x[4+1])*(x[4+1]) + (x[8+1])*(x[8+1])) - 1.6272983727*x[9+1] + 5.0652413175e-17*(5 - x[0+1])/sqrt((5 - x[0+1])*(5 - x[0+1]) + (x[4+1])*(x[4+1]) + (x[8+1])*(x[8+1]));

    return dx
end

function drone_10d()
    ## initial states
    X0 = Hyperrectangle(low=[-1,0,0,0,-1,0,0,0,-1,0], high=[1,0,0,0,1,0,0,0,1,0])

    ## initial-value problem
    prob = @ivp(x' = drone_10d!(x), dim: 10, x(0) ∈ X0)

    return prob
end

# ----------------------------------------
#  Case 1: smaller uncertainty
# ----------------------------------------
Wpos = 0.1
Wvel = 0.1
prob = drone_10d()
alg = TMJets(abstol=1e-7, orderT=5, orderQ=1, adaptive=false)

# Warm-up run
sol1 = solve(prob, tspan=Tspan, alg=alg);
solz1 = overapproximate(sol1, Zonotope);

using Plots

# Plots.plot(solz1, vars=(0, 1), linecolor="yellow", color=:yellow, alpha=0.8,
#     xlab="t", ylab="x")
# Plots.plot(solz1, vars=(0, 5), linecolor="green", color=:green, alpha=0.8,
#     xlab="t", ylab="y")
Plots.plot(solz1, vars=(0, 9), linecolor="blue", color=:blue, alpha=0.8,
    xlab="t", ylab="z")


text:
include("NaiveCFs.jl")
include("LocalCFs.jl")
include("OneVisionCFs.jl")


text:
struct Acronym
    word::String
    components::Vector{Union{String, Nothing}}
    nvals::Vector{Int}
end

struct AcronymFamily
    word::String
    n_max::Int
    components::Matrix{Vector{String}}
end

Base.show(io::IO, af::AcronymFamily) = print(io, typeof(af), "(\"", af.word, "\", ", af.n_max, ", ...)")

function _family_traces(af::AcronymFamily)
    max_n, l = size(af.components)
    a = [Vector{Tuple{Int, Int}}[] for _ in af.components]

    for i in l:-1:1, n in 1:max_n
        n > 1 && isempty(af.components[n, i]) && continue  # No components

        ind = (i, n)
        next = i + n

        if next == l + 1  # Last
            push!(a[n, i], [ind])

        elseif next <= l
            for nexttraces in a[:, next]
                for nexttrace in nexttraces
                    push!(a[n, i], vcat(ind, nexttrace))
                end
            end
        end

    end

    return vcat(a[:, 1]...)
end

function acronyms(af::AcronymFamily, trace::AbstractVector{NTuple{2, Int}})
    comp_lists = [af.components[n, i] for (i, n) in trace]
    nvals = last.(trace)
    comps_iter = Iterators.product((isempty(c) ? [nothing] : c for c in comp_lists)...)
    return (Acronym(af.word, collect(comp), nvals) for comp in comps_iter)
end

acronyms(af::AcronymFamily) = Iterators.flatten(acronyms(af, trace) for trace in _family_traces(af))

function generate_acronyms(word::String, max_n::Int, prefix_map::Dict)
    word = lowercase(word)
    l = length(word)

    components = Matrix{Vector{String}}(undef, max_n, l)

    for i in 1:l
        for j in 1:max_n
            if i + j - 1 > l
                components[j, i] = String[]

            else
                prefix = word[i:(i + j - 1)]
                components[j, i] = get(prefix_map, prefix, String[])
            end
        end
    end

    return AcronymFamily(word, max_n, components)
end


text:
using Test, AxisKeys

# These used to be special methods of AxisKeys.findindex,
# now handled by Base.findall(in(Interval), ...)
# and tested more carefully here:
# https://github.com/JuliaMath/IntervalSets.jl/pull/63

if VERSION >= v"1.2" # <(3) doesn't exist on 1.1, but Base.Fix2 is fine
@testset "unit ranges" begin

    for r in (Base.OneTo(5), 2:5)
        x = 2:4
        @test r[AxisKeys.findindex(x, r)] == x
        @test @inferred(AxisKeys.findindex(x, r)) isa AbstractUnitRange

        for x in -2:7

            @test AxisKeys.findindex(x, r) == findfirst(==(x), collect(r))

            for op in (isequal, Base.:(==), Base.:<, Base.:<=, Base.:>, Base.:>=)

                @test AxisKeys.findindex(op(x), r) == Base.findall(op(x), collect(r))
                @test AxisKeys.findindex(op(x), r) isa AbstractRange
                T = typeof(@inferred AxisKeys.findindex(op(x), r))
                T <: AbstractRange || @info "findindex($op($x), $r) isa $T, unit ranges"
            end
        end
    end

end
@testset "step ranges" begin

    for op in (Base.:<=, Base.:>=, Base.:<, Base.:>)

        for r in (0:3:9, 0.5:0.5:8)
        for x in [-0.1, 0, 0.2, 0.5, 0.7, 1.0, 1.2,  1.9, 2.0, 2.1, 7.9, 8, 8.1, 7.9, 8, 8.1]

            @test AxisKeys.findindex(op(x), r) == Base.findall(op(x), collect(r))
            T = typeof(AxisKeys.findindex(op(x), r))
            T <: AbstractRange || @info "findindex($op($x), $r) isa $T, step ranges"
        end
        end

        for r in ('b':'e', 'b':2:'f')
        for x in 'a':'g'

            @test AxisKeys.findindex(op(x), r) == Base.findall(op(x), collect(r))
            T = typeof(@inferred AxisKeys.findindex(op(x), r))
            T <: AbstractRange || @info "findindex($op($x), $r) isa $T, step ranges of Char"
        end
        end

    end

end
end # VERSION


text:
# https://www.hpl.hp.com/techreports/2003/HPL-2003-148.pdf
function ematchlist(e::EGraph, t::Array{Union{PatTerm, PatVar}} , v::Array{Id}, sub)
    Channel() do c
        if length(t) == 0
            put!(c, sub)
        else
            for sub1 in ematch(e, t[1], v[1], sub)
                for sub2 in ematchlist(e, t[2:end], v[2:end], sub1)
                    put!(c, sub2)
                end
            end
        end
    end
end
# sub should be a map from pattern variables to Id
function ematch(e::EGraph, t::PatVar, v::Id, sub)
    Channel() do c
        if haskey(sub, t)
            if find_root!(e, sub[t]) == find_root!(e, v)
                put!(c, sub)
            else
                pass
            end
        else
            put!(c,  Base.ImmutableDict(sub, t => find_root!(e, v)))
        end
    end
end

function ematch(e::EGraph, t::PatTerm, v::Id, sub)
    Channel() do c
        for n in e.classes[find_root!(e,v)].nodes
            if n.head == t.head
                for sub1 in ematchlist(e, t.args , n.args , sub)
                    put!(c,sub1)
                end
            end
        end
    end
end

function instantiate(e::EGraph, p::PatVar , sub)
    sub[p]
end

function instantiate(e::EGraph, p::PatTerm , sub)
    push!( e, Term(p.head, [ instantiate(e,a,sub) for a in p.args ] ))
end

struct Rule
    lhs::Pattern
    rhs::Pattern
end

function rewrite!(e::EGraph, r::Rule)
    matches = []
    EMPTY_DICT2 = Base.ImmutableDict{PatVar, Id}(PatVar(:____),  Id(-1))
    for (n, cls) in e.classes
        for sub in ematch(e, r.lhs, n, EMPTY_DICT2)
            push!( matches, ( instantiate(e, r.lhs ,sub)  , instantiate(e, r.rhs ,sub)))
        end
    end
    for (l,r) in matches
        union!(e,l,r)
    end
    rebuild!(e)
end


text:
include("./src/solver.jl")
using PyPlot

#=
solver.jl:
- Julia version: 0.6.2
- Author: jerrymei
- Date: 2018-04-27
=#

k(x, y)=0.6 * sin(2 * pi * x) + 0.2 * sin(4 * pi * y) + 0.1 * sin(10 * pi * x + 8 * pi * y)
N_mesh=512

x=linspace(0,1 - 1 / N_mesh,N_mesh)
y=linspace(0,1 - 1 / N_mesh,N_mesh)
sol=zeros(N_mesh,N_mesh)
for i=1:N_mesh
    for j=1:N_mesh
        sol[i,j]=k(x[i],y[j])
    end
end

f=get_f(sol,1.)

u,err=solver(f,1.,1e-6)
figure(1)
plot_surface(linspace(0,1 - 1 / N_mesh,N_mesh),linspace(0,1 - 1 / N_mesh,N_mesh),u)
println("Norm of the distance between Du and f ", vecnorm(get_f(u,1.) - f)/vecnorm(f))

u,err=multigrid_solver(f,1.,1e-6)
figure(1)
plot_surface(linspace(0,1 - 1 / N_mesh,N_mesh),linspace(0,1 - 1 / N_mesh,N_mesh),u)
println("Norm of the distance between Du and f ", vecnorm(get_f(u,1.) - f)/vecnorm(f))

max_stars_repo_name: JerryLingjieMei/NonlinearPDE

max_issues_repo_name: JerryLingjieMei/NonlinearPDE

max_forks_repo_name: JerryLingjieMei/NonlinearPDE


text:
include("init.jl")

@testset "MessageFederate test_message_federate_initialize" begin
    broker = createBroker()
    mFed, fedinfo = createMessageFederate()

    state = h.helicsFederateGetState(mFed)
    @test state == 0
    h.helicsFederateEnterExecutingMode(mFed)

    state = h.helicsFederateGetState(mFed)
    @test state == 2

    destroyFederate(mFed, fedinfo)
    destroyBroker(broker)
end

@testset "MessageFederate test_message_federate_endpoint_registration" begin
    broker = createBroker()
    mFed, fedinfo = createMessageFederate()

    epid1 = h.helicsFederateRegisterEndpoint(mFed, "ep1", "")
    epid2 = h.helicsFederateRegisterGlobalEndpoint(mFed, "ep2", "random")

    h.helicsFederateEnterExecutingMode(mFed)

    @test h.HELICS_STATE_EXECUTION == h.helicsFederateGetState(mFed)

    endpoint_name = h.helicsEndpointGetName(epid1)
    @test endpoint_name == "TestA Federate/ep1"

    endpoint_name = h.helicsEndpointGetName(epid2)
    @test endpoint_name == "ep2"

    endpoint_name = h.helicsEndpointGetType(epid1)
    @test endpoint_name == ""

    endpoint_name = h.helicsEndpointGetType(epid2)
    @test endpoint_name == "random"

    epid_b = h.helicsFederateGetEndpoint(mFed, "ep2")
    type = h.helicsEndpointGetType(epid_b)
    @test type == "random"

    epid_c = h.helicsFederateGetEndpointByIndex(mFed, 0)
    name = h.helicsEndpointGetName(epid_c)
    @test name == "TestA Federate/ep1"

    destroyFederate(mFed, fedinfo)
    destroyBroker(broker)
end

@testset "MessageFederate test_message_federate_send" begin
    broker = createBroker()
    mFed, fedinfo = createMessageFederate()

    epid1 = h.helicsFederateRegisterEndpoint(mFed, "ep1", "")
    epid2 = h.helicsFederateRegisterGlobalEndpoint(mFed, "ep2", "random")

    h.helicsFederateSetTimeProperty(mFed, h.HELICS_PROPERTY_TIME_DELTA, 1.0)
    h.helicsFederateEnterExecutingMode(mFed)

    data = "random-data"

    h.helicsEndpointSendEventRaw(epid1, "ep2", data, 1.0)

    granted_time = h.helicsFederateRequestTime(mFed, 2.0)
    @test granted_time == 1.0

    res = h.helicsFederateHasMessage(mFed)
    @test res == true

    res = h.helicsEndpointHasMessage(epid1)
    @test res == false

    res = h.helicsEndpointHasMessage(epid2)
    @test res == true

    message = h.helicsEndpointGetMessageObject(epid2)

    @test h.helicsMessageGetMessageID(message) == 55
    @test h.helicsMessageIsValid(message) == true
    @test h.helicsMessageGetString(message) == "random-data"
    @test h.helicsMessageGetRawDataSize(message) == 11
    @test h.helicsMessageGetOriginalDestination(message) == ""
    @test h.helicsMessageGetOriginalSource(message) == "TestA Federate/ep1"
    @test h.helicsMessageGetSource(message) == "TestA Federate/ep1"
    @test h.helicsMessageGetTime(message) == 1.0
    @test_broken false
    # h.helicsMessageGetRawData(message) crashes

    destroyFederate(mFed, fedinfo)
    destroyBroker(broker)
end

@testset "MessageFederate send_receive_2fed_multisend" begin

    broker = createBroker(2)
    mFed1, fedinfo1 = createMessageFederate(1, "A Federate")
    mFed2, fedinfo2 = createMessageFederate(1, "B Federate")

    epid1 = h.helicsFederateRegisterEndpoint(mFed1, "ep1", "")
    epid2 = h.helicsFederateRegisterGlobalEndpoint(mFed2, "ep2", "random")

    h.helicsEndpointSetOption(epid1, h.HELICS_HANDLE_OPTION_IGNORE_INTERRUPTS, true)

    h.helicsFederateSetTimeProperty(mFed1, h.HELICS_PROPERTY_TIME_DELTA, 1.0)
    h.helicsFederateSetTimeProperty(mFed2, h.HELICS_PROPERTY_TIME_DELTA, 1.0)

    h.helicsFederateEnterExecutingModeAsync(mFed1)
    h.helicsFederateEnterExecutingMode(mFed2)
    h.helicsFederateEnterExecutingModeComplete(mFed1)

    @test h.HELICS_STATE_EXECUTION == h.helicsFederateGetState(mFed1)
    @test h.HELICS_STATE_EXECUTION == h.helicsFederateGetState(mFed2)

    h.helicsEndpointSetDefaultDestination(epid1, "ep2");

    h.helicsEndpointSendMessageRaw(epid1, "", "a")
    h.helicsEndpointSendMessageRaw(epid1, "", "a")
    h.helicsEndpointSendMessageRaw(epid1, "", "a")

    h.helicsFederateRequestTimeAsync(mFed1, 1.0)
    granted_time = h.helicsFederateRequestTime(mFed2, 1.0)
    complete_time = h.helicsFederateRequestTimeComplete(mFed1)

    @test granted_time == 1.0
    @test complete_time == 1.0

    res = h.helicsEndpointPendingMessages(epid2)
    @test res == 3

    res = h.helicsFederatePendingMessages(mFed2)
    @test res == 3

    @test h.helicsEndpointGetDefaultDestination(epid1) == "ep2"

    # FIXME: Someday this will be implemented.
    @test_broken h.helicsEndpointGetOption(epid1, h.HELICS_HANDLE_OPTION_IGNORE_INTERRUPTS) == true

    destroyFederate(mFed1, fedinfo1)
    destroyFederate(mFed2, fedinfo2)
    destroyBroker(broker)

end

@testset "MessageFederate message_object_tests" begin

    broker = createBroker(1)
    mFed1, fedinfo1 = createMessageFederate(1, "test")

    epid1 = h.helicsFederateRegisterEndpoint(mFed1, "ep1", "")
    epid2 = h.helicsFederateRegisterGlobalEndpoint(mFed1, "ep2", "random");

    h.helicsFederateSetTimeProperty(mFed1, h.HELICS_PROPERTY_TIME_DELTA, 1.0)

    h.helicsFederateEnterExecutingMode(mFed1)

    @test h.helicsFederateGetState(mFed1) == h.HELICS_STATE_EXECUTION

    msg = h.helicsFederateCreateMessageObject(mFed1)
    h.helicsMessageSetDestination(msg, "ep2")
    h.helicsMessageGetDestination(msg) == "ep2"
    h.helicsMessageSetData(msg, repeat("a", 500))
    h.helicsMessageSetTime(msg, 0.0)

    h.helicsEndpointSendMessageObject(epid1, msg)
    time = h.helicsFederateRequestTime(mFed1, 1.0)
    @test time == 1.0

    @test h.helicsFederateHasMessage(mFed1) == true
    @test h.helicsEndpointHasMessage(epid1) == false
    @test h.helicsEndpointHasMessage(epid2) == true

    msg = h.helicsEndpointGetMessageObject(epid2)
    @test h.helicsMessageGetRawDataSize(msg) == 500
    # @show h.helicsMessageGetRawData(msg)
    @test_broken false
    # segfaults
    rdata = h.helicsMessageGetRawDataPointer(msg)
    @test Char(unsafe_load(Ptr{Cchar}(rdata), 245)) == 'a'

    h.helicsFederateFinalize(mFed1)

    @test h.helicsFederateGetState(mFed1) == h.HELICS_STATE_FINALIZE

    h.helicsMessageSetFlagOption(msg, 7, true);
    @test h.helicsMessageCheckFlag(msg, 7) == true
    h.helicsMessageClearFlags(msg);
    @test h.helicsMessageCheckFlag(msg, 7) == false

    h.helicsEndpointSetDefaultDestination(epid1, "ep2")

    destroyFederate(mFed1, fedinfo1)
    destroyBroker(broker)

end

@testset "MessageFederate timing tests" begin

    broker = createBroker(1)
    vFed1, fedinfo1 = createMessageFederate(1, "fed0")
    vFed2, fedinfo2 = createMessageFederate(1, "fed1")

    h.helicsFederateSetTimeProperty(vFed1, h.HELICS_PROPERTY_TIME_PERIOD, 0.1)
    h.helicsFederateSetTimeProperty(vFed2, h.HELICS_PROPERTY_TIME_PERIOD, 0.1)

    h.helicsFederateSetTimeProperty(vFed2, h.HELICS_PROPERTY_TIME_INPUT_DELAY, 0.1)

    h.helicsFederateSetFlagOption(vFed1, h.HELICS_FLAG_IGNORE_TIME_MISMATCH_WARNINGS, true)
    h.helicsFederateSetFlagOption(vFed2, h.HELICS_FLAG_IGNORE_TIME_MISMATCH_WARNINGS, true)

    ept1 = h.helicsFederateRegisterGlobalEndpoint(vFed1, "e1", "")
    h.helicsFederateRegisterGlobalEndpoint(vFed2, "e2", "")

    h.helicsFederateEnterExecutingModeAsync(vFed1)
    h.helicsFederateEnterExecutingMode(vFed2)
    h.helicsFederateEnterExecutingModeComplete(vFed1)
    h.helicsFederateRequestTimeAsync(vFed2, 2.0)
    gtime = h.helicsFederateRequestTime(vFed1, 1.0)
    # check that the request is only granted at the appropriate period
    @test gtime == 1.0

    @test h.helicsFederateGetIntegerProperty(vFed1, h.HELICS_PROPERTY_INT_CONSOLE_LOG_LEVEL) == -1
    @test h.helicsFederateGetIntegerProperty(vFed2, h.HELICS_PROPERTY_INT_CONSOLE_LOG_LEVEL) == -1

    @test h.helicsFederateGetFlagOption(vFed1, h.HELICS_FLAG_IGNORE_TIME_MISMATCH_WARNINGS) == true
    @test h.helicsFederateGetFlagOption(vFed2, h.HELICS_FLAG_IGNORE_TIME_MISMATCH_WARNINGS) == true

    h.helicsEndpointSendMessageRaw(ept1, "e2", "test1")
    h.helicsFederateRequestTimeAsync(vFed1, 1.9)
    gtime = h.helicsFederateRequestTimeComplete(vFed2)
    @test gtime >= 1.1 # the message should show up at the next available time point after the impact window
    h.helicsFederateRequestTimeAsync(vFed2, 2.0)
    gtime = h.helicsFederateRequestTimeComplete(vFed1)
    @test gtime >= 1.9

    tres = h.helicsFederateGetTimeProperty(vFed1, h.HELICS_PROPERTY_TIME_PERIOD)
    @test tres == 0.1

    # m = h.helicsEndpointGetMessageObject(ept1)
    # @show h.helicsMessageGetRawData(m)
    # TODO: null pointer received from C
    @test_broken false

    gtime = h.helicsFederateRequestTimeComplete(vFed2)
    @test gtime == 2.0
    h.helicsFederateFinalize(vFed1)
    h.helicsFederateFinalize(vFed2)

    destroyFederate(vFed1, fedinfo1)
    destroyFederate(vFed2, fedinfo2)
    destroyBroker(broker)

end


text:
before = """

foo

Note

This is a note.

bar

    """

after = """

foo

        Note

This is a note.

bar

    """

@test PlutoStaticHTML._convert_admonitions(before) == after

nb = Notebook([
    Cell("""
        md\"\"\"
        !!! note
            This is a note.
        \"\"\"
        """)
])
use_distributed = false
html, _ = notebook2html_helper(nb; use_distributed)
# This tests that there has been a hit on the `_convert_admonition` replacer.
@test contains(html, "admonition-header")


text:
"
`StopCondition` represents the set of stop conditions for the experiment to finish
if at least one of them has been fulfilled.

# Fields
- [`_maxEvaluations::Int64`]: maximum of evaluations permitted. An evaluation
    is made everytime the system computes an individual's fitness.
- [`_maxIterations::Int64`]: maximum of iterations permitted. An iteration is
    done in the main loop of the genetic algorithm, that is to say, everytime a
    new population replaces the old one.
- [`_maxIterNotImproving::Int64`]: maximum of iterations having the same best
    individual of the population.
- [`_maxTime::Float64`]: maximum run time permitted in seconds.
- [`_numEvaluations::Int64`]: current number of evaluations.
- [`_numIterations::Int64`]: current number of iterations.
- [`_numIterNotImproving::Int64`]: current number of iterations having the same
    best individual of the population.
- [`_bestFitness::Float64`]: fitness value of the best individual so far.
- [`_initialTime::Float64`]: starting time of the experiment in seconds.
"
mutable struct StopCondition
    _maxEvaluations::Int64
    _maxIterations::Int64
    _maxIterNotImproving::Int64
    _maxTime::Float64
    _numEvaluations::Int64
    _numIterations::Int64
    _numIterNotImproving::Int64
    _bestFitness::Float64
    _initialTime::Float64

    StopCondition() = new()
    function StopCondition(maxEvaluations, maxIterations,
                            maxIterNotImproving, maxTime,
                            numEvaluations, numIterations,
                            numIterNotImproving, initialTime)

        new(maxEvaluations, maxIterations,
            maxIterNotImproving, maxTime,
            numEvaluations, numIterations,
            numIterNotImproving, NaN, initialTime)
    end
end # struct

"""
    getMaxIterations(stopCondition::StopCondition)

Obtains the maximum of iterations permitted.
"""
getMaxIterations(stopCondition::StopCondition) = stopCondition._maxIterations
# function

"""
    getCurrentIteration(stopCondition::StopCondition)

Obtains the current iteration of the experiment.
"""
getCurrentIteration(stopCondition::StopCondition) = stopCondition._numIterations
# function

"""
    reached_(stopCondition::StopCondition)::Bool

Tells wether one of the stop conditions has been fulfilled or not.
"""
function reached_(stopCondition::StopCondition)::Bool

    reached = false

    if stopCondition._numEvaluations == stopCondition._maxEvaluations
        reached = true
    elseif stopCondition._numIterations == stopCondition._maxIterations
        reached = true
    elseif stopCondition._numIterNotImproving == stopCondition._maxIterNotImproving
        reached = true
    elseif (time() - stopCondition._initialTime) > stopCondition._maxTime
        reached = true
    end

    return reached
end # function

"""
    setTime(stopCondition::StopCondition, time::AbstractFloat)::Bool

Sets the starting time of the experiment.
"""
function setTime(stopCondition::StopCondition, time::AbstractFloat)
    stopCondition._initialTime = time
    nothing
end # function

"""
    notifyEvaluation(stopCondtition::StopCondition)

Adds one to the current number of evaluations.
"""
function notifyEvaluation(stopCondition::StopCondition)
    stopCondition._numEvaluations += 1
    nothing
end # function

"""
    notifyIteration(stopCondition::StopCondition)

Adds one to the current number of iterations.
"""
function notifyIteration(stopCondition::StopCondition)
    stopCondition._numIterations += 1
    nothing
end # function

"""
    notifyIterNotImproving(stopCondition::StopCondition)

Adds one to the current number of iterations having the same best individual
of the population.
"""
function notifyIterNotImproving(stopCondition::StopCondition)
    stopCondition._numIterNotImproving += 1
    nothing
end # function

"""
    resetIterNotImproving(stopCondition::StopCondition)

Resets the number of iterations having the same best individual of the population
back to 0.
"""
function resetIterNotImproving(stopCondition::StopCondition)
    stopCondition._numIterNotImproving = 0
    nothing
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
import ..UncertainDatasets.UncertainIndexDataset
import ..UncertainDatasets.AbstractUncertainIndexDataset

##############################
# Resample the entire dataset
##############################

"""
	resample(uvd::AbstractUncertainIndexDataset)

Draw a realisation of an uncertain index dataset according to the distributions
of the uncertain values comprising it.
"""
function resample(uvd::DT) where {DT <: AbstractUncertainIndexDataset}
	L = length(uvd)
	[resample(uvd.indices[i]) for i in 1:L]
end

"""
	resample(uv::AbstractUncertainIndexDataset, constraint::Vector{<:SamplingConstraint})

Draw a realisation of an uncertain index dataset according to the uncertain values
comprising it, while constraining the distributions furnishing the values in the dataset
according the provided sampling constraints.
"""
function resample(uv::DT, constraint::Vector{<:SamplingConstraint}) where {
		DT <: AbstractUncertainIndexDataset}
    [resample(uv.indices[i], constraint[i]) for i in 1:length(uv)]
end

function resample(uv::DT, constraint::SamplingConstraint) where {
		DT <: AbstractUncertainIndexDataset}
	[resample(uv.indices[i], constraint) for i in 1:length(uv)]
end

function resample(uv::DT, constraint::SamplingConstraint, n::Int) where {
		DT <: AbstractUncertainIndexDataset}
	[[resample(uv.indices[i], constraint) for i in 1:length(uv)] for k = 1:n]
end

function resample(uv::DT, constraint::Vector{<:SamplingConstraint}, n::Int) where {
		DT <: AbstractUncertainIndexDataset}
	[[resample(uv.indices[i], constraint[i]) for i in 1:length(uv)] for k = 1:n]
end

##############################
# Resample the entire dataset
##############################

"""
	resample_elwise(uvd::UncertainIndexDataset, n::Int)

Resample each element in `uvals` `n` times. The i-th entry in the returned
vector is a `n`-element vector consisting of `n` unique draws of `uvals[i]`.
"""
function resample_elwise(uvd::DT, n::Int) where {DT <: AbstractUncertainIndexDataset}
    [resample(uvd[i], n) for i = 1:length(uvd)]
end

"""
	resample(uvd::UncertainIndexDataset, n::Int)

Draw `n` realisations of an `UncertainIndexDataset` according to the
distributions of the `UncertainValue`s comprising it.
"""
function resample(uvd::DT, n::Int) where {DT <: AbstractUncertainIndexDataset}
	L = length(uvd)
	[[resample(uvd.indices[i]) for i in 1:L] for k in 1:n]
end

# ##########################################################################
# # Draw realisations of the uncertain dataset according to the distributions
# # of the uncertain values comprising it, according to sampling constraints
# ##########################################################################

# """
# 	resample(uv::UncertainIndexDataset, constraint::NoConstraint) -> Vector{Float64}

# Draw a realisation of an `UncertainIndexDataset` where each uncertain value
# is truncated according to `constraint`. In the case of `NoConstraint`,
# no trucation is performed and the whole distribution is sampled.
# """
# function resample(uv::DT, constraint::NoConstraint) where {
# 		DT <: AbstractUncertainIndexDataset}
# 	L = length(uv)
# 	[resample(uv.indices[i], constraint) for i in 1:L]
# end

# """
# 	resample(uv::UncertainIndexDataset, n::NoConstraint,
# 		n::Int) -> Vector{Vector{Float64}}

# Draw `n` realisation of an `UncertainIndexDataset` where each uncertain value
# is truncated according to `constraint`. In the case of `NoConstraint`,
# no trucation is performed and the whole distribution is sampled.
# """
# function resample(uv::DT, constraint::NoConstraint, n::Int) where {
# 		DT <: AbstractUncertainIndexDataset}
# 	L = length(uv)
# 	[[resample(uv.indices[i], constraint) for i in 1:L] for k = 1:n]
# end

# """
# 	resample(uv::UncertainIndexDataset,
# 		constraint::TruncateLowerQuantile) -> Vector{Float64}

# Draw a realisation of an `UncertainIndexDataset` where each uncertain value
# is truncated according to `constraint`. In the case of `TruncateLowerQuantile`,
# the supports of the distributions are truncated below at some quantile.
# """
# function resample(uv::DT, constraint::TruncateLowerQuantile) where {
# 		DT <: AbstractUncertainIndexDataset}
# 	L = length(uv)
# 	[resample(uv.indices[i], constraint) for i in 1:L]
# end

# """
# 	resample(uv::UncertainIndexDataset, constraint::TruncateLowerQuantile,
# 		n::Int) -> Vector{Vector{Float64}}

# Draw `n` realisation of an `UncertainIndexDataset` where each uncertain value
# is truncated according to `constraint`. In the case of `TruncateLowerQuantile`,
# the supports of the distributions are truncated below at some quantile.
# """
# function resample(uv::DT, constraint::TruncateLowerQuantile, n::Int) where {
# 		DT <: AbstractUncertainIndexDataset}
# 	L = length(uv)
# 	[[resample(uv.indices[i], constraint) for i in 1:L] for k = 1:n]
# end

# """
# 	resample(uv::UncertainIndexDataset,
# 		constraint::TruncateUpperQuantile) -> Vector{Float64}

# Draw a realisation of an `UncertainIndexDataset` where each uncertain value
# is truncated according to `constraint`. In the case of `TruncateLowerQuantile`,
# the supports of the distributions are truncated above at some quantile.
# """
# function resample(uv::DT, constraint::TruncateUpperQuantile) where {
# 		DT <: AbstractUncertainIndexDataset}
# 	L = length(uv)
# 	[resample(uv.indices[i], constraint) for i in 1:L]
# end

# """
# 	resample(uv::UncertainIndexDataset, constraint::TruncateUpperQuantile,
# 		n::Int) -> Vector{Vector{Float64}}

# Draw `n` realisation of an `UncertainIndexDataset` where each uncertain value
# is truncated according to `constraint`. In the case of `TruncateUpperQuantile`,
# the supports of the distributions are truncated above at some quantile.
# """
# function resample(uv::DT, constraint::TruncateUpperQuantile, n::Int) where {
# 		DT <: AbstractUncertainIndexDataset}
# 	L = length(uv)
# 	[[resample(uv.indices[i], constraint) for i in 1:L] for k = 1:n]
# end

# """
# 	resample(uv::UncertainIndexDataset,
# 		constraint::TruncateQuantiles) -> Vector{Float64}

# Draw a realisation of an `UncertainIndexDataset` where each uncertain value
# is truncated according to `constraint`. In the case of `TruncateLowerQuantile`,
# the supports of the distributions are truncated at some quantile range.
# """
# function resample(uv::UncertainIndexDataset, constraint::TruncateQuantiles) where {
# 		DT <: AbstractUncertainIndexDataset}
# 	L = length(uv)
# 	[resample(uv.indices[i], constraint) for i in 1:L]
# end

# """
# 	resample(uv::UncertainIndexDataset, n::Int,
# 		constraint::TruncateUpperQuantile) -> Vector{Vector{Float64}}

# Draw `n` realisation of an `UncertainIndexDataset` where each uncertain value
# is truncated according to `constraint`. In the case of `TruncateUpperQuantile`,
# the supports of the distributions are truncated above at some quantile.
# """
# function resample(uv::DT, constraint::TruncateQuantiles, n::Int) where {
# 		DT <: AbstractUncertainIndexDataset}
# 	L = length(uv)
# 	[[resample(uv.indices[i], constraint) for i in 1:L] for k = 1:n]
# end

# """
#     resample_elwise(uvd::UncertainIndexDataset, constraint::SamplingConstraint, n::Int)

# Resample each element in `uvals` `n` times. The i-th entry in the returned
# vector is a `n`-element vector consisting of `n` unique draws of `uvals[i]`, drawn
# after first truncating the support of `uvals[i]` according to the provided `constraint`.
# """
# function resample_elwise(uvd::DT, constraint::SamplingConstraint, n::Int) where {
# 		DT <: AbstractUncertainIndexDataset}
#     [resample(uvd[i], constraint, n) for i = 1:length(uvd)]
# end


text:
using YaoBase: @interface

parse_block(n::Int, x::Function) = x(n)

function parse_block(n::Int, x::AbstractBlock{N}) where N
    n == N || throw(ArgumentError("number of qubits does not match: $x"))
    return x
end

parse_block(n::Int, x::Pair{Int, <:AbstractBlock{N}}) where N = x

"""
    prewalk(f, src::AbstractBlock)

Walk the tree and call `f` once the node is visited.
"""
@interface function prewalk(f::Base.Callable, src::AbstractBlock)
    out = f(src)
    for each in subblocks(src)
        prewalk(f, each)
    end
    return out
end

"""
    postwalk(f, src::AbstractBlock)

Walk the tree and call `f` after the children are visited.
"""
@interface function postwalk(f::Base.Callable, src::AbstractBlock)
    for each in subblocks(src)
        postwalk(f, each)
    end
    return f(src)
end

@interface blockfilter!(f, v::Vector, blk::AbstractBlock) =
    postwalk(x -> f(x) ? push!(v, x) : v, blk)

@interface blockfilter(f, blk) = blockfilter!(f, [], blk)

"""
    collect_blocks(block_type, root)

Return a [`ChainBlock`](@ref) with all block of `block_type` in root.
"""
@interface collect_blocks(::Type{T}, x::AbstractBlock) where T <: AbstractBlock = blockfilter!(x->x isa T, T[], x)

#@interface expect(op::AbstractBlock, r::AbstractRegister) = r' * apply!(copy(r), op)

#expect(op::AbstractBlock, dm::DensityMatrix) = mapslices(x->sum(mat(op).*x)[], dm.state, dims=[1,2]) |> vec
expect(op::AbstractBlock, dm::DensityMatrix{1}) = sum(mat(op).*dropdims(dm.state, dims=3))

"""
    expect(op::AbstractBlock, reg::AbstractRegister{B}) -> Vector
    expect(op::AbstractBlock, dm::DensityMatrix{B}) -> Vector

expectation value of an operator.
"""
@interface function expect(op::AbstractBlock, dm::DensityMatrix{B}) where B
    mop = mat(op)
    [tr(view(dm.state,:,:,i)*mop) for i=1:B]
end

expect(op::AbstractBlock, reg::AbstractRegister{1}) = reg'*apply!(copy(reg), op)

function expect(op::AbstractBlock, reg::AbstractRegister{B}) where B
    ket = apply!(copy(reg), op)
    C = conj!(reshape(ket.state, :, B))
    A = reshape(reg.state, :, B)
    dropdims(sum(A.*C, dims=1), dims=1) |> conj
end

function expect(op::Sum, reg::AbstractRegister)
    sum(opi->expect(opi, reg), op)
end

function expect(op::Scale, reg::AbstractRegister)
    factor(op)*expect(content(op), reg)
end

expect(op::Sum, reg::AbstractRegister{1}) = invoke(expect, Tuple{Sum, AbstractRegister}, op, reg)
expect(op::Scale, reg::AbstractRegister{1}) = invoke(expect, Tuple{Scale, AbstractRegister}, op, reg)

for FUNC in [:measure!, :measure_collapseto!, :measure_remove!, :measure]
    @eval function YaoBase.$FUNC(op::AbstractBlock, reg::AbstractRegister, locs::AllLocs; kwargs...) where B
        $FUNC(eigen!(mat(op) |> Matrix), reg, locs; kwargs...)
    end
end


text:
@testset "rand" begin
# NOTE: tests should cover both pow2 and non-pow2 dims

# for rand_logn! and rand_poisson!
using AMDGPU.rocRAND

# in-place
for (f,T) in ((rand!, Float16),
              (rand!,Float32),
              (rand!,Cuint),
              (randn!,Float16),
              (randn!,Float32),
              (rand_logn!,Float32),
              (rand_poisson!,Cuint)),
    d in (2, (2,2), (2,2,2), 3, (3,3), (3,3,3))
    A = ROCArray{T}(undef, d)
    fill!(A, T(0))
    f(A)
    if f !== rand_poisson!
        @test !iszero(collect(A))
    end
end

# out-of-place, with implicit type
for (f,T) in ((AMDGPU.rand,Float32), (AMDGPU.randn,Float32),
              (AMDGPU.rand_logn,Float32), (AMDGPU.rand_poisson,Cuint),
              (rand,Float64), (randn,Float64)),
    args in ((2,), (2, 2), (3,), (3, 3))
    A = f(args...)
    @test eltype(A) == T
end

# out-of-place, with type specified
for (f,T) in ((AMDGPU.rand,Float32), (AMDGPU.randn,Float32), (AMDGPU.rand_logn,Float32),
              (AMDGPU.rand,Float64), (AMDGPU.randn,Float64), (AMDGPU.rand_logn,Float64),
              (AMDGPU.rand_poisson,Cuint),
              (rand,Float32), (randn,Float32),
              (rand,Float64), (randn,Float64)),
    args in ((T, 2), (T, 2, 2), (T, (2, 2)), (T, 3), (T, 3, 3), (T, (3, 3)))
    A = f(args...)
    @test eltype(A) == T
end

# FIXME: Int64 support
# GPUArrays rng
# unsupported types that always fall back to GPUArrays
for (f,T) in ((rand!,Int32),),
    d in (2, (2,2), (2,2,2), 3, (3,3), (3,3,3))
    A = ROCArray{T}(undef, d)
    fill!(A, T(0))
    f(AMDGPU.gpuarrays_rng(), A)
    @test !iszero(collect(A))
end
for (f,T) in ((AMDGPU.rand,Int32),),
    args in ((T, 2), (T, 2, 2), (T, (2, 2)), (T, 3), (T, 3, 3), (T, (3, 3)))
    A = f(args...)
    @test eltype(A) == T
end

@test_throws ErrorException randn!(ROCArray{Cuint}(undef, 10))
@test_throws ErrorException rand_logn!(ROCArray{Cuint}(undef, 10))
@test_throws ErrorException rand_poisson!(ROCArray{Float64}(undef, 10))

# seeding of both generators
AMDGPU.seed!()
AMDGPU.seed!(1)
## AMDGPU rocrand
AMDGPU.seed!(1)
a = AMDGPU.rand(Float32, 1)
AMDGPU.seed!(1)
b = AMDGPU.rand(Float32, 1)
# fixme: uncomment once mapreduce works
@test iszero(collect(a) - collect(b))
#@test all(a .== b)
## gpuarrays fallback
AMDGPU.seed!(1)
a = AMDGPU.rand(Int32, 1)
AMDGPU.seed!(1)
b = AMDGPU.rand(Int32, 1)
@test iszero(collect(a) - collect(b))
#@test all(a .== b)
end # testset


text:
module LowRankVortex

using Interpolations
using JLD
using LinearAlgebra
using PotentialFlow
import PotentialFlow.Motions: reset_velocity!
import PotentialFlow.Points
import PotentialFlow.Blobs
import PotentialFlow.Elements
using ProgressMeter
using Statistics
using TransportBasedInference
using NamedColors

include("vortex.jl")
include("convective_complexpotential.jl")
include("pressure.jl")
include("AD_pressure.jl")
include("analytical_jacobian_pressure.jl")
include("symmetric_pressure.jl")
include("symmetric_analytical_jacobian_pressure.jl")
include("symmetric_analytical_jacobian_pressure_freestream.jl")
include("forecast.jl")
include("generate_vortex.jl")
include("generate_patch.jl")
include("assimilation.jl")
include("lowrankassimilation.jl")
include("symmetric_assimilation.jl")
include("symmetric_lowrankassimilation.jl")
include("adaptive_symmetric_lowrankassimilation.jl")
include("plot_recipes.jl")

end # module


text:
using CTMCEnsemble
using Base.Test
using StatsBase

"""
    bench(sparsity=3, num=10000)

Find a good `maxiter` for power method.
`sparsity` controls generated distribution.
`num` is repeat number.
"""
function bench(sparsity=3, num=10000)
    newv() = normalize!(exp.(sparsity * randn(300)), 1)
    println("> sparsity=", sparsity, ", repeat=", num)
    stats = []
    for i = 1:num
        d = 600
        preds = [
            (newv(),     1:300),
            (newv(),   201:500),
            (newv(), ((401:700) - 1) % 600 + 1),
        ]
        i % 1000 == 0 && print(".")
        for maxiter = 0:8:32
            v_powermethod = powermethod(preds, maxiter=maxiter)
            v_ctmc        = ctmc(preds)
            diff = v_powermethod - v_ctmc
            same = indmax(v_powermethod) == indmax(v_ctmc)
            same || push!(stats, maxiter)
            # println(same ? "[O]" : "[X]"),
            #         "iter: ", maxiter,
            #         ", $d distinf = ", d * norm(diff, Inf),
            #         ", $d dist2 = ", d * norm(diff))
        end

    end
    m = collect(countmap(stats))
    print()
    println(join(join.(m[sortperm(getindex.(m, 1))], ["=>"]), "\n"))
    m[sortperm(getindex.(m, 1))]
end

bench(1)
bench(2)
bench(3)
bench(4)
bench(5)

# > sparsity=1, repeat=10000
# ..........0=>820
# 8=>146
# 16=>30
# 24=>7
# 32=>3
# > sparsity=2, repeat=10000
# ..........0=>1486
# 8=>240
# 16=>52
# 24=>7
# > sparsity=3, repeat=10000
# ..........0=>2180
# 8=>320
# 16=>86
# 24=>30
# 32=>10
# > sparsity=4, repeat=10000
# ..........0=>2465
# 8=>349
# 16=>115
# 24=>57
# 32=>33
# > sparsity=5, repeat=10000
# ..........0=>2567
# 8=>431
# 16=>187
# 24=>107
# 32=>82

println("love bufan")


text:
"""
    GDNetStats

Network statistics for `GDConfig`.
"""
mutable struct GDNetStats <: Statistics
    # gdfile stats
    gdfile_length::Int
    basis_size::Int
    deviation_size::Int
    hash_size::Int

    # network traffic
    uplink::Vector{Int}
    downlink::Vector{Int}

    num_requested_bases::Int
    num_unknown_bases::Int

    num_identical_hashes::Int
    num_identical_devs::Int

    GDNetStats(
        gdfile_length,
        basis_size,
        hash_size,
        deviation_size
    ) = new(
        gdfile_length,
        basis_size,
        hash_size,
        deviation_size,
        Vector{Int}(undef, 0),
        Vector{Int}(undef, 0),
        0,
        0,
        0,
        0,
    )
end

"""
    update_stats!(stats)

Updates the network traffic stats of a `GDConfig`.
"""
function update_stats!(stats::GDNetStats, round_num::Int)
    num_comm = round_num * STATS.base.num_clients_per_round
    gdfile_size = (stats.hash_size + stats.deviation_size) * stats.gdfile_length
    minimal_up_down = gdfile_size * num_comm

    overhead_down = stats.num_unknown_bases * stats.basis_size + stats.num_requested_bases * stats.hash_size
    overhead_up = stats.num_requested_bases * stats.basis_size + stats.num_unknown_bases * stats.hash_size

    gain_hash = (stats.hash_size-1) * stats.num_identical_hashes
    gain_dev = (stats.deviation_size-1) * stats.num_identical_devs
    gain_down = gain_hash + gain_dev

    downlink = minimal_up_down + overhead_down - gain_down
    uplink = minimal_up_down + overhead_up

    push!(stats.uplink, uplink)
    push!(stats.downlink, downlink)
end


text:
"""
		function layerbuilder(d::Int,k::Int,o::Int,n::Int,ftype::String,lastlayer::String = "",ltype::String = "Dense")
		create a chain with `n` layers of with `k` neurons with transfer function `ftype`.
		input and output dimension is `d` / `o`
		If lastlayer is no specified, all layers use the same function.
		If lastlayer is "linear", then the last layer is forced to be Dense.
		It is also possible to specify dimensions in a vector.
```juliadoctest
julia> FluxExtensions.layerbuilder(4,11,1,3,"relu")
Chain(Dense(4, 11, NNlib.relu), Dense(11, 11, NNlib.relu), Dense(11, 1, NNlib.relu))
julia> FluxExtensions.layerbuilder([4,11,11,1],"relu")
Chain(Dense(4, 11, NNlib.relu), Dense(11, 11, NNlib.relu), Dense(11, 1, NNlib.relu))
julia> FluxExtensions.layerbuilder(4,11,1,3,"relu","tanh")
Chain(Dense(4, 11, NNlib.relu), Dense(11, 11, NNlib.relu), Dense(11, 1, tanh))
# TODO fix this
julia> FluxExtensions.layerbuilder(4,11,1,3,"relu","tanh","ResDense")
Chain(ResDense(Dense(11, 11, NNlib.relu)), ResDense(Dense(11, 11, NNlib.relu)), ResDense(Dense(1, 1, tanh)))
# TODO fix this
julia> FluxExtensions.layerbuilder(4,11,1,3,"relu","linear","ResDense")
Chain(ResDense(Dense(11, 11, NNlib.relu)), ResDense(Dense(11, 11, NNlib.relu)), Dense(11, 1))
```
"""
layerbuilder(k::Vector{Int},l::Vector,f::Vector) = Flux.Chain(map(i -> i[1](i[3],i[4],i[2]),zip(l,f,k[1:end-1],k[2:end]))...)

layerbuilder(d::Int,k::Int,o::Int,n::Int, args...) =
	layerbuilder(vcat(d,fill(k,n-1)...,o), args...)

function layerbuilder(ks::Vector{Int},ftype::String,lastlayer::String = "",ltype::String = "Dense")
	ftype = (ftype == "linear") ? "identity" : ftype
	ls = Array{Any}(fill(eval(:($(Symbol(ltype)))),length(ks)-1))
	fs = Array{Any}(fill(eval(:($(Symbol(ftype)))),length(ks)-1))
	if !isempty(lastlayer)
		fs[end] = (lastlayer == "linear") ? identity : eval(:($(Symbol(lastlayer))))
		ls[end] = (lastlayer == "linear") ? Dense : ls[end]
	end
	layerbuilder(ks,ls,fs)
end


text:
type HiddenMarkovModel <: DynamicDiscreteModel
	#HMM specific field

	#DynamicDiscreteModel fields
	m::Array{Float64,4}			  	#the transition matrix given as m[x,y,x',y']
	mu::Array{Float64,2}  			#initial distribution (dx,dy)
	rho::Array{Float64,1}
	phi::Array{Float64,1}
	psi::Array{Float64,1}
end

function hmm(ab::Tuple,mu)
	a,b=ab
	dx,dy=size(mu)
	model=HiddenMarkovModel(Array(Float64,dx,dy,dx,dy),mu,Array(Float64,1),Array(Float64,dx),Array(Float64,dx))
	hmm2ddm!(model,a,b)
	model
end

#default initial distribution uniform
function hmm(ab::Tuple)
	dx,dy=size(ab[2])
	mu=fill(1/(dx*dy),dx,dy)
	hmm(ab,mu)
end

#ab=(a,b) parametrization
coef!(model::HiddenMarkovModel,ab::Tuple)=hmm2ddm!(model,ab[1],ab[2])

#canonical parametrization in case someone wants to apply mle() to it
coef!(model::HiddenMarkovModel,theta::Array{Float64,1})=coef!(model,theta2ab(theta,size(model.mu,1)))

theta2ab(theta::Array{Float64,1},dx::Int)=(z2q(theta[1:dx*(dx-1)]),z2q(theta[dx*(dx-1)+1:end],dx))

function dim(model::HiddenMarkovModel)
	dx,dy=size(model.mu)
	dx*(dx-1)+dx*(dy-1)
end

function em(model::HiddenMarkovModel,data)
	dx,dy=size(model.mu)
	w=zeros(dx,dy,dx,dy)
	estep(model,data,w)
	a=nsm(reshape(sum(w,(2,4)),dx,dx))
	b=nsm(reshape(sum(w,(1,2)),dx,dy))
	llk=loglikelihood(model,data,(a,b))
	println()
	println(" estimating hidden Markov model via Baum-Welch algorithm...")
	println("  log-likelihood: $(round(llk,4))")
	(a,b)
end

baumwelch=em


text:
@testset "KSP{$ST}" begin
  A = PETSc.Mat(ST, 3,3)
  A_julia = zeros(ST,3,3)
  for i=1:3
    for j=1:3
      pos = (i-1)*3 + j
      val = RC(Complex(pos, pos))
      A[i, j] = val
      A_julia[i,j] = val
     end
  end

  A[3,2] = RC(Complex(6,6))  # make A non singular
  A_julia[3,2] = RC(Complex(6,6))

  b = PETSc.Vec(ST, 3, vtype=PETSc.C.VECMPI)
  b_julia = zeros(ST, 3)
  b[3] = RC(Complex(1,1))
  b_julia[3] = RC(Complex(1,1))

  @testset "ksp GMRES solves" begin
      kspg = PETSc.KSP(A)
      x = kspg\b
      x_julia = A_julia\b_julia
      @test x ≈ x_julia

      x = A\b
      @test x ≈ x_julia

      x = KSPSolveTranspose(kspg, b)
      x_julia = A_julia.'\b_julia
      @test x ≈ x_julia

      x = KSPSolveTranspose(A, b)
      @test x ≈ x_julia

      pc   = PETSc.PC(ST,comm=comm(kspg),pc_type="jacobi")
      PETSc.chk(PETSc.C.PCSetOperators(pc.p,A.p,A.p))
      kspg = PETSc.KSP(pc)
      x = kspg\b
      x_julia = A_julia\b_julia

      @test x ≈ x_julia
      # fixme: test petscview without side effects, e.g.
      # via PetscViewerStringOpen
  end

  @testset "ksp BCGS solves" begin
      kspb = PETSc.KSP(A, ksp_type="bcgs")
      x = kspb\b
      x_julia = A_julia\b_julia
      @test x ≈ x_julia
  end
end


text:
using Test
using Openresty
using HTTP

function createconfig(workdir::String, sudo::Bool=false, log_to_files::Bool=true, daemon::Bool=false)
    me = ENV["USER"]
    user = sudo ? "user $me $me;" : ""
    error_log = string("error_log ", log_to_files ? "$workdir/logs/error.log" : "/dev/stderr", " debug;")
    access_log = string("access_log ", log_to_files ? "$workdir/logs/access.log" : "/dev/stdout", ";")
    daemon_mode = daemon ? "on" : "off";

    cfg = """$user
daemon $daemon_mode;
worker_processes  1;
$error_log
events {
    worker_connections  1024;
}
http {
    $access_log
    $error_log
    lua_package_path 'OPENRESTY_LUA_PACKAGE_PATH';
    lua_package_cpath 'OPENRESTY_LUA_PACKAGE_CPATH';
    include       mime.types;
    server {
        listen 8080;
        server_name localhost;
        root $workdir/html;
        index index.html index.htm index.nginx-debian.html;
    }
}
"""
    cfgfile = joinpath(workdir, "test.conf")

    open(cfgfile, "w") do f
        println(f, cfg)
    end

    cfgfile
end

function test_nginx_config()
    body = ""

    @info("waiting for nginx to come up")
    # try for 10 secs
    while isempty(body)
        sleep(2)
        try
            resp = HTTP.get("http://127.0.0.1:8080/")
            body = String(resp.body)
        catch ex
            @info("nginx not ready yet...")
        end
    end

    isempty(body) && error("nginx did not come up")

    @info("testing nginx response")
    @test occursin("Welcome to OpenResty", body)
    nothing
end

function count_requests(logcontents::String)
    nrequests = 0
    for line in readlines(IOBuffer(logcontents))
        if match(r"\"GET\s\/\sHTTP\/1\.1\" 200", line) !== nothing
            nrequests += 1
        end
    end
    nrequests
end

function test(; sudo::Bool=false, log_to_files::Bool=true, daemon::Bool=false)
    @info("starting tests", sudo, log_to_files, daemon)
    workdir = mktempdir()
    mkpath(workdir)
    cfgfile = createconfig(workdir, sudo, log_to_files, daemon)
    nginx = OpenrestyCtx(workdir; sudo=sudo)

    accesslog = nothing
    errorlog = nothing
    accesslog2 = nothing
    errorlog2 = nothing
    if !log_to_files
        accesslog = errorlog = PipeBuffer()
    end

    @info("setting up Openresty", workdir)

    # incorrect lua path should throw error
    @test nothing === setup(nginx, cfgfile; lua_package_path="~/lua/?.lua", lua_package_cpath="~/lua/?.so")
    @test isfile(Openresty.conffile(nginx))
    confstr = read(Openresty.conffile(nginx), String)
    @test occursin("~/lua/?.lua", confstr)
    @test occursin(Openresty.luapath, confstr)
    @test occursin("$(Openresty.luapath);~/lua/?.lua;;", confstr)
    @test occursin("~/lua/?.so", confstr)
    @test occursin(Openresty.luacpath, confstr)
    @test occursin("$(Openresty.luacpath);~/lua/?.so;;", confstr)

    @info("starting Openresty")
    start(nginx; accesslog=accesslog, errorlog=errorlog)
    sleep(2)
    @test isfile(Openresty.pidfile(nginx))
    @test isrunning(nginx)

    test_nginx_config()

    @info("restarting Openresty")
    if !log_to_files
        # we need to create new IOs because they may be closed when the process restarts
        accesslog2 = errorlog2 = PipeBuffer()
        flush(accesslog)
    end
    restart(nginx; delay_seconds=2, accesslog=accesslog2, errorlog=errorlog2)

    sleep(2)
    @test isfile(Openresty.pidfile(nginx))
    @test isrunning(nginx)
    @test nothing === reopen(nginx)
    @test nothing === reload(nginx)

    test_nginx_config()

    @info("stopping Openresty")
    stop(nginx)
    sleep(2)
    @test !isfile(Openresty.pidfile(nginx))
    @test !isrunning(nginx)
    if log_to_files
        @test isfile(Openresty.accesslogfile(nginx))
        @test isfile(Openresty.errorlogfile(nginx))
    else
        flush(accesslog2)
        logbytes = readavailable(accesslog)
        logbytes2 = readavailable(accesslog2)
        @test !isempty(logbytes)
        @test !isempty(logbytes2)
        strbytes = String(logbytes) * String(logbytes2)
        @test findfirst("start worker processes", strbytes) !== nothing
        @test findfirst("HTTP", strbytes) !== nothing

        # test that when daemon is off, we are able to capture access logs even after a reload
        # but in daemon mode we are not able to (because of https://trac.nginx.org/nginx/ticket/1715#ticket)
        if daemon
            @test count_requests(strbytes) == 1
        else
            @test count_requests(strbytes) == 2
        end
    end

    @test_throws Exception setup(nginx, cfgfile)
    @test nothing === setup(nginx, cfgfile; force=true)
    @test nothing === setup(nginx, cfgfile; force=true, lua_package_path=["~/lua/?.lua", "/a/different/path"], lua_package_cpath=["~/lua/?.so", "/a/different/cpath"])

    @info("cleaning up")
    rm(workdir; recursive=true, force=true)
    @info("done")

    nothing
end

function sudo_available()
    cmd = `sudo -n ls`
    try
        run(pipeline(`sudo -n ls`, stdout=devnull, stderr=devnull))
        return true
    catch
        return false
    end
end

for log_to_files in (true,false)
    for daemon in (true, false)
        test(; sudo=false, log_to_files=log_to_files, daemon=daemon)

        if sudo_available()
            test(; sudo=true, log_to_files=log_to_files)
        else
            @info("passwordless sudo not available, skipping sudo tests")
        end
    end
end


text:
function SimpleSyntheticHypergraph(n,m,K,pvals,rmin,rmax,cluster_sizes=ones(K),r_sizes=ones(rmax-rmin),cluster_prefs=ones(K))
    """
    Generate a synthetic hypergraph, store in incidence matrix
    Output several other things about hypergraph like edge sizes and degrees
    and edge list.

    This can be viewed as a simplified approximation to the HSBM, that is
    easier and faster to generate in practice.

    Procedure:
    * Set a fixed number of hyperedges m and clusters K (you can randomize this outside this function)
    * Set a distribution determining the relative proportion of hyperedges of each size from r_min to r_max
    * For each hyperedge, with probability p, assign it to one community at random,
        selecting a r nodes uniformly at random for a hyperedge of size r
    * Otherwise, with probability 1-p, select r nodes uniformly at random from across the entire hypergraph
    """

    # Normalize to get proportions/probabilities
    r_sizes = Weights(r_sizes/sum(r_sizes))
    cluster_sizes = Weights(cluster_sizes/sum(cluster_sizes))         # proportion of nodes in each cluster
    cluster_prefs = Weights(cluster_prefs/sum(cluster_prefs))        # relative proportion of interior hyperedge in each cluster

    # Generate the ground truth clustering of nodes
    # ground_truth = GenerateTruth(n,cluster_sizes)
    Clusters = Vector{Vector{Int64}}()
    for j = 1:K
        push!(Clusters,Vector{Int64}())
        @assert(n*cluster_sizes[j] > 2*rmax)
    end

    ground_truth = zeros(Int64,n)
    for i = 1:n
        c = sample(1:K,cluster_sizes)
        ground_truth[i] = c
        push!(Clusters[c],i)
    end

    U = Vector{Int64}()
    E = Vector{Int64}()
    EdgeList = Vector{Vector{Int64}}()
    E_lengths = zeros(Int64,m)
    for enum = 1:m
        e = SampleEdge(pvals,n,rmin,rmax,cluster_prefs,r_sizes,Clusters)
        push!(EdgeList,e)
        E_lengths[enum] = length(e)
        for node in e
            push!(U,node)
            push!(E,enum)
        end
    end
    He2n = SparseArrays.sparse(E,U,ones(length(U)),m,n)
    deg = vec(sum(He2n,dims = 1))
    return He2n, EdgeList, E_lengths, deg, ground_truth

end

function SampleEdge(pvals,n,rmin,rmax,cluster_prefs,r_sizes,Clusters)
    replace = false
    K = length(Clusters)
    r = sample(rmin:rmax,r_sizes)       # select a hyperedge size
    if rand(1)[1] < pvals[r-rmin+1]

        c = sample(1:K,cluster_prefs)   # select a random cluster to put the edge in
        return sort(sample(Clusters[c],r,replace = replace))
    else
        return sort(sample(1:n,r,replace = replace))
    end
end

max_stars_repo_name: nveldt/TestPkg

max_issues_repo_name: nveldt/TestPkg

max_forks_repo_name: nveldt/TestPkg


text:
module TestMap
    using Base.Test
    using NullableArrays

    # create m random arrays each with N dimensions and length dims[i] along
    # dimension i for i=1:N
    N = rand(2:5)
    dims = Int[ rand(3:8) for i in 1:N ]
    m = rand(4:6)
    n = rand(1:3)
    As = Array[ [rand(dims...) for i in 1:n] ; [rand(Int, dims...) for i in n+1:m] ]

    Ms = [ rand(Bool, dims...) for i in 1:m ]
    Xs = Array{NullableArray, 1}()
    for i in 1:m
        push!(Xs, NullableArray(As[i]))
    end
    Ys = Array{NullableArray, 1}()
    for i in 1:m
        push!(Ys, NullableArray(As[i], Ms[i]))
    end

    C = Array{Float64}(dims...)
    Z = NullableArray(Float64, dims...)

    R = map(|, Ms...)

    f(x...) = sum(x)
    dests = (C, Z)

    # 1 arg
    for (args, masks) in (
        ((As, Xs), fill((), m)), ((As, Ys), [ (Ms[i],) for i in 1:m ])
    )
        for i in 1:m
            # map!
            map!(f, args[1][i], args[1][i]) # map!(f, As[i], As[i])
            map!(f, args[2][i], args[2][i]) # map!(f, Xs[i], Xs[i])
            @test isequal(args[2][i], NullableArray(args[1][i], masks[i]...))
            # map
            A = map(f, args[1][i])
            X = map(f, args[2][i])
            @test isequal(X, NullableArray(A, masks[i]...))
        end
    end
    # 2 arg
    i, j = rand(1:m), rand(1:m)
    S = map(|, Ms[i], Ms[j])
    for (args, mask) in (
        ((As, Xs), ()), ((As, Ys), (S,))
    )
        # map!
        map!(f, dests[1], args[1][i], args[1][j])
        map!(f, dests[2], args[2][i], args[2][j])
        @test isequal(dests[2], NullableArray(dests[1], mask...))
        # map
        map(f, args[1][i], args[1][j])
        map(f, args[2][i], args[2][j])
        @test isequal(dests[2], NullableArray(dests[1], mask...))
    end
    # n arg
    for (args, mask) in (
        ((As, Xs), ()), ((As, Ys), (R,))
    )
        # map!
        map!(f, dests[1], args[1]...)
        map!(f, dests[2], args[2]...)
        @test isequal(dests[2], NullableArray(dests[1], mask...))
        # map
        map(f, args[1]...)
        map(f, args[2]...)
        @test isequal(dests[2], NullableArray(dests[1], mask...))
    end

    # test map over empty NullableArrays
    X = NullableArray(Int[])
    h1(x) = 5.0*x
    h2(x) = x
    h2(x...) = prod(x)

    Z1 = map(h1, X)
    Z2 = map(h2, X)
    @test isempty(Z1)
    @test isa(Z1, NullableArray{Float64})
    @test isempty(Z2)
    @test isa(Z2, NullableArray{Int})

    # if a function has no method for inner eltype of empty NullableArray,
    # result should be empty NullableArray{Any}() for consistency with generic map()
    h3(x::Float64...) = prod(x)
    Z3 = map(h3, X)
    @test isempty(Z3)
    @test isa(Z3, NullableArray{Any})
    Z3 = map(h3, X, X)
    @test isempty(Z3)
    @test isa(Z3, NullableArray{Any})
    Z3 = map(h3, X, X, X)
    @test isempty(Z3)
    @test isa(Z3, NullableArray{Any})

    # test map over all null NullableArray
    n = rand(10:100)
    Ys = [ NullableArray(rand(Int, n), fill(true, n)) for i in 1:rand(3:5) ]

    Z2 = map(h2, Ys[1])
    @test isequal(Z2, NullableArray(Int, n))
    @test isa(Z2, NullableArray{Int})
    Z2 = map(h2, Ys[1], Ys[2])
    @test isequal(Z2, NullableArray(Int, n))
    @test isa(Z2, NullableArray{Int})
    Z2 = map(h2, Ys...)
    @test isequal(Z2, NullableArray(Int, n))
    if VERSION >= v"0.6.0-dev.2544" # Commit known to work
        @test isa(Z2, NullableArray{Int})
    end
end # module TestMap


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

"""
    struct RenormalizedDensity <: AbstractDensity

Constructors:

* ```$(FUNCTIONNAME)(density::AbstractDensity, logrenormf::Real)```

A renormalized density derived from `density`, with

```julia
logdensityof(RenormalizedDensity(density, logrenormf))(v) ==
    logdensityof(density)(v) + logrenormf
```
"""
struct RenormalizedDensity{D<:AbstractDensity,T<:Real} <: AbstractDensity
    density::D
    logrenormf::T
end

Base.parent(density::RenormalizedDensity) = density.density

Base.:(==)(a::RenormalizedDensity, b::RenormalizedDensity) = a.density == b.density && a.logrenormf == b.logrenormf

var_bounds(density::RenormalizedDensity) = var_bounds(parent(density))

ValueShapes.varshape(density::RenormalizedDensity) = varshape(parent(density))

ValueShapes.unshaped(density::RenormalizedDensity) = RenormalizedDensity(unshaped(density.density), density.logrenormf)

(shape::AbstractValueShape)(density::RenormalizedDensity) = RenormalizedDensity(shape(density.density), density.logrenormf)

function Base.show(io::IO, d::RenormalizedDensity)
    print(io, Base.typename(typeof(d)).name, "(")
    show(io, d.density)
    print(io, ", ")
    show(io, d.logrenormf)
    print(io, ")")
end

function DensityInterface.logdensityof(density::RenormalizedDensity, v::Any)
    parent_logd = logdensityof(parent(density),v)
    R = float(typeof(parent_logd))
    convert(R, parent_logd + density.logrenormf)
end

function checked_logdensityof(density::RenormalizedDensity, v::Any)
    parent_logd = checked_logdensityof(parent(density),v)
    R = float(typeof(parent_logd))
    convert(R, parent_logd + density.logrenormf)
end

Distributions.sampler(density::RenormalizedDensity) = Distributions.sampler(parent(density))
bat_sampler(density::RenormalizedDensity) = bat_sampler(parent(density))

Statistics.cov(density::RenormalizedDensity) = cov(parent(density))

"""
    BAT.renormalize_density(density::AbstractDensity, logrenormf::Real)::AbstractDensity

*Experimental feature, not part of stable public API.*

Renormalies `density` with the logarithmic renormalization factor, so that

```julia
logdensityof(renormalize_density(density, logrenormf))(v) ==
    logdensityof(density)(v) + logrenormf
```
"""
function renormalize_density end
export renormalize_density

renormalize_density(density::Any, logrenormf::Real) = renormalize_density(convert(AbstractDensity, density), logrenormf)

renormalize_density(density::AbstractDensity, logrenormf::Real) = RenormalizedDensity(density, logrenormf)

renormalize_density(density::RenormalizedDensity, logrenormf::Real) = renormalize_density(parent(density), density.logrenormf + logrenormf)


text:
using Genie.Router
import GraphqlController

GraphqlController.register()

route("/") do
  serve_static_file("welcome.html")
end


text:
module TestAnalysis
using CustomTest
using Base.Test
using ControlSystems

@test_approx_eq stabregionPID("exp(-sqrt(s))", logspace(-3,3), kd=1)[2][1] -1.022356911142034

P = tf(1,[1,1])^4
gangoffourplot(P,tf(1))
ωp = 0.8
kp,ki,C = loopshapingPI(P,ωp,phasemargin=60, doplot=true)
@test_approx_eq kp 0.82274734724854
@test_approx_eq ki -0.45472580052281375


text:
module Williamsom2ThetaMethodFullNewtonTestsMPI

using  PartitionedArrays
using  Test
using  FillArrays
using  Gridap
using  GridapDistributed
using  GridapGeosciences
using  GridapPETSc
using  SparseArrays
using  MPI
using  SparseMatricesCSR

const PArrays=PartitionedArrays

include("../sequential/Williamson2InitialConditions.jl")

gamg = """
       -ksp_type gmres -ksp_rtol 1.0e-06 -ksp_atol 0.0 -ksp_gmres_restart 100
       -ksp_monitor -pc_type gamg -pc_gamg_type agg
       -mg_levels_esteig_ksp_type gmres -mg_coarse_sub_pc_type lu
       -mg_coarse_sub_pc_factor_mat_ordering_type nd -pc_gamg_process_eq_limit 50
       -pc_gamg_square_graph 9 pc_gamg_agg_nsmooths 1
       """

mumps = """
        -ksp_type preonly -ksp_error_if_not_converged true
        -pc_type lu -pc_factor_mat_solver_type mumps
        -mat_mumps_icntl_1 4
        -mat_mumps_icntl_14 100
        -mat_mumps_icntl_28 1
        -mat_mumps_icntl_29 2
        -mat_mumps_cntl_3 1.0e-6
        """

options = """
          -snes_type newtonls
          -snes_linesearch_type basic
          -snes_linesearch_damping 1.0
          -snes_rtol 1.0e-8
          -snes_atol 0.0
          -snes_monitor
          -snes_converged_reason
          -mm_ksp_type cg
          -mm_ksp_monitor
          -mm_ksp_rtol 1.0e-14
          -mm_pc_type jacobi
          """
function mysnessetup(snes)
  ksp      = Ref{GridapPETSc.PETSC.KSP}()
  pc       = Ref{GridapPETSc.PETSC.PC}()
  mumpsmat = Ref{GridapPETSc.PETSC.Mat}()
  @check_error_code GridapPETSc.PETSC.SNESSetFromOptions(snes[])
  @check_error_code GridapPETSc.PETSC.SNESGetKSP(snes[],ksp)
  @check_error_code GridapPETSc.PETSC.KSPView(ksp[],C_NULL)
  @check_error_code GridapPETSc.PETSC.KSPSetType(ksp[],GridapPETSc.PETSC.KSPPREONLY)
  @check_error_code GridapPETSc.PETSC.KSPGetPC(ksp[],pc)
  @check_error_code GridapPETSc.PETSC.PCSetType(pc[],GridapPETSc.PETSC.PCLU)
  @check_error_code GridapPETSc.PETSC.PCFactorSetMatSolverType(pc[],GridapPETSc.PETSC.MATSOLVERMUMPS)
  @check_error_code GridapPETSc.PETSC.PCFactorSetUpMatSolverType(pc[])
  @check_error_code GridapPETSc.PETSC.PCFactorGetMatrix(pc[],mumpsmat)
  @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[],  4, 1)
  # percentage increase in the estimated working space
  @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[],  14, 1000)
  @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 28, 2)
  @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 29, 2)
  @check_error_code GridapPETSc.PETSC.MatMumpsSetCntl(mumpsmat[], 3, 1.0e-6)
end

function set_ksp_mm(ksp)
  @check_error_code GridapPETSc.PETSC.KSPSetOptionsPrefix(ksp[],"mm_")
  @check_error_code GridapPETSc.PETSC.KSPSetFromOptions(ksp[])
  #@check_error_code GridapPETSc.PETSC.KSPView(ksp[],C_NULL)
end

function main(parts)
  GridapPETSc.with(args=split(options)) do
      if (PArrays.get_part_id(parts)==1)
        println(options)
      end

      # Solves the steady state Williamson2 test case for the shallow water equations on a sphere
      # of physical radius 6371220m. Involves a modified coriolis term that exactly balances
      # the potential gradient term to achieve a steady state
      # reference:
      # D. L. Williamson, J. B. Drake, J. J.HackRüdiger Jakob, P. N.Swarztrauber, (1992)
      # J Comp. Phys. 102 211-224

      l2_err_u = [0.011370921987771046 , 0.002991229234176333 ]
      l2_err_h = [0.005606685579166809, 0.001458107077200681 ]
      for matrix_type in [SparseMatrixCSR{0,GridapPETSc.PetscScalar,GridapPETSc.PetscInt},
                          SparseMatrixCSC{Float64,Int}]
          order=1
          degree=4
          θ=0.5
          for i in 1:1
            n      = i+1
            nstep  = 5*2^n
            Uc     = sqrt(g*H₀)
            dx     = 2.0*π*rₑ/(4*n)
            dt     = 0.25*dx/Uc
            println("timestep: ", dt)   # gravity wave time step
            T      = dt*nstep
            τ      = dt/2
            model  = CubedSphereDiscreteModel(parts, n; radius=rₑ)
            nls    = PETScNonlinearSolver(mysnessetup)
            mmls   = PETScLinearSolver(set_ksp_mm)

            function ts()
              shallow_water_theta_method_full_newton_time_stepper(
                nls, model, order, degree,
                h₀, u₀, f₀, topography,
                g, θ, T, nstep, τ;
                mass_matrix_solver=mmls,
                matrix_type=matrix_type,
                am_i_root=PArrays.get_part_id(parts)==1,
                write_solution=false,
                write_solution_freq=5,
                write_diagnostics=true,
                write_diagnostics_freq=1,
                dump_diagnostics_on_screen=true)
            end

            if (PArrays.get_part_id(parts)==1)
              @time hf, uf, _ = ts()
            else
              hf, uf, _ = ts()
            end

            Ω     = Triangulation(model)
            dΩ    = Measure(Ω, degree)
            hc    = CellField(h₀, Ω)
            e     = h₀-hf
            err_h = sqrt(sum(∫(e⋅e)*dΩ))/sqrt(sum(∫(hc⋅hc)*dΩ))
            uc    = CellField(u₀, Ω)
            e     = u₀-uf
            err_u = sqrt(sum(∫(e⋅e)*dΩ))/sqrt(sum(∫(uc⋅uc)*dΩ))
            if PArrays.get_part_id(parts)==1
              println("n=", n, ",\terr_u: ", err_u, ",\terr_h: ", err_h)
            end
            #@test abs(err_u - l2_err_u[i]) < 10.0^-12
            #@test abs(err_h - l2_err_h[i]) < 10.0^-12
          end
      end
  end
end
MPI.Init()
prun(main,mpi,MPI.Comm_size(MPI.COMM_WORLD))

end # module

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
"""
    orthogonalitydefect(B)

Find the orthogonality defect of the matrix B defined, for example,
on page 2 of Bennet

[Bennet](http://users.eecs.northwestern.edu/~hbennett/publications/bases.pdf)

# Examples
```jldoctest
julia> H= [1 2; 3 4];B,T=lll(H);

julia> [orthogonalitydefect(H) orthogonalitydefect(B)]
1×2 Array{Float64,2}:
 7.07107  1.0

```
"""
function orthogonalitydefect(B)
    N = size(B,2)
    prod = one(B[1,1])
    for n=1:N
        prod *= sqrt(B[:,n]'*B[:,n])
    end
    return prod/abs(det(B))
end

"""
    hermitefactor(B)

Find the Hermite factor of matrix B

# Examples
```jldoctest
julia> H= [1 2; 3 4];hermitefactor(H)
1.5811388300841898

```
"""
hermitefactor(B) = sqrt(B[:,1]'*B[:,1])/abs(det(B))

"""
    seysencond(B)

Seysen condition number as on, for example, page 3 of Bennet

[Bennet](http://users.eecs.northwestern.edu/~hbennett/publications/bases.pdf)

# Examples
```jldoctest
julia> H= [1 2; 3 4];seysencond(H)
2.8284271247461903

```
"""
function seysencond(B)
    N = size(B,2)
    Bi = inv(B)'
    quotients = zeros(typeof(Bi[1,1]),N)
    for n=1:N
        quotients[n] = sqrt(B[:,n]'*B[:,n])/sqrt(Bi[:,n]'*Bi[:,n])
    end
    return maximum(quotients)
end

"""
    islllreduced(B)

Determine if the matrix B is LLL reduced or not. See p 56 of Bremner for a
definition.

M. R. Bremner, "Lattice Basis Reduction: An Introduction to the LLL
 Algorithm and Its Applications" CRC Press, 2012.

# Examples
```jldoctest
julia> H= [1 2; 3 4];islllreduced(H)
false

julia> B,T=lll(H);islllreduced(B)
true

```
"""
function islllreduced(B)
# follows p 56 of Bremner
    Td = eltype(B)
    Xs,M = gso(B)
    isSizeReduced = all((M-I)[:]*2 .< one(Td))
    n = size(B,2)
    xiSq = zeros(float(Td),n)
    for i=1:n
        xiSq[i] = Xs[:,i]'*Xs[:,i]
    end
    deltas = zeros(float(Td),n-1)
    for i=2:n
        deltas[i-1] = xiSq[i]/xiSq[i-1]+M[i-1,i]^2
    end
    if isSizeReduced && minimum(deltas)>1/4
        return true
    else
        return false
    end
end

"""
    issizereduced(B)

Determine if the matrix B is size reduced or not.

# Examples
```jldoctest
julia> H= [1 2; 3 4];issizereduced(H)
false

julia> B,T = lll(H);issizereduced(B)
true

```
"""
function issizereduced(B)
    # to-do: replace gso with QR as on p 148 of LLL survey book
    Xs,M = gso(B)
    return all((M-I)[:]*2 .< one(eltype(B)))
end

"""
    Q,R = gso(H)

Find Gram-Schmidt orthogonal basis; H=Q*R
This is not quite a QR decomposition.

# Examples
```jldoctest
julia> H = [1 2; 3 4]; Q,R = gso(H)

```
"""
function gso(H::Matrix{Td}) where {Td}
# This is classic GS; see comparison between classic and modified in
# http://www.cis.upenn.edu/~cis610/Gram-Schmidt-Bjorck.pdf

# See also Fig 3.1 of "Lattice Basis Reduction: An Introduction to the LLL
# Algorithm and Its Applications" by Murray R. Bremner, CRC Press, 2012.

    Q = float.(H)            # Xs
    n = size(H,2)
    R = Matrix{float(Td)}(I,n,n) # M or μ

    for i=1:n
        for j=1:i-1
            R[j,i] = H[:,i]'*Q[:,j]/(Q[:,j]'*Q[:,j])
            Q[:,i]-= R[j,i]*Q[:,j]
        end
    end
    return Q,R
end

"""
    volume(B)

Volume of fundamental parallelepiped of a lattice with basis B.

# Examples
```jldoctest
julia> B = [1 2; 3 4]; volume(B)
1.9999999999999964

```
"""
volume(B) = sqrt(det(B'*B))


text:
using SDEs, Base.Test
include("setup.jl")
include("test_utils.jl")

Ncheck = 1000
tolerance = 2/sqrt(Ncheck) #CLT for comparing processes with different random number generators
id_tolerance = 1e-7 #For comparing things with same random number generators
θ=0.5
η=0.5

println(prob, ", dt:", t[2]-t[1], ", θ:",θ, ", η:",η)

EMPdev = check(prob, x0, t, SDEs.EMproto(), 1:Ncheck)
println("proto EM  :", EMPdev)

EMdev = check(prob, x0, t, EM(), 1:Ncheck)
println("EM        :", EMdev)

DeqEMdev = check(prob, x0, t, DEQEM(), 1:Ncheck)
println("Diffeq    :", DeqEMdev)

PCEPdev = check(prob, x0, t, SDEs.PCEproto(), 1:Ncheck; θ=θ, η=η)
println("proto PCE :", PCEPdev)

PCEdev = check(prob, x0, t, PCE(), 1:Ncheck; θ=θ, η=η)
println("PCE       :", PCEdev)

@test EMdev ≈ DeqEMdev atol=EMdev*tolerance
@test EMdev ≈ EMPdev atol=EMdev*id_tolerance
@test PCEPdev ≈ PCEdev atol=PCEdev*id_tolerance
@test PCEdev < EMdev

max_stars_repo_name: onoderat/SDEs

max_issues_repo_name: onoderat/SDEs

max_forks_repo_name: onoderat/SDEs


text:
"""
Analytic axisymmetric medium-size tokamak equilibrium in cartesian coordinates.
"""
module TokamakMediumCylindrical

    import ElectromagneticFields.AxisymmetricTokamakCylindrical

    export initial_conditions_barely_passing, initial_conditions_barely_trapped,
           initial_conditions_deeply_passing, initial_conditions_deeply_trapped

    export hamiltonian, toroidal_momentum

    AxisymmetricTokamakCylindrical.@code(2., 5., 2.) # inject magnetic field code

    initial_conditions_barely_passing() = ([2.5, 0., 0., 3.425E-1], (μ = 1E-2,)) # Δt=2.5, nt=50
    initial_conditions_barely_trapped() = ([2.5, 0., 0., 3.375E-1], (μ = 1E-2,)) # Δt=3.0, nt=100
    initial_conditions_deeply_passing() = ([2.5, 0., 0., 5E-1], (μ = 1E-2,))     # Δt=2.5, nt=25
    initial_conditions_deeply_trapped() = ([2.5, 0., 0., 1E-1], (μ = 1E-2,))     # Δt=5.0, nt=50

    μ_loop() = 1E-3
    μ_surface() = 1E-3

    function f_loop(s)
        R0 = 1.75
        Z0 = 0.0
        φ0 = 0.0
        u0 = 0.5
        rx = 0.1
        ry = 0.1

        Rs = R0 + rx*cos(2π*s)
        Zs = Z0 + ry*sin(2π*s)

        qs = [Rs, Zs, φ0, u0]

        return qs
    end

    function f_surface(s,t)
        R0 = 1.75
        Z0 = 0.0
        φ0 = 0.0
        u0 = 0.5
        r0 = 0.1

        Rt = R0 + r0*(s-0.5)
        Zt = Z0 + r0*(t-0.5)

        qt = [Rt, Zt, φ0, u0]

        return qt
    end

    include("guiding_center_4d_common.jl")
    include("guiding_center_4d_equations.jl")
    include("guiding_center_4d_loop.jl")
    include("guiding_center_4d_surface.jl")

    function toroidal_momentum(t,q)
        ϑ₃(t,q)
    end

    include("guiding_center_4d_diagnostics.jl")

end