text:
#=========================================================================================
    datastreams.jl

    This is the datastreams interface for Estuaries.  Basically just re-implements
    the interface for its source.
=========================================================================================#

#=========================================================================================
    <source interface>
=========================================================================================#
Data.schema(E::Estuary) = E.schema
schema(E::Estuary) = Data.schema(E)
schema!(E::Estuary) = (E.schema = Data.schema(E.src))
export schema, schema!

Data.isdone(E::Estuary, row, col) = Data.isdone(E.src, row, col)
Data.reference(E::Estuary) = Data.reference(E.src)

# size is implemented in the Base interface section in gradine.jl

Data.schema(E::Estuary, ::Type{Data.Field}) = Data.schema(E.src, Data.Field)

Data.streamtype(E::Estuary, ::Type{Data.Field}) = Data.streamtype(E.src, Data.Field)

# this for both regular and nullable types
function Data.streamfrom{T}(E::Estuary, ::Type{Data.Field}, ::Type{T}, row, col)
    Data.streamfrom(E.src, Data.Field, T, row, col)
end

Data.schema(E::Estuary, ::Type{Data.Column}) = Data.schema(E.src, Data.Column)

Data.streamtype(E::Estuary, ::Type{Data.Column}) = Data.streamtype(E.src, Data.Column)

# this for both regular and nullable types
function Data.streamfrom(E::Estuary, ::Type{Data.Column}, ::Type{T}, row, col) where T
    Data.streamfrom(E.src, Data.Column, T, row, col)
end
# TODO should fix other methods to pass 1 as argument
function Data.streamfrom(E::Estuary, ::Type{Data.Column}, ::Type{T}, col) where T
    Data.streamfrom(E, Data.Column, T, 1, col)
end
#=========================================================================================
    </source interface>
=========================================================================================#

#=========================================================================================
    <sink interface>
=========================================================================================#
Data.streamtypes(E::Estuary) = Data.streamptypes(sink)

# does this need special constructors???

function Data.streamto!(E::Estuary, ::Type{Data.Field}, val, row, col)
    Data.streamto!(E.sink, Data.Field, val, row, col)
end
function Data.streamto!(E::Estuary, ::Type{Data.Field}, val, row, col, schema::Data.Schema)
    Data.streamto!(E.sink, Data.Field, val, row, col, schema)
end

function Data.streamto!(E::Estuary, ::Type{Data.Column}, v, row, col)
    Data.streamto!(E.sink, Data.Column, v, row, col)
end
function Data.streamto!(E::Estuary, ::Type{Data.Column}, v, row, col, schema::Data.Schema)
    Data.streamto!(E.sink, Data.Column, v, row, col, schema)
end

Data.cleanup!(E::Estuary) = Data.cleanup!(E.sink)
Data.close!(E::Estuary) = Data.close!(E.sink)
#=========================================================================================
    </sink interface>
=========================================================================================#

#=========================================================================================
    <schema extensions>  consider adding to DataStreams.jl

    Note that usually objects with the DataStreams interface do not support operations
    that change the schema. That was probably the original intention.
=========================================================================================#
function altercolumn!{T}(schema::Data.Schema, idx::Integer, name::String, ::Type{T})
    schema.header[idx] = name
    schema.types[idx] = T
    schema
end

function insertcolumn!{T}(schema::Data.Schema, idx::Integer, name::String, ::Type{T})
    insert!(schema.header, idx, name)
    insert!(schema.types, idx, T)
    schema.cols += 1
    schema.index = Dict(n=>i for (i, n) ∈ enumerate(schema.header))
    schema
end

function appendcolumn!{T}(schema::Data.Schema, name::String, ::Type{T})
    push!(schema.header, name)
    push!(schema.types, T)
    schema.cols += 1
    schema.index[name] = schema.cols
    schema
end

function deletecolumn!(schema::Data.Schema, idx::Integer)
    name = splice!(schema.header, idx)
    deleteat!(schema.types, idx)
    schema.cols -= 1
    delete!(schema.index, name)
    schema
end

function deletecolumn!(schema::Data.Schema, name::String)
    idx = find(s -> s == name, schema.header)
    deleteat!(schema.header, idx)
    deleteat!(schema.types, idx)
    schema.cols -= 1
    delete!(schema.index, name)
    schema
end

# for now we leave on bounds checking
function Base.getindex(sch::Data.Schema, idx::Integer)
    if !(1 ≤ idx ≤ size(sch,2))
        throw(ArgumentError("Column index out of bounds."))
    end
    idx
end

Base.getindex(sch::Data.Schema, idx::Symbol) = sch[string(idx)]
Base.getindex(sch::Data.Schema, idx::AbstractVector{String}) = [sch[i] for i ∈ idx]
Base.getindex(sch::Data.Schema, idx::AbstractVector{<:Integer}) = [sch[i] for i ∈ idx]
Base.getindex(sch::Data.Schema, idx::AbstractVector{Symbol}) = [sch[string(i)] for i ∈ idx]

Base.haskey(sch::Data.Schema, idx::Integer) = (1 ≤ idx ≤ size(sch,2))
Base.haskey(sch::Data.Schema, idx::String) = haskey(sch.index, idx)
Base.haskey(sch::Data.Schema, idx::Symbol) = haskey(sch.index, string(idx))
#=========================================================================================
    </schema extensions>
=========================================================================================#


text:
using DSGE, OrderedCollections, CSV, DataFrames
using ClusterManagers, HDF5
# import DSGE: usual_settings!, usual_forecast!
using ModelConstructors, Dates, JLD, JLD2, SMC, StateSpaceRoutines, Nullables

include("../../../src/estimate/param_regimes.jl")
include("util_brookings.jl")

switch_all = false
switch_policy = false
switch_kappa = true
run_stored = true # Only set to true when on the Brookings_PC branch of DSGE.
## When using regime-switching parameters, set to false.

# Kalman Filter when run_stored is true
function kalman_stored(m,df,system)

