text:
"""
```
smc2(m::Union{AbstractDSGEModel,AbstractVARModel}, data::Matrix; verbose::Symbol, old_data::Matrix, kwargs...)
smc(m::Union{AbstractDSGEModel,AbstractVARModel},  data::DataFrame, kwargs...)
smc(m::Union{AbstractDSGEModel,AbstractVARModel},  kwargs...)
```

### Overview

These are wrapper functions to ensure simplicity of estimation of DSGE models while
navigating the DSGE package. Most keywords are passed as settings of the DSGE models
rather than keyword arguments, as in SMC.jl. See `SMC Settings` below.
There are also some additional
keyword arguments that are not passed as settings. Please see `Keyword Arguments` below.

### Arguments

- `m`: A model object, which stores parameter values, prior dists, bounds, and various
    other settings that will be referenced
- `data`: A matrix or dataframe containing the time series of the observables used in
    the calculation of the posterior/likelihood

### SMC Settings
For all of the following settings, unless otherwise noted, the names are the same
as the keyword argument for `SMC.smc`. The names are sometimes different from the
keyword argument for `SMC.smc` to distinguish the setting from MH estimation settings.

- `n_particles::Int`: Number of particles.
- `n_smc_blocks::Int`: Number of parameter blocks in mutation step (corresponds to kwarg `n_blocks` for `SMC.smc`).
- `n_mh_steps_smc::Int`: Number of Metropolis Hastings steps to attempt during the mutation step (corresponds to kwarg `n_mh_steps` for `SMC.smc`).
- `λ::S`: The 'bending coefficient' λ in Φ(n) = (n/N(Φ))^λ
- `n_Φ::Int`: Number of stages in the tempering schedule.
- `resampler_smc::Symbol`: Which resampling method to use (corresponds to kwarg `resampling_method` for `SMC.smc`).
    - `:systematic`: Will use sytematic resampling.
    - `:multinomial`: Will use multinomial resampling.
    - `:polyalgo`: Samples using a polyalgorithm.
- `threshold_ratio::S`: Threshold s.t. particles will be resampled when the population
    drops below threshold * N
- `step_size_smc::S`: Initial scaling factor for covariance of the particles. Controls size of steps in mutation step.
    This value will be adaptively set afterward to reach an accept rate of `target` (see kwarg below).
    This setting corresponds to the kwarg `c` for `SMC.smc`.
- `mixture_proportion::S = 1.0`: The mixture proportion for the mutation step's proposal distribution. See `?mvnormal_mixture_draw` for details.
    Note that a value of 0.9 has commonly been used in applications to DSGE models (see citations below).
    This setting corresponds to the kwarg `α` for `SMC.smc`.
- `target_accept::S`: The initial target acceptance rate for new particles during mutation
    (corresponds to the kwarg `target` for `SMC.smc`).
- `use_fixed_schedule::Bool`: Flag for whether or not to use a fixed tempering (ϕ) schedule.
- `adaptive_tempering_target_smc::S`: Coefficient of the sample size metric to be targeted when solving
    for an endogenous ϕ or 0.0 if using a fixed schedule (corresponds to the kwarg `tempering_target` for `SMC.smc`).
- `previous_data_vintage::String`: the old data vintage from which to start SMC when using a tempered update
    (corresponds to the kwarg `old_data_vintage` for `SMC.smc`).
- `smc_iteration::Int`: The iteration index for the number of times SMC has been run on the
     same data vintage. Primarily for numerical accuracy/testing purposes.

### Keyword Arguments
- `old_data::Matrix{Float64}` = []: A matrix containing the time series of observables of previous data
    (with `data` being the new data) for the purposes of a time tempered estimation
    (that is, using the posterior draws from a previous estimation as the initial set
    of draws for an estimation with new data). Running a bridge estimation
    requires both `old_data` and `old_cloud`.
- `old_cloud::Union{DSGE.ParticleCloud, DSGE.Cloud, SMC.Cloud} = DSGE.Cloud(m, 0)`: old Cloud object
    used to describe particles from a previous estimation with old data. Running a bridge estimation
    requires both `old_data` and `old_cloud`. If running a bridge estimation and
    no `old_cloud` is provided, then it will be loaded using
    the filepaths in `m`. If no `old_cloud` exists, then the bridge estimation will not run.
- `old_model::Union{AbstractDSGEModel, AbstractVARModel} = m`: model object from which we can build
    the old log-likelihood function for a time tempered SMC estimation. It should be possible
    to evaluate the old log-likelihood given `old_data` and the current draw of parameters.
    This may be nontrivial if, for example, *new* parameters have been added to `m` since the old
    estimation. In this case, `old_model` should include the *new* parameters but still return
    the old log-likelihood as the original estimation if given the same *old* parameters.
    By default, we assume the log-likelihood function has not changed
    and therefore coincides with the current one.
- `filestring_addl::Vector{String} = ""`: additional strings to be appended to the save file
    of estimation output.
- `save_intermediate::Bool = false`: Flag for whether one wants to save intermediate Cloud objects
- `intermediate_stage_increment::Int = 10`: Save Clouds at every increment
   (1 = each stage, 10 = every 10th stage, etc.). Useful if you are using a cluster with time
    limits because if you hit the time limit, then you can just
    start from an intermediate stage rather than start over.
- `continue_intermediate::Bool = false`: Flag to indicate whether one is continuing SMC from an
    intermediate stage.
- `intermediate_stage_start::Int = 0`: Intermediate stage at which one wishes to begin the estimation.
- `run_csminwel::Bool = true`: Set to true to run the csminwel algorithm to identify the true posterior mode
    (which may not exist) after completing an estimation. The mode identified by SMC is just
    the particle with the highest posterior value, but we do not check it is actually a mode (i.e.
    the Hessian is negative definite).
- `regime_switching::Bool = false`: Set to true if there are regime-switching parameters. Otherwise, not all the values of the
    regimes will be used or saved.
- `verbose`: Desired frequency of function progress messages printed to standard out.
	- `:none`: No status updates will be reported.
	- `:low`: Status updates for SMC initialization and recursion will be included.
	- `:high`: Status updates for every iteration of SMC is output, which includes
    the mean and standard deviation of each parameter draw after each iteration,
    as well as calculated acceptance rate, ESS, and number of times resampled.
- `log_prob_old_data::Float64 = 0.0`: Log p(\tilde y) which is the log marginal data density
    of the bridge estimation.

### Outputs

- `cloud`: The ParticleCloud object containing all of the information about the
    parameter values from the sample, their respective log-likelihoods, the ESS
    schedule, tempering schedule etc., which is saved in the saveroot.
"""
function smc2(m::Union{AbstractDSGEModel,AbstractVARModel}, data::Matrix{Float64};
              verbose::Symbol = :low,
              old_data::Matrix{Float64} = Matrix{Float64}(undef, size(data, 1), 0),
              old_cloud::Union{DSGE.ParticleCloud, DSGE.Cloud,
                               SMC.Cloud} = DSGE.ParticleCloud(m, 0),
              old_model::Union{AbstractDSGEModel, AbstractVARModel} = m,
              run_test::Bool = false,
              filestring_addl::Vector{String} = Vector{String}(),
              continue_intermediate::Bool = false, intermediate_stage_start::Int = 0,
              save_intermediate::Bool = false, intermediate_stage_increment::Int = 10,
              run_csminwel::Bool = true,
              regime_switching::Bool = false, log_prob_old_data::Float64 = 0.0)

    parallel    = get_setting(m, :use_parallel_workers)
    n_parts     = get_setting(m, :n_particles)
    n_blocks    = get_setting(m, :n_smc_blocks)
    n_mh_steps  = get_setting(m, :n_mh_steps_smc)
    old_vintage = get_setting(m, :previous_data_vintage)

    smc_iteration = get_setting(m, :smc_iteration)

    λ    = get_setting(m, :λ)
    n_Φ  = get_setting(m, :n_Φ)

    # Define tempering settings
    tempered_update_prior_weight = get_setting(m, :tempered_update_prior_weight)
    tempering_target             = get_setting(m, :adaptive_tempering_target_smc)
    use_fixed_schedule           = haskey(get_settings(m), :use_fixed_schedule) ?
        get_setting(m, :use_fixed_schedule) : tempering_target == 0.0

    # Print output for debugging?
    debug_assertion = haskey(get_settings(m), :debug_assertion) && get_setting(m, :debug_assertion)

    # Step 2 (Correction) settings
    resampling_method = get_setting(m, :resampler_smc)
    threshold_ratio   = get_setting(m, :resampling_threshold)
    threshold         = threshold_ratio * n_parts

    # Step 3 (Mutation) settings
    c      = get_setting(m, :step_size_smc)
    α      = get_setting(m, :mixture_proportion)
    target = get_setting(m, :target_accept)

    use_chand_recursion = get_setting(m, :use_chand_recursion)

    my_likelihood = if isa(m, AbstractDSGEModel)
        function _my_likelihood_dsge(parameters::ParameterVector, data::Matrix{Float64})::Float64
            update!(m, parameters)
            likelihood(m, data; sampler = false, catch_errors = true,
                       use_chand_recursion = use_chand_recursion, verbose = verbose)
        end
    else isa(m, AbstractVARModel)
        function _my_likelihood_var(parameters::ParameterVector, data::Matrix{Float64})::Float64
            update!(m, parameters)
            likelihood(m, data; sampler = false, catch_errors = true, verbose = verbose)
        end
    end

    my_old_likelihood = if isa(m, AbstractDSGEModel)
        function _my_old_likelihood_dsge(parameters::ParameterVector, data::Matrix{Float64})::Float64
            update!(old_model, parameters)
            likelihood(old_model, data; sampler = false, catch_errors = true,
                       use_chand_recursion = use_chand_recursion, verbose = verbose)
        end
    else isa(m, AbstractVARModel)
        function _my_old_likelihood_var(parameters::ParameterVector, data::Matrix{Float64})::Float64
            update!(old_model, parameters)
            likelihood(old_model, data; sampler = false, catch_errors = true, verbose = verbose)
        end
    end

    tempered_update = !isempty(old_data)

    # This step is purely for backwards compatibility purposes
    old_cloud_conv = isempty(old_cloud) ? SMC.Cloud(0,0) : SMC.Cloud(old_cloud)

    # Initialize Paths
    loadpath = ""
    if tempered_update
        if isempty(old_cloud)
            loadpath = rawpath(m, "estimate", "smc_cloud.jld2", filestring_addl)
            loadpath = replace(loadpath, r"vint=[0-9]{6}" => "vint=" * old_vintage)
        end
    elseif continue_intermediate
        loadpath = rawpath(m, "estimate", "smc_cloud", filestring_addl) *
            "_stage=$(intermediate_stage_start).jld2"
    end
    savepath = rawpath(m, "estimate", "smc_cloud.jld2", filestring_addl)
    particle_store_path = rawpath(m, "estimate", "smcsave.h5", filestring_addl)

    # Calls SMC package's generic SMC
    println("Calling SMC.jl's SMC estimation routine...")
    SMC.smc(my_likelihood, get_parameters(m), data;
            verbose      = verbose,
            testing      = m.testing,
            data_vintage = data_vintage(m),

            parallel   = parallel,
            n_parts    = n_parts,
            n_blocks   = n_blocks,
            n_mh_steps = n_mh_steps,

            λ = λ, n_Φ = n_Φ,

            resampling_method = resampling_method,
            threshold_ratio   = threshold_ratio,

            c = c, α = α, target = target,

            use_fixed_schedule = use_fixed_schedule,
            tempering_target   = tempering_target,

            old_data          = old_data,
            old_cloud         = old_cloud_conv,
            old_loglikelihood = my_old_likelihood,
            old_vintage       = old_vintage,
            smc_iteration     = smc_iteration,

            run_test = run_test,

            filestring_addl     = filestring_addl,
            loadpath            = loadpath,
            savepath            = savepath,
            particle_store_path = particle_store_path,

            continue_intermediate        = continue_intermediate,
            intermediate_stage_start     = intermediate_stage_start,
            save_intermediate            = save_intermediate,
            intermediate_stage_increment = intermediate_stage_increment,
	        tempered_update_prior_weight = tempered_update_prior_weight,

            regime_switching = regime_switching,
            debug_assertion = debug_assertion, log_prob_old_data = log_prob_old_data)

    if run_csminwel
        m <= Setting(:sampling_method, :SMC)
        update!(m, load_draws(m, :mode, filestring_addl = filestring_addl))
        out, H = optimize!(m, data)
        println("Saving to " * replace(savepath, "smc_cloud" => "paramsmode") * "...")
        jldopen(replace(savepath, "smc_cloud" => "paramsmode"), true, true, true, IOStream) do file
            write(file, "mode", out.minimizer)
        end
    end
end

function smc(m::Union{AbstractDSGEModel,AbstractVARModel}, data::DataFrame; verbose::Symbol = :low,
             old_data::Matrix{Float64} = Matrix{Float64}(undef, size(data, 1), 0),
             old_cloud::Union{DSGE.ParticleCloud, DSGE.Cloud,
                              SMC.Cloud} = DSGE.ParticleCloud(m, 0),
             old_model::Union{AbstractDSGEModel, AbstractVARModel} = m,
             filestring_addl::Vector{String} = Vector{String}(undef, 0),
             run_test::Bool = false,
             save_intermediate::Bool = false, intermediate_stage_increment::Int = 10,
             continue_intermediate::Bool = false, intermediate_stage_start::Int = 0,
             run_csminwel::Bool = true,
             regime_switching::Bool = false, log_prob_old_data::Float64 = 0.0)

    data_mat = df_to_matrix(m, data)
    return smc2(m, data_mat, verbose = verbose,
                old_data = old_data, old_cloud = old_cloud,
                old_model = old_model,
                filestring_addl = filestring_addl, run_test = run_test,
                save_intermediate = save_intermediate,
                intermediate_stage_increment = intermediate_stage_increment,
                continue_intermediate = continue_intermediate,
                intermediate_stage_start = intermediate_stage_start,
                run_csminwel = run_csminwel,
                regime_switching = regime_switching, log_prob_old_data = log_prob_old_data)
