text:
# ============================================================
# Main Program
# ============================================================

using ProgressMeter
using BenchmarkTools
using Langevin

ks, sol, flux, uq, simTime = Langevin.initialize("./config/cavity.txt");

dt = SKS.timestep(ks, uq, sol, simTime)

simTime = 0.0
residual = zeros(axes(ks.ib.primL))

KS = ks

@time for j = 1:KS.pSpace.ny
    for i = 1:KS.pSpace.nx+1
        un = KS.vSpace.u .* flux.n1[i, j][1] .+ KS.vSpace.v .* flux.n1[i, j][2]
        ut = KS.vSpace.v .* flux.n1[i, j][1] .- KS.vSpace.u .* flux.n1[i, j][2]

        for k in axes(sol.w[1, 1], 2)
            flux.fw1[i, j][:, k], flux.fh1[i, j][:, :, k], flux.fb1[i, j][:, :, k] =
                flux_kfvs(
                    sol.h[i-1, j][:, :, k] .+
                    0.5 .* KS.pSpace.dx[i-1, j] .* sol.sh[i-1, j][:, :, k, 1],
                    sol.b[i-1, j][:, :, k] .+
                    0.5 .* KS.pSpace.dx[i-1, j] .* sol.sb[i-1, j][:, :, k, 1],
                    sol.h[i, j][:, :, k] .-
                    0.5 .* KS.pSpace.dx[i, j] .* sol.sh[i, j][:, :, k, 1],
                    sol.b[i, j][:, :, k] .-
                    0.5 .* KS.pSpace.dx[i, j] .* sol.sb[i, j][:, :, k, 1],
                    un,
                    ut,
                    KS.vSpace.weights,
                    dt,
                    0.5 * (KS.pSpace.dy[i-1, j] + KS.pSpace.dy[i, j]),
                    sol.sh[i-1, j][:, :, k, 1],
                    sol.sb[i-1, j][:, :, k, 1],
                    sol.sh[i, j][:, :, k, 1],
                    sol.sb[i, j][:, :, k, 1],
                )
        end
    end
end

@time for j = 1:KS.pSpace.ny
    for i = 1:KS.pSpace.nx+1
        un = KS.vSpace.u .* flux.n1[i, j][1] .+ KS.vSpace.v .* flux.n1[i, j][2]
        ut = KS.vSpace.v .* flux.n1[i, j][1] .- KS.vSpace.u .* flux.n1[i, j][2]

        for k in axes(sol.w[1, 1], 2)
            flux_kfvs!(
                flux.fw1[i, j][:, k],
                flux.fh1[i, j][:, :, k],
                flux.fb1[i, j][:, :, k],
                sol.h[i-1, j][:, :, k] .+
                0.5 .* KS.pSpace.dx[i-1, j] .* sol.sh[i-1, j][:, :, k, 1],
                sol.b[i-1, j][:, :, k] .+
                0.5 .* KS.pSpace.dx[i-1, j] .* sol.sb[i-1, j][:, :, k, 1],
                sol.h[i, j][:, :, k] .-
                0.5 .* KS.pSpace.dx[i, j] .* sol.sh[i, j][:, :, k, 1],
                sol.b[i, j][:, :, k] .-
                0.5 .* KS.pSpace.dx[i, j] .* sol.sb[i, j][:, :, k, 1],
                un,
                ut,
                KS.vSpace.weights,
                dt,
                0.5 * (KS.pSpace.dy[i-1, j] + KS.pSpace.dy[i, j]),
                sol.sh[i-1, j][:, :, k, 1],
                sol.sb[i-1, j][:, :, k, 1],
                sol.sh[i, j][:, :, k, 1],
                sol.sb[i, j][:, :, k, 1],
            )
        end
    end
end


text:
#
# Returns the pearson coefficient of the fit
#
function pearson(X,Y,model,fit)
  ypred = similar(Y)
  ypred = model(X,fit.param)
  R = Statistics.cor(Y,ypred)
  return R
end

ext: jl


text:
@testset "Manifold Newton" begin
    function f_1(x::Vector)
        (x[1] - 5.0)^4
    end

    function g!_1(storage::Vector, x::Vector)
        storage[1] = 4.0 * (x[1] - 5.0)^3
    end

    function h!_1(storage::Matrix, x::Vector)
        storage[1, 1] = 12.0 * (x[1] - 5.0)^2
    end
    initial_x = [0.0]

    retract_flat!(x) = x
    method = ManifoldNewton(retract_flat!)
    Optim.optimize(NonDifferentiable(f_1, initial_x), [0.0], method)
    Optim.optimize(OnceDifferentiable(f_1, g!_1, initial_x), [0.0], method)

    options = Optim.Options(store_trace = false, show_trace = false,
                            extended_trace = true)
    results = Optim.optimize(f_1, g!_1, h!_1, [0.0], method, options)
    @test_throws ErrorException Optim.x_trace(results)
    @test Optim.g_converged(results)
    @test norm(Optim.minimizer(results) - [5.0]) < 0.01

    eta = 0.9

    function f_2(x::Vector)
      (1.0 / 2.0) * (x[1]^2 + eta * x[2]^2)
    end

    function g!_2(storage::Vector, x::Vector)
      storage[1] = x[1]
      storage[2] = eta * x[2]
    end

    function h!_2(storage::Matrix, x::Vector)
      storage[1, 1] = 1.0
      storage[1, 2] = 0.0
      storage[2, 1] = 0.0
      storage[2, 2] = eta
    end

    results = Optim.optimize(f_2, g!_2, h!_2, [127.0, 921.0], method)
    @test_throws ErrorException Optim.x_trace(results)
    @test Optim.g_converged(results)
    @test norm(Optim.minimizer(results) - [0.0, 0.0]) < 0.01
    @test summary(results) == "Manifold Newton's Method"
end


text:
using GLFW
using ModernGL
using Printf

## GLFW initialization
# set up GLFW key callbacks : press Esc to escape
function key_callback(window::GLFW.Window, key::GLFW.Key, scancode::Cint, action::GLFW.Action, mods::Cint)
	key == GLFW.KEY_ESCAPE && action == GLFW.PRESS && GLFW.SetWindowShouldClose(window, true)
end

# tell GLFW to run this function whenever the framebuffer size is changed
function framebuffer_size_callback(window::GLFW.Window, w::Cint, h::Cint)
	global fb_width = w
	global fb_height = h
	println("framebuffer size: ($w, $h)")
    return nothing
end

# tell GLFW to run this function whenever the window size is changed
function window_size_callback(window::GLFW.Window, w::Cint, h::Cint)
	global width = w
	global height = h
	println("window size: ($w, $h)")
    return nothing
end

# error callback
error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
GLFW.SetErrorCallback(error_callback)

# start OpenGL
function startgl(width, height)
	@static if Sys.isapple()
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, VERSION_MAJOR)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, VERSION_MINOR)
        GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
        GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE)
    else
        GLFW.DefaultWindowHints()
    end

    @info "starting GLFW ..."
    @info GLFW.GetVersionString()

	window = GLFW.CreateWindow(width, height, "Extended Init.")
	@assert window != C_NULL "could not open window with GLFW3."

	GLFW.SetKeyCallback(window, key_callback)
	GLFW.SetWindowSizeCallback(window, window_size_callback)
	GLFW.SetFramebufferSizeCallback(window, framebuffer_size_callback)
	GLFW.MakeContextCurrent(window)
	GLFW.WindowHint(GLFW.SAMPLES, 4)

	# get version info
    renderer = unsafe_string(glGetString(GL_RENDERER))
    version = unsafe_string(glGetString(GL_VERSION))
    @info "Renderder: $renderer"
    @info "OpenGL version supported: $version"

    return window
end

# _update_fps_counter
let previous = time()
    frame_count = 0
    global function updatefps(window::GLFW.Window)
        current = time()
        elapsed = current - previous
        if elapsed > 0.25
            previous = current
            fps = frame_count / elapsed
            GLFW.SetWindowTitle(window, @sprintf("opengl @ fps: %.2f", fps))
            frame_count = 0
        end
        frame_count += 1
    end
end

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_name: Gnimuc/Videre


text:
# Factor Graph OS type utilities
#  IIF methods should direclty detect extended types from user import
# of convert in their namespace

# FIXME, upgrade to AMP instead
KDE.getPoints(dfg::AbstractDFG, lbl::Symbol) = getBelief(dfg, lbl) |> getPoints

clampStringLength(st::AbstractString, len::Int=5) = st[1:minimum([len; length(st)])]

function clampBufferString(st::AbstractString, max::Int, len::Int=minimum([max,length(st)]))
  @assert 0 <= max "max must be greater or equal to zero"
  st = clampStringLength(st, len)
  for i in len:max-1  st *= " "; end
  return st
end

"""
    $SIGNATURES

Extract contiguous string of numbers at end of a label`::Symbol` -- e.g. `:x45_4` --> "4".
Returns `(string, suffix_substring)`

Related

[`incrSuffix`](@ref)
"""
function _getSuffix(lbl::Symbol; pattern::Regex=r"\d+")
  slbl = string(lbl)
  phrase_ = slbl |> reverse |> x->match(pattern,x).match
  slbl, reverse(phrase_)
end

"""
    $SIGNATURES

Utility for incrementing or decrementing suffix numbers in DFG variable labels, e.g.
```julia
incrSuffix(:x45_4)
# returns :x45_5

incrSuffix(:x45, +3)
# returns :x48

incrSuffix(:x45_4, -1)
# returns :x45_3
```

Notes
- Change `pattern::Regex=r"\\d+"` for alternative behaviour.
"""
function incrSuffix(lbl::Symbol, val::Integer=+1; pattern::Regex=r"\d+")
  slbl, phrase = _getSuffix(lbl, pattern=pattern)
  nint = phrase |> x->(parse(Int,x)+val)
  prefix = slbl[1:(end-length(phrase))]
  Symbol(prefix, nint)
end

"""
    $SIGNATURES
Get the CommonConvWrapper for this factor.
"""
_getCCW(gfnd::GenericFunctionNodeData) = gfnd.fnc
_getCCW(fct::DFGFactor) = getSolverData(fct) |> _getCCW
_getCCW(dfg::AbstractDFG, lbl::Symbol) = getFactor(dfg, lbl) |> _getCCW

DFG.getFactorType(ccw::CommonConvWrapper) = ccw.usrfnc!

_getZDim(ccw::CommonConvWrapper) = ccw.zDim
# TODO is MsgPrior piggy backing zdim on inferdim???
_getZDim(ccw::CommonConvWrapper{<:MsgPrior}) = length(ccw.usrfnc!.infoPerCoord) # ccw.usrfnc!.inferdim

_getZDim(fcd::GenericFunctionNodeData) = _getCCW(fcd) |> _getZDim
_getZDim(fct::DFGFactor) = _getCCW(fct) |> _getZDim

DFG.getDimension(fct::GenericFunctionNodeData) = _getZDim(fct)
DFG.getDimension(fct::DFGFactor) = _getZDim(fct)

"""
    $SIGNATURES

Return the manifold on which this ManifoldKernelDensity is defined.

DevNotes
- TODO currently ignores the .partial aspect (captured in parameter `L`)
"""
getManifold(mkd::ManifoldKernelDensity{M,B,Nothing}, asPartial::Bool=false) where {M,B} = mkd.manifold
getManifold(mkd::ManifoldKernelDensity{M,B,L}, asPartial::Bool=false) where {M,B,L <: AbstractVector} = asPartial ? mkd.manifold : getManifoldPartial(mkd.manifold, mkd._partial)

"""
    $TYPEDSIGNATURES

Return the number of dimensions this factor vertex `fc` influences.
"""
getFactorDim(w...) = getDimension(w...)
# getFactorDim(fcd::GenericFunctionNodeData) = isa(_getCCW(fcd).usrfnc!, MsgPrior) ? _getCCW(fcd).usrfnc!.inferdim : Int(_getCCW(fcd).zDim)
# getFactorDim(fc::DFGFactor) = getFactorDim(getSolverData(fc))
getFactorDim(fg::AbstractDFG, fctid::Symbol) = getFactorDim(getFactor(fg, fctid))

# function _getDimensionsPartial(ccw::CommonConvWrapper)
#   # @warn "_getDimensionsPartial not ready for use yet"
#   ccw.partialDims
# end
# _getDimensionsPartial(data::GenericFunctionNodeData) = _getCCW(data) |> _getDimensionsPartial
# _getDimensionsPartial(fct::DFGFactor) = _getDimensionsPartial(_getCCW(fct))
# _getDimensionsPartial(fg::AbstractDFG, lbl::Symbol) = _getDimensionsPartial(getFactor(fg, lbl))

# """
#     $SIGNATURES
# Get `.factormetadata` for each CPT in CCW for a specific factor in `fg`.
# """
# _getFMdThread(ccw::CommonConvWrapper,
#               thrid::Int=Threads.threadid()) = ccw.cpt[thrid].factormetadata
# #
# _getFMdThread(fc::Union{GenericFunctionNodeData,DFGFactor},
#               thrid::Int=Threads.threadid()) = _getFMdThread(_getCCW(fc), thrid)
# #
# _getFMdThread(dfg::AbstractDFG,
#               lbl::Symbol,
#               thrid::Int=Threads.threadid()) = _getFMdThread(_getCCW(dfg, lbl), thrid)
# #

# extend convenience function (Matrix or Vector{P})
function manikde!(variableType::Union{InstanceType{<:InferenceVariable}, InstanceType{<:AbstractFactor}},
                  pts::AbstractVector{P};
                  kw... ) where {P <: Union{<:AbstractArray,<:Number,<:ProductRepr,<:Manifolds.ArrayPartition} }
  #
  M = getManifold(variableType)
  infoPerCoord=ones(AMP.getNumberCoords(M, pts[1]))
  return AMP.manikde!(M, pts; infoPerCoord, kw...)
end

"""
    $SIGNATURES

Return params.N measurement samples for a factor in `<:AbstractDFG`.
"""
getMeasurements(dfg::AbstractDFG, fsym::Symbol, N::Int=getSolverParams(dfg).N) = sampleFactor(dfg, fsym, N)

