text:
module Pseudospectra
#=
Eigenvalue and Pseudospectrum Analysis for Julia

The Pseudospectra.jl package is a translation of EigTool, but no endorsement
or promotion by the authors of EigTool is implied.

This package is released under a BSD license, as described in the LICENSE file.

Julia code and supplements
Copyright (c) 2017-2019, Ralph Smith

Portions derived from EigTool:

 Copyright (c) 2002-2014, The Chancellor, Masters and Scholars
 of the University of Oxford, and the EigTool Developers. All rights reserved.
 EigTool is maintained on GitHub:  https://github.com/eigtool

SPDX-License-Identifier: BSD-3-Clause
License-Filename: LICENSES/BSD-3-Clause_Eigtool
=#

using ProgressMeter

using LinearAlgebra, SparseArrays, Arpack, Printf
using Requires

export new_matrix, driver!, spectralportrait
export psa_compute, psa_radius, psa_abscissa
export numerical_range, numerical_abscissa
export modeplot, mtxexpsplot, mtxpowersplot, isheadless, iscomputed
export PSAStruct, ArpackOptions, Portrait, GUIState

# Plotting packages should probably extend these:
export zoomin!, zoomout!

# Not exported, but may be used by plotting packages:
# vec2ax, expandlevels, isvalidax
# oneeigcond, psmode_inv_lanczos, transient_bestlb, set_method!

# Associated plotting packages should provide these, specialized on their
# own GUIState types:
# redrawcontour, surfplot, arnoldiplotter!, ewsplotter, plotmode,
# replzdlg, addmark

const smallσ = 1e-150

"""
by default, sparse matrices of this size or smaller are converted to full
for pseudospectra computation.
"""
const nmax4autofull = 200
"""
by default, iterative methods are used for computing pseudospectra of dense
matrices of this size or larger.
"""
const nmin4autoiter = 1600

const myname = "PSA"

include("types.jl")

# Placeholders for plot-specific code implemented elsewhere
function redrawcontour end
function surfplot end
function arnoldiplotter! end
function _portrait end

"""
    ewsplotter(gs::GUIState, ews::Vector, zoom)

plot eigenvalues

So we have something to look at while waiting for the compute engines.
"""
function ewsplotter end
function plotmode end
function replzdlg end

function addmark end
"""
    mtxexpsplot(ps_data,dt=0.1,nmax=50; gs::GUIState = defaultgs(), gradual=false)

plot the evolution of `∥e^(tA)∥`.

This is useful for analyzing linear initial value problems `∂x/∂t = Ax`.
"""
function mtxexpsplot(ps_data::PSAStruct, dt=0.1, nmax=50;
                     gs::GUIState=defaultgs(), kws...)
    mtxexpsplot(gs, ps_data, dt=dt, nmax=nmax; kws...)
end

"""
    mtxpowersplot(ps_data, nmax=50; gs::GUIState = defaultgs(), gradual=false)

plot norms of powers of a matrix `∥A^k∥`

This is useful for analyzing iterative linear algebra methods.
"""
function mtxpowersplot(ps_data::PSAStruct, nmax=50;
                       gs::GUIState = defaultgs(), kws...)
    mtxpowersplot(gs, ps_data, nmax=nmax; kws...)
end

function fillopts end
function isheadless end

include("utils.jl")
include("compute.jl")
include("xeigs.jl")
include("modes.jl")
include("abscissa.jl")
include("radius.jl")
include("numrange.jl")
include("transients.jl")
include("plotter.jl")
include("zooming.jl")

include("plots/PSAPlots.jl")

"""
    new_matrix(A::AbstractMatrix, opts::Dict{Symbol,Any}=()) -> ps_data

process a matrix into the auxiliary data structure used by Pseudospectra.

# Options
- `:direct::Bool`: force use of a direct algorithm?
- `:keep_sparse::Bool`: use sparse matrix code even if `A` is not large?
- `:real_matrix::Bool`: treat `A` as unitarily equivalent to a real matrix?
- `:verbosity::Int`: obvious
- `:eigA`: eigenvalues of `A`, if already known
- `:proj_lev`: projection level (see `psa_compute`)
- `:npts`: edge length of grid for computing and plotting pseudospectra
- `:arpack_opts::ArpackOptions`: (see type description)
- `:levels::Vector{Real}`: contour levels
- `:ax::Vector{Real}(4)`: bounding box for computation `[xmin,xmax,ymin,ymax]`
- `:scale_equal::Bool`: force isotropic axes for spectral portraits?
"""
function new_matrix(A::AbstractMatrix,
                    opts::Dict{Symbol,Any}=Dict{Symbol,Any}())
    m,n=size(A)
    (m >= n) || throw(ArgumentError(
        "Only square or tall rectangular matrices are supported."))
    (issparse(A) && (m != n)) && throw(ArgumentError(
        "Only square sparse matrices are supported."))
    (any(isnan.(A)) || any(isinf.(A))) && throw(ArgumentError(
        "Input matrix has infinite or invalid entries."))

    # flag for the M x (M-1) Hessenberg form
    # Presumably intended for the case where projection is done
    # by a Krylov scheme outside this package.
    AisHess = ((m == (n+1)) && all([x == 0 for x in tril(A,-2)]))

    # User may specify that A is unitarily equivalent to a real matrix
    # even if it is complex
    Aisreal = get(opts,:real_matrix, !(eltype(A) <: Complex))

    verbosity = get(opts,:verbosity,1)

    convert2full = issparse(A) & (n <= nmax4autofull) &
        !get(opts,:keep_sparse,false)
    if haskey(opts,:direct)
        direct = opts[:direct]
    else
        if issparse(A)
            direct = convert2full
        else
            direct = (n <= nmin4autoiter)
            if (verbosity > 0) && (n > nmin4autoiter)
                # Might merit a warning here since it is most likely not
                # expected so iteration options are probably inappropriate.
                # For now, unspec. iteration options already => a warning.
                println("defaulting to iterative for large dense mtx")
            end
        end
    end

    # TODO: check for correctness of
    # proj_lev, levels, ax, arpack stuff, etc. from opts

    # define as placeholder if not provided
    eigA = get(opts,:eigA,Vector{complex(eltype(A))}())

    input_unitary_mtx = get(opts,:unitary_mtx,I)
    proj_lev = get(opts,:proj_lev,Inf)
    npts = get(opts,:npts,setgridsize(n,24,80,!Aisreal))
    Aissquare = (m == n)

    local Tschur, U
    haveschur = false
    if Aissquare && !issparse(A) && direct
        # if small, delay is negligible
        (verbosity > 1) && (m > 100) &&
            println("Attempting initial decomposition...")
        # If square, dense, & direct, we prefer a Schur factorization.
        # Checking for schurfact! method should work,
        # but that's just asking for surprises. This should be robust.
        try
            if eltype(A) <: Complex
                F = schur(A)
            else
                # For some reason Julia devs think a real Schur decomp
                # should shadow the true (not real!) thing
                F = schur(A .+ zero(eltype(A))*im)
            end
            Tschur,U,eigA  = F.T,F.Z,F.values
            haveschur = true
        catch JE
            isa(JE,MethodError) || rethrow(JE)
        end
        if !haveschur && isempty(eigA)
            try
                eigA = eigvals(A)
            catch JE
                isa(JE,MethodError) || rethrow(JE)
                # ME: maybe trap algorithmic errors too; what are they?

                # Warning is needed here since it explains why we need axes
                # for the driver.
                @warn("Failed to compute eigenvalues; proceeding without.")
                if verbosity > 0
                    # If we display(JE) we get the whole damn matrix too
                    println("Exception was method error: ",JE.f,
                            " for ",typeof(A))
                end
            end
        end
        (verbosity > 1) && (m > 100) && println("...done.")
    end
    if haveschur
        ps_dict = Dict{Symbol,Any}(:Aisreal => Aisreal,
                                   :isHessenberg => AisHess,
                                   :schur_mtx => Tschur,
                                   :schur_unitary_mtx => U,
                                   :direct => true,
                                   :projection_on => true,
                                   :proj_lev => proj_lev,
                                   :ews => eigA)
        ps_data = PSAStruct(UpperTriangular(Tschur), input_unitary_mtx * U,
                            A, input_unitary_mtx, ps_dict)

    elseif issparse(A) || AisHess || !direct || Aissquare
        # sparse, Hessenberg, iterative, or needing special handling
        # CHECKME: previously seemed to force
        #   direct |= Aissquare
        #
        ps_dict = Dict{Symbol,Any}(:Aisreal => Aisreal,
                                   :isHessenberg => AisHess,
                                   :projection_on => false,
                                   :proj_lev => Inf,
                                   :ews => eigA)
        ps_data = PSAStruct(A, input_unitary_mtx, A, input_unitary_mtx,
                            ps_dict)
        if !direct
            if !haskey(opts,:arpack_opts)
                @info("setting default ARPACK options")
                ps_dict[:arpack_opts] = ArpackOptions{eltype(A)}()
            else
                isa(opts[:arpack_opts],ArpackOptions) || throw(
                    ArgumentError("type of opts[:arpack_options] must "
                                  * "be ArpackOptions"))
                ps_dict[:arpack_opts] = opts[:arpack_opts]
            end
        elseif issparse(A) && convert2full
            (verbosity > 0) &&
                println("converting to full for direct computation")
            Atmp = full(ps_data.matrix)
            try
                F  = schur(Atmp+complex(eltype(Atmp))(0))
                ps_dict[:schur_mtx] = F.T
                ps_dict[:schur_unitary_mtx] = F.Z
                ps_data.matrix = UpperTriangular(F.T)
                ps_data.unitary_mtx = ps_data.input_unitary_mtx * F.T
                ps_dict[:projection_on] = true
                ps_dict[:ews] = F.values
                ps_dict[:orig_ews] = copy(F.values)
            catch
                ps_data.matrix = Atmp
            end
        end
    else # dense, non-square (but not Hessenberg), and direct
        rfS,rfT = rect_fact(A)
        ps_dict = Dict{Symbol,Any}(:Aisreal => Aisreal,
                                   :isHessenberg => false,
                                   :projection_on => false,
                                   :proj_lev => Inf,
                                   :matrix2 => rfT,
                                   :ews => eigA)
        ps_data = PSAStruct(rfS,input_unitary_mtx,A,input_unitary_mtx,
                            ps_dict)
    end
    ps_dict[:orig_ews] = eigA
    ps_dict[:ew_estimates] = false

    if !ps_dict[:isHessenberg] && !isa(ps_data.unitary_mtx,UniformScaling)
        if size(ps_data.unitary_mtx,2) ∉ [m,1]
            ps_data.unitary_mtx = I
        end
    end

    lev = get(opts,:levels,zeros(0))
    zoom = Portrait(zeros(0),zeros(0),zeros(0,0),
                    npts, get(opts,:ax,zeros(0)),
                    LevelDesc(lev),
                    isempty(lev), proj_lev,
                    size(ps_data.matrix), false,
                    get(opts,:scale_equal,false))
    push!(ps_data.zoom_list,zoom)
    ps_data.zoom_pos = 1
    # save for use when returning to initial plot
    ps_dict[:init_opts] = deepcopy(zoom)
    ps_dict[:init_direct] = direct
    ps_dict[:direct] = direct
    ps_dict[:verbosity] = verbosity

    # DEVNOTE: if direct, upstream constructs axes and calls origplot/redraw
    return ps_data
end