data = df_to_matrix(m, df; cond_type = :none, in_sample = true)
start_date = max(date_presample_start(m), df[1, :date])
s_0 = Vector{Float64}(undef, 0)
P_0 = Matrix{Float64}(undef, 0, 0)
include_presample = true
outputs = [:loglh, :pred, :filt]
tol = 0.0

    T = size(data, 2)

    T = size(data, 2)
    if !isempty(data)
        # Calculate the number of periods since start date for each regime
        if haskey(m.settings, :n_hist_regimes)
            n_regime_periods = Vector{Int}(undef, get_setting(m, :n_hist_regimes) + get_setting(m, :n_cond_regimes))
            for k in 1:(get_setting(m, :n_hist_regimes)  + get_setting(m, :n_cond_regimes))
                n_regime_periods[k] = subtract_quarters(get_setting(m, :regime_dates)[k], start_date)
            end
        else
            n_regime_periods = Vector{Int}(undef, length(get_setting(m, :regime_dates)))
            for (k, v) in get_setting(m, :regime_dates)
                n_regime_periods[k] = subtract_quarters(v, start_date)
            end
        end

        if start_date < date_presample_start(m)
            error("Start date $start_date must be >= date_presample_start(m)")
        elseif 0 < subtract_quarters(date_zlb_start(m), start_date) < T
            regime_inds = Vector{UnitRange{Int64}}(undef, length(n_regime_periods) + 1)

            n_nozlb_periods = subtract_quarters(date_zlb_start(m), start_date) # number of periods since start date for start of ZLB

            # Get index of next regime after ZLB starts.
            # Note that it cannot be 1 b/c the first regime starts at the start date
            i_splice_zlb = findfirst(n_nozlb_periods .< n_regime_periods)

            if isnothing(i_splice_zlb)
                check_zlb_coincide, i_zlb_start = if n_nozlb_periods == n_regime_periods[end]
                    n_regime_periods[end], length(n_regime_periods)
                else
                    -1, 0
                end
            else
                i_zlb_start = i_splice_zlb - 1
                check_zlb_coincide = n_regime_periods[i_zlb_start]
            end
            if check_zlb_coincide == n_nozlb_periods
                # In this case, the post-ZLB coincides with a specific regime,
                # so we return the regime indices w/out splicing the ZLB in.
                return regime_indices(m, start_date, iterate_quarters(start_date, T - 1)), i_zlb_start, false
            else
                # Populate vector of regime indices
                if isnothing(i_splice_zlb) # post-ZLB is the last regime
                    for reg in 1:(length(n_regime_periods) - 1)
                        regime_inds[reg] = (n_regime_periods[reg] + 1):n_regime_periods[reg + 1]
                    end
                    regime_inds[end - 1] = (n_regime_periods[end] + 1):n_nozlb_periods
                    regime_inds[end]     = (n_nozlb_periods + 1):T

                    i_zlb_start = length(regime_inds)
                elseif i_splice_zlb > 2 # at least one full regime before ZLB starts
                    regime_inds[1] = 1:n_regime_periods[2]
                    for reg in 2:i_splice_zlb - 2 # if i_splice_zlb == 3, then this loop does not run
                        regime_inds[reg] = (n_regime_periods[reg] + 1):n_regime_periods[reg + 1]
                    end
                    regime_inds[i_splice_zlb - 1] = (n_regime_periods[i_splice_zlb - 1] + 1):n_nozlb_periods
                    regime_inds[i_splice_zlb]     = (n_nozlb_periods + 1):(n_regime_periods[i_splice_zlb])

                    i_zlb_start = i_splice_zlb

                    # Index reg + 1 b/c we have spliced pre- and post- ZLB regime in
                    for reg in i_splice_zlb:(length(n_regime_periods) - 1)
                        regime_inds[reg + 1] = (n_regime_periods[reg] + 1):n_regime_periods[reg + 1]
                    end
                    regime_inds[end] = (n_regime_periods[end] + 1):T
                else # post- ZLB regimes start in the very first regime specified by user
                    regime_inds[1] = 1:n_nozlb_periods
                    regime_inds[2] = (n_nozlb_periods + 1):n_regime_periods[2]

                    i_zlb_start = 2

                    # Index reg + 1 b/c spliced pre- and post-ZLB regime in
                    for reg in 2:(length(n_regime_periods) - 1)
                        regime_inds[reg + 1] = (n_regime_periods[reg] + 1):n_regime_periods[reg + 1]
                    end
                    regime_inds[end] = (n_regime_periods[end] + 1):T
                end
                splice_zlb_regime = true
            end
        else
            # This is the case that date_zlb_start <= start_date so the first regime is the post-ZLB regime (no pre-ZLB)
            # Since there is no ZLB regime to splice in, we return the regime_inds corresponding to the user-specified regimes
            regime_inds = Vector{UnitRange{Int64}}(undef, length(n_regime_periods))
            for reg in 1:(length(n_regime_periods) - 1)
                regime_inds[reg] = (n_regime_periods[reg] + 1):n_regime_periods[reg + 1]
            end
            regime_inds[end] = (n_regime_periods[end] + 1):T

            i_zlb_start = 1 # start immediately in post-ZLB

            splice_zlb_regime = date_zlb_start(m) == start_date ? false : true
        end
    else # Empty, so we ignore regime switching
        regime_inds = UnitRange{Int64}[1:T]
        i_zlb_start = 0
        splice_zlb_regime = false
    end
    # Remove regimes that only are for the forecast period (i.e. bigger than number of rows in dataframe)
    filter!(x -> last(x) <= T, regime_inds)

    # Partition sample into regimes, including the pre- and post-ZLB regimes.
    # regime_inds, i_zlb_start, splice_zlb_regime = zlb_plus_regime_indices(m, data, start_date)

    # Get system matrices for each regime
    n_reg = length(regime_inds)
    ind_zlb_start = i_zlb_start

    S = Float64

    if n_reg <= 0 || ind_zlb_start <= 0
        throw(DomainError())
    end

    if start_date < date_presample_start(m)
        error("Start date $start_date must be >= date_presample_start(m)")

        # TODO: This technically doesn't handle the case where the end_date of the sample
        # is before the start of the ZLB
    elseif date_presample_start(m) <= start_date <= date_zlb_start(m)
        TTTs = Vector{Matrix{S}}(undef, n_reg)
        RRRs = Vector{Matrix{S}}(undef, n_reg)
        CCCs = Vector{Vector{S}}(undef, n_reg)
        ZZs  = Vector{Matrix{S}}(undef, n_reg)
        DDs  = Vector{Vector{S}}(undef, n_reg)
        QQs  = Vector{Matrix{S}}(undef, n_reg)
        EEs  = Vector{Matrix{S}}(undef, n_reg)

        shock_inds = inds_shocks_no_ant(m)

        if splice_zlb_regime
            if ind_zlb_start == 2
                if false#n_mon_anticipated_shocks(m) == 0
                    QQ_preZLB = copy(system[1][:QQ]) # anticipated shocks are zero already
                else
                    QQ_preZLB                         = zeros(size(system[1][:QQ]))
                    QQ_preZLB[shock_inds, shock_inds] = system[1][:QQ][shock_inds, shock_inds]
                end
                QQ_postZLB = copy(system[1][:QQ])

                TTTs[1:2] = [system[1][:TTT], system[1][:TTT]]::Vector{Matrix{S}}
                RRRs[1:2] = [system[1][:RRR], system[1][:RRR]]::Vector{Matrix{S}}
                CCCs[1:2] = [system[1][:CCC], system[1][:CCC]]::Vector{Vector{S}}
                ZZs[1:2]  = [system[1][:ZZ],  system[1][:ZZ]]::Vector{Matrix{S}}
                DDs[1:2]  = [system[1][:DD],  system[1][:DD]]::Vector{Vector{S}}
                EEs[1:2]  = [system[1][:EE],  system[1][:EE]]::Vector{Matrix{S}}
                QQs[1]    = QQ_preZLB
                QQs[2]    = QQ_postZLB

                if n_reg >= 3
                    input = 2:(n_reg - 1) # minus 1 b/c spliced in a ZLB regime

                    TTTs[3:end] = [system[i][:TTT] for i in input]::Vector{Matrix{S}}
                    RRRs[3:end] = [system[i][:RRR] for i in input]::Vector{Matrix{S}}
                    CCCs[3:end] = [system[i][:CCC] for i in input]::Vector{Vector{S}}
                    ZZs[3:end]  = [system[i][:ZZ] for i in input]::Vector{Matrix{S}}
                    DDs[3:end]  = [system[i][:DD] for i in input]::Vector{Vector{S}}
                    QQs[3:end]  = [system[i][:QQ] for i in input]::Vector{Matrix{S}}
                    EEs[3:end]  = [system[i][:EE] for i in input]::Vector{Matrix{S}}
                end
            elseif ind_zlb_start >= 3
                # Populate matrices for regimes before the ZLB split
                pre_zlb_input = 1:(ind_zlb_start - 2) # minus 2 b/c ind_zlb_start - 1 will be the pre-ZLB regime
                splice_reg    = ind_zlb_start - 1

                TTTs[pre_zlb_input] = [system[i][:TTT] for i in pre_zlb_input]::Vector{Matrix{S}}
                RRRs[pre_zlb_input] = [system[i][:RRR] for i in pre_zlb_input]::Vector{Matrix{S}}
                CCCs[pre_zlb_input] = [system[i][:CCC] for i in pre_zlb_input]::Vector{Vector{S}}
                ZZs[pre_zlb_input]  = [system[i][:ZZ] for i in pre_zlb_input]::Vector{Matrix{S}}
                DDs[pre_zlb_input]  = [system[i][:DD] for i in pre_zlb_input]::Vector{Vector{S}}
                EEs[pre_zlb_input]  = [system[i][:EE] for i in pre_zlb_input]::Vector{Matrix{S}}

                if false#n_mon_anticipated_shocks(m) == 0
                    QQs[pre_zlb_input] = [system[i][:QQ] for i in pre_zlb_input]::Vector{Matrix{S}}
                else
                    for i in pre_zlb_input
                        QQ_preZLB                         = zeros(size(system[i][:QQ]))
                        QQ_preZLB[shock_inds, shock_inds] = system[i][:QQ][shock_inds, shock_inds]
                        QQs[i]                            = QQ_preZLB
                    end
                end

                # Regime in which the ZLB occurs
                if false#n_mon_anticipated_shocks(m) == 0
                    QQ_preZLB = copy(system[splice_reg, :QQ]) # anticipated shocks are zero already
                else
                    QQ_preZLB                         = zeros(size(system[splice_reg][:QQ]))
                    QQ_preZLB[shock_inds, shock_inds] = system[splice_reg][:QQ][shock_inds, shock_inds]
                end
                QQ_postZLB = copy(system[splice_reg][:QQ])
                TTTs[splice_reg:ind_zlb_start] = [system[splice_reg][:TTT], system[splice_reg][:TTT]]::Vector{Matrix{S}}
                RRRs[splice_reg:ind_zlb_start] = [system[splice_reg][:RRR], system[splice_reg][:RRR]]::Vector{Matrix{S}}
                CCCs[splice_reg:ind_zlb_start] = [system[splice_reg][:CCC], system[splice_reg][:CCC]]::Vector{Vector{S}}
                ZZs[splice_reg:ind_zlb_start]  = [system[splice_reg][:ZZ],  system[splice_reg][:ZZ]]::Vector{Matrix{S}}
                DDs[splice_reg:ind_zlb_start]  = [system[splice_reg][:DD],  system[splice_reg][:DD]]::Vector{Vector{S}}
                EEs[splice_reg:ind_zlb_start]  = [system[splice_reg][:EE],  system[splice_reg][:EE]]::Vector{Matrix{S}}
                QQs[splice_reg]     = QQ_preZLB
                QQs[ind_zlb_start] = QQ_postZLB

                # Handle regimes after the ZLB split
                if n_reg >= ind_zlb_start + 1
                    post_zlb_input = ind_zlb_start:(n_reg - 1) # minus 1 to both b/c spliced in the ZLB regime

                    TTTs[4:end] = [system[i][:TTT] for i in post_zlb_input]::Vector{Matrix{S}}
                    RRRs[4:end] = [system[i][:RRR] for i in post_zlb_input]::Vector{Matrix{S}}
                    CCCs[4:end] = [system[i][:CCC] for i in post_zlb_input]::Vector{Vector{S}}
                    ZZs[4:end]  = [system[i][:ZZ] for i in post_zlb_input]::Vector{Matrix{S}}
                    DDs[4:end]  = [system[i][:DD] for i in post_zlb_input]::Vector{Vector{S}}
                    QQs[4:end]  = [system[i][:QQ] for i in post_zlb_input]::Vector{Matrix{S}}
                    EEs[4:end]  = [system[i][:EE] for i in post_zlb_input]::Vector{Matrix{S}}
                end
             end
        else
            if ind_zlb_start > 1
                pre_zlb_input = 1:(ind_zlb_start - 1) # minus 1 b/c post-ZLB coincides with a regime

                TTTs[pre_zlb_input] = [system[i][:TTT] for i in pre_zlb_input]::Vector{Matrix{S}}
                RRRs[pre_zlb_input] = [system[i][:RRR] for i in pre_zlb_input]::Vector{Matrix{S}}
                CCCs[pre_zlb_input] = [system[i][:CCC] for i in pre_zlb_input]::Vector{Vector{S}}
                ZZs[pre_zlb_input]  = [system[i][:ZZ] for i in pre_zlb_input]::Vector{Matrix{S}}
                DDs[pre_zlb_input]  = [system[i][:DD] for i in pre_zlb_input]::Vector{Vector{S}}
                EEs[pre_zlb_input]  = [system[i][:EE] for i in pre_zlb_input]::Vector{Matrix{S}}

                if false#n_mon_anticipated_shocks(m) == 0
                    QQs[pre_zlb_input] = Matrix{S}[system[i][:QQ] for i in pre_zlb_input]::Vector{Matrix{S}}
                else
                    for i in pre_zlb_input
                        QQ_preZLB                         = zeros(size(system[i][:QQ]))
                        QQ_preZLB[shock_inds, shock_inds] = system[i][:QQ][shock_inds, shock_inds]
                        QQs[i]                            = QQ_preZLB
                    end
                end
            end

            # Handle regimes after the ZLB split
            post_zlb_input = ind_zlb_start:n_reg # no minus 1 b/c didn't splice in the regime

            TTTs[post_zlb_input] = [system[i][:TTT] for i in post_zlb_input]::Vector{Matrix{S}}
            RRRs[post_zlb_input] = [system[i][:RRR] for i in post_zlb_input]::Vector{Matrix{S}}
            CCCs[post_zlb_input] = [system[i][:CCC] for i in post_zlb_input]::Vector{Vector{S}}
            ZZs[post_zlb_input]  = [system[i][:ZZ] for i in post_zlb_input]::Vector{Matrix{S}}
            DDs[post_zlb_input]  = [system[i][:DD] for i in post_zlb_input]::Vector{Vector{S}}
            QQs[post_zlb_input]  = [system[i][:QQ] for i in post_zlb_input]::Vector{Matrix{S}}
            EEs[post_zlb_input]  = [system[i][:EE] for i in post_zlb_input]::Vector{Matrix{S}}
         end
    elseif date_zlb_start(m) < start_date
        # We always use post-ZLB matrices in this case,
        # hence we just return the matrices of the `RegimeSwitchingSystem`
        input = 1:n_reg
        TTTs = [system[i][:TTT] for i in input]::Vector{Matrix{S}}
        RRRs = [system[i][:RRR] for i in input]::Vector{Matrix{S}}
        CCCs = [system[i][:CCC] for i in input]::Vector{Vector{S}}
        ZZs  = [system[i][:ZZ] for i in input]::Vector{Matrix{S}}
        DDs  = [system[i][:DD] for i in input]::Vector{Vector{S}}
        QQs  = [system[i][:QQ] for i in input]::Vector{Matrix{S}}
        EEs  = [system[i][:EE] for i in input]::Vector{Matrix{S}}
    end
#=
    TTTs, RRRs, CCCs, QQs, ZZs, DDs, EEs = zlb_plus_regime_matrices(m, system, length(regime_inds),
                                                                    start_date;
                                                                    ind_zlb_start = i_zlb_start,
                                                                    splice_zlb_regime = splice_zlb_regime)
=#
    # If s_0 and P_0 provided, check that rows and columns corresponding to
    # anticipated shocks are zero in P_0
    if !isempty(s_0) && !isempty(P_0)
        ant_state_inds = setdiff(1:n_states_augmented(m), inds_states_no_ant(m))
        @assert all(x -> x == 0, P_0[:, ant_state_inds])
        @assert all(x -> x == 0, P_0[ant_state_inds, :])
    end

    # Specify number of presample periods if we don't want to include them in
    # the final results
    Nt0 = include_presample ? 0 : n_presample_periods(m)

    # Run Kalman filter, construct Kalman object, and return
    out = kalman_filter(regime_inds, data, TTTs, RRRs, CCCs, QQs,
                        ZZs, DDs, EEs, s_0, P_0; outputs = outputs,
                        Nt0 = Nt0, tol = tol)
return out[1]
end

# Initialize model objects
function model_init(;subspec="ss10", run_stored = false)
m = Model1002(subspec)
fp = dirname(@__FILE__)

if subspec == "ss22" || subspec == "ss24"
    standard_spec!(m, "191118", fcast_date = Date(2019, 12, 31))
    m <= Setting(:add_laborshare_measurement, false)
    m <= Setting(:hours_first_observable, false)

    m <= Setting(:date_regime2_start, Date(1990, 3, 31))
    m <= Setting(:date_regime2_start_text, "900331", true, "reg2start", "The text version to be saved of when regime 2 starts")
end

m <= Setting(:data_vintage, "191118")
m <= Setting(:saveroot, joinpath(fp, "output_data/"))
m <= Setting(:dataroot, joinpath(fp, "input_data/"))
m <= Setting(:date_forecast_start, Date(2019,12,31))

df = load_data(m, check_empty_columns = false)

m <= Setting(:n_particles, 500)#15000)
m <= Setting(:n_smc_blocks, 5)
m <= Setting(:n_mh_steps_smc, 1)
m <= Setting(:use_parallel_workers, false)
m <= Setting(:step_size_smc, 0.5)
m <= Setting(:resampler_smc, :systematic)
m <= Setting(:target_accept, 0.25)
m <= Setting(:mixture_proportion, 0.9)
m <= Setting(:use_fixed_schedule, false)
m <= Setting(:adaptive_tempering_target_smc, 0.97)
m <= Setting(:resampling_threshold, .5)
m <= Setting(:use_fixed_schedule, false)

