text:
###
### Sampler states
###

struct MH{space, P} <: InferenceAlgorithm
    proposals::P
end

function MH(space...)
    syms = Symbol[]

    prop_syms = Symbol[]
    props = AMH.Proposal[]

    check_support(dist) = insupport(dist, z)

    for s in space
        if s isa Symbol
            push!(syms, s)
        elseif s isa Pair || s isa Tuple
            push!(prop_syms, s[1])

            if s[2] isa AMH.Proposal
                push!(props, s[2])
            elseif s[2] isa Distribution
                push!(props, AMH.Proposal(AMH.Static(), s[2]))
            elseif s[2] isa Function
                push!(props, AMH.Proposal(AMH.Static(), s[2]))
            end
        end
    end

    proposals = NamedTuple{tuple(prop_syms...)}(tuple(props...))
    syms = vcat(syms, prop_syms)
    return MH{tuple(syms...), typeof(proposals)}(proposals)
end

alg_str(::Sampler{<:MH}) = "MH"

#################
# MH Transition #
#################

struct MHTransition{T, F<:AbstractFloat, M<:AMH.Transition}
    θ    :: T
    lp   :: F
    mh_trans :: M
end

function MHTransition(spl::Sampler{<:MH}, mh_trans::AMH.Transition)
    theta = tonamedtuple(spl.state.vi)
    return MHTransition(theta, mh_trans.lp, mh_trans)
end

#####################
# Utility functions #
#####################

"""
    set_namedtuple!(vi::VarInfo, nt::NamedTuple)

Places the values of a `NamedTuple` into the relevant places of a `VarInfo`.
"""
function set_namedtuple!(vi::VarInfo, nt::NamedTuple)
    for (n, vals) in pairs(nt)
        vns = vi.metadata[n].vns

        n_vns = length(vns)
        n_vals = length(vals)
        v_isarr = vals isa AbstractArray

        if v_isarr && n_vals == 1 && n_vns > 1
            for (vn, val) in zip(vns, vals[1])
                vi[vn] = val isa AbstractArray ? val : [val]
            end
        elseif v_isarr && n_vals > 1 && n_vns == 1
            vi[vns[1]] = vals
        elseif v_isarr && n_vals == 1 && n_vns == 1
            if vals[1] isa AbstractArray
                vi[vns[1]] = vals[1]
            else
                vi[vns[1]] = [vals[1]]
            end
        elseif !(v_isarr)
            vi[vns[1]] = [vals]
        else
            error("Cannot assign `NamedTuple` to `VarInfo`")
        end
    end
end

"""
    gen_logπ_mh(vi::VarInfo, spl::Sampler, model)

Generate a log density function -- this variant uses the
`set_namedtuple!` function to update the `VarInfo`.
"""
function gen_logπ_mh(spl::Sampler, model)
    function logπ(x)::Float64
        vi = spl.state.vi
        x_old, lj_old = vi[spl], getlogp(vi)
        # vi[spl] = x
        set_namedtuple!(vi, x)
        runmodel!(model, vi)
        lj = getlogp(vi)
        vi[spl] = x_old
        setlogp!(vi, lj_old)
        return lj
    end
    return logπ
end

"""
    dist_val_tuple(spl::Sampler{<:MH})

Returns two `NamedTuples`. The first `NamedTuple` has symbols as keys and distributions as values.
The second `NamedTuple` has model symbols as keys and their stored values as values.
"""
function dist_val_tuple(spl::Sampler{<:MH})
    vns = _getvns(spl.state.vi, spl)
    dt = _dist_tuple(spl.state.vi.metadata, spl.alg.proposals, spl.state.vi, vns)
    vt = _val_tuple(spl.state.vi.metadata, spl.state.vi, vns)
    return dt, vt
end

@generated function _val_tuple(metadata::NamedTuple, vi::VarInfo, vns::NamedTuple{names}) where {names}
    length(names) === 0 && return :(NamedTuple())
    expr = Expr(:tuple)
    map(names) do f
        push!(expr.args, Expr(:(=), f, :(
            length(metadata.$f.vns) == 1 ? getindex(vi, metadata.$f.vns)[1] : getindex.(Ref(vi), metadata.$f.vns)
        )))
    end
    return expr
end

@generated function _dist_tuple(
    metadata::NamedTuple,
    props::NamedTuple{propnames},
    vi::VarInfo,
    vns::NamedTuple{names}
) where {names, propnames}
    length(names) === 0 && return :(NamedTuple())
    expr = Expr(:tuple)
    map(names) do f
        if f in propnames
            # We've been given a custom proposal, use that instead.
            push!(expr.args, Expr(:(=), f, :(props.$f)))
        else
            # Otherwise, use the default proposal.
            push!(expr.args, Expr(:(=), f, :(AMH.Proposal(AMH.Static(), metadata.$f.dists[1]))))
        end
    end
    return expr
end

#################
# Sampler state #
#################

mutable struct MHState{V<:VarInfo} <: AbstractSamplerState
    vi :: V
    density_model :: AMH.DensityModel
end

###############################
# Static MH (from prior only) #
###############################

function Sampler(
    alg::MH,
    model::Model,
    s::Selector=Selector()
)
    # Set up info dict.
    info = Dict{Symbol, Any}()

    # Make a varinfo.
    vi = VarInfo(model)

    # Make a density model.
    dm = AMH.DensityModel(x -> 0.0)

    # Set up state struct.
    state = MHState(vi, dm)

    # Generate a sampler.
    spl = Sampler(alg, info, s, state)

    # Update the density model.
    spl.state.density_model = AMH.DensityModel(gen_logπ_mh(spl, model))

    return spl
end

function AbstractMCMC.sample_init!(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:MH},
    N::Integer;
    verbose::Bool=true,
    resume_from=nothing,
    kwargs...
)
    # Resume the sampler.
    set_resume!(spl; resume_from=resume_from, kwargs...)

    # Get `init_theta`
    initialize_parameters!(spl; verbose=verbose, kwargs...)
end

function AbstractMCMC.step!(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:MH},
    N::Integer,
    transition;
    kwargs...
)
    if spl.selector.rerun # Recompute joint in logp
        runmodel!(model, spl.state.vi)
    end

    # Retrieve distribution and value NamedTuples.
    dt, vt = dist_val_tuple(spl)

    # Create a sampler and the previous transition.
    mh_sampler = AMH.MetropolisHastings(dt)
    prev_trans = AMH.Transition(vt, getlogp(spl.state.vi))

    # Make a new transition.
    trans = AbstractMCMC.step!(rng, spl.state.density_model, mh_sampler, 1, prev_trans)

    # Update the values in the VarInfo.
    set_namedtuple!(spl.state.vi, trans.params)
    setlogp!(spl.state.vi, trans.lp)

    return Transition(spl)
end

####
#### Compiler interface, i.e. tilde operators.
####
function assume(
    spl::Sampler{<:MH},
    dist::Distribution,
    vn::VarName,
    vi::VarInfo
)
    updategid!(vi, vn, spl)
    r = vi[vn]
    return r, logpdf_with_trans(dist, r, istrans(vi, vn))
end

function dot_assume(
    spl::Sampler{<:MH},
    dist::MultivariateDistribution,
    vn::VarName,
    var::AbstractMatrix,
    vi::VarInfo,
)
    @assert dim(dist) == size(var, 1)
    getvn = i -> VarName(vn, vn.indexing * "[:,$i]")
    vns = getvn.(1:size(var, 2))
    updategid!.(Ref(vi), vns, Ref(spl))
    r = vi[vns]
    var .= r
    return var, sum(logpdf_with_trans(dist, r, istrans(vi, vns[1])))
end
function dot_assume(
    spl::Sampler{<:MH},
    dists::Union{Distribution, AbstractArray{<:Distribution}},
    vn::VarName,
    var::AbstractArray,
    vi::VarInfo,
)
    getvn = ind -> VarName(vn, vn.indexing * "[" * join(Tuple(ind), ",") * "]")
    vns = getvn.(CartesianIndices(var))
    updategid!.(Ref(vi), vns, Ref(spl))
    r = reshape(vi[vec(vns)], size(var))
    var .= r
    return var, sum(logpdf_with_trans.(dists, r, istrans(vi, vns[1])))
end

function observe(
    spl::Sampler{<:MH},
    d::Distribution,
    value,
    vi::VarInfo,
)
    return observe(SampleFromPrior(), d, value, vi)
end

function dot_observe(
    spl::Sampler{<:MH},
    ds::Union{Distribution, AbstractArray{<:Distribution}},
    value::AbstractArray,
    vi::VarInfo,
)
    return dot_observe(SampleFromPrior(), ds, value, vi)
end


text:
using Test
using Logging
using DataStructures
using Dates
using TimeSeries
using InfrastructureSystems
const IS = InfrastructureSystems
using PowerSystems
const PSY = PowerSystems
using PowerSystemCaseBuilder
const PSB = PowerSystemCaseBuilder

@testset "Test Serializtion/De-Serializtion" begin
    system_catelog = SystemCatalog(SYSTEM_CATALOG)
    for (type, descriptor_dict) in system_catelog.data
        for (name, descriptor) in descriptor_dict
            for forecasts in [true, false], reserves in [true, false]
                # build a new system from scratch
                sys = build_system(
                    type,
                    name;
                    add_forecasts = forecasts,
                    add_reserves = reserves,
                    force_build = true,
                )
                @test isa(sys, System)
                # build a new system from json
                @test PSB.is_serialized(name, forecasts, reserves)
                sys2 = build_system(
                    type,
                    name;
                    add_forecasts = forecasts,
                    add_reserves = reserves,
                )
                @test isa(sys2, System)

                PSB.clear_serialized_system(name)
                @test !PSB.is_serialized(name, forecasts, reserves)
            end
        end
    end
end

ext: jl

max_stars_repo_name: NREL-SIIP/PowerSystemBuilder

max_issues_repo_name: NREL-SIIP/PowerSystemBuilder

max_forks_repo_name: NREL-SIIP/PowerSystemBuilder


text:
#   Unit "timefrequency" of the FourierAnalysis Package for julia language
#
#   MIT License
#   Copyright (c) 2019, Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This unit implements time-frequency representations based
#   on analytic signal estimations using FFTW and DSP filters.

#### Time-Frequency Analytic Signal ####