"""
    new_matrix(A, opts::Dict{Symbol,Any}=()) -> ps_data

process a linear operator object into the auxiliary data structure used by
Pseudospectra.

There must be methods with `A` for `eltype`, `size`, and `mul!`.
It is up to the user to make sure that `mul!` is consistent with any
options passed to the iterative solver (see documentation for [`xeigs`](@ref)).
"""
function new_matrix(A, opts::Dict{Symbol,Any}=Dict{Symbol,Any}())
    # CHECKME: can A be anything other than a LinearMap here?
    m,n=size(A)
    (m == n) || throw(ArgumentError(
        "Only square linear operators are supported."))
    Aisreal = get(opts,:real_matrix, !(eltype(A) <: Complex))

    verbosity = get(opts,:verbosity,1)
    direct = false
    eigA = get(opts,:eigA,Vector{complex(eltype(A))}())

    input_unitary_mtx = get(opts,:unitary_mtx,I)
    npts = get(opts,:npts,setgridsize(n,24,80,!Aisreal))
    proj_lev = get(opts,:proj_lev,Inf)

    ps_dict = Dict{Symbol,Any}(:Aisreal => Aisreal,
                                   :isHessenberg => false,
                                   :projection_on => false,
                                   :proj_lev => Inf,
                                   :ews => eigA)
    ps_data = PSAStruct(A, input_unitary_mtx, A, input_unitary_mtx,
                        ps_dict)
    if !haskey(opts,:arpack_opts)
        @info("setting default ARPACK options")
        ps_dict[:arpack_opts] = ArpackOptions{eltype(A)}()
    else
        isa(opts[:arpack_opts],ArpackOptions) || throw(
            ArgumentError("type of opts[:arpack_options] must "
                          * "be ArpackOptions"))
        ps_dict[:arpack_opts] = opts[:arpack_opts]
    end
    lev = get(opts,:levels,zeros(0))
    zoom = Portrait(zeros(0),zeros(0),zeros(0,0),
                    npts, get(opts,:ax,zeros(0)),
                    LevelDesc(lev),
                    isempty(lev), proj_lev,
                    size(ps_data.matrix), false,
                    get(opts,:scale_equal,false))
    push!(ps_data.zoom_list,zoom)

    ps_dict[:orig_ews] = eigA
    ps_dict[:ew_estimates] = false
    ps_data.zoom_pos = 1
    # save for use when returning to initial plot
    ps_dict[:init_opts] = deepcopy(zoom)
    ps_dict[:init_direct] = direct
    ps_dict[:direct] = direct
    ps_dict[:verbosity] = verbosity

    return ps_data
end

# for verifying that tests cover intended cases
const logging_algo = Ref{Bool}(false)

"""
    driver!(ps_data, opts, gs=defaultgs(); revise_method=false)

Compute pseudospectra and plot a spectral portrait.

If using an iterative method to get some eigenvalues and a projection, invokes
that first.

# Arguments
- `ps_data::PSAStruct`: ingested matrix, as processed by `new_matrix`
- `gs::GUIState`: object handling graphical output
- `opts::Dict{Symbol,Any}`:
  - `:ax`, axis limits (overrides value stored in `ps_data`).
  - other options passed to `redrawcontour`, `arnoldiplotter!`

When revising a spectral portrait (`revise_method==true`), the following
entries in `opts` also apply:
 - `:arpack_opts::ArpackOptions`,
 - `:direct::Bool`.
"""
function driver!(ps_data::PSAStruct,
                 optsin::Dict{Symbol,Any}=Dict{Symbol,Any}(),
                 gs::GUIState=defaultgs();
                 myprintln=println, logger=:default, revise_method=false)
    # DEVNOTE: mostly corresponds to switch_redraw.m in EigTool
    opts = fillopts(gs,optsin)
    ps_dict = ps_data.ps_dict
    verbosity = get(ps_dict,:verbosity,1)

    # For changing from direct to iterative, or vice versa,
    if revise_method & haskey(opts,:direct)
        set_method!(ps_data, opts[:direct])
    end

    if revise_method & haskey(opts,:arpack_opts) & !ps_dict[:direct]
        ao = opts[:arpack_opts]
        if !isa(ao,ArpackOptions)
            @mywarn(logger,"invalid :arpack_opts option")
            return nothing
        end
        if haskey(ps_dict,:arpack_opts)
            pvalid = (ao == ps_dict[:arpack_opts])
        else
            pvalid = false
        end
        ps_dict[:proj_valid] = pvalid
        ps_dict[:arpack_opts] = ao
    end

    # if caller specifies ax, use it or bust.
    if haskey(opts,:ax)
        if isvalidax(opts[:ax])
            new_ax = opts[:ax]
        else
            @mywarn(logger,"opts[:ax] is not a valid bounding box")
            return nothing
        end
    else
        new_ax = zeros(0)
    end

    if ps_dict[:direct] || get(ps_dict,:proj_valid,false)
        # for iterative methods, we get here on reentrance with the projection
        n,m = size(ps_data.matrix)
        A = ps_data.matrix
        B = get(ps_dict,:matrix2,I)
        eigA = ps_dict[:ews]
        zoom = ps_data.zoom_list[ps_data.zoom_pos]

        if isempty(new_ax)
            if !isempty(eigA)
                # This sets the default domain for the typical case
                isempty(zoom.ax) && (zoom.ax = vec2ax(eigA))
                if !isheadless(gs)
                    # show eigenvalues while waiting
                    ewsplotter(gs, eigA, zoom)
                end
            else
                if isempty(zoom.ax)
                    @mywarn(logger,"bounding box must be specified")
                    return nothing
                end
            end
        else
            zoom.ax = new_ax
        end

        if haskey(opts, :npts)
            new_npts = opts[:npts]
            if isa(new_npts, Integer) && (new_npts > 2) && (new_npts < 2049)
                zoom.npts = new_npts
            else
                @mywarn(logger,"opts[:npts] is not a valid number of points")
                return nothing
            end
        end

        psa_opts = Dict{Symbol,Any}(:levels=>expandlevels(zoom.levels),
                                    :recompute_levels=>zoom.autolev,
                                    :proj_lev=>zoom.proj_lev,
                                    :scale_equal=>zoom.scale_equal,
                                    :real_matrix=>ps_dict[:Aisreal],
                                    :verbosity=>verbosity)
        ss = size(A)
        Z,x,y,t_levels,err,Tproj,eigAproj,algo = psa_compute(A,zoom.npts,
                                                             zoom.ax,
                                                             eigA,psa_opts,
                                                             B,
                                                        myprintln=myprintln,
                                                        logger=logger)

        # FIXME: handle projection properly
        ps_dict[:proj_ews] = eigAproj

        if err != 0
            @mywarn(logger,"PSA computation failed")
            # FIXME: reset GUI if any
            return nothing
        end

        (logging_algo[] | (verbosity > 1)) && println("algorithm: ",algo)
        zoom = ps_data.zoom_list[ps_data.zoom_pos]
        if zoom.autolev
            (verbosity > 1) && println("levels: $t_levels")
            zoom.levels = LevelDesc(t_levels)
        end
        zoom.x = x
        zoom.y = y
        zoom.Z = Z
        zoom.computed = true
        zoom.dims = size(ps_data.matrix)

        redrawcontour(gs, ps_data, opts)
    else
        # Iterative method (uses ARPACK):
        # This performs implicitly restart Arnoldi steps to
        # project on a Krylov subspace, yielding a Hessenberg matrix
        # with approximately the same spectral properties (locally) as A.

        ps_data.matrix = ps_data.input_matrix
        m,n = size(ps_data.matrix)

        ao = ps_dict[:arpack_opts]

        function xeigsproducer(chnl)
            local ews,H,V
            local nconv,niter,nmult
            try
                ews,v,nconv,niter,nmult,resid,H,V = xeigs(ps_data.matrix,I,chnl,
                                                      nev=ao.nev,ncv=ao.ncv,
                                                      which=ao.which,tol=ao.tol,
                                                      maxiter=ao.maxiter,
                                                      v0=ao.v0,
                                                      sigma=ao.sigma,
                                                      options=opts)

           catch JE
            # FIXME: post a dialog, reset GUI if any
                @warn("xeigs throws:")
                display(JE)
                println()
                stuff = (:failure,nothing)
                put!(chnl,stuff)
                return nothing
            end
        end

        local ews,H,V
        local nconv,niter,nmult
        old_ax = zeros(0)
        (verbosity > 1) && println("calling xeigs w/ $ao")
        chnl = Channel(xeigsproducer)
        xeigs_result = take!(chnl)
        ap_state = nothing
        while xeigs_result[1] ∉ [:finale,:failure]
            the_key,dispvec,the_str,the_ews,the_shifts = xeigs_result
            if !isheadless(gs)
                ap_state = arnoldiplotter!(gs,old_ax,opts,dispvec,
                                           the_str,the_ews, the_shifts,
                                           ap_state)
            end # if gs
            xeigs_result = take!(chnl)
        end
        if xeigs_result[1] == :failure
            @mywarn(logger,"xeigs failed")
            return nothing
        end
        the_key,ews,v,nconv,niter,nmult,resid,H,V = xeigs_result
        if verbosity > 0
            println("xeigs: $nconv of $(ao.nev) converged in $niter iters "
                    * "($nmult MxV)")
        end
        if verbosity > 1
            println("xeigs ews:")
            display(ews); println()
        end
        ews = filter(x->!isnan(x), ews)

        # We basically replace A with H, saving some projection information,
        # and proceed with the dense matrix algorithms.

        ps_dict[:ew_estimates] = true
        ps_dict[:proj_matrix] = H
        ps_data.matrix = H
        ps_dict[:isHessenberg] = true
        ps_dict[:proj_unitary_mtx] = ps_data.input_unitary_mtx * V
        ps_data.unitary_mtx = ps_dict[:proj_unitary_mtx]
        ps_dict[:proj_valid] = true
        ps_dict[:ews] = ews
        # reset zoom list
        ps_data.zoom_pos = 1
        resize!(ps_data.zoom_list,1)
        # CHECKME: do we need remove() here?
        ps_dict[:mode_markers] = []
        zoom = ps_data.zoom_list[1]

        if isempty(new_ax)
            origax = ps_dict[:init_opts].ax # init_opts is a Portrait!
            if !isempty(origax) && isvalidax(origax) # && !ps_dict[:init_direct]
                copy!(zoom.ax,origax)
            else
                # CHECKME: maybe use init_ews if available?
                zoom.ax = vec2ax(ews)
          # elseif gs.mainph != nothing
          #    println("using eigvals for axis limits")
          #    copy!(zoom.ax,getxylims(gs.mainph))
            end
        else
            zoom.ax = new_ax
        end
        zoom.autolev = ps_dict[:init_opts].autolev
        zoom.levels = deepcopy(ps_dict[:init_opts].levels)

        delete!(ps_dict,:proj_axes)
        delete!(ps_dict,:comp_proj_lev)
        origplot!(ps_data,opts,gs) # WARNING: reenters driver!()
        # caller must reset GUI if appropriate
    end
    nothing
end

function iscomputed(ps_data::PSAStruct, idx=ps_data.zoom_pos)
    ps_data.zoom_list[idx].computed
end

iscomputed(zoom::Portrait) = zoom.computed

"""
Make sure zoom list is ok, then redraw (unless `ax_only`).

Note: truncates zoom list, so use for a new problem or for a reset.
"""
function origplot!(ps_data::PSAStruct, opts, gs; ax_only = false)
    ps_data.zoom_pos = 1
    ps_dict = ps_data.ps_dict
    resize!(ps_data.zoom_list,1)
    zoom = ps_data.zoom_list[1]
    if isempty(zoom.ax)
        if isempty(get(ps_dict,:ews,[]))
            @warn("origplot called w/o preset axes or eigenvalues")
        else
            zoom.ax = vec2ax(ps_dict[:ews])
        end
    end
    ax_only || driver!(ps_data,opts,gs)
    nothing
end