# df_first = df[Date(1964, 3, 1) .<= df[:date]  .<= Date(1989, 12, 31), :]
# df_second = df[Date(1990, 3, 1) .<= df[:date], :]

# my_procs = addprocs_frbny(200)
# @everywhere using DSGE, OrderedCollections, DSGEModels, ModelConstructors
# @everywhere include("../../../src/estimate/param_regimes.jl")

# COMBINED
if run_stored
    m <= Setting(:period, "combined", true, "period", "period for this exercse (first or second)")
end
m <= Setting(:date_presample_start, DSGE.quartertodate("1964-Q1")) # need to reset presample
m <= Setting(:date_mainsample_start, DSGE.quartertodate("1964-Q3"))

# Set model regimes
# m <= Setting(:regime_switching, true)
# m <= Setting(:n_regimes, 2)
if !run_stored || subspec in ["ss22", "ss24"]
    m <= Setting(:regime_switching, true)
    m <= Setting(:n_regimes, 2)
    m <= Setting(:regime_dates, Dict(1 => Date(1964,3,31), 2 => Date(1990, 3, 31)))
end
if !run_stored
    setup_regime_switching_inds!(m)
end
return m, df
end

if switch_all

if run_stored
    m, df = model_init(run_stored = run_stored)

    cloud10_reg1 = JLD.load("m1002/ss10/estimate/raw/smc_cloud_npart=15000_period=r1_preZLB=false_reg2start=900331_vint=191118.jld2", "cloud").particles
    cloud10_reg2 = JLD.load("m1002/ss10/estimate/raw/smc_cloud_npart=15000_period=r2_preZLB=false_reg2start=900331_vint=191118.jld2", "cloud").particles

    inds = argmax(SMC.get_loglh(cloud10_reg1))
    para1 = cloud10_reg1[inds, 1:SMC.ind_para_end(size(cloud10_reg1, 2))]
    para2 = cloud10_reg2[inds, 1:SMC.ind_para_end(size(cloud10_reg2, 2))]

    for para in 1:length(m.parameters)
        m.parameters[para].value = para1[para]
    end

    df = df[Date(1964, 3, 1) .<= df[:date], :]
    df_first = df[Date(1964, 3, 1) .<= df[:date]  .<= Date(1989, 12, 31), :]
    df_second = df[Date(1990, 3, 1) .<= df[:date], :]
    system = compute_system(m)
    kal = DSGE.filter(m, df_first, system)
    loglh1 = kal[:loglh]
    # loglh1 = DSGE.forecast_one_draw(m, :mode, :none, [:forecastobs], para1, df_first,
    #                       return_loglh = true)

    old_system1 = compute_system(m)

    m <= Setting(:period, "second", true, "period", "period for this exercse (first or second)")
    m <= Setting(:date_presample_start, DSGE.quartertodate("1990-Q1")) # need to reset presample
    m <= Setting(:date_mainsample_start, DSGE.quartertodate("1990-Q3"))

    for para in 1:length(m.parameters)
        m.parameters[para].value = para2[para]
    end
    eq_old2 = eqcond(m)
    old_system2_1 = compute_system(m)

    m2_1 = deepcopy(m)

    system = compute_system(m)
    kal = DSGE.filter(m, df_second, system)
    loglh2 = kal[:loglh]

    # loglh2 = DSGE.forecast_one_draw(m, :mode, :none, [:forecastobs], para2, df_second,
    #                       return_loglh = true)
    m2 = deepcopy(m)
    old_system2 = compute_system(m)
end
m, df = model_init(subspec = "ss10")
# Add map from model regimes to parameter regimes
param_mat = repeat([1 2], length(m.parameters))
setup_param_regimes!(m, param_mat = param_mat)

# Set parameter values equal to Brookings cloud
cloud10_reg1 = JLD.load("m1002/ss10/estimate/raw/smc_cloud_npart=15000_period=r1_preZLB=false_reg2start=900331_vint=191118.jld2", "cloud").particles
cloud10_reg2 = JLD.load("m1002/ss10/estimate/raw/smc_cloud_npart=15000_period=r2_preZLB=false_reg2start=900331_vint=191118.jld2", "cloud").particles

inds = argmax(SMC.get_loglh(cloud10_reg1))
para1 = cloud10_reg1[inds, 1:SMC.ind_para_end(size(cloud10_reg1, 2))]
para2 = cloud10_reg2[inds, 1:SMC.ind_para_end(size(cloud10_reg2, 2))]

if !run_stored
    loglh1 = SMC.get_loglh(cloud10_reg1)[inds]
    loglh2 = SMC.get_loglh(cloud10_reg2)[inds]
end
# @show loglh1, loglh2

# create parameter regime values (initialized as
## being the same in both regimes)
paramed = copy(para1)
for p2 in 1:length(m.parameters)
    p = m.parameters[p2]
    for i in 1:2
        ModelConstructors.set_regime_prior!(p, i, p.prior)
        ModelConstructors.set_regime_fixed!(p, i, p.fixed)
        ModelConstructors.set_regime_valuebounds!(p, i, p.valuebounds)
    end
    ModelConstructors.set_regime_val!(p, 1, para1[p2])
    ModelConstructors.set_regime_val!(p, 2, para2[p2])

    append!(paramed, para2[p2])
end

new_system_1 = compute_system(m)
m_new = deepcopy(m)
# Call function forecast_one_draw with return_loglh = true
kal = DSGE.filter(m, df, new_system_1)
loglh = kal[:loglh]
# loglh = DSGE.forecast_one_draw(m, :mode, :none, [:forecastobs], paramed, df,
#                           return_loglh = true)

new_system = compute_system(m)

# Compare this value with loglh from cloud for that parameter
@show sum(loglh[1:104]) #≈ loglh1
@show sum(loglh[105:end]) #≈ loglh2: Earlier Skipping 1990-Q1 and Q2 as before mainsample
end

if switch_policy
#################################
# Only policy rule param change #
#################################
# Set parameter values equal to Brookings cloud
cloud24_reg = JLD.load("m1002/ss24/estimate/raw/smc_cloud_friday=true_npart=20000_period=full_incZLB_reg2start=900331_reg=2_vint=191118.jld2", "cloud").particles

inds = argmax(SMC.get_loglh(cloud24_reg))
para1 = cloud24_reg[inds, 1:SMC.ind_para_end(size(cloud24_reg, 2))]

if run_stored
m, df = model_init(subspec = "ss24", run_stored = true)

for para in 1:length(m.parameters)
    m.parameters[para].value = para1[para]
end

system = compute_system(m, regime_switching = true, n_regimes = 2)
loglh1 = kalman_stored(m,df,system)

#kal = filter(m, df, system)
#loglh1 = kal[:loglh]

JLD2.jldopen("output_data/loglh_policy.jld2", "w") do file
    write(file, "loglh", loglh1)  # alternatively, say "@write file A"
    write(file, "para", m.parameters)
end
else
    m, df = model_init(subspec = "ss10", run_stored = false)
    m <= Setting(:use_population_forecast, true)
#=
    m <= parameter(:ψ1, 1.3679, (1e-5, 10.), (1e-5, 10.00), ModelConstructors.Exponential(), Normal(1.5, 0.25), fixed=false,
                   description="ψ₁: Weight on inflation gap in monetary policy rule.",
                   tex_label="\\psi_1")

    m <= parameter(:ψ2, 0.0388, (-0.5, 0.5), (-0.5, 0.5), ModelConstructors.Untransformed(), Normal(0.12, 0.05), fixed=false,
                   description="ψ₂: Weight on output gap in monetary policy rule.",
                   tex_label="\\psi_2")

    m <= parameter(:ψ3, 0.2464, (-0.5, 0.5), (-0.5, 0.5), ModelConstructors.Untransformed(), Normal(0.12, 0.05), fixed=false,
                   description="ψ₃: Weight on rate of change of output gap in the monetary policy rule.",
                   tex_label="\\psi_3")

    m <= parameter(:ρ, 0.7126, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.75, 0.10), fixed=false,
                   description="ρ: The degree of inertia in the monetary policy rule.",
                   tex_label="\\rho_R")
=#
# create parameter regime values (initialized as
## being the same in both regimes)
Ψ_ind = 15#(1:length(m.parameters))[[m.parameters[i].key == :ψ1 for i in 1:length(m.parameters)]]
Ψ2_ind = 16
ψ3_ind = 17
ρ_ind = 20
ρ2_ind = 21
Ψ_r2_ind = 22
Ψ2_r2_ind = 23
Ψ3_r2_ind = 24

# Add map from model regimes to parameter regimes
param_mat = repeat([1 1], length(m.parameters))
for k in vcat(15:17,20)
    param_mat[k,:] = [1 2]
end
setup_param_regimes!(m, param_mat = param_mat)

# para2 = cloud10_reg2[inds, 1:SMC.ind_para_end(size(cloud10_reg2, 2))]

# loglh1 = SMC.get_loglh(cloud24_reg)[inds]
# loglh2 = SMC.get_loglh(cloud10_reg2)[inds]

# @show loglh1, loglh2

params = para1[vcat(1:20, 25:length(para1))]
for p2 in 1:length(para1[vcat(1:20, 25:length(para1))])
    p = m.parameters[p2]
    end_reg = p2 in [15,16,17,20] ? 2 : 1
    for i in 1:end_reg
        ModelConstructors.set_regime_prior!(p, i, p.prior)
        ModelConstructors.set_regime_fixed!(p, i, p.fixed)
        ModelConstructors.set_regime_valuebounds!(p, i, p.valuebounds)
    end
    if p2 > 20
        ModelConstructors.set_regime_val!(p, 1, para1[p2+4])
    else
        ModelConstructors.set_regime_val!(p, 1, para1[p2])
    end
    if end_reg == 1
        continue
    elseif p2 < 20
        ModelConstructors.set_regime_val!(p, 2, para1[p2+7])
        append!(params, para1[p2+7])
    else
        ModelConstructors.set_regime_val!(p, 2, para1[p2+1])
        append!(params, para1[p2+1])
    end
end

# Call function forecast_one_draw with return_loglh = true
system = compute_system(m)
kal = DSGE.filter(m, df, system)
loglh = kal[:loglh]
#=
loglh = DSGE.forecast_one_draw(m, :mode, :none, [:forecastobs], params, df,
                          return_loglh = true)
=#
# Compare this value with loglh from cloud for that parameter
@show sum(loglh)
# @show sum(loglh[1:104]) #≈ loglh1
# @show sum(loglh[105:end]) #≈ loglh2
end
end

if switch_kappa
#################################
### Only kappa_p param change ###
#################################
# Set parameter values equal to Brookings cloud
cloud2_reg = JLD.load("m1002/ss22/estimate/raw/smc_cloud_friday=true_npart=20000_period=full_incZLB_reg2start=900331_reg=2_vint=191118.jld2", "cloud").particles

inds = argmax(SMC.get_loglh(cloud2_reg))
para1 = cloud2_reg[inds, 1:SMC.ind_para_end(size(cloud2_reg, 2))]

if run_stored
m, df = model_init(subspec = "ss22", run_stored = true)

for para in 1:length(m.parameters)
    m.parameters[para].value = para1[para]
end

system = compute_system(m, regime_switching = true, n_regimes = 2)
loglh1 = kalman_stored(m, df, system)
# kal = DSGE.filter(m, df, system)
# loglh1 = kal[:loglh]

JLD2.jldopen("output_data/loglh_kappa.jld2", "w") do file
    write(file, "loglh", loglh1)  # alternatively, say "@write file A"
    write(file, "para", m.parameters)
end
else

m, df = model_init(subspec = "ss10")

m <= parameter(:ζ_p, 0.8940, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1), fixed=false,
                   description="ζ_p: The Calvo parameter (Regime 1). In every period, intermediate goods producers optimize prices with probability (1-ζ_p). With probability ζ_p, prices are adjusted according to a weighted average of the previous period's inflation (π_t1) and steady-state inflation (π_star).",
                   tex_label="\\zeta_p")

# create parameter regime values (initialized as
## being the same in both regimes)
ζ_p_ind = 2

# Add map from model regimes to parameter regimes
param_mat = repeat([1 1], length(m.parameters))
for k in [2]
    param_mat[k,:] = [1 2]
end
setup_param_regimes!(m, param_mat = param_mat)

# para2 = cloud10_reg2[inds, 1:SMC.ind_para_end(size(cloud10_reg2, 2))]