"""
```
(1)
function TFanalyticsignal(x         :: Vector{T},
                          sr        :: Int,
                          wl        :: Int          = 0,
                          bandwidth :: IntOrReal    = 2;
                      fmin          :: IntOrReal    = bandwidth,
                      fmax          :: IntOrReal    = sr÷2,
                      filtkind      :: FilterDesign = Butterworth(2),
                      nonlinear     :: Bool         = false,
                      fsmoothing    :: Smoother     = noSmoother,
                      tsmoothing    :: Smoother     = noSmoother,
                      planner       :: Planner      = getplanner,
                      ⏩           :: Bool         = true) where T<:Real

(2)
function TFanalyticsignal(𝐱         :: Vector{Vector{T}},
                      <same arguments as method (1)>

```
(1)

Given sampling rate `sr`, construct a [TFAnalyticSignal](@ref) object
from univariate data `x`, that is, a time-frequency representation of real
vector `x`.

Call three functions in series performing the three following operations:

i. pass the data throught a bank of pass-band filters calling function
[`filterbank`](@ref),

ii. compute the analytic signal calling function [`analyticsignal`](@ref)
(method (1) therein),

iii. If requested, smooth the analytic signal along the time and/or frequency
dimension calling function [`smooth`](@ref).

The arguments passed to each of these functions are:

|   filterbank   | analyticsignal |   smooth      |
|:--------------:|:--------------:|:-------------:|
|      `x`       |     `wl`       | `fsmoothing`  |
|     `sr`       |  `nonlinear`   | `tsmoothing`  |
|  `bandwidth`   |   `planner`    |               |
|  `filtkind`    |    `⏩`       |               |
|    `fmin`      |                |               |
|    `fmax`      |                |               |
|     `⏩`      |                |               |

Refer to the documentation of each function to learn the meaning
of each argument.

By default the `wl` argument is set to `length(x)`.

If `⏩` is true (default), filtering and computation fo the analytic signal
are multi-threaded across band-pass regions as long as the number of the regions
is at least twice the number of threads Julia is instructed to use.
See [Threads](@ref).

(2)

Given sampling rate `sr`, construct a
[TFAnalyticSignalVector](@ref) object from a vector of univariate data `𝐱`,
that is, from the time-frequency representations of the vectors in `𝐱`.

This method operates as method (1) with the following exceptions:

By default the `wl` argument is set to `length(x)` for each vectors in `𝐱`.
If another values is given, it will be used for all of them.

If `⏩` is true (default), the method is run in multi-threaded mode across the
vectors in `𝐱` as long as their number is at least twice the number of
threads Julia is instructed to use, otherwise this method attempts to run
each analytic signal estimation in multi-threaded mode like in method (1).
See [Threads](@ref).

If a `Planner` is not explicitly passed as an argument,
the FFTW plan is computed once and applied for all analytic signal
estimations.

**See**: [`filterbank`](@ref), [`analyticsignal`](@ref), [`smooth`](@ref).

**Examples**:
```
using Plots, FourierAnalysis

# generate some data
sr, t, bandwidth=128, 512, 2
h=taper(harris4, t)
x1=sinusoidal(10, 8, sr, t, 0)
x2=sinusoidal(10, 19, sr, t, 0)
x=Vector((x1+x2).*h.y+randn(t))
y1=sinusoidal(10, 6, sr, t, 0)
y2=sinusoidal(10, 16, sr, t, 0)
y=Vector((y1+y2).*h.y+randn(t))
plot([x, y])

# TFAnalyticSignal object for x (method (1))
Y = TFanalyticsignal(x, sr, sr*4; fmax=32)

# vector of TFAnalyticSignal objects for x and y (method (2))
𝒀 = TFanalyticsignal([x, y], sr, sr*4; fmax=32)

# (for shortening comments: TF=time-frequency, AS=analytic signal)

# mean AS in a TF region (8:12Hz and samples 1:128) for the first object
m=mean(𝒀[1], (8, 12), (1, 128)) # Output a complex number

# extract the AS in a TF region (8:12Hz and samples 1:128) for the first object
E=extract(𝒀[1], (8, 12), (1, 128)) # Output a complex matrix

# mean AS in a TF region (8:12Hz and samples 1:128) for the two objects
𝐦=mean(𝒀, (8, 12), (1, 128)) # Output a vector of complex numbers
# same computation without checking homogeneity of the two objects in 𝒀 (faster)
𝐦=mean(𝒀, (8, 12), (1, 128); check=false)

# extract the AS in a TF region (8:12Hz and samples 1:128) for the two objects
𝐄=extract(𝒀, (8, 12), (8, 12)) # Output a vector of complex matrices

# plot the real part of the AS of x (see unit recipes.jl)

# gather first useful attributes for the plot
using Plots.Measures
tfArgs=(right_margin = 2mm,
        top_margin = 2mm,
        xtickfont = font(10, "Times"),
        ytickfont = font(10, "Times"))

heatmap(tfAxes(Y)..., real(Y.y);
        c=:pu_or, tfArgs...)

# ...the imaginary part
heatmap(tfAxes(Y)..., imag(Y.y);
        c=:bluesreds, tfArgs...)

# ...the amplitude
heatmap(tfAxes(Y)..., amplitude(Y.y);
        c=:amp, tfArgs...)

# ...the amplitude of the AS smoothed in frequency and time
heatmap(tfAxes(Y)..., amplitude(smooth(hannSmoother, hannSmoother, Y).y);
        c=:fire, tfArgs...)

# ...the phase
heatmap(tfAxes(Y)..., phase(Y.y);
        c=:bluesreds, tfArgs...)

# or generate a TFAmplitude object from the AS
A=TFamplitude(Y)
# and plot it (with different colors)
heatmap(tfAxes(A)..., A.y;
        c=:fire, tfArgs...)

# generate a TFPhase object
ϴ=TFphase(Y)
# and plot it (with custom colors)
heatmap(tfAxes(ϴ)..., ϴ.y;
        c=:pu_or, tfArgs...)

# compute and plot phase in [0, 2π]
heatmap(tfAxes(Y)..., TFphase(Y; func=x->x+π).y;
        c=:amp, tfArgs...)

# compute and plot unwrapped phase
heatmap(tfAxes(Y)..., TFphase(Y; unwrapped=true).y;
        c=:amp, tfArgs...)

# smooth time-frequency AS: smooth frequency
Z=smooth(blackmanSmoother, noSmoother, Y)

# plot amplitude of smoothed analytic signal
heatmap(tfAxes(Z)..., amplitude(Z.y);
        c=:amp, tfArgs...)

# not equivalently (!), create an amplitude object and smooth it:
# in this case the amplitude is smoothed, not the AS
A=smooth(blackmanSmoother, noSmoother, TFamplitude(Y))
heatmap(tfAxes(A)..., A.y;
        c=:fire, tfArgs...)

# Smoothing raw phase estimates is unappropriate
# since the phase is a discontinous function, however it makes sense
# to smooth phase if the phase is unwrapped.
heatmap(tfAxes(Y)...,
        smooth(blackmanSmoother, noSmoother, TFphase(Y; unwrapped=true)).y;
        c=:amp, tfArgs...)

# smooth AS: smooth both frequency and time
E=smooth(blackmanSmoother, blackmanSmoother, Y)

# plot amplitude of smoothed analytic signal
heatmap(tfAxes(E)..., amplitude(E.y);
        c=:fire, tfArgs...)

# plot phase of smoothed analytic signal
heatmap(tfAxes(E)..., phase(E.y);
        c=:bluesreds, tfArgs...)

# not equivalently (!), create amplitude and phase objects and smooth them
A=smooth(blackmanSmoother, blackmanSmoother, TFamplitude(Y))
heatmap(tfAxes(A)..., A.y;
        c=:fire, tfArgs...)

ϴ=smooth(blackmanSmoother, blackmanSmoother, TFphase(Y, unwrapped=true))
heatmap(tfAxes(ϴ)..., ϴ.y;
        c=:pu_or, tfArgs...)

# smooth again
ϴ=smooth(blackmanSmoother, blackmanSmoother, ϴ)
heatmap(tfAxes(ϴ)..., ϴ.y;
        c=:pu_or, tfArgs...)
# and again ...

# create directly smoothed AS
Y=TFanalyticsignal(x, sr, t, bandwidth;
                   fmax=32,
                   fsmoothing=hannSmoother,
                   tsmoothing=hannSmoother)

# plot amplitude of smoothed analytic signal
heatmap(tfAxes(Y)..., amplitude(Y.y);
        c=:amp, tfArgs...)

# create directly smoothed Amplitude
A=TFamplitude(x, sr, t, bandwidth;
              fmax=32,
              fsmoothing=hannSmoother,
              tsmoothing=hannSmoother)

# plot smoothed amplitude
heatmap(tfAxes(A)..., A.y;
        c=:amp, tfArgs...)

# compute a TFAnalyticSignal object with non-linear AS
Y=TFanalyticsignal(x, sr, t, bandwidth; fmax=32, nonlinear=true)

# check that it is non-linear
Y.nonlinear

# check that the amplitude is now 1.0 everywhere
norm(amplitude(Y.y)-ones(eltype(Y.y), size(Y.y))) # must be zero

# plot non-linear phase
heatmap(tfAxes(Y)..., phase(Y.y);
        c=:bkr, tfArgs...)

# get the center frequencies of TFAmplitude object A
A.flabels

# extract the amplitude in a time-frequency region
extract(A, (2, 10), (1, 256)) # output a real matrix

# extract the amplitude in a time-frequency region at only one frequency
extract(A, 10, (1, 256)) # output a row vector

# extract the amplitude at one temporal sample at one frequency
extract(A, 10, 12) # or extract(A, 10.0, 12)

# extract amplitude at one temporal sample in a frequency region
extract(A, (10, 12), 12) # or extract(A, (10.0, 12.0), 12)

# extract amplitude at one temporal sample and all frequencies
extract(A, :, 12) # output a (column) vector

# compute the mean in a time-frequency region:
mean(A, (2, 10), (1, 256)) # output a real number
# is equivalent to (but may be less efficient than)
mean(extract(A, (2, 10), (1, 256)))

# using column sign for extracting all time samples
extract(A, (2, 10), :)

# This:
extract(A, :, :)
# is equivalent to this:
A.y
# but if you don't need to extract all frequencies,
# use the extract function to control what frequencies will be extracted:
# This
extract(A, (4, 5), 10)
# is not equivalent to this
A.y[4:5, 10]
# since the `extract` function finds the correct rows corresponding
# to the sought frequencies (in Hz), while A.y[4:5, 10]
# just returns the elements [4:5, 10] in the TF amplitude object

# Although the first center frequency in A is 2Hz, its
# band-pass region is 1-3Hz, therefore the frequency range 1:10 is accepted
mean(A, (1, 10), (1, 256))
# but this result in an error (see the REPL) since 0.5 is out of range:
mean(A, (0.5, 10), (1, 256))

# using a colon sign for time range
a=mean(A, (1, 10), :)
# using an integer for time range (indicates one specific sample)
a=mean(A, (1, 10), 16)

# using a colon sign for frequency range
a=mean(A, :, (1, 16))
# using a real number for frequency range
a=mean(A, 8.5, (1, 16))

# NB: the `extract` and `mean` functions work with the same syntax
#     for objects TFAnayticSignal, TFAmplitude and TFPhase.
```
"""
function TFanalyticsignal(x         :: Vector{T},
                          sr        :: Int,
                          wl        :: Int          = 0,
                          bandwidth :: IntOrReal    = 2;
                      fmin          :: IntOrReal    = bandwidth,
                      fmax          :: IntOrReal    = sr÷2,
                      filtkind      :: FilterDesign = Butterworth(2),
                      nonlinear     :: Bool         = false,
                      fsmoothing    :: Smoother     = noSmoother,
                      tsmoothing    :: Smoother     = noSmoother,
                      planner       :: Planner      = getplanner,
                      ⏩           :: Bool         = true) where T<:Real
    # NB  if `t` > 2^14 then `t` is set to 2^10 by `anaytic signal function`.
    f, F=filterbank(x, sr, bandwidth;
                filtkind=filtkind, fmin=fmin, fmax=fmax, ⏩=⏩)

    smooth(fsmoothing,
           tsmoothing,
           TFAnalyticSignal(Matrix(transpose(analyticsignal(F,
                                                            wl==0 ? length(x) : wl;
                                                        nonlinear=nonlinear,
                                                        planner=planner,
                                                        ⏩=⏩))),
                            bandwidth,
                            f,
                            nonlinear,
                            noSmoother,
                            noSmoother),
            "TFanalyticsignal")
end

