text:
# ------------------------------------------------------------------------------------------
# # Creating a dynamical system
#
# Topics:
# * What *is* a dynamical system in the context of **DynamicalSystems.jl**?
# * How can we define a discrete or continuous system?
# * What are the differences in handling large and small systems?
#     * Handy Dandy definition table!
# * Comment on DifferentialEquations.jl
# * Docstrings
#
# ---
#
# # Dynamical System
#
# A "dynamical system" is a law that describes how some variables should evolve in time.
# This law is described by the "equations of motion" function $\vec{f}$.
#
# There are two types of dynamical systems (in our case):
#
# 1. One is called a map, where time is a discrete quantity (like "steps" or "generations").
# The equations of motion then look something like $\vec{u}_{n+1} = \vec{f}(\vec{u}_n, p,
# n)$ where $n$ is an integer and $p$ are the parameters of the system.
#
# 2. The other type is called an Ordinary Differential Equation (ODE), where time is a
# continuous quantity. Then, the equations of motion give the time derivatives of the
# variables: $\frac{d\vec{u}}{dt} = \vec{f}(\vec{u}, p, t)$.
#
# In both cases $\vec{u}$ is the state of the system, a vector of the variables that define
# the system. For example, $\vec{u}$ may be defined as $(x, y, z)$, or as $(\theta, I)$, or
# as $(V, I, R, C_1, C_2)$, etc.
#
# ---
#
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Creating a simple discrete system
#
#
# For all intents and purposes, a `DynamicalSystem` is a Julia `struct` that contains all
# the information of a dynamical system:
# 1. Equations of motion function $\vec{f}$
# 3. Initial state.
# 4. Parameter container (if the equations of motion are parameterized).
# 2. *Optional* : Jacobian of the equations of motion.
#
#
# For simplicity let's focus on creating a simple discrete chaotic system, the [Hénon
# map](https://en.wikipedia.org/wiki/H%C3%A9non_map).
#
# The constructor we need is:
#
# ```julia
# DiscreteDynamicalSystem(eom, state, p)
# ```
#
# In order to construct a Hénon map, let's discuss the input arguments we need to pass to
# this constructor!
#
# #### 1. Equations of motion
# The first argument of the above constructor --`eom` -- is a **function** representing the
# equations of motion.
#
# Therefore, the first step in creating a `DynamicalSystem` is to define the equations of
# motion. For the Hénon map they have the form:
#
# $$
# \begin{aligned}
# x_{n+1} &= 1 - ax^2_n+y_n \\
# y_{n+1} & = bx_n
# \end{aligned}
# $$
#
# Now we must make a function out of them. There is some strictness when writing this
# function in Julia. Specifically, it can only be one of two forms: out-of-place (oop) or
# in-place (iip). Here is how to define it:
#
# * **oop** : The `eom` function **must** be in the form `eom(x, p, t) -> SVector`
#   which means that given a state `x::SVector` and some parameter container
#   `p` it returns a
# [`SVector`](http://juliaarrays.github.io/StaticArrays.jl/stable/pages/api.html#SVector-1)
#   (from the [StaticArrays](https://github.com/JuliaArrays/StaticArrays.jl) module)
#   containing the next state.
# * **iip** : The `eom` function **must** be in the form `eom!(xnew, x, p, t)`
#   which means that given a state `x::Vector` and some parameter container `p`,
#   it writes in-place the new state in `xnew`.
#
# We advise to use **oop** for systems with dimension < 11, and use **iip** otherwise.
#
# *If you are familiar with [DifferentialEquations.jl](http://docs.juliadiffeq.org/latest/),
# then notice that the equations of motion are defined in an identical manner*
#
# Because the Hénon map is only 2-dimensional, we follow the advice of the documentation and
# take advantage of the [`StaticArrays`](https://github.com/JuliaArrays/StaticArrays.jl)
# module.
# ------------------------------------------------------------------------------------------

using DynamicalSystems
h_eom(x, p, t) = SVector{2}(1.0 - p[1]*x[1]^2 + x[2], p[2]*x[1])

# ------------------------------------------------------------------------------------------
# * Remember: *both* `p` (for parameters) and `t` (for time) must be included in the
# equations of motion function, irrespective of whether they are used or not!
#
# #### 2. State
#
# The second argument for the `DiscreteDynamicalSystem` constructor is a `state`, which
# represents the initial condition for the system:
# ------------------------------------------------------------------------------------------

state = zeros(2) # doesn't matter if we use Vector of SVector for the `state`

# ------------------------------------------------------------------------------------------
# #### 3. Parameters for our EOM
#
# Then, the last argument, `p`, is simply a parameter container.
# ------------------------------------------------------------------------------------------

p = [1.4, 0.3] # p = [a, b] from the equations of motion

# ------------------------------------------------------------------------------------------
# These 3 things are enough to make a `DynamicalSystem`:
# ------------------------------------------------------------------------------------------

henon = DiscreteDynamicalSystem(h_eom, state, p)

# ------------------------------------------------------------------------------------------
# ## Getting a trajectory
#
# Now that we have created the system, the first (and most basic) thing to do is to simply
# plot its time evolution and see what it looks like.
#
# The `trajectory` function is a convenient tool that evolves the system and returns the
# output at equally spaced time intervals. The call signature is simply:
#
# ```julia
# trajectory(ds::DynamicalSystem, T [, u]; kwargs...)
# ```
# which evolves a system for total time `T`, optionally starting from a different state `u`.
# ------------------------------------------------------------------------------------------

# trajectory from initial condition
tr = trajectory(henon, 100000)

# trajectory from a different starting point
tr2 = trajectory(henon, 100000, 0.01rand(2))

using PyPlot
figure(figsize=(6,4))
plot(tr[:, 1], tr[:, 2], lw = 0.0, marker = "o", ms = 0.1, alpha = 0.5);
plot(tr2[:, 1], tr2[:, 2], lw = 0.0, marker = "o", ms = 0.1, alpha = 0.5);
xlabel("x"); ylabel("y");

# ------------------------------------------------------------------------------------------
# ## Crash-course on `Dataset`
#
# `trajectory` returns an object that is a `Dataset`:
# ------------------------------------------------------------------------------------------

tr = trajectory(henon, 100000)
println(typeof(tr))
println(summary(tr))

# ------------------------------------------------------------------------------------------
# `Dataset` instances handle most data in **DynamicalSystems.jl**. `Dataset` is a wrapper of
# a `Vector` of `SVector`s (statically sized vectors).
# ------------------------------------------------------------------------------------------

typeof(tr.data) # this is where the `Dataset` has the data

# ------------------------------------------------------------------------------------------
# When accessed with one index, a `Dataset` behaves as a vector of vectors
# ------------------------------------------------------------------------------------------

tr[1] # variables at first time point

tr[1:5]

# ------------------------------------------------------------------------------------------
# When accessed with two indices, a `Dataset` behaves like a matrix with each column being
# the timeseries of each dynamic variable
# ------------------------------------------------------------------------------------------

tr[:, 1] # timeseries of first variable

tr[1:56, 2] # time points 1:56 of second variable

tr[12, 1] # value of first variable at 12th timepoint

tr[1:10, 1:2] # using two ranges returns a `Dataset`

# ------------------------------------------------------------------------------------------
# # Adding a Jacobian
#
# Let's see `henon` again:
# ------------------------------------------------------------------------------------------

henon

# ------------------------------------------------------------------------------------------
# * The last line, "`jacobian:  ForwardDiff`"
# says that the Jacobian function of the equations of motion was computed automatically
# using the module
# [`ForwardDiff`](http://www.juliadiff.org/ForwardDiff.jl/stable/user/api.html).
# * The Jacobian function is a crucial component of a dynamical system, and that is why if
# it is not given, it is computed automatically.
#
#
#
# Even though the automatic computation is very efficient, the best possible performance
# will come if you pass a "hard-coded" jacobian:
# ------------------------------------------------------------------------------------------

h_jacobian(x, p, t) = @SMatrix [-2*p[1]*x[1] 1.0; p[2] 0.0]

# ------------------------------------------------------------------------------------------
# * Notice that for **out-of-place** systems, the Jacobian must also have the same form as
# the equations of motion, namely to return an `SMatrix`.
#
# Now, we can pass this Jacobian function to the `DiscreteDynamicalSystem` constructor as a
# 4th argument:
# ------------------------------------------------------------------------------------------

henon_with_jac = DiscreteDynamicalSystem(h_eom, state, p, h_jacobian)

# ------------------------------------------------------------------------------------------
# To see the difference in performance, let's call this Jacobian function
# ------------------------------------------------------------------------------------------

using BenchmarkTools
a = rand(SVector{2})
@btime $(henon_with_jac.jacobian)($a, $henon.prob.p, 0);
@btime $(henon.jacobian)($a, $henon.prob.p, 0);

# ------------------------------------------------------------------------------------------
# You can see that even though `ForwardDiff` is truly performant, the hard-coded version is
# much faster.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Continuous System
# * The process of creating a continuous system is identical to that of a discrete system,
# except that the constructor `ContinuousDynamicalSystem` is used instead.
# * In this case the function `eom` returns the time derivatives and not a "next state".
#
#
# ---
#
#
#
# * We will take the opportunity to show the process of using in-place equations of motion
# for a continuous system, which is aimed to be used for large systems (dimensionality $\ge$
# 10).
#
# * In addition, the system we will use (Henon-Heiles) does not have any parameters.
# $$
# \begin{aligned}
# \dot{x} &= p_x \\
# \dot{y} &= p_y \\
# \dot{p}_x &= -x -2 xy \\
# \dot{p}_y &= -y - (x^2 - y^2)
# \end{aligned}
# $$
# ------------------------------------------------------------------------------------------

# Henon-heiles famous system
# in-place form of equations of motion
# du is the derivatives vector, u is the state vector
function hheom!(du, u, p, t)
    du[1] = u[3]
    du[2] = u[4]
    du[3] = -u[1] - 2u[1]*u[2]
    du[4] = -u[2] - (u[1]^2 - u[2]^2)
    return nothing
end

# pass `nothing` as the parameters, because the system doesn't have any
hh = ContinuousDynamicalSystem(hheom!, [0, -0.25, 0.42081, 0], nothing)

# ------------------------------------------------------------------------------------------
# Great, now we can get a trajectory of this system, by solving using
# DifferentialEquations.jl
# ------------------------------------------------------------------------------------------

tr = trajectory(hh, 100.0, dt = 0.05)

figure(figsize = (6,4))
plot(tr[:, 1], tr[:, 2]);
xlabel("\$q_1\$"); ylabel("\$q_2\$");