# loglh1 = SMC.get_loglh(cloud2_reg)[inds]
# loglh2 = SMC.get_loglh(cloud10_reg2)[inds]

# @show loglh1, loglh2

params = para1[vcat(1:20, 22:length(para1))]
for p2 in 1:length(para1[vcat(1:20, 22:length(para1))])
    p = m.parameters[p2]
    end_reg = p2 in [2] ? 2 : 1
    for i in 1:end_reg
        ModelConstructors.set_regime_prior!(p, i, p.prior)
        ModelConstructors.set_regime_fixed!(p, i, p.fixed)
        ModelConstructors.set_regime_valuebounds!(p, i, p.valuebounds)
    end
    if p2 > 20
        ModelConstructors.set_regime_val!(p, 1, para1[p2+1])
    else
        ModelConstructors.set_regime_val!(p, 1, para1[p2])
    end
    if end_reg == 1
        continue
    else
        ModelConstructors.set_regime_val!(p, 2, para1[21])
        append!(params, para1[21])
    end
end

# Call function forecast_one_draw with return_loglh = true
system = compute_system(m)
kal = DSGE.filter(m, df, system)
loglh_kappa = kal[:loglh]
#=
loglh = DSGE.forecast_one_draw(m, :mode, :none, [:forecastobs], params, df,
                          return_loglh = true)
=#
# Compare this value with loglh from cloud for that parameter
@show sum(loglh_kappa)
# @show sum(loglh[1:104]) #≈ loglh1
# @show sum(loglh[105:end]) #≈ loglh2
end
end

#=
@show rawpath(m, "estimate", "smc_cloud", ["new_brookings"])
smc2(m, df_to_matrix(m, df), filestring_addl = ["new_brookings_fixed_sched"], regime_switching = true)

# FIRST
m <= Setting(:period, "first", true, "period", "period for this exercse (first or second)")
m <= Setting(:date_presample_start, DSGE.quartertodate("1964-Q1")) # need to reset presample
m <= Setting(:date_mainsample_start, DSGE.quartertodate("1964-Q3"))
smc2(m, df_to_matrix(m, df_first), filestring_addl = ["old_brookings_first"])

# SECOND
m <= Setting(:period, "second", true, "period", "period for this exercse (first or second)")
m <= Setting(:date_presample_start, DSGE.quartertodate("1990-Q1")) # need to reset presample
m <= Setting(:date_mainsample_start, DSGE.quartertodate("1990-Q3"))
smc2(m, df_to_matrix(m, df_second), filestring_addl = ["old_brookings_second"])
=#
# rmprocs(my_procs)

ext: jl


text:
using AdvancedPS
using AbstractMCMC
using Distributions
using Libtask
using Random
using Test

@testset "AdvancedPS.jl" begin
    @testset "Resampling tests" begin
        include("resampling.jl")
    end
    @testset "Container tests" begin
        include("container.jl")
    end
    @testset "SMC and PG tests" begin
        include("smc.jl")
    end
    @testset "RNG tests" begin
        include("rng.jl")
    end
end


text:
abstract type AbstractUninitializedReference{T} end
"""
This is a struct wrapping a RefValue rather than a mutable struct
because the ReverseDiff API wants the initialized parameters to be
Refs of some kind; it is simpler to wrap the first access in an
Uninitialized struct than a mutable struct.
"""
struct UninitializedRef{T} <: AbstractUninitializedReference{T}
    data::Base.RefValue{T}
end
@inline Base.getindex(uninit::UninitializedRef) = uninit.data[]
@inline Base.setindex!(uninit::UninitializedRef{T}, v::T) where {T} = uninit.data[] = v
struct UninitializedPtr{T} <: AbstractUninitializedReference{T}
    ptr::Ptr{T}
end
@inline Base.getindex(uninit::UninitializedPtr) = VectorizationBase.load(uninit.ptr)
@inline Base.setindex!(uninit::UninitializedPtr{T}, v::T) where {T} = VectorizationBase.store!(uninit.ptr, v)

function ∂getindex end
function ∂mul end
@inline alloc_adjoint(x::T) where {T<:Real} = Ref{T}()
@inline uninitialized(::Nothing) = nothing
@inline uninitialized(x::Base.RefValue) = UninitializedRef(x)
@inline function uninitialized(x::Pointer)
    # ptr = pointer(x)
    # SIMDPirates.lifetime_start!(ptr)
    UninitializedPtr(pointer(x))#ptr)
end
@inline initialized(x) = x
@inline initialized(x::UninitializedRef) = x.data
@inline initialized(x::UninitializedPtr) = Pointer(x.ptr)
isinitialized(::Any) = true
isinitialized(::Type{<:UninitializedRef}) = false
isinitialized(::Type{<:UninitializedPtr}) = false
@inline zero_initialize!(A::AbstractArray{T}) where {T} = fill!(A, zero(T))
@inline zero_initialize!(x::Base.RefValue{T}) where {T} = x[] = zero(T)
@inline zero_initialize!(x::Pointer{T}) where {T} = x[] = zero(T)

@inline alloc_adjoint(x::T) where {T} = Ref{T}()
@inline alloc_adjoint(sp::StackPointer, A::AbstractArray) = similar(sp, A)


text:
# Gauss-Kronrod quadrature rule, via the QuadGK package, since Genz-Malik
# rule does not handle the 1d case.  We will just use a fixed-order (7) G-K rule.

struct GaussKronrod{T<:AbstractFloat}
    x::Vector{T}
    w::Vector{T}
    wg::Vector{T}
end

# cache the Gauss-Kronrod rules so that we don't
# call QuadGK.kronrod every time.
const gkcache = Dict{Type, GaussKronrod}()

function GaussKronrod(::Type{T}) where {T<:AbstractFloat}
    haskey(gkcache, T) && return gkcache[T]::GaussKronrod{T}
    gkcache[T] = g = GaussKronrod{T}(QuadGK.kronrod(T,7)...)
    return g
end

# further speed up the common case of double precision (25% faster for a trivial integrand)
const gk_float64 = GaussKronrod{Float64}(QuadGK.kronrod(Float64,7)...)
GaussKronrod(::Type{Float64}) = gk_float64

countevals(g::GaussKronrod) = 17

function (g::GaussKronrod{T})(f, a_::SVector{1}, b_::SVector{1}, norm=vecnorm) where {T}
    a = a_[1]
    b = b_[1]
    c = (a+b)*T(0.5)
    Δ = (b-a)*T(0.5)

    fx⁰ = f(SVector(c))                # f(0)
    I = fx⁰ * g.w[end]
    I′ = fx⁰ * g.wg[end]
    @inbounds for i = 1:length(g.x)-1
        Δx = Δ * g.x[i]
        fx = f(SVector(c + Δx)) + f(SVector(c - Δx))
        I += fx * g.w[i]
        if iseven(i)
            I′ += fx * g.wg[i>>1]
        end
    end
    V = abs(Δ)
    I *= V
    I′ *= V
    return I, norm(I - I′), 1
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
    NoIntegerAlgorithm

Empty functor object for running an L-shaped algorithm without dealing with integer variables.

"""
struct NoIntegerAlgorithm <: AbstractIntegerAlgorithm
    skip_check::Bool
end

function initialize_integer_algorithm!(integer::NoIntegerAlgorithm, first_stage::JuMP.Model)
    # Sanity check
    if !integer.skip_check &&
       any(is_binary, all_decision_variables(first_stage, 1)) ||
       any(is_integer, all_decision_variables(first_stage, 1))
        @warn "First-stage has binary/integer decisions and no `IntegerStrategy` has been set. Procedure will fail if second-stage has binary/integer variables. Otherwise, the master_optimizer must be integer-capable."
    end
    return nothing
end

function initialize_integer_algorithm!(integer::NoIntegerAlgorithm, subproblem::SubProblem)
    if !integer.skip_check &&
       any(is_binary, all_decision_variables(subproblem.model, StochasticPrograms.stage(subproblem.model))) ||
       any(is_integer, all_decision_variables(subproblem.model, StochasticPrograms.stage(subproblem.model))) ||
       any(is_binary, all_auxiliary_variables(subproblem.model)) ||
       any(is_integer, all_auxiliary_variables(subproblem.model))
        error("Second-stage has binary/integer decisions. Rerun procedure with an `IntegerStrategy`.")
    end
    return nothing
end

function handle_integrality!(lshaped::AbstractLShaped, ::NoIntegerAlgorithm)
    # Ensure any binary/integer decisions are rounded
    for (i,dvar) in enumerate(all_decision_variables(lshaped.structure.first_stage, 1))
        if is_binary(dvar) || is_integer(dvar)
            lshaped.x[i] = round(lshaped.x[i])
        end
    end
    return nothing
end

function integer_variables(::NoIntegerAlgorithm)
    return MOI.VariableIndex[]
end

function check_optimality(::AbstractLShaped, ::NoIntegerAlgorithm)
    return true
end

function solve_subproblem(subproblem::SubProblem,
                          metadata,
                          ::NoFeasibilityAlgorithm,
                          ::NoIntegerAlgorithm,
                          x::AbstractVector)
    return solve_subproblem(subproblem, x)
end

function solve_subproblem(subproblem::SubProblem,
                          metadata,
                          feasibility_algorithm::FeasibilityCutsWorker,
                          ::NoIntegerAlgorithm,
                          x::AbstractVector)
    model = subproblem.optimizer
    if !prepared(feasibility_algorithm)
        prepare!(model, feasibility_algorithm)
    else
        activate!(model, feasibility_algorithm)
    end
    # Optimize auxiliary problem
    MOI.optimize!(model)
    # Sanity check that aux problem could be solved
    status = MOI.get(subproblem.optimizer, MOI.TerminationStatus())
    if !(status ∈ AcceptableTermination)
        error("Subproblem $(subproblem.id) was not solved properly during feasibility check, returned status code: $status")
    end
    sense = MOI.get(subproblem.optimizer, MOI.ObjectiveSense())
    correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
    w = correction * MOI.get(model, MOI.ObjectiveValue())
    # Ensure correction is available in master
    set_metadata!(metadata, subproblem.id, :correction, correction)
    # Check feasibility
    if w > sqrt(eps())
        # Subproblem is infeasible, create feasibility cut
        return FeasibilityCut(subproblem, x)
    end
    # Restore subproblem and solve as usual
    deactivate!(model, feasibility_algorithm)
    return solve_subproblem(subproblem, x)
end

"""
    IgnoreIntegers

Factory object for [`NoIntegerAlgorithm`](@ref). Passed by default to `integer_strategy` in `LShaped.Optimizer`.