function TFanalyticsignal(𝐱         :: Vector{Vector{T}},
                          sr        :: Int,
                          wl        :: Int,
                          bandwidth :: IntOrReal    = 2;
                      filtkind      :: FilterDesign = Butterworth(2),
                      fmin          :: IntOrReal    = bandwidth,
                      fmax          :: IntOrReal    = sr÷2,
                      nonlinear     :: Bool         = false,
                      fsmoothing    :: Smoother     = noSmoother,
                      tsmoothing    :: Smoother     = noSmoother,
                      planner       :: Planner      = getplanner,
                      ⏩           :: Bool         = true) where T<:Real
    # NB  if `t` > 2^14 then `t` is set to 2^10 by `anaytic signal function`.
    # differently from previous method, wl here must be explicitly provided
    k = length(𝐱)
    𝐙 = TFAnalyticSignalVector(undef, k)

    planner ≠ getplanner ? plan=planner : plan=Planner(plan_estimate, -1.0, wl, T, true)

    function TF!(i::Int, threaded::Bool)
        f, F = filterbank(𝐱[i],
                          sr,
                          bandwidth;
                      filtkind=filtkind,
                      fmin=fmin,
                      fmax=fmax,
                      ⏩=threaded)

        𝐙[i] = TFAnalyticSignal(Matrix(transpose(analyticsignal(F,
                                                                wl;
                                                            nonlinear=nonlinear,
                                                            planner=plan,
                                                            ⏩=threaded))),
                                bandwidth,
                                f,
                                nonlinear,
                                noSmoother,
                                noSmoother)
    end

    _thread(⏩, k) ? (@threads for i=1:k TF!(i, false) end) : (for i=1:k TF!(i, true) end)
    return smooth(fsmoothing, tsmoothing, 𝐙, "TFanalyticsignal")
end

# ++++++++++++++++++++  Show override  +++++++++++++++++++ # (REPL output)
function Base.show(io::IO, ::MIME{Symbol("text/plain")}, Y::TFAnalyticSignal)
    r=size(Y.y, 1)
    c=size(Y.y, 2)
    l=length(Y.flabels)
    println(io, titleFont, "▤ TFAnayticSignal type; $r freq. x $c samples")
    #println(io, "□  □    □      □        □           □", defaultFont)
    println(io, separatorFont, "⭒  ⭒    ⭒      ⭒        ⭒           ⭒", defaultFont)
    println(io, "non-linear(.nonlinear): $(Y.nonlinear)")
    println(io, "freq. sm.(.fsmoothing): ", string(Y.fsmoothing))
    println(io, "time  sm.(.tsmoothing): ", string(Y.tsmoothing))
    println(io, "bandwidth             : $(Y.bandwidth) Hz")
    println(io, "freq. lab.  (.flabels): $(l)-", typeof(Y.flabels))
    println(io, "data              (.y): $(r)x$(c)-", typeof(Y.y))
    r≠l && @warn "number of frequency labels does not match the data matrix size" l r
end

function Base.show(io::IO, ::MIME{Symbol("text/plain")}, 𝐘::TFAnalyticSignalVector)
    println(io, titleFont, "▤ ⋯ ▤ TFAnalyticSignalVector Type")
    println(io, separatorFont, "⭒  ⭒    ⭒      ⭒        ⭒           ⭒", defaultFont)
    println(io, "$(length(𝐘))-element Vector{TFAnalyticSignal}")
end
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

#### Time-Frequency Amplitude ####

#############################################################################
# Generic constructors of TFAmplitude objects:
##############################################
# Enable construction giving only `y`, `bandwidth`, `flabels`, `fsmoothing`
# and `tsmoothing`. `func` is set to `identity`.
TFAmplitude(y, bandwidth, flabels, fsmoothing, tsmoothing) =
    TFAmplitude(y, bandwidth, flabels, fsmoothing, tsmoothing, identity)
#
# As above, but setting by default also both
# `fsmoothing` and `tsmoothing` to `noSmoother`.
TFAmplitude(y, bandwidth, flabels) =
    TFAmplitude(y, bandwidth, flabels, noSmoother, noSmoother, identity)

#############################################################################
# Constructors from Data
##############################################
"""
```
(1)
function TFamplitude(Z::TFAnalyticSignal;
                func::Function=identity)

(2)
function TFamplitude(𝐙::TFAnalyticSignalVector;
                func::Function=identity)

(3)
function TFamplitude(x         :: Vector{T},
                     sr        :: Int,
                     wl        :: Int,
                     bandwidth :: IntOrReal = 2;
                func       :: Function     = identity,
                filtkind   :: FilterDesign = Butterworth(2),
                fmin       :: IntOrReal    = bandwidth,
                fmax       :: IntOrReal    = sr÷2,
                fsmoothing :: Smoother     = noSmoother,
                tsmoothing :: Smoother     = noSmoother,
                planner    :: Planner      = getplanner,
                ⏩        :: Bool         = true) where T<:Real =

(4)
function TFamplitude(𝐱 :: Vector{Vector{T}},
                < same arguments as method (3) >
```

(1)

Construct a [TFAmplitude](@ref) object computing the amplitude of
[TFAnalyticSignal](@ref) object `Z`.
Optional keyword argument `func` is a function to be applied
element-wise to the data matrix of the output. By default,
the `identity` (do nothing) function is applied.

(2)

Construct a [TFAmplitudeVector](@ref) object from a
[TFAnalyticSignalVector](@ref) object executing method (1)
for all [TFAnalyticSignal](@ref) objects in `𝐙`

(3)

Call [`TFanalyticsignal`](@ref) to obtain the time-frequency
analytic signal of real signal vector `x` and construct a [TFAmplitude](@ref)
object holding the time-frequency amplitude (the mudulus, often
referred to as the *envelope*) of `x`.

All arguments are used for regulating the estimation of the analytic signal,
with the exception of `func`, `fsmoothing` and `fsmoothing`:

`func` is an optional function to be applied to the amplitude data matrix output.

In order to estimate the analytic signal in the time-frequency domain
this function calls the [`TFanalyticsignal`](@ref) constructor
(method (1) therein), with both `fsmoothing` and `tsmoothing` arguments
set to `noSmoother`.
Arguments `fsmoothing` and `fsmoothing` are then used to smooth the amplitude.

In order to obtain amplitude estimations on smoothed analytic signal instead,
create a [TFAnalyticSignal](@ref) object passing a
[Smoother](@ref) to the [`TFanalyticsignal`](@ref)
constructor and then use method (1) to obtain the amplitude.
Such amplitude estimation can be further smoothed
using the [`smooth`](@ref) function, as shown in the examples.

For the meaning of all other arguments, which are passed to function
[`TFanalyticsignal`](@ref), see the documentation therein.

(4)

Construct a [TFAmplitudeVector](@ref) object from a
vector of real signal vectors `𝐱`, executing method (3) for all of them.
In order to estimate the time-frequency analytic signal for a
vector of signals, method (2) of [`TFanalyticsignal`](@ref) is called.

**See**: [`TFanalyticsignal`](@ref), [TFAmplitude](@ref).

**See also**: [`amplitude`](@ref).

**Examples**: see the examples of [`TFanalyticsignal`](@ref).
"""
TFamplitude(Z::TFAnalyticSignal; func::Function=identity) =
    TFAmplitude(amplitude(Z.y; func=func), Z.bandwidth, Z.flabels, Z.fsmoothing, Z.tsmoothing, func)

TFamplitude(𝐙::TFAnalyticSignalVector; func::Function=identity) =
    TFAmplitudeVector([TFAmplitude(amplitude(Z.y; func=func), Z.bandwidth, Z.flabels, Z.fsmoothing, Z.tsmoothing, func) for Z ∈ 𝐙])

TFamplitude(𝐱         :: Union{Vector{T}, Vector{Vector{T}}},
            sr        :: Int,
            wl        :: Int,
            bandwidth :: IntOrReal    = 2;
        func          :: Function     = identity,
        filtkind      :: FilterDesign = Butterworth(2),
        fmin          :: IntOrReal    = bandwidth,
        fmax          :: IntOrReal    = sr÷2,
        fsmoothing    :: Smoother     = noSmoother,
        tsmoothing    :: Smoother     = noSmoother,
        planner       :: Planner      = getplanner,
        ⏩           :: Bool         = true) where T<:Real =

    smooth(fsmoothing,
           tsmoothing,
           TFamplitude(TFanalyticsignal(𝐱,
                                        sr,
                                        wl,
                                        bandwidth;
                                    filtkind  = filtkind,
                                    fmin      = fmin,
                                    fmax      = fmax,
                                    planner   = planner,
                                    ⏩       = ⏩);
                        func=func),
           "TFamplitude")

# ++++++++++++++++++++  Show override  +++++++++++++++++++ # (REPL output)
function Base.show(io::IO, ::MIME{Symbol("text/plain")}, Y::TFAmplitude)
   r=size(Y.y, 1)
   c=size(Y.y, 2)
   l=length(Y.flabels)
   println(io, titleFont, "▤ TFAmplitude type; $r freq. x $c samples")
   #println(io, "□  □    □      □        □           □", defaultFont)
   println(io, separatorFont, "⭒  ⭒    ⭒      ⭒        ⭒           ⭒", defaultFont)
   println(io, "freq. sm.(.fsmoothing): ", string(Y.fsmoothing))
   println(io, "time  sm.(.tsmoothing): ", string(Y.tsmoothing))
   println(io, "bandwidth             : $(Y.bandwidth) Hz")
   println(io, "function     (.func)  : ", string(Y.func))
   println(io, "freq. lab.(.flabels)  : $(l)-", typeof(Y.flabels))
   println(io, "data            (.y)  : $(r)x$(c)-", typeof(Y.y))
   r≠l && @warn "number of frequency labels does not match the data matrix size" l r
end

function Base.show(io::IO, ::MIME{Symbol("text/plain")}, 𝐘::TFAmplitudeVector)
   println(io, titleFont, "▤ ⋯ ▤ TFAmplitudeVector Type")
   println(io, separatorFont, "⭒  ⭒    ⭒      ⭒        ⭒           ⭒", defaultFont)
   println(io, "$(length(𝐘))-element Vector{TFAmplitude}")
end
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

#### Time-Frequency Phase ####

#############################################################################
# Generic constructors of TFAmplitude objects:
##############################################
#
# Enable construction giving only `y`, `bandwidth`, `flabels`, `fsmoothing`
# and `tsmoothing`. `unwrapped` is set to false and `func` is set to `identity`.
TFPhase(y, bandwidth, flabels, nonlinear, fsmoothing, tsmoothing) =
    TFPhase(y, bandwidth, flabels, nonlinear, fsmoothing, tsmoothing, false, identity)
#
# As above, but setting by default also `nonlinear` to true and both
# `fsmoothing` and `tsmoothing` to `noSmoother`.
TFPhase(y, bandwidth, flabels) =
    TFPhase(y, bandwidth, flabels, true, noSmoother, noSmoother, false, identity)