# ------------------------------------------------------------------------------------------
# # Handy Dandy definition table
#
# Depending on whether your system is small or large, you want to use out-of-place or in-
# place equations of motion. The Jacobian function (if you wish to provide it) must also be
# of the same form.
#
# Here is a handy table that summarizes what the definitions should look like:
#
# |          System Type         |    equations of motion    |            Jacobian
# |
# |:----------------------------:|:-------------------------:|:-----------------------------
# -:|
# | in-place (big systems)       | `eom!(du, u, p, t)`       | `jacobian!(J, u, p, t)`
# |
# | out-of-place (small systems) | `eom(u, p, t) -> SVector` | `jacobian(u, p, t) ->
# SMatrix` |
#
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Time-evolution of systems & DifferentialEquations.jl
# Discrete systems are evolved using internal algorithms. However, all time evolution of
# continuous systems is done through the
# [DifferentialEquations.jl](http://docs.juliadiffeq.org/latest/) library. In fact,
# `trajectory` for continuous systems simply wraps `solve` with some extra arguments.
#
# Keep in mind that by default all continuous systems are evolved using the solver `Vern9()`
# (9th order Verner solver) with tolerances `abstol = reltol = 1e-9`. This is especially
# important when one wants to compare benchmarks with different packages.
#
# It is almost certain that if you use **DynamicalSystems.jl** you want to use also
# DifferentialEquations.jl, due to the huge list of available features.
#
#
# ## When to use **DynamicalSystems.jl**?
# At this point in the tutorial you might be thinking:
# > How is DynamicalSystems.jl different from DifferentialEquations.jl? Seems the same to me
# so far...
#
# That's because we haven't seen any special features of **DynamicalSystems.jl** yet!
# `trajectory` is nothing more than a "convenient shortcut".
#
# The [contents page](https://juliadynamics.github.io/DynamicalSystems.jl/latest/#contents)
# of the documentation does a good job describing what is possible through
# **DynamicalSystems.jl**.
#
# In addition, for usage with DifferentialEquations.jl, we can create specialized
# integrators that evolve the system and the tangent space, or many states of the system in
# parallel (at *exactly* same times):
# * `tangent_integrator`
# * `parallel_integrator`
#
# These functions work for both continuous and discrete systems and also work regardless of
# whether the system is in-place (large), out-of-place (small) or auto-differentiated.
# Special attention has also been given to the performance of the integrators.
#
# Also, keep in mind that regardless of whether you use a `DynamicalSystem` to produce a
# timeseries or not, more than half of the **DynamicalSystems.jl** library is using
# numerical data as an input. Most of the time this numerical data is expected in the form
# of an `AbstractDataset` instance, which is what `trajectory` returns.
#
# ## Using DifferentialEquations.jl from a `DynamicalSystem`
#
# It is *very* likely that other features of DifferentialEquations.jl will be useful to
# someone using **DynamicalSystems.jl**. However, you can still use DifferentialEquations.jl
# *after* you have defined a continuous dynamical system, because the field `prob` gives an
# `ODEProblem`:
# ------------------------------------------------------------------------------------------

hh.prob

# ------------------------------------------------------------------------------------------
# * **please be careful when using this problem directly, because as you can see the
# `tspan[end]` field is `Inf`!!!**
#
#
# * One final comment: using Callbacks is *not* possible with **DynamicalSystems.jl**,
# because the the equations of motion function has to be assumed differentiable
# "everywhere".
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Docstrings
# ------------------------------------------------------------------------------------------

?DynamicalSystem

?trajectory

?Dataset

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/introduction-to-dynamicalsystems.jl/1. Creating a dynamical system.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using BinDeps
using Compat.Libdl

@BinDeps.setup

fmmlib2d = library_dependency("fmmlib2d")

provides(Sources, URI("https://github.com/ludvigak/fmmlib2d/archive/master.zip"),
         fmmlib2d, unpacked_dir = "fmmlib2d-master")

fmmsrcdir = joinpath(BinDeps.srcdir(fmmlib2d), "fmmlib2d-master/src")
libname = "fmmlib2d." * Libdl.dlext
libfile = joinpath(BinDeps.libdir(fmmlib2d),libname)
buildfile = joinpath(fmmsrcdir, libname)

fsrcs = "hfmm2dpart.f hfmm2drouts.f d2tstrcr_omp.f d2mtreeplot.f h2dterms.f helmrouts2d.f cdjseval2d.f hank103.f prini.f cfmm2dpart.f zfmm2dpart.f rfmm2dpart.f lfmm2dpart.f lfmm2drouts.f l2dterms.f laprouts2d.f"

fsrcs = split(fsrcs, " ")
buildcmd = `gfortran -O2 -fPIC -shared $fsrcs -o $libname`

provides(BuildProcess,
         (@build_steps begin
          GetSources(fmmlib2d)
          @build_steps begin
          ChangeDirectory(fmmsrcdir)
          FileRule(libfile, @build_steps begin
                   buildcmd
                   CreateDirectory(libdir(fmmlib2d))
                   `cp $libname $libfile`
                   end)
          end
          end),
         fmmlib2d)

@BinDeps.install Dict(:fmmlib2d => :fmmlib2d)

ext: jl


text:
abstract type Alignable end

"""
Used to specify space that is occupied in a grid. Like 1:1|1:1 for the first square,
or 2:3|1:4 for a rect over the 2nd and 3rd row and the first four columns.
"""
struct Span
    rows::UnitRange{Int64}
    cols::UnitRange{Int64}
end

"""
An object that can be aligned that also specifies how much space it occupies in
a grid via its span.
"""
struct SpannedAlignable{T <: Alignable}
    al::T
    sp::Span
end

"""
    side_indices(c::SpannedAlignable)::RowCols{Int}

Indices of the rows / cols for each side
"""
function side_indices(c::SpannedAlignable)
    return RowCols(
        c.sp.cols.start,
        c.sp.cols.stop,
        c.sp.rows.start,
        c.sp.rows.stop,
    )
end

"""
These functions tell whether an object in a grid touches the left, top, etc. border
of the grid. This means that it is relevant for the grid's own protrusion on that side.
"""
ismostin(sp::SpannedAlignable, grid, ::Left) = sp.sp.cols.start == 1
ismostin(sp::SpannedAlignable, grid, ::Right) = sp.sp.cols.stop == grid.ncols
ismostin(sp::SpannedAlignable, grid, ::Bottom) = sp.sp.rows.stop == grid.nrows
ismostin(sp::SpannedAlignable, grid, ::Top) = sp.sp.cols.start == 1

isleftmostin(sp::SpannedAlignable, grid) = ismostin(sp, grid, Left())
isrightmostin(sp::SpannedAlignable, grid) = ismostin(sp, grid, Right())
isbottommostin(sp::SpannedAlignable, grid) = ismostin(sp, grid, Bottom())
istopmostin(sp::SpannedAlignable, grid) = ismostin(sp, grid, Top())

struct SolvedAxisLayout <: Alignable
    inner::BBox
    outer::BBox
    axis::LayoutedAxis
end

struct AxisLayout <: Alignable
    decorations::BBox
    axis::LayoutedAxis
end

struct SolvedGridLayout <: Alignable
    bbox::BBox
    content::Vector{SpannedAlignable}
    nrows::Int
    ncols::Int
    grid::RowCols{Vector{Float64}}
end

struct GridLayout <: Alignable
    content::Vector{SpannedAlignable}
    nrows::Int
    ncols::Int
    colratios::Vector{Float64}
    rowratios::Vector{Float64}
    colgapfraction::Float64
    rowgapfraction::Float64
end

"""
All the protrusion functions calculate how much stuff "sticks out" of a layoutable object.
This is so that collisions are avoided, while what is actually aligned is the
"important" edges of the layout objects.
"""
leftprotrusion(x) = protrusion(x, Left())
rightprotrusion(x) = protrusion(x, Right())
bottomprotrusion(x) = protrusion(x, Bottom())
topprotrusion(x) = protrusion(x, Top())

protrusion(u::AxisLayout, side::Side) = u.decorations[side]
protrusion(sp::SpannedAlignable, side::Side) = protrusion(sp.al, side)

function protrusion(gl::GridLayout, side::Side)
    return mapreduce(max, gl.content, init = 0.0) do elem
        # we use only objects that stick out on this side
        # And from those we use the maximum protrusion
        ismostin(elem, gl, side) ? protrusion(elem, side) : 0.0
    end
end

protrusion(s::SolvedAxisLayout, ::Left) = left(s.inner) - left(s.outer)
function protrusion(s::SolvedAxisLayout, side::Side)
    return s.outer[side] - s.inner[side]
end

"""
This function solves a grid layout such that the "important lines" fit exactly
into a given bounding box. This means that the protrusions of all objects inside
the grid are not taken into account. This is needed if the grid is itself placed
inside another grid.
"""
function solve(gl::GridLayout, bbox::BBox)

    # first determine how big the protrusions on each side of all columns and rows are
    maxgrid = RowCols(gl.ncols, gl.nrows)
    # go through all the layout objects placed in the grid
    for c in gl.content
        idx_rect = side_indices(c)
        mapsides(idx_rect, maxgrid) do side, idx, grid
            grid[idx] = max(grid[idx], protrusion(c.al, side))
        end
    end
    # compute what size the gaps between rows and columns need to be
    colgaps = maxgrid.lefts[2:end] .+ maxgrid.rights[1:end-1]
    rowgaps = maxgrid.tops[2:end] .+ maxgrid.bottoms[1:end-1]

    # determine the biggest gap
    # using the biggest gap size for all gaps will make the layout more even, but one
    # could make this aspect customizable, because it might waste space
    maxcolgap = maximum(colgaps)
    maxrowgap = maximum(rowgaps)

    # determine the vertical and horizontal space needed just for the gaps
    # again, the gaps are what the protrusions stick into, so they are not actually "empty"
    # depending on what sticks out of the plots
    sumcolgaps = maxcolgap * (gl.ncols - 1)
    sumrowgaps = maxrowgap * (gl.nrows - 1)

    # compute what space remains for the inner parts of the plots
    remaininghorizontalspace = width(bbox) - sumcolgaps
    remainingverticalspace = height(bbox) - sumrowgaps

    # compute how much gap to add, in case e.g. labels are too close together
    # this is given as a fraction of the space used for the inner parts of the plots
    # so far, but maybe this should just be an absolute pixel value so it doesn't change
    # when resizing the window
    addedcolgap = gl.colgapfraction * remaininghorizontalspace
    addedrowgap = gl.rowgapfraction * remainingverticalspace

    # compute the actual space available for the rows and columns (plots without protrusions)
    spaceforcolumns = remaininghorizontalspace - addedcolgap * (gl.ncols - 1)
    spaceforrows = remainingverticalspace - addedrowgap * (gl.nrows - 1)

    # compute the column widths and row heights using the specified row and column ratios
    colwidths = gl.colratios ./ sum(gl.colratios) .* spaceforcolumns
    rowheights = gl.rowratios ./ sum(gl.rowratios) .* spaceforrows

    # this is the vertical / horizontal space between the inner lines of all plots
    colgap = maxcolgap + addedcolgap
    rowgap = maxrowgap + addedrowgap

    # compute the x values for all left and right column boundaries
    xleftcols = [left(bbox) + sum(colwidths[1:i-1]) + (i - 1) * colgap for i in 1:gl.ncols]
    xrightcols = xleftcols .+ colwidths

    # compute the y values for all top and bottom row boundaries
    ytoprows = [top(bbox) - sum(rowheights[1:i-1]) - (i - 1) * rowgap for i in 1:gl.nrows]
    ybottomrows = ytoprows .- rowheights

    # now we can solve the content thats inside the grid because we know where each
    # column and row is placed, how wide it is, etc.
    # note that what we did at the top was determine the protrusions of all grid content,
    # but we know the protrusions before we know how much space each plot actually has
    # because the protrusions should be static (like tick labels etc don't change size with the plot)

    gridboxes = RowCols(
        xleftcols, xrightcols,
        ytoprows, ybottomrows
    )
    solvedcontent = map(gl.content) do c
        idx_rect = side_indices(c)
        bbox_cell = mapsides(idx_rect, gridboxes) do side, idx, gridside
            gridside[idx]
        end
        solved = solve(c.al, bbox_cell)
        return SpannedAlignable(solved, c.sp)
    end
    # return a solved grid layout in which all objects are also solved layout objects
    return SolvedGridLayout(
        bbox, solvedcontent,
        gl.nrows, gl.ncols,
        gridboxes
    )
end