"""
struct IgnoreIntegers <: AbstractIntegerStrategy
    skip_check::Bool

    IgnoreIntegers(skip_check::Bool) = new(skip_check)
end

IgnoreIntegers(; skip_check = false) = IgnoreIntegers(skip_check)

function master(ignore::IgnoreIntegers, ::Type{T}) where T <: AbstractFloat
    return NoIntegerAlgorithm(ignore.skip_check)
end

function worker(ignore::IgnoreIntegers, ::Type{T}) where T <: AbstractFloat
    return NoIntegerAlgorithm(ignore.skip_check)
end
function worker_type(::IgnoreIntegers)
    return NoIntegerAlgorithm
end


text:
@testset "Triangular" begin

    for f in [LowerTriangular, UpperTriangular]
        let P = 10, rng = MersenneTwister(123456), N = 10

            sc = ()->randn(rng)
            mPP = ()->abs.(randn(rng, P, P))
            dPP = ()->Diagonal(abs.(randn(rng, P)))
            tPP = ()->f(abs.(randn(rng, P, P)))
            m0 = randn!(rng, similar(mPP()))

            # Construction.
            @test check_errs(N, unary_ȲD(f)..., tPP, mPP, mPP)
            @test check_errs(N, unary_ȲD_inplace(f, m0)..., tPP, mPP, mPP)

            # Determinant.
            @test check_errs(N, unary_ȲD(det)..., sc, tPP, tPP)
            @test check_errs(N, unary_ȲD_inplace(det, tPP())..., sc, tPP, tPP)
            @test check_errs(N, unary_ȲD_inplace(det, dPP())..., sc, tPP, tPP)

            # Log Determinant.
            @test check_errs(N, unary_ȲD(logdet)..., sc, tPP, tPP)
            @test check_errs(N, unary_ȲD_inplace(logdet, tPP())..., sc, tPP, tPP)
            @test check_errs(N, unary_ȲD_inplace(logdet, dPP())..., sc, tPP, tPP)
        end
    end
end


text:
"""
Gauss Process
"""
struct GaussProcess{MT, KT, VT}
	mean::MT
	kernel::KT
	lnoise::VT
end
@functor GaussProcess

function Base.show(io::IO, gp::GaussProcess)
	print(io, "GaussProcess(", gp.mean, ", ", gp.kernel, ", ", "lnoise=", gp.lnoise[1], ")")
end

"""
Construct multivariate normal distribution from GP, facilitate sampling
"""
function MvNormal(gp::GaussProcess, x::AbstractArray)
	n = size(x, 2)
	μ = reshape(gp.mean(x), n)
	K = gp.kernel(x) + Diagonal(exp(2*gp.lnoise[1])*ones(n))
	MvNormal(μ, K)
end

"""
loss function, negative log likelihood
"""
# loss function
function negloglik(gp::GaussProcess, x, y; λ=1e-6)
	μ = reshape(gp.mean(x), size(y))
	Σ = gp.kernel(x, λ=λ) + Diagonal(exp(2*gp.lnoise[1])*ones(size(x, 2)))

  d = length(y)
  # L = BackwardsLinalg.cholesky(Σ)
	C = cholesky(Σ)
	L = C.L
	ȳ = y .- μ
	z = L \ ȳ
  0.5*dot(z, z) + (d/2.0)*log(2π) + logdet(L)
end

"""
inference from fitted GP
"""
# inference
function predict(gp::GaussProcess, x, x_old, y_old; λ=1e-6)
	μn = ndims(y_old)>1 ? gp.mean(x) : reshape(gp.mean(x), size(x, 2))
	μo = reshape(gp.mean(x_old), size(y_old))
	Σ_oo = gp.kernel(x_old, λ=λ) + Diagonal(exp(2*gp.lnoise[1])*ones(size(x_old, 2)))
	Σ_no = gp.kernel(x, x_old)
	Σ_nn = gp.kernel(x, λ=λ)

	C = cholesky(Σ_oo)
	L = C.L
	ȳo = y_old .- μo
	zo = L \ ȳo
	α = L' \ zo
	k = L \ Σ_no'
	k̄ = L' \ k

	yn = μn .+ reshape(Σ_no*α, size(μn))
	Σ̄_nn = Σ_nn .- Σ_no*k̄

	yn, diag(Σ̄_nn)
end

"""
dispatch parameters (Array) to a ML model, facilitate use of Optim.jl
"""
function nparams(ps)
	nps = []
	for p in ps
		push!(nps, length(p))
	end
	nps
end

function flatten_params(ps)
	psarray = []
	for p in ps
		push!(psarray, vec(p))
	end
	vcat(psarray...)
end

function dispatch!(model, x::AbstractVector)
	ps = params(model)
	nps = nparams(ps)

	loc = 1
	i = 1
	for p in ps
		p .= reshape(x[loc:loc+nps[i]-1], size(p))
		loc += nps[i]
		i += 1
	end

	model
end

max_issues_repo_name: HamletWantToCode/GPFlux

max_forks_repo_name: HamletWantToCode/GPFlux


text:
###############
# API methods #
###############

"""
    ForwardDiff.derivative(f, x::Real)

Return `df/dx` evaluated at `x`, assuming `f` is called as `f(x)`.

This method assumes that `isa(f(x), Union{Real,AbstractArray})`.
"""
@inline function derivative(f::F, x::R) where {F,R<:Real}
    T = typeof(Tag(f, R))
    return extract_derivative(T, f(Dual{T}(x, one(x))))
end

"""
    ForwardDiff.derivative(f!, y::AbstractArray, x::Real, cfg::DerivativeConfig = DerivativeConfig(f!, y, x), check=Val{true}())

Return `df!/dx` evaluated at `x`, assuming `f!` is called as `f!(y, x)` where the result is
stored in `y`.

Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care.
"""
@inline function derivative(f!, y::AbstractArray, x::Real,
                            cfg::DerivativeConfig{T} = DerivativeConfig(f!, y, x), ::Val{CHK}=Val{true}()) where {T, CHK}
    CHK && checktag(T, f!, x)
    ydual = cfg.duals
    seed!(ydual, y)
    f!(ydual, Dual{T}(x, one(x)))
    map!(value, y, ydual)
    return extract_derivative(T, ydual)
end

"""
    ForwardDiff.derivative!(result::Union{AbstractArray,DiffResult}, f, x::Real)

Compute `df/dx` evaluated at `x` and store the result(s) in `result`, assuming `f` is called
as `f(x)`.

This method assumes that `isa(f(x), Union{Real,AbstractArray})`.
"""
@inline function derivative!(result::Union{AbstractArray,DiffResult},
                             f::F, x::R) where {F,R<:Real}
    T = typeof(Tag(f, R))
    ydual = f(Dual{T}(x, one(x)))
    result = extract_value!(T, result, ydual)
    result = extract_derivative!(T, result, ydual)
    return result
end

"""
    ForwardDiff.derivative!(result::Union{AbstractArray,DiffResult}, f!, y::AbstractArray, x::Real, cfg::DerivativeConfig = DerivativeConfig(f!, y, x), check=Val{true}())