"""
possibly change from direct to iterative method or vice versa
"""
function set_method!(ps_data::PSAStruct, todirect::Bool)
    # this is the part of switch_method which pertains to ps_data
    ps_dict = ps_data.ps_dict
    (todirect == ps_dict[:direct]) && return
    m,n = size(ps_data.input_matrix)
    if todirect
        if haskey(ps_dict,:schur_matrix)
            ps_data.matrix = ps_dict[:schur_matrix]
            ps_data.ews = ps_dict[:orig_ews]
            ps_dict[:ew_estimates] = false
        elseif m==n && issparse(ps_data.input_matrix)
            ps_data.matrix = ps_data.input_matrix
        end
        if haskey(ps_dict,:schur_unitary_mtx)
            ps_data.unitary_mtx  = ps_data.input_unitary_mtx *
                ps_data.schur_unitary_mtx
        else
            ps_data.unitary_mtx = ps_data.input_unitary_mtx
        end
        ps_dict[:proj_valid] = false
        # if reverting to a square matrix, no longer have ARPACK projection
        ss = size(ps_data.matrix)
        (ss[1]==ss[2]) && (ps_dict[:isHessenberg] = false)
    else # switch to iterative
        (m == n) || throw(ArgumentError("Iterative method not implemented "
                                        * "for rectangular matrices"))
        # apparently that's all we need for now
    end
    ps_dict[:direct] = todirect
end

"""
    spectralportrait(A::AbstractMatrix; npts=100) => Plots object

compute pseudospectra of matrix `A` and display as a spectral portrait.

Pseudospectra are computed on a grid of `npts` by `npts` points in
the complex plane, including a neighborhood of the spectrum.
Contour levels are `log10(ϵ)` where `ϵ` is the inverse resolvent norm.
This is a convenience wrapper for simple cases; see the Pseudospectra
package documentation for more elaborate interfaces.
"""
function spectralportrait(A0 :: AbstractMatrix; npts=100)
    if _currentplotter[] == :undef
        setpsplotter()
    end
    local ps_data
    try
        ps_data = new_matrix(A0)
    catch JE
        @warn "The spectralportrait function only works for simple cases."
        rethrow(JE)
    end
    n,m = size(ps_data.matrix)
    A = ps_data.matrix
    ps_dict = ps_data.ps_dict
    B = get(ps_dict,:matrix2,I)
    eigA = ps_dict[:ews]
    zoom = ps_data.zoom_list[ps_data.zoom_pos]
    isempty(zoom.ax) && (zoom.ax = vec2ax(eigA))
    psa_opts = _basic_psa_opts(zoom,ps_dict)
    ss = size(A)
    Z,xs,ys,t_levels,err,Tproj,eigAproj,algo = psa_compute(A,npts,
                                                             zoom.ax,
                                                             eigA,psa_opts,
                                                             B)
    return _portrait(xs,ys,Z,eigA)
end

_basic_psa_opts(zoom,ps_dict) = Dict{Symbol,Any}(
    :levels=>expandlevels(zoom.levels),
    :recompute_levels=>zoom.autolev,
    :proj_lev=>zoom.proj_lev,
    :scale_equal=>zoom.scale_equal,
    :real_matrix=>ps_dict[:Aisreal],
    :verbosity=>0)

################################################################
# FIXME: until we think of a better way to handle this:
include("../examples/demo_mtx.jl")

function __init__()
    @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" link_pyplot()
    @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" link_plots()
    # GLMakie="e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
    @require AbstractPlotting = "537997a7-5e4e-5d89-9595-2241ea00577e" link_makie()
end

end # module


text:
#-------------------------------------------------------------------------------
#caching function for gaussian quadrature of multiple streams over the azimuth

@memoize function streamnodes(n::Int64)::NTuple{2,Vector{Float64}}
    #gauss-legendre quadrature points and weights in [-1,1]
    x, w = gausslegendre(n)
    #map to θ ∈ [0,2π] and precompute trig factors
    𝓂 = Vector{Float64}(undef,n)
    𝒲 = Vector{Float64}(undef,n)
    for i ∈ 1:n
        #azimuth angle of stream
        θᵢ = (π/2)*(x[i] + 1)/2
        #mapped weight
        wᵢ = (π/2)*w[i]/2
        #factor of relevance to optical depth and such
        𝓂[i] = 1/cos(θᵢ)
        #weight with precomputed trig factors for hemispheric geometry
        𝒲[i] = 2π*wᵢ*cos(θᵢ)*sin(θᵢ)
    end
    return 𝓂, 𝒲
end

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

copyprofile(x::MoistAdiabat) = copy(x)
copyprofile(x::AtmosphericProfile) = copy(x)
copyprofile(x::OneDimensionalInterpolator) = copy(x)
copyprofile(x) = x
copyprofiles(x...) = map(copyprofile, x)

#-------------------------------------------------------------------------------
#little types to dispatch on numerical cores/methods

export Radau, Discretized

abstract type AbstractNumericalCore end

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

struct Radau <: AbstractNumericalCore
    nstream::Int64
    tol::Float64
end

function Radau(; nstream::Int=5, tol::AbstractFloat=1e-5)
    Radau(nstream, tol)
end

function Base.show(io::IO, r::Radau)
    print(io, "Radau(nstream=$(r.nstream), tol=$(r.tol))")
end

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

struct Discretized <: AbstractNumericalCore
    nstream::Int64
    nlobatto::Int64
end

function Discretized(; nstream::Int=5, nlobatto::Int=2)
    Discretized(nstream, nlobatto)
end

function Base.show(io::IO, d::Discretized)
    print(io, "Discretized(nstream=$(d.nstream), nlobatto=$(d.nlobatto))")
end

#-------------------------------------------------------------------------------
#a simple container of arrays needed for whole atmosphere characterization

export FluxPack

struct FluxPack{T<:Real}
    τ::Matrix{T}
    M⁺::Matrix{T}
    M⁻::Matrix{T}
    F⁺::Vector{T}
    F⁻::Vector{T}
    Fnet::Vector{T}
end

function Base.show(io::IO, M::FluxPack{T}) where {T}
    np, nν = size(M.M⁺)
    print(io, "size ($np x $nν) FluxPack{$T}\n")
    print("  TOA:\n")
    print(io, "    outgoing radiation = $(M.F⁺[1]) W/m^2\n")
    print(io, "    incoming radiation = $(M.F⁻[1]) W/m^2\n")
    print("  surface:\n")
    print(io, "    outgoing radiation = $(M.F⁺[end]) W/m^2\n")
    print(io, "    incoming radiation = $(M.F⁻[end]) W/m^2\n")
end

function FluxPack(np::Int, nν::Int, T::Type=Float64)
    FluxPack(
        zeros(T, np - 1, nν),
        zeros(T, np, nν),
        zeros(T, np, nν),
        zeros(T, np),
        zeros(T, np),
        zeros(T, np)
    )
end

function FluxPack(P::AbstractVector, ν::AbstractVector, T::Type=Float64)
    FluxPack(length(P), length(ν), T)
end

Base.size(F::FluxPack) = size(F.M⁺)

for op ∈ (:+, :-, :\, :*)
    @eval begin
        Base.$op(A::FluxPack, B::FluxPack) = FluxPack(
            broadcast($op, A.τ, B.τ),
            broadcast($op, A.M⁺, B.M⁺),
            broadcast($op, A.M⁻, B.M⁻),
            broadcast($op, A.F⁺, B.F⁺),
            broadcast($op, A.F⁻, B.F⁻),
            broadcast($op, A.Fnet, B.Fnet)
        )
    end
end

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

function ∫F!(F⁺, F⁻, M⁺, M⁻, ν)::Nothing
    n, m = size(M⁻)
    @assert size(M⁺) == size(M⁻)
    @assert length(ν) == m
    @assert length(F⁻) == length(F⁺) == n
    @inbounds for i ∈ 1:n
        Mᵢ⁺ = @view M⁺[i,:]
        Mᵢ⁻ = @view M⁻[i,:]
        F⁺[i] = trapz(ν, Mᵢ⁺)
        F⁻[i] = trapz(ν, Mᵢ⁻)
    end
    return nothing
end

ext: jl


text:
using Pkg
Pkg.Registry.add(RegistrySpec(url = "https://github.com/JuliaRegistries/General"))
Pkg.update()
Pkg.add(Pkg.PackageSpec(;name="Merly", version="1.0.2"))
using Merly

max_stars_repo_name: neomatrixcode/merly-app

max_issues_repo_name: neomatrixcode/merly-app

max_forks_repo_name: neomatrixcode/merly-app


text:
randn_tensor(dims, nonneg::Bool) = nonneg ? abs.(randn(dims)) : randn(dims)

rand_candecomp(r::Int64, dims::NTuple{N, Int};
               lambdas_nonneg::Bool=false, factors_nonneg::Bool=false) where N =
    CANDECOMP(ntuple(i -> randn_tensor((dims[i], r), factors_nonneg), N),
              randn_tensor(r, lambdas_nonneg))

rand_kruskal3(r::Int64, dims::NTuple{N, Int}, nonnegative::Bool) where N =
    compose(rand_candecomp(r, dims, lambdas_nonneg=nonnegative, factors_nonneg=nonnegative))

rand_tucker(core_dims::NTuple{N, Int}, dims::NTuple{N, Int};
            core_nonneg::Bool=false, factors_nonneg::Bool=false) where N =
    Tucker(ntuple(i -> randn_tensor((dims[i], core_dims[i]), factors_nonneg), N),
           randn_tensor(core_dims, core_nonneg))

function add_noise(tnsr::Array{T,N}, sn_ratio = 0.6, nonnegative::Bool = false) where {T <: Number, N}
    tnsr_noise = randn(size(tnsr)...)
    if nonnegative
        tnsr_noise .= max.(tnsr_noise, 0.0)
    end
    noise_scale = 10^(-sn_ratio/0.2)*norm(tnsr)/norm(tnsr_noise)
    return tnsr .+ noise_scale .* tnsr_noise
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions to render the numbers in text backend.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

function _render_number_text(
    number::Rational;
    compact::Bool = true
)
    if compact
        aux = abs(number.num)
        num_str = ""

        while aux ≥ 1
            i = aux % 10
            num_str = _EXPONENTS[i + 1] * num_str
            aux = floor(aux / 10)
        end

        aux = abs(number.den)
        den_str = ""

        while aux ≥ 1
            i = aux % 10
            den_str = _SUBNUMBERS[i + 1] * den_str
            aux = floor(aux / 10)
        end

        number_str = (number.num ≥ 0 ? "" : "- ") * num_str * "/" * den_str

        return number_str
    else
        num_str = sprint(print, abs(number.num))
        den_str = sprint(print, abs(number.den))

        num_width = textwidth(num_str)
        den_width = textwidth(den_str)

        fraction_width = max(num_width, den_width)

        # Align the numbers in the center.
        Δ = div(fraction_width - num_width, 2)
        num_str = " "^Δ * num_str

        Δ = div(fraction_width - den_width, 2)
        den_str = " "^Δ * den_str

        # Check if the fraction is negative. In this case, we must change the
        # alignemnt.
        is_neg = number.num < 0

        num_str  = (!is_neg ? "" : "  ") * num_str * "\n"
        num_str *= (!is_neg ? "" : "- ") * "—"^fraction_width * "\n"
        num_str *= (!is_neg ? "" : "  ") * den_str

        return num_str
    end
end

function _render_number_text(number::Irrational; kwargs...)
    return string(number)
end