#############################################################################
# Constructors from Data
##############################################
"""
```
(1)
function TFphase(Z :: TFAnalyticSignal;
            func      :: Function = identity,
            unwrapped :: Bool     = false)

(2)
function TFphase(𝐙 :: TFAnalyticSignalVector;
            func      ::Function = identity,
            unwrapped ::Bool     = false)

(3)
function TFphase(x         :: Vector{T},
                 sr        :: Int,
                 wl        :: Int,
                 bandwidth :: IntOrReal = 2;
           unwrapped  :: Bool         = false,
           func       :: Function     = identity,
           filtkind   :: FilterDesign = Butterworth(2),
           fmin       :: IntOrReal    = bandwidth,
           fmax       :: IntOrReal    = sr÷2,
           nonlinear  :: Bool         = false,
           fsmoothing :: Smoother     = noSmoother,
           tsmoothing :: Smoother     = noSmoother,
           planner    :: Planner      = getplanner,
           ⏩        :: Bool         = true) where T<:Real

(4)
function TFphase(𝐱 :: Vector{Vector{T}},
                < same arguments as method (3) >
```

(1)

Construct a [TFPhase](@ref) object computing the phase of
[TFAnalyticSignal](@ref) object `Z`. By default the phase is
represented in ``[−π, π]``.

If optional keyword argument `unwrapped` is true (false by defalut),
the phase is unwrapped, that is, it holds the cumulative sum of the phase
along the time dimension once this is represented in ``[0, 2π]``.

Optional keyword argument `func` is a function to be applied
element-wise to the data matrix of the output. By default,
the `identity` (do nothing) function is applied. If `unwrapped`
is true, the function is applied on the unwrapped phase.

(2)

Construct a [TFPhaseVector](@ref) object from a
[TFAnalyticSignalVector](@ref) object executing method (1)
for all [TFAnalyticSignal](@ref) objects in `𝐙`

(3)

Call [`TFanalyticsignal`](@ref) to obtain the time-frequency
analytic signal of real signal vector `x` and construct a [TFPhase](@ref)
object holding the time-frequency phase (argument) of `x`.

All arguments are used for regulating the estimation of the analytic signal,
with the exception of `unwrapped`, `func`, `fsmoothing` and `fsmoothing`.

`unwrapped` has the same meaning as in method (1) and (2).

`func` is an optional function to be applied to the phase data matrix output.
If `unwrapped` is true, it is applied to the unwrapped phase.

In order to estimate the analytic signal in the time-frequency domain
this function calls the [`TFanalyticsignal`](@ref) constructor
(method (1) therein),
with both `fsmoothing` and `tsmoothing` arguments set to `noSmoother`.
`fsmoothing` and `fsmoothing` are then used to smooth the phase if `unwrapped`
is true.

In order to obtain phase estimations on smoothed analytic signal instead,
create a [TFAnalyticSignal](@ref) object passing a
[Smoother](@ref) to the [`TFanalyticsignal`](@ref)
constructor and then use method (1) to obtain the phase.

For the meaning of all other arguments, which are passed to function
[`TFanalyticsignal`](@ref), see the documentation therein.

(4)

Construct a [TFPhaseVector](@ref) object from a
vector of real signal vectors `𝐱`, executing method (3) for all of them.
In order to estimate the time-frequency analytic signal for a
vector of signals, method (2) of [`TFanalyticsignal`](@ref) is called.

**See**: [`TFanalyticsignal`](@ref), [TFPhase](@ref).

**See also**: [`phase`](@ref), [`unwrapPhase`](@ref), [`polar`](@ref).

**Examples**: see the examples of [`TFanalyticsignal`](@ref).
"""
TFphase(Z::TFAnalyticSignal; func::Function=identity, unwrapped::Bool=false) =
    TFPhase(phase(Z.y; unwrapdims=2*unwrapped, func=func), Z.bandwidth, Z.flabels, Z.nonlinear, Z.fsmoothing, Z.tsmoothing, unwrapped, func)

TFphase(𝐙::TFAnalyticSignalVector; func::Function=identity, unwrapped::Bool=false) =
    TFPhaseVector([TFphase(Z; func=func, unwrapped=unwrapped) for Z in 𝐙])

function TFphase(𝐱         :: Union{Vector{T}, Vector{Vector{T}}},
                 sr        :: Int,
                 wl        :: Int,
                 bandwidth :: IntOrReal    = 2;
           unwrapped  :: Bool         = false,
           func       :: Function     = identity,
           filtkind   :: FilterDesign = Butterworth(2),
           fmin       :: IntOrReal    = bandwidth,
           fmax       :: IntOrReal    = sr÷2,
           nonlinear  :: Bool         = false,
           fsmoothing :: Smoother     = noSmoother,
           tsmoothing :: Smoother     = noSmoother,
           planner    :: Planner      = getplanner,
           ⏩        :: Bool         = true) where T<:Real

   TFP = TFphase(TFanalyticsignal(𝐱,
                                 sr,
                                 wl,
                                 bandwidth;
                                 filtkind  = filtkind,
                                 fmin      = fmin,
                                 fmax      = fmax,
                                 nonlinear = nonlinear,
                                 planner   = planner,
                                 ⏩       = ⏩);
                         unwrapped = unwrapped,
                         func      = func)

   !unwrapped && (fsmoothing≠noSmoother || tsmoothing≠noSmoother) &&
   @warn 📌*", call to TFPhase constructor; smoothing is meaningful only on unwrapped phase. Smoothing has been disabled.)" unwrapped

   return unwrapped ? smooth(fsmoothing,
                             tsmoothing,
                             TFP,
                             "TFphase") :
                      TFP
end

# ++++++++++++++++++++  Show override  +++++++++++++++++++ # (REPL output)
function Base.show(io::IO, ::MIME{Symbol("text/plain")}, Y::TFPhase)
    r=size(Y.y, 1)
    c=size(Y.y, 2)
    l=length(Y.flabels)
    println(io, titleFont, "▤ TFPhase type; $r freq. x $c samples")
    #println(io, "□  □    □      □        □           □", defaultFont)
    println(io, separatorFont, "⭒  ⭒    ⭒      ⭒        ⭒           ⭒", defaultFont)
    println(io, "non-linear(.nonlinear): $(Y.nonlinear)")
    println(io, "freq. sm.(.fsmoothing): ", string(Y.fsmoothing))
    println(io, "time  sm.(.tsmoothing): ", string(Y.tsmoothing))
    println(io, "bandwidth             : $(Y.bandwidth) Hz")
    println(io, "unwrapped             : $(Y.unwrapped)")
    println(io, "function     (.func)  : ", string(Y.func))
    println(io, "freq. lab.(.flabels)  : $(l)-", typeof(Y.flabels))
    println(io, "data            (.y)  : $(r)x$(c)-", typeof(Y.y))
    r≠l && @warn "number of frequency labels does not match the data matrix size" l r
end

function Base.show(io::IO, ::MIME{Symbol("text/plain")}, 𝐘::TFPhaseVector)
    println(io, titleFont, "▤ ⋯ ▤ TFPhaseVector Type")
    println(io, separatorFont, "⭒  ⭒    ⭒      ⭒        ⭒           ⭒", defaultFont)
    println(io, "$(length(𝐘))-element Vector{TFPhase}")
end
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #


text:
# Abstract QLearning Policy
abstract type QPolicy <: Reinforce.AbstractPolicy end
save_policy(policy::QPolicy, filename="model_checkpoint.bson") = BSON.bson(filename, Dict(:policy => policy))
load_policy(filename="model_checkpoint.bson") = BSON.load(filename)[:policy]
Reinforce.reset!(policy::QPolicy) = nothing

# Individual Policies are in the policies folder for easier organization
include("policies/EpGreedyPolicy.jl")
include("policies/DeepQPolicy.jl")
include("policies/DoubleDeepQPolicy.jl")
include("policies/DuelingDouble_DeepQPolicy.jl")


text:
using MathOptInterface
const MOI = MathOptInterface
const CI = MOI.ConstraintIndex
const VI = MOI.VariableIndex

const MOIU = MOI.Utilities

"""
    ADMMIterations()

The number of ADMM iterations completed during the solve.
"""
struct ADMMIterations <: MOI.AbstractModelAttribute
end

MOI.is_set_by_optimize(::ADMMIterations) = true

mutable struct MOISolution
    ret_val::Int
    raw_status::String
    primal::Vector{Float64}
    dual::Vector{Float64}
    slack::Vector{Float64}
    objective_value::Float64
    dual_objective_value::Float64
    objective_constant::Float64
    solve_time_sec::Float64
    iterations::Int
end
MOISolution() = MOISolution(0, # SCS_UNFINISHED
                            "", Float64[], Float64[], Float64[], NaN, NaN, NaN,
                            0.0, 0)

# Used to build the data with allocate-load during `copy_to`.
# When `optimize!` is called, a the data is passed to SCS
# using `SCS_solve` and the `ModelData` struct is discarded
mutable struct ModelData
    m::Int # Number of rows/constraints
    n::Int # Number of cols/variables
    I::Vector{Int} # List of rows
    J::Vector{Int} # List of cols
    V::Vector{Float64} # List of coefficients
    b::Vector{Float64} # constants
    objective_constant::Float64 # The objective is min c'x + objective_constant
    c::Vector{Float64}
end

# This is tied to SCS's internal representation
mutable struct ConeData
    f::Int # number of linear equality constraints
    l::Int # length of LP cone
    q::Int # length of SOC cone
    qa::Vector{Int} # array of second-order cone constraints
    s::Int # length of SD cone
    sa::Vector{Int} # array of semi-definite constraints
    ep::Int # number of primal exponential cone triples
    ed::Int # number of dual exponential cone triples
    p::Vector{Float64} # array of power cone params
    nrows::Dict{Int, Int} # The number of rows of each vector sets, this is used by `constrrows` to recover the number of rows used by a constraint when getting `ConstraintPrimal` or `ConstraintDual`
    function ConeData()
        new(0, 0, 0, Int[], 0, Int[], 0, 0, Float64[], Dict{Int, Int}())
    end
end

mutable struct Optimizer <: MOI.AbstractOptimizer
    cone::ConeData
    maxsense::Bool
    data::Union{Nothing, ModelData} # only non-Void between MOI.copy_to and MOI.optimize!
    sol::MOISolution
    silent::Bool
    options::Dict{Symbol, Any}
    function Optimizer(; kwargs...)
        optimizer = new(ConeData(), false, nothing, MOISolution(), false,
                        Dict{Symbol, Any}())
        for (key, value) in kwargs
            MOI.set(optimizer, MOI.RawParameter(key), value)
        end
        return optimizer
    end

end

MOI.get(::Optimizer, ::MOI.SolverName) = "SCS"

function MOI.set(optimizer::Optimizer, param::MOI.RawParameter, value)
    optimizer.options[param.name] = value
end
function MOI.get(optimizer::Optimizer, param::MOI.RawParameter)
    return optimizer.options[param.name]
end

MOI.supports(::Optimizer, ::MOI.Silent) = true
function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool)
    optimizer.silent = value
end
MOI.get(optimizer::Optimizer, ::MOI.Silent) = optimizer.silent

function MOI.is_empty(optimizer::Optimizer)
    !optimizer.maxsense && optimizer.data === nothing
end
function MOI.empty!(optimizer::Optimizer)
    optimizer.maxsense = false
    optimizer.data = nothing # It should already be nothing except if an error is thrown inside copy_to
    optimizer.sol.ret_val = 0
end

MOIU.supports_allocate_load(::Optimizer, copy_names::Bool) = !copy_names

function MOI.supports(::Optimizer,
                      ::Union{MOI.ObjectiveSense,
                              MOI.ObjectiveFunction{MOI.SingleVariable},
                              MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}})
    return true
end

function MOI.supports(::Optimizer, ::MOI.VariablePrimalStart,
                      ::Type{MOI.VariableIndex})
    return true
end

function MOI.supports(::Optimizer,
                      ::Union{MOI.ConstraintPrimalStart,
                              MOI.ConstraintDualStart},
                      ::Type{<:MOI.ConstraintIndex})
    return true