Compute `df!/dx` evaluated at `x` and store the result(s) in `result`, assuming `f!` is
called as `f!(y, x)` where the result is stored in `y`.

Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care.
"""
@inline function derivative!(result::Union{AbstractArray,DiffResult},
                             f!, y::AbstractArray, x::Real,
                             cfg::DerivativeConfig{T} = DerivativeConfig(f!, y, x), ::Val{CHK}=Val{true}()) where {T, CHK}
    CHK && checktag(T, f!, x)
    ydual = cfg.duals
    seed!(ydual, y)
    f!(ydual, Dual{T}(x, one(x)))
    result = extract_value!(T, result, y, ydual)
    result = extract_derivative!(T, result, ydual)
    return result
end

derivative(f, x::AbstractArray) = throw(DimensionMismatch("derivative(f, x) expects that x is a real number. Perhaps you meant gradient(f, x)?"))

#####################
# result extraction #
#####################

# non-mutating #
#--------------#

@inline extract_derivative(::Type{T}, y::Dual) where {T}          = partials(T, y, 1)
@inline extract_derivative(::Type{T}, y::Real) where {T}          = zero(y)
@inline extract_derivative(::Type{T}, y::AbstractArray) where {T} = map(d -> extract_derivative(T,d), y)

# mutating #
#----------#

extract_derivative!(::Type{T}, result::AbstractArray, y::AbstractArray) where {T} =
    map!(d -> extract_derivative(T,d), result, y)
extract_derivative!(::Type{T}, result::DiffResult, y) where {T} =
    DiffResults.derivative!(d -> extract_derivative(T,d), result, y)


text:
@testset "Test Compiler" begin
    @testset "Integer Arithmetic" begin
        for (input, expected_constants, expected_instructions) in [
            (
                "1; 2",
                [1, 2],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpPop),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpPop),
                ],
            ),
            (
                "1 + 2",
                [1, 2],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpAdd),
                    m.make(m.OpPop),
                ],
            ),
            (
                "1 - 2",
                [1, 2],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpSub),
                    m.make(m.OpPop),
                ],
            ),
            (
                "1 * 2",
                [1, 2],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpMul),
                    m.make(m.OpPop),
                ],
            ),
            (
                "1 / 2",
                [1, 2],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpDiv),
                    m.make(m.OpPop),
                ],
            ),
            ("-1", [1], [m.make(m.OpConstant, 0), m.make(m.OpMinus), m.make(m.OpPop)]),
        ]
            run_compiler_tests(input, expected_constants, expected_instructions)
        end
    end

    @testset "Boolean Expressions" begin
        for (input, expected_constants, expected_instructions) in [
            ("true", [], [m.make(m.OpTrue), m.make(m.OpPop)]),
            ("false", [], [m.make(m.OpFalse), m.make(m.OpPop)]),
            (
                "1 == 2",
                [1, 2],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpEqual),
                    m.make(m.OpPop),
                ],
            ),
            (
                "1 != 2",
                [1, 2],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpNotEqual),
                    m.make(m.OpPop),
                ],
            ),
            (
                "1 < 2",
                [1, 2],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpLessThan),
                    m.make(m.OpPop),
                ],
            ),
            (
                "1 > 2",
                [1, 2],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpGreaterThan),
                    m.make(m.OpPop),
                ],
            ),
            (
                "true == false",
                [],
                [m.make(m.OpTrue), m.make(m.OpFalse), m.make(m.OpEqual), m.make(m.OpPop)],
            ),
            (
                "true != false",
                [],
                [
                    m.make(m.OpTrue),
                    m.make(m.OpFalse),
                    m.make(m.OpNotEqual),
                    m.make(m.OpPop),
                ],
            ),
            ("!true", [], [m.make(m.OpTrue), m.make(m.OpBang), m.make(m.OpPop)]),
        ]
            run_compiler_tests(input, expected_constants, expected_instructions)
        end
    end

    @testset "Conditionals" begin
        for (input, expected_constants, expected_instructions) in [
            (
                "if (true) {};",
                [],
                [
                    m.make(m.OpTrue),
                    m.make(m.OpJumpNotTruthy, 8),
                    m.make(m.OpNull),
                    m.make(m.OpJump, 9),
                    m.make(m.OpNull),
                    m.make(m.OpPop),
                ],
            ),
            (
                "if (true) { 10 }; 3333",
                [10, 3333],
                [
                    m.make(m.OpTrue),
                    m.make(m.OpJumpNotTruthy, 10),
                    m.make(m.OpConstant, 0),
                    m.make(m.OpJump, 11),
                    m.make(m.OpNull),
                    m.make(m.OpPop),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpPop),
                ],
            ),
            (
                "if (true) { 10 } else { 20 }; 3333",
                [10, 20, 3333],
                [
                    m.make(m.OpTrue),
                    m.make(m.OpJumpNotTruthy, 10),
                    m.make(m.OpConstant, 0),
                    m.make(m.OpJump, 13),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpPop),
                    m.make(m.OpConstant, 2),
                    m.make(m.OpPop),
                ],
            ),
        ]
            run_compiler_tests(input, expected_constants, expected_instructions)
        end
    end

    @testset "Global Let Statements" begin
        for (input, expected_constants, expected_instructions) in [
            (
                "let one = 1;\nlet two = 2;\n",
                [1, 2],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpSetGlobal, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpSetGlobal, 1),
                ],
            ),
            (
                "let one = 1;\none;\n",
                [1],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpSetGlobal, 0),
                    m.make(m.OpGetGlobal, 0),
                    m.make(m.OpPop),
                ],
            ),
            (
                "let one = 1;\nlet two = one;\ntwo;\n",
                [1],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpSetGlobal, 0),
                    m.make(m.OpGetGlobal, 0),
                    m.make(m.OpSetGlobal, 1),
                    m.make(m.OpGetGlobal, 1),
                    m.make(m.OpPop),
                ],
            ),
        ]
            run_compiler_tests(input, expected_constants, expected_instructions)
        end
    end

    @testset "Global Let Statements" begin
        for (input, expected_constants, expected_instructions) in [
            (
                "let num = 55;\nfn() { num }",
                [55, vcat(m.make(m.OpGetGlobal, 0), m.make(m.OpReturnValue))],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpSetGlobal, 0),
                    m.make(m.OpClosure, 1, 0),
                    m.make(m.OpPop),
                ],
            ),
            (
                "fn() {\nlet num = 55;\nnum\n}",
                [
                    55,
                    vcat(
                        m.make(m.OpConstant, 0),
                        m.make(m.OpSetLocal, 0),
                        m.make(m.OpGetLocal, 0),
                        m.make(m.OpReturnValue),
                    ),
                ],
                [m.make(m.OpClosure, 1, 0), m.make(m.OpPop)],
            ),
            (
                """
                fn() {
                  let a = 55;
                  let b = 77;
                  a + b
                }
                """,
                [
                    55,
                    77,
                    vcat(
                        m.make(m.OpConstant, 0),
                        m.make(m.OpSetLocal, 0),
                        m.make(m.OpConstant, 1),
                        m.make(m.OpSetLocal, 1),
                        m.make(m.OpGetLocal, 0),
                        m.make(m.OpGetLocal, 1),
                        m.make(m.OpAdd),
                        m.make(m.OpReturnValue),
                    ),
                ],
                [m.make(m.OpClosure, 2, 0), m.make(m.OpPop)],
            ),
        ]
            run_compiler_tests(input, expected_constants, expected_instructions)
        end
    end

    @testset "String Expressions" begin
        for (input, expected_constants, expected_instructions) in [
            ("\"monkey\"", ["monkey"], [m.make(m.OpConstant, 0), m.make(m.OpPop)]),
            (
                "\"mon\" + \"key\"",
                ["mon", "key"],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpAdd),
                    m.make(m.OpPop),
                ],
            ),
        ]
            run_compiler_tests(input, expected_constants, expected_instructions)
        end
    end

    @testset "Array Literals" begin
        for (input, expected_constants, expected_instructions) in [
            ("[]", [], [m.make(m.OpArray, 0), m.make(m.OpPop)]),
            (
                "[1, 2, 3]",
                [1, 2, 3],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpConstant, 2),
                    m.make(m.OpArray, 3),
                    m.make(m.OpPop),
                ],
            ),
            (
                "[1 + 2, 3 - 4, 5 * 6]",
                [1, 2, 3, 4, 5, 6],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpAdd),
                    m.make(m.OpConstant, 2),
                    m.make(m.OpConstant, 3),
                    m.make(m.OpSub),
                    m.make(m.OpConstant, 4),
                    m.make(m.OpConstant, 5),
                    m.make(m.OpMul),
                    m.make(m.OpArray, 3),
                    m.make(m.OpPop),
                ],
            ),
        ]
            run_compiler_tests(input, expected_constants, expected_instructions)
        end
    end

    @testset "Hash Literals" begin
        for (input, expected_constants, expected_instructions) in [
            ("{}", [], [m.make(m.OpHash, 0), m.make(m.OpPop)]),
            (
                "{1: 2, 3: 4, 5: 6}",
                [1, 2, 3, 4, 5, 6],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpConstant, 2),
                    m.make(m.OpConstant, 3),
                    m.make(m.OpConstant, 4),
                    m.make(m.OpConstant, 5),
                    m.make(m.OpHash, 6),
                    m.make(m.OpPop),
                ],
            ),
            (
                "{1: 2 + 3, 4: 5 * 6}",
                [1, 2, 3, 4, 5, 6],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpConstant, 2),
                    m.make(m.OpAdd),
                    m.make(m.OpConstant, 3),
                    m.make(m.OpConstant, 4),
                    m.make(m.OpConstant, 5),
                    m.make(m.OpMul),
                    m.make(m.OpHash, 4),
                    m.make(m.OpPop),
                ],
            ),
        ]
            run_compiler_tests(input, expected_constants, expected_instructions)
        end
    end

    @testset "Index Expressions" begin
        for (input, expected_constants, expected_instructions) in [
            (
                "[1, 2, 3][1 + 1]",
                [1, 2, 3, 1, 1],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpConstant, 2),
                    m.make(m.OpArray, 3),
                    m.make(m.OpConstant, 3),
                    m.make(m.OpConstant, 4),
                    m.make(m.OpAdd),
                    m.make(m.OpIndex),
                    m.make(m.OpPop),
                ],
            ),
            (
                "{1: 2}[2 - 1]",
                [1, 2, 2, 1],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpHash, 2),
                    m.make(m.OpConstant, 2),
                    m.make(m.OpConstant, 3),
                    m.make(m.OpSub),
                    m.make(m.OpIndex),
                    m.make(m.OpPop),
                ],
            ),
        ]
            run_compiler_tests(input, expected_constants, expected_instructions)
        end
    end

    @testset "Functions" begin
        for (input, expected_constants, expected_instructions) in [
            (
                "fn() { return 5 + 10 }",
                [
                    5,
                    10,
                    vcat(
                        m.make(m.OpConstant, 0),
                        m.make(m.OpConstant, 1),
                        m.make(m.OpAdd),
                        m.make(m.OpReturnValue),
                    ),
                ],
                [m.make(m.OpClosure, 2, 0), m.make(m.OpPop)],
            ),
            (
                "fn() { 5 + 10 }",
                [
                    5,
                    10,
                    vcat(
                        m.make(m.OpConstant, 0),
                        m.make(m.OpConstant, 1),
                        m.make(m.OpAdd),
                        m.make(m.OpReturnValue),
                    ),
                ],
                [m.make(m.OpClosure, 2, 0), m.make(m.OpPop)],
            ),
            (
                "fn() { 1; 2 }",
                [
                    1,
                    2,
                    vcat(
                        m.make(m.OpConstant, 0),
                        m.make(m.OpPop),
                        m.make(m.OpConstant, 1),
                        m.make(m.OpReturnValue),
                    ),
                ],
                [m.make(m.OpClosure, 2, 0), m.make(m.OpPop)],
            ),
            (
                "fn() { }",
                [m.make(m.OpReturn)],
                [m.make(m.OpClosure, 0, 0), m.make(m.OpPop)],
            ),
        ]
            run_compiler_tests(input, expected_constants, expected_instructions)
        end
    end

    @testset "Compiler Scopes" begin
        c = m.Compiler()
        g = c.symbol_table
        @test length(c.scopes) == 1

        m.emit!(c, m.OpMul)
        m.enter_scope!(c)
        @test length(c.scopes) == 2

        m.emit!(c, m.OpSub)
        @test length(c) == 1

        last = m.last_instruction(c)
        @test last.op == m.OpSub
        @test c.symbol_table.outer == g

        m.leave_scope!(c)
        @test length(c.scopes) == 1
        @test c.symbol_table == g
        @test isnothing(c.symbol_table.outer)

        m.emit!(c, m.OpAdd)
        @test length(c) == 2

        last = m.last_instruction(c)
        @test last.op == m.OpAdd

        previous = m.prev_instruction(c)
        @test previous.op == m.OpMul
    end

    @testset "Function Calls" begin
        for (input, expected_constants, expected_instructions) in [
            (
                "fn() { 24 }()",
                [24, vcat(m.make(m.OpConstant, 0), m.make(m.OpReturnValue))],
                [m.make(m.OpClosure, 1, 0), m.make(m.OpCall, 0), m.make(m.OpPop)],
            ),
            (
                "let noArg = fn() { 24 }; noArg()",
                [24, vcat(m.make(m.OpConstant, 0), m.make(m.OpReturnValue))],
                [
                    m.make(m.OpClosure, 1, 0),
                    m.make(m.OpSetGlobal, 0),
                    m.make(m.OpGetGlobal, 0),
                    m.make(m.OpCall, 0),
                    m.make(m.OpPop),
                ],
            ),
            (
                "let oneArg = fn(a) { }; oneArg(24)",
                [m.make(m.OpReturn), 24],
                [
                    m.make(m.OpClosure, 0, 0),
                    m.make(m.OpSetGlobal, 0),
                    m.make(m.OpGetGlobal, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpCall, 1),
                    m.make(m.OpPop),
                ],
            ),
            (
                "let manyArg = fn(a, b, c) { }; manyArg(24, 25, 26)",
                [m.make(m.OpReturn), 24, 25, 26],
                [
                    m.make(m.OpClosure, 0, 0),
                    m.make(m.OpSetGlobal, 0),
                    m.make(m.OpGetGlobal, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpConstant, 2),
                    m.make(m.OpConstant, 3),
                    m.make(m.OpCall, 3),
                    m.make(m.OpPop),
                ],
            ),
            (
                "let oneArg = fn(a) { a }; oneArg(24)",
                [vcat(m.make(m.OpGetLocal, 0), m.make(m.OpReturnValue)), 24],
                [
                    m.make(m.OpClosure, 0, 0),
                    m.make(m.OpSetGlobal, 0),
                    m.make(m.OpGetGlobal, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpCall, 1),
                    m.make(m.OpPop),
                ],
            ),
            (
                "let manyArg = fn(a, b, c) { a; b; c; }; manyArg(24, 25, 26)",
                [
                    vcat(
                        m.make(m.OpGetLocal, 0),
                        m.make(m.OpPop),
                        m.make(m.OpGetLocal, 1),
                        m.make(m.OpPop),
                        m.make(m.OpGetLocal, 2),
                        m.make(m.OpReturnValue),
                    ),
                    24,
                    25,
                    26,
                ],
                [
                    m.make(m.OpClosure, 0, 0),
                    m.make(m.OpSetGlobal, 0),
                    m.make(m.OpGetGlobal, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpConstant, 2),
                    m.make(m.OpConstant, 3),
                    m.make(m.OpCall, 3),
                    m.make(m.OpPop),
                ],
            ),
        ]
            run_compiler_tests(input, expected_constants, expected_instructions)
        end
    end

    @testset "Builtin Functions" begin
        for (input, expected_constants, expected_instructions) in [
            (
                "len([]); push([], 1)",
                [1],
                [
                    m.make(m.OpGetBuiltin, 0),
                    m.make(m.OpArray, 0),
                    m.make(m.OpCall, 1),
                    m.make(m.OpPop),
                    m.make(m.OpGetBuiltin, 4),
                    m.make(m.OpArray, 0),
                    m.make(m.OpConstant, 0),
                    m.make(m.OpCall, 2),
                    m.make(m.OpPop),
                ],
            ),
            (
                "fn() { len([]) }",
                [
                    vcat(
                        m.make(m.OpGetBuiltin, 0),
                        m.make(m.OpArray, 0),
                        m.make(m.OpCall, 1),
                        m.make(m.OpReturnValue),
                    ),
                ],
                [m.make(m.OpClosure, 0, 0), m.make(m.OpPop)],
            ),
        ]
            run_compiler_tests(input, expected_constants, expected_instructions)
        end
    end

    @testset "Closures" begin
        for (input, expected_constants, expected_instructions) in [
            (
                """
                fn(a) {
                    fn(b) {
                        a + b
                    }
                }
                """,
                [
                    vcat(
                        m.make(m.OpGetFree, 0),
                        m.make(m.OpGetLocal, 0),
                        m.make(m.OpAdd),
                        m.make(m.OpReturnValue),
                    ),
                    vcat(
                        m.make(m.OpGetLocal, 0),
                        m.make(m.OpClosure, 0, 1),
                        m.make(m.OpReturnValue),
                    ),
                ],
                [m.make(m.OpClosure, 1, 0), m.make(m.OpPop)],
            ),
            (
                """
                fn(a) {
                    fn(b) {
                        fn(c) {
                            a + b + c
                        }
                    }
                }
                """,
                [
                    vcat(
                        m.make(m.OpGetFree, 0),
                        m.make(m.OpGetFree, 1),
                        m.make(m.OpAdd),
                        m.make(m.OpGetLocal, 0),
                        m.make(m.OpAdd),
                        m.make(m.OpReturnValue),
                    ),
                    vcat(
                        m.make(m.OpGetFree, 0),
                        m.make(m.OpGetLocal, 0),
                        m.make(m.OpClosure, 0, 2),
                        m.make(m.OpReturnValue),
                    ),
                    vcat(
                        m.make(m.OpGetLocal, 0),
                        m.make(m.OpClosure, 1, 1),
                        m.make(m.OpReturnValue),
                    ),
                ],
                [m.make(m.OpClosure, 2, 0), m.make(m.OpPop)],
            ),
            (
                """
                let global = 55;
                fn() {
                    let a = 66;
                    fn() {
                        a = a + 2;
                        let b = 77;
                        fn() {
                            global = global + 4;
                            let c = 88;
                            global + a + b + c
                        }
                    }
                }
                """,
                [
                    55,
                    66,
                    2,
                    77,
                    4,
                    88,
                    vcat(
                        m.make(m.OpGetGlobal, 0),
                        m.make(m.OpConstant, 4),
                        m.make(m.OpAdd),
                        m.make(m.OpSetGlobal, 0),
                        m.make(m.OpConstant, 5),
                        m.make(m.OpSetLocal, 0),
                        m.make(m.OpGetGlobal, 0),
                        m.make(m.OpGetFree, 0),
                        m.make(m.OpAdd),
                        m.make(m.OpGetFree, 1),
                        m.make(m.OpAdd),
                        m.make(m.OpGetLocal, 0),
                        m.make(m.OpAdd),
                        m.make(m.OpReturnValue),
                    ),
                    vcat(
                        m.make(m.OpGetFree, 0),
                        m.make(m.OpConstant, 2),
                        m.make(m.OpAdd),
                        m.make(m.OpSetFree, 0),
                        m.make(m.OpConstant, 3),
                        m.make(m.OpSetLocal, 0),
                        m.make(m.OpGetFree, 0),
                        m.make(m.OpGetLocal, 0),
                        m.make(m.OpClosure, 6, 2),
                        m.make(m.OpReturnValue),
                    ),
                    vcat(
                        m.make(m.OpConstant, 1),
                        m.make(m.OpSetLocal, 0),
                        m.make(m.OpGetLocal, 0),
                        m.make(m.OpClosure, 7, 1),
                        m.make(m.OpReturnValue),
                    ),
                ],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpSetGlobal, 0),
                    m.make(m.OpClosure, 8, 0),
                    m.make(m.OpPop),
                ],
            ),
            (
                """
                let i = 2;
                while (i > 0) {
                    i = i - 1;
                    let j = 3;
                    let f = fn() {
                        i + j;
                    }
                    while (j > 0) {
                        j = j - 1;
                        puts(f() + j);
                    }
                }
                """,
                [
                    2,
                    0,
                    1,
                    3,
                    vcat(
                        m.make(m.OpGetGlobal, 0),
                        m.make(m.OpGetFree, 0),
                        m.make(m.OpAdd),
                        m.make(m.OpReturnValue),
                    ),
                    0,
                    1,
                    vcat(
                        m.make(m.OpGetOuter, 1, 1, 0),
                        m.make(m.OpConstant, 6),
                        m.make(m.OpSub),
                        m.make(m.OpSetOuter, 1, 1, 0),
                        m.make(m.OpGetBuiltin, 5),
                        m.make(m.OpGetOuter, 1, 1, 1),
                        m.make(m.OpCall, 0),
                        m.make(m.OpGetOuter, 1, 1, 0),
                        m.make(m.OpAdd),
                        m.make(m.OpCall, 1),
                        m.make(m.OpPop),
                        m.make(m.OpContinue),
                    ),
                    vcat(
                        m.make(m.OpGetGlobal, 0),
                        m.make(m.OpConstant, 2),
                        m.make(m.OpSub),
                        m.make(m.OpSetGlobal, 0),
                        m.make(m.OpConstant, 3),
                        m.make(m.OpSetLocal, 0),
                        m.make(m.OpGetLocal, 0),
                        m.make(m.OpClosure, 4, 1),
                        m.make(m.OpSetLocal, 1),
                        m.make(m.OpGetLocal, 0),
                        m.make(m.OpConstant, 5),
                        m.make(m.OpGreaterThan),
                        m.make(m.OpJumpNotTruthy, 44),
                        m.make(m.OpClosure, 7, 0),
                        m.make(m.OpCall, 0),
                        m.make(m.OpJumpNotTruthy, 44),
                        m.make(m.OpJump, 23),
                        m.make(m.OpNull),
                        m.make(m.OpPop),
                        m.make(m.OpContinue),
                    ),
                ],
                [
                    m.make(m.OpConstant, 0),
                    m.make(m.OpSetGlobal, 0),
                    m.make(m.OpGetGlobal, 0),
                    m.make(m.OpConstant, 1),
                    m.make(m.OpGreaterThan),
                    m.make(m.OpJumpNotTruthy, 28),
                    m.make(m.OpClosure, 8, 0),
                    m.make(m.OpCall, 0),
                    m.make(m.OpJumpNotTruthy, 28),
                    m.make(m.OpJump, 6),
                    m.make(m.OpNull),
                    m.make(m.OpPop),
                ],
            ),
            (
                """
                while (true) {
                    break;
                    continue;
                }
                """,
                [vcat(m.make(m.OpBreak), m.make(m.OpContinue), m.make(m.OpContinue))],
                [
                    m.make(m.OpTrue),
                    m.make(m.OpJumpNotTruthy, 16),
                    m.make(m.OpClosure, 0, 0),
                    m.make(m.OpCall, 0),
                    m.make(m.OpJumpNotTruthy, 16),
                    m.make(m.OpJump, 0),
                    m.make(m.OpNull),
                    m.make(m.OpPop),
                ],
            ),
        ]
            run_compiler_tests(input, expected_constants, expected_instructions)
        end
    end

    @testset "Recursive Functions" begin
        for (input, expected_constants, expected_instructions) in [
            (
                """
                let countDown = fn(x) { countDown(x - 1); };
                countDown(1);
                """,
                [
                    1,
                    vcat(
                        m.make(m.OpCurrentClosure),
                        m.make(m.OpGetLocal, 0),
                        m.make(m.OpConstant, 0),
                        m.make(m.OpSub),
                        m.make(m.OpCall, 1),
                        m.make(m.OpReturnValue),
                    ),
                    1,
                ],
                [
                    m.make(m.OpClosure, 1, 0),
                    m.make(m.OpSetGlobal, 0),
                    m.make(m.OpGetGlobal, 0),
                    m.make(m.OpConstant, 2),
                    m.make(m.OpCall, 1),
                    m.make(m.OpPop),
                ],
            ),
            (
                """
                let wrapper = fn() {
                    let countDown = fn(x) { countDown(x - 1); };
                    countDown(1);
                }
                wrapper();
                """,
                [
                    1,
                    vcat(
                        m.make(m.OpCurrentClosure),
                        m.make(m.OpGetLocal, 0),
                        m.make(m.OpConstant, 0),
                        m.make(m.OpSub),
                        m.make(m.OpCall, 1),
                        m.make(m.OpReturnValue),
                    ),
                    1,
                    vcat(
                        m.make(m.OpClosure, 1, 0),
                        m.make(m.OpSetLocal, 0),
                        m.make(m.OpGetLocal, 0),
                        m.make(m.OpConstant, 2),
                        m.make(m.OpCall, 1),
                        m.make(m.OpReturnValue),
                    ),
                ],
                [
                    m.make(m.OpClosure, 3, 0),
                    m.make(m.OpSetGlobal, 0),
                    m.make(m.OpGetGlobal, 0),
                    m.make(m.OpCall, 0),
                    m.make(m.OpPop),
                ],
            ),
        ]
            run_compiler_tests(input, expected_constants, expected_instructions)
        end
    end

    @testset "Error handling" begin
        @test_throws ErrorException("unknown operator: &") begin
            c = m.Compiler()
            m.compile!(
                c,
                m.InfixExpression(
                    m.Token(m.PLUS, "+"),
                    m.IntegerLiteral(m.Token(m.INT, "2"), 2),
                    "&",
                    m.IntegerLiteral(m.Token(m.INT, "2"), 2),
                ),
            )
        end

        @test_throws ErrorException("unknown operator: +") begin
            c = m.Compiler()
            m.compile!(
                c,
                m.PrefixExpression(
                    m.Token(m.PLUS, "+"),
                    "+",
                    m.IntegerLiteral(m.Token(m.INT, "2"), 2),
                ),
            )
        end
    end
end


text:
export PeriodicTorsion

"""
    PeriodicTorsion(; periodicities, phases, ks, proper)