end

function smc(m::Union{AbstractDSGEModel,AbstractVARModel}; verbose::Symbol = :low,
             old_data::Matrix{Float64} = Matrix{Float64}(undef, size(data, 1), 0),
             old_cloud::Union{DSGE.ParticleCloud, DSGE.Cloud,
                              SMC.Cloud} = DSGE.ParticleCloud(m, 0),
             old_model::Union{AbstractDSGEModel, AbstractVARModel} = m,
             filestring_addl::Vector{String} = Vector{String}(undef, 0),
             run_test::Bool = false,
             save_intermediate::Bool = false, intermediate_stage_increment::Int = 10,
             continue_intermediate::Bool = false, intermediate_stage_start::Int = 0,
             run_csminwel::Bool = true,
             regime_switching::Bool = false, log_prob_old_data::Float64 = 0.0)

    data = load_data(m)
    data_mat = df_to_matrix(m, data)
    return smc2(m, data_mat, verbose = verbose,
                old_data = old_data, old_cloud = old_cloud,
                old_model = old_model,
                filestring_addl = filestring_addl, run_test = run_test,
                save_intermediate = save_intermediate,
                intermediate_stage_increment = intermediate_stage_increment,
                continue_intermediate = continue_intermediate,
                intermediate_stage_start = intermediate_stage_start,
                run_csminwel = run_csminwel,
                regime_switching = regime_switching, log_prob_old_data = log_prob_old_data)
end

function isempty(c::ParticleCloud)
    length(c.particles) == 0
end

function cloud_isempty(c::ParticleCloud)
    length(c.particles) == 0
end

function get_cloud(m::AbstractDSGEModel; filepath::String = rawpath(m, "estimate", "smc_cloud.jld2"))
    return load(filepath, "cloud")
end

"""
```
function vector_particles_to_cloud(m::AbstractDSGEModel, particles::Vector{Particles})
```
Converts a ParticleCloud (old) to Cloud (new).
"""
function vector_particles_to_cloud(m::AbstractDSGEModel, particles::Vector{Particle})
    cloud = Cloud(m, length(particles))
    N = size(cloud.particles, 2)
    for i in 1:size(cloud.particles,1)
        cloud.particles[i, 1:length(m.parameters)] = particles[i].value
        cloud.particles[i, ind_loglh(N)] = particles[i].loglh
        update!(m, particles[i].value)
        cloud.particles[i, ind_logprior(N)] = prior(m)
        cloud.particles[i, ind_old_loglh(N)] = particles[i].old_loglh
        cloud.particles[i, ind_accept(N)] = particles[i].accept
        cloud.particles[i, ind_weight(N)] = particles[i].weight
    end
    return cloud
end


text:
# Mentoring note: CircularDeque below is sufficient to complete the exercise without bonus
# tasks. The lengthier CircularBuffer afterward completes the bonus tasks and is optimized
# for performance.

module CircularDeques

"""
    CircularDeque{T}(n)

Empty circular deque with capacity for `n` elements of type `T`.
"""
mutable struct CircularDeque{T}
    length::Int  # number of inserted items
    # Can use firstindex instead. When fully optimized, tracking lastindex leads to a faster
    # push! whereas firstindex leads to a faster pushfirst!. Alternatively, instead of
    # tracking length and one end, can track both ends for symmetric performance.
    lastindex::Int  # index in data of the last (i.e., most-recently push!ed) item
    data::Vector{T}

    CircularDeque{T}(n) where {T} = new(0, n, Vector{T}(undef, n))
end

Base.length(cd::CircularDeque) = cd.length
Base.isempty(cd::CircularDeque) = length(cd) == 0
capacity(cd::CircularDeque) = length(cd.data)
isfull(cd::CircularDeque) = length(cd) == capacity(cd)

Base.empty!(cd::CircularDeque) = (cd.length = 0; cd)

function Base.push!(cd::CircularDeque, item; overwrite::Bool=false)
    overwrite || !isfull(cd) || throw(BoundsError(cd, length(cd) + 1))
    cd.lastindex = mod1(cd.lastindex + 1, capacity(cd))
    cd.data[cd.lastindex] = item
    cd.length = min(capacity(cd), cd.length + 1)
    return cd
end

function Base.popfirst!(cd::CircularDeque)
    !isempty(cd) || throw(BoundsError(cd, 0))
    i = mod1(cd.lastindex - length(cd) + 1, capacity(cd))
    cd.length -= 1
    return cd.data[i]
end

end  # module CircularDeques

using .CircularDeques

module CircularBuffers

using Base: @propagate_inbounds
using Printf

export CircularBuffer, capacity, isfull

"""
    CircularBuffer{T}(n) <: AbstractVector{T}

Empty circular buffer with capacity for `n` elements of type `T`.

`CicularBuffer`s are optimized for inserting elements at the end¸ so `push!` performs better
(~1.5x faster) than `pushfirst!`. Deleting elements at either end (`pop!`, and `popfirst!`)
as well as `getindex` and `setindex!` on random elements are all efficient.
"""
mutable struct CircularBuffer{T} <: AbstractVector{T}
    length::Int  # number of inserted items
    lastindex::Int  # index in data of the last (i.e., most-recently push!ed) item
    data::Vector{T}

    CircularBuffer{T}(n) where {T} = new(0, n, Vector{T}(undef, n))
end

# Plain arrays are IndexLinear, but AbstractArrays are IndexCartesian by default for
# compatibility reasons that are irrelevant here. IndexLinear provides speedups in some
# specialized code in Base.
Base.IndexStyle(::Type{<:CircularBuffer}) = IndexLinear()

Base.size(cb::CircularBuffer) = (cb.length,)
Base.empty!(cb::CircularBuffer) = (cb.length = 0; cb)

"""
    capacity(cb::CircularBuffer) -> Int

Return the maximum number of elements that `cb` can store.
"""
capacity(cb::CircularBuffer) = length(cb.data)

"""
    isfull(cb::CircularBuffer) -> Bool

Determine whether `cb` is full.
"""
isfull(cb::CircularBuffer) = length(cb) == capacity(cb)

"""
    _dataindex(cb::CircularBuffer, i::Integer) -> Int

Given an index `i` of `cb`, return the corresponding index `j` of `cb.data` such that `cb[i]
== cb.data[j]`.
"""
@inline function _dataindex(cb::CircularBuffer, i::Integer)
    @boundscheck checkbounds(cb, i)
    j = cb.lastindex - length(cb) + i
    return j < 1 ? j + length(cb.data) : j  # optimized vs mod1
end

@propagate_inbounds function Base.getindex(cb::CircularBuffer, i::Integer)
    j = _dataindex(cb, i)
    @inbounds return cb.data[j]
end

@propagate_inbounds function Base.setindex!(cb::CircularBuffer, item, i::Integer)
    j = _dataindex(cb, i)
    @inbounds return cb.data[j] = item
end

"""
    push!(cb::CircularBuffer, items...; overwrite=false) -> cb

Insert each item from `items` to the end of `cb`, then return `cb`. `overwrite` determines
whether items at the beginning of `cb` can be removed to provide space for new items. If
`overwrite=false` and `push!` attempts to insert an item into a full `CircularBuffer`, then
throw a `BoundsError`.
"""
@inline function Base.push!(cb::CircularBuffer, item; overwrite::Bool=false)
    full = isfull(cb)
    @boundscheck overwrite || !full || throw(BoundsError(cb, length(cb) + 1))
    i = cb.lastindex
    cb.lastindex = (i == capacity(cb) ? 1 : i + 1)
    @inbounds cb.data[cb.lastindex] = item
    !full && (cb.length += 1)
    return cb
end

@propagate_inbounds Base.push!(cb::CircularBuffer, itr...; kw...) = append!(cb, itr; kw...)

@propagate_inbounds function Base.append!(cb::CircularBuffer, itr; kw...)
    for item in itr
        push!(cb, item; kw...)
    end
    return cb
end

"""
    pushfirst!(cb::CircularBuffer, item; overwrite=false) -> cb

Like `push!`, but insert `item` to the *beginning* of `cb`. See [`push!`](@ref).

Note that `pushfirst!(cb::CircularBuffer, items...)` and
`prepend!(cb::CircularBuffer, items)` are overridden to throw `MethodError` due to ambiguity
about the order of insertion.
"""
@inline function Base.pushfirst!(cb::CircularBuffer, item; overwrite::Bool=false)
    full = isfull(cb)
    @boundscheck overwrite || !full || throw(BoundsError(cb, 0))
    if full
        i = cb.lastindex
        cb.lastindex = (i == 1 ? capacity(cb) : i - 1)
    else
        cb.length += 1
    end
    @inbounds cb[1] = item
    return cb
end

Base.pushfirst!(cb::CircularBuffer, args...) = throw(MethodError(pushfirst!, (cb, args...)))
Base.prepend!(cb::CircularBuffer, args...) = throw(MethodError(pushfirst!, (cb, args...)))

@propagate_inbounds function Base.popfirst!(cb::CircularBuffer)
    out = first(cb)
    cb.length -= 1
    return out
end

@propagate_inbounds function Base.pop!(cb::CircularBuffer)
    out = last(cb)
    cb.length -= 1
    i = cb.lastindex
    cb.lastindex = (i == 1 ? capacity(cb) : i - 1)
    return out
end

end  # module CircularBuffers

using .CircularBuffers

enable_bonus_tests = true

ext: jl


text:
ENV["GKSwstype"] = "100"
ENV["PLOTS_TEST"] = "true"
using Documenter, Wavelets, WaveletsExt