end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{<:MOI.VectorAffineFunction{Float64}},
    ::Type{<:Union{MOI.Zeros, MOI.Nonnegatives, MOI.SecondOrderCone,
                   MOI.ExponentialCone, MOI.DualExponentialCone,
                   MOI.PositiveSemidefiniteConeTriangle,
                   MOI.PowerCone, MOI.DualPowerCone}})
    return true
end

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; kws...)
    return MOIU.automatic_copy_to(dest, src; kws...)
end

using SparseArrays

# Computes cone dimensions
function constroffset(cone::ConeData,
                      ci::CI{<:MOI.AbstractFunction, MOI.Zeros})
    return ci.value
end
#_allocate_constraint: Allocate indices for the constraint `f`-in-`s`
# using information in `cone` and then update `cone`
function _allocate_constraint(cone::ConeData, f, s::MOI.Zeros)
    ci = cone.f
    cone.f += MOI.dimension(s)
    return ci
end
function constroffset(cone::ConeData,
                      ci::CI{<:MOI.AbstractFunction, MOI.Nonnegatives})
    return cone.f + ci.value
end
function _allocate_constraint(cone::ConeData, f, s::MOI.Nonnegatives)
    ci = cone.l
    cone.l += MOI.dimension(s)
    return ci
end
function constroffset(cone::ConeData,
                      ci::CI{<:MOI.AbstractFunction, MOI.SecondOrderCone})
    return cone.f + cone.l + ci.value
end
function _allocate_constraint(cone::ConeData, f, s::MOI.SecondOrderCone)
    push!(cone.qa, s.dimension)
    ci = cone.q
    cone.q += MOI.dimension(s)
    return ci
end
function constroffset(cone::ConeData,
                      ci::CI{<:MOI.AbstractFunction,
                             MOI.PositiveSemidefiniteConeTriangle})
    return cone.f + cone.l + cone.q + ci.value
end
function _allocate_constraint(cone::ConeData, f,
                              s::MOI.PositiveSemidefiniteConeTriangle)
    push!(cone.sa, s.side_dimension)
    ci = cone.s
    cone.s += MOI.dimension(s)
    return ci
end
function constroffset(cone::ConeData,
                      ci::CI{<:MOI.AbstractFunction, MOI.ExponentialCone})
    return cone.f + cone.l + cone.q + cone.s + ci.value
end
function _allocate_constraint(cone::ConeData, f, s::MOI.ExponentialCone)
    ci = 3cone.ep
    cone.ep += 1
    return ci
end
function constroffset(cone::ConeData,
                      ci::CI{<:MOI.AbstractFunction, MOI.DualExponentialCone})
    return cone.f + cone.l + cone.q + cone.s + 3cone.ep + ci.value
end
function _allocate_constraint(cone::ConeData, f, s::MOI.DualExponentialCone)
    ci = 3cone.ed
    cone.ed += 1
    return ci
end
function constroffset(cone::ConeData,
                      ci::CI{<:MOI.AbstractFunction, <:MOI.PowerCone})
    return cone.f + cone.l + cone.q + cone.s + 3cone.ep + 3cone.ed + ci.value
end
function _allocate_constraint(cone::ConeData, f, s::MOI.PowerCone)
    ci = length(cone.p)
    push!(cone.p, s.exponent)
    return ci
end
function constroffset(cone::ConeData,
                      ci::CI{<:MOI.AbstractFunction, <:MOI.DualPowerCone})
    return cone.f + cone.l + cone.q + cone.s + 3cone.ep + 3cone.ed + ci.value
end
function _allocate_constraint(cone::ConeData, f, s::MOI.DualPowerCone)
    ci = length(cone.p)
    # SCS' convention: dual cones have a negative exponent.
    push!(cone.p, -s.exponent)
    return ci
end
function constroffset(optimizer::Optimizer, ci::CI)
    return constroffset(optimizer.cone, ci::CI)
end
function MOIU.allocate_constraint(optimizer::Optimizer, f::F, s::S) where {F <: MOI.AbstractFunction, S <: MOI.AbstractSet}
    return CI{F, S}(_allocate_constraint(optimizer.cone, f, s))
end

# Vectorized length for matrix dimension n
sympackedlen(n) = div(n*(n+1), 2)
# Matrix dimension for vectorized length n
sympackeddim(n) = div(isqrt(1+8n) - 1, 2)
trimap(i::Integer, j::Integer) = i < j ? trimap(j, i) : div((i-1)*i, 2) + j
trimapL(i::Integer, j::Integer, n::Integer) = i < j ? trimapL(j, i, n) : i + div((2n-j) * (j-1), 2)
function _sympackedto(x, n, mapfrom, mapto)
    @assert length(x) == sympackedlen(n)
    y = similar(x)
    for i in 1:n, j in 1:i
        y[mapto(i, j)] = x[mapfrom(i, j)]
    end
    y
end
sympackedLtoU(x, n=sympackeddim(length(x))) = _sympackedto(x, n, (i, j) -> trimapL(i, j, n), trimap)
sympackedUtoL(x, n=sympackeddim(length(x))) = _sympackedto(x, n, trimap, (i, j) -> trimapL(i, j, n))

function sympackedUtoLidx(x::AbstractVector{<:Integer}, n)
    y = similar(x)
    map = sympackedLtoU(1:sympackedlen(n), n)
    for i in eachindex(y)
        y[i] = map[x[i]]
    end
    y
end

# Scale coefficients depending on rows index
# rows: List of row indices
# coef: List of corresponding coefficients
# d: dimension of set
# rev: if true, we unscale instead (e.g. divide by √2 instead of multiply for PSD cone)
function _scalecoef(rows, coef, d, rev)
    scaling = rev ? 1 / √2 : 1 * √2
    output = copy(coef)
    diagidx = BitSet()
    for i in 1:d
        push!(diagidx, trimap(i, i))
    end
    for i in 1:length(output)
        if !(rows[i] in diagidx)
            output[i] *= scaling
        end
    end
    return output
end
# Unscale the coefficients in `coef` with respective rows in `rows` for a set `s`
scalecoef(rows, coef, s) = _scalecoef(rows, coef, MOI.dimension(s), false)
# Unscale the coefficients in `coef` with respective rows in `rows` for a set of type `S` with dimension `d`
unscalecoef(rows, coef, d) = _scalecoef(rows, coef, sympackeddim(d), true)

output_index(t::MOI.VectorAffineTerm) = t.output_index
variable_index_value(t::MOI.ScalarAffineTerm) = t.variable_index.value
variable_index_value(t::MOI.VectorAffineTerm) = variable_index_value(t.scalar_term)
coefficient(t::MOI.ScalarAffineTerm) = t.coefficient
coefficient(t::MOI.VectorAffineTerm) = coefficient(t.scalar_term)
# constrrows: Recover the number of rows used by each constraint.
# When, the set is available, simply use MOI.dimension
constrrows(s::MOI.AbstractVectorSet) = 1:MOI.dimension(s)
# When only the index is available, use the `optimizer.ncone.nrows` field
constrrows(optimizer::Optimizer, ci::CI{<:MOI.AbstractVectorFunction, <:MOI.AbstractVectorSet}) = 1:optimizer.cone.nrows[constroffset(optimizer, ci)]

orderval(val, s) = val
function orderval(val, s::MOI.PositiveSemidefiniteConeTriangle)
    sympackedUtoL(val, s.side_dimension)
end
orderidx(idx, s) = idx
function orderidx(idx, s::MOI.PositiveSemidefiniteConeTriangle)
    sympackedUtoLidx(idx, s.side_dimension)
end
function MOIU.load_constraint(optimizer::Optimizer, ci::MOI.ConstraintIndex, f::MOI.VectorAffineFunction, s::MOI.AbstractVectorSet)
    A = sparse(output_index.(f.terms), variable_index_value.(f.terms), coefficient.(f.terms))
    # sparse combines duplicates with + but does not remove zeros created so we call dropzeros!
    dropzeros!(A)
    I, J, V = findnz(A)
    offset = constroffset(optimizer, ci)
    rows = constrrows(s)
    optimizer.cone.nrows[offset] = length(rows)
    i = offset .+ rows
    b = f.constants
    if s isa MOI.PositiveSemidefiniteConeTriangle
        # FIXME shouldn't scale before ?
        b = sympackedUtoL(b, s.side_dimension)
        b = scalecoef(rows, b, s)
        V = scalecoef(I, V, s)
        I = sympackedUtoLidx(I, s.side_dimension)
    end
    # The SCS format is b - Ax ∈ cone
    optimizer.data.b[i] = b
    append!(optimizer.data.I, offset .+ I)
    append!(optimizer.data.J, J)
    append!(optimizer.data.V, -V)
end

function MOIU.allocate_variables(optimizer::Optimizer, nvars::Integer)
    optimizer.cone = ConeData()
    VI.(1:nvars)
end

function MOIU.load_variables(optimizer::Optimizer, nvars::Integer)
    cone = optimizer.cone
    m = cone.f + cone.l + cone.q + cone.s + 3cone.ep + 3cone.ed + 3length(cone.p)
    I = Int[]
    J = Int[]
    V = Float64[]
    b = zeros(m)
    c = zeros(nvars)
    optimizer.data = ModelData(m, nvars, I, J, V, b, 0., c)
    # `optimizer.sol` contains the result of the previous optimization.
    # It is used as a warm start if its length is the same, e.g.
    # probably because no variable and/or constraint has been added.
    if length(optimizer.sol.primal) != nvars
        optimizer.sol.primal = zeros(nvars)
    end
    @assert length(optimizer.sol.dual) == length(optimizer.sol.slack)
    if length(optimizer.sol.dual) != m
        optimizer.sol.dual = zeros(m)
        optimizer.sol.slack = zeros(m)
    end
end

function MOIU.allocate(::Optimizer, ::MOI.VariablePrimalStart,
                       ::MOI.VariableIndex, ::Union{Nothing, Float64})
end
function MOIU.allocate(::Optimizer, ::MOI.ConstraintPrimalStart,
                       ::MOI.ConstraintIndex,
                       ::Union{Nothing, AbstractVector{Float64}})
end
function MOIU.allocate(::Optimizer, ::MOI.ConstraintDualStart,
                       ::MOI.ConstraintIndex,
                       ::Union{Nothing, AbstractVector{Float64}})