A periodic torsion angle between four atoms.
"""
struct PeriodicTorsion{N, T, E} <: SpecificInteraction
    periodicities::NTuple{N, Int}
    phases::NTuple{N, T}
    ks::NTuple{N, E}
    proper::Bool
end

function PeriodicTorsion(; periodicities, phases, ks, proper::Bool=true,
                            n_terms=length(periodicities))
    T, E = eltype(phases), eltype(ks)
    if n_terms > length(periodicities)
        n_to_add = n_terms - length(periodicities)
        periodicities_pad = vcat(collect(periodicities), ones(Int, n_to_add))
        phases_pad = vcat(collect(phases), zeros(T, n_to_add))
        ks_pad = vcat(collect(ks), zeros(E, n_to_add))
    else
        periodicities_pad, phases_pad, ks_pad = periodicities, phases, ks
    end
    PeriodicTorsion{n_terms, T, E}(tuple(periodicities_pad...), tuple(phases_pad...),
                                    tuple(ks_pad...), proper)
end

@inline @inbounds function force(d::PeriodicTorsion, coords_i, coords_j, coords_k,
                                    coords_l, box_size)
    ab = vector(coords_i, coords_j, box_size)
    bc = vector(coords_j, coords_k, box_size)
    cd = vector(coords_k, coords_l, box_size)
    cross_ab_bc = ab × bc
    cross_bc_cd = bc × cd
    bc_norm = norm(bc)
    θ = atan(
        ustrip(dot(cross_ab_bc × cross_bc_cd, bc / bc_norm)),
        ustrip(dot(cross_ab_bc, cross_bc_cd)),
    )
    fs = sum(zip(d.periodicities, d.phases, d.ks)) do (periodicity, phase, k)
        dEdθ = -k * periodicity * sin((periodicity * θ) - phase)
        fi =  dEdθ * bc_norm * cross_ab_bc / dot(cross_ab_bc, cross_ab_bc)
        fl = -dEdθ * bc_norm * cross_bc_cd / dot(cross_bc_cd, cross_bc_cd)
        v = (dot(-ab, bc) / bc_norm^2) * fi - (dot(-cd, bc) / bc_norm^2) * fl
        fj =  v - fi
        fk = -v - fl
        return SpecificForce4Atoms(fi, fj, fk, fl)
    end
    return fs
end

@inline @inbounds function potential_energy(d::PeriodicTorsion, coords_i, coords_j, coords_k,
                                            coords_l, box_size)
    θ = torsion_angle(coords_i, coords_j, coords_k, coords_l, box_size)
    E = sum(zip(d.periodicities, d.phases, d.ks)) do (periodicity, phase, k)
        return k + k * cos((periodicity * θ) - phase)
    end
    return E
end


text:
using NestedTuples
using Test

using NestedTuples: with, TypelevelExpr

@testset "NestedTuples.jl" begin
    x = (a = (a = :a, b = :b), q = (l = :l, u = :u))

    @test NestedTuples.exprify(x; rename=false) == ([:a, :b, :l, :u], :((a = (a = a, b = b), q = (l = l, u = u))))

    @test NestedTuples.flatten(x) == (:a, :b, :l, :u)

    @test NestedTuples.keysort(x) == (a = (a = :a, b = :b), q = (l = :l, u = :u))

    @test NestedTuples.leaf_setter(x)(1,2,3,4) == (a = (a = 1, b = 2), q = (l = 3, u = 4))

    @test NestedTuples.schema(x) == (a = (a = Symbol, b = Symbol), q = (l = Symbol, u = Symbol))

    @test NestedTuples.@with((x=1, y=2), x+y) == 3

    @test let nt = (x=1, y=(a=2, b=3))
        @with(nt, (x + @with(y, a + b))) == 6
    end

    @test keysort(((b=1,a=2),(d=3,c=4))) == ((a = 2, b = 1), (c = 4, d = 3))

    @test keysort(lazymerge((b=1,a=2),(d=3,c=4))) == lazymerge((a = 2, b = 1), (c = 4, d = 3))

    @test keysort(NamedTuple()) == NamedTuple()

    @test with(@__MODULE__, (x=1,), (y=2,), TypelevelExpr(:(x+y))) == 3

    @test keysort(3) == 3

    @test convert(NamedTuple, lazymerge((a = 1, b = 2, c = 3), (b = 4, d = 5))) == (a = 1, c = 3, d = 5, b = 4)
end


text:
"""
---
Generates a correlated random rayleigh sequence of size N, parametrized by the doppler frequency fd and the sampling frqeuency fs.
The method uses IFFT approach described in [1]
# --- Syntax
α		= rayleighChan(nbPath,N,fs,fd,seed=-1)
# --- Input parameters
- nbPath : Number of path to be generated (>0) [Int]
- N	: Size of desired output (size of FFT) [Union{Int,Float64}]
- fs	: Sampling frequency [Union{Int,Float64}]
- fd	: doppler frequency [Union{Int,Float64}]
- seed: Seed for random generation ([Int], default =-1)
# --- Output parameters
- α	: Complex rayleigh coefficient [Array{Complex{Float64},N}]
# --- References
[1] D. J. Young and N. C. Beaulieu, "The generation of correlated Rayleigh random variates by inverse discrete Fourier transform," in IEEE Transactions on Communications