function _render_number_text(
    number::Number;
    always_print_base::Bool = false,
    multiplication_sign::Char = '×',
    significand_format::String = "%g",
    show_base::Bool = true,
    show_significand::Bool = true,
    new_decimal_base::Union{Nothing, Number} = nothing
)
    # Get the significand and the base.
    if new_decimal_base !== nothing
        significand, base = _get_significand_and_base(number, new_decimal_base)
    else
        significand, base = _get_significand_and_base(number)
    end

    # Significand
    # ==========================================================================

    significand_str = ""

    if show_significand
        fmt = Printf.Format(significand_format)
        significand_str = Printf.format(fmt, significand)
    end

    # Base
    # ==========================================================================

    # If `base` is 0, then only show it if the user wants.
    base_str = ""

    if show_base && ((base != 0) || always_print_base)
        # Create the string representation.
        aux = abs(base)
        exponent_str = ""

        while aux ≥ 1
            i = aux % 10
            exponent_str = _EXPONENTS[i + 1] * exponent_str
            aux = floor(aux / 10)
        end

        # If `exponent_str` is empty, then the base is 0.
        if always_print_base
            exponent_str = "⁰"
        end

        if base < 0
            exponent_str = "⁻" * exponent_str
        end

        if show_significand
            base_str *= " " * multiplication_sign * " "
        end

        base_str *= "10" * exponent_str
    end

    # Output
    # ==========================================================================

    number_str = significand_str * base_str

    return number_str
end


text:
function Γ!(ˍ₋out, ˍ₋arg1; )
    let α = @inbounds(ˍ₋arg1[1]), β = @inbounds(ˍ₋arg1[2]), ρ = @inbounds(ˍ₋arg1[3]), δ = @inbounds(ˍ₋arg1[4]), σ = @inbounds(ˍ₋arg1[5]), Ω_1 = @inbounds(ˍ₋arg1[6])
        @inbounds begin
                ˍ₋out[1] = σ
                nothing
            end
    end
end

function Ω!(ˍ₋out, ˍ₋arg1; )
    let α = @inbounds(ˍ₋arg1[1]), β = @inbounds(ˍ₋arg1[2]), ρ = @inbounds(ˍ₋arg1[3]), δ = @inbounds(ˍ₋arg1[4]), σ = @inbounds(ˍ₋arg1[5]), Ω_1 = @inbounds(ˍ₋arg1[6])
        @inbounds begin
                ˍ₋out[1] = Ω_1
                ˍ₋out[2] = Ω_1
                nothing
            end
    end
end

function H̄!(ˍ₋out, ˍ₋arg1, ˍ₋arg2; )
    let c = @inbounds(ˍ₋arg1[1]), q = @inbounds(ˍ₋arg1[2]), i = @inbounds(ˍ₋arg1[3]), k = @inbounds(ˍ₋arg1[4]), z = @inbounds(ˍ₋arg1[5]), α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                ˍ₋out[1] = (+)((/)(1, c), (/)((*)(-1, β, (+)(1, (*)(-1, δ), (*)(α, (^)(k, (+)(-1, α)), (exp)(z)))), c))
                ˍ₋out[2] = (+)(c, k, (*)(-1, q), (*)(-1, k, (+)(1, (*)(-1, δ))))
                ˍ₋out[3] = (+)(q, (*)(-1, (^)(k, α), (exp)(z)))
                ˍ₋out[4] = (+)(z, (*)(-1, z, ρ))
                ˍ₋out[5] = (+)(i, (*)(-1, k), (*)(k, (+)(1, (*)(-1, δ))))
                nothing
            end
    end
end

function H̄_w!(ˍ₋out, ˍ₋arg1, ˍ₋arg2; )
    let c = @inbounds(ˍ₋arg1[1]), q = @inbounds(ˍ₋arg1[2]), i = @inbounds(ˍ₋arg1[3]), k = @inbounds(ˍ₋arg1[4]), z = @inbounds(ˍ₋arg1[5]), α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                ˍ₋out[1] = (+)((/)(-1, (^)(c, 2)), (/)((*)(β, (+)(1, (*)(-1, δ), (*)(α, (^)(k, (+)(-1, α)), (exp)(z)))), (^)(c, 2)))
                ˍ₋out[2] = 1
                ˍ₋out[7] = -1
                ˍ₋out[8] = 1
                ˍ₋out[15] = 1
                ˍ₋out[16] = (/)((*)(-1, α, β, (^)(k, (+)(-2, α)), (+)(-1, α), (exp)(z)), c)
                ˍ₋out[17] = δ
                ˍ₋out[18] = (*)(-1, α, (^)(k, (+)(-1, α)), (exp)(z))
                ˍ₋out[20] = (*)(-1, δ)
                ˍ₋out[21] = (/)((*)(-1, α, β, (^)(k, (+)(-1, α)), (exp)(z)), c)
                ˍ₋out[23] = (*)(-1, (^)(k, α), (exp)(z))
                ˍ₋out[24] = (+)(1, (*)(-1, ρ))
                nothing
            end
    end
end

function ȳ_iv!(ˍ₋out, ˍ₋arg1; )
    let α = @inbounds(ˍ₋arg1[1]), β = @inbounds(ˍ₋arg1[2]), ρ = @inbounds(ˍ₋arg1[3]), δ = @inbounds(ˍ₋arg1[4]), σ = @inbounds(ˍ₋arg1[5]), Ω_1 = @inbounds(ˍ₋arg1[6])
        @inbounds begin
                ˍ₋out[1] = (+)((^)((/)((+)(-1, δ, (/)(1, β)), α), (/)(α, (+)(-1, α))), (*)(-1, δ, (^)((/)((+)(-1, δ, (/)(1, β)), α), (/)(1, (+)(-1, α)))))
                ˍ₋out[2] = (^)((/)((+)(-1, δ, (/)(1, β)), α), (/)(α, (+)(-1, α)))
                ˍ₋out[3] = (*)(δ, (^)((/)((+)(-1, δ, (/)(1, β)), α), (/)(1, (+)(-1, α))))
                nothing
            end
    end
end

function x̄_iv!(ˍ₋out, ˍ₋arg1; )
    let α = @inbounds(ˍ₋arg1[1]), β = @inbounds(ˍ₋arg1[2]), ρ = @inbounds(ˍ₋arg1[3]), δ = @inbounds(ˍ₋arg1[4]), σ = @inbounds(ˍ₋arg1[5]), Ω_1 = @inbounds(ˍ₋arg1[6])
        @inbounds begin
                ˍ₋out[1] = (^)((/)((+)(-1, δ, (/)(1, β)), α), (/)(1, (+)(-1, α)))
                nothing
            end
    end
end

function ȳ!(ˍ₋out, ˍ₋arg1; )
    let α = @inbounds(ˍ₋arg1[1]), β = @inbounds(ˍ₋arg1[2]), ρ = @inbounds(ˍ₋arg1[3]), δ = @inbounds(ˍ₋arg1[4]), σ = @inbounds(ˍ₋arg1[5]), Ω_1 = @inbounds(ˍ₋arg1[6])
        @inbounds begin
                ˍ₋out[1] = (+)((^)((/)((+)(-1, δ, (/)(1, β)), α), (/)(α, (+)(-1, α))), (*)(-1, δ, (^)((/)((+)(-1, δ, (/)(1, β)), α), (/)(1, (+)(-1, α)))))
                ˍ₋out[2] = (^)((/)((+)(-1, δ, (/)(1, β)), α), (/)(α, (+)(-1, α)))
                ˍ₋out[3] = (*)(δ, (^)((/)((+)(-1, δ, (/)(1, β)), α), (/)(1, (+)(-1, α))))
                nothing
            end
    end
end

function x̄!(ˍ₋out, ˍ₋arg1; )
    let α = @inbounds(ˍ₋arg1[1]), β = @inbounds(ˍ₋arg1[2]), ρ = @inbounds(ˍ₋arg1[3]), δ = @inbounds(ˍ₋arg1[4]), σ = @inbounds(ˍ₋arg1[5]), Ω_1 = @inbounds(ˍ₋arg1[6])
        @inbounds begin
                ˍ₋out[1] = (^)((/)((+)(-1, δ, (/)(1, β)), α), (/)(1, (+)(-1, α)))
                nothing
            end
    end
end