"""
    $SIGNATURES

Get the folder location where debug and solver information is recorded for a particular factor graph.
"""
getLogPath(opt::SolverParams) = opt.logpath
getLogPath(dfg::AbstractDFG) = getSolverParams(dfg) |> getLogPath

"""
    $SIGNATURES

Append `str` onto factor graph log path as convenience function.
"""
joinLogPath(opt::SolverParams, str...) = joinpath(getLogPath(opt), str...)
joinLogPath(dfg::AbstractDFG, str...) = joinLogPath(getSolverParams(dfg), str...)

"""
    $(SIGNATURES)

Set variable(s) `sym` of factor graph to be marginalized -- i.e. not be updated by inference computation.
"""
function setfreeze!(dfg::AbstractDFG, sym::Symbol)
  if !isInitialized(dfg, sym)
    @warn "Vertex $(sym) is not initialized, and won't be frozen at this time."
    return nothing
  end
  vert = DFG.getVariable(dfg, sym)
  data = getSolverData(vert)
  data.ismargin = true
  nothing
end
function setfreeze!(dfg::AbstractDFG, syms::Vector{Symbol})
  for sym in syms
    setfreeze!(dfg, sym)
  end
end

"""
    $(SIGNATURES)

Freeze nodes that are older than the quasi fixed-lag length defined by `fg.qfl`, according to `fg.fifo` ordering.

Future:
- Allow different freezing strategies beyond fifo.
"""
function fifoFreeze!(dfg::AbstractDFG)
  if DFG.getSolverParams(dfg).qfl == 0
    @warn "Quasi fixed-lag is enabled but QFL horizon is zero. Please set a valid window with FactoGraph.qfl"
  end

  # the fifo history
  tofreeze = DFG.getAddHistory(dfg)[1:(end-DFG.getSolverParams(dfg).qfl)]

  # check that the variable to freeze exists fix issue #966
  filter!(v->exists(dfg,v), tofreeze)

  if length(tofreeze) == 0
      @info "[fifoFreeze] QFL - no nodes to freeze."
      return nothing
  end
  @info "[fifoFreeze] QFL - Freezing nodes $(tofreeze[1]) -> $(tofreeze[end])."
  setfreeze!(dfg, tofreeze)
  nothing
end

DFG.getPoint(typ::InferenceVariable, w...;kw...) = getPoint(typeof(typ), w...;kw...)
DFG.getCoordinates(typ::InferenceVariable, w...;kw...) = getCoordinates(typeof(typ), w...;kw...)

# WIP
# _getMeasurementRepresentation(::AbstractPrior, coord::AbstractVector{<:Number}) =

"""
    $SIGNATURES

Get the ParametricPointEstimates---based on full marginal belief estimates---of a variable in the distributed factor graph.

DevNotes
- TODO update for manifold subgroups.
- TODO standardize after AMP3D

Related

[`getVariablePPE`](@ref), [`setVariablePosteriorEstimates!`](@ref), [`getVariablePPE!`](@ref)
"""
function calcPPE( var::DFGVariable,
                  varType::InferenceVariable=getVariableType(var);
                  ppeType::Type{<:MeanMaxPPE}=MeanMaxPPE,
                  solveKey::Symbol=:default,
                  ppeKey::Symbol=solveKey,
                  timestamp=now()  )
  #
  P = getBelief(var, solveKey)
  maniDef = convert(MB.AbstractManifold, varType)
  manis = convert(Tuple, maniDef) # LEGACY, TODO REMOVE
  ops = buildHybridManifoldCallbacks(manis)
  Pme = calcMean(P)  # getKDEMean(P) #, addop=ops[1], diffop=ops[2]

  # returns coordinates at identify
  Pma = getKDEMax(P, addop=ops[1], diffop=ops[2])
  # calculate point

  ## TODO make PPE only use getCoordinates for now (IIF v0.25)
  Pme_ = getCoordinates(varType,Pme)
  # Pma_ = getCoordinates(M,Pme)

  # suggested, max, mean, current time
  # TODO, poor constructor argument assumptions on `ppeType`
  ppeType(ppeKey, Pme_, Pma, Pme_, timestamp)
end

# calcPPE(var::DFGVariable; method::Type{<:AbstractPointParametricEst}=MeanMaxPPE, solveKey::Symbol=:default) = calcPPE(var, getVariableType(var), method=method, solveKey=solveKey)

"""
    $TYPEDSIGNATURES

Calculate new Parametric Point Estimates for a given variable.

Notes
- Different methods are possible, currently [`MeanMaxPPE`](@ref) `<: AbstractPointParametricEst`.

Aliases
- `calcVariablePPE`

Related

[`setPPE!`](@ref)
"""
function calcPPE( dfg::AbstractDFG,
                  label::Symbol;
                  solveKey::Symbol=:default,
                  ppeType::Type{<:AbstractPointParametricEst}=MeanMaxPPE )
  #
  var = getVariable(dfg, label)
  calcPPE(var, getVariableType(var), ppeType=ppeType, solveKey=solveKey)
end

const calcVariablePPE = calcPPE

function setThreadModel!( fgl::AbstractDFG;
                          model=IIF.SingleThreaded )
  #
  for (key, id) in fgl.fIDs
    _getCCW(fgl, key).threadmodel = model
  end
  nothing
end

"""
    $SIGNATURES

Return bool on whether a certain factor has user defined multihypothesis.

Related

[`getMultihypoDistribution`](@ref)
"""
isMultihypo(fct::DFGFactor) = isa(_getCCW(fct).hypotheses, Distribution)

"""
    $SIGNATURES

Return the categorical distributed used for multihypothesis selection in a factor.

Related

isMultihypo
"""
getMultihypoDistribution(fct::DFGFactor) = _getCCW(fct).hypotheses

"""
    $SIGNATURES

Free all variables from marginalization.
"""
function dontMarginalizeVariablesAll!(fgl::AbstractDFG)
  fgl.solverParams.isfixedlag = false
  fgl.solverParams.qfl = (2^(Sys.WORD_SIZE-1)-1)
  fgl.solverParams.limitfixeddown = false
  for sym in ls(fgl)
    setMarginalized!(fgl, sym, false)
  end
  nothing
end

"""
    $SIGNATURES

Free all variables from marginalization.

Related

dontMarginalizeVariablesAll!
"""
function unfreezeVariablesAll!(fgl::AbstractDFG)
  dontMarginalizeVariablesAll!(fgl)
end

# WIP
# function resetSolvableAllExcept!(dfg::AbstractDFG,
#                                   fltr::NothingUnion{Regex}=nothing)
#   #
#   unfreezeVariablesAll!(dfg)
# end

"""
    $SIGNATURES

Reset initialization flag on all variables in `::AbstractDFG`.

Notes
- Numerical values remain, but inference will overwrite since init flags are now `false`.
"""
function resetVariableAllInitializations!(fgl::AbstractDFG)
  vsyms = ls(fgl)
  for sym in vsyms
    setVariableInitialized!(getVariable(fgl, sym), :false)
  end
  nothing
end

"""
    $SIGNATURES

Enable defaults for fixed-lag-like operation by using smart message passing on the tree.

Notes:
- These are only default settings, and can be modified in each use case scenario.
- Default does not update downsolve through to leaves of the tree.
"""
function defaultFixedLagOnTree!(dfg::AbstractDFG,
                                len::Int=30;
                                limitfixeddown::Bool=true )
  #
  getSolverParams(dfg).isfixedlag = true
  getSolverParams(dfg).qfl = len
  getSolverParams(dfg).limitfixeddown = limitfixeddown
  getSolverParams(dfg)
end

"""
    $SIGNATURES

Return `::Tuple` with matching variable ID symbols and `Suggested` PPE values.

Related

getVariablePPE
"""
function getPPESuggestedAll(dfg::AbstractDFG,
                            regexFilter::Union{Nothing, Regex}=nothing )
  #
  # get values
  vsyms = listVariables(dfg, regexFilter) |> sortDFG
  slamPPE = map(x->getVariablePPE(dfg, x).suggested, vsyms)
  # sizes to convert to matrix
  rumax = zeros(Int, 2)
  for ppe in slamPPE
    rumax[2] = length(ppe)
    rumax[1] = maximum(rumax)
  end

  # populate with values
  XYT = zeros(length(slamPPE),rumax[1])
  for i in 1:length(slamPPE)
    XYT[i,1:length(slamPPE[i])] = slamPPE[i]
  end
  return (vsyms, XYT)
end