makedocs(
    sitename = "WaveletsExt.jl",
    format = Documenter.HTML(),
    authors = "Zeng Fung Liew, Shozen Dan",
    clean = true,
    pages = Any[
        "Home" => "index.md",
        "Manual" => Any[
            "Transforms" => "manual/transforms.md",
            "Best Basis" => "manual/bestbasis.md",
            "Denoising" => "manual/denoising.md",
            "Local Discriminant Basis" => "manual/localdiscriminantbasis.md"
        ],
        "API" => Any[
            "DWT" => "api/dwt.md",
            "ACWT" => "api/acwt.md",
            "SWT" => "api/swt.md",
            "TIWT" => "api/tiwt.md",
            "Best Basis" => "api/bestbasis.md",
            "Denoising" => "api/denoising.md",
            "LDB" => "api/ldb.md",
            "SIWPD" => "api/siwpd.md",
            "Utils" => "api/utils.md",
            "Visualizations" => "api/visualizations.md",
        ]
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/UCD4IDS/WaveletsExt.jl.git"
)


text:
using PrettyLogging
using Base.Test

L = Logger(STDOUT, 0, "", 1)

logspecial(L, "starting tests", 1)
logmsg(L, "butts")
loghold(L, "processing [")
for i = 1:10
    loghold(L,".")
    L2 = indent(L)
    logmsg(L2, "testing interwoven indentations $i")
end
loghold(L, "]")

logmsg(L, "Done!")


text:
export FuseLowlevelOps
export FUSE_INO_ROOT, RENAME_NOREPLACE, RENAME_EXCHANGE, XATTR_CREATE, XATTR_REPLACE
export FuseConnInfo, FuseFileInfo, FuseEntryParam, FuseCmdlineArgs, FuseCmdlineOpts, FuseReq, FuseIno, FuseMode
export FuseBufvec, FuseBuf, FuseBufFlags, FuseBufCopyFlags
export Cstat, Cflock, Timespec

export FUSE_SET_ATTR_MODE, FUSE_SET_ATTR_UID, FUSE_SET_ATTR_GID, FUSE_SET_ATTR_SIZE
export FUSE_SET_ATTR_ATIME, FUSE_SET_ATTR_MTIME, FUSE_SET_ATTR_ATIME_NOW, FUSE_SET_ATTR_MTIME_NOW, FUSE_SET_ATTR_CTIME

const CFu = Ptr{Cvoid}

"""
    FuseLoopConfig <: Layout

Layout template as input for `fuse_session_loop_mt`.

Named `struct fuse_loop_config` in `libfuse3`.
"""
mutable struct FuseLoopConfig
    clone_fd::Cint
    max_idle_threads::Cuint
end

"""
    FuseCmdlineArgs <: Layout

Layout template for commandline arguments as to be passed to fuse_parse_cmdline

The length of vector pointed to by `argv` is stored in field `argc`.
Field `allocated` must be set to `0`, if the data is Julia generated.

If generated by Julia, use `CStructGuarded{FuseCmdlineArgs}`.
"""
struct FuseCmdlineArgs <: Layout
    argc::Cint
    argv::Ptr{LVarVector{Cstring, (x) -> x.argc + 1}}
    allocated::Cint
end

"""
    FuseCmdlineOpts <: Layout

Layout template for output argument of `fuse_parse_cmdline`.

Has to be allocated by Julia caller, by `CStructGuarded{FuseCmdlineOpts}`.

Field `mountpoint` has to be copied immediately ofter the call, because
the `*char` pointer stored in there is unsafe.
"""
struct LFuseCmdlineOpts <: Layout
    singlethread::Cint
    foreground::Cint
    debug::Cint
    nodefault_subtype::Cint
    mountpoint::Cstring
    show_version::Cint
    show_help::Cint
    clone_fd::Cint
    max_idle_threads::Cuint
end
struct FuseCmdlineOpts
    singlethread::Bool
    foreground::Bool
    debug::Bool
    show_version::Bool
    show_help::Bool
    clone_fd::Bool
    nodefault_subtype::Int
    mountpoint::String
    max_idle_threads::Int
end

"""
    FuseSession <: Layout

Opaque object as returned by `fuse_session_new` and used by subsequent
fuse session management calls. (see `main_loop`)
"""
struct FuseSession <: Layout
end

"""
    FuseLowlevelOps

This structure mainly defines the names and sequence of callback functions,
which can be used by the `libfuse3` lowlevel interface.

Named `struct fuse_lowlevel_ops` in `Fuse3` library.
"""
struct FuseLowlevelOps
    init::CFu
    destroy::CFu
    lookup::CFu
    forget::CFu
    getattr::CFu
    setattr::CFu
    readlink::CFu
    mknod::CFu
    mkdir::CFu
    unlink::CFu
    rmdir::CFu
    symlink::CFu
    rename::CFu
    link::CFu
    open::CFu
    read::CFu
    write::CFu
    flush::CFu
    release::CFu
    fsync::CFu
    opendir::CFu
    readdir::CFu
    releasedir::CFu
    fsyncdir::CFu
    statfs::CFu
    setxattr::CFu
    getxattr::CFu
    listxattr::CFu
    removexattr::CFu
    access::CFu
    create::CFu
    getlk::CFu
    setlk::CFu
    bmap::CFu
    ioctl::CFu
    poll::CFu
    write_buf::CFu
    retrieve_reply::CFu
    forget_multi::CFu
    flock::CFu
    fallocate::CFu
    readdirplus::CFu
    copy_file_range::CFu
    lseek::CFu
end
const F_SIZE = sizeof(FuseLowlevelOps) ÷ sizeof(CFu
)

const FuseIno = UInt64
const FuseMode = UInt32
const FuseDev = UInt64

const Cuid_t = UInt32
const Cgid_t = UInt32
const Coff_t = Csize_t
const Coff64_t = UInt64
const Cuint64_t = UInt64
const Cpid_t = Cint
const Cfsblkcnt_t = UInt64
const Cfsfilcnt_t = UInt64

const FUSE_INO_ROOT = FuseIno(1)

struct FuseBufFlags
    flag::Cint
end
const FUSE_BUF_IS_FD = FuseBufFlags(1 << 1)
const FUSE_BUF_FD_SEEK = FuseBufFlags(1 << 2)
const FUSE_BUF_FD_RETRY = FuseBufFlags(1 << 3)

struct FuseBufCopyFlags
    flag::Cint
end
const FUSE_BUF_NO_SPLICE = FuseBufCopyFlags(1 << 1)
const FUSE_BUF_FORCE_SPLICE = FuseBufCopyFlags(1 << 2)
const FUSE_BUF_SPLICE_MOVE = FuseBufCopyFlags(1 << 3)
const FUSE_BUF_SPLICE_NONBLOCK = FuseBufCopyFlags(1 << 4)

"""
    FuseReq

Opaque structure to identify the request. Contains the pointer as obtained by fuselib.

All typical fuse callback functions use that as their first argument.
"""
struct FuseReq
    pointer::Ptr{Nothing}
end

"""
    Timespec <: Layout

C-structure as contained in `Cstat` (`struct stat`).
"""
struct Timespec <: Layout
    seconds::Int64
    nanoseconds::Int64
end

"""
    FuseCtx <: Layout

C-structure returned by `fuse_req_ctx(req)` containg `uid`, `gid`, `pid`, and `umask`
of the process, which is responsible for the request.

The data is only valid while the request `req` is being processed within one of the
callback bunctions.
"""
struct FuseCtx <: Layout
    uid::Cuid_t
    gid::Cgid_t
    pid::Cpid_t
    umask::FuseMode
end

"""
    Cstat <: Layout

C-structure containing inode data. Named `struct stat` in std C library.
"""
struct Cstat <: Layout
    dev     :: UInt64
    ino     :: UInt64
    nlink   :: UInt64
    mode    :: FuseMode
    uid     :: Cuid_t
    gid     :: Cgid_t
    pad0    :: UInt32
    rdev    :: UInt64
    size    :: Int64
    blksize :: Int64
    blocks  :: Int64
    atime   :: Timespec
    mtime   :: Timespec
    ctime   :: Timespec
end

"""
    FuseEntryParam <: Layout

Fuse wrapper type for a `struct stat`, adding some extra fields.
Contains an embedded field `attr::Cstat` with `struct stat` data.
"""
struct FuseEntryParam <: Layout
    ino::FuseIno
    generation::UInt64
    attr::Cstat
    attr_timeout::Cdouble
    entry_timeout::Cdouble
end

# Capability bits for 'fuse_conn_info.capable' and 'fuse_conn_info.want'
const FUSE_CAP_ASYNC_READ = Cuint(1 << 0)
const FUSE_CAP_POSIX_LOCKS = Cuint(1 << 1)
const FUSE_CAP_ATOMIC_O_TRUNC = Cuint(1 << 3)
const FUSE_CAP_EXPORT_SUPPORT = Cuint(1 << 4)
const FUSE_CAP_DONT_MASK = Cuint(1 << 6)
const FUSE_CAP_SPLICE_WRITE = Cuint(1 << 7)
const FUSE_CAP_SPLICE_MOVE = Cuint(1 << 8)
const FUSE_CAP_SPLICE_READ = Cuint(1 << 9)
const FUSE_CAP_FLOCK_LOCKS = Cuint(1 << 10)
const FUSE_CAP_IOCTL_DIR = Cuint(1 << 11)
const FUSE_CAP_AUTO_INVAL_DATA = Cuint(1 << 12)
const FUSE_CAP_READDIRPLUS = Cuint(1 << 13)
const FUSE_CAP_READDIRPLUS_AUTO = Cuint(1 << 14)
const FUSE_CAP_ASYNC_DIO = Cuint(1 << 15)
const FUSE_CAP_WRITEBACK_CACHE = Cuint(1 << 16)
const FUSE_CAP_NO_OPEN_SUPPORT = Cuint(1 << 17)
const FUSE_CAP_PARALLEL_DIROPS = Cuint(1 << 18)
const FUSE_CAP_POSIX_ACL = Cuint(1 << 19)
const FUSE_CAP_HANDLE_KILLPRIV = Cuint(1 << 20)
const FUSE_CAP_CACHE_SYMLINKS = Cuint(1 << 23)
const FUSE_CAP_NO_OPENDIR_SUPPORT = Cuint(1 << 24)
const FUSE_CAP_EXPLICIT_INVAL_DATA = Cuint(1 << 25)

"""
    FuseConnInfo <: Layout

"""
struct FuseConnInfo <: Layout
    proto_major::Cuint
    proto_minor::Cuint
    max_write::Cuint
    max_read::Cuint
    max_readahead::Cuint
    capable::Cuint
    want::Cuint
    max_background::Cuint
    congestion_threshold::Cuint
    time_gran::Cuint
    # reserved::LFixedVector{Cuint,22}
end

# bit masks for 2nd field of FuseFileInfo
const FUSE_FI_WRITEPAGE = Cuint(1 << 0)
const FUSE_FI_DIRECT_IO = Cuint(1 << 1)
const FUSE_FI_KEEP_CACHE = Cuint(1 << 2)
const FUSE_FI_FLUSH = Cuint(1 << 3)
const FUSE_FI_NONSEEKABLE = Cuint(1 << 4)
const FUSE_FI_CACHE_READDIR = Cuint(1 << 5)

"""
    FuseFileInfo <: Layout

Layout of C-structure used by many fuse callbacks.

The C-structure is allocated and populated by the `libfuse3` C-library
and passed as a pointer to the callback functions as an argument
`fi::Ptr{FuseFileInfo}`, which is made accessible to Julia by
the `CStruct(fi)` conversion.
"""
struct FuseFileInfo <: Layout
    flags::Cint
    bits::Cuint
    fh::UInt64
    lock_owner::UInt64
    poll_events::UInt32
end

"""
    Cflock <: Layout

C-structure for `flock` callback. Named `struct flock` in std C library.
"""
struct Cflock <: Layout
    type::Cshort
    whence::Cshort
    start::Coff64_t
    len::Coff64_t
    pid::Cpid_t
end

"""
    Ciovec <: Layout

C-structure for `ioctl` callback. Named `struct iovec` in std C library.
"""
struct Ciovec <: Layout
    base::Ptr{Cvoid}
    len::Csize_t
end

"""
    Cstatfs <: Layout

C-structure for `statfs` callback. Named `struct statfs` in std C library.
"""
struct Cstatvfs <: Layout
    bsize::Culong
    frsize::Culong
    blocks::Cfsblkcnt_t
    bfree::Cfsblkcnt_t
    bavail::Cfsblkcnt_t
    files::Cfsfilcnt_t
    ffree::Cfsfilcnt_t
    favail::Cfsfilcnt_t
    fsid::Clong
    flag::Culong
    namemax::Culong
    __spare::NTuple{6,Cint}
end

struct FusePollHandle <: Layout
end

"""
    FuseBuf <: Layout

Flexible size buffer containing a pointer to variable size data and a size.

The buffer object itself has fixed size.
Used in conjunction with `FuseBufvec`.
"""
struct FuseBuf <: Layout
    size::Csize_t
    flags::FuseBufFlags
    mem::Ptr{LVarVector{UInt8, (x) -> x.size}}
    fd::Cint
    pos::Coff_t
end

"""
    FuseBufvec <: Layout

Variable sized vector of buffers of type `FuseBuf`.
"""
struct FuseBufvec <: Layout
    count::Csize_t
    idx::Csize_t
    off::Csize_t
    buf::LVarVector{FuseBuf, (x) -> x.count}
end
struct FuseForgetData
end

# used for setattr function
const FUSE_SET_ATTR_MODE      = Cuint(1 << 0)
const FUSE_SET_ATTR_UID       = Cuint(1 << 1)
const FUSE_SET_ATTR_GID       = Cuint(1 << 2)
const FUSE_SET_ATTR_SIZE      = Cuint(1 << 3)
const FUSE_SET_ATTR_ATIME     = Cuint(1 << 4)
const FUSE_SET_ATTR_MTIME     = Cuint(1 << 5)
const FUSE_SET_ATTR_ATIME_NOW = Cuint(1 << 7)
const FUSE_SET_ATTR_MTIME_NOW = Cuint(1 << 8)
const FUSE_SET_ATTR_CTIME     = Cuint(1 << 10)

# used for ioctl function
const FUSE_IOCTL_COMPAT = Cuint(1 << 0)
const FUSE_IOCTL_UNRESTRICTED = Cuint(1 << 1)
const FUSE_IOCTL_RETRY = Cuint(1 << 2)
const FUSE_IOCTL_DIR = Cuint(1 << 4)
const FUSE_IOCTL_MAX_IOV = 256

# for replace function
const RENAME_NOREPLACE = (1 << 0)
const RENAME_EXCHANGE = (1 << 1)

# for setxattr function
const XATTR_CREATE = 1
const XATTR_REPLACE = 2


text:
__precompile__(true)
module CLArrays

using GPUArrays
using GPUArrays: LocalMemory
using OpenCL
using Transpiler
import Transpiler: cli

function context end

include("memory.jl")
include("array.jl")
include("ondevice.jl")
include("device.jl")
include("context.jl")
include("compilation.jl")
include("3rdparty.jl")

export CLArray, gpu_call, is_gpu, is_cpu

end # module


text:
export AbstractStage,
       PreEpisodeStage,
       PRE_EPISODE_STAGE,
       PostEpisodeStage,
       POST_EPISODE_STAGE,
       PreActStage,
       PRE_ACT_STAGE,
       PostActStage,
       POST_ACT_STAGE,
       AbstractHook,
       ComposedHook,
       EmptyHook,
       StepsPerEpisode,
       RewardsPerEpisode,
       TotalRewardPerEpisode,
       CumulativeReward,
       TimePerStep

abstract type AbstractStage end

struct PreEpisodeStage <: AbstractStage end
struct PostEpisodeStage <: AbstractStage end
struct PreActStage <: AbstractStage end
struct PostActStage <: AbstractStage end

const PRE_EPISODE_STAGE = PreEpisodeStage()
const POST_EPISODE_STAGE = PostEpisodeStage()
const PRE_ACT_STAGE = PreActStage()
const POST_ACT_STAGE = PostActStage()

"""
A hook is called at different stage duiring a [`run`](@ref). One can inject customized runtime logic in it.
"""
abstract type AbstractHook end

(hook::AbstractHook)(::T, agent, env, obs) where {T<:AbstractStage} = nothing

# https://github.com/JuliaLang/julia/issues/14919
# function (f::Function)(stage::T, args...;kw...) where T<: AbstractStage end

#####
# ComposedHook
#####

"""
    ComposedHook(hooks...)

Compose different hooks into a single hook.
"""
struct ComposedHook{T<:Tuple} <: AbstractHook
    hooks::T
    ComposedHook(hooks...) = new{typeof(hooks)}(hooks)
end

function (hook::ComposedHook)(stage::AbstractStage, args...; kw...)
    for h in hook.hooks
        h(stage, args...; kw...)
    end
end

Base.getindex(hook::ComposedHook, inds...) = getindex(hook.hooks, inds...)

#####
# EmptyHook
#####

"""
Do nothing
"""
struct EmptyHook <: AbstractHook end

const EMPTY_HOOK = EmptyHook()

#####
# StepsPerEpisode
#####

"""
    StepsPerEpisode(; steps = Int[], count = 0, tag = "TRAINING")

Store steps of each episode in the field of `steps`.
"""
mutable struct StepsPerEpisode <: AbstractHook
    steps::Vector{Int}
    count::Int
    tag::String
end

function StepsPerEpisode(; steps = Int[], count = 0, tag = "TRAINING")
    StepsPerEpisode(steps, count, tag)
end

function (hook::StepsPerEpisode)(::PostActStage, agent, env, action_obs)
    hook.count += 1
end

function (hook::StepsPerEpisode)(::PostEpisodeStage, agent, env, obs)
    push!(hook.steps, hook.count)
    hook.count = 0
    @debug hook.tag STEPS_PER_EPISODE = hook.steps[end]
end

#####
# RewardsPerEpisode
#####

"""
    RewardsPerEpisode(; rewards = Vector{Vector{Float64}}(), tag = "TRAINING")

Store each reward of each step in every episode in the field of `rewards`.
"""
mutable struct RewardsPerEpisode <: AbstractHook
    rewards::Vector{Vector{Float64}}
    tag::String
end

function RewardsPerEpisode(; rewards = Vector{Vector{Float64}}(), tag = "TRAINING")
    RewardsPerEpisode(rewards, tag)
end

function (hook::RewardsPerEpisode)(::PreEpisodeStage, agent, env, obs)
    push!(hook.rewards, [])
end

function (hook::RewardsPerEpisode)(::PostActStage, agent, env, action_obs)
    action, obs = action_obs
    push!(hook.rewards[end], get_reward(obs))
end

function (hook::RewardsPerEpisode)(::PostEpisodeStage, agent, env, obs)
    @debug hook.tag REWARDS_PER_EPISODE = hook.rewards[end]
end

#####
# TotalRewardPerEpisode
#####

"""
    TotalRewardPerEpisode(; rewards = Float64[], reward = 0.0, tag = "TRAINING")

Store the total rewards of each episode in the field of `rewards`.
"""
mutable struct TotalRewardPerEpisode <: AbstractHook
    rewards::Vector{Float64}
    reward::Float64
    tag::String
end

function TotalRewardPerEpisode(; rewards = Float64[], reward = 0.0, tag = "TRAINING")
    TotalRewardPerEpisode(rewards, reward, tag)
end

function (hook::TotalRewardPerEpisode)(::PostActStage, agent, env, action_obs)
    action, obs = action_obs
    hook.reward += get_reward(obs)
end

function (hook::TotalRewardPerEpisode)(::PostEpisodeStage, agent, env, obs)
    push!(hook.rewards, hook.reward)
    hook.reward = 0
    @debug hook.tag REWARD_PER_EPISODE = hook.rewards[end]
end

#####
# CumulativeReward
#####

"""
    CumulativeReward(rewards::Vector{Float64} = [0.0], tag::String = "TRAINING")

Store cumulative rewards since the beginning to the field of `rewards`.
"""
Base.@kwdef struct CumulativeReward <: AbstractHook
    rewards::Vector{Float64} = [0.0]
    tag::String = "TRAINING"
end

function (hook::CumulativeReward)(::PostActStage, agent, env, action_obs)
    action, obs = action_obs
    push!(hook.rewards, get_reward(obs) + hook.rewards[end])
    @debug hook.tag CUMULATIVE_REWARD = hook.rewards[end]
end

#####
# TimePerStep
#####

"""
    TimePerStep(;max_steps=100)
    TimePerStep(times::CircularArrayBuffer{Float64, 1}, t::UInt64)

Store time cost of the latest `max_steps` in the `times` field.
"""
mutable struct TimePerStep <: AbstractHook
    times::CircularArrayBuffer{Float64, 1}
    t::UInt64
end

TimePerStep(;max_steps=100) = TimePerStep(CircularArrayBuffer{Float64}(max_steps), time_ns())

function (hook::TimePerStep)(::PostActStage, agent, env, obs_action)
    push!(hook.times, (time_ns() - hook.t)/1e9)
    hook.t = time_ns()
end


text:
@testset "predicates" verbose = true begin
    inc = [1.0, 2.0, 3.0, 3.0, 4.0]
    sinc = [1.0, 2.0, 3.0, 4.0, 5.0]
    ninc = [1.0, 2.0, 3.0, 2.0, 4.0]
    @test is_increasing(inc)
    @test !is_increasing(inc; strict=true)
    @test is_increasing(sinc; strict=true)
    @test !is_increasing(ninc)
end


text:
module DiffEqsWrapperTests

using Test
using Gridap
using Gridap.ODEs
using Gridap.ODEs.ODETools
using Gridap.ODEs.TransientFETools
using Gridap.ODEs.DiffEqWrappers

# using DifferentialEquations
# using Sundials
using Gridap.Algebra: NewtonRaphsonSolver
using Base.Iterators

# FE problem (heat eq) using Gridap
function fe_problem(u, n)

  f(t) = x -> ∂t(u)(x, t) - Δ(u(t))(x)

  domain = (0, 1, 0, 1)
  partition = (n, n)
  model = CartesianDiscreteModel(domain, partition)

  order = 1

  reffe = ReferenceFE(lagrangian,Float64,order)
  V0 = FESpace(
    model,
    reffe,
    conformity = :H1,
    dirichlet_tags = "boundary",
  )
  U = TransientTrialFESpace(V0, u)

  Ω = Triangulation(model)
  degree = 2 * order
  dΩ = Measure(Ω, degree)

  a(u, v) = ∫( ∇(v) ⋅ ∇(u) )dΩ
  b(v, t) = ∫( v * f(t) )dΩ
  m(u, v) = ∫( v * u )dΩ

  res(t, u, v) = a(u, v) + m(∂t(u), v) - b(v, t)
  jac(t, u, du, v) = a(du, v)
  jac_t(t, u, dut, v) = m(dut, v)

  op = TransientFEOperator(res, jac, jac_t, U, V0)

  U0 = U(0.0)
  uh0 = interpolate_everywhere(u(0.0), U0)
  u0 = get_free_dof_values(uh0)

  return op, u0

end

# Solving the heat equation using Gridap.ODEs and DiffEqs
tspan = (0.0, 1.0)

u(x, t) = t
u(t) = x -> u(x, t)

# ISSUE 1: When I choose n > 2, even though the problem that we will solve is
# linear, the Sundials solvers seems to have convergence issues in the nonlinear
# solver (?). Ut returns errors
# [IDAS ERROR]  IDACalcIC Newton/Linesearch algorithm failed to converge.
# ISSUE 2: When I pass `jac_prototype` the code gets stuck

n = 3 # cells per dim (2D)
op, u0 = fe_problem(u,n)

# Some checks
res!, jac!, mass!, stif! = diffeq_wrappers(op)
J = prototype_jacobian(op,u0)
r = copy(u0)
θ = 1.0; t0 = 0.0; tF = 1.0; dt = 0.1; tθ = 1.0; dtθ = dt*θ
res!(r, u0, u0, nothing, tθ)
jac!(J, u0, u0, nothing, (1 / dtθ), tθ)

K = prototype_jacobian(op,u0)
M = prototype_jacobian(op,u0)
stif!(K, u0, u0, nothing, tθ)
mass!(M, u0, u0, nothing, tθ)
# Here you have the mass matrix M

@test (1/dtθ)*M+K ≈ J

# To explore the Sundials solver options, e.g., BE with fixed time step dtd
f_iip = DAEFunction{true}(res!; jac = jac!)#, jac_prototype=J)
# jac_prototype is the way to pass my pre-allocated jacobian matrix
prob_iip = DAEProblem{true}(f_iip, u0, u0, tspan, differential_vars = [true,true,true,true])
# When I pass `jac_prototype` the code get stuck here:
# sol_iip = Sundials.solve(prob_iip, IDA(), reltol = 1e-8, abstol = 1e-8)
# @show sol_iip.u

# or iterator version
# integ = init(prob_iip, IDA(), reltol = 1e-8, abstol = 1e-8)
# step!(integ)

# Show using integrators as iterators
# for i in take(integ, 100)
  # @show integ.u
# end

end # module

# FUTURE WORK: Check other options, not only Sundials

# ISSUE: Future work, add own (non)linear solver.
# Let us assume that we just want to consider a fixed point algorithm
# and we consider an implicit time integration of a nonlinear PDE.
# Our solvers are efficient since they re-use cache among
# nonlinear steps (e.g., symbolic factorization, arrays for numerical factorization)
# and for the transient case in our library, we can also reuse all this
# between time steps. Could we attain something like this using
# DifferentialEquations/Sundials?
# @ChrisRackauckas suggests to take a look at:
# https://docs.sciml.ai/latest/tutorials/advanced_ode_example/ shows swapping out linear solvers.
# https://docs.sciml.ai/latest/features/linear_nonlinear/ is all of the extra details.

# Try to pass solver too
# ls = LUSolver()
# ls_cache = nothing
# x = copy(u0)
# solve!(x,J,r,ls_cache)
#


text:
module NondimLinearOneHalfLayerQG

export
  Problem,
  set_q!,
  updatevars!,
  get_bickley,
  get_σ,
  streamfunctionfrompv!,
  energies

using
  FFTW,
  Reexport,
  DocStringExtensions

@reexport using FourierFlows

using LinearAlgebra: mul!, ldiv!
using FourierFlows: parsevalsum

"""
    Problem(dev::Device=CPU();
                nx = 256,
                ny = nx,
                Lx = 2.0,
                Ly = Lx,
                x0 = 0.0,
                 F = 25.0,
                 U = nothing,
               Uyy = nothing,
                 β = 0.0,
                dt = 2e-3,
           stepper = "RK4",
            effort = FFTW.MEASURE,
  aliased_fraction = 0,
                 T = Float64)
Construct a two-dimensional nondimensional linear two-layer QG `problem` on device `dev`.
Keyword arguments
=================
  - `dev`: (required) `CPU()` or `GPU()`; computer architecture used to time-step `problem`.
  - `nx`: Number of grid points in ``x``-domain.
  - `ny`: Number of grid points in ``y``-domain.
  - `Lx`: Extent of the ``x``-domain.
  - `Ly`: Extent of the ``y``-domain.
  - `x0`: Origin of the ``x``-domain.
  - `F`: Inverse Burger number.
  - `U`: Imposed zonal flow U(y).
  - `Uyy`: Meridional relative vorticity PV gradient.
  - `β`: Meridional linear planetary PV gradient.
  - `dt`: Time-step.
  - `stepper`: Time-stepping method.
  - `effort`: FFTW effort.
  - `aliased_fraction`: the fraction of high-wavenumbers that are zero-ed out by `dealias!()`.
  - `T`: `Float32` or `Float64`; floating point type used for `problem` data.
"""
function Problem(dev = CPU();
    # Numerical parameters
                  nx = 256,
                  ny = nx,
                  Lx = 2.0,
                  Ly = Lx,
                  x0 = 0.0,
    # Physical parameters
                   F = 25.0,     # inverse Burger number
                   U = nothing,  # imposed zonal flow U(y)
                 Uyy = nothing,  # meridional relative vorticity PV gradient
                   β = 0.0,      # meridional linear planetary PV gradient
    # Timestepper and equation options
                  dt = 2e-3,
             stepper = "RK4",
              effort = FFTW.MEASURE,
  # Float type and dealiasing
    aliased_fraction = 0.0,
                   T = Float64)

    # The grid.
    grid = TwoDGrid(dev, nx, Lx, ny, Ly; x0=x0, aliased_fraction=aliased_fraction, T=T, effort=effort)

    # No mean flow case.
    U === nothing && (U = zeros(dev, T, (1)); Uyy = zeros(dev, T, (1)))

    params = Params(F, U, Uyy, β, grid, dev=dev)
    vars = Vars(dev, grid)
    equation = Equation(dev, params, grid)

    FourierFlows.Problem(equation, stepper, dt, grid, vars, params, dev)
end

abstract type NondimLinearOneHalfLayerQGParams <: AbstractParams end

"""
    Params{T, Aphys2D}(F, U, Uyy, β, grid, Qy)
A struct containing the parameters for the NondimLinearOneHalfLayerQG problem. Included are:
$(TYPEDFIELDS)
"""
struct Params{T, Aphys2D} <: NondimLinearOneHalfLayerQGParams
    # prescribed params
    "inverse Burger number"
         F :: T
    "background zonal flow U(y)"
         U :: Aphys2D
    "background PV gradient due to zonal flow's lateral shear"
       Uyy :: Aphys2D
     "planetary vorticity gradient (in the y direction)"
         β :: T

    # derived param
    "background potential vorticity gradient in y-direction due to β and U"
        Qy :: Aphys2D
end

function convert_U_to_U3D(dev, grid, U::AbstractArray{TU, 1}) where TU
  T = eltype(grid)
  if length(U) == 1
    U_2D = zeros(dev, T, (1, 2))
    U_2D[:] = U
    U_2D = repeat(U_2D, outer=(grid.ny, 1))
  else
    U_2D = zeros(dev, T, (grid.ny, 1))
    U_2D[:] = U
  end
  U_3D = zeros(dev, T, (1, grid.ny))
  @views U_3D[1, :, :] = U_2D
  return U_3D
end

function convert_U_to_U3D(dev, grid, U::AbstractArray{TU, 2}) where TU
  T = eltype(grid)
  U_3D = zeros(dev, T, (1, grid.ny))
  @views U_3D[1, :, :] = U
  return U_3D
end

function Params(F, U, Uyy, β, grid; dev::Device=CPU()) where TU
    T = eltype(grid)
    A = ArrayType(dev)

    U = convert_U_to_U3D(dev, grid, U)
    Uyy = convert_U_to_U3D(dev, grid, Uyy)

    ny, nx = grid.ny, grid.nx
    Krsq = grid.Krsq

    Qy = zeros(dev, T, (nx, ny))
    @. Qy = F*U - Uyy + β

    return Params(F, U, Uyy, β, Qy)
end

"""
    Equation(params::NondimLinearOneHalfLayerQGParams, grid)
Return the `equation` for a nondimensional one and a half-layer QG problem with `params` and `grid`. Linear operator
``L`` is empty in this case (no hypo-/hyer-viscosity, and β terms are absorbed in the background PV gradient term.)
Nonlinear term is computed via `calcN!` function.
"""
function Equation(dev, params, grid)
  T = eltype(grid)
  L = zeros(dev, T, (grid.nkr, grid.nl))

  return FourierFlows.Equation(L, calcN!, grid)
end

abstract type NondimLinearOneHalfLayerQGVars <: AbstractVars end

"""
    Vars{Aphys, Atrans, F, P}(q, ψ, u, v, qh, ψh, uh, vh, prevsol)
The variables for NondimLinearOneHalfLayer QG:
$(FIELDS)
"""
struct Vars{Aphys, Atrans, P} <: NondimLinearOneHalfLayerQGVars
    "potential vorticity"
        q :: Aphys
    "streamfunction"
        ψ :: Aphys
    "x-component of velocity"
        u :: Aphys
    "y-component of velocity"
        v :: Aphys
    "Fourier transform of potential vorticity"
       qh :: Atrans
    "Fourier transform of streamfunction"
       ψh :: Atrans
    "Fourier transform of x-component of velocity"
       uh :: Atrans
    "Fourier transform of y-component of velocity"
       vh :: Atrans
    "`sol` at previous time-step"
  prevsol :: P
end

"""
    Vars(dev, grid)
Return the `vars` for for the NondimLinearOneHalfLayerQG problem on device `dev` and
with `grid`.
"""
function Vars(::Dev, grid::AbstractGrid) where Dev
  T = eltype(grid)

  @devzeros Dev T (grid.nx, grid.ny) q ψ u v
  @devzeros Dev Complex{T} (grid.nkr, grid.nl) qh ψh uh vh prevsol

  return Vars(q, ψ, u, v, qh, ψh, uh, vh, prevsol)
end

"""
    calcN!(N, sol, t, clock, vars, params, grid)
Calculate the nonlinear term, that is the advection term,
```math
N = - \\widehat{𝖩(ψ, q+η)}.
```
"""
function calcN!(N, sol, t, clock, vars, params, grid)
    @. vars.qh = sol
    streamfunctionfrompv!(vars.ψh, vars.qh, params, grid)
    @. vars.uh = -im * grid.l  * vars.ψh
    @. vars.vh =  im * grid.kr * vars.ψh

    qx, qxh = deepcopy(vars.q), deepcopy(vars.qh)
    @. qxh = im * grid.kr * vars.qh

    ldiv!(vars.q, grid.rfftplan, deepcopy(vars.qh))
    ldiv!(vars.u, grid.rfftplan, deepcopy(vars.uh))
    ldiv!(vars.v, grid.rfftplan, deepcopy(vars.vh))
    ldiv!(qx, grid.rfftplan, deepcopy(qxh))

    # use vars.u, vars.v, vars.uh, vars.vh as scratch variables
    rhs1, rhs1h = deepcopy(vars.u), deepcopy(vars.uh)
    rhs2, rhs2h = deepcopy(vars.v), deepcopy(vars.vh)
    rhs1 = @. params.U*qx
    rhs2 = @. vars.v*params.Qy
    mul!(rhs1h, grid.rfftplan, deepcopy(rhs1)) # \hat{U*qx}
    mul!(rhs2h, grid.rfftplan, deepcopy(rhs2)) # \hat{v*Qy}

    @. N = - (rhs1h + rhs2h)
    return nothing
end

"""
    streamfunctionfrompv!(ψh, qh, params::NondimLinearOneHalfLayerQGParams, grid)
Invert the PV to obtain the Fourier transform of the streamfunction `ψh`.
"""
function streamfunctionfrompv!(ψh, qh, params::NondimLinearOneHalfLayerQGParams, grid)
  @. ψh  = - qh / (grid.Krsq + params.F)

  return nothing
end

# ----------------
# Helper functions
# ----------------

"""
    updatevars!(vars, params, grid, sol)
    updatevars!(prob)
Update all problem variables using `sol`.
"""
function updatevars!(vars, params, grid, sol)
  @. vars.qh = sol
  streamfunctionfrompv!(vars.ψh, vars.qh, params, grid)
  @. vars.uh = -im * grid.l  * vars.ψh
  @. vars.vh =  im * grid.kr * vars.ψh

  ldiv!(vars.q, grid.rfftplan, deepcopy(vars.qh))
  ldiv!(vars.ψ, grid.rfftplan, deepcopy(vars.ψh))
  ldiv!(vars.u, grid.rfftplan, deepcopy(vars.uh))
  ldiv!(vars.v, grid.rfftplan, deepcopy(vars.vh))

  return nothing
end

updatevars!(prob) = updatevars!(prob.vars, prob.params, prob.grid, prob.sol)

"""
    set_q!(sol, params, vars, grid, q)
    set_q!(prob, q)
Set the solution `prob.sol` as the transform of `q` and update variables.
"""
function set_q!(sol, params, vars, grid, q)
  A = typeof(vars.q)
  mul!(vars.qh, grid.rfftplan, A(q))

  vars.qh[1, 1] = 0
  @. sol = vars.qh

  updatevars!(vars, params, grid, sol)

  return nothing
end

set_q!(prob, q) = set_q!(prob.sol, prob.params, prob.vars, prob.grid, q)

"""
    get_bickley(grid, U₀, Lj)
Get a Bickley jet of half-width `Lj` and amplitude `U₀` on grid `grid`.
"""
function get_bickley(grid, U₀, Lj)
    T = eltype(grid)
    y, ny = grid.y, grid.ny
    U, Uyy = zeros(T, ny), zeros(T, ny)

    @. U = U₀*sech(y/Lj)^2
    @. Uyy = (U₀*2/Lj^2)*U*(3*tanh(y/Lj)^2 - 1)

    return U, Uyy
end

"""
    get_σ(t, KE, PE)
Estimate the growth rate of instabilities from a least-squares fit to energy
histories. `KE` is the kinetic energy and `PE` is the potential energy.
"""
function get_σ(t, KE, PE)
  f = findlast(diff(KE + PE).<-0.5) + 1 # Get first time index from the last renormalization cycle.
  t, KE = t[f:end], KE[f:end]
  n = size(t)[1]
  d = Matrix(reshape(log.(KE), (1, n)))
  gm = Matrix(reshape(t, (1, n)))
  Gm = Matrix([ones(n, 1) gm'])
  GmT = Gm'
  mv = inv(GmT*Gm)*(GmT*d')
  σ = mv[2]

  return σ
end

"""
    energies(vars, params, grid, sol)
    energies(prob)
Return the kinetic energy KE, and the potential energy PE.
The kinetic energy is
```math
𝖪𝖤 = \\int \\frac{1}{2} |{\\bf ∇} ψ|^2 dx dy = \\frac{1}{2} \\sum_{𝐤} |𝐤|² |ψ̂|²,
```
while the potential energy that corresponds to the interface (i.e., the interface between the active and inert layers) is
```math
𝖯𝖤 = \\int \\frac{1}{2} F ψ^2 dx dy = \\frac{1}{2} F \\sum_{𝐤} |ψ|².
```
"""
function energies(vars, params::NondimLinearOneHalfLayerQGParams, grid, sol)
  KE, PE = zeros(1), zeros(1)

  @. vars.qh = sol
  streamfunctionfrompv!(vars.ψh, vars.qh, params, grid)

  abs²∇𝐮h = vars.uh        # use vars.uh as scratch variable
  @. abs²∇𝐮h = grid.Krsq * abs2(vars.ψh)

  KE[1] = @views 0.5 * parsevalsum(abs²∇𝐮h, grid)
  PE[1] = @views 0.5 * params.F * parsevalsum(abs2.(vars.ψh), grid)

  return KE, PE
end

energies(prob) = energies(prob.vars, prob.params, prob.grid, prob.sol)

end # module

max_stars_repo_name: apaloczy/RoughOceanQGEnergetics

max_issues_repo_name: apaloczy/RoughOceanQGEnergetics

max_forks_repo_name: apaloczy/RoughOceanQGEnergetics


text:
# this is going to hold the main updating functions e.g. upduate activation etc.
function update_activation!(t, layer, nn, matrices, parameters)

	vt, vr, C, a, b, c, d, k, σ, δt = get_parameters(update_activation!,parameters,layer)
	v, sp, spt, rec, I = get_matrices(update_activation!, matrices, layer)

	ξ = generate_noise(σ,nn[layer])

	update_voltage!(v, vr, vt, rec, I, ξ, C, k, δt)
	update_recovery!(v, vr, rec, a, b, δt)
	update_spikes!(v, sp, spt, t, vt, rec, c, d)
end

function generate_noise(σ,nn)

	(rand(nn) .- 0.5) .* σ
end

function update_voltage!(v, vr, vt, rec, I, ξ, C, k, δt)

	δv = (k .* (v .- vr) .* (v .- vt) .- rec .+ I .+ ξ) ./ C

	v .= v .+ δv .* δt
end

function update_recovery!(v, vr, rec, a, b, δt)

	δrec = a .* (b .* (v ./ vr) .- rec)

	rec .= rec + δrec .* δt
end

function update_spikes!(v, sp, spt, t, vt, rec, c, d)

	# check spikes
	sp .= v .> vt
	spt[sp] .= t
	v .= v .* (.! sp) .+ c .* sp # turn all into c if above resting
	rec .= rec .+ (.! sp) .* d
end

function calc_input!(l, m, p; rev=0)

	input, w, ACh, v = get_matrices(calc_input!, m, l)

	l==2 && winner_takes_all(ACh,p.winners[l])
	input .= dropdims(sum( @.(w * ACh * (rev - v)'), dims=1), dims=1) #dropdims(sum(w .* ACh .* (rev .- v)', dims = 1), dims=1)

	l==1 && winner_takes_all(input, p.winners[l])
end

function divisive_normalise!(input, inmax)
	if sum(input) == 0
		return
	end
	input .= inmax * input ./ sqrt(sum(input .^ 2))
end

function subtractive_normalise(input, λ)

	input = - λ * sum(input) / length(input)
end

function winner_takes_all(input, numwinners)
	if length(input) == numwinners
		return
	end
	anti_assign(input, maxk(input,numwinners), 0.0)
end

function maxk!(ix, a, k; initialized=false)
    partialsortperm!(ix, a, 1:k, rev=true, initialized=initialized)
    return ix
end

function maxk(a,k)

	partialsortperm(a,1:k,rev=true)
end

function anti_assign(array, indices, assign_otherwise)

	for i in 1:length(array)
		array[i] = i in indices ? array[i] : assign_otherwise
	end
end


text:
user_dir = joinpath(dirname(@__FILE__), "..", "myMatrixDepot")
dirdata = joinpath(dirname(@__FILE__),"..","data")
 if isdir(user_dir)
     rm(user_dir, recursive = true)
end
if isdir(string(dirdata, '/', "uf"))
    rm(string(dirdata, '/', "uf"), recursive = true)
end
if isdir(string(dirdata, '/', "mm"))
    rm(string(dirdata, '/', "mm"), recursive = true)
end
workspace()
using MatrixDepot
using Base.Test


text:
using LinearAlgebra: norm
using PSSFSS
using PSSFSS.FillZY: fillz, filly
using PSSFSS.Elements: s₁s₂2β₁β₂
using PSSFSS.Sheets: Sheet
using PSSFSS.RWG: setup_rwg
using Unitful: ustrip
using Logging: Error, ConsoleLogger, default_metafmt, global_logger
using Test

testlogger = ConsoleLogger(stderr, Error,
                       meta_formatter=default_metafmt, show_limited=true,
                       right_justify=0)
oldlogger = global_logger(testlogger)

@testset "fillz" begin

    FGHz = 2.0; wl_inch = 11.80285 / FGHz
    inch2m = 2.54/100
    wl_m = wl_inch * inch2m
    k0 = 2π / wl_m
    strata = [ Layer()
               Layer(ϵᵣ=3.4, tanδ=0.003, width=1mil)
               rectstrip(Nx=20, Ny=20, Lx=20, Ly=50, Px=100, Py=100, units=mil)
               Layer(ϵᵣ=3.4, tanδ=0.003, width=1mil)
             ]

    layers = [s for s in strata if s isa Layer]

    s = findfirst(x -> x isa Sheet, strata) - 1 # location of Sheet
    metal = [x for x in strata if x isa Sheet][1]
    β₁, β₂ = metal.β₁, metal.β₂
    ufact = 0.5
    units_per_meter = ustrip(Float64, metal.units, 1u"m")
    u = ufact * max(norm(β₁), norm(β₂)) * units_per_meter

    rwgdat = setup_rwg(metal)
    ψ₁ = ψ₂ = 0.0

    zmat = fillz(k0,u,layers,s,ψ₁,ψ₂,metal,rwgdat)

    @test abs(zmat[20,30] - (0.8990238712012797 - 354.1682770720877im)) < 0.005
    @test abs(zmat[100,100] - (110.95156005861463 - 37529.820828933574im)) < 0.01

end

@testset "filly" begin

    FGHz = 2.0; wl_inch = 11.80285 / FGHz
    inch2m = 2.54/100
    wl_m = wl_inch * inch2m
    k0 = 2π / wl_m
    strata = [ Layer()
              Layer(ϵᵣ=3.4, tanδ=0.003, width=1mil)
              rectstrip(class='M', Nx=20, Ny=20, Lx=20, Ly=50, Px=100, Py=100, units=mil)
              Layer(ϵᵣ=3.4, tanδ=0.003, width=1mil)
            ]

    layers = [s for s in strata if s isa Layer]

    s = findfirst(x -> x isa Sheet, strata) - 1 # location of Sheet
    apert = [x for x in strata if x isa Sheet][1]
    β₁, β₂ = apert.β₁, apert.β₂
    ufact = 0.5
    units_per_meter = ustrip(Float64, apert.units, 1u"m")
    u = ufact * max(norm(β₁), norm(β₂)) * units_per_meter

    rwgdat = setup_rwg(apert)
    ψ₁ = ψ₂ = 0.0

    ymat = filly(k0,u,layers,s,ψ₁,ψ₂,apert,rwgdat)

    @test abs(ymat[20,30] - (1.2168396602563192e-6 - 0.02476973884018213im)) < 0.0001
    @test abs(ymat[100,100] - (2.2127855040521283e-6 - 3.5131557304457406im)) < 0.0001

end

global_logger(oldlogger)

nothing


text:
struct SimplePar <: AbstractComponent
    par::Parameter
    SimplePar(val::Number) = new(Parameter(val))
end

prepare!(comp::SimplePar, domain::AbstractDomain) = [NaN]

function evaluate!(buffer, comp::SimplePar, domain::AbstractDomain,
                   par)
    buffer[1] = par
end


text:
using SafeTestsets

@time @safetestset "neighborhoods" begin include("neighborhoods.jl") end
@time @safetestset "mask" begin include("mask.jl") end
@time @safetestset "framework" begin include("framework.jl") end
@time @safetestset "outputs" begin include("outputs.jl") end
@time @safetestset "utils" begin include("utils.jl") end
@time @safetestset "integration" begin include("integration.jl") end
@time @safetestset "image" begin include("image.jl") end
@time @safetestset "multi" begin include("multi.jl") end


text:
using Random, BenchmarkTools, SNOPT7
using PGFPlots
const PGF = PGFPlots

# Car escape
T = Float64

# options
max_con_viol = 1.0e-8
verbose=false

opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,live_plotting=:off)

opts_al = AugmentedLagrangianSolverOptions{T}(verbose=verbose,
    opts_uncon=opts_ilqr,
    cost_tolerance=1.0e-6,
    cost_tolerance_intermediate=1.0e-2,
    constraint_tolerance=max_con_viol,
    penalty_scaling=50.,
    penalty_initial=10.)

opts_pn = ProjectedNewtonSolverOptions{T}(verbose=verbose,
    feasibility_tolerance=max_con_viol,
    solve_type=:feasible)

opts_altro = ALTROSolverOptions{T}(verbose=verbose,
    opts_al=opts_al,
    R_inf=1.0e-1,
    resolve_feasible_problem=false,
    opts_pn=opts_pn,
    projected_newton=true,
    projected_newton_tolerance=1.0e-3);

opts_ipopt = DIRCOLSolverOptions{T}(verbose=verbose,
    nlp=Ipopt.Optimizer(),
    feasibility_tolerance=max_con_viol)

opts_snopt = DIRCOLSolverOptions{T}(verbose=verbose,
    nlp=SNOPT7.Optimizer(),
    feasibility_tolerance=max_con_viol)

x0 = Problems.car_escape.x0
xf = Problems.car_escape.xf

# ALTRO

prob_altro = copy(Problems.car_escape)
@time p1, s1 = solve(prob_altro, opts_altro)
@benchmark p1, s1 = solve($prob_altro, $opts_altro)
max_violation_direct(p1)
Problems.plot_escape(p1.X,x0,xf)

# DIRCOL w/ Ipopt
prob_ipopt = update_problem(copy(Problems.car_escape),model=Dynamics.car) # get continuous time model
@time p2, s2 = solve(prob_ipopt, opts_ipopt)
@benchmark p2, s2 = solve($prob_ipopt, $opts_ipopt)
max_violation_direct(p2)
Problems.plot_escape(p2.X,x0,xf)

# DIRCOL w/ SNOPT
prob_snopt = update_problem(copy(Problems.car_escape),model=Dynamics.car) # get continuous time model
@time p3, s3 = solve(prob_snopt, opts_snopt)
@benchmark p3, s3 = solve($prob_snopt, $opts_snopt)
max_violation_direct(p3)
Problems.plot_escape(p3.X,x0,xf)

# AL-iLQR
opts_ilqr = iLQRSolverOptions{T}(verbose=true,live_plotting=:off,iterations=150)

opts_al = AugmentedLagrangianSolverOptions{T}(verbose=true,
    opts_uncon=opts_ilqr,
    cost_tolerance=1.0e-6,
    cost_tolerance_intermediate=1.0e-2,
    constraint_tolerance=max_con_viol,
    penalty_scaling=50.,
    penalty_initial=10.,
    iterations=9)

prob_altro = copy(Problems.car_escape)
prob_altro.X .*= NaN

@time p4, s4 = solve(prob_altro, opts_al)
max_violation(p4)
Problems.plot_escape(p4.X,x0,xf)
x = [x[1] for x in p4.X]
y = [x[2] for x in p4.X]
PGF.Plots.Linear(x,y);

###############################################
#            Create PGF Plot                  #
###############################################

# Plot the walls
color_wall = "gray"
style = "color=$color_wall, fill=$color_wall"
p = [PGF.Plots.Circle(circle..., style=style) for circle in Problems.circles_escape]

# Plot the trajectories
t0 = PGF.Plots.Linear(Problems.X0_escape[1,:],Problems.X0_escape[2,:],
    legendentry="initial guess",
    mark="none",
    style="color=black, dashed, thick")

t1 = trajectory_plot(p1, mark="none", legendentry="ALTRO", style="very thick, color=$col_altro, dashed")
t2 = trajectory_plot(p1, mark="none", legendentry="Ipopt", style="very thick, color=$col_ipopt")
t3 = trajectory_plot(p3, mark="none", legendentry="SNOPT", style="very thick, color=$col_snopt");
t4 = trajectory_plot(p4, mark="none", legendentry="AL-iLQR", style="very thick, color=cyan");

# Plot the initial and final points
goal = ([x0[1], xf[1]],
        [x0[2], xf[2]])

z = ["start", "end"]
sc = "{start={mark=*,yellow, scale=1.5, mark options={fill=yellow}},end={mark=square*,red,scale=1.5, mark options={fill=red}}}"
g = PGF.Plots.Scatter(goal[1], goal[2], z, scatterClasses=sc)

# Plot the whole thing
a = Axis([p;t0; t3; t2; t4; t1; g],
    xmin=-1, ymin=-1, xmax=11, ymax=8,
    axisEqualImage=true,
    legendPos="north west",
    hideAxis=true)

# Save to tikz format
# paper = "/home/taylor/Documents/research/ALTRO_paper/images"
save(joinpath(paper,"escape_traj.tikz"), a, include_preamble=false)

# Max constraint plot
t_pn = s1.stats[:time_al]
t_span_al = range(0,stop=s1.stats[:time_al],length=s1.solver_al.stats[:iterations])
t_span_pn = range(t_pn,stop=s1.stats[:time],length=s1.solver_pn.stats[:iterations]+1)
t_span = [t_span_al;t_span_pn[2:end]]
c_span = [s1.solver_al.stats[:c_max]...,s1.solver_pn.stats[:c_max]...]

s_range = [1,(2:1:length(s3.stats[:iter_time])-1)...,length(s3.stats[:iter_time])]
i_range = [1,(2:1:length(s2.stats[:iter_time])-1)...,length(s2.stats[:iter_time])]
p = plot(t_pn*ones(100),range(1.0e-9,stop=1.0,length=100),color=:red,linestyle=:dash,label="Projected Newton",width=2)
p = plot!(s3.stats[:iter_time][s_range],s3.stats[:c_max][s_range],marker=:circle,yscale=:log10,ylim=[1.0e-9,1.0],color=:green,label="SNOPT")
p = plot!(s2.stats[:iter_time][i_range],s2.stats[:c_max][i_range],marker=:circle,yscale=:log10,ylim=[1.0e-9,1.0],color=:blue,label="Ipopt")
p = plot!(t_span,c_span,title="Car Escape c_max",xlabel="time (s)",marker=:circle,color=:orange,width=2,yscale=:log10,ylim=[1.0e-9,1.0],label="ALTRO")

savefig(p,joinpath(pwd(),"examples/IROS_2019/car_escape_c_max.png"))

#PGFPlots version
w1 = PGF.Plots.Linear(s3.stats[:iter_time][s_range],s3.stats[:c_max][s_range],mark="none", legendentry="DIRCOL-S",style="very thick, color=green")
w2 = PGF.Plots.Linear(s2.stats[:iter_time][i_range],s2.stats[:c_max][i_range], mark="none",legendentry="DIRCOL-I",style="very thick, color=cyan")
w3 = PGF.Plots.Linear(t_span,c_span, mark="none", legendentry="ALTRO",style="very thick, color=orange!80!black")
w4 = PGF.Plots.Linear(t_pn*ones(100),range(1.0e-9,stop=1.0,length=100),legendentry="Projected Newton",mark="none",style="very thick, color=red, dashed")

a = Axis([w4;w1;w2;w3],
    xmin=0., ymin=1e-9, xmax=45, ymax=1.0,
    legendPos="north east",
    ymode="log",
    hideAxis=false,
    xlabel="time (s)")

paper = "/home/taylor/Documents/research/ALTRO_paper/images"
save(joinpath(paper,"escape_c_max.tikz"), a, include_preamble=false)

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
using PhaseFromInterferograms
using Test

@testset "PhaseFromInterferograms.jl" begin
    # Write your tests here.
end


text:
module Sim
include("./name2idx/parameters.jl")
include("./name2idx/species.jl")
include("./set_model.jl")
include("./observable.jl")
include("./load_csv.jl")

using .C
using .V

using Sundials
using SteadyStateDiffEq

# Options for ODE solver
const ABSTOL = 1e-8
const RELTOL = 1e-8

normalization = Dict{String,Dict{}}()

const dt = 1.0
const t = collect(0.0:dt:120.0)

const conditions = [
    "MCF7_EGF", "MCF7_HRG",
    "BT474_EGF", "BT474_HRG",
    "MDAMB231_EGF", "MDAMB231_HRG",
    "SKBR3_EGF", "SKBR3_HRG"
]

simulations = Array{Float64,3}(
    undef, length(observables), length(t), length(conditions)
)

function solveode(
        f::Function,
        u0::Vector{Float64},
        t::Vector{Float64},
        p::Vector{Float64})::Union{ODESolution{},Nothing}
    local sol::ODESolution{}, is_successful::Bool
    try
        prob = ODEProblem(f, u0, (t[1], t[end]), p)
        sol = solve(
            prob,CVODE_BDF(),
            abstol=ABSTOL,reltol=RELTOL,saveat=dt,dtmin=eps(),verbose=false
        )
        is_successful = ifelse(sol.retcode === :Success, true, false)
    catch
        is_successful = false
    finally
        if !is_successful
            GC.gc()
        end
    end
    return is_successful ? sol : nothing
end

function get_steady_state(
        f::Function,
        u0::Vector{Float64},
        p::Vector{Float64})::Vector{Float64}
    local sol::SteadyStateSolution{}, is_successful::Bool
    try
        prob = ODEProblem(diffeq, u0, (0.0, Inf), p)
        prob = SteadyStateProblem(prob)
        sol = solve(
            prob,
            DynamicSS(
                CVODE_BDF();abstol=ABSTOL,reltol=RELTOL
            ),
            dt=dt,verbose=false
        )
        is_successful = ifelse(sol.retcode === :Success, true, false)
    catch
        is_successful = false
    finally
        if !is_successful
            GC.gc()
        end
    end
    return is_successful ? sol.u : []
end

function simulate!(p::Vector{Float64}, u0::Vector{Float64})::Union{Bool,Nothing}
    p_in_mcf7::Vector{Float64} = copy(p)
    u0_in_mcf7::Vector{Float64} = copy(u0)
    # add ligand
    @inbounds for (i, condition) in enumerate(conditions)
        if condition == "MCF7_EGF"
            p .= p_in_mcf7
            u0 .= u0_in_mcf7
            # u0 = get_steady_state(diffeq,u0,p)
            if isempty(u0)
                return false
            else
                u0[V.EGF] = 10.0
                u0[V.HRG] = 0.0
            end
        elseif condition == "MCF7_HRG"
            p .= p_in_mcf7
            u0 .= u0_in_mcf7
            # u0 = get_steady_state(diffeq,u0,p)
            if isempty(u0)
                return false
            else
                u0[V.EGF] = 0.0
                u0[V.HRG] = 10.0
            end
        elseif condition == "BT474_EGF"
            p .= p_in_mcf7
            u0 .= u0_in_mcf7
            (p, u0) = scale_using_mcf7!(p, u0, "BT474_BREAST")
            # u0 = get_steady_state(diffeq,u0,p)
            if isempty(u0)
                return false
            else
                u0[V.EGF] = 10.0
                u0[V.HRG] = 0.0
            end
        elseif condition == "BT474_HRG"
            p .= p_in_mcf7
            u0 .= u0_in_mcf7
            (p, u0) = scale_using_mcf7!(p, u0, "BT474_BREAST")
            # u0 = get_steady_state(diffeq,u0,p)
            if isempty(u0)
                return false
            else
                u0[V.EGF] = 0.0
                u0[V.HRG] = 10.0
            end
        elseif condition == "MDAMB231_EGF"
            p .= p_in_mcf7
            u0 .= u0_in_mcf7
            (p, u0) = scale_using_mcf7!(p, u0, "MDAMB231_BREAST")
            # u0 = get_steady_state(diffeq,u0,p)
            if isempty(u0)
                return false
            else
                u0[V.EGF] = 10.0
                u0[V.HRG] = 0.0
            end
        elseif condition == "MDAMB231_HRG"
            p .= p_in_mcf7
            u0 .= u0_in_mcf7
            (p, u0) = scale_using_mcf7!(p, u0, "MDAMB231_BREAST")
            # u0 = get_steady_state(diffeq,u0,p)
            if isempty(u0)
                return false
            else
                u0[V.EGF] = 0.0
                u0[V.HRG] = 10.0
            end
        elseif condition == "SKBR3_EGF"
            p .= p_in_mcf7
            u0 .= u0_in_mcf7
            (p, u0) = scale_using_mcf7!(p, u0, "SKBR3_BREAST")
            # u0 = get_steady_state(diffeq,u0,p)
            if isempty(u0)
                return false
            else
                u0[V.EGF] = 10.0
                u0[V.HRG] = 0.0
            end
        elseif condition == "SKBR3_HRG"
            p .= p_in_mcf7
            u0 .= u0_in_mcf7
            (p, u0) = scale_using_mcf7!(p, u0, "SKBR3_BREAST")
            # u0 = get_steady_state(diffeq,u0,p)
            if isempty(u0)
                return false
            else
                u0[V.EGF] = 0.0
                u0[V.HRG] = 10.0
            end
        end
        sol = solveode(diffeq, u0, t, p)
        if sol === nothing
            return false
        else
            @inbounds @simd for j in eachindex(t)
                simulations[observables_index("Phosphorylated_Akt"),j,i] = (
                    sol.u[j][V.AktP]
                )
                simulations[observables_index("Phosphorylated_ERK"),j,i] = (
                    sol.u[j][V.ERKP] + sol.u[j][V.ERKPP]
                )
                simulations[observables_index("Phosphorylated_c-Myc"),j,i] = (
                    sol.u[j][V.cMycP] + sol.u[j][V.cMycPP]
                )
            end
        end
    end
end
end # module


text:
println("===============================================")
println("Testing AFM_processing.jl")
println("===============================================")

function _contact_threshold_f()

    f = [1, 1.5, 2, 2.5, 3];
    d = [2.0, 4, 6, 8, 10];
    index = RHEOSafm.contact_threshold(f, d, (threshold = 2,))

    index == 3
end

@test _contact_threshold_f()

function _contact_hertz_f()
    # d = collect(0:0.1:10);
    # R = 1e-3
    # E = 10

    # F = 4/3 *E *R^0.5 .*d.^1.5;
    # d = collect(0:0.1:10.5)
    # F = vcat(zeros(5),F)
    # index = RHEOSafm.contact_hertz(F,d, (R = R,))
    input_file = joinpath(@__DIR__, "testdata", "AFM_contact_test.txt")
    data = importJPK(input_file, AFM(150e-9), sections = ["extend"]);
    data_contact = contact_point(data, AFM(150e-9), "hertz", (R = 150e-9, s = 0.5));
    length_contact = length(data_contact.t);

    length_contact == 127
end

@test _contact_hertz_f()


text:
function sodium_init()
    ccall((:sodium_init, libsodium), Cint, ())
end

function crypto_box_sealbytes()
    ccall((:crypto_box_sealbytes, libsodium), Csize_t, ())
end

function crypto_box_seal(c, m, mlen, pk)
    ccall(
        (:crypto_box_seal, libsodium),
        Cint,
        (Ptr{Cuchar}, Ptr{Cuchar}, Culonglong, Ptr{Cuchar}),
        c,
        m,
        mlen,
        pk,
    )
end


text:
# -----------------------------------------------
# Generic Regression and Solver routines
# -----------------------------------------------

# using Optim
# include("types.jl")

# #########   GENERIC SOLVER   #########

function solve{T<:Number}(X::Array{T}, Y::Array{T}, regtype::AbstractRegressionType)
    # Test for compatibility and correctedness of input:
    if issingletonarray(X)
        error("X is a singleton Array - Misspecified input data")
    elseif issingletonarray(Y)
        error("Y is a singleton Array - Misspecified input data")
    elseif !(size(X)[1] == size(Y)[1] || size(X)[2] == size(Y)[2])
        error("X and Y seem not to have matching dimension - Misspecified input data")
    #elseif size(X)[2] == size(Y)[2]
    #    X = X'
    #    Y = Y'
    end
    # Dispatch to correct regression function:
    if issubtype(typeof(regtype), LinearRegressionType)
        linregress(X, Y, regtype)
    elseif issubtype(typeof(regtype), NonLinearRegressionType)
        nlinfit(X, Y, regtype)
    end
end


text:
module DiscreteWrapper

using JavaCall

const AbstractIntegerDistribution = @jimport org.hipparchus.distribution.discrete.AbstractIntegerDistribution
const BinomialDistribution = @jimport org.hipparchus.distribution.discrete.BinomialDistribution
const Class = @jimport java.lang.Class
const EnumeratedIntegerDistribution = @jimport org.hipparchus.distribution.discrete.EnumeratedIntegerDistribution
const GeometricDistribution = @jimport org.hipparchus.distribution.discrete.GeometricDistribution
const HypergeometricDistribution = @jimport org.hipparchus.distribution.discrete.HypergeometricDistribution
const IntegerDistribution = @jimport org.hipparchus.distribution.IntegerDistribution
const JString = @jimport java.lang.String
const List = @jimport java.util.List
const Object = @jimport java.lang.Object
const PascalDistribution = @jimport org.hipparchus.distribution.discrete.PascalDistribution
const PoissonDistribution = @jimport org.hipparchus.distribution.discrete.PoissonDistribution
const UniformIntegerDistribution = @jimport org.hipparchus.distribution.discrete.UniformIntegerDistribution
const ZipfDistribution = @jimport org.hipparchus.distribution.discrete.ZipfDistribution

export AbstractIntegerDistribution
export BinomialDistribution
export EnumeratedIntegerDistribution
export GeometricDistribution
export HypergeometricDistribution
export PascalDistribution
export PoissonDistribution
export SaddlePointExpansion
export UniformIntegerDistribution
export ZipfDistribution
export cumulative_probability
export equals
export get_class
export get_exponent
export get_mean
export get_number_of_elements
export get_number_of_successes
export get_number_of_trials
export get_numerical_mean
export get_numerical_variance
export get_pmf
export get_population_size
export get_probability_of_success
export get_sample_size
export get_support_lower_bound
export get_support_upper_bound
export hash_code
export inverse_cumulative_probability
export is_support_connected
export log_probability
export normal_approximate_probability
export notify
export notify_all
export probability
export to_string
export upper_cumulative_probability
export wait

include("abstract_integer_distribution.jl")
include("binomial_distribution.jl")
include("enumerated_integer_distribution.jl")
include("geometric_distribution.jl")
include("hypergeometric_distribution.jl")
include("pascal_distribution.jl")
include("poisson_distribution.jl")
include("saddle_point_expansion.jl")
include("uniform_integer_distribution.jl")
include("zipf_distribution.jl")

end


text:
module TestTransformedTargetModel

using MLJBase
using Test
using ..Models
using StableRNGs
const MMI = MLJBase.MLJModelInterface

rng = StableRNG(698790187)

atom = DeterministicConstantRegressor()
p_atom = ConstantRegressor()
whitener = UnivariateStandardizer()

@testset "constructor and clean!" begin

    @test_logs TransformedTargetModel(atom,
                                              target=UnivariateStandardizer)
    model = @test_logs TransformedTargetModel(atom, target=whitener)
    @test model.model == atom
    @test model.inverse == nothing
    @test model.target == whitener
    @test model.cache

    model = @test_logs TransformedTargetModel(model=atom, target=whitener)
    @test model.model == atom
    @test model.inverse == nothing

    model = @test_logs(
        TransformedTargetModel(atom, target=whitener, inverse=identity))
    @test model.model == atom
    @test model.inverse == identity

    @test_throws(MLJBase.ERR_MODEL_UNSPECIFIED,
                 TransformedTargetModel(target=whitener))

    @test_throws(MLJBase.ERR_TOO_MANY_ARGUMENTS,
                 TransformedTargetModel(atom, whitener))

    @test_throws(MLJBase.err_tt_unsupported(whitener),
                 TransformedTargetModel(whitener, target=whitener))

    model = @test_logs((:warn, MLJBase.WARN_IDENTITY_INVERSE),
                       TransformedTargetModel(p_atom, target=whitener))
    @test model.inverse == identity

    model = @test_logs((:warn, MLJBase.WARN_MISSING_INVERSE),
                       TransformedTargetModel(atom, target=y->log.(y)))
    @test model.inverse == identity

end

# a function for transforming and it's inverse:
f(y) = sin.(y)
g(y) = asin.(y)

# implemented as a static model:
mutable struct Nonlinear <: Static
    λ::Float64 # ignored
end
Nonlinear(; λ=1) = Nonlinear(λ)
MMI.transform(::Nonlinear, _, y) = f(y)
MMI.inverse_transform(::Nonlinear, _, z) = g(z)

# some data:
X, _ = make_regression(5)
z = rand(rng, 5)
y = asin.(z)
@test sin.(y) ≈ z

# average of y on untransformed scale:
avg = mean(y)

# average of y on the non-linear scale defined by f:
avg_nonlinear = g(mean(f(y))) # = g(mean(z))

@testset begin "fit and predict"

    # Remember that `atom` predicts the constant mean of the training
    # target on all new observations. Let's
    # check it's expected behaviour before wrapping:
    fr, _, _ = MMI.fit(atom, 0, X, y)
    @test predict(atom, fr, X) ≈ fill(avg, 5)

    # Test wrapping using f and g:
    model = TransformedTargetModel(atom, target=f, inverse=g)
    fr1, _, _ = MMI.fit(model, 0, X, y)
    @test predict(model, fr1, X) ≈ fill(avg_nonlinear, 5)

    # Test wrapping using a `Static` transformer:
    model = TransformedTargetModel(atom, target=Nonlinear())
    fr1, _, _ = MMI.fit(model, 0, X, y)
    @test predict(model, fr1, X) ≈ fill(avg_nonlinear, 5)

    # Test wrapping using a non-static `Unsupervised` model:
    model = TransformedTargetModel(atom, target=whitener)
    fr1, _, _ = MMI.fit(model, 0, X, y)
    @test predict(model, fr1, X) ≈ fill(avg, 5) # whitener is linear

    # Test with `inverse=identity`:
    model = TransformedTargetModel(atom, target=Nonlinear(), inverse=identity)
    fr1, _, _ = MMI.fit(model, 0, X, y)
    @test predict(model, fr1, X) ≈ fill(mean(z), 5)

    # Test a probablistic model:
    model = TransformedTargetModel(p_atom, target=whitener, inverse=identity)
    fr1, _, _ = MMI.fit(model, 0, X, y)
    yhat = predict(model, fr1, X)
    @test isapprox(first(yhat).μ, 0, atol=1e-15)
end

MMI.iteration_parameter(::Type{DeterministicConstantRegressor}) = :n

@testset "traits" begin
    model = TransformedTargetModel(atom, target=Nonlinear())
    @test input_scitype(model) == input_scitype(atom)
    @test target_scitype(model) == target_scitype(atom)
    @test is_wrapper(model)
    @test iteration_parameter(model) == :(model.n)
    @test package_name(model) == "MLJBase"
    @test occursin("2229d", package_uuid(model))
    @test package_license(model) == "MIT"
    @test package_url(model) == "https://github.com/JuliaAI/MLJBase.jl"
end

@testset "integration" begin
    model = TransformedTargetModel(atom, target=Nonlinear())
    mach = machine(model, X, y)
    fit!(mach, verbosity=0)
    @test predict(mach, X) ≈ fill(avg_nonlinear, 5)
    @test issubset([:model, :target], keys(report(mach)))
    @test issubset([:model, :target], keys(fitted_params(mach)))
    @test fitted_params(mach).model.fitresult ≈ mean(z)
end

mutable struct FooModel46 <: Deterministic
    epochs
end

MMI.fit(::FooModel46, verbosity, X, y) =
    nothing, nothing, (training_losses=ones(length(y)),)
MMI.predict(::FooModel46, fitresult, Xnew) = ones(nrows(Xnew))
MMI.supports_training_losses(::Type{<:FooModel46}) = true
MMI.iteration_parameter(::Type{<:FooModel46}) = :epochs
MMI.training_losses(::FooModel46, report) = report.training_losses

X = rand(5)
y = rand(5)

@testset "training_losses" begin
    atom = FooModel46(10)
    model = TransformedTargetModel(atom, target=Nonlinear())
    @test supports_training_losses(model)
    @test iteration_parameter(model) == :(model.epochs)
    mach = machine(model, X, y)
    fit!(mach, verbosity=0)
    rep = report(mach)
    @test rep.model.training_losses == ones(5)
    @test training_losses(mach) == ones(5)
end

end

true


text:
using GCode
using Test

@testset "GCode.jl" begin
    # Write your tests here.
end


text:
module ITensorGaussianMPS

using ITensors
using ITensors.NDTensors
using LinearAlgebra

import LinearAlgebra: Givens

export slater_determinant_to_mps,
  slater_determinant_to_gmps, hopping_hamiltonian, slater_determinant_matrix

include("gmps.jl")

end


text:
module noise_models
    using Symplectic
    using Cubature: hquadrature

    type SingleAngleModel
        n::Int64
        d_x::Int64
        d_z::Int64
        rate::Float64
        theta::Array{Float64, 1}
        grains::Float64
        S::Array{Float64, 1}
        M::Array{Float64, 1}
        p_theta::Array{Float64, 1}
        theta_hat::Float64
        M_hat::Int64
        theta_true::Float64
        M_true::Int64
        sampling_rates::Dict
        update_functions::Dict

        function SingleAngleModel(p::Float64, grains::Number)
            this = new();
            this.n = 15;
            this.d_x = 3;
            this.d_z = 1;
            this.rate = p;
            this.theta = linspace(0, pi/2, grains);
            this.grains = convert(Float64, grains);
            this.S = linspace(-pi, pi, 4*grains+1);
            this.M = 0.5*(this.S[1:end-1] + this.S[2:end]);
            this.p_theta = log(ones(convert(Int64, grains)) ./ grains);
            this.theta_hat = pi/4;
            this.M_hat = floor(grains/2);
            this.theta_true = (0.05+0.9*rand())*pi/2;
            this.M_true = floor(this.theta_true/(pi/2)*grains);

            # SAMPLING RATES
            P_XZ = zeros(this.n+1, this.n+1, length(this.M));
            COSES = cos(this.M - this.theta_true).^2;
            SINES = sin(this.M - this.theta_true).^2;
            RUNNING_COSES = ones(length(COSES));
            RUNNING_SINES = ones(length(SINES));
            COSESK = ones(this.n+1, length(COSES));
            SINESK = ones(this.n+1, length(COSES));

            tic();
            for k=0:this.n
                COSESK[k+1, :] = deepcopy(RUNNING_COSES);
                SINESK[k+1, :] = deepcopy(RUNNING_SINES);
                RUNNING_COSES = RUNNING_COSES .* COSES;
                RUNNING_SINES = RUNNING_SINES .* SINES;
            end
            t = toq();
            # print("marginal arrays time: $t\n");

            tic();
            BINOMS = zeros(this.n+1, this.n+1);
            for j=0:this.n
                for k=0:(this.n-j)
                    BINOMS[j+1, k+1] = binomial(this.n, j) * binomial(this.n-j, k);
                end
            end
            t = toq();
            # print("binomial array time: $t\n");

            tic()
            RATES = ones(this.n+1, this.n+1);
            for j=0:this.n
                for k=0:(this.n-j)
                    RATES[j+1, k+1] = (this.rate)^(j+k) * (1-this.rate)^(this.n - j - k);
                end
            end
            t = toq();
            # print("rates array time: $t\n");

            COEFFICIENTS = BINOMS .* RATES;

            tic();
            for j=0:1:this.n
                for k=0:(this.n-j)
                    P_XZ[j+1, k+1, :] = COEFFICIENTS[j+1, k+1] * COSESK[j+1, :] .* SINESK[k+1, :];
                end
            end
            t = toq();
            # print("joint arrays time: $t\n");

            tic();
            this.sampling_rates = Dict("W_X>0" => reshape(sum(P_XZ[2:end, :,  :], [1, 2]), length(this.M)),
                                       "W_Z>0" => reshape(sum(P_XZ[:, 2:end, :], [1, 2]), length(this.M)),
                                       "W_X>d_x" => reshape(sum(P_XZ[this.d_x+2:end, :, :], [1, 2]), length(this.M)),
                                       "W_X=k|d_x>=W_X>0" => reshape(sum(P_XZ[2:this.d_x+1, :, :], [2]) ./ sum(P_XZ[2:this.d_x+1, :, :], [1, 2]), this.d_x, length(this.M)),
                                       "W_Z>d_z|W_Z>0" => reshape(sum(P_XZ[:, this.d_z+2:end, :], [1, 2]) ./ sum(P_XZ[:, 2:end, :], [1, 2]), length(this.M)),
                                       "uncorrectable" => reshape(sum(P_XZ[1:this.d_x+1, this.d_z+2:end, :], [1, 2])
                                                           + sum(P_XZ[this.d_x+2:end, 1:this.d_z+1, :], [1, 2])
                                                           + sum(P_XZ[this.d_x+2:end, this.d_z+2:end, :], [1, 2]), length(this.M))
                                    )
            t = toq();
            # print("summations time: $t\n");

            # UPDATE FUNCTIONS
            # Z=1, X=3, p^4 (5x/128 + 1/64*sin(2x) - 1/128 *sin(4x) - 1/192*sin(6x) - 1/1024*sin(8x)) * (1-p)^11
            # Z=1, X=2, p^3 (x/16 + 1/64*sin(2x) - 1/64*sin(4x) - 1/192*sin(6x)) * (1-p)^12
            # Z=1, X=1, p^2 (x/8 - 1/32*sin(4x)) * (1-p)^13
            # Z=0, X=3, p^3 (5x/16 + 15/64*sin(2x) + 3/64*sin(4x) + 1/192*sin(6x)) * (1-p)^12
            # Z=0, X=2, p^2 (3x/8 + 1/4*sin(2x) + 1/32*sin(4x)) * (1-p)^13
            # Z=0, X=1, p^1 (x/2 + 1/4*sin(2x)) * (1-p)^14

            tic();
            p = this.rate;
            PL = Dict("Z=1, X=0" => x -> p^1 * (x/2 - 1/4*sin(2*x)) * (1-p)^14,
                      # "Z=1, X=1" => x -> p^2 * (x/8 - 1/32*sin(4*x)) * (1-p)^13,
                      # "Z=1, X=2" => x -> p^3 * (x/16 + 1/64*sin(2*x) - 1/64*sin(4*x) - 1/192*sin(6*x)) * (1-p)^12,
                      # "Z=1, X=3" => x -> p^4 * (5*x/128 + 1/64*sin(2*x) - 1/128 *sin(4*x) - 1/192*sin(6*x) - 1/1024*sin(8*x)) * (1-p)^11,
                      "Z=0, X=1" => x -> p^1 * (x/2 + 1/4*sin(2*x)) * (1-p)^14,
                      # "Z=0, X=2" => x -> p^2 * (3*x/8 + 1/4*sin(2*x) + 1/32*sin(4*x)) * (1-p)^13,
                      # "Z=0, X=3" => x -> p^3 * (5*x/16 + 15/64*sin(2*x) + 3/64*sin(4*x) + 1/192*sin(6*x)) * (1-p)^12
                    )

            this.update_functions = Dict("Z=1, X=0" => log(PL["Z=1, X=0"](this.S[2:end]) - PL["Z=1, X=0"](this.S[1:end-1])),
                                         # "Z=1, X=1" => PL["Z=1, X=1"](this.S[2:end]) - PL["Z=1, X=1"](this.S[1:end-1]),
                                         # "Z=1, X=2" => PL["Z=1, X=2"](this.S[2:end]) - PL["Z=1, X=2"](this.S[1:end-1]),
                                         # "Z=1, X=3" => PL["Z=1, X=3"](this.S[2:end]) - PL["Z=1, X=3"](this.S[1:end-1]),
                                         "Z=0, X=1" => log(PL["Z=0, X=1"](this.S[2:end]) - PL["Z=0, X=1"](this.S[1:end-1])),
                                         # "Z=0, X=2" => PL["Z=0, X=2"](this.S[2:end]) - PL["Z=0, X=2"](this.S[1:end-1]),
                                         # "Z=0, X=3" => PL["Z=0, X=3"](this.S[2:end]) - PL["Z=0, X=3"](this.S[1:end-1])
                                        )

            update_functions_profile = toq();
            # print("Update functions time: $update_functions_profile\n")
            return this
        end
    end
    export SingleAngleModel

    function rand_error(model::SingleAngleModel, n::Int64, dt::Float64)
        error_types = (rand(1,n) .< sin(dt)^2);
        error_sites = (rand(1,n) .< model.rate);
        symp = Symp((error_sites & error_types), (error_sites & ~error_types));
        return symp
    end
    export rand_error

    function update(model::SingleAngleModel, err::Symp, theta_est::Float64)
        w_x = convert(Float64, sum(err.x));
        w_z = convert(Float64, sum(err.z));
        c_t = cos(model.theta - theta_est).^(2*w_z);
        s_t = sin(model.theta - theta_est).^(2*w_x);
        p_t_given_err = c_t .* s_t .* model.p_theta;
        p_t_given_err = p_t_given_err / sum(p_t_given_err);
        model.p_theta = p_t_given_err;
        return model
    end
    export update

end


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

macro stdout_silent(f)
	ex = quote
		stdouttemp = stdout
		redirect_stdout()
		try
			$f
		catch e
			throw(e)
		finally
			redirect_stdout(stdouttemp)
		end
	end
	return esc(ex)
end


text:
function set_qp1test_solutions!(solver)
    # SOLVE 1
    LQOI.set_solution!(solver,
        variable_primal = [0.571429, 0.428572, 0.857143],
        constraint_primal = [4.0, 1.0],
        variable_dual = [0.0, 0.0, 0.0],
        constraint_dual = [0.714286, 0.857144]
    )
    return
end

function set_qp2test_solutions!(solver)
    # SOLVE 1
    LQOI.set_solution!(solver,
        variable_primal = [0.571429, 0.428572, 0.857143],
        constraint_primal = [4.0, 1.0],
        variable_dual = [0.0, 0.0, 0.0],
        constraint_dual = [0.714286, 0.857144]
    )
    # SOLVE 2
    LQOI.set_solution!(solver,
        variable_primal = [0.571429, 0.428572, 0.857143],
        constraint_primal = [4.0, 1.0],
        variable_dual = [-0.0, -0.0, -0.0],
        constraint_dual = [-1.42857, -1.71429]
    )
    return
end

function set_qp3test_solutions!(solver)
    # SOLVE 1
    LQOI.set_solution!(solver,
        variable_primal = [0.250002, 0.749998],
        constraint_primal = [1.0],
        variable_dual = [6.97141e-6, 4.06505e-8],
        constraint_dual = [2.75]
    )
    # SOLVE 2
    LQOI.set_solution!(solver,
        variable_primal = [1.0, 0.0],
        constraint_primal = [1.0],
        variable_dual = [-0.0, -1.0],
        constraint_dual = [2.0]
    )
    return
end

function set_qcp1test_solutions!(solver)
    # SOLVE 1
    LQOI.set_solution!(solver,
        variable_primal = [0.500003, 1.75],
        constraint_primal = [1.25, 2.25],
        quadratic_primal = [2.0],
        variable_dual = [-0.0, -0.0],
        constraint_dual = [4.48497e-8, 5.48002e-8],
        quadratic_dual = [1.0]
    )
    return
end

function set_qcp2test_solutions!(solver)
    # SOLVE 1
    LQOI.set_solution!(solver,
        variable_primal = [1.41421],
        constraint_primal = Float64[],
        quadratic_primal = [2.0],
        variable_dual = [-0.0],
        constraint_dual = Float64[],
        quadratic_dual = [0.353553]
    )
    return
end

function set_qcp3test_solutions!(solver)
    # SOLVE 1
    LQOI.set_solution!(solver,
        variable_primal = [1.41421],
        constraint_primal = Float64[],
        quadratic_primal = [2.0],
        variable_dual = [0.0],
        constraint_dual = Float64[],
        quadratic_dual = [-0.353553]
    )
    return
end

function set_socp1test_solutions!(solver)
    # SOLVE 1
    LQOI.set_solution!(solver,
        variable_primal = [0.5, 0.5, 0.707107],
        constraint_primal = [1.0],
        quadratic_primal = [-9.55708e-13],
        variable_dual = [-3.72382e-10, -3.72383e-10, 2.4911e-10],
        constraint_dual = [0.707107],
        quadratic_dual = [-1.41421]
    )
    return
end


text:
module HyperPA

using FBD
using BenchmarkTools
using Juno

function source()
    joinpath(@__DIR__, "hyper_pa")
end
function data()
    joinpath(@__DIR__, "..", "data", "hyper_pa")
end

function params(name, nodes)
    edgesize_distribution = FBD.read_probabilities(joinpath(source(), "size distribution", name * " size distribution.txt"))
    degree_distribution = FBD.read_probabilities(joinpath(source(), "simplex per node", name * "-simplices-per-node-distribution.txt"))
    max_edgesize = findlast(edgesize_distribution.accept .> 0)
    degree_distribution, edgesize_distribution, max_edgesize, nodes
end

"""
    jl_with_io(name="DAWN", nodes=3029)

Load probabilities from disk, generate a graph via FBD, save the graph to disk,
and return it.
"""
function jl_with_io(name="DAWN", nodes=3029)
    ps = params(name, nodes)
    graph = hyper_pa(ps...)
    FBD.write_graph(joinpath(data(), "julia", name * ".txt"), graph)
    graph
end

"""
    jl_benchmark(name="DAWN", nodes=3029)

Benchmark the FBD implementation (skip disk IO) and return runtime in seconds
"""
function jl_benchmark(name="DAWN", nodes=3029)
    ps = params(name, nodes)
    m = @benchmark graph = hyper_pa($ps...)
    display(m)
    time(median(m))/1e9
end

"""
    jl_profile(name="DAWN", nodes=3029, runs=20)

Profile the FBD implementation (skip disk IO) and display the results in Juno's GUI
"""
function jl_profile(name="DAWN", nodes=3029, runs=20)
    ps = params(name, nodes)
    Juno.@profiler for i in 1:runs; graph = hyper_pa(ps...); end
end

"""
    external(name="DAWN", nodes=200)

Run the external implementation.
"""
function external(name="DAWN", nodes=200)
    printstyled("expected runtime on the order of $(round(Integer, (nodes/200)^2)) minute(s)\n", color=Base.warn_color())
    output_directory = joinpath(data(), "external")
    @elapsed cd(source()) do
        run(`python3 hyper_preferential_attachment.py --name=$name --file_name=$name --num_nodes=$nodes --simplex_per_node_directory='simplex per node' --size_distribution_directory='size distribution' --output_directory=$output_directory`)
    end
end

#= at default size:
x@X Generator % time python3 hyper_preferential_attachment.py --name=DAWN --file_name=DAWN --num_nodes=3029 --simplex_per_node_directory='simplex per node' --size_distribution_directory='size distribution' --output_directory=output_directory
done with DAWN
python3 hyper_preferential_attachment.py --name=DAWN --file_name=DAWN      22469.39s user 1842.45s system 93% cpu 7:14:08.42 total
x@X Generator %
=#

end