end
function MOIU.allocate(optimizer::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    optimizer.maxsense = sense == MOI.MAX_SENSE
end
function MOIU.allocate(::Optimizer, ::MOI.ObjectiveFunction,
                       ::MOI.Union{MOI.SingleVariable,
                                   MOI.ScalarAffineFunction{Float64}})
end

function MOIU.load(::Optimizer, ::MOI.VariablePrimalStart,
                   ::MOI.VariableIndex, ::Nothing)
end
function MOIU.load(optimizer::Optimizer, ::MOI.VariablePrimalStart,
                   vi::MOI.VariableIndex, value::Float64)
    optimizer.sol.primal[vi.value] = value
end
function MOIU.load(::Optimizer, ::MOI.ConstraintPrimalStart,
                   ::MOI.ConstraintIndex, ::Nothing)
end
function MOIU.load(optimizer::Optimizer, ::MOI.ConstraintPrimalStart,
                   ci::MOI.ConstraintIndex, value)
    offset = constroffset(optimizer, ci)
    rows = constrrows(optimizer, ci)
    optimizer.sol.slack[offset .+ rows] .= value
end
function MOIU.load(::Optimizer, ::MOI.ConstraintDualStart,
                   ::MOI.ConstraintIndex, ::Nothing)
end
function MOIU.load(optimizer::Optimizer, ::MOI.ConstraintDualStart,
                   ci::MOI.ConstraintIndex, value)
    offset = constroffset(optimizer, ci)
    rows = constrrows(optimizer, ci)
    optimizer.sol.dual[offset .+ rows] .= value
end
function MOIU.load(::Optimizer, ::MOI.ObjectiveSense, ::MOI.OptimizationSense)
end
function MOIU.load(optimizer::Optimizer, ::MOI.ObjectiveFunction,
                   f::MOI.SingleVariable)
    MOIU.load(optimizer,
              MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
              MOI.ScalarAffineFunction{Float64}(f))
end
function MOIU.load(optimizer::Optimizer, ::MOI.ObjectiveFunction,
                   f::MOI.ScalarAffineFunction)
    c0 = Vector(sparsevec(variable_index_value.(f.terms), coefficient.(f.terms),
                          optimizer.data.n))
    optimizer.data.objective_constant = f.constant
    optimizer.data.c = optimizer.maxsense ? -c0 : c0
    return nothing
end

function MOI.optimize!(optimizer::Optimizer)
    cone = optimizer.cone
    m = optimizer.data.m
    n = optimizer.data.n
    A = sparse(optimizer.data.I, optimizer.data.J, optimizer.data.V)
    b = optimizer.data.b
    objective_constant = optimizer.data.objective_constant
    c = optimizer.data.c
    optimizer.data = nothing # Allows GC to free optimizer.data before A is loaded to SCS

    options = optimizer.options
    if optimizer.silent
        options = copy(options)
        options[:verbose] = 0
    end

    linear_solver, options = sanitize_SCS_options(options)

    sol = SCS_solve(linear_solver, m, n, A, b, c,
                    cone.f, cone.l, cone.qa, cone.sa, cone.ep, cone.ed, cone.p,
                    optimizer.sol.primal, optimizer.sol.dual,
                    optimizer.sol.slack; options...)

    ret_val = sol.ret_val
    primal = sol.x
    dual = sol.y
    slack = sol.s
    objective_value = (optimizer.maxsense ? -1 : 1) * sol.info.pobj
    dual_objective_value = (optimizer.maxsense ? -1 : 1) * sol.info.dobj
    solve_time = (sol.info.setupTime + sol.info.solveTime) / 1000
    optimizer.sol = MOISolution(ret_val, raw_status(sol.info), primal, dual,
                                slack, objective_value, dual_objective_value,
                                objective_constant, solve_time, sol.info.iter)
end

function MOI.get(optimizer::Optimizer, ::MOI.SolveTime)
    return optimizer.sol.solve_time_sec
end
function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
    return optimizer.sol.raw_status
end

function MOI.get(optimizer::Optimizer, ::ADMMIterations)
    return optimizer.sol.iterations
end

# Implements getter for result value and statuses
# SCS returns one of the following integers:
# -7 SCS_INFEASIBLE_INACCURATE
# -6 SCS_UNBOUNDED_INACCURATE
# -5 SCS_SIGINT
# -4 SCS_FAILED
# -3 SCS_INDETERMINATE
# -2 SCS_INFEASIBLE  : primal infeasible, dual unbounded
# -1 SCS_UNBOUNDED   : primal unbounded, dual infeasible
#  0 SCS_UNFINISHED  : never returned, used as placeholder
#  1 SCS_SOLVED
#  2 SCS_SOLVED_INACCURATE
function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus)
    s = optimizer.sol.ret_val
    @assert -7 <= s <= 2
    if s == -7
        return MOI.ALMOST_INFEASIBLE
    elseif s == -6
        return MOI.ALMOST_DUAL_INFEASIBLE
    elseif s == 2
        return MOI.ALMOST_OPTIMAL
    elseif s == -5
        return MOI.INTERRUPTED
    elseif s == -4
        return MOI.NUMERICAL_ERROR
    elseif s == -3
        return MOI.SLOW_PROGRESS
    elseif s == -2
        return MOI.INFEASIBLE
    elseif s == -1
        return MOI.DUAL_INFEASIBLE
    elseif s == 1
        return MOI.OPTIMAL
    else
        @assert s == 0
        return MOI.OPTIMIZE_NOT_CALLED
    end
end

function MOI.get(optimizer::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(optimizer, attr)
    value = optimizer.sol.objective_value
    if !MOIU.is_ray(MOI.get(optimizer, MOI.PrimalStatus()))
        value += optimizer.sol.objective_constant
    end
    return value
end
function MOI.get(optimizer::Optimizer, attr::MOI.DualObjectiveValue)
    MOI.check_result_index_bounds(optimizer, attr)
    value = optimizer.sol.dual_objective_value
    if !MOIU.is_ray(MOI.get(optimizer, MOI.DualStatus()))
        value += optimizer.sol.objective_constant
    end
    return value
end

function MOI.get(optimizer::Optimizer, attr::MOI.PrimalStatus)
    if attr.N > MOI.get(optimizer, MOI.ResultCount())
        return MOI.NO_SOLUTION
    end
    s = optimizer.sol.ret_val
    if s in (-3, 1, 2)
        MOI.FEASIBLE_POINT
    elseif s in (-6, -1)
        MOI.INFEASIBILITY_CERTIFICATE
    else
        MOI.INFEASIBLE_POINT
    end
end
function MOI.get(optimizer::Optimizer, attr::MOI.VariablePrimal, vi::VI)
    MOI.check_result_index_bounds(optimizer, attr)
    optimizer.sol.primal[vi.value]
end
function MOI.get(optimizer::Optimizer, attr::MOI.VariablePrimal, vi::Vector{VI})
    return MOI.get.(optimizer, attr, vi)
end
function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintPrimal,
                 ci::CI{<:MOI.AbstractFunction, S}) where S <: MOI.AbstractSet
    MOI.check_result_index_bounds(optimizer, attr)
    offset = constroffset(optimizer, ci)
    rows = constrrows(optimizer, ci)
    primal = optimizer.sol.slack[offset .+ rows]
    if S == MOI.PositiveSemidefiniteConeTriangle
        primal = sympackedLtoU(primal)
        primal = unscalecoef(rows, primal, length(rows))
    end
    return primal
end

function MOI.get(optimizer::Optimizer, attr::MOI.DualStatus)
    if attr.N > MOI.get(optimizer, MOI.ResultCount())
        return MOI.NO_SOLUTION
    end
    s = optimizer.sol.ret_val
    if s in (-3, 1, 2)
        MOI.FEASIBLE_POINT
    elseif s in (-7, -2)
        MOI.INFEASIBILITY_CERTIFICATE
    else
        MOI.INFEASIBLE_POINT
    end
end
function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintDual,
                 ci::CI{<:MOI.AbstractFunction, S}) where S <: MOI.AbstractSet
    MOI.check_result_index_bounds(optimizer, attr)
    offset = constroffset(optimizer, ci)
    rows = constrrows(optimizer, ci)
    dual = optimizer.sol.dual[offset .+ rows]
    if S == MOI.PositiveSemidefiniteConeTriangle
        dual = sympackedLtoU(dual)
        dual = unscalecoef(rows, dual, length(rows))
    end
    return dual
end

MOI.get(optimizer::Optimizer, ::MOI.ResultCount) = 1


text:
#!/usr/bin/env -S julia -iq
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 4., Zadanie 5.
#

include("interpolation.jl")
using .Interpolation

f = x -> exp(x)
a = 0.0
b = 1.0
n = 5

if length(ARGS) > 0
    ns = ARGS[1]
    n = parse(Int64, ns)
end

if length(ARGS) > 1
    fn = ARGS[2]
    if fn == "x"
        f = x -> x^2 * sin(x)
        a = -1.0
        b = 1.0
    end
end

Interpolation.rysujNnfx(f, a, b, n)

println("end the interactive session to close all plots")

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
"""

`MTH229`: helper functions for using `Julia` with MTH229

This module does two things:

* Install other useful packages with one command (`Roots`,  `ForwardDiff`, `QuadGK`, `SpecialFunctions`, ...) and re-exports their methods.

* Add a number of helper functions.

The helper functions include:

- `secant(f, a, b)`: return a function giving the secant line between ``(a,f(a))`` and ``(b,f(b))``.

- `tangent(f, c)`:  return a function giving the tangent line to ``f(x)`` at the point ``(c,f(c))``.

- `bisection(f, a, b)`: A simple implementation of the bisection
  method. The interval ``[a,b]`` should be a bracketing interval. This
  function makes an illustrative graphic. For real use of the bisection method, the `fzero(f,
  a, b)` function, from the `Roots` package, should be used.

- `'`: As in `f'`. Overloads `adjoint` allowing the derivative of   a function to be found as with math notation: `f'`.  The notation can be used for higher-order derivatives too: `f''`, `f'''`, ... This uses automatic differentiation from the `ForwardDiff` package.

- `plotif(f, g, a, b)`: Plot the function `f` over the interval `[a,b]` and color differently where ``g(x) > 0`` over ``[a,b]``. By passing in `f` for `g` shows where `f` is positive on `[a,b]`; passing in `f'` shows where `f` is increasing on `[a,b]`; and passing in `f''` shows where `f` is concave up on `[a,b]`.

- `sign_chart(f, a, b)`: shows a *signchart* of `f` by numerically identifying the zero crossings or infinities of `f` over `[a,b]` (assuming `f(a)` and `f(b)` are non zero, then checking the sign between these values. Calling `sign_chart(f', a, b)` is useful for the first-derivative test and `sign_chart(f'', a, b)` for the second-derivative test.

- `fisheye(f)` returns the composition `atan ∘ f ∘ tan`, which can be useful to find zeros of a function over the entire range of real numbers.

- `riemann(f, a, b, n; method="right")` An implementation of Riemann sums. The method can be "right" or "left" for Riemann sums, or "trapezoid" or "simpsons" for related approximations.

This package provides some plotting routines for `Plots`, `SimplePlots`, and `Makie`. For the latter two, some "recipes" for plotting functions and symbolic directions; and for all some convenience methods.
"""
module MTH229

# msg = """
# Loading the `MTH229` package.

# * Run the command `?MTH229` for a short description.

# """

# @info msg

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
    @eval Base.Experimental.@optlevel 1
end

using Reexport
@reexport using Roots
@reexport using SpecialFunctions
@reexport using SymPy
@reexport using QuadGK
@reexport using LinearAlgebra
@reexport using ForwardDiff
using PlotUtils

using Requires

function __init__()
     @require SimplePlots="307c2aad-90be-4152-b348-f51955fac6ce" include("simpleplots.jl")
     @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("plots.jl")
     @require AbstractPlotting="537997a7-5e4e-5d89-9595-2241ea00577e" include("makie.jl")
 end

## start it up

"""
    mth229(dirnm)
Entry point to install projects and start notebook
"""
function mth229(dirnm=homedir())
    eval(:(using IJulia))
    eval(:(import ZipFile))

    if !isfile(joinpath(dirnm, "01-calculator.ipynb"))
        @warn "installing projects in $dirnm"
        zf = "https://www.github.com/mth229/229-projects/archive/master.zip"
        zarchive = ZipFile.Reader(download(zf))
        !isdir(dirnm) && mkdir(dirnm)
        cd(dirnm)

        @show zarchive.files
        for f in zarchive.files
            nm = basename(f.name)
            occursin("ipynb", nm) || continue
            @info "installing $nm"
            open(nm, "w") do io
                write(io, read(f, String))
            end
        end
    end

    Base.invokelatest(IJulia.notebook)