"""
    $SIGNATURES

Find and return a `::Tuple` of variables and distances to `loc::Vector{<:Real}`.

Related

findVariablesNearTimestamp
"""
function findVariablesNear( dfg::AbstractDFG,
                            loc::Vector{<:Real},
                            regexFilter::Union{Nothing, Regex}=nothing;
                            number::Int=3  )
  #

  xy = getPPESuggestedAll(dfg, regexFilter)
  dist = sum( (xy[2][:,1:length(loc)] .- loc').^2, dims=2) |> vec
  prm = (dist |> sortperm)[1:number]
  return (xy[1][prm], sqrt.(dist[prm]))
end

"""
    $SIGNATURES

Convenience wrapper to `DFG.loadDFG!` taking only one argument, the file name, to load a DFG object in standard format.
"""
loadDFG(filename::AbstractString) = loadDFG!(initfg(), filename)

"""
    $SIGNATURES

Find all factors that go `from` variable to any other complete variable set within `between`.

Notes
- Developed for downsolve in CSM, expanding the cliqSubFg to include all frontal factors.
"""
function findFactorsBetweenFrom(dfg::G,
                                between::Vector{Symbol},
                                from::Symbol ) where {G <: AbstractDFG}
  # get all associated factors
  allfcts = ls(dfg, from)

  # remove candidates with neighbors outside between with mask
  mask = ones(Bool, length(allfcts))
  i = 0
  for fct in allfcts
    i += 1
    # check if immediate neighbors are all in the `between` list
    immnei = ls(dfg, fct)
    if length(immnei) != length(intersect(immnei, between))
      mask[i] = false
    end
  end

  # return only masked factors
  return allfcts[mask]
end

"""
    $SIGNATURES

Return list of factors which depend only on variables in variable list in factor
graph -- i.e. among variables.

Notes
-----
* `unused::Bool=true` will disregard factors already used -- i.e. disregard where `potentialused=true`
"""
function getFactorsAmongVariablesOnly(dfg::G,
                                      varlist::Vector{Symbol};
                                      unused::Bool=true  ) where G <: AbstractDFG
  # collect all factors attached to variables
  prefcts = Symbol[]
  for var in varlist
    union!(prefcts, DFG.ls(dfg, var))
  end

  almostfcts = Symbol[]
  if unused
    # now check if those factors have already been added
    for fct in prefcts
      vert = DFG.getFactor(dfg, fct)
      if !getSolverData(vert).potentialused
        push!(almostfcts, fct)
      end
    end
  else
    almostfcts = prefcts
  end

  # Select factors that have all variables in this clique var list
  usefcts = Symbol[]
  for fct in almostfcts
    if length(setdiff(DFG.getNeighbors(dfg, fct), varlist)) == 0
      push!(usefcts, fct)
    end
  end

  return usefcts
end

"""
    $SIGNATURES

Calculate new and then set PPE estimates for variable from some distributed factor graph.

DevNotes
- TODO solve key might be needed if one only wants to update one
- TODO consider a more fiting name.
- guess it would make sense that :default=>variableNodeData, goes with :default=>MeanMaxPPE

Aliases
- `setVariablePosteriorEstimates!`

DevNotes:

JT - TODO if subfg is in the cloud or from another fg it has to be updated
it feels like a waste to update the whole variable for one field.
currently i could find mergeUpdateVariableSolverData()
might be handy to use a setter such as updatePointParametricEst(dfg, variable, solverkey)
This might also not be the correct place, if it is uncomment:
````
if (subfg <: InMemoryDFGTypes)
  updateVariable!(subfg, var)
end
```

Related

[`calcPPE`](@ref), getVariablePPE, (updatePPE! ?)
"""
function setPPE!( variable::DFGVariable,
                  solveKey::Symbol = :default,
                  ppeType::Type{T} = MeanMaxPPE,
                  newPPEVal::T = calcPPE(variable, ppeType=ppeType, solveKey=solveKey) ) where {T <: AbstractPointParametricEst}
  #
  # vnd = getSolverData(variable, solveKey)

  #TODO in the future one can perhaps populate other solver data types here by looking at the typeof ppeDict entries
  getPPEDict(variable)[solveKey] = newPPEVal

  return variable
end

function setPPE!( subfg::AbstractDFG,
                  label::Symbol,
                  solveKey::Symbol = :default,
                  ppeType::Type{T} = MeanMaxPPE,
                  newPPEVal::NothingUnion{T} = nothing )  where {T <: AbstractPointParametricEst}
  #
  variable = getVariable(subfg,label)
  # slight optimization to avoid double variable lookup (should be optimized out during code lowering)
  newppe = newPPEVal !== nothing ? newPPEVal : calcPPE(variable, solveKey=solveKey, ppeType=ppeType)
  setPPE!(variable, solveKey, ppeType, newppe)
end

const setVariablePosteriorEstimates! = setPPE!

## ============================================================================
# Starting integration with Manifolds.jl, via ApproxManifoldProducts.jl first
## ============================================================================

"""
    $SIGNATURES
Fetch and unpack JSON dictionary stored as a data blob.
"""
function fetchDataJSON(dfg::AbstractDFG, varsym::Symbol, lbl::Symbol)
  gde,rawData = getData(dfg, varsym, lbl)
  if gde.mimeType == "application/json/octet-stream"
    JSON2.read(IOBuffer(rawData))
  else
    error("Unknown JSON Blob format $(gde.mimeType)")
  end
end

#


text:
@testset "incomplete objective" begin
    function myfun!(F, J, x)
           if !(F == nothing)
               F[1] = (x[1]+3)*(x[2]^3-7)+18
               F[2] = sin(x[2]*exp(x[1])-1)
           end

           if !(J == nothing)
               J[1, 1] = x[2]^3-7
               J[1, 2] = 3*x[2]^2*(x[1]+3)
               u = exp(x[1])*cos(x[2]*exp(x[1])-1)
               J[2, 1] = x[2]*u
               J[2, 2] = u
           end
    end

    r = nlsolve(only_fj!(myfun!), [0.5, 0.5])
    @test norm(r.zero - [0.0, 1.0], Inf) < 1e-8

    function myfun(x)
        F = similar(x)
        F[1] = (x[1]+3)*(x[2]^3-7)+18
        F[2] = sin(x[2]*exp(x[1])-1)

        J = NLSolversBase.alloc_DF(x, F)
        J[1, 1] = x[2]^3-7
        J[1, 2] = 3*x[2]^2*(x[1]+3)
        u = exp(x[1])*cos(x[2]*exp(x[1])-1)
        J[2, 1] = x[2]*u
        J[2, 2] = u

        F, J
    end
    r = nlsolve(only_fj(myfun), [0.5, 0.5])
    @test norm(r.zero - [0.0, 1.0], Inf) < 1e-8
end


text:
using NorvigsTrieSpellchecker
using BenchmarkTools

function testset(lines)
    return [(right, wrong)
            for (right, wrongs) in [split(line, ":") for line in lines]
            for wrong in split(wrongs)]
end

function spelltest(t, tests; verbose=false, silent=false)
    start = time()
    good, unknown = 0, 0
    n = 0
    for (right, wrong) in tests
        n += 1
        w = correction(t, wrong)
        good += (w == right)
        if w != right
            unknown += !haskey(t, right)
            # unknown += !(right in keys(WORDS))
            if verbose
                println("correction($wrong) => $(w) ($(WORDS[w])); expected $(right) ($(WORDS[right]))")
            end
        end
    end
    dt = time() - start
    if !silent
        println("$(round((good / n), digits=3)) of $(n) correct ($(round((unknown / n), digits=3)) unknown) at $(round((n / dt), digits=3)) words per second")
    end
end

const TRIE = load_text("big.txt")

@time spelltest(TRIE, testset(readlines(open("spell-testset1.txt"))))
@time spelltest(TRIE, testset(readlines(open("spell-testset2.txt"))))

ts1 = testset(readlines(open("spell-testset1.txt")))
ts2 = testset(readlines(open("spell-testset2.txt")))

b1 = @benchmark spelltest($TRIE, $ts1, silent = true)
b2 = @benchmark spelltest($TRIE, $ts2, silent = true)

println("Benchmark testset1, average speed: ", round(length(ts1)/(sum(b1.times)/length(b1.times)/10^6/10^3), digits=3), " words per second")
println("Benchmark testset2, average speed: ", round(length(ts2)/(sum(b2.times)/length(b2.times)/10^6/10^3), digits=3), " words per second")


text:
@testset "test_simplification" begin
## SMINREAL ##
G = ss([-5 0 0 0; 0 -1 -2.5 0; 0 4 0 0; 0 0 0 -6], [2 0; 0 1; 0 0; 0 2],
       [0 3 0 0; -2 0 0 1], [0 0; 1 0])
@test sminreal(G) == G
@test sminreal(G[1, 1]) == ss(0)
@test sminreal(G[1, 2]) == ss([-1 -2.5; 4 0], [1; 0], [3 0], [0])
@test sminreal(G[2, 1]) == ss([-5], [2], [-2], [1])
@test sminreal(G[2, 2]) == ss([-6], [2], [1], [0])

using ControlSystems.DemoSystems: resonant
R = resonant()*resonant()
@test sminreal(R) == R # https://github.com/JuliaControl/ControlSystems.jl/issues/409

# https://github.com/JuliaControl/ControlSystems.jl/issues/475
A = [
 0.0  0.0   0.0;
 1.0  0.0  -0.00578297;
 0.0  1.0   0.0;
]
B = [1.0; 0.0; 0.0]
C = [0.0 0.0 1.0]
D = [0]

ss_sys = ss(A,B,C,D)
@test sminreal(ss_sys) == ss_sys

## MINREAL ##

s = tf("s")
P = [1/(s+1) 2/(s+3); 1/(s+1) 1/(s+1)]
sys = ss(P)
sysmin = minreal(sys)

@test size(sysmin.A,1) == 3 # Test that the reduction of sys worked

@test hinfnorm(sys - sysmin)[1] < 1e-15 # And that the answer is correct

sysr, Σ = balreal(sys-sysmin)
@test maximum(Σ) < sqrt(eps())

@test all(sigma(sys-sysmin, [0.0, 1.0, 2.0])[1] .< 1e-15)  # Previously crashed because of zero dimensions in tzeros

t = 0:0.1:10
y1,x1 = step(sys,t)[[1,3]]
y2,x2 = step(sysmin,t)[[1,3]]
@test sum(abs2,y1.-y2) < 1e-6 # Test that the output from the two systems are the same

## Baltrunc

    sys = let
        _A = [-2.467030078096955 -1.7579086436230489 0.2642212618287617 -1.8139967744112642 0.9581654870157803 -1.1573282803258818 -2.409554294298775 1.2006959533714947 0.2661208318398341 0.9999177692959037 -0.3558865668962797 0.9416230396873052 -0.04719968633569187 1.9048861053182367 0.04105479772871481; 1.7913628319734 -2.762572481036251 1.1195768845996525 -0.5685580615144388 0.6726614103668934 1.0083399722116977 0.5812702064616616 -0.7829110624242633 -0.5724789585031453 0.2739159931739683 0.2757716078817949 -1.1772430269658822 -0.6381809965413997 -0.5567081736321586 0.3625746395689792; -0.005143392441108892 -1.7086991517875587 -1.717585106651333 1.5009119207685369 1.8738081845828842 -1.0939749683914888 -0.5264010960315267 0.22497823448076093 1.246044372509282 -0.5303445465251139 -0.5548425351067188 0.33884707949495035 0.5101563763162115 0.20539506245514677 0.26057781689650356; 0.604556580615853 1.2138308974005823 -2.059380841591183 -1.191665990395876 1.150735331516496 0.8110103414198591 -0.154085845809975 1.1119072467025846 0.11430616747450531 1.2819431451075372 0.921519933726758 0.4697140300022479 0.005096755657399703 -0.6758568389946026 0.23258764043708394; -0.3279280635059737 -0.3899618604207743 0.12657097470583745 0.8299061544706142 -2.6661102102915484 -0.33139258711550823 -0.012925000635212677 -0.42098301162071433 1.5144416025925511 -0.796229814982033 -0.0657636569172887 1.2193694935643702 -1.700730366738557 -0.480614226127818 0.5646449944635483; 0.2450793035368965 -1.8805339894110746 0.30667982435893143 0.9479179427233335 -1.2301806768500736 -3.136669168734015 -0.2481694826668052 0.380557579742639 -1.416193374954648 0.6987044103237035 -0.8853533466891377 -1.030250851383474 1.4525483702515944 0.36288849878688056 0.385057654662975; -2.0022001785391392 -0.2835451807794339 -0.26028704624052934 0.07217513465675776 0.08414471225474451 0.024542390758292842 -2.6763138550617667 0.011591697680752557 1.7133530919443034 0.17602929220838714 -0.6726695861811227 1.5441809892864893 -0.00015597605070589597 1.301000782058375 0.03551676106571688; -0.16249013289239514 0.13198122357867412 1.3236770560172366 0.3301510019655505 0.45179991512285006 0.07486588236983925 -1.6288494827036406 -1.5392724756838398 -0.6217318765063355 0.2858777881749128 -1.4166215418299406 -0.35825226416886896 1.8697114631330003 -0.3059163232549132 1.934355556146854; -0.3097986537718366 1.9026121572125902 -0.6714075159339921 -0.2504951318869757 -1.0960812903814519 0.551662724714743 -0.5206166054513874 -0.10057595046577751 -4.293450503414131 0.24751394898682735 0.5853194271595069 1.2157066833567851 -0.5437491323644809 0.03154589360470189 -1.3601989491821653; -1.8758967387335015 -0.08090769024414653 1.703446780809511e-5 -0.15226980080851651 -0.4474269369804122 -1.0627549839203736 1.1462306659249069 0.4825852860273927 1.1227319809099574 -2.5566746127781825 0.2744599619933535 1.514041120509779 0.8781499781073471 -0.7005275854271636 0.8540288603061493; 1.3102799390114126 2.0268412603985713 -0.26445369495117343 0.1673277757620624 0.33995168379334895 1.3424350975481887 -0.5262869773852599 1.772445455267548 0.7427203416097469 -0.6894548775913741 -4.648846973081957 0.6160480706048977 0.27693416312465485 -0.5843844597444985 -0.05046410580231136; -0.2667927151895562 -0.8394567537703295 0.7356804516400021 0.021352734365649183 1.6248465860949257 1.0375631106350516 -0.44557649819583367 -1.318671094466728 0.1198337174336245 0.028769936034846715 1.1651225420369067 -2.5802262184585376 -0.000961942820959872 0.3531297512430618 0.8801838685743447; -0.2340108920478918 -0.47020830048023965 0.0905085276240309 -1.159187923078835 -0.32263237957704327 -0.20398562485896923 -1.3882045935474685 0.6799676640930128 0.34600863452034686 -0.8950627077509554 0.37192323444628994 -0.17983996993004395 -3.254198004959321 -0.12498721862068969 -0.7674453121406706; -0.48890022859302995 0.7856693831160491 0.5530727968750053 0.28679439003867857 -0.353796629578632 1.3881939385941933 0.04935762134562992 -0.6880042763129538 1.3414828049193512 -1.6887109143288497 -0.846493928570519 -2.8372692994687916 -0.8703006934200097 -2.0491003542049913 -0.04580815825737827; 0.14776326405901308 -1.0181952707299127 0.17409898161459436 0.5045755253054325 1.2756549793646714 0.0513146326521474 1.9399237561334002 -1.3740316671021375 -1.6331290401813572 -0.4569927755409647 -0.6655757376846796 -0.11273958657543114 0.17478370341836297 0.9495603386771154 -2.069076828262977]
        _B = [-0.9890109448948233; 0.94489034588656; -0.6141178465083617; 0.16165201706201165; 0.3159060979974231; 0.48695887811075705; 1.4986210850273796; 0.20368746793706913; 0.30197970526128687; -0.8327748143726169; 0.5412418165171653; -0.6001356417667288; 0.9563775528650693; 3.0913194636365113; -0.34466316244242773]
        _C = [0.3132381087641251 -0.2983781324325293 -0.27755322605976035 0.722187317484002 0.10756083023764713 -0.8432926363233019 -1.2942110080251952 0.08989673121535467 -1.6845175332374567 -2.535783822800354 -0.6469091057919979 0.6817108268640213 -1.4395236987685638 -0.3003320643423357 0.49320084051761887]
        _D = [0.0]
        ss(_A, _B, _C, _D)
    end
    sysr,Σ = baltrunc(sys, n=10)
    @test sysr.nx == 10
    @test norm(sys-sysr) < 1e-5

    sysr,_ = baltrunc(sys)
    @test norm(sys-sysr) < 0.02

    sysr,Σ = baltrunc(sys, n=3, residual=true)
    @test sysr.nx == 3
    @test dcgain(sysr)[] ≈ dcgain(sys)[] rtol=sqrt(eps())

    sys = c2d(sys, 0.01)
    sysr,Σ = baltrunc(sys, n=3, residual=true)
    @test sysr.nx == 3
    @test dcgain(sysr)[] ≈ dcgain(sys)[] rtol=sqrt(eps())

    ## Large random system
    errors = map(1:3) do _
        G = ssrand(1,1,300)
        Gr,_ = balreal(G)
        norm(G-Gr)
    end
    @test maximum(errors) < 5e-7

end


text:
module LinearShiftRegisters

export LSR_period, LCG_period

function _repeated_substring(source)
	len = length(source)

	# had to account for weird edge case in which length 2 vectors always returns themselves
	if len < 3
		len == 1 && return nothing

		s1, s2 = source
		if len == 2 && s1 == s2
			return s1
		else
			return nothing
		end
	end

	# Check candidate strings
	for i in 1:(len ÷ 2 + 1)
		repeat_count, remainder = divrem(len, i)

		# Check for no leftovers characters, and equality when repeated
		if remainder == 0 && source == repeat(source[1:i], repeat_count)
			return source[1:i]#, repeat_count
		end
	end

	return nothing
end

function _LSR_period(const_vec::AbstractVector, init_vec::AbstractVector)
	sequence_vec = reverse(init_vec)

	res = nothing
	while isnothing(res)
		new_val = mod(sum(map(i -> const_vec[i] * sequence_vec[end - (i - 1)], 1:4)), 2)
		push!(sequence_vec, new_val)
		res = _repeated_substring(sequence_vec)
	end

	return res
end

@doc raw"""
```julia
LSR_period(const_str::AbstractVector{T}, init_str::AbstractVector{R}) where {T, R <: Integer} -> AbstractVector
LSR_period(const_str::AbstractString, init_str::AbstractString) -> String
```

Given a starting binary string (`init_str`) and a constant mutating string (`const_str`), `LSR_period` will find the repeating sequence generated from these.

!!! note

    As this method looks for recurring sequences, it is possible that this function does not terminate—or maybe it does terminate, but I haven't waited long enough to find out.

---

### Description

A (linear) shift register with feedback is an arrangement of registers in a row, each register being capable of holding either the digit 1 (on) or 0 (off).  Suppose the system has ``n`` registers, ``R_1,\ldots,R_n``, and ``X_i(t)`` denotes the content of register ``R_i`` at time ``t``.  Initially, let the system be given

Given

### Examples

```julia
julia> LSR_period("1001", "0100")
"001000111101011"

julia> LSR_period("1001", "0101")
"101011001000111"
```
"""
function LSR_period(const_vec::AbstractVector{T}, init_vec::AbstractVector{R}) where {T, R <: Integer}
	return _LSR_period(const_vec, init_vec)
end
function LSR_period(const_str::AbstractString, init_str::AbstractString)
	const_vec = parse.(Int8, collect(const_str))
	init_vec = parse.(Int8, collect(init_str))

	return join(_LSR_period(const_vec, init_vec))
end

function _LCG_period(m::Integer, a::Integer, b::Integer, x0::T) where {T <: Integer}
	sequence_vec = T[x0]

	new_x, res = x0, nothing
	while isnothing(res)
		new_x =  mod(a * new_x + b, m)
		push!(sequence_vec, new_x)
		res = _repeated_substring(sequence_vec)
	end

	return res
end

@doc raw"""
```julia
LCG_period(m::Integer, a::Integer, b::Integer, x0::Integer) -> String
```

Given `m`, `a`, `b`, and `x0`, finds the period generated by these parameters from a Linear Congruential Generator, and returns this sequence as a string.

!!! note

    As this method looks for recurring sequences, it is possible that this function does not terminate—or maybe it does terminate, but I haven't waited long enough to find out.

---

### Description

Choose m, a positive integer.  Choose a, b (mod m) integers.

The Linear Congruential Generator is given to us by

```math
x_1 \equiv ax_0+b\mod m\\
x_2 \equiv ax_1+b\mod m\\
\vdots\\
x_{n + 1} \equiv ax_n + b \mod m\\
\vdots
```

Eventually, ``x_n = x_i`` for some ``n``, ``i`` (let's suppose ``n>i`` and ``i`` is not necessarily ``0``).  Then, clearly,

```math
x_{n+1} = x_{i+1},\ldots.
```

So

```math
x_i, x_{i+1},\dots, x_{n-1}
```

Will be repeated forever.

From a cryptography point of view, LCG is too vulnerable (although all pseudorandom number generators are vulnerable), but it is still useful for other perposes that require pseudorandom numbers.

For essentially a more sophisticated version of LCG in the sense that it is based on recurring linear relations, see `LSR_period`.

---

### Examples

```julia
julia> println(LinearShiftRegisters._LCG_period(9, 5, 2, 0))
[0, 2, 3, 8, 6, 5]

julia> LCG_period(9, 5, 2, 0)
"023865"
```
"""
function LCG_period(m::Integer, a::Integer, b::Integer, x0::Integer)
	return join(_LCG_period(m, a, b, x0))
end

end # end module


text:
using IMFData
using DataFrames
using DataFramesMeta
using CSV

datasets = get_imf_datasets()
# Use @where macro from DataFramesMeta to filter based on dataset name
ds_ifs = @where(datasets, occursin.("IFS", :dataset_id))
ds_dot = @where(datasets, occursin.("DOT", :dataset_id))

ifs_structure  = get_imf_datastructure("IFS")
# Search for GDP indicators
ifs_indicators = ifs_structure["Parameter Values"]["CL_INDICATOR_IFS"]
gdp_indicators = @where(ifs_indicators,
                 occursin.("Gross Domestic Product", :description),
                 occursin.("Domestic Currency", :description))
# CSV.write("ifs_gdp_indicators.csv", gdp_indicators; delim='\t')

indic = "NGDP_SA_XDC"
area  = "US"
data_available = get_ifs_data(area, indic, "Q", 1947, 2016)
data_not_available = get_ifs_data(area, indic, "M", 2000, 2015)
data_not_defined = get_ifs_data(area, "NGDP_SA", "Q", 2000, 2005)


text:
@testset "998.maximum-binary-tree-ii.jl" begin
    @test insert_into_max_tree(TreeNode(4, TreeNode(1), TreeNode(3, TreeNode(2), nothing)),
                               5) ==
          TreeNode(5, TreeNode(4, TreeNode(1), TreeNode(3, TreeNode(2), nothing)))
    @test insert_into_max_tree(TreeNode(5, TreeNode(2, nothing, TreeNode(1)), TreeNode(4)),
                               3) ==
          TreeNode(5, TreeNode(2, nothing, TreeNode(1)), TreeNode(4, nothing, TreeNode(3)))
    @test insert_into_max_tree(TreeNode(5, TreeNode(2, nothing, TreeNode(1)), TreeNode(3)),
                               4) ==
          TreeNode(5, TreeNode(2, nothing, TreeNode(1)), TreeNode(4, TreeNode(3), nothing))
end


text:
using LinearAlgebra
using FileIO, Colors, GeometryBasics
using ReferenceTests: loadasset, RNG
using Makie: Record, volume

@cell "Image on Geometry (Moon)" begin
    moon = loadasset("moon.png")
    fig, ax, meshplot = mesh(Sphere(Point3f(0), 1f0), color=moon, shading=false, show_axis=false)
    update_cam!(ax.scene, Vec3f(-2, 2, 2), Vec3f(0))
    fig
end

@cell "Image on Geometry (Earth)" begin
    earth = loadasset("earth.png")
    m = uv_mesh(Tesselation(Sphere(Point3f(0), 1f0), 60))
    mesh(m, color=earth, shading=false)
end

@cell "Orthographic Camera" begin
    function colormesh((geometry, color))
        mesh1 = normal_mesh(geometry)
        npoints = length(GeometryBasics.coordinates(mesh1))
        return GeometryBasics.pointmeta(mesh1; color=fill(color, npoints))
    end
    # create an array of differently colored boxes in the direction of the 3 axes
    x = Vec3f(0); baselen = 0.2f0; dirlen = 1f0
    rectangles = [
        (Rect(Vec3f(x), Vec3f(dirlen, baselen, baselen)), RGBAf(1, 0, 0, 1)),
        (Rect(Vec3f(x), Vec3f(baselen, dirlen, baselen)), RGBAf(0, 1, 0, 1)),
        (Rect(Vec3f(x), Vec3f(baselen, baselen, dirlen)), RGBAf(0, 0, 1, 1))
    ]

    meshes = map(colormesh, rectangles)
    fig, ax, meshplot = mesh(merge(meshes))
    scene = ax.scene
    center!(scene)
    cam = cameracontrols(scene)
    dir = widths(data_limits(scene)) ./ 2.
    dir_scaled = Vec3f(
        dir[1] * scene.transformation.scale[][1],
        0.0,
        dir[3] * scene.transformation.scale[][2],
    )
    cam.upvector[] = (0.0, 0.0, 1.0)
    cam.lookat[] = minimum(data_limits(scene)) + dir_scaled
    cam.eyeposition[] = (cam.lookat[][1], cam.lookat[][2] + 6.3, cam.lookat[][3])
    cam.attributes[:projectiontype][] = Makie.Orthographic
    cam.zoom_mult[] = 0.61f0
    update_cam!(scene, cam)
    fig
end

@cell "Volume Function" begin
    volume(RNG.rand(32, 32, 32), algorithm=:mip)
end

@cell "Textured Mesh" begin
    catmesh = loadasset("cat.obj")
    mesh(catmesh, color=loadasset("diffusemap.png"))
end

@cell "Load Mesh" begin
    mesh(loadasset("cat.obj"))
end

@cell "Colored Mesh" begin
    x = [0, 1, 2, 0]
    y = [0, 0, 1, 2]
    z = [0, 2, 0, 1]
    color = [:red, :green, :blue, :yellow]
    i = [0, 0, 0, 1]
    j = [1, 2, 3, 2]
    k = [2, 3, 1, 3]
    # indices interpreted as triangles (every 3 sequential indices)
    indices = [1, 2, 3,   1, 3, 4,   1, 4, 2,   2, 3, 4]
    mesh(x, y, z, indices, color=color)
end

@cell "Wireframe of a Mesh" begin
    wireframe(loadasset("cat.obj"))
end

@cell "Wireframe of Sphere" begin
    wireframe(Sphere(Point3f(0), 1f0))
end

@cell "Wireframe of a Surface" begin
    function xy_data(x, y)
        r = sqrt(x^2 + y^2)
        r == 0.0 ? 1f0 : (sin(r) / r)
    end
    N = 30
    lspace = range(-10, stop=10, length=N)
    z = Float32[xy_data(x, y) for x in lspace, y in lspace]
    r = range(0, stop=3, length=N)
    wireframe(r, r, z)
end

@cell "Surface with image" begin
    N = 30
    function xy_data(x, y)
        r = sqrt(x^2 + y^2)
        r == 0.0 ? 1f0 : (sin(r) / r)
    end
    xrange = range(-2, stop=2, length=N)
    surf_func(i) = [Float32(xy_data(x * i, y * i)) for x = xrange, y = xrange]
    surface(
        xrange, xrange, surf_func(10),
        color=RNG.rand(RGBAf, 124, 124)
    )
end

@cell "Meshscatter Function" begin
    large_sphere = Sphere(Point3f(0), 1f0)
    positions = decompose(Point3f, large_sphere)
    colS = [RGBAf(RNG.rand(), RNG.rand(), RNG.rand(), 1.0) for i = 1:length(positions)]
    sizesS = [RNG.rand(Point3f) .* 0.05f0 for i = 1:length(positions)]
    meshscatter(positions, color=colS, markersize=sizesS)
end

@cell "scatter" begin
    scatter(RNG.rand(20), RNG.rand(20), markersize=10)
end

@cell "Marker sizes" begin
    scatter(RNG.rand(20), RNG.rand(20), markersize=RNG.rand(20) .* 20, color=to_colormap(:Spectral, 20))
end

@cell "Record Video" begin
    f(t, v, s) = (sin(v + t) * s, cos(v + t) * s, (cos(v + t) + sin(v)) * s)
    t = Observable(0.0) # create a life signal
    limits = Rect3f(Vec3f(-1.5, -1.5, -3), Vec3f(3, 3, 6))
    fig, ax, p1 = meshscatter(lift(t -> f.(t, range(0, stop=2pi, length=50), 1), t), markersize=0.05)
    p2 = meshscatter!(ax, lift(t -> f.(t * 2.0, range(0, stop=2pi, length=50), 1.5), t), markersize=0.05)

    linepoints = lift(p1[1], p2[1]) do pos1, pos2
        map((a, b) -> (a => b), pos1, pos2)
    end

    linesegments!(ax, linepoints, linestyle=:dot)

    Record(fig, 1:2) do i
        t[] = i / 10
    end
end

@cell "3D Contour with 2D contour slices" begin
    function test(x, y, z)
        xy = [x, y, z]
        ((xy') * Matrix(I, 3, 3) * xy) / 20
    end
    x = range(-2pi, stop=2pi, length=100)
    # c[4] == fourth argument of the above plotting command
    fig, ax, c = contour(x, x, x, test, levels=6, alpha=0.3, transparency=true)

    xm, ym, zm = minimum(data_limits(c))
    contour!(ax, x, x, map(v -> v[1, :, :], c[4]), transformation=(:xy, zm), linewidth=2)
    heatmap!(ax, x, x, map(v -> v[:, 1, :], c[4]), transformation=(:xz, ym))
    contour!(ax, x, x, map(v -> v[:, :, 1], c[4]), fillrange=true, transformation=(:yz, xm))
    # reorder plots for transparency
    ax.scene.plots[:] = ax.scene.plots[[1, 3, 4, 5, 2]]
    fig
end

@cell "Contour3d" begin
    function xy_data(x, y)
        r = sqrt(x * x + y * y)
        r == 0.0 ? 1f0 : (sin(r) / r)
    end
    r = range(-1, stop=1, length=100)
    contour3d(r, r, (x, y) -> xy_data(10x, 10y), levels=20, linewidth=3)
end

@cell "Arrows 3D" begin
    function SphericalToCartesian(r::T, θ::T, ϕ::T) where T <: AbstractArray
        x = @.r * sin(θ) * cos(ϕ)
        y = @.r * sin(θ) * sin(ϕ)
        z = @.r * cos(θ)
        Point3f.(x, y, z)
    end
    n = 100^2 # number of points to generate
    r = ones(n);
    θ = acos.(1 .- 2 .* RNG.rand(n))
    φ = 2π * RNG.rand(n)
    pts = SphericalToCartesian(r, θ, φ)
    arrows(pts, (normalize.(pts) .* 0.1f0), arrowsize=0.02, linecolor=:green, arrowcolor=:darkblue)
end

@cell "Image on Surface Sphere" begin
    n = 20
    θ = [0;(0.5:n - 0.5) / n;1]
    φ = [(0:2n - 2) * 2 / (2n - 1);2]
    x = [cospi(φ) * sinpi(θ) for θ in θ, φ in φ]
    y = [sinpi(φ) * sinpi(θ) for θ in θ, φ in φ]
    z = [cospi(θ) for θ in θ, φ in φ]
    RNG.rand([-1f0, 1f0], 3)
    pts = vec(Point3f.(x, y, z))
    f, ax, p = surface(x, y, z, color=Makie.logo(), transparency=true)
end

@cell "Arrows on Sphere" begin
    n = 20
    f   = (x, y, z) -> x * exp(cos(y) * z)
    ∇f  = (x, y, z) -> Point3f(exp(cos(y) * z), -sin(y) * z * x * exp(cos(y) * z), x * cos(y) * exp(cos(y) * z))
    ∇ˢf = (x, y, z) -> ∇f(x, y, z) - Point3f(x, y, z) * dot(Point3f(x, y, z), ∇f(x, y, z))

    θ = [0;(0.5:n - 0.5) / n;1]
    φ = [(0:2n - 2) * 2 / (2n - 1);2]
    x = [cospi(φ) * sinpi(θ) for θ in θ, φ in φ]
    y = [sinpi(φ) * sinpi(θ) for θ in θ, φ in φ]
    z = [cospi(θ) for θ in θ, φ in φ]

    pts = vec(Point3f.(x, y, z))
    ∇ˢF = vec(∇ˢf.(x, y, z)) .* 0.1f0
    surface(x, y, z)
    arrows!(
        pts, ∇ˢF,
        arrowsize=0.03, linecolor=(:white, 0.6), linewidth=0.03
    )
    current_figure()
end

@cell "surface + contour3d" begin
    vx = -1:0.01:1
    vy = -1:0.01:1

    f(x, y) = (sin(x * 10) + cos(y * 10)) / 4
    fig = Figure()
    ax1 = fig[1, 1] = Axis(fig, title = "surface")
    ax2 = fig[1, 2] = Axis(fig, title = "contour3d")
    surface!(ax1, vx, vy, f)
    contour3d!(ax2, vx, vy, (x, y) -> f(x, y), levels=15, linewidth=3)
    fig
end

@cell "FEM mesh 3D" begin
    cat = loadasset("cat.obj")
    vertices = decompose(Point3f, cat)
    faces = decompose(TriangleFace{Int}, cat)
    coordinates = [vertices[i][j] for i = 1:length(vertices), j = 1:3]
    connectivity = [faces[i][j] for i = 1:length(faces), j = 1:3]
    mesh(
        coordinates, connectivity,
        color=RNG.rand(length(vertices))
    )
end

@cell "OldAxis + Surface" begin
    vx = -1:0.01:1
    vy = -1:0.01:1

    f(x, y) = (sin(x * 10) + cos(y * 10)) / 4
    scene = Scene(resolution=(500, 500), camera=cam3d!)
    # One way to style the axis is to pass a nested dictionary / named tuple to it.
    psurf = surface!(scene, vx, vy, f)
    axis3d!(scene, frame = (linewidth = 2.0,))
    center!(scene)
    # One can also directly get the axis object and manipulate it
    axis = scene[OldAxis] # get axis

    # You can access nested attributes likes this:
    axis[:names, :axisnames] = ("\\bf{ℜ}[u]", "\\bf{𝕴}[u]", " OK\n\\bf{δ}\n γ")
    tstyle = axis[:names] # or just get the nested attributes and work directly with them

    tstyle[:textsize] = 10
    tstyle[:textcolor] = (:red, :green, :black)
    tstyle[:font] = "helvetica"

    psurf[:colormap] = :RdYlBu
    wh = widths(scene)
    t = text!(
        campixel(scene),
        "Multipole Representation of first resonances of U-238",
        position=(wh[1] / 2.0, wh[2] - 20.0),
        align=(:center,  :center),
        textsize=20,
        font="helvetica"
    )
    c = lines!(scene, Circle(Point2f(0.1, 0.5), 0.1f0), color=:red, offset=Vec3f(0, 0, 1))
    psurf.converted[3][] = f.(vx .+ 0.5, (vy .+ 0.5)')
    scene
end

@cell "Fluctuation 3D" begin
    # define points/edges
    perturbfactor = 4e1
    N = 3; nbfacese = 30; radius = 0.02

    large_sphere = Sphere(Point3f(0), 1f0)
    positions = decompose(Point3f, large_sphere, 30)
    np = length(positions)
    pts = [positions[k][l] for k = 1:length(positions), l = 1:3]
    pts = vcat(pts, 1.1 .* pts + RNG.randn(size(pts)) / perturbfactor) # light position influence ?
    edges = hcat(collect(1:np), collect(1:np) .+ np)
    ne = size(edges, 1); np = size(pts, 1)
    cylinder = Cylinder(Point3f(0), Point3f(0, 0, 1.0), 1f0)
    # define markers meshes
    meshC = normal_mesh(Tesselation(cylinder, nbfacese))
    meshS = normal_mesh(Tesselation(large_sphere, 20))
    # define colors, markersizes and rotations
    pG = [Point3f(pts[k, 1], pts[k, 2], pts[k, 3]) for k = 1:np]
    lengthsC = sqrt.(sum((pts[edges[:,1], :] .- pts[edges[:, 2], :]).^2, dims=2))
    sizesC = [Vec3f(radius, radius, lengthsC[i]) for i = 1:ne]
    sizesC = [Vec3f(1) for i = 1:ne]
    colorsp = [RGBA{Float32}(RNG.rand(), RNG.rand(), RNG.rand(), 1.0) for i = 1:np]
    colorsC = [(colorsp[edges[i, 1]] .+ colorsp[edges[i, 2]]) / 2.0 for i = 1:ne]
    sizesC = [Vec3f(radius, radius, lengthsC[i]) for i = 1:ne]
    Qlist = zeros(ne, 4)
    for k = 1:ne
        ct = Cylinder(
            Point3f(pts[edges[k, 1], 1], pts[edges[k, 1], 2], pts[edges[k, 1], 3]),
            Point3f(pts[edges[k, 2], 1], pts[edges[k, 2], 2], pts[edges[k, 2], 3]),
            1f0
        )
        Q = GeometryBasics.rotation(ct)
        r = 0.5 * sqrt(1 .+ Q[1, 1] .+ Q[2, 2] .+ Q[3, 3]); Qlist[k, 4] = r
        Qlist[k, 1] = (Q[3, 2] .- Q[2, 3]) / (4 .* r)
        Qlist[k, 2] = (Q[1, 3] .- Q[3, 1]) / (4 .* r)
        Qlist[k, 3] = (Q[2, 1] .- Q[1, 2]) / (4 .* r)
    end

    rotationsC = [Vec4f(Qlist[i, 1], Qlist[i, 2], Qlist[i, 3], Qlist[i, 4]) for i = 1:ne]
    # plot
    fig, ax, meshplot = meshscatter(
        pG[edges[:, 1]],
        color=colorsC, marker=meshC,
        markersize=sizesC,  rotations=rotationsC,
    )
    meshscatter!(
        ax, pG,
        color=colorsp, marker=meshS, markersize=radius,
    )
    fig
end

@cell "Connected Sphere" begin
    large_sphere = Sphere(Point3f(0), 1f0)
    positions = decompose(Point3f, large_sphere)
    linepos = view(positions, RNG.rand(1:length(positions), 1000))
    fig, ax, lineplot = lines(linepos, linewidth=0.1, color=:black, transparency=true)
    scatter!(
        ax, positions, markersize=10,
        strokewidth=2, strokecolor=:white,
        color=RGBAf(0.9, 0.2, 0.4, 0.3), transparency=true,
    )
    fig
end

@cell "image scatter" begin
    scatter(
        1:10, 1:10, RNG.rand(10, 10) .* 10,
        rotations=normalize.(RNG.rand(Quaternionf, 10 * 10)),
        markersize=20,
        # can also be an array of images for each point
        # need to be the same size for best performance, though
        marker=Makie.logo()
    )
end

@cell "Simple meshscatter" begin
    large_sphere = Sphere(Point3f(0), 1f0)
    positions = decompose(Point3f, large_sphere)
    meshscatter(positions, color=RGBAf(0.9, 0.2, 0.4, 1), markersize=0.05)
end

@cell "Animated surface and wireframe" begin
    function xy_data(x, y)
        r = sqrt(x^2 + y^2)
        r == 0.0 ? 1f0 : (sin(r) / r)
    end

    xrange = range(-2, stop=2, length=50)
    surf_func(i) = [Float32(xy_data(x * i, y * i)) for x = xrange, y = xrange]
    z = surf_func(20)
    fig, ax, surf = surface(xrange, xrange, z)

    wf = wireframe!(ax, xrange, xrange, lift(x -> x .+ 1.0, surf[3]),
        linewidth=2f0, color=lift(x -> to_colormap(x)[5], surf[:colormap])
    )
    Record(fig, range(5, stop=40, length=3)) do i
        surf[3] = surf_func(i)
    end
end

@cell "Normals of a Cat" begin
    x = loadasset("cat.obj")
    mesh(x, color=:black)
    pos = map(decompose(Point3f, x), GeometryBasics.normals(x)) do p, n
        p => p .+ Point(normalize(n) .* 0.05f0)
    end
    linesegments!(pos, color=:blue)
    current_figure()
end

@cell "Sphere Mesh" begin
    mesh(Sphere(Point3f(0), 1f0), color=:blue)
end

@cell "Stars" begin
    stars = 100_000
    scene = Scene(backgroundcolor=:black, camera=cam2d!)
    scatter!(
        scene,
        map(i -> (RNG.randn(Point3f) .- 0.5) .* 10, 1:stars),
        color=RNG.rand(stars),
        colormap=[(:white, 0.4), (:blue, 0.4), (:yellow, 0.4)], strokewidth=0,
        markersize=RNG.rand(range(2, stop=8, length=100), stars),
    )
    update_cam!(scene, Rect3f(Vec3f(-5), Vec3f(10)))
    scene
end

@cell "Unicode Marker" begin
    scatter(Point3f[(1, 0, 0), (0, 1, 0), (0, 0, 1)], marker=[:x, :circle, :cross],
            markersize=35)
end

@cell "Merged color Mesh" begin
    function colormesh((geometry, color))
        mesh1 = normal_mesh(geometry)
        npoints = length(GeometryBasics.coordinates(mesh1))
        return GeometryBasics.pointmeta(mesh1; color=fill(color, npoints))
    end
    # create an array of differently colored boxes in the direction of the 3 axes
    x = Vec3f(0); baselen = 0.2f0; dirlen = 1f0
    rectangles = [
        (Rect(Vec3f(x), Vec3f(dirlen, baselen, baselen)), RGBAf(1, 0, 0, 1)),
        (Rect(Vec3f(x), Vec3f(baselen, dirlen, baselen)), RGBAf(0, 1, 0, 1)),
        (Rect(Vec3f(x), Vec3f(baselen, baselen, dirlen)), RGBAf(0, 0, 1, 1))
    ]

    meshes = map(colormesh, rectangles)
    mesh(merge(meshes))
end

@cell "Line GIF" begin
    us = range(0, stop=1, length=100)
    f, ax, p = linesegments(Rect3f(Vec3f(0, -1, 0), Vec3f(1, 2, 2)))
    p = lines!(ax, us, sin.(us), zeros(100), linewidth=3, transparency=true)
    lineplots = [p]
    Makie.translate!(p, 0, 0, 0)
    colors = to_colormap(:RdYlBu)
    N = 5
    Record(f, 1:N) do i
        t = i/(N/5)
        if length(lineplots) < 20
            p = lines!(
                ax,
                us, sin.(us .+ t), zeros(100),
                color=colors[length(lineplots)],
                linewidth=3
            )
            pushfirst!(lineplots, p)
            translate!(p, 0, 0, 0)
        else
            lineplots = circshift(lineplots, 1)
            lp = first(lineplots)
            lp[2] = sin.(us .+ t)
            translate!(lp, 0, 0, 0)
        end
        for lp in Iterators.drop(lineplots, 1)
            z = translation(lp)[][3]
            translate!(lp, 0, 0, z + 0.1)
        end
    end
end

@cell "Surface + wireframe + contour" begin
    N = 51
    x = range(-2, stop=2, length=N)
    y = x
    z = (-x .* exp.(-x.^2 .- (y').^2)) .* 4
    fig, ax, surfaceplot = surface(x, y, z)
    xm, ym, zm = minimum(data_limits(ax.scene))
    contour!(ax, x, y, z, levels=15, linewidth=2, transformation=(:xy, zm))
    wireframe!(ax, x, y, z, transparency=true, color=(:black, 0.1))
    center!(ax.scene) # center the Scene on the display
    fig
end

let
    struct FitzhughNagumo{T}
        ϵ::T
        s::T
        γ::T
        β::T
    end

    @cell "Streamplot 3D" begin
        P = FitzhughNagumo(0.1, 0.0, 1.5, 0.8)
        f(x, P::FitzhughNagumo) = Point3f(
            (x[1] - x[2] - x[1]^3 + P.s) / P.ϵ,
            P.γ * x[2] - x[2] + P.β,
            P.γ * x[1] - x[3] - P.β,
        )
        f(x) = f(x, P)
        streamplot(f, -1.5..1.5, -1.5..1.5, -1.5..1.5, colormap=:magma, gridsize=(10, 10), arrow_size=0.1, transparency=true)
    end
end

@cell "Volume on black background" begin
    r = LinRange(-3, 3, 100);  # our value range

    ρ(x, y, z) = exp(-(abs(x))) # function (charge density)

    fig, ax, pl = volume(
        r, r, r,          # coordinates to plot on
        ρ,                # charge density (functions as colorant)
        algorithm=:mip  # maximum-intensity-projection
    )
    ax.scene[OldAxis].names.textcolor = :gray # let axis labels be seen on dark background
    fig.scene.backgroundcolor[] = to_color(:black)
    fig
end

@cell "Depth Shift" begin
    # Up to some artifacts from fxaa the left side should be blue and the right red.
    fig = Figure(resolution = (800, 400))

    prim = Rect3(Point3f(0), Vec3f(1))
    ps  = RNG.rand(Point3f, 10) .+ Point3f(0, 0, 1)
    mat = RNG.rand(4, 4)
    A   = RNG.rand(4,4,4)

    # This generates two sets of plots each on two axis. Both axes have one set
    # without depth_shift (0f0, red) and one at ∓10eps(1f0) (blue, left/right axis).
    # A negative shift should push the plot in the foreground, positive in the background.
    for (i, _shift) in enumerate((-10eps(1f0), 10eps(1f0)))
        ax = LScene(fig[1, i], show_axis = false)

        for (color, shift) in zip((:red, :blue), (0f0, _shift))
            mesh!(ax, prim, color = color, depth_shift = shift)
            lines!(ax, ps, color = color, depth_shift = shift)
            linesegments!(ax, ps .+ Point3f(-1, 1, 0), color = color, depth_shift = shift)
            scatter!(ax, ps, color = color, markersize=10, depth_shift = shift)
            text!(ax, "Test", position = Point3f(0, 1, 1.1), color = color, depth_shift = shift)
            surface!(ax, -1..0, 1..2, mat, colormap = (color, color), depth_shift = shift)
            meshscatter!(ax, ps .+ Point3f(-1, 1, 0), color = color, depth_shift = shift)
            # # left side in axis
            heatmap!(ax, 0..1, 0..1, mat, colormap = (color, color), depth_shift = shift)
            # # right side in axis
            image!(ax, -1..0, 1..2, mat, colormap = (color, color), depth_shift = shift)
            p = volume!(ax, A, colormap = (:white, color), depth_shift = shift)
            translate!(p, -1, 0, 0)
            scale!(p, 0.25, 0.25, 0.25)
        end

        center!(ax.scene)
    end
    fig
end

@cell "Order Independent Transparency" begin
    # top row (yellow, cyan, magenta) contains stacks with the same alpha value
    # bottom row (red, green, blue) contains stacks with varying alpha values
    fig = Figure()
    ax = LScene(fig[1, 1])
    r = Rect2f(-1, -1, 2, 2)
    for x in (0, 1)
        for (i, a) in enumerate((0.25, 0.5, 0.75, 1.0))
            ps = [Point3f(a, (0.15 + 0.01y)*(2x-1) , 0.2y) for y in 1:8]
            if x == 0
                cs = [RGBAf(1, 0, 0, 0.75), RGBAf(0, 1, 0, 0.5), RGBAf(0, 0, 1, 0.25)]
            elseif x == 1
                cs = [RGBAf(1, x, 0, a), RGBAf(0, 1, x, a), RGBAf(x, 0, 1, a)]
            end
            idxs = [1, 2, 3, 2, 1, 3, 1, 2, 1, 2, 3][i:7+i]
            meshscatter!(
                ax, ps, marker = r,
                color = cs[idxs], transparency = true
            )
        end
    end
    cam = cameracontrols(ax.scene)
    cam.attributes.fov[] = 22f0
    update_cam!(ax.scene, cam, Vec3f(0.625, 0, 3.5), Vec3f(0.625, 0, 0), Vec3f(0, 1, 0))
    fig
end

@cell "space 3D" begin
    fig = Figure()
    for ax in [LScene(fig[1, 1]), Axis3(fig[1, 2])]
        mesh!(ax, Rect3(Point3f(-10), Vec3f(20)), color = :orange)
        mesh!(ax, Rect2f(0.8, 0.1, 0.1, 0.8), space = :relative, color = :blue, shading = false)
        linesegments!(ax, Rect2f(-0.5, -0.5, 1, 1), space = :clip, color = :cyan, linewidth = 5)
        text!(ax, "Clip Space", position = Point2f(0, 0.52), align = (:center, :bottom), space = :clip)
        image!(ax, 0..40, 0..800, [x for x in range(0, 1, length=40), _ in 1:10], space = :pixel)
    end
    fig
end


text:
# This file contains a function to generate time series

"""
    @array(array::AbstractString)

Wrap an array for zero-based indexing.

# Usage

```jldoctest; output = false
using OffsetArrays
using Scats.Internal.Input: InputStruct
using Scats.Internal.Gen: @array

file, _ = mktemp()
Input = InputStruct()

Input.example(file)
Input(file)

N₋₁ = Input.N - 1
@array("Input.t") === OffsetArray(Input.t, 0:N₋₁)

# output

true
```
"""
macro array(array::AbstractString)
    return esc(Meta.parse("OffsetArray($array, 0:N₋₁)"))
end

"""
    gen!(Gen::GenStruct, Input::InputStruct)

Generate time series for an instance of [`InputStruct`](@ref) using generator
parameters from an instance of [`GenStruct`](@ref).

# Usage

```jldoctest; output = false
using Scats
s = Scats.API()
file, _ = mktemp()
s.Gen.example(file)
s.gen!()

# output

```
"""
function gen!(Gen::GenStruct, Input::InputStruct)

    @unpack "Gen" "N" "Δt" "q" "α" "β" "r" "A" "ν" "ϕ" "γ"

    # Save some values in input
    @pack "Input" "N" "Δt" "q"

    # Calculate standard deviation
    σ = √(sum(A .* A) / (2 * γ))

    # Calculate auxiliary variable
    N₋₁ = N - 1

    # Check sizes of arrays
    if size(Input.t, 1) != N
        Input.t = Vector{RT}(undef, N)
    end
    if size(Input.x, 1) != N
        Input.x = Vector{RT}(undef, N)
    end

    # Wrap arrays from input
    ta = @array "Input.t"
    xa = @array "Input.x"

    # Generate random numbers array
    rng = MersenneTwister()
    rand = @array "randn(rng, RT, N)"

    # Generate time series
    for k in 0:N₋₁

        # Calculate current time moment
        t = Δt * k

        # Write it in the array
        ta[k] = t

        # Calculate the linear trend part
        c = α + β * t

        # Add the harmonics part
        for l in 1:r
            c = c + A[l] * cos(2 * π * ν[l] * t - ϕ[l])
        end

        # Add the noise part
        c = c + rand[k]

        # Write a result in the array
        xa[k] = c

    end

end


text:
"""
	EconFormula(formula::Formula)

	EconFormula is a struct which is composed of various formulas used in
	the construction of a regression model for econometrics.

	It parses a formula with the following syntax:
		response ~ exogenous + (endogenous ~ instruments) ~ absorb

	Returns a struct which holds
		representation::String
		exogenous::Formula
		endogenous::Formula
		instruments::Formula
		absorb::Formula
		clusters::Formula
"""
mutable struct EconFormula
	representation::String
	exogenous::Formula
	endogenous::Formula
	instruments::Formula
	absorb::Formula
	clusters::Formula
	function EconFormula(formula::Formula)
		Representation, Exogenous, Endogenous, Instruments, Absorb, Cluster = formulaparser(formula)
		new(Representation, Exogenous, Endogenous, Instruments, Absorb, Cluster)
	end
end
endoinst(obj::Any) = nothing
function endoinst(obj::Expr)
	if (obj.head == :call) & (obj.args[1] == :(~))
		return (obj.args[2], obj.args[3])
	end
end
absorb(obj::Any) = nothing
function absorb(obj::Expr)
	if (obj.head == :(=)) & (obj.args[1] == :absorb)
		return obj.args[2]
	end
end
cluster(obj::Any) = nothing
function cluster(obj::Expr)
	if (obj.head == :(=)) & (obj.args[1] == :cluster)
		return obj.args[2]
	end
end
function formulaparser(formula::Formula)
	formula = deepcopy(formula)
	Representation = deepcopy(formula)
	response = formula.lhs
	formula = formula.rhs
	if isa(formula, Symbol)
		Representation = string(Representation)
		Exogenous = Formula(response, formula)
		Endogenous = Formula(response, 0)
		Instruments = Formula(response, 0)
		Absorb = Formula(response, 0)
		Cluster = Formula(response, 0)
	else
		EndoInst = endoinst.(formula.args)
		Absorb = absorb.(formula.args)
		Cluster = cluster.(formula.args)
		Representation.rhs.args =
			Representation.rhs.args[(EndoInst .== nothing) .&
			(Absorb .== nothing) .& (Cluster .== nothing)]
		if length(Representation.rhs.args) == 2
			setfield!(Representation, :rhs, Representation.rhs.args[2])
		end
		Representation = string(Representation)
		formula.args = formula.args[(EndoInst .== nothing) .&
			(Absorb .== nothing) .& (EndoInst .== nothing) .& (Cluster .== nothing)]
		if (formula.args[1] == :(+)) & (length(formula.args) == 2)
			Exogenous = Formula(response, formula.args[2])
		else
			Exogenous = Formula(response, formula)
		end
		filter!(elem -> elem != nothing, EndoInst)
		filter!(elem -> elem != nothing, Absorb)
		filter!(elem -> elem != nothing, Cluster)
		if isempty(EndoInst)
			Endogenous = Formula(response, 0)
			Instruments = Formula(response, 0)
		else
			EndoInst = EndoInst[1]
			Endogenous = Formula(response, EndoInst[1])
			Instruments = Formula(response, EndoInst[2])
		end
		if isempty(Absorb)
			Absorb = Formula(response, 0)
		else
			Absorb = Formula(response, Absorb[1])
		end
		if isempty(Cluster)
			Cluster = Formula(response, 0)
		else
			Cluster = Formula(response, Cluster[1])
		end
	end
	return (Representation, Exogenous, Endogenous, Instruments, Absorb, Cluster)
end


text:
import Base.start, Base.done, Base.next, Base.put!, Base.take!
export ctask

type CTask
    name::String
    inq::SyncObjRef
    outq::SyncObjRef
    rr::RemoteRef               # Holds the task
    CTask(name) = new(name, channel(), channel(), RemoteRef())
end

register_ct(name::String, ct::CTask) = put!(MessageUtils.dnsref(), (:ctask, name), ct)
deregister_ct(name::String) = take!(MessageUtils.dnsref(), (:ctask, name))
fetch_ct(name::String; timeout=0.001) = fetch(MessageUtils.dnsref(), (:ctask, name); timeout=timeout)

# Should only be executed on the node where the task is expected to be started
function ctwrap(ct, f)
    task_local_storage(:CTASK, ct)
    try
        f()

    finally
        if (ct.name != "")
             deregister_ct(ct.name)   # Handle error condition where entry is not found
        end
    end
    nothing
end

function ctask_start(f, name)
    ct = CTask(name)
    if (name != "")
        register_ct(name, ct)
    end
    t=Task(()->ctwrap(ct,f))
    put!(ct.rr, t)
    schedule(t)
    ct
end

function ctask(f::Function; pid=myid(), name="")
    if pid == myid()
        ctask_start(f, name)
    else
        remotecall_fetch(pid, () -> ctask_start(f, name))
    end
end

start(ct::CTask) = nothing
done(ct::CTask, state) = false
next(ct::CTask, state) = take!(ct.inq)

ctask_self() = task_local_storage(:CTASK)

# Sending a message to a task
(|>)(msg::Tuple, ctname::String) = put!(ctname, msg)
(|>)(msg::Tuple, ct::CTask) = put!(ct, msg)

put!(msg::Tuple; timeout=0.0) = (put!(ctask_self().outq, msg; timeout=timeout); nothing)
put!(ct::CTask, msg::Tuple; timeout=0.0) = (put!(ct.inq, msg; timeout=timeout); nothing)
put!(ctname::String, msg::Tuple; timeout=0.0) = put!(fetch_ct(ctname), msg; timeout=timeout)

# Pulling a message from a task
take!(; timeout=0.0) = take!(ctask_self().inq; timeout=timeout)
take!(ct::CTask; timeout=0.0) = take!(ct.outq; timeout=timeout)
take!(ctname::String; timeout=0.0) = take!(fetch_ct(ctname); timeout=timeout)


text:
# Copyright (c) 2020-2021, Bank of Canada
# All rights reserved.

try    # ============================================

using LightXML

if !haskey(Base.ENV, "FAME")
    error("FAME environment variable is not set!")
end

# The CHLI help file containing the table of codes
help_file = joinpath(Base.ENV["FAME"], "help", "chli", "chli_status_codes.htm")
help_file = Base.Filesystem.realpath(help_file)

function do_row(el_tr)
    td_vals = map(content, get_elements_by_tagname(el_tr, "td"))
    @assert length(td_vals) == 3
    return (strip(td_vals[1]), parse(Int32, td_vals[2]), strip(td_vals[3]))
end

find_table(els::Array) = isempty(els) ? nothing : (z=filter(!isequal(nothing), find_table.(els)); length(z) == 0 ? nothing : length(z) == 1 ? z[1] : z)
find_table(el::XMLElement) = name(el) == "table" ? el : find_table( collect(child_elements(el)) )

# main loop
open("./FAMEMessages.jl", "w") do f
    println(f, "\n# This file is autogenerted.  Do not edit.")
    println(f, "\n\nchli_help_file = \"", help_file, "\"")
    # Read and parse the help file
    xdoc = parse_file(help_file)
    table = find_table(root(xdoc))
    @assert table !== nothing && isa(table, XMLElement)
    table_rows = get_elements_by_tagname(table, "tr")
    table_rows = table_rows[2:end]
    # Print codes and messages to file
    println(f, "\nchli_status_description = Dict{Int32, String}(")
    for (hval, code, msg) in map(do_row, table_rows)
        println(f, "    $code => \"$msg\",")
    end
    println(f, ")\n\n")
    for (hval, code, msg) in map(do_row, table_rows)
        println(f, "const $hval = Int32($code)")
    end
    println(f)
end

catch e # ============================================
    @error "$(sprint(showerror, e))"
    isfile("./FAMEMessages.jl") && rm("./FAMEMessages.jl")
    # use a stub - these are the status codes we explicitly use in our code.
    open("./FAMEMessages.jl", "w") do f
        println(f, "\n# This file is autogenerted.  Do not edit.")
        println(f, "\n\nchli_help_file = ", "not found")
        println(f, "\nchli_status_description = Dict{Int32, String}(")
        println(f, "    0 => ", "Success.")
        println(f, "    13 => ", "The given object does not exist.")
        println(f, "    67 => ", "Bad option.")
        println(f, ")\n\n")
        println(f, "global const HSUCC = Int32(0)")
        println(f, "global const HNOOBJ = Int32(13)")
        println(f, "global const HBOPT = Int32(67)")
        println(f)
    end
end


text:
#**************************************************************************************
# calculations.jl
# =============== part of the GeoEfficiency.jl package.
#
# this file contains all the required function to calculate the Geometrical Efficiency.
#
#**************************************************************************************

const relativeError = 0.0001			#set the global relative precession of the Geometrical Efficiency Calculation

"""

	GeoEff_Pnt(Detector::CylDetector, aPnt::Point)

return the Geometrical Efficiency for a point source `aPnt` on front of the cylindrical detector `Detector` face.

`Throw` an Error if the point is out of cylinderical detector the face.

This is the base function that all other function call directly or indirectly to calculate Geometrical Efficiency.
"""
function GeoEff_Pnt(Detector::CylDetector, aPnt::Point)

	function MaxPhi(theta::Float64 )
		side = aPnt.Height * sin(theta)
		return clamp((aPnt.Rho^2 + side^2 - Detector.CryRadius^2 )/ side / aPnt.Rho /2.0, -1.0, 1.0) |> acos
	end # function

	func(theta::Float64 ) = MaxPhi(theta) * sin(theta)

	if 0.0 == aPnt.Rho
		strt = 0.0
		fine = atan2(Detector.CryRadius , aPnt.Height)
		return quadgk(sin, strt, fine, reltol = relativeError)[1]

	else
		strt = 0.0
		transtion = atan2(Detector.CryRadius - aPnt.Rho, aPnt.Height)
		fine = atan2(Detector.CryRadius + aPnt.Rho, aPnt.Height)
		if transtion >= 0.0

			return quadgk(sin, strt, transtion, reltol = relativeError)[1] + quadgk(func, transtion, fine, reltol = relativeError)[1] / pi

		else
			Error("Point off-axis: out of the Detector face, this case is not implemented yet")

		end #if

	end #if
end #function

"""

	GeoEff_Disk(Detector::CylDetector, SurfacePnt::Point, SrcRadius::Real)

return the Geometrical Efficiency for a disk source. The disk center is `SurfacePnt` and its radius is `SrcRadius` on front of the cylindrical detector `Detector` face.

`Throw` an Error if the disk is out of cylindrical detector the face.
"""
function GeoEff_Disk(Detector::CylDetector, SurfacePnt::Point, SrcRadius::Real)
	integrand(xRho) = xRho * GeoEff_Pnt(Detector, setRho!(SurfacePnt, xRho))
	return  quadgk(integrand , 0, SrcRadius, reltol = relativeError)[1] / SrcRadius^2
end #function

"""

	geoEff(Detector::CylDetector, aSurfacePnt::Point, SrcRadius::Real = 0.0, SrcLength::Real = 0.0)

return the Geometrical Efficiency for a source (point, disk or cylinder) with the cylindrical detector `Detector`.

example:-

	newDet = CylDetector()
	eff = geoEff(newDet, aSurfacePNT, SrcRadius, SrcLength)

`aSurfacePNT`: a surface point.
if `SrcRadius` = `SrcLength` = `0`; the method returns the Geometrical Efficiency at this point.

`SrcRadius`: Radius of the source.
if `srcHieght` = 0; the method returns Geometrical Efficiency for disc of Radius = `SrcRadius` and
its center is defined by the `aSurfacePNT`.

`srcHieght`:  the height of an upright cylinder source having a base like described above.

>`Throw` an Error if the source location is inappropriate.

`Note please`
\n- `aSurfacePnt`: point height is consider to be measured from the detector face surface.
\n*****
"""
function geoEff(Detector::CylDetector, aSurfacePnt::Point, SrcRadius::Real = 0.0, SrcLength::Real = 0.0)
	pnt::Point = deepcopy(aSurfacePnt)
	@assert Detector.CryRadius > SrcRadius		" Source Raduis: Expected less than 'Detector Raduis=$(Detector.CryRadius)', get $SrcRadius."

    if 0.0 == SrcRadius
		@assert Detector.CryRadius  > pnt.Rho	" Point off-axis: Expected less than 'Detector Raduis=$(Detector.CryRadius)', get $(pnt.Rho)."
        return GeoEff_Pnt(Detector, pnt)/2            	#Point source

	elseif 0.0 == SrcLength								#Disk source
        return GeoEff_Disk(Detector, pnt, SrcRadius)

	else												#Cylinderical source
        integrand(xH) = GeoEff_Disk(Detector, setHeight!(pnt, xH), SrcRadius)
		return quadgk(integrand , pnt.Height, pnt.Height + SrcLength, reltol = relativeError)[1] / SrcLength

	end #if
end #function

"""

	geoEff(Detector::RadiationDetector = detectorFactory())

return the Geometrical Efficiency of the given detector or if no detector is supplied it ask for a detector from the `console`. Also prompt the user to input a source via the `console`.

>`Throw` an Error if the source location is inappropriate.
\n*****
"""
function geoEff(Detector::RadiationDetector = detectorFactory())
	geoEff(Detector, source()...)
end #function

"""

	geoEff(Detector::BoreDetector, aCenterPnt::Point, SrcRadius::Real = 0.0, SrcLength::Real = 0.0)

return the Geometrical Efficiency for the given source (point , disk or cylinder) with the Bore-Hole detector `Detector`.

example:-

	newDet = BoreDetector()
	efff = geoEff(newDet, aCenterPNT, SrcRadius, SrcLength)

`aCenterPNT`: a center point represent the anchoring point of the source.
if `SrcRadius` = `SrcLength` = `0`; the method returns the Geometrical Efficiency at the anchoring point.

`SrcRadius`: Raduis of the source.
if SrcLength = 0;  the method returns Geometrical Efficiency for disc of Radius = `SrcRadius`
and its center is defined by the `aCenterPNT`.

`SrcLength`: the height of an upright cylinder source having a base like described above.

>`Throw` an Error if the source location is inappropriate.

`Note please`
\n-  `aCenterPNT` : point `height` is consider to be measured from the detector middle, +ve value are above the detector center while -ve are below.
\n*****
"""
function geoEff(Detector::BoreDetector, aCenterPnt::Point, SrcRadius::Real = 0.0, SrcLength::Real = 0.0)

	HeightWup = aCenterPnt.Height - Detector.CryLength/2.0
	HeightWdown = aCenterPnt.Height + Detector.CryLength/2.0
	if HeightWdown < 0.0
		if HeightWup + SrcLength < 0.0 		#invert the source.
			return geoEff(Detector, Point(aCenterPnt.Height - Detector.CryLength, aCenterPnt.Rho), SrcRadius, SrcLength)

		else # the source span the detector and emerges from both sides, split the source into two sources.
			#res = (1 - 2 * geoEff(detin, Point(0.0), SrcRadius, SrcLength))* Detector.CryLength /SrcLength
			res = geoEff(Detector, Point(0.0), SrcRadius, -aCenterPnt.Height )* (-aCenterPnt.Height /SrcLength)
			res += geoEff(Detector, Point(0.0), SrcRadius, SrcLength + aCenterPnt.Height )* (1.0 + aCenterPnt.Height /SrcLength)
			return res

		end
	end

	pntWup::Point = deepcopy(aCenterPnt);
	setHeight(aCenterPnt, abs(HeightWup));  #0.0 == SrcRadius && setRho!(pntWup, 0.0)

	pntWdown::Point = deepcopy(aCenterPnt);
	setHeight!(pntWdown, abs(HeightWdown)); #0.0 == SrcRadius && setRho!(pntWdown, 0.0)

	detin::CylDetector = CylDetector(Detector.HoleRadius)
	detout::CylDetector = CylDetector(Detector.CryRadius)

	if HeightWup >= 0.0						# the source as a whole out of detector
		res = geoEff(detout, pntWup, SrcRadius, SrcLength)[1]
		res -= geoEff(detin, pntWdown, SrcRadius, SrcLength)[1]

	elseif HeightWup + SrcLength < 0.0 		# the source as a whole in the detector
		res = 1 - geoEff(detin, setHeight(pntWup,abs(HeightWup + SrcLength)), SrcRadius, SrcLength)[1]
		res -= geoEff(detin, pntWdown, SrcRadius, SrcLength)[1]

	else # elseif SrcLength > 0.0
		res = (1.0 - geoEff(detin, Point(0.0), SrcRadius, -HeightWup)[1] )* -HeightWup/SrcLength
		res += geoEff(detout, Point(0.0), SrcRadius, HeightWup + SrcLength)[1] * (1.0 + HeightWup/SrcLength)
		return res

	#=else
		return 1.0 - geoEff(detin, setHeight!(pnt, -Height), SrcRadius)[1]
	else
		res = 1 - quadgk(xH -> GeoEff_Disk(detin, setHeight!(pnt, xH), SrcRadius), 0.0, -pnt.Height, reltol = relativeError)[1]
		res = res + quadgk(xH -> GeoEff_Disk(detout, setHeight!(pntWup, xH), SrcRadius), 0.0, pntWup.Height , reltol = relativeError)[1]
			=#
	end #if

	return res
end #function

"""

	geoEff(Detector::WellDetector, aWellPnt::Point, SrcRadius::Real = 0.0, SrcLength::Real = 0.0)

return the Geometrical Efficiency for the given source (point, disk or cylinder) with the Well-Type detector `Detector`.

Example:-

	newDet = WellDetector()
	eff = geoEff(newDet, aWellPNT, SrcRadius, SrcLength)

`aWellPNT`: a Well point represent the anchoring point of the source.
if `SrcRadius` = `SrcLength` = `0`; the method returns the Geometrical Efficiency at the anchoring point.

`SrcRadius`: Raduis of the source.
if SrcLength = 0;  the method returns Geometrical Efficiency for disc of Radius = `SrcRadius`
and its center is defined by the `aWellPNT`.

SrcLength:  the height of upright cylinder source having a base like described above.

>`Throw` an Error if the source location is inappropriate.

`Note Please`
\n- `aWellPNT` : point `height` is considered to be measured from the detector hole surface.
\n*****
"""
function geoEff(Detector::WellDetector, aWellPnt::Point, SrcRadius::Real = 0.0, SrcLength::Real = 0.0)
	pnt::Point = deepcopy(aWellPnt)
	Height = pnt.Height - Detector.HoleDepth

	detin::CylDetector = CylDetector(Detector.HoleRadius)
	detout::CylDetector = CylDetector(Detector.CryRadius)
	setHeight!(pnt, Height); #0.0 == SrcRadius && setRho!(pnt, 0.0)

	if Height > 0.0							# the source as a whole out of the detector
		return geoEff(detout, setHeight!(pnt, Height), SrcRadius, SrcLength)[1]

	elseif Height + SrcLength < 0.0 		# the source as a whole inside of the detector
		return 1.0 - geoEff(detin, setHeight!(pnt, -(Height + SrcLength) ), SrcRadius, SrcLength)[1]

	elseif SrcLength > 0.0
		res = (1.0 - geoEff(detin, Point(0.0), SrcRadius, -Height)[1] )* -Height/SrcLength
		res += geoEff(detout, Point(0.0), SrcRadius, Height + SrcLength)[1] * (1.0 + Height/SrcLength)
		return res

	else
		return 1.0 - geoEff(detin, setHeight!(pnt, -Height), SrcRadius)[1]

	end #if
end #function


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

using Test

using QuantumCircuits
using QuantumCircuits.QCircuits.Circuit: getCode, toQiskit
using QuantumCircuits.QCircuits.Math

n = 7
qr = QuantumRegister(n)
qc = QCircuit(qr)
qc.h(qr)
qc.barrier()
qc.x(qr)

@test length(getCode(qc)) == 2*n+1

eqr = QuantumRegister(n)
eqc = QCircuit(eqr)
eqc.h(eqr)
eqc.x(eqr)

toQiskit(qc)
@test unitary_error(tomatrix(qc), tomatrix(eqc)) < 1e-16
@test unitary_error(tomatrix(inv(qc)), tomatrix(inv(eqc))) < 1e-16
@test unitary_error(tomatrix(decompose(qc)), tomatrix(decompose(eqc))) < 1e-16
@test unitary_error(tomatrix(simplify(qc)), tomatrix(simplify(eqc))) < 1e-16


text:
using LCMCore
using CaesarLCMTypes

function typed_callback(channel::String, msg::example_t)
    @show channel
    @show msg
end
lcm=LCM()
subscribe(lcm, "TEST_CHANNEL", typed_callback, example_t)

while true
    handle(lcm)
end


text:
abstract type TrustRegionUpdater end
struct TrustRegion{M,SP,D}
    scheme::M
    spsolve::SP
    Δupdate::D
end
summary(tr::TrustRegion) = "$(summary(modelscheme(tr))) with $(summary(algorithm(tr)))"
function initial_preconditioner(approach::TrustRegion, x)
    nothing
end
"""
  BTR() <: TrustRegionUpdater

Basic trust region updater following, and named after [CGT].
"""
struct BTR{T}
    Δmin::T
end
TrustRegion(; deltamin = nothing) = TrustRegion(Newton(), NTR(), BTR(deltamin))
TrustRegion(m, sp = NTR(); deltamin = nothing) = TrustRegion(m, sp, BTR(deltamin))
modelscheme(tr::TrustRegion) = tr.scheme
algorithm(tr::TrustRegion) = tr.spsolve

solve(problem::OptimizationProblem, x0, scheme, options::OptimizationOptions) =
    solve(problem, (x0, nothing), TrustRegion(scheme, NWI()), options)

solve(
    problem::OptimizationProblem{<:Any,<:Nothing,<:Nothing,<:Nothing},
    x0,
    approach::TrustRegion,
    options::OptimizationOptions,
) = solve(problem, (x0, nothing), approach, options)

solve(
    problem::OptimizationProblem{<:Any,<:Nothing,<:Nothing,<:Nothing},
    s0::Tuple,
    approach::TrustRegion,
    options::OptimizationOptions,
) = solve(problem, s0, approach, options)

function solve(
    problem::OptimizationProblem,
    x0,
    scheme::Newton,
    options::OptimizationOptions,
)
    solve(problem, (x0, nothing), TrustRegion(scheme, NTR()), options)
end
function solve(
    problem::OptimizationProblem,
    x0,
    approach::TrustRegion,
    options::OptimizationOptions,
)
    solve(problem, (x0, nothing), approach, options)
end
include("optimize/inplace_loop.jl")
include("optimize/outofplace_loop.jl")

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
function blocking_test_prep(obs=ℝ.([1.0, 1.2, 0.8, 1.3, 2.0]),
                            tt=[0.0, 1.0, 1.5, 2.3, 4.0],
                            knots=collect(1:length(obs)-2)[1:1:end],
                            change_pt_buffer=100)
    θ₀ = [10.0, -8.0, 25.0, 0.0, 3.0]
    P˟ = BSI.FitzhughDiffusion(θ₀...)
    P̃ = [BSI.FitzhughDiffusionAux(θ₀..., t₀, u[1], T, v[1]) for (t₀,T,u,v)
         in zip(tt[1:end-1], tt[2:end], obs[1:end-1], obs[2:end])]
    L = @SMatrix [1. 0.]
    Σdiagel = 10^(-10)
    Σ = @SMatrix [Σdiagel]

    Ls = [L for _ in P̃]
    Σs = [Σ for _ in P̃]
    τ(t₀,T) = (x) ->  t₀ + (x-t₀) * (2-(x-t₀)/(T-t₀))
    m = length(obs) - 1
    P = Array{BSI.ContinuousTimeProcess,1}(undef,m)
    dt = 1/50
    for i in m:-1:1
        num_pts = Int64(ceil((tt[i+1]-tt[i])/dt))+1
        t = τ(tt[i], tt[i+1]).( range(tt[i], stop=tt[i+1], length=num_pts) )
        P[i] = ( (i==m) ? BSI.GuidPropBridge(Float64, t, P˟, P̃[i], Ls[i], obs[i+1], Σs[i];
                                         change_pt=BSI.NoChangePt(change_pt_buffer),
                                         solver=BSI.Vern7()) :
                          BSI.GuidPropBridge(Float64, t, P˟, P̃[i], Ls[i], obs[i+1], Σs[i],
                                         P[i+1].H[1], P[i+1].Hν[1], P[i+1].c[1];
                                         change_pt=BSI.NoChangePt(change_pt_buffer),
                                         solver=BSI.Vern7()) )
    end

    T = SArray{Tuple{2},Float64,1,2}
    blocking_params = (knots, 10^(-7), BSI.SimpleChangePt(change_pt_buffer))
    BSI.ChequeredBlocking(blocking_params..., P)
end

@testset "blocking object" begin
    obs = ℝ.([1.0, 1.2, 0.8, 1.3, 2.0])
    tt = [0.0, 1.0, 1.5, 2.3, 4.0]
    L = @SMatrix [1. 0.]
    ϵ = 10^(-7)
    Σ = @SMatrix [10^(-10)]

    𝔅 = blocking_test_prep(obs, tt)

    @testset "validity of initial set-up" begin
        @test 𝔅.idx == 1
        @test 𝔅.knots[1] == [1, 3]
        @test 𝔅.knots[2] == [2]
        @test 𝔅.blocks[1] == [[1], [2, 3], [4]]
        @test 𝔅.blocks[2] == [[1, 2], [3, 4]]
        @test 𝔅.change_pts[1] == [BSI.SimpleChangePt(100), BSI.NoChangePt(100), BSI.SimpleChangePt(100), BSI.NoChangePt(100)]
        @test 𝔅.change_pts[2] == [BSI.NoChangePt(100), BSI.SimpleChangePt(100), BSI.NoChangePt(100), BSI.NoChangePt(100)]
        @test 𝔅.vs == obs[2:end]
        @test 𝔅.Ls[1] == [I, L, I, L]
        @test 𝔅.Ls[2] == [L, I, L, L]
        @test 𝔅.Σs[1] == [I*ϵ, Σ, I*ϵ, Σ]
        @test 𝔅.Σs[2] == [Σ, I*ϵ, Σ, Σ]
    end

    θ = [10.0, -8.0, 15.0, 0.0, 3.0]
    # check function `next_set_of_blocks` (requires setting up Workspace)
    # NOTE it doesn't really check blocking object, it's more about Workspace
end


text:
context("Lines") do
    context("linesegment intersect") do
        a = LineSegment(Point2f0(0,0), Point2f0(1,0))
        b = LineSegment(Point2f0(0.5,0.5), Point2f0(0.5,-0.5))
        @fact intersects(a,b) --> (true, Point2f0(0.5,0.0))

        a = LineSegment(Point2f0(0,0), Point2f0(0.499,0))
        b = LineSegment(Point2f0(0.5,0.5), Point2f0(0.5,-0.5))
        @fact intersects(a,b) --> (false, zero(Point2f0))
    end
    context("poly self intersections") do
        p1 = Point2f0[(x, sin(x)) for x in 0:.5:6]
        p2 = Point2f0[(p1[i][1], 0) for i=length(p1):-1:1]
        points = vcat(p1, p2)
        inds, intersects = self_intersections(points)
        @fact inds --> [7,19]
        @fact intersects --> Point2f0[(3.1434426,0.0)]
        polys = split_intersections(points)
        @fact length(polys) --> 2
        @fact length(polys[1]) --> 12
        @fact length(polys[2]) --> 15

        u = linspace(0, 2pi, 200)
        points = Point2f0[(sin(x), sin(2x)) for x in u]

        inds, intersects = self_intersections(points)
        @fact inds --> [1,100]
        @fact intersects --> Point2f0[(0.0,0.0)]
        polys = split_intersections(points)
        @fact length(polys) --> 2
        @fact length(polys[1]) --> 99
        @fact length(polys[2]) --> 102
    end
end


text:
"""
This file is for creating the structure of Weight matrix W in the paper by using Flux.
This way, we may be able to use NN in the future work.
"""

#import Flux: functor, @functor
#using ForwardDiff, Zygote

"""
PredictionMat(in::Integer, out::Integer)

Creates a Matrix with parameter `W`.

y = W * x

The input `x` must be a vector of length `in`, or a batch of vectors represented
as an `in × N` matrix. The out `y` will be a vector or batch of length `out`.

```julia
julia> d = PredictionMat(5, 2)
PredictionMat(5, 2)

julia> d(rand(5))
Tracked 2-element Array{Float32,1}:
0.47773355f0
-0.8844353f0
"""
struct PredictionMat{S}
    W::S
end

"""
glorot_uniform borrowed from Flux.jl https://github.com/FluxML/Flux.jl/blob/master/src/utils.jl

duplicated to avoid adding Flux package to this project
"""
nfan() = 1, 1 # fan_in, fan_out
nfan(n) = 1, n # A vector is treated as a n×1 matrix
nfan(n_out, n_in) = n_in, n_out # In case of Dense kernels: arranged as matrices
nfan(dims...) = prod(dims[1:end-2]) .* (dims[end-1], dims[end]) # In case of convolution kernels
glorot_uniform(dims...) = (rand(Float32, dims...) .- 0.5f0) .* sqrt(24.0f0 / sum(nfan(dims...)))

function PredictionMat(in::Integer, out::Integer; initW = glorot_uniform)
    return PredictionMat(initW(out, in))
end

#@functor PredictionMat

function (a::PredictionMat)(x::AbstractArray)
    W = a.W
    W*x
end

function Base.show(io::IO, l::PredictionMat)
    print(io, "PredictionMat(", size(l.W, 2), ", ", size(l.W, 1))
    print(io, ")")
end

(a::PredictionMat{W})(x::AbstractArray{T}) where {T <: Union{Float32,Float64}, W <: AbstractArray{T}} =
invoke(a, Tuple{AbstractArray}, x)

(a::PredictionMat{W})(x::AbstractArray{<:AbstractFloat}) where {T <: Union{Float32,Float64}, W <: AbstractArray{T}} =
a(T.(x))


text:
# matrix and transformations

"""
    getmatrix()

Get the current matrix. Returns an array of six float64 numbers:

- xx component of the affine transformation

- yx component of the affine transformation

- xy component of the affine transformation

- yy component of the affine transformation

- x0 translation component of the affine transformation

- y0 translation component of the affine transformation

Some basic matrix transforms:

- translate

`transform([1, 0, 0, 1, dx, dy])` shifts by `dx`, `dy`

- scale

`transform([fx 0 0 fy 0 0])` scales by `fx` and `fy`

- rotate

`transform([cos(a), -sin(a), sin(a), cos(a), 0, 0])` rotates around to `a` radians

    rotate around O: [c -s s c 0 0]

- shear

`transform([1 0 a 1 0 0])` shears in x direction by `a`
-  shear in y: [1  B 0  1 0 0]

- x-skew

`transform([1, 0, tan(a), 1, 0, 0])` skews in x by `a`

- y-skew

`transform([1, tan(a), 0, 1, 0, 0])` skews in y by `a`

- flip

`transform([fx, 0, 0, fy, centerx * (1 - fx), centery * (fy-1)])` flips with center at `centerx`/`centery`

- reflect

`transform([1 0 0 -1 0 0])` reflects in xaxis

`transform([-1 0 0 1 0 0])` reflects in yaxis

When a drawing is first created, the matrix looks like this:

    getmatrix() = [1.0, 0.0, 0.0, 1.0, 0.0, 0.0]

When the origin is moved to 400/400, it looks like this:

    getmatrix() = [1.0, 0.0, 0.0, 1.0, 400.0, 400.0]

To reset the matrix to the original:

    setmatrix([1.0, 0.0, 0.0, 1.0, 0.0, 0.0])

"""
function getmatrix()
    gm = Cairo.get_matrix(get_current_cr())
    return([gm.xx, gm.yx, gm.xy, gm.yy, gm.x0, gm.y0])
end

"""
    setmatrix(m::Array)

Change the current matrix to matrix `m`. Use `getmatrix()` to get the current matrix.
"""
function setmatrix(m::Array)
    if eltype(m) != Float64
        m = map(Float64, m)
    end
    # some matrices make Cairo freak out and need reset. Not sure what the rules are yet…
    if length(m) < 6
        throw("didn't like that matrix $m: not enough values")
    elseif count(!iszero, m) == 0
        throw("didn't like that matrix $m: too many zeroes")
    else
        cm = Cairo.CairoMatrix(m[1], m[2], m[3], m[4], m[5], m[6])
        Cairo.set_matrix(get_current_cr(), cm)
    end
end

"""
    transform(a::Array)

Modify the current matrix by multiplying it by matrix `a`.

For example, to skew the current state by 45 degrees in x and move by 20 in y direction:

    transform([1, 0, tand(45), 1, 0, 20])

Use `getmatrix()` to get the current matrix.
"""
function transform(a::Array)
    b = Cairo.get_matrix(get_current_cr())
    setmatrix([
        (a[1] * b.xx)  + a[2]  * b.xy,             # xx
        (a[1] * b.yx)  + a[2]  * b.yy,             # yx
        (a[3] * b.xx)  + a[4]  * b.xy,             # xy
        (a[3] * b.yx)  + a[4]  * b.yy,             # yy
        (a[5] * b.xx)  + (a[6] * b.xy) + b.x0,     # x0
        (a[5] * b.yx)  + (a[6] * b.yy) + b.y0      # y0
    ])
end

"""
    rotationmatrix(a)

Return a 3x3 Julia matrix that will apply a rotation through `a` radians.
"""
function rotationmatrix(a)
    return ([cos(a)  -sin(a)  0.0 ;
             sin(a)   cos(a)  0.0 ;
             0.0         0.0  1.0 ])
end

"""
    translationmatrix(x, y)

Return a 3x3 Julia matrix that will apply a translation in `x` and `y`.
"""
function translationmatrix(x, y)
    return ([1.0   0.0   x ;
             0.0   1.0   y ;
             0.0   0.0   1.0 ])
end

"""
    scalingmatrix(sx, sy)

Return a 3x3 Julia matrix that will apply a scaling by `sx` and `sy`.
"""
function scalingmatrix(sx, sy)
    return ([sx 0.0  0.0 ;
           0.0   sy  0.0 ;
           0.0  0.0  1.0])
end

"""
    cairotojuliamatrix(c)

Return a 3x3 Julia matrix that's the equivalent of the six-element matrix in `c`.
"""
function cairotojuliamatrix(c::Array)
    return [c[1] c[3] c[5] ; c[2] c[4] c[6] ; 0.0 0.0 1.0]
end

"""
    juliatocairomatrix(c)

Return a six-element matrix that's the equivalent of the 3x3 Julia matrix in `c`.
"""
function juliatocairomatrix(c::Matrix)
    return [c[1] c[2] c[4] c[5] c[7] c[8]]
end

"""
    getrotation(R::Matrix)
    getrotation()

Get the rotation of a Julia 3x3 matrix, or the current Luxor rotation.

```math
\\begin{bmatrix}
a & b & tx \\\\
c & d & ty \\\\
0 & 0 & 1  \\\\
\\end{bmatrix}
```

The rotation angle is `atan(-b, a)` or `atan(c, d)`.
"""
function getrotation(R::Matrix)
    # t = atan(-R[4], R[1]) # should be the same as:
    t = atan(R[2], R[5])
    return mod2pi(t)
end

function getrotation()
    getrotation(cairotojuliamatrix(getmatrix()))
end

"""
    getscale(R::Matrix)
    getscale()

Get the current scale of a Julia 3x3 matrix, or the current Luxor scale.

Returns a tuple of x and y values.
"""
function getscale(R::Matrix)
    sx = hypot(R[1], R[2])
    sy = hypot(R[4], R[5])
    return (sx, sy)
end

function getscale()
    getscale(cairotojuliamatrix(getmatrix()))
end

"""
    gettranslation(R::Matrix)
    gettranslation()

Get the current translation of a Julia 3x3 matrix, or the current Luxor translation.

Returns a tuple of x and y values.
"""
function gettranslation(R::Matrix)
    return (R[7], R[8])
end

function gettranslation()
    gettranslation(cairotojuliamatrix(getmatrix()))
end


text:
export girth, girth_cycle

"""
`girth_cycle(G)` returns a shorest cycle of `G` as an array listing
the vertices on that cycle, or an empty array if `G` is acyclic.

*Warning*: This implementation is quite inefficient.
"""
function girth_cycle(G::SimpleGraph{T}) where {T}
    if cache_check(G, :girth_cycle)
        return cache_recall(G, :girth_cycle)
    end
    best_path = T[]
    if is_acyclic(G)
        cache_save(G, :girth_cycle, best_path)
        return best_path
    end
    best = NV(G) + 1

    GG = deepcopy(G)

    for e in elist(GG)
        u, v = e
        SimpleGraphs.delete!(GG, u, v)
        P = find_path(GG, u, v)
        add!(GG, u, v)
        nP = length(P)
        if 0 < nP < best
            best = nP
            best_path = P
        end
    end
    cache_save(G, :girth_cycle, best_path)
    return best_path
end

"""
`girth(G)` computes the length of a shortest cycle in `G` or returns `0`
if `G` is acyclic.

**Warning**: This implementation is quite inefficient.
"""
function girth(G::SimpleGraph)
    if cache_check(G, :girth)
        return cache_recall(G, :girth)
    end
    g = length(girth_cycle(G))
    cache_save(G, :girth, g)
    return g
end


text:
using CairoMakie
using Distributed
using Distributions

proc_number = 4;

if nprocs() < proc_number
        addprocs(proc_number - nprocs())
end

@everywhere include("Calculation/calc_settings.jl")
# x-axis of dI/dV line plot in Bohr radii
X = d1[1] .* xs
## Plotting
#To plot simulated dI/dV line cut along with experimental data for charged SV
function spectral_cut()
        res = @showprogress pmap(x -> spectral_bulk(ω, Location(x, 0), s, xs)) #Computation, Get the results
        fig = Figure(resolution = (800, 300))
        ax = fig[1, 1] = Axis(
                fig,
                xlabel = "x/nm",
                ylabel = "Spectral function(A(x))",
                title = "Spectral cut",
                xlabelpadding = 0,
                ylabelpadding = 0,
                xlabelsize = 16,
                ylabelsize = 16,
                xticklabelsize = 14,
                yticklabelsize = 14,
                titlesize = 18,
                xticklabelfont = "serif-roman",
                yticklabelfont = "serif-roman",
                xlabelfont = "serif-italic",
                ylabelfont = "serif-roman",
        )
        qft_data = lines!(X_charged * a0 / 10, res, label = "A single local potential, ω = 0.60eV")
        axislegend()
        fig
end

@time spectral_cut()


text:
using Documenter, ShiftedArrays

DocMeta.setdocmeta!(ShiftedArrays, :DocTestSetup, :(using ShiftedArrays); recursive=true)

makedocs(
    # options
    modules = [ShiftedArrays],
    sitename = "ShiftedArrays.jl",
    format = Documenter.HTML(),
    pages = Any[
        "Introduction" => "index.md",
        "API" => "api.md",
    ]
)

# Deploy built documentation from Travis.
# =======================================

deploydocs(
    # options
    repo = "github.com/piever/ShiftedArrays.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)