"""
This function solves a grid layout so that it fits exactly inside a bounding box.
Exactly means that the protrusions of all other objects inside this grid layout
also have to fit into the bounding box. This is needed if the grid is the outermost
object in the layout, the bounding box would then be the scene boundary.
"""
function outersolve(gl::GridLayout, bbox::BBox)
    maxgrid = RowCols(gl.ncols, gl.nrows)
    for c in gl.content
        idx_rect = side_indices(c)
        mapsides(idx_rect, maxgrid) do side, idx, grid
            grid[idx] = max(grid[idx], protrusion(c.al, side))
        end
    end

    topprot = maxgrid.tops[1]
    bottomprot = maxgrid.bottoms[end]
    leftprot = maxgrid.lefts[1]
    rightprot = maxgrid.rights[end]

    colgaps = maxgrid.lefts[2:end] .+ maxgrid.rights[1:end-1]
    rowgaps = maxgrid.tops[2:end] .+ maxgrid.bottoms[1:end-1]

    maxcolgap = gl.ncols <= 1 ? 0 : maximum(colgaps)
    maxrowgap = gl.nrows <= 1 ? 0 : maximum(rowgaps)

    sumcolgaps = maxcolgap * (gl.ncols - 1)
    sumrowgaps = maxrowgap * (gl.nrows - 1)

    remaininghorizontalspace = width(bbox) - sumcolgaps - leftprot - rightprot
    remainingverticalspace = height(bbox) - sumrowgaps - topprot - bottomprot

    addedcolgap = gl.colgapfraction * remaininghorizontalspace
    addedrowgap = gl.rowgapfraction * remainingverticalspace

    spaceforcolumns = remaininghorizontalspace - addedcolgap * (gl.ncols - 1)
    spaceforrows = remainingverticalspace - addedrowgap * (gl.nrows - 1)

    colwidths = gl.colratios ./ sum(gl.colratios) .* spaceforcolumns
    rowheights = gl.rowratios ./ sum(gl.rowratios) .* spaceforrows

    colgap = maxcolgap + addedcolgap
    rowgap = maxrowgap + addedrowgap

    xleftcols = map(1:gl.ncols) do i
        left(bbox) + leftprot + sum(colwidths[1:i-1]) + (i - 1) * colgap
    end
    xrightcols = xleftcols .+ colwidths

    ytoprows = map(1:gl.nrows) do i
        top(bbox) - topprot - sum(rowheights[1:i-1]) - (i - 1) * rowgap
    end

    ybottomrows = ytoprows .- rowheights
    gridboxes = RowCols(
        xleftcols, xrightcols,
        ytoprows, ybottomrows
    )
    solvedcontent = map(gl.content) do c
        idx_rect = side_indices(c)
        bbox_cell = mapsides(idx_rect, gridboxes) do side, idx, gridside
            gridside[idx]
        end
        solved = solve(c.al, bbox_cell)
        return SpannedAlignable(solved, c.sp)
    end
    # solvedcontent = solve.(gl.content)
    SolvedGridLayout(
        bbox, solvedcontent, gl.nrows, gl.ncols,
        gridboxes
    )
end

function solve(ua::AxisLayout, innerbbox)
    bbox = mapsides(innerbbox, ua.decorations) do side, iside, decside
        op = side isa Union{Left, Top} ? (-) : (+)
        return op(iside, decside)
    end
    SolvedAxisLayout(innerbbox, bbox, ua.axis)
end

const Indexables = Union{UnitRange, Int, Colon}

"""
This function allows indexing syntax to add a layout object to a grid.
You can do:

grid[1, 1] = obj
grid[1, :] = obj
grid[1:3, 2:5] = obj

and all combinations of the above
"""
function Base.setindex!(g, a::Alignable, rows::Indexables, cols::Indexables)
    if rows isa Int
        rows = rows:rows
    elseif rows isa Colon
        rows = 1:g.nrows
    end
    if cols isa Int
        cols = cols:cols
    elseif cols isa Colon
        cols = 1:g.ncols
    end

    if !((1 <= rows.start <= g.nrows) || (1 <= rows.stop <= g.nrows))
        error("invalid row span $rows for grid with $(g.nrows) rows")
    end
    if !((1 <= cols.start <= g.ncols) || (1 <= cols.stop <= g.ncols))
        error("invalid col span $cols for grid with $(g.ncols) columns")
    end

    push!(g.content, SpannedAlignable(a, Span(rows, cols)))
end


text:
using OrdinaryDiffEq
using GLMakie

#=
 # 1. simplest case

 ̈x =  -sin(x) - λ ̇x

=#
dampedpendulum(du, u, p, t) = @. -sin(u) - p.λ * du

du0 = [0.]
u0  = [1.0]
tspan = (0.0, 100.0)
p = (; λ = 0.1)

prob = SecondOrderODEProblem(dampedpendulum, du0, u0, tspan, p)
sol = solve(prob, VelocityVerlet())

# other solvers, e.g. DPRKN6

begin
    f = Figure()

    ax = Axis(f[1,1], title = "Position")
    lines!(sol.t, sol[1,:])

    ax = Axis(f[1,2], title = "Velocity")
    lines!(sol.t, sol[2,:])

    Label(f[0,1:2], "Damped pendulum", textsize = 24)

    mkpath("plots")
    save("plots/01_damped_pendulum.png", f)
    f
end

max_stars_repo_name: SteffenPL/Julia-for-mathematicians

max_issues_repo_name: SteffenPL/Julia-for-mathematicians

max_forks_repo_name: SteffenPL/Julia-for-mathematicians


text:
module CorrelationFunctions
using JSON: JSON, parse
import Images

include("utility.jl")

module Directional
using StatsBase: fit, Histogram, mean
using LinearAlgebra: normalize, norm
using Images: Kernel, imgradients, feature_transform,
    distance_transform, label_components
using PrettyTables: pretty_table
using Base.Iterators
using FFTW: plan_rfft, plan_irfft, fft
using CircularArrays: CircularArray
import ..Plane, ..Torus, ..Topology

include("directional/directions.jl")
include("directional/corrdata.jl")
include("directional/slicer.jl")
include("directional/indicator.jl")
include("directional/l2.jl")
include("directional/s2.jl")
include("directional/c2.jl")
include("directional/surface.jl")
include("directional/chord-length.jl")
# FIXME: not actually directional, but uses the same set of functions
include("directional/pore-size.jl")

export l2, s2, c2,
    surfsurf, surfvoid, lowfreq_energy_ratio, chord_length, pore_size,
    AbstractIndicator, SeparableIndicator, InseparableIndicator, S2FTPlans,
    CorrelationData, default_directions, directions, unit_length

# These are exported for their docstrings.
export direction1Dp, direction2Dp, direction3Dp
end # Directional

module Map
using LinearAlgebra: norm
using CUDA: CuArray, cu, @cuda, blockIdx, blockDim, threadIdx, gridDim
using CUDA.CUFFT: fft!, ifft!
using Images: imgradients, KernelFactors, label_components
using Interpolations: interpolate, Gridded, Linear, extrapolate, Periodic
import ..Plane, ..Torus

include("map/result.jl")
include("map/general_map.jl")
include("map/algorithms.jl")
include("map/iterators.jl")
include("map/l2_map.jl")
include("map/s2_map.jl")
include("map/c2_map.jl")
include("map/ss_map.jl")
include("map/sv_map.jl")

export l2, s2, c2, surfsurf, surfvoid,
    dir_from_map, restore_full_map, mean_dir
end # Map

using .Directional
export read_cuboid, pore_size, lowfreq_energy_ratio,
    Directional, Map
end # module


text:
register(DataDep(
    "pudl-raw-ferc714",
    """
    Dataset: PUDL Raw FERC Form 714
    Website: https://zenodo.org/record/4127101
    Author: Gosnell, Christina, Selvans, Zane
    Date of Publication: October 24, 2020
    License:

    Raw Federal Energy Regulatory Commission (FERC) Form 714 Data, archived from https://www.ferc.gov/industries-data/electric/general-information/electric-industry-forms/form-no-714-annual-electric/data

    Electric transmitting utilities operating balancing authority areas and planning areas with annual peak demand over 200MW are required to file Form 714 with the Federal Energy Regulatory Commission (FERC), reporting balancing authority area generation, actual and scheduled inter-balancing authority area power transfers, and net energy for load, summer-winter generation peaks and system lambda.

    $PUDL_DESCRIPTION
    """,
    [
        "https://zenodo.org/api/files/f083d5ce-b7b2-48a0-8896-925f6f0d0b24/datapackage.json",
        "https://zenodo.org/api/files/f083d5ce-b7b2-48a0-8896-925f6f0d0b24/ferc714.zip",
    ],
    "78f8b660810de6b18dbe91de07ff7a976b0d04964063ec9cb3f3215bfe5269f0",
))


text:
"""
$(TYPEDEF)
Optimizer for finding a sparse basis vector in a subspace based on [this paper](https://arxiv.org/pdf/1412.4659.pdf).
It solves the following problem
```math
\\argmin_{x} \\|x\\|_0 ~s.t.~Ax= 0
```
# Fields
$(FIELDS)

# Example
```julia
ADM()
ADM(λ = 0.1)
```

## Note
While useable for implicit problems, a better choice in general is given by the
`ImplicitOptimizer` which tends to be more robust.
"""
mutable struct ADM{T} <: AbstractSubspaceOptimizer{T}
    """Sparsity threshold"""
    λ::T

    function ADM(threshold::T = 1e-1) where T
        @assert all(one(eltype(T)) .> threshold .> zero(eltype(T))) "Threshold must be positive definite and less than 1"

        return new{T}(threshold)
    end
end

Base.summary(::ADM) = "ADM"