end
export mth229

###
export tangent, secant, D, grad, sign_chart, fisheye, rangeclamp
export lim
export bisection, newton
export riemann, fubini
export uvec, xs_ys, unzip, parametric_grid
export e

@info "This package defines `e = exp(1)`"
"e is not Base.MathConstants.ℯ, rather `exp(1)` so that it plays nicely with ForwardDiff"
e = exp(1)

"""
    rangeclamp(f, hi=20, lo=-hi)

Replace large values of `f` with `NaN`; useful for plotting with vertical asymptotes.
"""
rangeclamp(f, hi=20, lo=-hi; replacement=NaN) = x -> lo < f(x) < hi ? f(x) : replacement

# give a warning
@info "This package overloads the `'` operation for derivatives. This may cause issues with linear algebra usage"

 " f'(x) will find the derivative of `f` using Automatic Differentation from the `ForwardDiff` package "
Base.adjoint(f::Function) = x -> ForwardDiff.derivative(f, float(x))
function D(f, n::Int=1)
    n < 0 && throw(ArgumentError("n must be non-negative"))
    n == 0 && return f
    n == 1 && return x -> ForwardDiff.derivative(f, float(x))
    D(D(f), n-1)
end
grad(f) = (x, xs...) -> ForwardDiff.gradient(f, vcat(x, xs...))
# could aslo wrap as function ForwardDiff.jacobian, ForwardDiff.hessian

"""
   sign_chart(f, a, b; atol=1e-4)

Create a sign chart for `f` over `(a,b)`. Returns a tuple with an identified zero or vertical asymptote and the corresponding sign change. The tolerance is used to disambiguate numerically found values.

# Example

```
julia> sign_chart(x -> x/(x-1)^2, -5, 5)
2-element Vector{NamedTuple{(:∞0, :sign_change), Tuple{Float64, String}}}:
 (∞0 = 0.0, sign_change = "- → +")
 (∞0 = 1.0000000000000002, sign_change = "+ → +")
```

"""
function sign_chart(f, a, b; atol=1e-6)
    pm(x) = x < 0 ? "-" : x > 0 ? "+" : "0"
    summarize(f,cp,d) = (DNE_0_∞=cp, sign_change=pm(f(cp-d)) * " → " * pm(f(cp+d)))

    if Roots._is_f_approx_0(f(a),a, eps(), eps()) ||
        Roots._is_f_approx_0(f(b), b, eps(), eps())
        return "Sorry, the endpoints must not be zeros for the function"
    end

    zs = find_zeros(f, a, b)
    pts = vcat(a, zs, b)
    for (u,v) ∈ zip(pts[1:end-1], pts[2:end])
        zs′ = find_zeros(x -> 1/f(x), u, v)
        for z′ ∈ zs′
            flag = false
            for z ∈ zs
                if isapprox(z′, z, atol=atol)
                    flag = true
                    break
                end
            end
            !flag && push!(zs, z′)
        end
    end

    if isempty(zs)
	fc = f(a + (b-a)/2)
	return "No sign change, always " * (fc > 0 ? "positive" : iszero(fc) ? "zero" : "negative")
    end

    sort!(zs)
    m,M = extrema(zs)
    d = min((m-a)/2, (b-M)/2)
    if length(zs) > 1
        d′ = minimum(diff(zs))/2
        d = min(d, d′ )
    end
    summarize.(f, zs, d)
end

"""
Returns a function describing the tangent line to the graph of f at x=c.

Example. Where does the tangent line intersect the y axis?
```
f(x) = sin(x)
tl(x) = tangent(f, pi/4)(x)  # or tl = tangent(f, pi/3) to use a non-generic function
tl(0)
```

Uses the automatic derivative of `f` to find the slope of the tangent line at `x=c`.

"""
tangent(f,c) = x -> f(c) + f'(c) * (x-c)

"""
Returns a function describing the secant line to the graph of f at x=a and x=b.

Example. Where does the secant line intersect the y axis?
```
f(x) = sin(x)
a, b = pi/4, pi/3
sl(x) = secant(f, a, b)(x)  # or sl = sl(f, a, b) to use a non-generic function
sl(0)
```

"""
secant(f, a, b) = x -> f(a) + (f(b) - f(a)) / (b-a) * (x - a)

"""

`lim(f, c, n, dir="+")`: means to generate numeric table of values of `f` as `h` gets close to `c`.

Example:
```
f(x) = sin(x) / x
lim(f, 0)
```
"""
function lim(f::Function, c::Real; n::Int=6, dir="+")
	 hs = [(1/10)^i for i in 1:n] # close to 0
	 if dir == "+"
	   xs = c .+ hs
	 else
	   xs = c .- hs
	 end
	 ys = map(f, xs)
	 [xs ys]
end

"""

Simple implementation of the bisection method.

Example:

```julia
bisection(sin, 3, 4)
f(x) = x^5 - x^4 - x^3 - x^2 - x - 1
a = bisection(f, 1, 2)
f(a)
```

The display shows a simple graphic illustrating the method's division for the first few steps.

An easier-to-understand alternative to `Roots.find_zero(f, (a,b), Bisection())`.

"""
function bisection(f::Function, a, b)
    a,b = sort([a,b])

    if f(a) * f(b) > 0
        error("[a,b] is not a bracket. A bracket means f(a) and f(b) have different signs!")
    end

    M = a + (b-a) / 2

    i, j = 0, 64
    ss = fill("#", 65)
    ss[i+1]="a"; ss[j+1]="b"
    println("")
    println(join(ss))
    flag = true

    while a < M < b
        if flag && j-i == 1
            ss = fill(" ", 65)
            ss[j:(j+1)] .= "⋮"
            println(join(ss))
            println("")
            flag = false
        end

        if f(M) == 0.0
            println("... exact answer found ...")
	    break
        end
        ## update step
	if f(a) * f(M) < 0
	    a, b = a, M

            if flag
                j = div(i + j, 2)
            end

	else
	    a, b = M, b

            if flag
                i = div(i + j, 2)
            end

	end

        if flag
            ss = fill(".", 65)
            ss[i+1]="a"; ss[j+1]="b"; ss[(i+2):j] .= "#"
            println(join(ss))
        end

        M = a + (b-a) / 2
    end
    M
end