function ȳ_p!(ˍ₋out, ::Val{:α}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                ˍ₋out[1] = (+)((/)((*)(-1, α, (/)((+)(-1, δ, (/)(1, β)), (^)(α, 2)), (^)((/)((+)(-1, δ, (/)(1, β)), α), (+)(-1, (/)(α, (+)(-1, α))))), (+)(-1, α)), (*)(-1, δ, (+)((/)((*)(-1, (/)((+)(-1, δ, (/)(1, β)), (^)(α, 2)), (^)((/)((+)(-1, δ, (/)(1, β)), α), (+)(-1, (/)(1, (+)(-1, α))))), (+)(-1, α)), (*)(-1, (/)(1, (^)((+)(-1, α), 2)), (^)((/)((+)(-1, δ, (/)(1, β)), α), (/)(1, (+)(-1, α))), (log)((/)((+)(-1, δ, (/)(1, β)), α))))), (*)((+)((/)(1, (+)(-1, α)), (/)((*)(-1, α), (^)((+)(-1, α), 2))), (^)((/)((+)(-1, δ, (/)(1, β)), α), (/)(α, (+)(-1, α))), (log)((/)((+)(-1, δ, (/)(1, β)), α))))
                ˍ₋out[2] = (+)((/)((*)(-1, α, (/)((+)(-1, δ, (/)(1, β)), (^)(α, 2)), (^)((/)((+)(-1, δ, (/)(1, β)), α), (+)(-1, (/)(α, (+)(-1, α))))), (+)(-1, α)), (*)((+)((/)(1, (+)(-1, α)), (/)((*)(-1, α), (^)((+)(-1, α), 2))), (^)((/)((+)(-1, δ, (/)(1, β)), α), (/)(α, (+)(-1, α))), (log)((/)((+)(-1, δ, (/)(1, β)), α))))
                ˍ₋out[3] = (*)(δ, (+)((/)((*)(-1, (/)((+)(-1, δ, (/)(1, β)), (^)(α, 2)), (^)((/)((+)(-1, δ, (/)(1, β)), α), (+)(-1, (/)(1, (+)(-1, α))))), (+)(-1, α)), (*)(-1, (/)(1, (^)((+)(-1, α), 2)), (^)((/)((+)(-1, δ, (/)(1, β)), α), (/)(1, (+)(-1, α))), (log)((/)((+)(-1, δ, (/)(1, β)), α)))))
                nothing
            end
    end
end

function ȳ_p!(ˍ₋out, ::Val{:ρ}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                nothing
            end
    end
end

function ȳ_p!(ˍ₋out, ::Val{:σ}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                nothing
            end
    end
end

function ȳ_p!(ˍ₋out, ::Val{:Ω_1}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                nothing
            end
    end
end

function ȳ_p!(ˍ₋out, ::Val{:δ}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                ˍ₋out[1] = (+)((/)((*)(α, (^)((/)((+)(-1, δ, (/)(1, β)), α), (+)(-1, (/)(α, (+)(-1, α))))), (*)(α, (+)(-1, α))), (*)(-1, (^)((/)((+)(-1, δ, (/)(1, β)), α), (/)(1, (+)(-1, α)))), (/)((*)(-1, δ, (^)((/)((+)(-1, δ, (/)(1, β)), α), (+)(-1, (/)(1, (+)(-1, α))))), (*)(α, (+)(-1, α))))
                ˍ₋out[2] = (/)((*)(α, (^)((/)((+)(-1, δ, (/)(1, β)), α), (+)(-1, (/)(α, (+)(-1, α))))), (*)(α, (+)(-1, α)))
                ˍ₋out[3] = (+)((/)((*)(δ, (^)((/)((+)(-1, δ, (/)(1, β)), α), (+)(-1, (/)(1, (+)(-1, α))))), (*)(α, (+)(-1, α))), (^)((/)((+)(-1, δ, (/)(1, β)), α), (/)(1, (+)(-1, α))))
                nothing
            end
    end
end

function ȳ_p!(ˍ₋out, ::Val{:β}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                ˍ₋out[1] = (+)((/)((*)(δ, (/)(1, (^)(β, 2)), (^)((/)((+)(-1, δ, (/)(1, β)), α), (+)(-1, (/)(1, (+)(-1, α))))), (*)(α, (+)(-1, α))), (/)((*)(-1, α, (/)(1, (^)(β, 2)), (^)((/)((+)(-1, δ, (/)(1, β)), α), (+)(-1, (/)(α, (+)(-1, α))))), (*)(α, (+)(-1, α))))
                ˍ₋out[2] = (/)((*)(-1, α, (/)(1, (^)(β, 2)), (^)((/)((+)(-1, δ, (/)(1, β)), α), (+)(-1, (/)(α, (+)(-1, α))))), (*)(α, (+)(-1, α)))
                ˍ₋out[3] = (/)((*)(-1, δ, (/)(1, (^)(β, 2)), (^)((/)((+)(-1, δ, (/)(1, β)), α), (+)(-1, (/)(1, (+)(-1, α))))), (*)(α, (+)(-1, α)))
                nothing
            end
    end
end

function x̄_p!(ˍ₋out, ::Val{:α}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                ˍ₋out[1] = (+)((/)((*)(-1, (/)((+)(-1, δ, (/)(1, β)), (^)(α, 2)), (^)((/)((+)(-1, δ, (/)(1, β)), α), (+)(-1, (/)(1, (+)(-1, α))))), (+)(-1, α)), (*)(-1, (/)(1, (^)((+)(-1, α), 2)), (^)((/)((+)(-1, δ, (/)(1, β)), α), (/)(1, (+)(-1, α))), (log)((/)((+)(-1, δ, (/)(1, β)), α))))
                nothing
            end
    end
end

function x̄_p!(ˍ₋out, ::Val{:ρ}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                nothing
            end
    end
end

function x̄_p!(ˍ₋out, ::Val{:σ}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                nothing
            end
    end
end

function x̄_p!(ˍ₋out, ::Val{:Ω_1}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                nothing
            end
    end
end

function x̄_p!(ˍ₋out, ::Val{:δ}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                ˍ₋out[1] = (/)((^)((/)((+)(-1, δ, (/)(1, β)), α), (+)(-1, (/)(1, (+)(-1, α)))), (*)(α, (+)(-1, α)))
                nothing
            end
    end
end

function x̄_p!(ˍ₋out, ::Val{:β}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                ˍ₋out[1] = (/)((*)(-1, (/)(1, (^)(β, 2)), (^)((/)((+)(-1, δ, (/)(1, β)), α), (+)(-1, (/)(1, (+)(-1, α))))), (*)(α, (+)(-1, α)))
                nothing
            end
    end
end

const steady_state! = nothing

function Γ_p!(ˍ₋out, ::Val{:α}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                nothing
            end
    end
end

function Γ_p!(ˍ₋out, ::Val{:ρ}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                nothing
            end
    end
end

function Γ_p!(ˍ₋out, ::Val{:σ}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                ˍ₋out[1] = 1
                nothing
            end
    end
end

function Γ_p!(ˍ₋out, ::Val{:Ω_1}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                nothing
            end
    end
end

function Γ_p!(ˍ₋out, ::Val{:δ}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                nothing
            end
    end
end

function Γ_p!(ˍ₋out, ::Val{:β}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                nothing
            end
    end
end

function Ω_p!(ˍ₋out, ::Val{:α}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                nothing
            end
    end
end

function Ω_p!(ˍ₋out, ::Val{:ρ}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                nothing
            end
    end
end

function Ω_p!(ˍ₋out, ::Val{:σ}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                nothing
            end
    end
end

function Ω_p!(ˍ₋out, ::Val{:Ω_1}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                ˍ₋out[1] = 1
                ˍ₋out[2] = 1
                nothing
            end
    end
end

function Ω_p!(ˍ₋out, ::Val{:δ}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                nothing
            end
    end
end

function Ω_p!(ˍ₋out, ::Val{:β}, ˍ₋arg2; )
    let α = @inbounds(ˍ₋arg2[1]), β = @inbounds(ˍ₋arg2[2]), ρ = @inbounds(ˍ₋arg2[3]), δ = @inbounds(ˍ₋arg2[4]), σ = @inbounds(ˍ₋arg2[5]), Ω_1 = @inbounds(ˍ₋arg2[6])
        @inbounds begin
                nothing
            end
    end
end

max_stars_repo_name: HighDimensionalEconLab/HMCExamples

max_issues_repo_name: HighDimensionalEconLab/HMCExamples

max_forks_repo_name: HighDimensionalEconLab/HMCExamples


text:
using SolitonDistribution
using StatsBase, Statistics, Distributions
using Test

@testset "SolitonDistribution.jl" begin
    K, M, δ, atol = 100, 60, 0.2, 0
    Ω = Soliton(K, M, δ, atol)
    @test pdf(Ω, M) > pdf(Ω, M-1)
    @test pdf(Ω, M) > pdf(Ω, M+1)
    @test cumsum(pdf.(Ω, 1:K)) ≈ cdf.(Ω, 1:K)
    @test cdf(Ω, 0) ≈ 0
    @test cdf(Ω, K) ≈ 1
    @test mean(Ω) ≈ 7.448826535558562
    @test var(Ω) ≈ 178.91717915136957
    @test minimum(Ω) == 1
    @test maximum(Ω) == K
    @test quantile(Ω, 0) == 1
    @test quantile(Ω, 1) == K
    @test insupport(Ω, 1) && insupport(Ω, 2) && insupport(Ω, K)
    @test !insupport(Ω, 0) && !insupport(Ω, 2.1) && !insupport(Ω, K + 1)

    K, M, δ, atol = 100, 60, 0.2, 1e-3
    Ω = Soliton(K, M, δ, atol)
    ds = [d for d in 1:K if pdf(Ω, d) > 0]
    @test all(pdf.(Ω, ds) .> atol)
    @test cdf(Ω, 0) ≈ 0
    @test cdf(Ω, K) ≈ 1
    @test minimum(Ω) == 1
    @test maximum(Ω) == K
    @test quantile(Ω, 0) == 1
    @test quantile(Ω, 1) == M
end

max_stars_repo_name: severinson/SolitonDistribution

max_issues_repo_name: severinson/SolitonDistribution

max_forks_repo_name: severinson/SolitonDistribution


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct FlowLogInformation <: SwaggerModel
    targetResourceId::Any # spec type: Union{ Nothing, String } # spec name: targetResourceId
    properties::Any # spec type: Union{ Nothing, FlowLogProperties } # spec name: properties

    function FlowLogInformation(;targetResourceId=nothing, properties=nothing)
        o = new()
        validate_property(FlowLogInformation, Symbol("targetResourceId"), targetResourceId)
        setfield!(o, Symbol("targetResourceId"), targetResourceId)
        validate_property(FlowLogInformation, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type FlowLogInformation

const _property_map_FlowLogInformation = Dict{Symbol,Symbol}(Symbol("targetResourceId")=>Symbol("targetResourceId"), Symbol("properties")=>Symbol("properties"))
const _property_types_FlowLogInformation = Dict{Symbol,String}(Symbol("targetResourceId")=>"String", Symbol("properties")=>"FlowLogProperties")
Base.propertynames(::Type{ FlowLogInformation }) = collect(keys(_property_map_FlowLogInformation))
Swagger.property_type(::Type{ FlowLogInformation }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_FlowLogInformation[name]))}
Swagger.field_name(::Type{ FlowLogInformation }, property_name::Symbol) =  _property_map_FlowLogInformation[property_name]

function check_required(o::FlowLogInformation)
    (getproperty(o, Symbol("targetResourceId")) === nothing) && (return false)
    (getproperty(o, Symbol("properties")) === nothing) && (return false)
    true
end

function validate_property(::Type{ FlowLogInformation }, name::Symbol, val)
end


text:
"""An Analytic1DFilter has only positive frequencies. Its Fourier-domain support
ranges between posfirst and (posfirst+length(pos)-1)<N/2."""
immutable Analytic1DFilter{T} <: AbstractFilter{T,FourierDomain{1}}
    pos::Vector{T}
    posfirst::Int
end

"""A Coanalytic1DFilter has only negative frequencies. Its Fourier-domain
support ranges between (neglast-length(neg)+1)>(-N/2) and neglast."""
immutable Coanalytic1DFilter{T} <: AbstractFilter{T,FourierDomain{1}}
    neg::Vector{T}
    neglast::Int
end

"""A FourierSymmetric1DFilter has a real impulse respnse in the spatial domain,
hence a symmetric response in the Fourier domain. It consists of a leg of
nonzero frequencies and a ""zero"" coefficient."""
immutable FourierSymmetric1DFilter{T} <: AbstractFilter{T,FourierDomain{1}}
    leg::Vector{T}
    zero::T
end

"""A FullResolution1DFilter spans the entire spectrum. Its is defined in the
Fourier domain as a vector of size N."""
immutable FullResolution1DFilter{T} <: AbstractFilter{T,FourierDomain{1}}
    coeff::Vector{T}
end

"""A Vanishing1DFilter has a limited Fourier-domain support, yet spanning both
in positive and negative frequencies. Moreover, it is null for ω=0 and
at the midpoint ω=N/2. Thus, it is defined as the combination of an
Analytic1DFilter (positive frequencies) and a Coanalytic1DFilter (negative
frequencies)."""
immutable Vanishing1DFilter{T} <: AbstractFilter{T,FourierDomain{1}}
    an::Analytic1DFilter
    coan::Coanalytic1DFilter{T}
end

"""A VanishingWithMidpoint1DFilter has a limited Fourier-domain support, yet
spanning both in positive and negative frequencies. It is null for ω=0 but
nonzero at the midpoint ω=N/2. Thus, it is defined as the combination of
an Analytic1DFilter (positive frequencies), a Coanalytic1DFilter (negative
frequencies), and a midpoint (frequency ω=N/2)."""
immutable VanishingWithMidpoint1DFilter{T} <: AbstractFilter{T,FourierDomain{1}}
    an::Analytic1DFilter{T}
    coan::Coanalytic1DFilter{T}
    midpoint::T
end

"""Given an array `y` and a filter bank specification `spec`, the
constructor `AbstractFilter` truncates `y` according to the magnitude
threshold `spec.ɛ` and returns the appropriate subtype of `AbstractFilter`
to represent the truncated array. The output is cast with the supertype
`AbstractFilter` so that the output of this function can be stored
efficiently into an array of abstract objects.
For one-dimensional filters in the Fourier domain, the output type may be:
* `Analytic1DFilter` if `y` is non-negligible only for positive frequencies
* `Coanalytic1DFilter` if `y` is non-negligible only for negative frequencies
* `Vanishing1DFilter` if the support of `y` encompasses a portion of
  positive frequencies and a portion of negative frequencies, excluding
  midpoint.
* `VanishingWithMidpoint1DFilter` if the support of `y` encompasses a portion
  of positive frequencies and a portion of negative frequencies, including
  midpoint.
* `FullResolution1DFilter` if `y` is non-negligible along the whole spectrum"""
function AbstractFilter{T}(y::Vector{T}, spec::AbstractSpec{T,FourierDomain{1}})
    supertype = AbstractFilter{T,FourierDomain{1}}
    N = 1 << spec.log2_size
    halfN = N >> 1
    ɛ2 = T(spec.ɛ * spec.ɛ)
    y2 = abs2.(y)
    negbools = y2[2:halfN] .> ɛ2
    negfirst, neglast = findfirst(negbools), findlast(negbools)
    posbools = y2[(2+halfN):end] .> ɛ2
    posfirst, poslast = findfirst(posbools), findlast(posbools)
    hasmidpoint = y2[1] .> ɛ2
    if hasmidpoint
        (negfirst == 1) && (neglast == (halfN-1)) &&
            (posfirst == 1) && (poslast == (halfN-1)) &&
            return FullResolution1DFilter(fftshift(y))::supertype
        midpoint = y[1]
        if (neglast == 0)
            coan = Coanalytic1DFilter(zeros(T, 1), -halfN + 1)
        else
            coan = Coanalytic1DFilter(y[(1+negfirst):(1+neglast)],
                neglast - halfN)
        end
        if (poslast == 0)
            an = Analytic1DFilter(zeros(T, 1), halfN - 1)
        else
            an = Analytic1DFilter(y[(1+halfN+posfirst):(1+halfN+poslast)],
                posfirst)
        end
        return VanishingWithMidpoint1DFilter(an, coan, midpoint)::supertype
    end
    if (negfirst == 0)
        return Analytic1DFilter(y[(1+halfN+posfirst):(1+halfN+poslast)],
            posfirst)::supertype
    end
    if (posfirst == 0)
        return Coanalytic1DFilter(y[(1+negfirst):(1+neglast)],
            neglast)::supertype
    end
    an = Analytic1DFilter(y[(1+halfN+posfirst):(1+halfN+poslast)], posfirst)
    coan = Coanalytic1DFilter(y[(1+negfirst):(1+neglast)], neglast - halfN)
    return Vanishing1DFilter(an, coan)::supertype
end

# multiplication operator with scalar *
Base.:*{T}(ψ::Analytic1DFilter{T}, b::Number) =
    Analytic1DFilter{T}(b * ψ.pos, ψ.posfirst)
Base.:*{T}(ψ::Coanalytic1DFilter{T}, b::Number) =
    Coanalytic1DFilter{T}(b * ψ.neg, ψ.neglast)
Base.:*{T}(ψ::FullResolution1DFilter{T}, b::Number) =
    FullResolution1DFilter{T}(b * ψ.coeff)
Base.:*{T}(ψ::FourierSymmetric1DFilter{T}, b::Number) =
    FourierSymmetric1DFilter{T}(b * ψ.leg, T(b) * ψ.zero)
Base.:*{T}(ψ::Vanishing1DFilter{T}, b::Number) =
    Vanishing1DFilter{T}(b * ψ.an, b * ψ.coan)
Base.:*{T}(ψ::VanishingWithMidpoint1DFilter{T}, b::Number) =
    VanishingWithMidpoint1DFilter{T}(b * ψ.an, b * ψ.coan, T(b) * ψ.midpoint)

"""Given a narrowband filter `ψ` and its corresponding filter bank
specification `spec`, returns the base-2 logarithm of the resampling factor
such that the passband of `ψ` is not aliased. This resmpling factor is
expressed in fraction of signal length ; thus, its base-2 logarithm is
non-positive. The result does not account for manual oversampling."""
critical_log2_sampling(ψ::Analytic1DFilter, spec::AbstractSpec) =
    1 + nextpow2_exponent(ψ.posfirst + length(ψ.pos) - 1) - spec.log2_size
critical_log2_sampling(ψ::Coanalytic1DFilter, spec::AbstractSpec) =
    1 + nextpow2_exponent(ψ.neglast - length(ψ.neg) + 1) - spec.log2_size
critical_log2_sampling(ψ::FullResolution1DFilter, spec::AbstractSpec) = 0
critical_log2_sampling(ψ::FourierSymmetric1DFilter, spec::AbstractSpec) =
    1 + nextpow2_exponent(length(ψ.leg)) - spec.log2_size
critical_log2_sampling(ψ::Vanishing1DFilter, spec::AbstractSpec) = max(
    critical_log2_sampling(ψ.an, spec), critical_log2_sampling(ψ.coan, spec))
critical_log2_sampling(ψ::VanishingWithMidpoint1DFilter, spec::AbstractSpec) = 0

# indexing between -N/2 and N/2-1
function Base.getindex{T}(ψ::Analytic1DFilter{T}, i::Integer)
    i < ψ.posfirst && return zero(T)
    i > (ψ.posfirst + length(ψ.pos) - 1) && return zero(T)
    @inbounds return ψ.pos[1 - ψ.posfirst + i]
end
function Base.getindex{T}(ψ::Analytic1DFilter{T}, I::UnitRange{Int64})
    I.stop < ψ.posfirst && return zeros(T, length(I))
    start = max(I.start, ψ.posfirst)
    stop = min(I.stop, ψ.posfirst + length(ψ.pos) - 1)
    return T[
        zeros(T, max(start - I.start, 0)) ;
        ψ.pos[1 - ψ.posfirst + (start:stop)] ;
        zeros(T, max(I.stop - max(I.start - 1, stop), 0)) ]
end
function Base.getindex{T}(ψ::Coanalytic1DFilter{T}, i::Integer)
    i<(ψ.neglast - length(ψ.neg) + 1) && return zero(T)
    i>ψ.neglast && return zero(T)
    return ψ.neg[i - ψ.neglast + end]
end
function Base.getindex{T}(ψ::Coanalytic1DFilter{T}, I::UnitRange{Int64})
    I.start > ψ.neglast && return zeros(T, length(I))
    start = max(I.start, ψ.neglast - length(ψ.neg) + 1)
    stop = min(I.stop, ψ.neglast)
    return T[
        zeros(T, max(min(start, I.stop + 1) - I.start, 0));
        ψ.neg[end - ψ.neglast + (start:stop)];
        zeros(T, max(I.stop - max(I.start - 1, stop), 0)) ]
end
function Base.getindex{T}(ϕ::FourierSymmetric1DFilter{T}, i::Integer)
    i==0 && return ϕ.zero
    abs(i)>length(ϕ.leg) && return zero(T)
    return ϕ.leg[abs(i)]
end
function Base.getindex{T}(ϕ::FourierSymmetric1DFilter{T}, I::UnitRange{Int64})
    T[ ϕ[ω] for ω in I ]
end
function Base.getindex{T}(ψ::FullResolution1DFilter{T}, i::Integer)
    N = length(ψ.coeff)
    halfN = N >> 1
    if i < -halfN
        return zero(T)
    elseif i < 0
        return ψ.coeff[1+i+N]
    elseif i < halfN
        return ψ.coeff[1+i]
    else
        return zero(T)
    end
end
function Base.getindex{T}(ψ::FullResolution1DFilter{T}, I::UnitRange{Int64})
    N = length(ψ.coeff)
    halfN = N >> 1
    I.start > halfN && return zeros(T, length(I))
    I.stop < -halfN && return zeros(T, length(I))
    return T[
        zeros(T, max(-halfN-I.start, 0));
        ψ.coeff[1 + N + (max(I.start, -halfN):min(I.stop, -1))];
        ψ.coeff[1 + (max(I.start, 0):min(I.stop, halfN-1))];
        zeros(T, max(I.stop-halfN+1, 0))
    ]
end
function Base.getindex{T}(ψ::Vanishing1DFilter{T}, i::Integer)
    return (i > 0 ? ψ.an[i] : ψ.coan[i])
end
function Base.getindex{T}(ψ::Vanishing1DFilter{T}, I::UnitRange{Int64})
    return T[
        ψ.coan[min(0, I.start):min(-1, I.stop)] ;
        zeros(T, Int((I.start < 0) && (I.stop > 0))) ;
        ψ.an[max(1, I.start):max(0, I.stop)] ]
end
function Base.getindex(ψ::VanishingWithMidpoint1DFilter, i::Integer)
    halfN = ψ.an.posfirst + length(ψ.an.pos)
    i==(-halfN) && return ψ.midpoint
    return (i > 0 ? ψ.an[i] : ψ.coan[i])
end
function Base.getindex{T}(ψ::VanishingWithMidpoint1DFilter{T},
                          I::UnitRange{Int64})
    halfN = ψ.an.posfirst + length(ψ.an.pos)
    output = T[
        ψ.coan[min(0, I.start):min(-1, I.stop)] ;
        zeros(T, Int((I.start < 0) && (I.stop > 0))) ;
        ψ.an[max(1, I.start):max(0, I.stop)] ]
    if I.start <= -halfN
        output[-halfN-I.start+1] = ψ.midpoint
    end
    return output
end

"""Adds the squared magnitude of a Fourier-domain wavelet `ψ` to an accumulator
`lp` (squared Littlewood-Paley sum)."""
function littlewoodpaleyadd!(lp::Vector, ψ::Analytic1DFilter)
    @inbounds for ω in eachindex(ψ.pos)
        @fastmath lp[ψ.posfirst + ω] += abs2(ψ.pos[ω])
    end
end
function littlewoodpaleyadd!(lp::Vector, ψ::Coanalytic1DFilter)
    @inbounds for ω in eachindex(ψ.neg)
        @fastmath lp[1 + length(lp) + ψ.neglast - length(ψ.neg) + ω] +=
            abs2(ψ.neg[ω])
    end
end
function littlewoodpaleyadd!(lp::Vector, ψ::FullResolution1DFilter)
    @inbounds for ω in eachindex(ψ.coeff)
        @fastmath lp[ω] += abs2(ψ.coeff[ω])
    end
end
function littlewoodpaleyadd!(lp::Vector, ψ::FourierSymmetric1DFilter)
    @inbounds for ω in eachindex(ψ.leg)
        @fastmath lp[1 + ω] += abs2(ψ.leg[ω])
        @fastmath lp[length(lp) + 1 - ω] += abs2(ψ.leg[ω])
    end
    @fastmath lp[1 + 0] += abs2(ψ.zero)
end
function littlewoodpaleyadd!(lp::Vector, ψ::Vanishing1DFilter)
    littlewoodpaleyadd!(lp, ψ.an)
    littlewoodpaleyadd!(lp, ψ.coan)
end
function littlewoodpaleyadd!(lp::Vector, ψ::VanishingWithMidpoint1DFilter)
    littlewoodpaleyadd!(lp, ψ.an)
    littlewoodpaleyadd!(lp, ψ.coan)
    @fastmath lp[1 + (length(lp)>>1)] += abs2(ψ.midpoint)
end

"""Returns the maximum Fourier-domain absolute value of a filter."""
Base.maximum(ψ::Analytic1DFilter) = sqrt(maximum(abs2.(ψ.pos)))
Base.maximum(ψ::Coanalytic1DFilter) = sqrt(maximum(abs2.(ψ.neg)))
Base.maximum(ψ::FullResolution1DFilter) = sqrt(maximum(abs2.(ψ.coeff)))
Base.maximum(ψ::FourierSymmetric1DFilter) =
    sqrt(max(maximum(abs2.(ψ.leg)), abs2(ψ.zero)))
Base.maximum(ψ::Vanishing1DFilter) = max(maximum(ψ.an), maximum(ψ.coan))
Base.maximum(ψ::VanishingWithMidpoint1DFilter) =
    max(maximum(ψ.an), maximum(ψ.coan), abs(ψ.midpoint))

"""The exponent of the smallest power of two not less than input integer `n`.
Returns `0` for `n==0`, and returns `-nextpow2_exponent(-n)` for negative
arguments."""
nextpow2_exponent(n::Unsigned) = (sizeof(n)<<3)-leading_zeros(n-1)
function nextpow2_exponent(n::Integer)
    if n == 0
        return oftype(n, 0)
    elseif n < 0
        return oftype(n, -nextpow2_exponent(unsigned(-n)))
    else
        return oftype(n, nextpow2_exponent(unsigned(n)))
    end
end

"""Renormalizes an array of Fourier-domain wavelets `ψs` by the square root
of the maximum value of its Littlewood-Paley sum `lp`. This ensures approximate
an energy conservation property for the subsequent wavelet transform operator.
The Littlewood-Paley sum `lp` is defined, for each frequency `ω`, as the sum of
wavelet energies (squared magnitudes)."""
function renormalize!{T<:Number,G<:LineGroups}(
        ϕ::FourierSymmetric1DFilter{T},
        ψs::Array{AbstractFilter{T,FourierDomain{1}},3},
        spec::AbstractSpec{T,FourierDomain{1},G})
    lp = zeros(real(T), 1<<spec.log2_size)
    littlewoodpaleyadd!(lp, ϕ)
    for idψ in eachindex(ψs[1, :, :]) littlewoodpaleyadd!(lp, ψs[idψ]); end
    get_n_orientations(spec.pointgroup)==1 && symmetrize!(lp)
    ψs .*= inv(sqrt(lp[1+end>>1]))
    return scale!(lp, inv(lp[1+end>>1]))
end

"""Reverses the coefficients of a Fourier-domain 1D filter `ψ` to yield a
""mirror filter"" whose center frequency is of opposite sign. This corresponds
to a conjugation in the time domain."""
spin(ψ::Analytic1DFilter) = Coanalytic1DFilter(reverse(ψ.pos), -ψ.posfirst)
spin(ψ::Coanalytic1DFilter) = Analytic1DFilter(reverse(ψ.neg), -ψ.neglast)
spin(ψ::FullResolution1DFilter) = FullResolution1DFilter(reverse(ψ.coeff))
spin(ψ::Vanishing1DFilter) = Vanishing1DFilter(spin(ψ.coan), spin(ψ.an))
spin(ψ::VanishingWithMidpoint1DFilter) =
    VanishingWithMidpoint1DFilter(spin(ψ.coan), spin(ψ.an), ψ.midpoint)

"""Returns the symmetric part of a input vector. This function is called by
`renormalize!` to symmetrize the Littlewood-Paley sum of a one-dimensional
filter bank."""
function symmetrize!(lp::Vector)
    N = length(lp)
    for ω in 1:(N>>1 - 1)
        halfsum = 0.5 * (lp[1 + ω] + lp[1 + N - ω])
        lp[1 + ω] = halfsum
        lp[1 + N - ω] = halfsum
    end
end


text:
@testset "linalg" begin
    @testset "dot" begin
        @testset "Vector" begin
            M = 3
            x, y = randn(M), randn(M)
            ẋ, ẏ = randn(M), randn(M)
            x̄, ȳ = randn(M), randn(M)
            frule_test(dot, (x, ẋ), (y, ẏ))
            rrule_test(dot, randn(), (x, x̄), (y, ȳ))
        end
        @testset "Matrix" begin
            M, N = 3, 4
            x, y = randn(M, N), randn(M, N)
            ẋ, ẏ = randn(M, N), randn(M, N)
            x̄, ȳ = randn(M, N), randn(M, N)
            frule_test(dot, (x, ẋ), (y, ẏ))
            rrule_test(dot, randn(), (x, x̄), (y, ȳ))
        end
        @testset "Array{T, 3}" begin
            M, N, P = 3, 4, 5
            x, y = randn(M, N, P), randn(M, N, P)
            ẋ, ẏ = randn(M, N, P), randn(M, N, P)
            x̄, ȳ = randn(M, N, P), randn(M, N, P)
            frule_test(dot, (x, ẋ), (y, ẏ))
            rrule_test(dot, randn(), (x, x̄), (y, ȳ))
        end
    end
    @testset "inv" begin
        N = 3
        B = generate_well_conditioned_matrix(N)
        frule_test(inv, (B, randn(N, N)))
        rrule_test(inv, randn(N, N), (B, randn(N, N)))
    end
    @testset "det" begin
        N = 3
        B = generate_well_conditioned_matrix(N)
        frule_test(det, (B, randn(N, N)))
        rrule_test(det, randn(), (B, randn(N, N)))
    end
    @testset "logdet" begin
        N = 3
        B = generate_well_conditioned_matrix(N)
        frule_test(logdet, (B, randn(N, N)))
        rrule_test(logdet, randn(), (B, randn(N, N)))
    end
    @testset "tr" begin
        N = 4
        frule_test(tr, (randn(N, N), randn(N, N)))
        rrule_test(tr, randn(), (randn(N, N), randn(N, N)))
    end
    @testset "*" begin
        dims = [3,4,5]
        for n in dims, m in dims, p in dims
            n > 3 && n == m == p && continue  # don't need to test square case multiple times
            A = randn(m, n)
            B = randn(n, p)
            Ȳ = randn(m, p)
            rrule_test(*, Ȳ, (A, randn(m, n)), (B, randn(n, p)))
        end
    end
    @testset "$f" for f in [/, \]
        @testset "Matrix" begin
            for n in 3:5, m in 3:5
                A = randn(m, n)
                B = randn(m, n)
                Ȳ = randn(size(f(A, B)))
                rrule_test(f, Ȳ, (A, randn(m, n)), (B, randn(m, n)))
            end
        end
        @testset "Vector" begin
            x = randn(10)
            y = randn(10)
            ȳ = randn(size(f(x, y))...)
            rrule_test(f, ȳ, (x, randn(10)), (y, randn(10)))
        end
        if f == (/)
            @testset "$T on the RHS" for T in (Diagonal, UpperTriangular, LowerTriangular)
                RHS = T(randn(T == Diagonal ? 10 : (10, 10)))
                Y = randn(5, 10)
                Ȳ = randn(size(f(Y, RHS))...)
                rrule_test(f, Ȳ, (Y, randn(size(Y))), (RHS, randn(size(RHS))))
            end
        else
            @testset "$T on LHS" for T in (Diagonal, UpperTriangular, LowerTriangular)
                LHS = T(randn(T == Diagonal ? 10 : (10, 10)))
                y = randn(10)
                ȳ = randn(size(f(LHS, y))...)
                rrule_test(f, ȳ, (LHS, randn(size(LHS))), (y, randn(10)))
                Y = randn(10, 10)
                Ȳ = randn(10, 10)
                rrule_test(f, Ȳ, (LHS, randn(size(LHS))), (Y, randn(size(Y))))
            end
            @testset "Matrix $f Vector" begin
                X = randn(10, 4)
                y = randn(10)
                ȳ = randn(size(f(X, y))...)
                rrule_test(f, ȳ, (X, randn(size(X))), (y, randn(10)))
            end
            @testset "Vector $f Matrix" begin
                x = randn(10)
                Y = randn(10, 4)
                ȳ = randn(size(f(x, Y))...)
                rrule_test(f, ȳ, (x, randn(size(x))), (Y, randn(size(Y))))
            end
        end
    end
    @testset "norm" begin
        for dims in [(), (5,), (3, 2), (7, 3, 2)]
            A = randn(dims...)
            p = randn()
            ȳ = randn()
            rrule_test(norm, ȳ, (A, randn(dims...)), (p, randn()))
        end
    end
end


text:
using SparseArrays: sparse
using MatrixNetworks: scomponents

export adjacency
"""
    adjacency(cells::AbstractVector{Cell{N}}) where N

Find an adjacency matrix for the graph with `cells` as vertices and such that two cells are
considered adjacent if they share a child. Additionally return a dictionary mapping pairs
of adjacent cells to their shared child.
"""
function adjacency(cells::AbstractVector{Cell{N}}) where N
    function common_face_index(c1, c2)
        for (i, c) in enumerate(c1.children)
            if c2 in keys(c.parents)
                return i
            end
        end
        return 0
    end
    faces = Dict{Tuple{Cell{N}, Cell{N}}, Cell{N}}()
    row_inds, col_inds, vals = Int[], Int[], Int[]
    for (i1, i2) in combinations(1:length(cells), 2)
        c1, c2 = cells[i1], cells[i2]
        i = common_face_index(c1, c2)
        if i > 0
            push!(row_inds, i1); push!(col_inds, i2); push!(vals, 1)
            push!(row_inds, i2); push!(col_inds, i1); push!(vals, 1)
            faces[(c1, c2)] = faces[(c2, c1)] = c1.children[i]
        end
    end
    return sparse(row_inds, col_inds, vals, length(cells), length(cells)), faces
end

export connected_components
"""
    connected_components(cells::AbstractVector{Cell{N}}) where N

Find the connected components of the given cells using the adjacency matrix given by
`adjacency`.
"""
function connected_components(cells::AbstractVector{Cell{N}}) where N
    adj, _ = adjacency(cells)
    cc = scomponents(adj)
    components = [Cell{N}[] for i in 1:cc.number]
    for (i, c) in enumerate(cells)
        push!(components[cc.map[i]], c)
    end
    return components
end

"""
    connected_components(comp::CellComplex)

Find the connected components of the highest dimensional cells of a cell complex.
"""
connected_components(comp::CellComplex) = connected_components(comp.cells[end])

export boundary
"""
    boundary(cells::AbstractVector{Cell{N}}) where N

Find the boundary cells, defined as those child cells with exactly one parent.
"""
function boundary(cells::AbstractVector{Cell{N}}) where N
    children = unique(vcat([c.children for c in cells]...))
    return filter(c -> length(c.parents) == 1, children)
end

"""
    boundary(comp::CellComplex)

Find the boundary of the highest dimensional cells of `comp` and create a cell complex with
the boundary and all of its descendants.
"""
boundary(comp::CellComplex) = CellComplex(boundary(comp.cells[end]))

export boundary_components_connected
"""
    boundary_components_connected(cells::AbstractVector{Cell{N}}) where N

Find the connected components of the boundary assuming the cells are connected. Return the
interior boundaries and the exterior boundary.
"""
function boundary_components_connected(cells::AbstractVector{Cell{N}}) where N
    boundary_comps = connected_components(boundary(cells))
    # identify the exterior boundary by finding which component contains
    # a point with maximal first coordinate
    exterior_ind = argmax([maximum([maximum(p.coords[1] for p in f.points)
        for f in component]) for component in boundary_comps])
    interior_inds = collect(1:length(boundary_comps))
    deleteat!(interior_inds, exterior_ind)
    return boundary_comps[interior_inds], boundary_comps[exterior_ind]
end

"""
    boundary_components_connected(comp::CellComplex)

Find the connected components of the boundary of the highest dimensional cells, assuming
the cells are connected. For each component of the boundary, create a cell complex with the
component boundary and all of its descendants. Return the cell complexes for the interior
boundaries and the exterior boundary.
"""
function boundary_components_connected(comp::CellComplex)
    interiors, exterior = boundary_components_connected(comp.cells[end])
    return map(CellComplex, interiors), CellComplex(exterior)
end

export boundary_components
"""
    boundary_components(cells::AbstractVector{Cell{N}}) where N

Find the connected components of the boundary. Return the interior boundaries and the
exterior boundary.
"""
function boundary_components(cells::AbstractVector{Cell{N}}) where N
    interior_boundaries, exterior_boundaries = Vector{Cell{N}}[], Vector{Cell{N}}[]
    for component in connected_components(cells)
        ibs, e = boundary_components_connected(component)
        append!(interior_boundaries, ibs)
        push!(exterior_boundaries, e)
    end
    return interior_boundaries, exterior_boundaries
end

"""
    boundary_components(comp::CellComplex)

Find the connected components of the boundary of the highest dimensional cells. For each
component of the boundary, create a cell complex with the component boundary and all of its
descendants. Return the cell complexes for the interior boundaries and the exterior
boundary.
"""
function boundary_components(comp::CellComplex)
    interiors, exteriors = boundary_components(comp.cells[end])
    return map(CellComplex, interiors), map(CellComplex, exteriors)
end


text:
# initialize hidden and cell arrays
function initstate(atype, hidden, batchsize)
    state = Array(Any, 2)
    state[1] = zeros(batchsize, hidden)
    state[2] = zeros(batchsize, hidden)
    return map(s->convert(atype,s), state)
end

# initialize all weights of decoder network
function initweights(o::Dict)
    w = Dict()
    srnn, w["rnn"] = rnninit(o[:embed], o[:hidden])
    w["wsoft"] = o[:winit]*randn(o[:vocabsize], o[:hidden])
    w["bsoft"] = zeros(o[:vocabsize], 1)
    w["vemb"] = o[:winit]*randn(o[:embed], o[:visual][1])
    w["wemb"] = o[:winit]*randn(o[:embed], o[:vocabsize])
    return convert_weight(o[:atype], w), srnn
end

max_stars_repo_name: ilkerkesen/kusrp

max_issues_repo_name: ilkerkesen/kusrp

max_forks_repo_name: ilkerkesen/kusrp


text:
# Client for FedDCD with linear local model
mutable struct FedDCDClient{T1<:Int64, T2<:Float64, T3<:SparseMatrixCSC{Float64, Int64}, T4<:Matrix{Float64}, T5<:Vector{Int64}, T6<:Function} <: AbstractClient
    id::T1                                  # client index
    Xtrain::T3                              # training data
    XtrainT::T3                             # Row copy
    Ytrain::T5                              # training label
    W::T4                                   # (model) primal variable
    y::T4                                   # (model) dual variable
    η::T2                                   # learning rate
    λ::T2                                   # L2 regularization parameter
    oracle!::T6                             # gradient oracle
    function FedDCDClient(id::Int64, Xtrain::SparseMatrixCSC{Float64, Int64}, Ytrain::Vector{Int64}, config::Dict{String, Real}, oracle!::Function)
        numClasses = config["num_classes"]
        λ = config["lambda"]
        η = config["learning_rate"]
        d = size(Xtrain, 2)
        W = zeros(Float64, d, numClasses)
        y = zeros(Float64, d, numClasses)
        XtrainT = copy(Xtrain')
        new{Int64, Float64, SparseMatrixCSC{Float64, Int64}, Matrix{Float64}, Vector{Int64}, Function}(id, Xtrain, XtrainT, Ytrain, W, y, η, λ, oracle!)
    end
end

# Model update on local device
function update!(
    client::FedDCDClient
)
    client.oracle!(client.Xtrain, client.XtrainT, client.Ytrain, client.λ, client.W, client.y, client.η)
end

# Get objective value
function getObjValue(
    client::FedDCDClient
)
    objValue = obj(client.Xtrain, client.Ytrain, client.W, client.λ)
    return objValue
end

# Client for FedDCD with multilayer perceptron local model
mutable struct FedDCDClientNN{T1<:Int64, T2<:Float64, T3<:SparseMatrixCSC{Float64, Int64}, T4<:Zygote.Params, T5<:Flux.OneHotArray, T6<:Function, T7<:Flux.Chain} <: AbstractClient
    id::T1                                  # client index
    Xtrain::T3                              # training data
    XtrainT::T3                             # Row copy
    Ytrain::T5                              # training label
    W::T7                                   # (model) primal variable
    y::T4                                   # (model) dual variable
    λ::T2                                   # regularization parameter
    oracle!::T6                             # gradient oracle
    function FedDCDClientNN(id::Int64, Xtrain::SparseMatrixCSC{Float64, Int64}, Ytrain::Vector{Int64}, dim::Int64, config::Dict{String, Real}, oracle!::Function)
        numClasses = config["num_classes"]
        λ = config["lambda"]
        Ytrain = Flux.onehotbatch(Ytrain, 1:numClasses)
        W = Chain( Dense(780, dim, relu, bias=false), Dense(dim, 10, bias=false), NNlib.softmax)
        y = deepcopy(params(W))
        for j = 1:length(y)
            fill!(y[j], 0.0)
        end
        XtrainT = copy(Xtrain')
        new{Int64, Float64, SparseMatrixCSC{Float64, Int64}, Zygote.Params, Flux.OneHotArray, Function, Flux.Chain}(id, Xtrain, XtrainT, Ytrain, W, y, λ, oracle!)
    end
end

# Model update on local device
function update!(
    client::FedDCDClientNN
)
    # @printf("Client %d locally update\n", client.id)
    client.oracle!(client.XtrainT, client.Ytrain, client.W, client.y, client.λ)
end

# Client for AccFedDCD with linear local model
mutable struct AccFedDCDClient{T1<:Int64, T2<:Float64, T3<:SparseMatrixCSC{Float64, Int64}, T4<:Matrix{Float64}, T5<:Vector{Int64}, T6<:Function} <: AbstractClient
    id::T1                                  # client index
    Xtrain::T3                              # training data
    XtrainT::T3                             # Row copy
    Ytrain::T5                              # training label
    W::T4                                   # (model) primal variable
    y::T4                                   # (model) dual variable 1
    z::T4                                   # (model) dual variable 2
    v::T4                                   # (model) dual variable 3
    u::T4                                   # (model) dual variable 4
    η::T2                                   # learning rate
    r::T2                                   # participation rate
    λ::T2                                   # L2 regularization parameter
    κ::T2                                   # conditional number
    a::T2                                   # hyper parameter
    b::T2                                   # hyper parameter
    oracle!::T6                             # gradient oracle
    function AccFedDCDClient(id::Int64, Xtrain::SparseMatrixCSC{Float64, Int64}, Ytrain::Vector{Int64}, config::Dict{String, Real}, oracle!::Function)
        numClasses = config["num_classes"]
        λ = config["lambda"]
        η = config["learning_rate"]
        r = config["participation_rate"]
        d = size(Xtrain, 2)
        W = zeros(Float64, d, numClasses)
        y = zeros(Float64, d, numClasses)
        z = zeros(Float64, d, numClasses)
        v = zeros(Float64, d, numClasses)
        u = zeros(Float64, d, numClasses)
        XtrainT = copy(Xtrain')
        κ = λ
        a = sqrt(κ) / (1/r + sqrt(κ))
        b = a*κ*r^2
        new{Int64, Float64, SparseMatrixCSC{Float64, Int64}, Matrix{Float64}, Vector{Int64}, Function}(id, Xtrain, XtrainT, Ytrain, W, y, z, v, u, η, r, λ, κ, a, b, oracle!)
    end
end

# Model updates on local device
function updateW!(client::AccFedDCDClient)
    # @printf("Client %d locally update\n", client.id)
    client.oracle!(client.Xtrain, client.XtrainT, client.Ytrain, client.λ, client.W, client.v, client.η)
end

function updatev!(client::AccFedDCDClient)
    a = client.a
    client.v .= ((1-a)*client.y + a*client.z)
    client.y .= copy(client.v)
end

function updateu!(client::AccFedDCDClient)
    a = client.a; b = client.b
    θ1 = a^2 / (a^2 + b); θ2 = b / (a^2 + b)
    client.u .= (θ1*client.z + θ2*client.v)
    client.z .= copy(client.u)
end

# Client for AccFedDCD with multilayer perceptron local model
mutable struct AccFedDCDClientNN{T1<:Int64, T2<:Float64, T3<:SparseMatrixCSC{Float64, Int64}, T4<:Zygote.Params, T5<:Flux.OneHotArray, T6<:Function, T7<:Flux.Chain} <: AbstractClient
    id::T1                                  # client index
    Xtrain::T3                              # training data
    XtrainT::T3                             # Row copy
    Ytrain::T5                              # training label
    W::T7                                   # (model) primal variable
    y::T4                                   # (model) dual variable 1
    z::T4                                   # (model) dual variable 2
    v::T4                                   # (model) dual variable 3
    u::T4                                   # (model) dual variable 4
    λ::T2                                   # regularization parameter
    r::T2                                   # participation rate
    κ::T2                                   # conditional number
    a::T2                                   # hyper parameter
    b::T2                                   # hyper parameter
    oracle!::T6                             # gradient oracle
    function AccFedDCDClientNN(id::Int64, Xtrain::SparseMatrixCSC{Float64, Int64}, Ytrain::Vector{Int64}, dim::Int64, config::Dict{String, Real}, oracle!::Function)
        numClasses = config["num_classes"]
        λ = config["lambda"]
        r = config["participation_rate"]
        Ytrain = Flux.onehotbatch(Ytrain, 1:numClasses)
        W = Chain( Dense(780, dim, relu, bias=false), Dense(dim, 10, bias=false), NNlib.softmax)
        y = deepcopy(params(W))
        z = deepcopy(params(W))
        v = deepcopy(params(W))
        u = deepcopy(params(W))
        for j = 1:length(y)
            fill!(y[j], 0.0)
            fill!(z[j], 0.0)
            fill!(v[j], 0.0)
            fill!(u[j], 0.0)
        end
        XtrainT = copy(Xtrain')
        κ = 30*λ
        a = sqrt(κ) / (1/r + sqrt(κ))
        b = a*κ*r^2
        new{Int64, Float64, SparseMatrixCSC{Float64, Int64}, Zygote.Params, Flux.OneHotArray, Function, Flux.Chain}(id, Xtrain, XtrainT, Ytrain, W, y, z, v, u, λ, r, κ, a, b, oracle!)
    end
end

# Model updates on local device
function updateW!(client::AccFedDCDClientNN)
    client.oracle!(client.XtrainT, client.Ytrain, client.W, client.v, client.λ)
end

function updatev!(client::AccFedDCDClientNN)
    a = client.a
    for j = 1:length(client.v)
        client.v[j] .= (1-a)*client.y[j] + a*client.z[j]
        client.y[j] .= client.v[j]
    end
end

function updateu!(client::AccFedDCDClientNN)
    a = client.a; b = client.b
    θ1 = a^2 / (a^2 + b); θ2 = b / (a^2 + b)
    for j = 1:length(client.u)
        client.u[j] .= θ1*client.z[j] + θ2*client.v[j]
        client.z[j] .= client.u[j]
    end
end


text:
@cell "Layouting" begin
    p1 = scatter(RNG.rand(10), markersize=1)
    p2 = lines(RNG.rand(10), RNG.rand(10))
    p3 = surface(0..1, 0..1, RNG.rand(100, 100))
    p4 = heatmap(RNG.rand(100, 100))
    x = 0:0.1:10
    p5 = lines(0:0.1:10, sin.(x))
    pscene = vbox(
        hbox(p1, p2),
        p3,
        hbox(p4, p5, sizes=[0.7, 0.3]),
        sizes=[0.2, 0.6, 0.2]
    )
end

@cell "Comparing contours, image, surfaces and heatmaps" begin
    N = 20
    x = LinRange(-0.3, 1, N)
    y = LinRange(-1, 0.5, N)
    z = x .* y'
    hbox(
        vbox(
            contour(x, y, z, levels=20, linewidth=3),
            contour(x, y, z, levels=0, linewidth=0, fillrange=true),
            heatmap(x, y, z),
        ),
        vbox(
            image(x, y, z, colormap=:viridis),
            surface(x, y, fill(0f0, N, N), color=z, shading=false),
            image(-0.3..1, -1..0.5, AbstractPlotting.logo())
        )
    )
end


text:
using Documenter
using GraphPlot
push!(ENV,"DFG_USE_CGDFG"=>"true")
using DistributedFactorGraphs

makedocs(
    modules = [DistributedFactorGraphs],
    format = Documenter.HTML(),
    sitename = "DistributedFactorGraphs.jl",
    pages = Any[
        "Home" => "index.md",
        "Getting Started" => [
            "DFG Data Structures" => "DataStructure.md",
            "Building Graphs" => "BuildingGraphs.md",
            "Using Graph Elements" => "GraphData.md",
            "Drawing Graphs" => "DrawingGraphs.md",
            "Quick API Reference" => "ref_api.md"
        ],
        # "DistributedFactorGraph API's" => [
        #     "LightDFGs" => "apis/LightGraphs.md",
        #     "CloudGraphsDFGs" => "apis/CloudGraphs.md",
        # ],
        "Function Reference" => "func_ref.md"
    ]
    # html_prettyurls = !("local" in ARGS),
    )

deploydocs(
    repo   = "github.com/JuliaRobotics/DistributedFactorGraphs.jl.git",
    target = "build",
    # deps   = Deps.pip("mkdocs", "python-markdown-math")
)