function (opt::ADM{T})(X, A, Y, λ::V = first(opt.λ);
    maxiter::Int64 = maximum(size(A)), abstol::V = eps(eltype(T)), progress = nothing,
    rtol::V = zero(eltype(T)),
    atol::V = convert(eltype(T), 0.99),
    f::Function = F(opt), g::Function = G(opt))  where {T, V}

    n,m = size(A)
    ny, my = size(Y)
    nx, mx = size(X)
    nq, mq = 0,0

    # Closure for the pareto function
    fg(x, A, y) = (g∘f)(x, A, y)
    fg(x, A) = (g∘f)(x,A)

    # Init all variables
    R = SoftThreshold()

    xzero = zero(eltype(T))
    obj = xzero
    sparsity = xzero
    conv_measure = xzero

    iters = 0
    converged = false

    max_ind = 0

    _progress = isa(progress, Progress)
    initial_prog = 0
    N = nullspace(A[:,:], atol = atol)
    Q = deepcopy(N)
    map(x->normalize!(x), eachrow(Q))
    Q_ = deepcopy(N)
    x = N'Q

    @views while (iters < maxiter) && !converged
        iters += 1

        R(x, N'Q, λ)
        mul!(Q, N, x)
        map(x->normalize!(x), eachrow(Q))

        conv_measure = norm(Q .- Q_, 2)

        if _progress
            sparsity, obj = f(Q,A,λ)

            ProgressMeter.next!(
            progress;
            showvalues = [
                (:Threshold, λ), (:Objective, obj), (:Sparsity, sparsity),
                (:Convergence, conv_measure),
                ]
                )
        end

        if conv_measure < abstol
            converged = true
        else
            Q_ .= Q
        end
    end

    clip_by_threshold!(Q, λ)

    # Reduce the solution size to linear independent columns
    @views Q = linear_independent_columns(Q, rtol)

    # Indicate if already used
    _included = zeros(Bool, my, size(Q, 2))

    @views for i in 1:my, j in 1:size(Q, 2)
        # Check, if already included
        any(_included[:, j]) && continue
        if @views evaluate_pareto!(X[:, i], Q[:,j] , fg, A, λ)
            _included[i,j] = true
        end
    end

    if _progress
        sparsity, obj = f(X,A,λ)

        ProgressMeter.update!(
            progress,
            initial_prog + maxiter -1
        )

        ProgressMeter.next!(
        progress;
        showvalues = [
            (:Threshold, λ), (:Objective, obj), (:Sparsity, sparsity),
            ]
        )
    end

    if rank(X'X) < my
        @warn "$opt @ $λ has found illconditioned equations. Vary the threshold or relative tolerance."
    end

    return
end


text:
# This mechanism allows us to import a package list.
# Source: https://discourse.julialang.org/t/building-a-dockerfile-with-packages/37272/2

using Pkg
pkg"add Colors"
pkg"add CSV"
pkg"add CSVFiles"
pkg"add DataFrames"
pkg"add DelimitedFiles"
pkg"add GMT"
pkg"add GR"
pkg"add IJulia"
pkg"add ImageView"
pkg"add ImageInTerminal"
pkg"add FileIO"
pkg"add ImageFiltering"
pkg"add Images"
pkg"add ImageFeatures"
pkg"add Interact"
pkg"add LaTeXStrings"
pkg"add LinearAlgebra"
pkg"add Plots"
pkg"add PyPlot"
pkg"add Random"
pkg"add GtkReactive"
pkg"add Roots"
pkg"add SymEngine"
pkg"add Statistics"
pkg"add TestImages"
pkg"add Plotly"
pkg"add Distributions"
pkg"add ProgressBars"
pkg"add https://github.com/VMLS-book/VMLS.jl"
pkg"add https://github.com/korsbo/Latexify.jl"
pkg"add WebIO"

max_stars_repo_name: jgwerner/umich-stacks

max_issues_repo_name: jgwerner/umich-stacks

max_forks_repo_name: jgwerner/umich-stacks


text:
using PropertyFiles
using Test

@testset "PropertyFiles tests" begin
    @testset "basic tests" begin
        include("basictests.jl")
    end

    @testset "overload tests" begin
        include("overloadtests.jl")
    end

    @testset "load-store tests" begin
        include("storetests.jl")
    end
end


text:
"""
```
const DenseSolutionArray=Matrix
```

Dense storage of solution
"""
const DenseSolutionArray=Matrix

"""
$(SIGNATURES)

Get degree of freedom number
"""
dof(a::DenseSolutionArray{Tv}, ispec::Integer, K::Integer) where Tv = (K-1)*size(a,1)+ispec

unknown_indices(a::DenseSolutionArray{Tv}) where Tv=LinearIndices(a)

"""
$(SIGNATURES)

Array of values in solution array.
"""
values(a::DenseSolutionArray{Tv}) where Tv = vec(a)

#
# Accessors for node-dof based loops
#
_firstnodedof(U::DenseSolutionArray{Tv},K) where Tv = (K-1)*size(U,1)+1
_lastnodedof(U::DenseSolutionArray{Tv},K) where Tv = K*size(U,1)
_spec(U::DenseSolutionArray{Tv},idof,K) where Tv =   idof-(K-1)*size(U,1)
_add(U::DenseSolutionArray{Tv},idof,val) where Tv=U[CartesianIndices(U)[idof]]+=val

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


text:
"""
    check_print(vec::Vector{Symbol}, ps::Vector{Symbol})

Check whether something should get printed. Only if the log_levels option
include the necessary level
"""
function check_print(log_levels::Vector{Symbol}, necessary::Vector{Symbol})
    for v in log_levels
        if v in necessary
            return true
        end
    end
    return false
end

function print_info(m::JuniperProblem)
    println("#Variables: ", m.num_var)
    println("#IntBinVar: ", m.num_disc_var)
    println("Obj Sense: ", m.obj_sense)
    println()
end

function print_dict(d)
    longest_key_name = maximum([length(string(key)) for key in keys(d)])+2
    for key in keys(d)
        skey = string(key)
        pkey = skey*repeat(" ", longest_key_name-length(skey))
        if d[key] == nothing
            println(pkey, ": NA")
        else
            println(pkey, ": ",d[key])
        end
    end
end

function get_non_default_options(options)
    defaults = Juniper.get_default_options()
    non_defaults = Dict{Symbol,Any}()
    for fname in fieldnames(SolverOptions)
        # TODO: Better printing of nl_solver/mip_solver name when SolverName exists

        # doesn't work for arrays but the only array atm is log_levels
        # and the default doesn't include :Options therefore !== should work...
        if getfield(options,fname) !== getfield(defaults,fname)
            non_defaults[fname] = getfield(options,fname)
        end
    end
    return non_defaults
end

function print_options(m::JuniperProblem; all=true)
    if all
        println(m.options)
    else
        print_dict(get_non_default_options(m.options))
    end
    println()
end

function print_fp_table(mip_obj,nlp_obj,t, fields, field_chars, catol)
    ln, arr = get_fp_table(mip_obj,nlp_obj,t, fields, field_chars, catol)
    println(ln)
end

function print_final_timing(time_bnb_solve::Float64, time_obj::TimeObj)
    println("BnB time: ", round(time_bnb_solve; digits=2))
    println("% solve child time: ", round((time_obj.solve_leaves_get_idx+time_obj.solve_leaves_branch)/time_bnb_solve*100; digits=1))
    println("Solve node time get idx: ", round(time_obj.solve_leaves_get_idx; digits=2))
    println("Solve node time branch: ", round(time_obj.solve_leaves_branch; digits=2))
    println("Branch time: ", round(time_obj.branch; digits=2))
    println("Get idx time: ", round(time_obj.get_idx; digits=2))
    println("Upd gains time: ", round(time_obj.upd_gains; digits=2))
end


text:
module Gunrock
    include("config.jl")
	include("bfs.jl")
	include("cc.jl")
	include("sssp.jl")
	include("pagerank.jl")
	include("bc.jl")
	export bfs, cc, sssp, bc, pagerank
end


text:
AccFree_MTiXsens_sync_perturb=Dict()
Xsensor_reduce_sync_perturb=Dict()
IBS_sync_perturb =Dict()

for experiment=1:size(filenameMTiXsens,1)
    #AccFree_MTiXsens_sync_perturb[filenameMTiXsens[experiment,1]]=Dict()
    #for n in sensorXsensNames
    #       AccFree_MTiXsens_sync_perturb[filenameMTiXsens[experiment,1]][n]=AccFree_MTiXsens_sync[filenameMTiXsens[experiment,1]][n][sync_[experiment].MTiXsens:sync_[experiment].MTiXsens+step_max_perturb,:]
    #end

    AccFree_MTiXsens_sync_perturb[filenameMTiXsens[experiment,1]]=Dict()
    for n in sensorXsensNames
        are = AccFree_MTiXsens_sync[filenameMTiXsens[experiment,1]][n][:,"PacketCounter"]
        ligne = findall(x -> x == sync_.MTiXsens[experiment], are)
        AccFree_MTiXsens_sync_perturb[filenameMTiXsens[experiment,1]][n]=
            AccFree_MTiXsens_sync[filenameMTiXsens[experiment,1]][n][
                ligne[1]:min(ligne[1]+step_max_perturb,size(AccFree_MTiXsens_sync[filenameMTiXsens[experiment,1]][n][ligne[1]:end,:],1)),
                :]
    end
end

for experiment=1:size(filenameXsensor,1)
    Xsensor_reduce_sync_perturb[filenameXsensor[experiment,1]]=
        Xsensor_reduce_sync[filenameXsensor[experiment,1]][
            sync_.Xsensor[experiment]:min(sync_.Xsensor[experiment]+convert(Int,step_max_perturb÷(Ts_Xsensor/Ts_MTiXsens)),size(Xsensor_reduce_sync[filenameXsensor[experiment,1]][sync_.Xsensor[experiment]:end,:],1)),
            :]
end

for experiment=1:size(filenameIBS,1)
    IBS_sync_perturb[filenameIBS[experiment,1]] =
        IBS_sync[filenameIBS[experiment,1]][
            sync_.IBS[experiment]:min(sync_.IBS[experiment]+convert(Int,step_max_perturb÷(Ts_IBS/Ts_MTiXsens)),size(IBS_sync[filenameIBS[experiment,1]][sync_.IBS[experiment]:end,:],1)),
            :]
end

println("Synchro perturb data finished")


text:
using Documenter, Parameters
makedocs(
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    sitename = "Parameters.jl",
    # pages also make the side-bar
    pages = Any[
        "index.md",
        "manual.md",
        "api.md"]
)

deploydocs(
    repo = "github.com/mauro3/Parameters.jl.git"
)


text:
#!/usr/bin/env julia
using DocOpt

doc = """
Pope Benhmark and Live Test
Usage:
  benchmark.jl
  benchmark.jl [--nsamples=<ns>] [--nchannels=<nc>] [--cps=<cps>] [--nozmq] [--runtime_s=<rt>]
  benchmark.jl -h | --help

Options:
  --nsamples=<ns>    Samples per pulse. [default: 1000]
  --nchannels=<nc>   Number of channels to read and write. [default: 240]
  --cps=<cps>        Average counts per second per channel. [default: 100]
  --runtime_s=<rt>   Run for roughly this long, in seconds. [default: 30]
  --nozmq            Include this flag to not use the ZMQ port feature.

Help exposition:
  Launch a 2nd process to write LJH files, using an exponential distribution to determine time between
  writing records. Use the 1st process to analyze those LJH files using Pope in realtime.

  If you have succesfully run `Pkg.test("Pope") from within Julia, this should work:
  ./benchmark.jl
"""

arguments = docopt(doc, version=v"0.0.1")
nsamples = parse(Int,arguments["--nsamples"])
nchannels = parse(Int,arguments["--nchannels"])
cps = parse(Float64,arguments["--cps"])
runtime_s = parse(Float64,arguments["--runtime_s"])
nozmq = arguments["--nozmq"]

println("Starting Pope livetest/benchmark.")
if nprocs() == 1
  addprocs(1)
  println("added 1 proc")
end

@everywhere begin
using ProgressMeter
using Pope: LJH, LJH
using HDF5
using Distributions: Distribution, Exponential
"    gen_ljh_files(dt=9.6e-6, npre=200, nsamp=1000,channels=1:2:480,dname=tempdir(); version=\"2.2.0\")
generate ljhs files (1 per channel) in directory `dname`, with the specified parameters
return a Vector{LJHFile} of the created files, with both read and write intents"
function gen_ljh_files(dt=9.6e-6, npre=200, nsamp=1000,channels=1:2:480,dname=joinpath(tempdir(),randstring(12)); version="2.2.0")
  !isdir(dname) && mkdir(dname)
  basename = last(split(dname,'/'))
  ljhs = LJH.LJHFile[]
  for ch in channels
    fname = joinpath(dname, basename*"_chan$ch.ljh")
    ljh=LJH.create(fname,dt, npre, nsamp; version=version, channel=ch)
    push!(ljhs,ljh)
  end
  ljhs
end
function launch_stocastic_writer{T}(endchannel, ljh::LJH.LJHFile{LJH.LJH_22,T}, d::Distribution, data=zeros(UInt16,ljh.record_nsamples))
  @schedule begin
    i=0
    while !isready(endchannel)
      sleep(rand(d))
      write(ljh,data,i+=1,round(Int, 1e6*time()))
    end
    close(ljh)
  end
end

function launch_stocastic_writers(endchannel, ljhs::Vector{LJH.LJHFile}, d::Distribution, data)
  for ljh in ljhs
    launch_stocastic_writer(endchannel, ljh, d, data)
  end
end

ulimit() = a=parse(Int,readstring(`bash -c "ulimit -n"`))

"    launch_writer_other_process(;d=Exponential(0.01),data=zeros(UInt16,1000),dt=9.6e-6, npre=200, nsamp=length(data),channels=1:2:480,dname=joinpath(tempdir(),randstring(12)), version=\"2.2.0\")
Opens one LJH file per channel in `channels` and starts writing pulse records with `data`
"
function launch_writer_other_process(;d=Exponential(0.01),data=zeros(UInt16,1000),dt=9.6e-6, npre=200, nsamp=length(data),channels=1:2:480,dname=joinpath(tempdir(),randstring(12)), version="2.2.0")
  nprocs() >= 2 || error("nprocs needs to be 2 or greater, try starting julia with `julia -p 1`")
  ulimit() >= 50+length(channels) || error("open file limit too low, run `ulimit -n 1000` before opening julia")
  endchannel = RemoteChannel(1)
  ljhmadechannel = RemoteChannel(1)
  s=@spawn begin
    ljhs = gen_ljh_files(dt, npre, nsamp, channels, dname; version=version)
    localendchannel = Channel(1)
    launch_stocastic_writers(localendchannel, ljhs, d, data)
    put!(ljhmadechannel,true)
    @schedule begin
      # checking a RemoteChannel causes lots of CPU usage in both tasks, avoid checking it alot
      wait(endchannel)
      put!(localendchannel,true)
    end
    @show s
  end
  return dname, endchannel, ljhmadechannel, s
end
end #@everywhere
data = collect(UInt16(0):UInt16(nsamples-1))
filter_values = zeros(nsamples-1)
filter_at = zeros(nsamples-1)
npresamples = 200
average_pulse_peak_index = 250
shift_threshold = 5
pretrigger_rms_cuts = postpeak_deriv_cuts = [0.0,0.0]
analyzer_pk_string = "manually created in benchmark.jl"
channels = 1:2:2*nchannels
frametime = 9.6e-6
λ = 1/cps # mean time (seconds) between pulses
d = Exponential(λ)
dname = joinpath(tempdir(),randstring(12))
outputname = tempname()

println("reading/analyzing process pid $(getpid())")
wait(@spawn println("writing process pid (probably) $(getpid())"))
println("writing and analzing $(length(channels)) channels")
println("$(1/λ) counts per second per channel")
println("$(length(data)) samples per pulse, 2 bytes/sample")
println("output file (wrriten by 1st process) $outputname")
println("ljh files (written by 2nd process) in $dname")
println("check activity monitor or top for CPU/memory usage")
println("temporary files are probably delted after you reboot, but you may want to check on your platform")

println("Starting writing")
dname,endchannel,ljhmadechannel,s=launch_writer_other_process(d=d,dt=frametime,channels=channels,npre=npresamples,data=data, dname=dname)
wait(ljhmadechannel) # makes sure all the LJH files are created before moving on
println("Writing in progress")

nozmq || Pope.init_for_zmqdatasink(Pope.ZMQ_PORT,verbose=true)
println("Starting analyzing")
h5 = h5open(outputname,"w", "libver_bounds", (HDF5.H5F_LIBVER_LATEST, HDF5.H5F_LIBVER_LATEST))
readers = Pope.Readers()
fname=""
for channel in channels
  fname = LJH.fnames(dname, channel)
  analyzer = Pope.MassCompatibleAnalysisFeb2017(filter_values, filter_at, npresamples, nsamples, average_pulse_peak_index, frametime, shift_threshold, pretrigger_rms_cuts, postpeak_deriv_cuts, analyzer_pk_string)
  if nozmq
    product_writer = Pope.make_buffered_hdf5_writer(h5, channel, analyzer)
  else
    product_writer = Pope.make_buffered_hdf5_and_zmq_multisink(h5, channel, analyzer)
  end
  reader = Pope.make_reader(fname, analyzer, product_writer)
  push!(readers, reader)
end
schedule(readers)
println("Analyzing in progress")
println("Analyzing for $runtime_s seconds")
wait(@schedule begin
  tstart = time()
  tdiff = time()-tstart
  runtime_ms = round(Int, 1000*runtime_s)
  p = Progress(runtime_ms,1,"livetest/benchmark: ")
  while tdiff<runtime_s
    ProgressMeter.update!(p,round(Int,1000*(tdiff)))
    sleep(min(1,tdiff))
    tdiff = time()-tstart
  end
  ProgressMeter.update!(p,runtime_ms)
  put!(endchannel,true)
  println("writing stopped") end);
sleep(3) # make sure ljh files are all fully written, I get errors without this
Pope.stop(readers)
wait(readers)
println("analyzing stopped")

function check_values(channel, h5)
  fname = LJH.fnames(dname, channel)
  ljh = LJH.LJHFile(fname)
  # @assert 80<length(ljh)/runtime_s<120  # the efault value in launch_writer_other_process has 100 cps
  g = h5["chan$channel"]
  for name in names(g)
    if name == "calculated_cuts"
      continue
    end
    length(g[name])==length(ljh) || @show name, length(g[name]), length(ljh), fname
    @assert length(g[name])==length(ljh)
  end
  @assert(read(g["rowcount"])==collect(1:length(ljh)))
end

println("sanity checking analysis output")
for channel in channels
  check_values(channel,h5)
end


text:
function lineposterior(; time=time, flux=flux, σ2=σ2, idx=idx, maxiter = 10, maxinit = 10, maxrandom = 10, seed = 1, show_trace = true, usercheck = true)

    logs, predictTest = multigp(time, flux, σ2, idx; maxiter = maxinit,
                                                   maxrandom = maxrandom,
                                                        seed = seed,
                                                  show_trace = show_trace)

    if usercheck
        for i in unique(idx)

            figure()
            cla()
            title(@sprintf("filter index %d", i))
            xtest = collect(LinRange(minimum(time), maximum(time), 200))
            μpred, Σpred = predictTest(xtest, i)
            PyPlot.plot(time[findall(idx.==i)], flux[findall(idx.==i)], "ko")
            PyPlot.plot(xtest, μpred, "r-")
            PyPlot.fill_between(xtest, μpred - sqrt.(diag(Σpred)), μpred + sqrt.(diag(Σpred)), color="r", alpha=0.1)

        end
    end

    @show logs

    # gpmodelVI(time, flux, σ2, idx; logs=logs, maxiterml=15_000, maxitervi=maxiter, S=200, seed=seed)

    gpmodelVI(time, flux, σ2, idx; logs=logs, maxiter = maxiter, seed = seed, show_trace = show_trace)

end


text:
"""
Default Speed Limits in Kilometers Per Hour in urban areas
"""
const SPEED_ROADS_URBAN = Dict{Int,Float64}(
    1 => 100,    # Motorway
    2 => 90,    # Trunk
    3 => 90,    # Primary
    4 => 70,    # Secondary
    5 => 50,    # Tertiary
    6 => 40,    # Residential/Unclassified
    7 => 20,     # Service
    8 => 10)     # Living street
"""
Default Speed Limits in Kilometers Per Hour in rural areas
"""
const SPEED_ROADS_RURAL = Dict{Int,Float64}(
    1 => 100,
    2 => 90,
    3 => 90,
    4 => 70,
    5 => 50,
    6 => 40,
    7 => 40,
    8 => 30)


text:
using LinearAlgebra,Distributions,DataFrames,StatsBase,GLM,LsqFit,Optim,BlackBoxOptim,HypothesisTests,Colors,Images,StatsModels,CategoricalArrays,
ImageFiltering,SpecialFunctions,DSP,HCubature,Combinatorics,DataStructures,ANOVA,StatsFuns,Trapz,CircStats,ImageSegmentation,ProgressMeter,PyCall,
Dierckx,BandedMatrices,OffsetArrays
import Base: vec,range
import StatsBase: predict

include("NeuroDataType.jl")
include("Stats.jl")
include("Spike.jl")
include("LFP.jl")
include("Image.jl")
include("Condition.jl")
include("2P.jl")

vec(x::RGBA)=[x.r,x.g,x.b,x.alpha]
anscombe(x) = 2*sqrt(x+(3/8))

"`Gaussian` function"
gaussianf(x;a=1,μ=0,σ=1) = a*exp(-0.5((x-μ)/σ)^2)
function gaussianf(x,y;a=1,μ₁=0,σ₁=1,μ₂=0,σ₂=1,θ=0)
    sinv,cosv = sincos(θ)
    x₀ = x-μ₁
    y₀ = y-μ₂
    x′ = cosv * x₀ + sinv * y₀
    y′ = cosv * y₀ - sinv * x₀
    a*exp(-0.5((x′/σ₁)^2 + (y′/σ₂)^2))
end

"""
`von Mises` function [^1]

```math
f(α) =  βe^{κ(cos(n(α - μ)) - 1)}
```

[^1]

Swindale, N.V. (1998). Orientation tuning curves: empirical description and estimation of parameters. Biol Cybern 78, 45–56.

- β: amplitude at μ
- μ: angle of peak
- κ: width parameter
- n: frequency parameter
"""
vmf(α;β=1,μ=0,κ=1,n=1) = β*exp(κ*(cos(n*(α-μ))-1))
"""
`Generalized von Mises` function [^1]

```math
f(α) =  βe^{κ₁(cos(α - μ₁) - 1) + κ₂(cos2(α - μ₂) - 1)}
```

[^1]

Gatto, R., and Jammalamadaka, S.R. (2007). The generalized von Mises distribution. Statistical Methodology 4, 341–353.
"""
gvmf(α;β=1,μ₁=0,κ₁=1,μ₂=0,κ₂=1) = β*exp(κ₁*cos(α-μ₁) + κ₂*cos(2(α-μ₂)))

"""
`Difference of Gaussians` function
"""
dogf(x;aₑ=2,μₑ=0,σₑ=1,aᵢ=1,μᵢ=0,σᵢ=2) = aₑ*exp(-0.5((x-μₑ)/σₑ)^2) - aᵢ*exp(-0.5((x-μᵢ)/σᵢ)^2)
function dogf(x,y;aₑ=2,μₑ₁=0,σₑ₁=1,μₑ₂=0,σₑ₂=1,θₑ=0,aᵢ=1,μᵢ₁=0,σᵢ₁=2,μᵢ₂=0,σᵢ₂=2,θᵢ=0)
    sinvₑ,cosvₑ = sincos(θₑ)
    xₑ₀ = x-μₑ₁
    yₑ₀ = y-μₑ₂
    xₑ′ = cosvₑ * xₑ₀ + sinvₑ * yₑ₀
    yₑ′ = cosvₑ * yₑ₀ - sinvₑ * xₑ₀
    sinvᵢ,cosvᵢ = sincos(θᵢ)
    xᵢ₀ = x-μᵢ₁
    yᵢ₀ = y-μᵢ₂
    xᵢ′ = cosvᵢ * xᵢ₀ + sinvᵢ * yᵢ₀
    yᵢ′ = cosvᵢ * yᵢ₀ - sinvᵢ * xᵢ₀
    aₑ*exp(-0.5((xₑ′/σₑ₁)^2 + (yₑ′/σₑ₂)^2)) - aᵢ*exp(-0.5((xᵢ′/σᵢ₁)^2 + (yᵢ′/σᵢ₂)^2))
end

"""
`sin` grating function

- μ: x offset
- f: Frequency in cycle/unit_x
- phase: Phase of a cycle in [0, 1] scale
"""
gratingf(x;μ=0, f=1, phase=0) = sin(2π * (f * (x-μ) + phase))

"""
2D `sin` grating function

- μ₁: x offset
- μ₂: y offset
- θ: Orientation in radius, 0 is -, increase counter-clock wise
- f: Frequency in cycle/unit_x/y
- phase: Phase of a cycle in [0, 1] scale
"""
function gratingf(x,y;μ₁=0,μ₂=0,θ=0,f=1,phase=0)
    sinθ,cosθ = sincos(θ)
    y′ = cosθ * (y-μ₂) - sinθ * (x-μ₁)
    sin(2π * (f * y′ + phase))
end

"""
`cas` function defined as ``cas(x) = cos(x) + sin(x)``

- f: Frequency in cycle/unit_x
- phase: Phase of a cycle in [0, 1] scale
- isnorm: scale `cas` in [-√2, √2] to [-1, 1]
"""
function cas(x;f=1, phase=0, isnorm::Bool=true)
    r = sum(sincos(2π * (f * x + phase)))
    if isnorm
        r /=sqrt(2)
    end
    return r
end

"""
2D `cas` function defined as ``cas(x+y) = cos(x+y) + sin(x+y)``

- kx: Frequency in cycle/unit_x
- ky: Frequency in cycle/unit_y
- phase: Phase of a cycle in [0, 1] scale
- isnorm: scale `cas` in [-√2, √2] to [-1, 1]
"""
function cas(x,y;kx=1,ky=1, phase=0, isnorm::Bool=true)
    r = sum(sincos(2π * (kx * x + ky * y + phase)))
    isnorm && (r /=sqrt(2))
    return r
end

"""
`cas` phase to `sin` phase, phase is in [0, 1] scale

!!! note
    ``cas(x) = √2 sin(x + π/4)``
"""
cas2sin(phase) = phase + 0.125

"""
2D `cas` to 2D sin `gratingf`
"""
function cas2sin(kx,ky,phase)
    θ = atan(ky,kx) - π/2
    f = sqrt(kx*kx + ky*ky)
    return (θ=θ,f=f,phase=phase + 0.125)
end

"""
`sin` phase to `cas` phase, phase is in [0, 1] scale

!!! note
    ``cas(x) = √2 sin(x + π/4)``
"""
sin2cas(phase) = phase - 0.125

"""
2D sin `gratingf` to 2D `cas`
"""
function sin2cas(θ,f,phase)
    sinθ,cosθ = sincos(θ + π/2)
    kx = cosθ*f
    ky = sinθ*f
    return (kx=kx,ky=ky,phase=phase - 0.125)
end

"`Gabor` function"
gaborf(x;a=1,μ=0,σ=1,f=1,phase=0) = a*exp(-0.5((x-μ)/σ)^2)*sin(2π*(f*(x-μ)+phase))
function gaborf(x,y;a=1,μ₁=0,σ₁=1,μ₂=0,σ₂=1,θ=0,f=1,phase=0)
    sinv,cosv = sincos(θ)
    x₀ = x-μ₁
    y₀ = y-μ₂
    x′ = cosv * x₀ + sinv * y₀
    y′ = cosv * y₀ - sinv * x₀
    a*exp(-0.5((x′/σ₁)^2 + (y′/σ₂)^2)) * sin(2π*(f * y′ + phase))
end

"Binary mask for gabor envelope"
function gaborenvelopemask(x,y;fσ=2.5,μ₁=0,σ₁=1,μ₂=0,σ₂=1,θ=0)
    sinv,cosv = sincos(θ)
    x₀ = x-μ₁
    y₀ = y-μ₂
    x′ = cosv * x₀ + sinv * y₀
    y′ = cosv * y₀ - sinv * x₀
    (x′/fσ/σ₁)^2 + (y′/fσ/σ₂)^2 <= 1 ? true : false
end
"Binary mask for concentric circular dog envelope"
dogenvelopemask(x,y;fσ=2.5,μ₁=0,σₑ₁=1,μ₂=0,rσᵢₑ=2) = edogenvelopemask(x,y;fσ,μ₁,σₑ₁,rσ₂₁=1,μ₂,rσᵢₑ,θ=0)
"Binary mask for concentric orientated elliptical dog envelope"
function edogenvelopemask(x,y;fσ=2.5,μ₁=0,σₑ₁=1,rσ₂₁=1,μ₂=0,rσᵢₑ=2,θ=0)
    sinv,cosv = sincos(θ)
    x₀ = x-μ₁
    y₀ = y-μ₂
    x′ = cosv * x₀ + sinv * y₀
    y′ = cosv * y₀ - sinv * x₀
    σ₁ = max(σₑ₁,rσᵢₑ*σₑ₁)
    σ₂ = σ₁ * rσ₂₁
    (x′/fσ/σ₁)^2 + (y′/fσ/σ₂)^2 <= 1 ? true : false
end

"Fit 1D model to data"
function fitmodel(model,x,y)
    lb,ub = extrema(y)
    bm = (lb+ub)/2
    br = (ub-lb)/2
    alb,aub = abs.((lb,ub))
    ab = max(alb,aub)

    xlb,xub = extrema(x)
    xbm = (xlb+xub)/2
    xalb,xaub = abs.((xlb,xub))
    xab = max(xalb,xaub)

    rlt = fun = missing
    if model == :vmn2
        fun = (x,p) -> vmf.(x,β=p[1],μ=p[2],κ=p[3],n=2)
        ofun = (p;x=x,y=y) -> sum((y.-fun(x,p)).^2)

        ub=[1.8ab,   prevfloat(float(π)),   200]
        lb=[0.2ab,            0,              0]
        p0=[ab,               0,              1]
    elseif model == :gvm
        # fun = (x,p) -> gvmf.(x,β=p[1],μ₁=p[2],κ₁=p[3],μ₂=p[4],κ₂=p[5])
        fun = gvmff
        ofun = (p;x=x,y=y) -> sum((y.-fun(x,p)).^2)

        ub=[1.8ab,   prevfloat(float(2π)),   40,    prevfloat(float(π)),     40]
        lb=[nextfloat(0.0),   0,              0,             0,               0]
        p0=[ab,               0,              1,             0,               1]
    elseif model == :dog
        fun = (x,p) -> dogf.(x,aₑ=p[1],μₑ=p[2],σₑ=p[3],aᵢ=p[4],μᵢ=p[5],σᵢ=p[6])
        ofun = (p;x=x,y=y) -> sum((y.-fun(x,p)).^2)

        ub=[1.8ab,   10xab,            10xab,    1.8ab,     10xab,              10xab]
        lb=[0,      -10xab,   nextfloat(0.0),        0,    -10xab,     nextfloat(0.0)]
        p0=[ab,        0,                  1,       ab,         0,                  1]
    elseif model == :sfdog
        # fun = (x,p) -> dogf.(x,aₑ=p[1],μₑ=p[2],σₑ=p[3],aᵢ=p[4],μᵢ=p[5],σᵢ=p[6]) .+ p[7]
        fun = sfdogff
        ofun = (p;x=x,y=y) -> sum((y.-fun(x,p)).^2)

        ub=[1.5ab,    10,             10,          1.5ab,       10,                10,         bm+br/3]
        lb=[0,         0,      nextfloat(0.0),       0,          0,        nextfloat(0.0),        0]
        p0=[ab,     0.5xab,         0.5xab,          0,        0.5xab,           0.5xab,       bm-br/3]
    elseif model == :sfgaussian
        # fun = (x,p) -> gaussianf.(log2.(x),a=p[1],μ=p[2],σ=p[3]) .+ p[4]
        fun = sfgaussianff
        ofun = (p;x=x,y=y) -> sum((y.-fun(x,p)).^2)

        # limit gaussian center ~[-8 10], sigma ~[1.4 8]
        ub=[1.5ab,     3.3,          3,       bm]
        lb=[0,          -3,        0.5,        0]
        p0=[ab,          0,          1,        0]
    end
    if !ismissing(fun)
        # ofit = optimize(ofun,lb,ub,p0,SAMIN(rt=0.92),Optim.Options(iterations=220000))
        # param=ofit.minimizer
        ofit = bboptimize(ofun,p0;SearchRange=collect(zip(lb,ub)),Method=:adaptive_de_rand_1_bin_radiuslimited,MaxSteps=200000)
        param = best_candidate(ofit)

        rlt = (;model,fun,param, goodnessoffit(y,fun(x,param),k=length(param))...)
    end
    return rlt
end
gvmff(x,p) = gvmf.(x,β=p[1],μ₁=p[2],κ₁=p[3],μ₂=p[4],κ₂=p[5])
sfdogff(x,p) = dogf.(x,aₑ=p[1],μₑ=p[2],σₑ=p[3],aᵢ=p[4],μᵢ=p[5],σᵢ=p[6]) .+ p[7]
sfgaussianff(x,p) = gaussianf.(log2.(x),a=p[1],μ=p[2],σ=p[3]) .+ p[4]

"Fit 2D model to image"
function fitmodel2(model,data::Matrix,ppu;w=0.5)
    rspx = (size(data).-1)./2
    radii = rspx./ppu
    x = (mapreduce(i->[i[2] -i[1]],vcat,CartesianIndices(data)) .+ [-(rspx[2]+1) (rspx[1]+1)])/ppu
    y = vec(data)

    # try estimate solution
    roi = peakroi(localcontrast(data,round(Int,w*ppu)))
    alb,aub = abs.(extrema(data[roi.i]))
    ab = max(alb,aub)
    r = roi.radius/ppu
    c = [roi.center[2] - (rspx[2]+1), -roi.center[1] + (rspx[1]+1)]/ppu

    rlt = fun = missing
    if model == :dog
        if aub >= alb
            ai = 5alb
            ae = aub + ai
        else
            ae = 5aub
            ai = alb + ae
        end
        # concentric circular dog
        # fun = (x,y,p) -> dogf.(x,y,aₑ=p[1],μₑ₁=p[2],σₑ₁=p[3],μₑ₂=p[4],σₑ₂=p[3],θₑ=0,aᵢ=p[5],μᵢ₁=p[2],σᵢ₁=p[6]*p[3],μᵢ₂=p[4],σᵢ₂=p[6]*p[3],θᵢ=0)
        # mfun = (x,y,p) -> dogenvelopemask.(x,y;fσ=2.5,μ₁=p[2],μ₂=p[4],σₑ₁=p[3],rσᵢₑ=p[6])
        fun = dogff
        mfun = dogfmf
        ofun = (p;x=x,y=y) -> @views sum((y.-fun(x[:,1],x[:,2],p)).^2)
        ub=[5ae,    0.5r+c[1],    0.9r,    0.5r+c[2],     5ai,    4]
        lb=[0,     -0.5r+c[1],    0.1r,   -0.5r+c[2],     0,      0.25]
        p0=[ae,      c[1],        0.3r,     c[2],         ai,     1]
    elseif model == :edog
        if aub >= alb
            ai = 5alb
            ae = aub + ai
        else
            ae = 5aub
            ai = alb + ae
        end
        # concentric orientated elliptical dog
        # fun = (x,y,p) -> dogf.(x,y,aₑ=p[1],μₑ₁=p[2],σₑ₁=p[3],μₑ₂=p[4],σₑ₂=p[5]*p[3],θₑ=p[6],aᵢ=p[7],μᵢ₁=p[2],σᵢ₁=p[8]*p[3],μᵢ₂=p[4],σᵢ₂=p[5]*p[8]*p[3],θᵢ=p[6])
        # mfun = (x,y,p) -> edogenvelopemask.(x,y;fσ=2.5,μ₁=p[2],μ₂=p[4],σₑ₁=p[3],rσ₂₁=p[5],rσᵢₑ=p[8],θ=p[6])
        fun = edogff
        mfun = edogfmf
        ofun = (p;x=x,y=y) -> @views sum((y.-fun(x[:,1],x[:,2],p)).^2)
        ub=[5ae,    0.5r+c[1],    0.9r,    0.5r+c[2],       1,       prevfloat(float(π)),    5ai,     4]
        lb=[0,     -0.5r+c[1],    0.1r,   -0.5r+c[2],      0.5,              0,               0,      0.25]
        p0=[ae,      c[1],        0.3r,     c[2],           1,               0,               ai,     1]
    elseif model == :gabor
        # fun = (x,y,p) -> gaborf.(x,y,a=p[1],μ₁=p[2],σ₁=p[3],μ₂=p[4],σ₂=p[5]*p[3],θ=p[6],f=p[7],phase=p[8])
        # mfun = (x,y,p) -> gaborenvelopemask.(x,y;fσ=2.5,μ₁=p[2],μ₂=p[4],σ₁=p[3],σ₂=p[5]*p[3],θ=p[6])
        fun = gaborff
        mfun = gaborfmf
        ofun = (p;x=x,y=y) -> @views sum((y.-fun(x[:,1],x[:,2],p)).^2)

        ori,sf = f1orisf(powerspectrum2(data,ppu)...)
        ub=[5ab,   0.5r+c[1],   0.9r,    0.5r+c[2],    6.0,     prevfloat(float(π)),     12,            prevfloat(1.0)]
        lb=[0,    -0.5r+c[1],   0.1r,   -0.5r+c[2],    0.2,             0,              0.05,                  0]
        p0=[ab,    c[1],        0.3r,     c[2],          1,            ori,     clamp(sf,lb[7],ub[7]),       0.5]
    end
    if !ismissing(fun)
        # ofit = optimize(ofun,lb,ub,p0,SAMIN(rt=0.92),Optim.Options(iterations=220000))
        # param = ofit.minimizer
        ofit = bboptimize(ofun,p0;SearchRange=collect(zip(lb,ub)),Method=:adaptive_de_rand_1_bin_radiuslimited,MaxSteps=200000)
        param = best_candidate(ofit)

        @views rlt = (;model,fun,mfun,param,radii, goodnessoffit(y,fun(x[:,1],x[:,2],param),k=length(param))...)
    end
    return rlt
end
# workaround since JLD2 can not save/load anonymous functions yet
dogff(x,y,p) = dogf.(x,y,aₑ=p[1],μₑ₁=p[2],σₑ₁=p[3],μₑ₂=p[4],σₑ₂=p[3],θₑ=0,aᵢ=p[5],μᵢ₁=p[2],σᵢ₁=p[6]*p[3],μᵢ₂=p[4],σᵢ₂=p[6]*p[3],θᵢ=0)
edogff(x,y,p) = dogf.(x,y,aₑ=p[1],μₑ₁=p[2],σₑ₁=p[3],μₑ₂=p[4],σₑ₂=p[5]*p[3],θₑ=p[6],aᵢ=p[7],μᵢ₁=p[2],σᵢ₁=p[8]*p[3],μᵢ₂=p[4],σᵢ₂=p[5]*p[8]*p[3],θᵢ=p[6])
gaborff(x,y,p) = gaborf.(x,y,a=p[1],μ₁=p[2],σ₁=p[3],μ₂=p[4],σ₂=p[5]*p[3],θ=p[6],f=p[7],phase=p[8])
dogfmf(x,y,p) = dogenvelopemask.(x,y;fσ=2.5,μ₁=p[2],μ₂=p[4],σₑ₁=p[3],rσᵢₑ=p[6])
edogfmf(x,y,p) = edogenvelopemask.(x,y;fσ=2.5,μ₁=p[2],μ₂=p[4],σₑ₁=p[3],rσ₂₁=p[5],rσᵢₑ=p[8],θ=p[6])
gaborfmf(x,y,p) = gaborenvelopemask.(x,y;fσ=2.5,μ₁=p[2],μ₂=p[4],σ₁=p[3],σ₂=p[5]*p[3],θ=p[6])

predict(fit,x) = fit.fun(x,fit.param)
function predict(fit,x,y;xygrid=true,yflip=false)
    if xygrid
        z = [fit.fun(i,j,fit.param) for j in y, i in x]
        yflip && (z=reverse(z,dims=1))
    else
        z = fit.fun(x,y,fit.param)
    end
    z
end

"""
Goodness of Fit Metrics:

- r: Pearson Correlation Coefficient
- mae: Mean Absolute Error
- rmse: Root Mean Squared Error
- rae: Relative Absolute Error
- rse: Relative Squared Error
- r2: R Squared
- adjr2: Adjusted-R²
- s: Residual Standard Error
- aic: Akaike Information Criterion
- bic: Bayesian Information Criterion

1. y: responses
2. ŷ: model predictions

- n: sample size
- e: errors(y - ŷ)
- k: number of predictors
- df: degree of freedom(n - k - 1)
"""
function goodnessoffit(y,ŷ;n = length(y),e = y .- ŷ,k=missing,df = n-k-1)
    r = cor(y,ŷ)
    ae = abs.(e)
    e2 = e.^2
    ssᵣ = sum(e2)
    mae = mean(ae)
    rmse = sqrt(mean(e2))
    ydm = y .- mean(y)
    ssₜ = sum(ydm.^2)
    rae = sum(ae)/sum(abs.(ydm))
    fvu = ssᵣ/ssₜ
    rse = sqrt(fvu)
    r2 = 1 - fvu
    s = ssᵣ/df
    s = s < 0 ? missing : sqrt(s)
    adjr2 = 1 - fvu*(n-1)/df
    aic = n*log(ssᵣ) + 2k
    bic = n*log(ssᵣ/n) + k*log(n) # under iid normal error assumption
    (;r,mae,rmse,rae,rse,r2,adjr2,s,aic,bic)
end

function searchclosest(v,vs;start::Integer=1,step::Integer=1,circ=false)
    n=length(vs);ssign = sign(vs[start]-v)
    i = start
    for _ in 1:n
        if !circ
            i<1 && return -Inf
            n<i && return Inf
        end
        sign(vs[i]-v) != ssign && return i
        i += step
        if circ
            i<1 && (i+=n)
            n<i && (i-=n)
        end
    end
    Inf
end

"""
left and right half width of `v` relative to `y[start]`, -Inf/Inf when no `v` is found.

- start: index of `y` which is the center of the width
- v: value on which width is cutoff
- circ: whether `y` is defined on circular domain and width can wrap around
- x: domain of `y`, return width when `x` provided, otherwise return cutoff indices
"""
function halfwidth(y;start=argmax(y),v=y[start]/2,circ=false,x=nothing)
    li = searchclosest(v,y;start,step=-1,circ)
    ri = searchclosest(v,y;start,step=1,circ)
    if isnothing(x)
        return li,ri
    else
        if circ
            lw = isinf(li) ? li : li<=start ? abs(x[start]-x[li]) : abs(x[start]-x[1])+abs(x[end]-x[li])
            rw = isinf(ri) ? ri : ri>=start ? abs(x[ri]-x[start]) : abs(x[ri]-x[1])+abs(x[end]-x[start])
        else
            lw = isinf(li) ? li : abs(x[start]-x[li])
            rw = isinf(ri) ? ri : abs(x[ri]-x[start])
        end
        return lw,rw
    end
end

circtuningfeature(mfit;od=[π,0.5π],fn=:a,x = 0:0.002:2π) = circtuningfeature(x,predict(mfit,x);od,fn) # 0.002rad ≈ 0.11deg
"""
Properties of Circular Tuning

    - Prefered Angle with Peak Response
    - Selectivity Index
        - version 1: (ResponsePeak - ResponseOpposing)/ResponsePeak
        - version 2: (ResponsePeak - ResponseOpposing)/(ResponsePeak + ResponseOpposing)
    - Half Width at Half Peak-to-Trough

1. x: angles in radius
2. y: responses
- od: opposing angle distance to prefered angle, e.g. π for DSI, 0.5π for OSI
- fn: factor name
"""
function circtuningfeature(x,y;od=[π,0.5π],fn=:a)
    maxr,maxi = findmax(y)
    minr,mini = findmin(y)
    maxx = x[maxi]
    ox = maxx.+od
    _,oi = findclosestangle(x,ox)
    or = y[oi]

    si1 = 1 .- or./maxr
    si2 = (maxr .- or)./(maxr .+ or)
    hw = halfwidth(y,start=maxi,v=(maxr+minr)/2,circ=true,x=x)

    (;Symbol(:p,fn)=>rad2deg(mod2pi(maxx)),Symbol(fn,:hw)=>rad2deg.(hw),Symbol(fn,:si1)=>si1,Symbol(fn,:si2)=>si2,Symbol(fn,:od)=>od)
end

sftuningfeature(mfit;x = 0:0.003:10) = sftuningfeature(x,predict(mfit,x))
"""
Properties of Spatial Frequency Tuning

    - Prefered Spatial Frequency with Peak Response
    - Half Width at Half Peak-to-Trough
    - Freq Passing Type {A:All Pass, H:High Pass, L:Low Pass, B:Band Pass}
    - Bandwidth ``log2(H_cut/L_cut)``
    - Passwidth at Half Peak-to-Trough constrained by low freq lim and high freq lim

1. x: sf in cycle/degree
2. y: responses
"""
function sftuningfeature(x,y;low=minimum(x),high=maximum(x))
    maxr,maxi = findmax(y)
    minr,mini = findmin(y)
    maxx = x[maxi]

    hw = halfwidth(y,start=maxi,v=(maxr+minr)/2,circ=false,x=x)
    pt = all(isinf.(hw)) ? 'A' : isinf(hw[1]) ? 'L' : isinf(hw[2]) ? 'H' : 'B'
    bw = log2((maxx+hw[2])/(maxx-hw[1]))
    pw = pt == 'A' ? high-low : pt == 'L' ? maxx-low+hw[2] : pt == 'H' ? high-maxx+hw[1] : sum(hw)

    (;psf=maxx,sfhw=hw,sftype=pt,sfbw=bw,sfpw=pw)
end

"""
Tuning properties of factor response

1. fl: factor levels
2. fr: factor responses for each level

    Angle, Orientation and Direction follow the same convention such that 0 is -/→, then increase counter-clock wise.

    For cases where Orientation and Direction are interlocked(drifting grating):
        - when Orientation is -(0), then Direction is ↑(90)
        - when Direction is →(0), then Orientation is |(-90)
"""
function factorresponsefeature(fl,fr;fm=mean.(fr),factor=:Ori,isfit::Bool=true)
    i = .!ismissing.(fr)
    fl = fl[i];fr=fr[i];fm=fm[i]
    ls = mapreduce((l,r)->fill(l,length(r)),append!,fl,fr)
    rs = mapreduce(deepcopy,append!,fr)

    if factor in [:Ori,:Ori_Final]
        θ = mod2pi.(deg2rad.(ls))
        α = mod2pi.(deg2rad.(fl))
        d = mean(diff(sort(unique(α)))) # angle spacing
        up, = circ_otest(α,w=fm) # Omnibus test for non-uniformity
        # for orientation
        oα = mod.(α,π)
        ocv, = circ_var(2oα,w=fm,d=2d)
        ocm, = circ_mean(2oα,w=fm)
        ocm = rad2deg(mod2pi(ocm)/2)
        # for direction
        dcv, = circ_var(α;w=fm,d)
        dcm, = circ_mean(α,w=fm)
        dcm = rad2deg(mod2pi(dcm+0.5π))
        maxr,maxi = findmax(fm)
        maxl = fl[maxi]

        fit = ()
        if isfit
            try
                mfit = fitmodel(:gvm,θ,rs) # fit Generalized von Mises
                fit = (;circtuningfeature(mfit,od=[π,0.5π],fn=:o)...,mfit)
            catch
                display.(stacktrace(catch_backtrace()))
            end
        end

        return (;up,ocv,ocm,dcv,dcm,max=maxl=>maxr,fit)
    elseif factor == :Dir
        θ = deg2rad.(fl)
        d = mean(diff(sort(unique(θ)))) # angle spacing
        # for orientation
        oθ = mod.(θ.-0.5π,π)
        om, = circ_mean(2oθ,w=fr)
        oo = rad2deg(mod(angle(om),2π)/2)
        ocv = circ_var(2oθ,w=fr,d=2d)
        # for direction
        dm, = circ_mean(θ,w=fr)
        od = rad2deg(mod(angle(dm),2π))
        dcv = circ_var(θ;w=fr,d)
        # fit Generalized von Mises for direction
        fit = ()
        if isfit
            try
                mfit = fitmodel(:gvm,θ,fr)
                fit = (circtuningfeature(mfit,od=π,fn=:d)...,gvm=mfit)
            catch
            end
            # fit von Mises for orientation
            try
                mfit = fitmodel(:vmn2, θ.-0.5π,fr)
                fit =(fit...,circtuningfeature(mfit, od=0.5π,fn=:o)...,vmn2=mfit)
            catch
            end
        end

        return (;dm,od,dcv,om,oo,ocv,fit)
    elseif factor == :SpatialFreq
        up = PyOnewayANOVA.anova_oneway(fr,use_var="unequal").pvalue
        msf = 2^(sum(rs.*log2.(ls))/sum(rs)) # weighted average
        maxr,maxi = findmax(fm)
        maxl = fl[maxi]

        fit = ()
        if isfit
            try
                # mfit = fitmodel(:sfdog,ls,rs) # fit Difference of Gaussians
                mfit = fitmodel(:sfgaussian,ls,rs) # fit Gaussian of logarithm sf
                fit = (;sftuningfeature(mfit)...,mfit)
            catch
                display.(stacktrace(catch_backtrace()))
            end
        end

        return (;up,msf,max=maxl=>maxr,fit)
    elseif factor == :ColorID
        # transform colorId to hue angle
        ucid = sort(unique(fl))
        hstep = 2pi/length(ucid)
        ha = map(l->hstep*(findfirst(c->c==l,ucid)-1),fl)
        # for hue direction
        hm, = circ_mean(ha,w=fr)
        oh = mod(rad2deg(angle(hm)),360)
        hcv = circ_var(ha,w=fr,d=hstep)
        maxi = argmax(fr)
        maxh = rad2deg(ha[maxi])
        maxr = fr[maxi]

        fit = ()
        if isfit
            # fit Generalized von Mises for hue
            try
                mfit = fitmodel(:gvm,θ,fr)
                fit = (circtuningfeature(mfit,od=π,fn=:h)...,gvm=mfit)
            catch
            end
            # fit von Mises for hue axis
            try
                mfit = fitmodel(:vmn2,θ,fr)
                fit = (fit...,circtuningfeature(mfit,od=0.5π,fn=:ha)...,vmn2=mfit)
            catch
            end
        end

        return (;ham,oha,hacv,hm,oh,hcv,maxh,maxr,fit)
    elseif factor in [:HueAngle,:Angle]
        θ = mod2pi.(deg2rad.(ls))
        α = mod2pi.(deg2rad.(fl))
        d = mean(diff(sort(unique(α)))) # angle spacing
        up, = circ_otest(α,w=fm) # Omnibus test for non-uniformity
        # for axis
        aα = mod.(α,π)
        acv, = circ_var(2aα,w=fm,d=2d)
        acm, = circ_mean(2aα,w=fm)
        acm = rad2deg(mod2pi(acm)/2)
        # for angle
        cv, = circ_var(α;w=fm,d)
        cm, = circ_mean(α,w=fm)
        cm = rad2deg(mod2pi(cm))
        maxr,maxi = findmax(fm)
        maxl = fl[maxi]

        fit = ()
        if isfit
            try
                mfit = fitmodel(:gvm,θ,rs) # fit Generalized von Mises
                fit = (;circtuningfeature(mfit,od=[π,0.5π],fn=:a)...,mfit)
            catch
                display.(stacktrace(catch_backtrace()))
            end
        end

        return (;up,acv,acm,cv,cm,max=maxl=>maxr,fit)
    else
        return ()
    end
end

"""
Spike Triggered Average of Images

1. x: Matrix where each row is one image
2. y: Vector of image response

- norm: normalization factor, default no normalization.
        it could be ``sum(y)`` if y is number of spike or spike rate, then STA would be spiking probability.
- whiten: whiten factor, default no whiten.
        it could be ``(xᵀx)⁻¹`` or inverse of covariance matrix, that decorrelate STA.
"""
function sta(x::AbstractMatrix,y::AbstractVector;norm=nothing,whiten=nothing)
    r = x'*y
    isnothing(norm) || (r/=norm)
    isnothing(whiten) || (r=whiten*r)

    # r = x'*x\r
    # r=length(y)*inv(cov(x,dims=1))*r

    return r
end

function psthsts(xs::Vector,binedges::Vector,c;israte::Bool=true,normfun=nothing)
    m,se,x = psth(xs,binedges,israte=israte,normfun=normfun)
    df = DataFrame(x=x,m=m,se=se,c=fill(c,length(x)))
end
function psthsts(xs::Vector,binedges::Vector,cond::DataFrame;israte::Bool=true,normfun=nothing)
    fs = finalfactor(cond)
    vcat([psth(xs[r[:i]],binedges,condstring(r,fs),israte=israte,normfun=normfun) for r in eachrow(cond)]...)
end
function psthsts(xs::Vector,binedges::Vector,ctc::DataFrame,factor;israte::Bool=true,normfun=nothing)
    vf = intersect(names(ctc),factor)
    isempty(vf) && error("No Valid Factor Found.")
    psth(xs,binedges,condin(ctc[:,vf]),israte=israte,normfun=normfun)
end
function psth(ds::DataFrame,binedges::Vector,conds::Vector;normfun=nothing,spike=:spike,isse::Bool=true)
    is,ss = findcond(ds,conds)
    df = psth(map(x->ds[spike][x],is),binedges,ss,normfun=normfun)
    if isse
        df[:ymin] = df[:y]-df[:ysd]./sqrt(df[:n])
        df[:ymax] = df[:y]+df[:ysd]./sqrt(df[:n])
    end
    return df,ss
end
function psthstss(xss::Vector,binedges::Vector,conds;normfun=nothing)
    n = length(xss)
    n!=length(conds) && error("Length of xss and conds don't match.")
    dfs = [psth(xss[i],binedges,conds[i],normfun=normfun) for i=1:n]
    return cat(1,dfs)
end
function unitdensity(pos;w=ones(length(pos)),spacerange=extrema(pos),bw=0.01(last(spacerange)-first(spacerange)),
                    step=bw/2,r=nothing,wfun=sum,s=nothing)
    hbw = bw/2
    y = first(spacerange):step:last(spacerange)
    n = [wfun(w[i-hbw .<=pos.< i+hbw]) for i in y]
    if !isnothing(r)
        n = n/(bw*π*r^2)
    end
    i = isnan.(n) .| isinf.(n)
    if any(i)
        n = Spline1D(y[.!i],n[.!i],k=3,bc="extrapolate")(y)
    end
    if !isnothing(s)
        g = gaussianf.(-5:5,σ=s)
        g ./= sum(g)
        n = imfilter(n,centered(g),ImageFiltering.Fill(0))
    end
    return (;n,y)
end
function spacepsth(unitpsth,unitposition;w=ones(size(unitposition,1)),spacerange=extrema(unitposition),
                    bw=0.01(last(spacerange)-first(spacerange)),step=bw/2)
    hbw = bw/2
    x = unitpsth[1].x
    y = first(spacerange):step:last(spacerange)

    n = zeros(length(y))
    psth = zeros(length(y),length(x))
    for i in eachindex(y)
        @views j = y[i]-hbw .<=unitposition[:,2].< y[i]+hbw
        n[i] = sum(w[j])
        if n[i] > 0
            @views psth[i,:] = mapreduce(u->u.m,.+,unitpsth[j])/n[i]
        end
    end
    return (;psth,x,y,n)
end

"""
Shift(shuffle) corrected, normalized(coincidence/spike), condition/trial-averaged Cross-Correlogram of binary spike trains.
(Bair, W., Zohary, E., and Newsome, W.T. (2001). Correlated Firing in Macaque Visual Area MT: Time Scales and Relationship to Behavior. J. Neurosci. 21, 1676–1697.)
"""
function correlogram(bst1,bst2;lag=nothing,isnorm=true,shiftcorrection=true,condis=nothing)
    if !isnothing(condis)
        cccg=[];x=[]
        for ci in condis
            ccg,x = correlogram(bst1[:,ci],bst2[:,ci],lag=lag,isnorm=isnorm,shiftcorrection=shiftcorrection,condis=nothing)
            push!(cccg,ccg)
        end
        ccg=dropdims(mean(hcat(cccg...),dims=2),dims=2)
        return ccg,x
    end
    n,nepoch = size(bst1)
    lag = floor(Int,isnothing(lag) ? min(n-1, 10*log10(n)) : lag)
    x = -lag:lag;xn=2lag+1
    cc = Array{Float64}(undef,xn,nepoch)
    for k in 1:nepoch
        cc[:,k]=crosscov(bst1[:,k],bst2[:,k],x,demean=false)*n
    end
    ccg = dropdims(mean(cc,dims=2),dims=2)
    if isnorm
        λ1 = mean(mean(bst1,dims=1))
        λ2 = mean(mean(bst2,dims=1))
        gmsr = sqrt(λ1*λ2)
        Θ = n.-abs.(x)
        normfactor = 1 ./ Θ ./ gmsr
    end
    if shiftcorrection
        psth1 = dropdims(mean(bst1,dims=2),dims=2)
        psth2 = dropdims(mean(bst2,dims=2),dims=2)
        s = crosscov(psth1,psth2,x,demean=false)*n
        shiftccg = (nepoch*s .- ccg)/(nepoch-1)
        if isnorm
            ccg .*= normfactor
            shiftccg .*= normfactor
        end
        ccg .-= shiftccg
    elseif isnorm
        ccg .*= normfactor
    end
    ccg,x
end
function circuitestimate(unitbinspike;lag=nothing,maxprojlag=3,minepoch=5,minspike=10,esdfactor=5,isdfactor=3.5,unitid=[],condis=nothing)
    nunit=length(unitbinspike)
    n,nepoch = size(unitbinspike[1])
    lag = floor(Int,isnothing(lag) ? min(n-1, 10*log10(n)) : lag)
    x = -lag:lag;xn=2lag+1

    isenoughspike = (i,j;minepoch=5,minspike=10) -> begin
        vsi = sum(i,dims=1)[:] .>= minspike
        vsj = sum(j,dims=1)[:] .>= minspike
        count(vsi) >= minepoch && count(vsj) >= minepoch
    end

    ccgs=[];ccgis=[];projs=[];eunits=[];iunits=[];projweights=[]
    for (i,j) in combinations(1:nunit,2)
        if isnothing(condis)
            !isenoughspike(unitbinspike[i],unitbinspike[j],minepoch=minepoch,minspike=minspike) && continue
            vcondis=nothing
        else
            vci = map(ci->isenoughspike(unitbinspike[i][:,ci],unitbinspike[j][:,ci],minepoch=minepoch,minspike=minspike),condis)
            all(.!vci) && continue
            vcondis = condis[vci]
        end

        ccg,_ = correlogram(unitbinspike[i],unitbinspike[j],lag=lag,condis=vcondis)
        ps,es,is,pws = projectionfromcorrelogram(ccg,i,j,maxprojlag=maxprojlag,esdfactor=esdfactor,isdfactor=isdfactor)
        if !isempty(ps)
            push!(ccgs,ccg);push!(ccgis,(i,j))
            append!(projs,ps);append!(eunits,es);append!(iunits,is);append!(projweights,pws)
        end
    end
    unique!(eunits);unique!(iunits)
    if length(unitid)==nunit
        map!(t->(unitid[t[1]],unitid[t[2]]),ccgis,ccgis)
        map!(t->(unitid[t[1]],unitid[t[2]]),projs,projs)
        map!(t->unitid[t],eunits,eunits)
        map!(t->unitid[t],iunits,iunits)
    end
    return ccgs,x,ccgis,projs,eunits,iunits,projweights
end
function projectionfromcorrelogram(cc,i,j;maxprojlag=3,minbaselag=maxprojlag+1,esdfactor=5,isdfactor=5)
    midi = Int((length(cc)+1)/2)
    base = [cc[midi+minbaselag:end];cc[1:midi-minbaselag]]
    bm,bsd = mean_and_std(base);hl = bm + esdfactor*bsd;ll = bm - isdfactor*bsd
    forwardlags = midi+1:midi+maxprojlag
    backwardlags = midi-maxprojlag:midi-1
    fcc = cc[forwardlags]
    bcc = cc[backwardlags]
    ps=[];ei=[];ii=[];pws=[]

    if ll <= cc[midi] <= hl
        if any(fcc .> hl)
            push!(ps,(i,j));push!(ei,i);push!(pws,(maximum(fcc)-bm)/bsd)
        elseif any(fcc .< ll)
            push!(ps,(i,j));push!(ii,i);push!(pws,(minimum(fcc)-bm)/bsd)
        end
        if any(bcc .> hl)
            push!(ps,(j,i));push!(ei,j);push!(pws,(maximum(bcc)-bm)/bsd)
        elseif any(bcc .< ll)
            push!(ps,(j,i));push!(ii,j);push!(pws,(minimum(bcc)-bm)/bsd)
        end
    end
    return ps,ei,ii,pws
end

"Check Layer Boundaries"
function checklayer!(ls::Dict;ln=["1", "2", "3", "23", "4A", "4B", "4AB", "4Cα", "4Cβ", "4C", "5", "6", "56", "WM"])
    n = length(ln)
    for i in 1:n-1
        if haskey(ls,ln[i])
            for j in (i+1):n
                if haskey(ls,ln[j])
                    ls[ln[i]][1] = ls[ln[j]][2]
                    break
                end
            end
        end
    end
    return ls
end

"""
Try to locate cell layer.

1. y coordinate of cell postion
2. layer definition
"""
function assignlayer(y,layer::Dict)
    l = missing
    for k in keys(layer)
        if layer[k][1] <= y < layer[k][2]
            l=k;break
        end
    end
    return l
end

"Check circuit consistency and remove duplicates"
function checkcircuit(projs,eunits,iunits,projweights)
    ivu = intersect(eunits,iunits)
    veunits = setdiff(eunits,ivu)
    viunits = setdiff(iunits,ivu)
    ivp = map(p->p[1] in ivu,projs)
    vprojs = projs[.!ivp]
    vprojweights = projweights[.!ivp]
    ui = indexin(unique(vprojs),vprojs)
    return vprojs[ui],veunits,viunits,vprojweights[ui]
end