# ---
"""
function rayleighChan(nbPath,sizeFFT,fs,fd,randSeed=-1)
	# ----------------------------------------------------
	# --- Parameters init
	# ----------------------------------------------------
	# --- FFT size and FFT parameters
	# Core number
	#FFTW.set_num_threads(4);
	# FFT size
	# Forcing power of 2 for FFT computation
	N	= nextpow(2,sizeFFT);
	# --- Compute normalized doppler frequency
	normFd	= fd / fs;
	# --- Deduce grid-based index associated to normFd
	# This is denoted as km in [1]
	km		= Int(floor(normFd * N));
	# Defined as F_M[k] in [1]
	fM	= zeros(Complex{Float64},N);
	if km == 0
		# --- Target doppler frequency cannot be
		# obtained with frequency grid
		# Residual doppler frequency
		# We set the spectrum to have a constant output
		fM[1] = 1;
	else
		# ----------------------------------------------------
		# --- Filter definition
		# ----------------------------------------------------
		# --- Init filter
		# First area
		for k = 0:km-1
			fM[1+km] = sqrt( 1/ (2*sqrt(1 - (k/(normFd*N)^2))));
		end
		# @ km +1
		fM[km+1]       = sqrt( km/2 *( pi/2 - atan( (km-1)/sqrt(2km-1))));
		# Second area
		fM[km+2:N-km]  .= 0;
		# @ N - km +1
		fM[N - km + 1] = fM[km+1];
		# Third area
		fM[ N-km+2:N]  = fM[km:-1:2];
	end
	# ----------------------------------------------------
	# --- Casting output
	# ----------------------------------------------------
	# --- Seed
	if randSeed != -1
		# ---
		Random.seed!(randSeed);
	end
	# --- Creating random sequence
	# We need to have unitary variance @output of system model
	# /2 for I/Q equivalent power
	# Σ | F[k] | ^2 -> Normalisation wrt doppler filter
	variance  = sqrt(1/2/sum(abs2.(fM)));
	seqPhase  = variance* randn(N,nbPath) .* fM;	  # --- In phase
	seqQuad	  = variance* randn(N,nbPath) .* fM;	  # --- Quadrature
	α 		  = sizeFFT * ifft(seqPhase + 1im* seqQuad,1);
	# --- Creating random sequence
	return α
end


text:
include("operations/save_basic_features_instances.jl")
include("operations/save_features_instances.jl")
include("operations/save_single_features_instances.jl")
include("operations/serialize_instances.jl")
include("operations/generate_decompositions.jl")
include("operations/merge_decompositions.jl")
include("operations/combine_decompositions.jl")
include("operations/check_sanity.jl")
include("operations/delete_duplicates.jl")
include("operations/export_matpowerm.jl")
include("operations/load_matctr_instances.jl")
include("operations/solve_decompositions.jl")

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
include("types.jl")
include("shared.jl")
include("util.jl")
include("adv_cg.jl")
include("adv_bfgs.jl")
include("adv_sgd.jl")

## set

# alg = :bfgs
alg = :sgd

# sgd setting
step = 0.1
use_adagrad = true

verbose = false

### prepare data

dname = "glass"
D_all = readcsv("data-example/" * dname * ".csv")
id_train = readcsv("data-example/" * dname * ".train")
id_test = readcsv("data-example/" * dname * ".test")

id_train = round(Int64, id_train)
id_test = round(Int64, id_test)

println(dname)

### Cross Validation, using first split
## First stage

id_tr = vec(id_train[1,:])
id_ts = vec(id_test[1,:])
X_train = D_all[id_tr,1:end-1]
y_train = round(Int, D_all[id_tr, end])

X_test = D_all[id_ts,1:end-1]
y_test = round(Int, D_all[id_ts, end])

X_train, mean_vector, std_vector = standardize(X_train)
X_test = standardize(X_test, mean_vector, std_vector)

Cs =  [2.0^i for i=0:3:12]
ncs = length(Cs)

# fold
n_train = size(X_train, 1)
n_test = size(X_test, 1)
kf = 5

# k folds
folds = k_fold(n_train, kf)

loss_list = zeros(ncs)
loss01_list = zeros(ncs)

# The first stage of CV
idx = randperm(n_train)
X_train = X_train[idx,:]
y_train = y_train[idx]

for i = 1:ncs

  println(i, " | Adversarial | C = ", string(Cs[i]))

  losses = zeros(n_train)
  losses1 = zeros(n_train)
  # k fold
  for j = 1:kf
    # prepare training and validation
    id_tr = vcat(folds[[1:j-1; j+1:end]]...)
    id_val = folds[j]

    X_tr = X_train[id_tr, :];  y_tr = y_train[id_tr]
    X_val = X_train[id_val, :];  y_val = y_train[id_val]

    print("    ",j, "-th fold : ")
    tol = Cs[i] * length(y_tr) * 1e-6
    if alg == :bfgs
      @time model = train_adv_bfgs(X_tr, y_tr, Cs[i], ftol = tol, grtol = tol, max_iter=1000, show_trace = verbose)
    elseif alg == :sgd
      @time model = train_adv_sgd(X_tr, y_tr, Cs[i], ftol = tol, step = step, use_adagrad = use_adagrad,
                    grtol = tol, max_iter=10000, show_trace = verbose)
    end

    _, ls, _, ls1, _, _ = test_adv(model, X_val, y_val)

    losses[id_val] = ls
    losses1[id_val] = ls1

  end

  loss_list[i] = mean(losses)
  loss01_list[i] = mean(losses1)
  # println("loss : ", string(mean(losses)))
  println("  => loss01 : ", string(mean(losses1)))
  println()

end

ind_max= indmin(loss01_list)
C0 = Cs[ind_max]
Cs =  [C0*2.0^(i-3) for i=1:5]
ncs = length(Cs)

## Second stage
idx = randperm(n_train)
X_train = X_train[idx,:]
y_train = y_train[idx]

for i = 1:ncs

  println(i, " | Adversarial | C = ", string(Cs[i]))

  losses = zeros(n_train)
  losses1 = zeros(n_train)
  # k fold
  for j = 1:kf
    # prepare training and validation
    id_tr = vcat(folds[[1:j-1; j+1:end]]...)
    id_val = folds[j]

    X_tr = X_train[id_tr, :];  y_tr = y_train[id_tr]
    X_val = X_train[id_val, :];  y_val = y_train[id_val]

    print("    ",j, "-th fold : ")
    tol = Cs[i] * length(y_tr) * 1e-6
    if alg == :bfgs
      @time model = train_adv_bfgs(X_tr, y_tr, Cs[i], ftol = tol, grtol = tol, max_iter=1000, show_trace = verbose)
    elseif alg == :sgd
      @time model = train_adv_sgd(X_tr, y_tr, Cs[i], ftol = tol, step = step, use_adagrad = use_adagrad,
                    grtol = tol, max_iter=10000, show_trace = verbose)
    end

    _, ls, _, ls1, _, _ = test_adv(model, X_val, y_val)
    losses[id_val] = ls
    losses1[id_val] = ls1
    #println("loss : ", string(ls))
    #println("loss01 : ", string(ls))

  end

  loss_list[i] = mean(losses)
  loss01_list[i] = mean(losses1)
  # println("loss : ", string(mean(losses)))
  println("  => loss01 : ", string(mean(losses1)))
  println()

end

ind_max= indmin(loss01_list)
C_best = Cs[ind_max]

### Evaluation

n_split = size(id_train, 1)

v_model = Vector{ClassificationModel}()
v_result = Vector{Tuple}()
v_acc = zeros(n_split)

for i = 1:n_split
  # standardize
  id_tr = vec(id_train[i,:])
  id_ts = vec(id_test[i,:])
  X_train = D_all[id_tr,1:end-1]
  y_train = round(Int, D_all[id_tr, end])

  X_test = D_all[id_ts,1:end-1]
  y_test = round(Int, D_all[id_ts, end])

  X_train, mean_vector, std_vector = standardize(X_train)
  X_test = standardize(X_test, mean_vector, std_vector)

  #train and test
  tol = C_best * length(y_train) * 1e-6
  if alg == :bfgs
    @time model = train_adv_bfgs(X_train, y_train, C_best, ftol = tol, grtol = tol, max_iter=1000, show_trace = verbose)
  elseif alg == :sgd
    @time model = train_adv_sgd(X_train, y_train, C_best, ftol = tol, step = step, use_adagrad = use_adagrad,
                  grtol = tol, max_iter=10000, show_trace = verbose)
  end

  result = test_adv(model, X_test, y_test)
  loss01 = result[3]
  acc = 1.0 - loss01

  println("accuracy : ", acc)

  push!(v_model, model)
  push!(v_result, result)
  v_acc[i] = acc
end

println(dname)
println("mean accuracy : ", mean(v_acc))
println("std accuracy : ", std(v_acc))

max_stars_repo_name: Koprvhdix/adversarial-multiclass

max_issues_repo_name: Koprvhdix/adversarial-multiclass

max_forks_repo_name: Koprvhdix/adversarial-multiclass


text:
module SumOfSquares

using LinearAlgebra

import Reexport

import MutableArithmetics
const MA = MutableArithmetics

# MultivariatePolynomials extension

import MultivariatePolynomials
const MP = MultivariatePolynomials
Reexport.@reexport using MultivariateBases
# @set assumes that `SemialgebraicSets` is defined
Reexport.@reexport using SemialgebraicSets
Reexport.@reexport using MultivariateMoments

include("gram_matrix.jl")
include("sosdec.jl")

using PolyJuMP
abstract type SOSLikeCone <: PolyJuMP.PolynomialSet end
Base.broadcastable(cone::SOSLikeCone) = Ref(cone)
# FIXME Maybe replace by PSDLike MOI sets

function matrix_cone_type end

export Certificate
include("Certificate/Certificate.jl")
using .Certificate: Sparsity, SignSymmetry, ChordalCompletion, ClusterCompletion, Symmetry
export Sparsity, SignSymmetry, ChordalCompletion, ClusterCompletion
export Symmetry

# TODO remove in SumOfSquares v1.0
@deprecate NoSparsity Sparsity.NoPattern
@deprecate VariableSparsity Sparsity.Variable
@deprecate MonomialSparsity Sparsity.Monomial

include("rand.jl")

# MOI extension

using MathOptInterface
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

Reexport.@reexport using PolyJuMP

include("attributes.jl")
include("psd2x2.jl")
include("diagonally_dominant.jl")
include("sos_polynomial.jl")
include("copositive_inner.jl")

# Bridges
include("Bridges/Bridges.jl")

# JuMP extension

Reexport.@reexport using JuMP

include("utilities.jl")
include("variable.jl")
include("constraint.jl")

function setdefaults!(data::PolyJuMP.Data)
    PolyJuMP.setdefault!(data, PolyJuMP.NonNegPoly, SOSCone)
    PolyJuMP.setdefault!(data, PolyJuMP.PosDefPolyMatrix, SOSMatrixCone)
    return
end

export SOSModel
function SOSModel(args...; kwargs...)
    model = Model(args...; kwargs...)
    setpolymodule!(model, SumOfSquares)
    model
end

end # module

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
using Test

@testset "TDD.jl" begin
    function neighbor(n)
        a=reshape(graph[n],(1,length(graph[n])))
        b=hcat(n,a)
        println(sort(b,dims=2))
    end
    graph = [
        [2,3],
        [1],
        [1],
        [5],
        [4]]
    @test neighbor(1)==[1 2 3]
    @test neighbor(4)==[4 5]

    graph = [
        [2],
        [1],
        [1],
        [5],
        [4]]
    @test neighbor(1)==[1 2]
    @test neighbor(3)==[1 3]

    graph = [
        [2],
        [1],
        [1],
        [],
        [4]]
    @test neighbor(4)==[4]
end


text:
push!(LOAD_PATH,"../src/")

using Documenter, NURBS

makedocs(
	format = :html,
	sitename = "NURBS.jl",
	assets = ["assets/nurbs.css", "assets/logo.png"],
	pages = [
		"Home" => "index.md",
		"Getting Started" => "gettingStarted.md",
		"Curves" => [
			"Bezier Curves" => "bezierC.md",
			"B-Spline Curves" => "splineC.md",
			"Rational B-Spline Curves" => "rSplineC.md"
		],
		"Surfaces" => [
			"Bezier Surfaces" => "bezierS.md",
			"B-Spline Surfaces" => "splineS.md",
			"Rational B-Spline Surfaces" => "rSplineS.md"
		],
		"Authors" => "authors.md",
		"Project Conclusions" => [
			"About Parallelization" => "parallelization-report.md",
			"Further Developments" => "ToDo.md"
		],
	]
)