newton(f, fp, x0; kwargs...) = Roots.find_zero((f,fp), x0, Roots.Newton(); kwargs...)
newton(f, x0; kwargs...) = newton(f, f', x0; kwargs...)

"""
    fisheye(f)

Transform `f` defined on `(-∞, ∞)` to a new function whose domain is in `(-π/2, π/2)` and range is within `(-π/2, π/2)`. Useful for finding all zeros over the real line. For example

```
f(x) = 1 + 100x^2 - x^3
fzeros(f, -100, 100) # empty just misses the zero found with:
fzeros(fisheye(f), -pi/2, pi/2) .|> tan  # finds 100.19469143521222, not perfect but easy to get
```

By Gunter Fuchs.
"""
fisheye(f) = atan ∘ f ∘ tan

##
## --------------------------------------------------
##

"""
riemann: compute Riemann sum approximations to a definite integral. As well, implement trapezoid and Simpson's rule.

Example:
```
f(x) = exp(x^2)
riemann(f, 0, 1, 1000)   # default right-Riemann sums
riemann(f, 0, 1, 1000, method="left")       # left sums
riemann(f, 0, 1, 1000, method="trapezoid")  # use trapezoid rule
riemann(f, 0, 1, 1000, method="simpsons")   # use Simpson's rule
```

"""
function riemann(f::Function, a::Real, b::Real, n::Int; method="right")
  if method == "right"
     meth = (f,l,r) -> f(r) * (r-l)
  elseif method == "left"
     meth= (f,l,r) -> f(l) * (r-l)
  elseif method == "trapezoid"
     meth = (f,l,r) -> (1/2) * (f(l) + f(r)) * (r-l)
  elseif method == "simpsons"
     meth = (f,l,r) -> (1/6) * (f(l) + 4*(f((l+r)/2)) + f(r)) * (r-l)
  end

  xs = a .+ (0:n) * (b-a)/n
  as = [meth(f, l, r) for (l,r) in zip(xs[1:end-1], xs[2:end])]
  sum(as)
end

#######
## simplified multivariable integrals

# limits of integration
endpoints(ys,x) = ((f,x) -> isa(f, Function) ? f(x...) : f).(ys, Ref(x))
# avoid specialization in quadgk
struct FWrapper
    f
end
(F::FWrapper)(x) = F.f(x)

"""
fubini(f, dy, dx)
fubini(f, dz, dy, dx)

Computes numeric integral of `f` over region specified by `dz`, `dy`, `dx`. These are a tuple of values of numbers or univariate functions depending on the value of the term on the right (`dy` can depend on `dx`s value).

*Much* slower than `hcubature` from the `HCubature` package, as it refines flat areas too many times, allocates too much, etc. But does allow a more flexible specification of the region to integrate over, as `hcubature` requires box-like regions.

```
f(x,y,z) = x * y^2 * z^3
fubini(f, (0,1), (0,2), (0,3))  # int_0^3 int_0^2 int_0^1 f(x,y,z) dz dy dx
g(v) = f(v...)
hcubature(g, (0,0,0), (3,2,1))  # same. Not order switched

# triangular like region
fubini(f, (0, y->y), (0, x->x), (0,3))
```
"""
fubini(f, dx)     = quadgk(FWrapper(f), dx...)[1]
fubini(f, ys, xs) = fubini(x -> fubini(y -> f(x,y), endpoints(ys, x)), xs)
fubini(f, zs, ys, xs) = fubini(x ->
    fubini(y ->
        fubini(z -> f(x,y,z),
            endpoints(zs, (x,y))),
        endpoints(ys,x)),
    xs)

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

uvec(x) = x / norm(x)

"""
    `unzip(vs)`
    `unzip(v1, v2, ...)`
    `unzip(r::Function, a, b)`
    `unzip(r::Function, a, b, n)`

Take a vector of points described by vectors (as returned by, say
`r(t)=[sin(t),cos(t)], r.([1,2,3])`, and return a tuple of collected x
values, y values, and optionally z values.

If the argument is specified as a comma separated collection of vectors, then these are combined and passed along.

If the argument is a function and two end point, then the  points are chosen by `PlotUtils.adaptedgrid`.

If the argument is a function, two endpoints and a number of points, then  `n` points are chosen evenly spaced over `[a,b]`.

This is useful for plotting when the data is more conveniently
represented in terms of vectors, but the plotting interface requires the x and y values collected.

Examples:
```
using SimplePlots
r(t) = [sin(t), cos(t)]
rp(t) = [cos(t), -sin(t)]
plot(unzip(r, 0, 2pi)...)  # calls plot(xs, ys)

t0, t1 = pi/6, pi/4

p, v = r(t0), rp(t0)
plot!(unzip(p, p+v)...)  # connect p to p+v with line

p, v = r(t1), rp(t1)
quiver!(unzip([p])..., quiver=unzip([v]))
```

Based on `unzip` from the `Plots` package.
"""
unzip(vs) = Tuple([vs[j][i] for j in eachindex(vs)] for i in eachindex(vs[1]))
unzip(v,vs...) = unzip([v, vs...])
unzip(r::Function, a, b, n) = unzip(r.(range(a, stop=b, length=n)))
# return (xs, f.(xs)) or (f₁(xs), f₂(xs), ...)
function unzip(f::Function, a, b)
    n = length(f(a))
    if n == 1
        return PlotUtils.adapted_grid(f, (a,b))
    else
        xsys = [PlotUtils.adapted_grid(x->f(x)[i], (a,b)) for i ∈ 1:n]
        xs = sort(vcat([xsys[i][1] for i ∈ 1:n]...))
        return unzip(f.(xs))
    end
end

# unzip(a) = map(x -> getfield.(a, x), fieldnames(eltype(a)))

# alternate, should deprecate
xs_ys(vs) = (A=hcat(vs...); Tuple([A[i,:] for i in eachindex(vs[1])]))
xs_ys(v,vs...) = xs_ys([v, vs...])
xs_ys(r::Function, a, b, n=100) = xs_ys(r.(range(a, stop=b, length=n)))

"""
    parametric_grid(us, vs, r)

Create matrices for `xs`, `ys`, `zs` from `r(u,v) = [x(u,v), y(u,v), z(u,v)]`

Used to plot parametrically defined surfaces.
"""
function parametric_grid(us, vs, r)
    unzip(r.(us, vs'))
end

# for plotif. This identifies a vector of colors
function identify_colors(g, xs, colors=(:red, :blue, :black))
    F = (a,b) -> begin
        ga,gb=g(a),g(b)
        ga * gb < 0 && return nothing
        ga >= 0 && return true
        return false
    end
    find_colors(F, xs, colors)
end

# F(a,b) returns true, false, or nothing
function find_colors(F, xs, colors=(:red, :blue, :black))
    n = length(xs)
    cols = repeat([colors[1]], n)
    for i in 1:n-1
        a,b = xs[i], xs[i+1]
        val = F(a,b)
        if val == nothing
            cols[i] = colors[3]
        elseif val
            cols[i] = colors[1]
        else
            cols[i] = colors[2]
        end
    end
    cols[end] = cols[end-1]
    cols
end

end


text:
"""
Ansillary is a package for interacting with ANSI compatible terminals.

Ansillary aims to support only the commonly supported capabilities of ANSI terminals, i.e. completely replacing ncurses is not a goal of Ansillary. Explicitly, if a capability is not available in one of the following terminals then it will not be supported in Ansillary:

* LibVTE terminals.
* LibVTerm terminals.
* The Linux console.
* Konsole.
* XTerm.

The reason for this is simple, reliable capability detection is harder than writing a package that wraps ncurses (or similar libraries such as termbox). Do note, however, that if a terminal does not support a capability but the lack of that capability does not cause issues for the user then it will be supported. For example, bracketed paste mode is supported despite the Linux console not actually having a paste feature because the Linux console still parses and accepts the required escape sequences even though they have no effect. Whereas scrolling is _not_ supported because the Linux console does not support scrolling and not having something scroll is something that the user will definitely notice.

The package is currently split into three modules: [`Cursor`](@ref) for controlling the position and visibilty of the cursor, [`Inputs`](@ref) for reading input events from a terminal, and [`Screen`](@ref) for controlling what is displayed.
"""
module Ansillary

using REPL.Terminals: TTYTerminal

export Cursor, Inputs, Screen, TTYTerminal

const TERMINAL = Ref{TTYTerminal}()

function __init__()
    TERMINAL[] =
        TTYTerminal(get(ENV, "TERM", Sys.iswindows() ? "" : "dumb"), stdin, stdout, stderr)
end

# Needs to be at top because `Cursor` needs it.
include("Inputs.jl")

include("Cursor.jl")
include("Screen.jl")

end # module


text:
module TensorNetworkStatesBackend

    include("./TensorNetworkStatesFunctions.jl")
    include("../backend_Models.jl")
    export execute_backend
    using .TensorNetworkStatesFunctions
    using .BackendModels

    function execute_backend(n_qregs::Int, gates, model, meas_all::Bool)
        _mps = execute_operation(n_qregs, gates, model)
        output_state = Vector{ComplexF64}(undef, 0)
        for i=0:2^n_qregs-1
            res = restore(_mps, i) # MPS で表された状態から、| bin(i) > の係数を計算する関数
            push!(output_state, res)
        end
        info = mps_size(_mps)
        return (output_state, info)
    end

    function execute_backend(n_qregs::Int, gates, model, measure_basis::Array{Int, 1})
        if length(measure_basis) != n_qregs
            error("Invalid measure basis.")
        end
        _mps = execute_operation(n_qregs, gates, model)
        meas_idx = 0
        for (i, val) in enumerate(measure_basis)
            if val == 1
                meas_idx += 2^(n_qregs - i)
            elseif val > 1
                error("Invalid measure basis.")
            end
        end
        res = restore(_mps, meas_idx)
        info = mps_size(_mps)
        return (res, info)
    end

    function execute_operation(n_qregs::Int, gates, model)
        if model.zero_state == true
            mps = init_zero_MPS(n_qregs)
        else
            if length(model.init_state) != 2^n_qregs
                error("Initial state length must be equal to 2^(n_qubits).")
            else
                mps = init_MPS(model.init_state, model.param)
            end
        end

        backend_gates = call_backend_gates(gates)
        _backend_gates = copy(backend_gates)
        parsed_gates = parse_gates(_backend_gates)
        _mps = copy(mps)

        for gate in parsed_gates
            _mps = apply!(gate, _mps, model.param)
        end
        return _mps
    end

    function call_backend_gates(gates)
        backend_gates = []
        for gate in gates
            name = gate._name
            fields = fieldnames(typeof(gate))
            field_names = []
            for i in fields
                push!(field_names, string(i))
            end
            given_property = [getproperty(gate, property) for property in fields]
            code = gate._name * "("
            for prop in given_property[1:end-1]
                code *= string(prop)
                code *= ", "
            end
            code *= ")"
            new_gate = eval(Meta.parse(code))
            push!(backend_gates, new_gate)
        end
        return backend_gates
    end

    function parse_gates(backend_gates)
        parsed_gates =[]
        for gate in backend_gates
            _gates = gate_parser(gate)
            for _gate in _gates
                push!(parsed_gates, _gate)
            end
        end
        return parsed_gates
    end

end

max_stars_repo_name: ryuNagai/blueqat-julia

max_issues_repo_name: ryuNagai/blueqat-julia

max_forks_repo_name: ryuNagai/blueqat-julia


text:
using ReinforcementLearning
using StableRNGs
using Flux
using Flux.Losses

function RL.Experiment(
    ::Val{:JuliaRL},
    ::Val{:BCQD},
    ::Val{:CartPole},
    type::AbstractString;
    save_dir = nothing,
    seed = 123,
)
    rng = StableRNG(seed)
    UPDATE_FREQ = 1
    env = CartPoleEnv(; T = Float32, rng = rng)
    ns, na = length(state(env)), length(action_space(env))

    trajectory_num = 10000
    dataset_size = 10000
    batch_size = 64

    init = glorot_uniform(rng)

    create_actor_critic() = ActorCritic(
        actor = Chain(
            Dense(ns, 128, relu; init),
            Dense(128, 128, relu; init),
            Dense(128, na; init),
        ),
        critic = Chain(
            Dense(ns, 128, relu; init),
            Dense(128, 128, relu; init),
            Dense(128, na; init),
        ),
        optimizer = ADAM(1e-3),
    )

    agent = Agent(
        policy = OfflinePolicy(
            learner = BCQDLearner(
                approximator = create_actor_critic() |> cpu,
                target_approximator = create_actor_critic() |> cpu,
                γ = 0.99f0,
                τ = 0.01f0,
                θ = 1f-2,
                threshold = 0.3f0,
                batch_size = batch_size,
                update_freq = UPDATE_FREQ,
            ),
            dataset = gen_JuliaRL_dataset(:BasicDQN, :CartPole, type; dataset_size = dataset_size),
            continuous = false,
            batch_size = batch_size,
        ),
        trajectory = CircularArraySARTTrajectory(
            capacity = 1000,
            state = Vector{Float32} => (ns,),
        ),
    )

    stop_condition = StopAfterStep(trajectory_num, is_show_progress=!haskey(ENV, "CI"))
    hook = TotalRewardPerEpisode()
    Experiment(agent, env, stop_condition, hook, "BCQD <-> CartPole ($type dataset)")
end

using Plots
pyplot() #hide
ex = E`JuliaRL_BCQD_CartPole(medium)`
run(ex)
plot(ex.hook.rewards)
savefig("assets/JuliaRL_BCQD_CartPole_medium.png") #hide

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

"""
    juldate(date::DateTime) -> reduced_julia_days

### Purpose ###

Convert from calendar to Reduced Julian Days.

### Explanation ###

Julian Day Number is a count of days elapsed since Greenwich mean noon on 1
January 4713 B.C.  Julian Days are the number of Julian days followed by the
fraction of the day elapsed since the preceding noon.

This function takes the given `date` and returns the number of Julian calendar
days since epoch `1858-11-16T12:00:00` (Reduced Julian Days = Julian Days -
2400000).

### Argument ###

* `date`: date in Julian Calendar, UTC standard.  It can be either e single date
  or an array of dates.  Each element can be given in `DateTime` type or
  anything that can be converted to that type.  In the case of vectorial input,
  each element is considered as a date, so you cannot provide a date by parts.

### Output ###

The number of Reduced Julian Days is returned.  If `date` is an array, an array
of the same length as `date` is returned.

### Example ###

Get number of Reduced Julian Days at 2016-03-20T15:24:00.

``` julia
juldate(DateTime(2016, 03, 20, 15, 24))
# => 57468.14166666667
juldate(2016, 03, 20, 15, 24)
# => 57468.14166666667
juldate("2016-03-20T15:24")
# => 57468.14166666667
```

### Notes ###

Julian Calendar is assumed, thus before `1582-10-15T00:00:00` this function is
*not* the inverse of `daycnv`.  For the conversion proleptic Gregorian date to
number of Julian days, use `jdcnv`, which is the inverse of `daycnv`.

Code of this function is based on IDL Astronomy User's Library.
"""
# Before 1582-10-15 Dates.datetime2julian uses proleptic Gregorian Calendar,
# instead AstroLib's juldate uses Julian Calendar.  In addition, after that day
# using Dates.datetime2julian(date)-2.4e6 would give results not precise enough.
# For all these reasons, we use here the same algorithm as AstroLib's juldate.
function juldate(dt::DateTime)
    year, month, day = Dates.yearmonthday(dt)
    hours, minutes, seconds = Dates.hour(dt), Dates.minute(dt), Dates.second(dt)

    if year < 0
        year += 1
    else
        year == 0 && error("There is no year zero in Julian Calendar")
    end

    day = day + (hours + minutes/60.0 + seconds/3600.0)/24.0

    # Do not include leap year in January and February.
    if month < 3
        year -= 1
        month += 12
    end

    jd = fld(year, 4.0) + 365.0*(year - 1860.0) + floor(30.6001*(month + 1.0)) +
         day - 105.5

    # Adjust for Gregorian Calendar, started on 1582-10-15 (= RJD -100830.5).
    if jd > -100830.5
        a = div(year, 100)
        jd += 2.0 - a + fld(a, 4.0)
    end
    return jd
end

juldate(dt...) = juldate(DateTime(dt...))

function juldate{T<:Any}(dt::AbstractArray{T})
    jd = similar(dt, Float64)
    for i in eachindex(dt)
        jd[i] = juldate(dt[i])
    end
    return jd
end