text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

export HDF5DataLayer

type HDF5DataEnsemble <: DataEnsemble
    name         :: Symbol
    neurons      :: Array{DataNeuron}
    train_id     :: Cint
    test_id      :: Cint
    train_epoch  :: Int
    test_epoch   :: Int
    net_subgroup :: Cint
    connections  :: Vector{Connection}
    phase        :: Phase
end

function HDF5DataEnsemble(name::Symbol, neurons::Array{DataNeuron}, train_id::Cint, test_id::Cint)
    HDF5DataEnsemble(name, neurons, train_id, test_id, 1, 1, 1, Connection[], TrainTest)
end

@eval function HDF5DataEnsemble(net::Net, train_id::Cint, test_id::Cint, target::Symbol)
    if target == :data
        ndim = ccall((:get_data_ndim, $libIO), Cint, (Cint,), train_id)
        _shape = ccall((:get_data_shape, $libIO), Ptr{Cint}, (Cint,), train_id)
    else
        ndim = ccall((:get_label_ndim, $libIO), Cint, (Cint,), train_id)
        _shape = ccall((:get_label_shape, $libIO), Ptr{Cint}, (Cint,), train_id)
    end
    # first index is batch so skip, reverse c order
    shape = pointer_to_array(_shape, ndim)[end:-1:2]
    neurons = Array(DataNeuron, shape...)
    for i = 1:length(neurons)
        neurons[i] = DataNeuron(0.0f0)
    end
    ens = HDF5DataEnsemble(target, neurons, train_id, test_id)
    add_ensemble(net, ens)
    ens
end

@eval function init(ens::HDF5DataEnsemble, net::Net)
    arr = Array(Float32, size(ens)..., net.batch_size)
    set_buffer(net, symbol(ens.name,:value), arr; _copy=false)
    set_buffer(net, symbol(ens.name,:∇), zeros(arr); _copy=false)
    if ens.name == :data
        ccall((:set_data_pointer, $libIO), Void, (Cint, Ptr{Float32},), ens.train_id, arr)
        ccall((:set_data_pointer, $libIO), Void, (Cint, Ptr{Float32},), ens.test_id, arr)
    else
        ccall((:set_label_pointer, $libIO), Void, (Cint, Ptr{Float32},), ens.train_id, arr)
        ccall((:set_label_pointer, $libIO), Void, (Cint, Ptr{Float32},), ens.test_id, arr)
    end
end

@eval function forward{N}(ens::HDF5DataEnsemble, data::Array{Float32,N}, net::Net, phase::Phase)
    if phase == Train
        id = ens.train_id
    else
        id = ens.test_id
    end
    epoch = ccall((:get_epoch, $libIO), Cint, (Cint,), id)
    if phase == Train
        net.train_epoch = epoch
    else
        net.test_epoch = epoch
    end
    if ens.name == :data
        _shape = ccall((:get_next_batch, $libIO), Void, (Cint,), id)
    end
end

function parse_hdf5_source(source::AbstractString)
    open(source, "r") do s
        strip(readlines(s)[1])
    end
end

@eval function HDF5DataLayer(net::Net, train_data_source::AbstractString,
                       test_data_source::AbstractString;
                       shuffle=true, scale=1.0f0)
    batch_size = net.batch_size
    @assert(batch_size > 0, "Data Layer batch_size must be greater than 0")
    train_data_source = parse_hdf5_source(train_data_source)
    test_data_source = parse_hdf5_source(test_data_source)
    train_id = ccall((:init_dataset, $libIO), Cint, (Cint, Ptr{UInt8}, Cuchar, Cuchar, Cuchar), batch_size, train_data_source, shuffle, LATTE_MPI, false)
    test_id = ccall((:init_dataset, $libIO), Cint, (Cint, Ptr{UInt8}, Cuchar, Cuchar, Cuchar), batch_size, test_data_source, false, LATTE_MPI, true)
    HDF5DataEnsemble(net, train_id, test_id, :data), HDF5DataEnsemble(net, train_id, test_id, :label)
end


text:
# negative of obj, since the original formulation is max_\theta
# whereas the optimizer only support minimization

# all the opt in the metric needs to be done in cpu

# objectve: cuda vector
function ap_objective(ps::CuArrays.CuVector, y::CuArrays.CuVector, pm::PerformanceMetric; args...)
    psc = CuArrays.collect(ps)    # to cpu
    yc = CuArrays.collect(y)      # to cpu

    obj, _ = objective(pm, psc, yc; args...)
    obj = obj + dot(psc, yc)
    return -obj
end

# custom gradient
@adjoint function ap_objective(ps::CuArrays.CuVector, y::CuArrays.CuVector, pm::PerformanceMetric; args...)
    psc = CuArrays.collect(ps)    # to cpu
    yc = CuArrays.collect(y)      # to cpu

    obj, q = objective(pm, psc, yc; args...)
    obj = obj + dot(psc, yc)

    grad = CuArrays.cu(q - yc)   # to cuda
    return -obj, Δ -> (Δ * grad, nothing, nothing)
end

ext: jl


text:
# # [Multirate Runge-Kutta Timestepping](@id Multirate-RK-Timestepping)

# In this tutorial, we shall explore the use of Multirate Runge-Kutta
# methods for the solution of nonautonomous (or non time-invariant) equations.
# For our model problem, we shall reuse the acoustic wave test in the GCM
# configuration. See its [code](@ref Acoustic-Wave-Configuration)
# for details on the model and parameters. For the purposes of this tutorial,
# we will only run the experiment for a total of 3600 simulation seconds.
# Details on this test case can be found in Sec. 4.3 of [Giraldo2013](@cite).

using ClimateMachine
const clima_dir = dirname(dirname(pathof(ClimateMachine)));
include(joinpath(
    clima_dir,
    "tutorials",
    "Numerics",
    "TimeStepping",
    "tutorial_acousticwave_config.jl",
))

FT = Float64;

# The typical context for Multirate splitting methods is given by problems
# in which the tendency ``\mathcal{T}`` is assumed to have single parts that
# operate on different time rates (such as a slow time scale and a fast time scale).
# A general form is given by
#
# ``
# 	\dot{\boldsymbol{q}} = \mathcal{T}(\boldsymbol{q}) \equiv
# 	{\mathcal{T}}_{f}(\boldsymbol{q}) + {\mathcal{T}}_{s}(\boldsymbol{q}),
# ``
#
# where the right-hand side has been split into a "fast" component ``{\mathcal{T}_{f}}``,
# and a "slow" component ``{\mathcal{T}_{s}}``.

# Referencing the canonical form introduced in [Time integration](@ref
# Time-integration), both ``{\mathcal{T}_{f}}`` and ``{\mathcal{T}_{s}}``
# could be discretized either explicitly or implicitly, hence, they could
# belong to either ``\mathcal{F}(t, \boldsymbol{q})`` or ``\mathcal{G}(t, \boldsymbol{q})``
# term.
#
# For a given time-step size ``\Delta t``, the two-rate method in [Schlegel2009](@cite)
# is summarized as the following:
#
# ```math
# \begin{align}
#     \boldsymbol{Q}_1 &= \boldsymbol{q}(t_n), \\
#     \boldsymbol{r}_{i} &= \sum_{j=1}^{i-1}\tilde{a}^O_{ij} {\mathcal{T}_{s}}(\boldsymbol{Q}_{j}), \\
#     \boldsymbol{w}_{i,1} &= \boldsymbol{Q}_{i-1},\\
#     \boldsymbol{w}_{i,k} &= \boldsymbol{w}_{i,k-1} + \Delta t \tilde{c}_i^O \sum_{j=1}^{k-1}\tilde{a}^I_{k,j}
#     \left(\frac{1}{\tilde{c}_i^O}\boldsymbol{r}_i + {\mathcal{T}_{f}}(\boldsymbol{w}_{i,j})\right),\\
#     & \quad\quad i = 2, \cdots, s^O + 1 \text{ and } k = 2, \cdots, s^I + 1,\nonumber \\
#     \boldsymbol{Q}_i &= \boldsymbol{w}_{i,s^I + 1}
# \end{align}
# ```
#
# where the tilde parameters denote increments per RK stage:
#
# ```math
# \begin{align}
#     \tilde{a}_{ij} &= \begin{cases}
#         a_{i,j} - a_{i-1, j} & \text{if } i < s + 1 \\
#         b_j - a_{s,j} & \text{if } i = s + 1
#     \end{cases},\\
#     \tilde{c}_{i} &= \begin{cases}
#         c_{i} - c_{i-1} & \text{if } i < s + 1 \\
#         1 - c_{s} & \text{if } i = s + 1
#     \end{cases},
# \end{align}
# ```
#
# where the coefficients ``a``, ``b``, and ``c`` correspond to the Butcher
# tableau for a given RK method. The superscripts ``O`` and ``I`` denote the
# *outer* (slow) and *inner* (fast) components of the multirate method
# respectively. Thus, tilde coefficients should be associated with the RK
# method indicated by the superscripts. In other words, the RK methods
# for the slow ``{\mathcal{T}_{s}}`` and fast
# ``{\mathcal{T}_{f}}`` components have Butcher tables given by:
#
# ```math
# \begin{align}
#     \begin{array}{c|c}
#     \boldsymbol{c}_{O} &\boldsymbol{A}_{O} \\
#     \hline
#     & \boldsymbol{b}_O^T
#     \end{array}, \quad
#     \begin{array}{c|c}
#     \boldsymbol{c}_{I} &\boldsymbol{A}_{I} \\
#     \hline
#     & \boldsymbol{b}_I^T
#     \end{array},
# \end{align}
# ```
#
# where ``\boldsymbol{A}_O = \lbrace a_{i,j}^O\rbrace``, ``\boldsymbol{b}_O = \lbrace b_i^O \rbrace``, and
# ``c_O = \lbrace c_i^O \rbrace`` (similarly for ``\boldsymbol{A}_I``, ``\boldsymbol{b}_I``, and ``\boldsymbol{c}_I``).
# The method described here is for an explicit RK outer method with ``s`` stages.
# More details can be found in [Schlegel2012](@cite).

# The acoustic wave test case used in this tutorial represents a global-scale
# problem with inertia-gravity waves traveling around the entire planet.
# It has a hydrostatically balanced initial state that is given a pressure
# perturbation.
# This initial pressure perturbation causes an acoustic wave to travel to
# the antipode, coalesce, and return to the initial position. The exact solution
# of this test case is simple in that the (linear) acoustic theory allows one
# to verify the analytic speed of sound based on the thermodynamics variables.
# The initial condition is defined as a hydrostatically balanced atmosphere
# with background (reference) potential temperature.

ode_solver = ClimateMachine.MultirateSolverType(
    splitting_type = ClimateMachine.HEVISplitting(),
    slow_method = LSRK54CarpenterKennedy,
    fast_method = ARK2GiraldoKellyConstantinescu,
    implicit_solver_adjustable = true,
    timestep_ratio = 100,
)

timeend = FT(3600)
CFL = FT(5)
cfl_direction = HorizontalDirection()
run_acousticwave(ode_solver, CFL, cfl_direction, timeend);

# The interested reader can explore the combination of different slow and
# fast methods for Multirate solvers, consulting the ones available in
# `ClimateMachine.jl`, such as the
# [`Low-Storage-Runge-Kutta-methods`](@ref ClimateMachine.ODESolvers.LowStorageRungeKutta2N),
# [`Strong-Stability-Preserving-RungeKutta-methods`](@ref ClimateMachine.ODESolvers.StrongStabilityPreservingRungeKutta),
# and [`Additive-Runge-Kutta-methods`](@ref ClimateMachine.ODESolvers.AdditiveRungeKutta).

# ## References
# - [Giraldo2013](@cite)
# - [Schlegel2009](@cite)
# - [Schlegel2012](@cite)


text:
using LinearAlgebra
using DataFrames
using CSV
using Plots

obs = CSV.File("../../DATA/ovariancancer_obs.csv", header = false) |> Tables.matrix
grp = CSV.File("../../DATA/ovariancancer_grp.csv", header = false) |> Tables.matrix

U, S, V = svd(obs)

p1 = plot(S, yaxis = :log, legend = false)
p2 = plot(cumsum(S)./sum(S), legend = false)

cancer = [[],[],[]]
normal = [[],[],[]]

for i = 1:size(obs)[1]
    x = V[:,1]'*obs[i,:]
    y = V[:,2]'*obs[i,:]
    z = V[:,3]'*obs[i,:]
    if (grp[i] == "Cancer")
        append!.(cancer, [x, y, z])
    else
        append!.(normal, [x, y, z])
    end
end

p3 = scatter(cancer..., mc = :red, label = "Cancer")
scatter!(normal..., mc = :blue, label = "Normal")


text:
# Line geometry connects (x, y) coordinates with lines.

struct LineGeometry <: Gadfly.GeometryElement
    default_statistic::Gadfly.StatisticElement

    # Do not reorder points along the x-axis.
    preserve_order::Bool

    order::Int

    tag::Symbol
end

function LineGeometry(default_statistic=Gadfly.Stat.identity();
                      preserve_order=false, order=2, tag=empty_tag)
    LineGeometry(default_statistic, preserve_order, order, tag)
end

"""
    Geom.line[(; preserve_order=false, order=2)]

Draw a line connecting the `x` and `y` coordinates.  Optionally plot multiple
lines according to the `group`, `color` or `linestyle` aesthetics.  `order` controls whether
the lines(s) are underneath or on top of other forms.

Set `preserve_order=true` to *not* sort the points according to their
position along the x axis, or use the equivalent [`Geom.path`](@ref) alias.
"""
const line = LineGeometry

### why would one ever want to set preserve_order to false here
"""
    Geom.contours[(; levels=15, samples=150, preserve_order=true)]

Draw contour lines of the 2D function, matrix, or DataFrame in the `z`
aesthetic.  This geometry is equivalent to [`Geom.line`](@ref) with
[`Stat.contour`](@ref); see the latter for more information.
"""
function contour(; levels=15, samples=150, preserve_order=true)
    return LineGeometry(Gadfly.Stat.contour(levels=levels, samples=samples),
                                            preserve_order=preserve_order)
end

# Only allowing identity statistic in paths b/c I don't think any
# any of the others will work with `preserve_order=true` right now
"""
    Geom.path

Draw lines between `x` and `y` points in the order they are listed.  This
geometry is equivalent to [`Geom.line`](@ref) with `preserve_order=true`.
"""
path() = LineGeometry(preserve_order=true)

"""
    Geom.density[(; bandwidth=-Inf)]

Draw a line showing the density estimate of the `x` aesthetic.
This geometry is equivalent to [`Geom.line`](@ref) with
[`Stat.density`](@ref); see the latter for more information.
"""
density(; bandwidth::Real=-Inf) =
    LineGeometry(Gadfly.Stat.density(bandwidth=bandwidth))

"""
    Geom.density2d[(; bandwidth=(-Inf,-Inf), levels=15)]

Draw a set of contours showing the density estimate of the `x` and `y`
aesthetics.  This geometry is equivalent to [`Geom.line`](@ref) with
[`Stat.density2d`](@ref); see the latter for more information.
"""
density2d(; bandwidth::Tuple{Real,Real}=(-Inf,-Inf), levels=15) =
    LineGeometry(Gadfly.Stat.density2d(bandwidth=bandwidth, levels=levels); preserve_order=true)

"""
    Geom.smooth[(; method:loess, smoothing=0.75)]

Plot a smooth function estimated from the line described by `x` and `y`
aesthetics.  Optionally group by `color` and plot multiple independent smooth
lines.  This geometry is equivalent to [`Geom.line`](@ref) with
[`Stat.smooth`](@ref); see the latter for more information.
"""
smooth(; method::Symbol=:loess, smoothing::Float64=0.75) =
    LineGeometry(Gadfly.Stat.smooth(method=method, smoothing=smoothing),
    order=5)

"""
    Geom.step[(; direction=:hv)]

Connect points described by the `x` and `y` aesthetics using a stepwise
function.  Optionally group by `color` or `group`.  This geometry is equivalent
to [`Geom.line`](@ref) with [`Stat.step`](@ref); see the latter for more
information.
"""
step(; direction::Symbol=:hv) = LineGeometry(Gadfly.Stat.step(direction=direction))

default_statistic(geom::LineGeometry) = geom.default_statistic

element_aesthetics(::LineGeometry) = [:x, :y, :color, :group, :linestyle]

function Gadfly.Geom.render(geom::LineGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.line", aes, :x, :y)
    Gadfly.assert_aesthetics_equal_length("Geom.line", aes, :x, :y)

    default_aes = Gadfly.Aesthetics()
    default_aes.group = IndirectArray([1])
    default_aes.color = [theme.default_color]
    default_aes.linestyle = theme.line_style[1:1]
    aes = inherit(aes, default_aes)

    # Render lines, using multivariate groupings:
    XT, YT = eltype(aes.x), eltype(aes.y)
    GT, CT, LST = Int, eltype(aes.color), eltype(aes.linestyle)
    groups = collect(Tuple{GT, CT, LST}, Compose.cyclezip(aes.group, aes.color, aes.linestyle))
    ugroups = unique(groups)
    nugroups = length(ugroups)
    # Recycle groups
    (1 .< length(groups) .< length(aes.x))  && (groups = [b for (a, b) in zip(aes.x, cycle(groups))])

    # Point order
    aes_x, aes_y, zgroups  = if !geom.preserve_order
        p = sortperm(aes.x)
        aes.x[p], aes.y[p], (nugroups==1 ? ugroups : groups[p])
    else
        aes.x, aes.y, (nugroups==1 ? ugroups : groups)
    end

    gs = Vector{GT}(undef, nugroups)
    cs = Vector{CT}(undef, nugroups)
    lss = Vector{LST}(undef, nugroups)
    lines = Vector{Vector{Tuple{XT,YT}}}(undef, nugroups)
    linestyle_palette_length = length(theme.line_style)
    if nugroups==1
        gs[1], cs[1], lss[1] = zgroups[1]
        lines[1] = collect(Tuple{XT, YT}, zip(aes_x, aes_y))
    elseif nugroups>1
        for (k,g) in enumerate(ugroups)
            i = zgroups.==[g]
            gs[k], cs[k], lss[k] = g
            lines[k] = collect(Tuple{XT,YT}, zip(aes_x[i], aes_y[i]))
        end
    end

    linestyles = Gadfly.get_stroke_vector.(LST<:Int ?
         theme.line_style[mod1.(lss, linestyle_palette_length)] : lss)

    classes = svg_color_class_from_label.(aes.color_label(cs))
    ctx = context(order=geom.order)
    ctx = compose!(ctx, (context(), Compose.line(lines, geom.tag),
            stroke(cs), strokedash(linestyles),
            svgclass(classes)), svgclass("geometry"))

    return compose!(ctx, fill(nothing), linewidth(theme.line_width))
end


text:
export switchings

mutable struct DiscreteSwitchingSequence{S <: HybridSystem, MT <: AbstractMatrix, TT}
    s::S
    A::MT
    seq::Vector{TT}
    len::Int
end

function Base.show(io::IO, s::DiscreteSwitchingSequence)
    println(io, "Discrete switching sequence of length $(s.len):")
    println(io, s.seq)
end

HybridSystems.source(s::HybridSystem, seq::DiscreteSwitchingSequence) = source(s, seq.seq[1])
HybridSystems.target(s::HybridSystem, seq::DiscreteSwitchingSequence) = target(s, seq.seq[end])
# Short circuit for unconstrained system
HybridSystems.source(::HybridSystem{OneStateAutomaton}, ::DiscreteSwitchingSequence) = 1
HybridSystems.target(::HybridSystem{OneStateAutomaton}, ::DiscreteSwitchingSequence) = 1

function DiscreteSwitchingSequence(s::HybridSystem, A::AbstractMatrix, seq::Vector)
    DiscreteSwitchingSequence(s, A, seq, length(seq))
end
switchingsequence(s::HybridSystem, A::AbstractMatrix, seq::Vector) = DiscreteSwitchingSequence(s, A, seq)
using LinearAlgebra
function switchingsequence(s::HybridSystem, len::Int=0, v::Int=1)
    n = statedim(s, v)
    A = Matrix(1.0LinearAlgebra.I, n, n)
    DiscreteSwitchingSequence(s, A, Vector{transitiontype(s)}(undef, len), 0)
end

function Base.prepend!(s::DiscreteSwitchingSequence, other::DiscreteSwitchingSequence)
    s.A = s.A * other.A
    if s.len < length(s.seq)
        if s.len + other.len <= length(s.seq)
            s.seq[other.len .+ (1:s.len)] = s.seq[1:s.len] # Cannot use view here
            s.seq[1:other.len] = @view other.seq[1:other.len]
        else
            len = length(s.seq)
            off = max(0, len - other.len)
            append!(s.seq, @view s.seq[(off+1):s.len])
            s.seq[other.len .+ (1:off)] = s.seq[1:off] # Cannot use view here
            if len < other.len
                s.seq[:] = @view other.seq[1:len]
                append!(s.seq, @view other.seq[(len+1):other.len])
            else
                s.seq[1:other.len] = @view other.seq[1:other.len]
            end
        end
    else
        prepend!(s.seq, other.seq)
    end
    s.len += other.len
end

function Base.append!(s::DiscreteSwitchingSequence, other::DiscreteSwitchingSequence)
    s.A = other.A * s.A
    if s.len < length(s.seq)
        if s.len + other.len <= length(s.seq)
            s.seq[(s.len+1):(s.len+other.len)] = @view other.seq[1:other.len]
        else
            off = length(s.seq) - s.len
            s.seq[(s.len+1):end] = @view other.seq[1:off]
            append!(s.seq, @view other.seq[(off+1):other.len])
        end
    else
        append!(s.seq, other.seq)
    end
    s.len += other.len
end

struct SwitchingIterator{S<:HybridSystem}
    s::S
    k::Int        # length of sequence
    v0::Int       # starting mode
    forward::Bool # Is the sequence going forward or backward
end
Base.IteratorSize(::SwitchingIterator) = Base.SizeUnknown()

"""
    switchings(s::HybridSystem, k::Int, v0::Int, forward::Bool=true)

Iterates over all the `forward` switching of length `k` starting at mode `v0`.
"""
function switchings(s::HybridSystem, k::Int, v0::Int, forward::Bool=true)
    SwitchingIterator(s, k, v0, forward)
end

struct SwitchingIteratorState{IT, ST, MT, ET}
    # modeit[i] is a list of all the possible transitions for the (i-1)th mode
    modeit::Vector{IT}
    # modest[i] is the ith state of iterator modeit[i]
    modest::Vector{ST}
    As::Vector{MT}
    seq::Vector{ET}
    function SwitchingIteratorState{IT, ST, MT, ET}(k::Int) where {IT, ST, MT, ET}
        new{IT, ST, MT, ET}(Vector{IT}(undef, k), Vector{ST}(undef, k),
                            Vector{MT}(undef, k), Vector{ET}(undef, k))
    end
end

function io_transitions(s::HybridSystem, mode, forward::Bool)
    if forward
        return out_transitions(s, mode)
    else
        return in_transitions(s, mode)
    end
end
function io_state(s::HybridSystem, t, forward::Bool)
    if forward
        return target(s, t)
    else
        return source(s, t)
    end
end

function Base.iterate(it::SwitchingIterator)
    ts = io_transitions(it.s, it.v0, it.forward)
    IT = typeof(ts)
    state_item = iterate(ts)
    if state_item === nothing
        return nothing
    end
    ST = typeof(state_item[2])
    MT = typeof(first(it.s.resetmaps).A)
    ET = transitiontype(it.s)
    st = SwitchingIteratorState{IT, ST, MT, ET}(it.k)
    return complete_switching(it, st, it.forward ? 1 : it.k)
end
function Base.iterate(it::SwitchingIterator, st::SwitchingIteratorState)
    next_switching(it, st, it.forward ? it.k : 1)
end
function cur_mode(it::SwitchingIterator, st::SwitchingIteratorState, i::Int)
    j = i + (it.forward ? -1 : 1)
    if j <= 0 || j > it.k
        return it.v0
    else
        return io_state(it.s, st.seq[j], it.forward)
    end
end
function prev_matrix(it::SwitchingIterator, st::SwitchingIteratorState, i::Int)
    j = i + (it.forward ? -1 : 1)
    if j <= 0 || j > it.k
        return I
    else
        return st.As[j]
    end
end
function process_item_state(it::SwitchingIterator, st::SwitchingIteratorState, i::Int, item_state::Nothing)
    inc = it.forward ? 1 : -1
    return next_switching(it, st, i - inc)
end
function process_item_state(it::SwitchingIterator, st::SwitchingIteratorState, i::Int, item_state)
    st.modest[i] = item_state[2]
    st.seq[i] = item_state[1]
    B = resetmap(it.s, st.seq[i]).A
    A = prev_matrix(it, st, i)
    st.As[i] = it.forward ? B * A : A * B
    return complete_switching(it, st, i + (it.forward ? 1 : -1))
end
function next_switching(it::SwitchingIterator, st::SwitchingIteratorState, i::Int)
    inc = it.forward ? 1 : -1
    if i <= 0 || i > it.k
        return nothing
    else
        item_state = iterate(st.modeit[i], st.modest[i])
        process_item_state(it, st, i, item_state)
    end
end
function complete_switching(it::SwitchingIterator, st::SwitchingIteratorState, i::Int)
    inc = it.forward ? 1 : -1
    if i <= 0 || i > it.k
        return switchingsequence(it.s, st.As[i - inc], copy(st.seq)), st
    else
        st.modeit[i] = io_transitions(it.s, cur_mode(it, st, i), it.forward)
        item_state = iterate(st.modeit[i])
        process_item_state(it, st, i, item_state)
    end
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

contraction analysis example adapted from
"Stability and robustness analysis of nonlinear systems via contraction metrics
and SOS programming" by Aylward, E.M., Parrilo, P.A. and Slotine, J.J.E
=#

import DynamicPolynomials
const DP = DynamicPolynomials
import Combinatorics
import PolyJuMP
import SumOfSquares

struct ContractionJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    beta::Real
    M_deg::Int
    delta::Real
    use_matrixwsos::Bool # use wsos matrix cone, else PSD formulation
    is_feas::Bool
end

function build(inst::ContractionJuMP{T}) where {T <: Float64}
    delta = inst.delta
    n = 2
    dom = PolyUtils.FreeDomain{T}(n)

    M_halfdeg = div(inst.M_deg + 1, 2)
    (U_M, pts_M, Ps_M) = PolyUtils.interpolate(dom, M_halfdeg)
    lagrange_polys = get_lagrange_polys(pts_M, 2 * M_halfdeg)
    x = DP.variables(lagrange_polys)

    # dynamics according to the Moore-Greitzer model
    dx1dt = -x[2] - 1.5 * x[1]^2 - 0.5 * x[1]^3
    dx2dt = 3 * x[1] - x[2]
    dynamics = [dx1dt; dx2dt]

    model = JuMP.Model()
    JuMP.@variable(model, polys[1:3], PolyJuMP.Poly(
        PolyJuMP.MultivariateBases.FixedPolynomialBasis(lagrange_polys)))

    M = [polys[1] polys[2]; polys[2] polys[3]]
    dMdt = [dot(DP.differentiate(M[i, j], x), dynamics) for i in 1:n, j in 1:n]
    dfdx = DP.differentiate(dynamics, x)'
    Mdfdx = [sum(M[i, k] * dfdx[k, j] for k in 1:n) for i in 1:n, j in 1:n]
    R = Mdfdx + Mdfdx' + dMdt + inst.beta * M

    if inst.use_matrixwsos
        deg_R = maximum(DP.maxdegree.(R))
        d_R = div(deg_R + 1, 2)
        (U_R, pts_R, Ps_R) = PolyUtils.interpolate(dom, d_R)
        M_gap = [M[i, j](pts_M[u, :]) - (i == j ? delta : 0.0)
            for i in 1:n for j in 1:i for u in 1:U_M]
        R_gap = [-R[i, j](pts_R[u, :]) - (i == j ? delta : 0.0)
            for i in 1:n for j in 1:i for u in 1:U_R]

        wsosmatT = Hypatia.WSOSInterpPosSemidefTriCone{T}
        rt2 = sqrt(2)
        M_scal = Cones.scale_svec!(M_gap, rt2, incr = U_M)
        R_scal = Cones.scale_svec!(R_gap, rt2, incr = U_R)
        JuMP.@constraint(model, M_scal in wsosmatT(n, U_M, Ps_M))
        JuMP.@constraint(model, R_scal in wsosmatT(n, U_R, Ps_R))
    else
        PolyJuMP.setpolymodule!(model, SumOfSquares)
        JuMP.@constraint(model, M - Matrix(delta * I, n, n) in JuMP.PSDCone())
        JuMP.@constraint(model, -R - Matrix(delta * I, n, n) in JuMP.PSDCone())
    end

    return model
end

function test_extra(inst::ContractionJuMP{T}, model::JuMP.Model) where T
    stat = JuMP.termination_status(model)
    @test stat == (inst.is_feas ? MOI.OPTIMAL : MOI.INFEASIBLE)
    return
end

# TODO below functions used for recovering lagrange polys are not numerically
# stable for high degree (since they essentially calculate 2^degree);
# may be unnecessary if already doing a QR of the U*U Vandermonde,
# so better to replace them with linear algebra without using DP
function get_lagrange_polys(pts::Matrix{T}, deg::Int) where {T <: Real}
    if deg > 8
        @warn("get_lagrange_polys is not numerically stable for large degree")
    end
    (U, n) = size(pts)
    DP.@polyvar x[1:n]
    basis = get_chebyshev_polys(x, deg)
    @assert length(basis) == U
    vand_inv = inv([basis[j](x => view(pts, i, :)) for i in 1:U, j in 1:U])
    lagrange_polys = [dot(view(vand_inv, :, i), basis) for i in 1:U]
    return lagrange_polys
end

# returns the multivariate Chebyshev polynomials in x up to degree deg
function get_chebyshev_polys(x::Vector{DP.PolyVar{true}}, deg::Int)
    if deg > 8
        @warn("get_chebyshev_polys is not numerically stable for large degree")
    end
    n = length(x)
    u = get_chebyshev_univ(x, deg)
    V = Vector{DP.Polynomial{true, Int}}(undef, PolyUtils.get_L(n, deg))
    V[1] = DP.Monomial(1)
    col = 1
    for t in 1:deg, xp in Combinatorics.multiexponents(n, t)
        col += 1
        V[col] = u[1][xp[1] + 1]
        for j in 2:n
            V[col] *= u[j][xp[j] + 1]
        end
    end
    return V
end

function get_chebyshev_univ(monovec::Vector{DP.PolyVar{true}}, deg::Int)
    if deg > 8
        @warn("get_chebyshev_univ is not numerically stable for large degree")
    end
    n = length(monovec)
    u = Vector{Vector}(undef, n)
    for j in 1:n
        uj = u[j] = Vector{DP.Polynomial{true, Int}}(undef, deg + 1)
        uj[1] = DP.Monomial(1)
        uj[2] = monovec[j]
        for t in 3:(deg + 1)
            uj[t] = 2 * uj[2] * uj[t - 1] - uj[t - 2]
        end
    end
    return u
end


text:
__precompile__()

module GDAL

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("GDAL not installed properly, run Pkg.build(\"GDAL\"), restart Julia and try again")
end
include(depsjl_path)

module C
    import GDAL: libgdal

    include("C/misc.jl")
    include("C/common.jl")
    include("C/cpl_conv.jl")
    include("C/cpl_error.jl")
    include("C/cpl_minixml.jl")
    include("C/cpl_progress.jl")
    include("C/cpl_virtualmem.jl")
    include("C/cpl_vsi.jl")
    include("C/gdal.jl")
    include("C/gdal_alg.jl")
    include("C/gdal_vrt.jl")
    include("C/gdal_utils.jl")
    include("C/gdalwarper.jl")
    include("C/ogr_api.jl")
    include("C/ogr_core.jl")
    include("C/ogr_srs_api.jl")
end

include("C/misc.jl")
include("types.jl")
include("common.jl")
include("gdal_h.jl")
include("gdal_alg.jl")
include("gdal_vrt.jl")
include("gdal_utils.jl")
include("ogr_api.jl")
include("ogr_srs_api.jl")
include("cpl_error.jl")

include("error.jl")

function __init__()
    h = Libdl.dlopen(libgdal)
    finfo = Libdl.dlsym(h, :GDALVersionInfo)
    fconf = Libdl.dlsym(h, :CPLSetConfigOption)
    ferrh = Libdl.dlsym(h, :CPLSetErrorHandler)

    # Always check your dependencies from `deps.jl`
    check_deps()

    # register custom error handler
    funcptr = cfunction(gdaljl_errorhandler, Ptr{Void}, (CPLErr, Cint, Cstring))
    ccall(ferrh, Ptr{Void}, (Ptr{Void},), funcptr)

    # get GDAL version number
    versionstring = unsafe_string(ccall(finfo, Cstring, (Cstring,), "RELEASE_NAME"))
    global const GDALVERSION = convert(VersionNumber, versionstring)

    # set the GDAL_DATA variable
    global const GDAL_DATA = abspath(@__DIR__, "..", "deps", "usr", "share", "gdal")
    ENV["GDAL_DATA"] = GDAL_DATA
    ccall(fconf, Void, (Cstring, Cstring), "GDAL_DATA", GDAL_DATA)
end

"""
Load a null-terminated list of strings

That is it expects a `StringList`, in the sense of the CPL functions,
as a null-terminated array of strings.
"""
function unsafe_loadstringlist(ptr::Ptr{Cstring})
    strings = Vector{String}()
    (ptr == C_NULL) && return strings
    i = 1
    cstring = unsafe_load(ptr, i)
    while cstring != C_NULL
        push!(strings, unsafe_string(cstring))
        i += 1
        cstring = unsafe_load(ptr, i)
    end
    strings
end

end # module


text:
global charNum = parse(Int, ARGS[1])
##
using Pkg
Pkg.activate(".")
Pkg.instantiate()
##

using CSV
using DataFrames
using Pipe
using Statistics
using ProgressMeter
using LinearAlgebra
using StatsFuns
using MCPhylo
using Random
Random.seed!(6046749401947475242)

##

wals = CSV.read("../data/charMtx.csv", DataFrame)

famFreqs = combine(groupby(wals, :glot_fam), nrow)
sort!(famFreqs, :nrow, rev=true)

families = filter(x -> x.nrow>1, famFreqs).glot_fam
isolates = filter(x -> x.nrow==1, famFreqs).glot_fam

lineages = vcat(families, isolates)

isoDict = @pipe wals |>
      filter(x -> x.glot_fam ∈ isolates, _) |>
      zip(_.glot_fam, _.longname) |>
      Dict

##

d = CSV.read("../data/fpairMtx.csv", DataFrame)
d = innerjoin(d, wals[:,[:longname, :glot_fam]], on=:taxon => :longname)
##

function myvar(x)
    if length(x) == 1
        return 0.
    end
    var(x, corrected=false)
end

##

fm2trees = Dict()
@showprogress for fm in families
      fm2trees[fm] = MCPhylo.ParseNewick("../data/posteriorTrees/$fm.posterior.tree")
end

for fm in isolates
      fm2trees[fm] = repeat([Node(isoDict[fm])], 1000)
end

##

offDiag = [CartesianIndex(x,y) for y in 1:4 for x in 1:4 if x!=y];
taxa = d.taxon;

defined = findall(d[:,charNum+1] .!= "-")

##

nSims = 1000

sim = read("modelFitting/output/universal_$(lpad(charNum, 2, "0")).jls", ModelChains)

samples = vcat([sim.value[:,:,k] for k in sim.chains]...)

##
stats = []
@showprogress for g in 1:nSims
    i = rand(1:size(samples, 1))
    treeI = ceil(Int,invlogit(samples[i,1])*1000)

    simChar = Dict{String,Int}()

    fullSrates = exp.(samples[i,2:(end-1)])
    q = zeros(4, 4)
    q[offDiag] = fullSrates
    for i = 1:4
        q[i, i] -= sum(q[i, :])
    end

    D, U = eigen(q)
    Uinv = inv(U)

    eq = real.(Uinv[4, :])
    eq /= sum(eq)

    for fmI = 1:length(lineages)
        fmI
        tree = fm2trees[lineages[fmI]][treeI]
        preo = pre_order(tree)
        tStates = Vector{Int}(undef, length(preo))
        tStates[tree.num] = rand(Categorical(eq))
        for nd in preo[2:end]
            mState = tStates[nd.mother.num]
            p = real.(U * diagm(exp.(nd.inc_length * D)) * Uinv)
            dState = rand(Categorical(p[mState, :]))
            tStates[nd.num] = dState
        end
        for nd in preo
            if nd.nchild == 0
                simChar[nd.name] = tStates[nd.num]
            end
        end
    end

    x = [simChar[l] for l in taxa]
    x_bin = DataFrame(hcat([x .== i for i = 1:4]...))
    x_bin[!, :glot_fam] = d.glot_fam
    x_bin = x_bin[defined,:]

    tv = mapslices(myvar, Array(x_bin[:, 1:4]), dims = 1) |> sum
    lmv = @pipe x_bin |>
          groupby(_, :glot_fam) |>
          combine(x -> sum(mapslices(myvar, Array(x[:, 1:4]), dims = 1)), _) |>
          mean(_[:, 2])

    vc = @pipe x_bin |>
          groupby(_, :glot_fam) |>
          combine(x -> mapslices(mean, Array(x[:, 1:4]), dims = 1), _) |>
          mapslices(myvar, Array(_[:, 2:end]), dims = 1) |>
          sum

    push!(stats, [tv, lmv, vc])
end

##
uSim = DataFrame(Tuple.(stats))
rename!(uSim, [:totalVariance, :lineagewiseVariance, :crossLineageVariance])

##

uSim[:,:method] .= "simulated"
uSim[:,:model] .= "universal"

##
sim = read("modelFitting/output/lineage_$(lpad(charNum, 2, "0")).jls", ModelChains)

samples = vcat([sim.value[:,:,k] for k in sim.chains]...)

##

stats = []
@showprogress for g in 1:nSims
    i = rand(1:size(samples, 1))
    treeI = ceil(Int,invlogit(samples[i,1])*1000)
    simChar = Dict{String,Int}()
    fullSrates = reshape(exp.(samples[i,2:(end-1)]), 12, :)

    for fmI in 1:length(lineages)
        q = zeros(4,4)
        q[offDiag] = fullSrates[:,fmI]
        for i in 1:4
            q[i,i] -= sum(q[i,:])
        end
        D, U = eigen(q)
        Uinv = inv(U)

        eq = abs.(Uinv[4,:])
        eq /= sum(eq)
        tree = fm2trees[lineages[fmI]][treeI]
        preo = pre_order(tree)
        tStates = Vector{Int}(undef, length(preo))
        tStates[tree.num] = rand(Categorical(eq/sum(eq)))
        for nd in preo[2:end]
              mState = tStates[nd.mother.num]
              p = real.(U * diagm(exp.(nd.inc_length * D)) * Uinv)
              ps = abs.(p[mState,:]/sum(p[mState,:]))
              ps /= sum(ps)
              dState = rand(Categorical(ps))
              tStates[nd.num] = dState
        end
        for nd in preo
            if nd.nchild == 0
                simChar[nd.name] = tStates[nd.num]
            end
        end
    end

    x = [simChar[l] for l in taxa]
    x_bin = DataFrame(hcat([x .== i for i in 1:4]...))
    x_bin[!,:glot_fam] = d.glot_fam
    x_bin = x_bin[defined,:]

    tv = mapslices(myvar, Array(x_bin[:,1:4]), dims=1) |> sum
    lmv = @pipe x_bin |> groupby(_, :glot_fam) |>
        combine(x -> sum(mapslices(myvar, Array(x[:,1:4]), dims=1)), _) |>
        mean(_[:,2])

    vc = @pipe x_bin |> groupby(_, :glot_fam) |>
        combine(x -> mapslices(mean, Array(x[:,1:4]), dims=1), _) |>
    mapslices(myvar, Array(_[:,2:end]), dims=1) |> sum

    push!(stats, [tv, lmv, vc])

end
##
lSim = DataFrame(Tuple.(stats))
rename!(lSim, names(uSim)[1:3])

lSim[:,:method] .= "simulated"
lSim[:,:model] .= "lineageSpecific"

##

try
    mkdir("ppp")
catch e
end

CSV.write("ppp/$(lpad(charNum,2,"0")).csv", vcat(uSim, lSim))

max_stars_repo_name: erathorn/phylogeneticTypology

max_issues_repo_name: erathorn/phylogeneticTypology

max_forks_repo_name: erathorn/phylogeneticTypology


text:
using OrthogonalPolynomialsQuasi, ContinuumArrays, DifferentialEquations, Plots

###
# Heat
###

S = JacobiWeight(1.0,1.0) .* Jacobi(1.0,1.0)
D = Derivative(axes(S,1))
Δ = -((D*S)'*(D*S))
M = S'S

function evolution(u, (M,A), t)
    M\(A*u)
end
n = 50
u0 = [[1,2,3]; zeros(n-3)]
prob = ODEProblem(evolution,u0,(0.0,3.0),(cholesky(Symmetric(M[1:n,1:n])),Δ[1:n,1:n]))
@time u = solve(prob,TRBDF2()); u(1.0)

###
# Heat natural
###
L = LinearSpline(range(-1,1;length=2))
P = apply(hcat,L,S)
ApplyArray(hcat,Legendre() \ (D*L), Legendre() \ (D*S))

(D*L)'*(D*L)

(D*S)'*(D*S)

((D*S)'*(D*L)).args[3]
(D*L)'*(D*S)
H = (D*L).args[1]
P = Legendre()
(P\H)' * P'P * (P \ (D*S))

###
# Beam
###

# [u, u_t]_t = [0 I; M] * [u, u_t]

S = JacobiWeight(2.0,2.0) .* Jacobi(2.0,2.0)
D = Derivative(axes(S,1))
Δ² = (D^2*S)'*(D^2*S)
M = S'S
function wave(uv, (M,A), t)
    n = size(M,1)
    u,v = uv[1:n],uv[n+1:end]
    [v; M\(A*u)]
end
n = 30
u0 = [1; 2; 3; zeros(n-3)]
prob = ODEProblem(wave,[u0; zeros(n)],(0.0,3.0),(cholesky(Symmetric(M[1:n,1:n])),-(Δ²)[1:n,1:n]))
@time u = solve(prob,TRBDF2()); u(1.0)

g = range(-1,1;length=200)
V = S[g,1:n]
@gif for t in 0.0:0.01:3.0
    plot(g, V*u(t)[1:n]; ylims=(-5,5))
end

prob = ODEProblem(wave,[u0; zeros(n)],(0.0,3.0),(qr(M),-(Δ²)))


text:
@testset "Tensor-CUR" begin
Random.seed!(12345)

@testset "Small case" begin
    T = rand_kruskal3(2, (20, 40, 60), false)
    @testset "slab axis: $i" for i in 1:3
        @time cur3 = tensorcur3(T, 3, 15, i)
        @test sum(cur3.Cweight) == 3
        @test sum(cur3.Rweight) == 15
        @test all(err -> err < 1e-5, cur3.error)
    end
end

@testset "Large case without reconstruction" begin
    T = rand_kruskal3(2, (100, 120, 60), false)
    @time factors = tensorcur3(T, 10, 200, compute_u=false)
    @test sum(factors.Cweight) == 10
    @test sum(factors.Rweight) == 200
    @test size(factors.U) == (0, 0)
    @test size(factors.error) == (0,)
end

end


text:
export TransposeOperator

type TransposeOperator{T<:Number,B<:Operator} <: Operator{T}
    op::B
end

TransposeOperator{T<:Number}(B::Operator{T})=TransposeOperator{T,typeof(B)}(B)

Base.convert{T}(::Type{Operator{T}},A::TransposeOperator)=TransposeOperator(convert(Operator{T},A.op))

domainspace(P::TransposeOperator)=rangespace(P.op)
rangespace(P::TransposeOperator)=domainspace(P.op)

domain(P::TransposeOperator)=domain(P.op)

bandinds(P::TransposeOperator)=-bandinds(P.op)[end],-bandinds(P.op)[1]

getindex(P::TransposeOperator,k::Integer,j::Integer) =
    P.op[j,k]

function Base.convert{T,TO<:TransposeOperator}(::Type{BandedMatrix},S::SubOperator{T,TO})
    kr,jr=parentindexes(S)
    BandedMatrix(view(parent(S).op,jr,kr)).'
end

Base.transpose(A::Operator)=TransposeOperator(A)
Base.ctranspose{T<:Real}(A::Operator{T})=TransposeOperator(A)


text:
export Circle
export Square
export Rectangle
export AbstractDisk
export AbstractQuadrilateral
export AbstractRectangle
export AbstractUniformDisk

using Formatting
using LinearAlgebra
using NLopt
using RecipesBase
using StaticArrays

abstract type AbstractDisk <: AbstractShape{2,1} end
abstract type AbstractUniformDisk <: AbstractDisk end
abstract type AbstractQuadrilateral <: AbstractShape{2,4} end
abstract type AbstractRectangle <: AbstractQuadrilateral end

################################################################################
# CIRCLE
"""
    Circle <: AbstractDisk

Fields: `radius`, `origin`, `ϕ`, `sinϕ`, `cosϕ`, `corners`, `models`
"""
struct Circle <: AbstractUniformDisk
    radius::Float64
    origin::Point{2,Cartesian}
    ϕ::Float64
    sinϕ::Float64
    cosϕ::Float64
    corners::Tuple{}
    models::Array{Opt,1}

    function Circle(radius::Real, origin::Point{2}, ϕ::Real)
        snϕ, csϕ = sincos(ϕ)
        c = new(radius,origin,ϕ,snϕ,csϕ)
        models = build_models(c)
        return new(radius,origin,ϕ,snϕ,csϕ,(),models)
    end
end

"""
    Circle(radius, origin::Point{2}; ϕ=0, reference=:center) -> circle
"""
function Circle(radius::Real, origin::Point{2}; ϕ::Real=0, reference::Symbol=:center)
    if Base.sym_in(reference, (:bottom, :Bottom, :b, :B))
        origin = origin + rotate( Point( 0, radius), cos(ϕ), sin(ϕ))
    elseif Base.sym_in(reference,(:top, :Top, :T, :t))
        origin = origin + rotate( Point( 0,-radius), cos(ϕ), sin(ϕ))
    elseif Base.sym_in(reference,(:left, :Left, :L, :l))
        origin = origin + rotate( Point( radius, 0), cos(ϕ), sin(ϕ))
    elseif Base.sym_in(reference,(:right, :Right, :R, :r))
        origin = origin + rotate( Point(-radius, 0), cos(ϕ), sin(ϕ))
    end
    return Circle(radius,origin,ϕ)
end

"""
    Circle(radius, x0, y0; ϕ=0, reference=:center) -> circle
"""
Circle(radius::Real, x0::Real, y0::Real; kwargs...) = Circle(radius, Point(x0,y0); kwargs...)

"""
    (::Circle)(p::Point{2}) --> is_in_circle::Bool
"""
(c::Circle)(p::Point{2}) = norm(p-c.origin) < c.radius-eps(c.radius)

"""
    perimeter(t,shape,side) -> x,y,dx/dt,dy/dt
"""
function perimeter(t::Number,c::Circle,side::Int)
    sn,cs = sincos(t)
    x = c.origin.x + c.radius*cs
    y = c.origin.y + c.radius*sn
    ẋ = -c.radius*sn
    ẏ =  c.radius*cs
    return x,y,ẋ,ẏ
end

function Base.getproperty(c::Circle, sym::Symbol)
    if Base.sym_in(sym,(:R,:r,:rad))
        return getfield(c,:radius)
    else
        return getfield(c,sym)
    end
end

function Base.propertynames(::Circle,private=false)
    if private
        return fieldnames(Circle)
    else
        return (:radius, :origin)
    end
end

import ..PRINTED_COLOR_NUMBER
import ..PRINTED_COLOR_DARK

function Base.show(io::IO, circle::Circle)
    printstyled(io, "Circle", color=PRINTED_COLOR_DARK)
    print(io, "(radius: ")
    printstyled(io, fmt("2.2f",circle.radius), color=PRINTED_COLOR_NUMBER)
    print(io, ", center: ", circle.origin, ")")
end

################################################################################
# SQUARE
"""
    Square <: AbstractQuadrilateral <: AbstractShape{2,4}

Fields: `a`, `origin`, `ϕ`, `sinϕ`, `cosϕ`, `corners`, `models`
"""
struct Square <: AbstractRectangle
    a::Float64
    origin::Point{2,Cartesian}
    ϕ::Float64
    sinϕ::Float64
    cosϕ::Float64
    corners::NTuple{4,Point{2,Cartesian}}
    models::Array{Opt,1}

    function Square(a::Real, origin::Point{2}, ϕ::Number)
        s = new(a,origin,ϕ,sin(ϕ),cos(ϕ))
        models = build_models(s)
        c1 = rotate(s.origin+Point( s.a/2, s.a/2), s)
        c2 = rotate(s.origin+Point(-s.a/2, s.a/2), s)
        c3 = rotate(s.origin+Point(-s.a/2,-s.a/2), s)
        c4 = rotate(s.origin+Point( s.a/2,-s.a/2), s)
        return new(a,origin,ϕ,sin(ϕ),cos(ϕ),(c1,c2,c3,c4),models)
    end
end

"""
    Square(a, origin::Point; ϕ=0, reference=:center) -> square

`a` length of sides
"""
function Square(a::Real, origin::Point{2}; ϕ::Real=0, reference::Symbol=:center)
    if Base.sym_in(reference,(:bottom, :Bottom, :b, :B))
        origin += rotate(Point( 0, a/2),cos(ϕ),sin(ϕ))
    elseif Base.sym_in(reference,(:top, :Top, :T, :t))
        origin += rotate(Point( 0,-a/2),cos(ϕ),sin(ϕ))
    elseif Base.sym_in(reference,(:left, :Left, :L, :l))
        origin += rotate(Point( a/2, 0),cos(ϕ),sin(ϕ))
    elseif Base.sym_in(reference,(:right, :Right, :R, :r))
        origin += rotate(Point(-a/2, 0),cos(ϕ),sin(ϕ))
    elseif Base.sym_in(reference,(:topright, :TopRight, :TR, :tr, :ne, :NE))
        origin += rotate(Point(-a/2,-a/2),cos(ϕ),sin(ϕ))
    elseif Base.sym_in(reference,(:topleft, :TopLeft, :TL, :tl, :nw, :NW))
        origin += rotate(Point( a/2,-a/2),cos(ϕ),sin(ϕ))
    elseif Base.sym_in(reference,(:bottomleft, :BottomLeft, :BL, :bl, :sw, :SW))
        origin += rotate(Point( a/2, a/2),cos(ϕ),sin(ϕ))
    elseif Base.sym_in(reference,(:bottomright, :BottomRight, :BR, :br, :se, :SE))
        origin += rotate(Point(-a/2, a/2),cos(ϕ),sin(ϕ))
    end
    return Square(a,origin,ϕ)
end

"""
    Square(a,x0,y0; ϕ=0, reference=:center)
"""
Square(a::Real, x0::Real, y0::Real; kwargs...) = Square(a,Point(x0,y0); kwargs...)

"""
    (::Square)(p::Point) -> is_in_square::Bool
"""
function (s::Square)(p::Point{2})
    purot = unrotate(p, s)
    p = purot - s.origin
    return -s.a/2+eps(s.a/2) < p.x < s.a/2-eps(s.a/2)  &&  -s.a/2+eps(s.a/2) < p.y < s.a/2-eps(s.a/2)
end

function perimeter(t::Number,s::Square,side::Int)
    if side==1
        x = -s.a/2
        y = -s.a/2 + t*s.a
        ẋ = 0
        ẏ = s.a
    elseif side==2
        x =  s.a/2
        y = -s.a/2 + t*s.a
        ẋ = 0
        ẏ = s.a
    elseif side==3
        x = -s.a/2 + t*s.a
        y = -s.a/2
        ẋ = s.a
        ẏ = 0
    elseif side==4
        x = -s.a/2 + t*s.a
        y =  s.a/2
        ẋ = s.a
        ẏ = 0
    end
    x += s.origin.x
    y += s.origin.y
    p = rotate(Point(x,y),s)
    ṗ = rotate(Point(ẋ,ẏ),s.cosϕ,s.sinϕ)
    return p.x, p.y, ṗ.x, ṗ.y
end

function Base.getproperty(s::Square, sym::Symbol)
    if Base.sym_in(sym,(:A,:B,:b))
        return getfield(s,:a)
    else
        return getfield(s,sym)
    end
end

function Base.propertynames(::Square,private=false)
    if private
        return fieldnames(Square)
    else
        return (:a, :origin, :ϕ)
    end
end

import ..PRINTED_COLOR_NUMBER
import ..PRINTED_COLOR_DARK

function Base.show(io::IO, square::Square)
    printstyled(io, "Square", color=PRINTED_COLOR_DARK)
    print(io, "(a: ")
    printstyled(io, fmt("2.2f",square.a), color=PRINTED_COLOR_NUMBER)
    print(io, ", ϕ: ")
    printstyled(io, fmt("3.1f",180square.ϕ/π),"°", color=PRINTED_COLOR_NUMBER)
    print(io, ", center: ", square.origin, ")")
end

################################################################################
# RECTANGLE
"""
    Rectangle <: AbstractQuadrilateral <: AbstractShape{2,4}

Fields: `a`, `b`, `origin`, `ϕ`, `sinϕ`, `cosϕ`, `corners`, `models`
"""
struct Rectangle <: AbstractRectangle
    a::Float64
    b::Float64
    origin::Point{2,Cartesian}
    ϕ::Float64
    sinϕ::Float64
    cosϕ::Float64
    corners::NTuple{4,Point{2,Cartesian}}
    models::Array{Opt,1}

    function Rectangle(a::Real, b::Real, origin::Point{2}, ϕ::Real)
        r = new(a,b,origin,ϕ,sin(ϕ),cos(ϕ))
        models = build_models(r)
        c1 = rotate(r.origin+Point( r.a/2, r.b/2), r)
        c2 = rotate(r.origin+Point(-r.a/2, r.b/2), r)
        c3 = rotate(r.origin+Point(-r.a/2,-r.b/2), r)
        c4 = rotate(r.origin+Point( r.a/2,-r.b/2), r)
        return new(a,b,origin,ϕ,sin(ϕ),cos(ϕ),(c1,c2,c3,c4),models)
    end
end

"""
    Rectangle(a, b, origin::Point; ϕ=0, reference=:center) -> r

`a`, `b` length of sides, `ϕ` is angle of rotation relative to `reference`
"""
function Rectangle(a::Real, b::Real, origin::Point{2}; ϕ::Real=0, reference::Symbol=:center)
    if Base.sym_in(reference,(:bottom, :Bottom, :b, :B))
        origin += rotate(Point( 0, b/2), cos(ϕ), sin(ϕ))
    elseif Base.sym_in(reference,(:top, :Top, :T, :t))
        origin += rotate(Point( 0,-b/2), cos(ϕ), sin(ϕ))
    elseif Base.sym_in(reference,(:left, :Left, :L, :l))
        origin += rotate(Point( a/2, 0), cos(ϕ), sin(ϕ))
    elseif Base.sym_in(reference,(:right, :Right, :R, :r))
        origin += rotate(Point(-a/2, 0), cos(ϕ), sin(ϕ))
    elseif Base.sym_in(reference,(:topright, :TopRight, :TR, :tr, :ne, :NE))
        origin += rotate(Point(-a/2,-b/2), cos(ϕ), sin(ϕ))
    elseif Base.sym_in(reference,(:topleft, :TopLeft, :TL, :tl, :nw, :NW))
        origin += rotate(Point( a/2,-b/2), cos(ϕ), sin(ϕ))
    elseif Base.sym_in(reference,(:bottomleft, :BottomLeft, :BL, :bl, :sw, :SW))
        origin += rotate(Point( a/2, b/2), cos(ϕ), sin(ϕ))
    elseif Base.sym_in(reference,(:bottomright, :BottomRight, :BR, :br, :se, :SE))
        origin += rotate(Point(-a/2, b/2), cos(ϕ), sin(ϕ))
    end
    return Rectangle(a,b,origin,ϕ)
end

"""
    Rectangle(a, b, x0, y0; ϕ=0, reference=:center)
"""
Rectangle(a::Real, b::Real, x0::Real, y0::Real; kwargs...) = Rectangle(a,b,Point(x0,y0); kwargs...)

"""
    Rectangle(s::Square) -> r
"""
Rectangle(s::Square) = Rectangle(s.a,s.a,s.origin; ϕ = s.ϕ)

"""
    (::Rectangle)(p::Point) -> is_in_rectangle::Bool
"""
function (r::Rectangle)(p::Point{2})
    purot = unrotate(p, r)
    p = purot - r.origin
    return -r.a/2+eps(r.a/2) < p.x < r.a/2-eps(r.a/2)  &&  -r.b/2+eps(r.b/2) < p.y < r.b/2-eps(r.b/2)
end

function perimeter(t::Number,r::Rectangle,side::Int)
    if side==1
        x = -r.a/2
        y = -r.b/2 + t*r.b
        ẋ = 0
        ẏ = r.b
    elseif side==2
        x =  r.a/2
        y = -r.b/2 + t*r.b
        ẋ = 0
        ẏ = r.b
    elseif side==3
        x = -r.a/2 + t*r.a
        y = -r.b/2
        ẋ = r.a
        ẏ = 0
    elseif side==4
        x = -r.a/2 + t*r.a
        y =  r.b/2
        ẋ = r.a
        ẏ = 0
    end
    x += r.origin.x
    y += r.origin.y
    p = rotate(Point(x,y),r)
    ṗ = rotate(Point(ẋ,ẏ),r.cosϕ,r.sinϕ)
    return p.x, p.y, ṗ.x, ṗ.y
end

function Base.getproperty(r::Rectangle, sym::Symbol)
    if Base.sym_in(sym,(:A,))
        return getfield(r,:a)
    elseif Base.sym_in(sym,(:B,))
        return getfield(r,:b)
    else
        return getfield(r,sym)
    end
end

function Base.propertynames(::Rectangle,private=false)
    if private
        return fieldnames(Rectangle)
    else
        return (:a, :b, :origin, :ϕ)
    end
end

import ..PRINTED_COLOR_NUMBER
import ..PRINTED_COLOR_DARK

function Base.show(io::IO, rect::Rectangle)
    printstyled(io, "Rectangle", color=PRINTED_COLOR_DARK)
    print(io, "(a: ")
    printstyled(io, fmt("2.2f", rect.a), color=PRINTED_COLOR_NUMBER)
    print(io, ", b: ")
    printstyled(io, fmt("2.2f", rect.b), color=PRINTED_COLOR_NUMBER)
    print(io, ", ϕ: ")
    printstyled(io, fmt("3.1f", 180rect.ϕ/π),"°", color=PRINTED_COLOR_NUMBER)
    print(io, ", center: ", rect.origin, ")")
end

################################################################################
# AUXILLIARIES

import ..NL_NORMAL_ALGORITHM
import ..NL_NORMAL_XTOL
import ..NL_NORMAL_FTOL
import ..NL_NORMAL_MAXEVAL

"""
    rotate(x,y,cosϕ,sinϕ) -> xrot, yrot

rotate by `ϕ` about origin

----------------
    rotate(x,y,shape) -> xrot, yrot

rotate about `shape`'s center by `shape`'s angle
"""
rotate(p::Point{2},cosϕ::Real,sinϕ::Real) = iszero(sinϕ) ? p : Point(cosϕ*p.x-sinϕ*p.y, sinϕ*p.x+cosϕ*p.y)
rotate(p::Point{2},s::AbstractShape{2}) = s.origin + rotate(p-s.origin,s.cosϕ,s.sinϕ)

function rotate(p::AbstractArray,cosϕ::Real,sinϕ::Real)
    prot = Array{Point{2}}(undef,size(p))
    rotate!(prot,p,cosϕ,sinϕ)
    return prot
end
function rotate(p::AbstractArray,s::AbstractShape{2})
    prot = Array{Point{2}}(undef,size(p))
    rotate!(prot,p.-s.origin,s.cosϕ,s.sinϕ)
    return prot .+ s.origin
end

function rotate!(prot::AbstractArray,p::Point{2},cosϕ::Real,sinϕ::Real)
    for i ∈ eachindex(p) prot[i] = rotate(p[i],cosϕ,sinϕ) end
    return nothing
end

"""
    unrotate(x,y,cosϕ,sinϕ) -> xrot, yrot

rotate by `-ϕ` about origin

----------------
    unrotate(x,y,shape) -> xrot, yrot

rotate about `shape`'s center by the negative of `shape`'s angle
"""
unrotate(p::Point{2},cosϕ::Real,sinϕ::Real) = rotate(p,cosϕ,-sinϕ)
unrotate(p::Point{2},s::AbstractShape{2}) = s.origin + unrotate(p-s.origin,s.cosϕ,s.sinϕ)
unrotate(p::AbstractArray,s::AbstractShape{2}) = s.origin .+ unrotate(p.-s.origin,s.cosϕ,s.sinϕ)

################################################################################
# UNSYMMETRIC AUXILLIARIES

function build_models(s::AbstractShape{N}) where N
    models = Array{Opt}(undef,N)
    for i ∈ eachindex(models)
        models[i] = Opt(NL_NORMAL_ALGORITHM, 1)
        models[i].xtol_rel = NL_NORMAL_XTOL
        models[i].ftol_rel = NL_NORMAL_FTOL
        models[i].maxeval = NL_NORMAL_MAXEVAL
        if typeof(s)<:AbstractQuadrilateral
            models[i].lower_bounds = 0
            models[i].upper_bounds = 1
        end
    end
    return models
end

function model_objective(v::Vector,grad::Vector,X::Number,Y::Number,s::AbstractShape,side::Int)
    x,y,ẋ,ẏ = perimeter(v[1],s,side)
    if length(grad) > 0
        grad[1] = 2(x-X)*ẋ + 2(y-Y)*ẏ
    end
    return (x-X)^2 + (y-Y)^2
end

"""
    normal_distance(shape,x,y) -> normal, tangent, distance, side

Get the `normal` vector from the surface of `shape` to a point `x,y`,
and the associated `tangent` vector and distance from surface `side`.

The returned arguments are sorted in ascending order by `distance`, so that, e.g.,
`normal[1]` is the normal vector associated with the closest side.
"""
function normal_distance(s::TS,x::Number,y::Number) where TS<:AbstractShape{2}
    n = Vector{Vector{Float64}}(undef,length(s.models))
    d = Vector{Float64}(undef,length(s.models))
    sides = 1:length(s.models)
    for i ∈ eachindex(s.models)
        model = s.models[i]
        model.min_objective = (v,g)->model_objective(v,g,x,y,s,i)
        minf,mint,ret = optimize(model,[t0_shape(x,y,s)])
        X,Y = perimeter(mint[1],s,i)
        n[i] = [x-X, y-Y]
        n[i] = n[i]/norm(n[i])
        d[i] = sqrt(minf)
    end
    perm = sortperm(d)
    n, d = n[perm], d[perm]
    NORMAL = Vector{Vector{Float64}}(undef,length(s.models))
    TANGENT = Vector{Vector{Float64}}(undef,length(s.models))
    for i ∈ eachindex(NORMAL)
        nx, ny = n[i][1],n[i][2]
        NORMAL[i] = [nx, ny]
        TANGENT[i] = [ny,-nx]
    end
    return NORMAL, TANGENT, d, sides[perm]
end

t0_shape(x::Number,y::Number,s::AbstractQuadrilateral) = .5
t0_shape(x::Number,y::Number,s::AbstractDisk) = atan(y-s.origin.y,x-s.origin.x)

################################################################################
# Plotting

import ..SHAPE_COLOR
import ..SHAPE_FILL_ALPHA

# plotting for all disks
@recipe function f(d::AbstractDisk)
    ϕ = LinRange(0,2π,201)
    x = Vector{Float64}(undef,length(ϕ))
    y = Vector{Float64}(undef,length(ϕ))
    for i ∈ eachindex(ϕ) x[i],y[i] = perimeter(ϕ[i],d,1) end
    alpha --> 0
    seriestype --> :path
    fillcolor --> SHAPE_COLOR
    fillrange --> d.origin.y
    fillalpha --> SHAPE_FILL_ALPHA
    aspect_ratio --> 1
    legend-->false
    (x,y)
end

# plot all quadrilaterals
@recipe function f(q::AbstractQuadrilateral)
    t = LinRange(0,1,201)
    x1,y1 = Array{Float64}(undef,length(t)), Array{Float64}(undef,length(t))
    x2,y2 = Array{Float64}(undef,length(t)), Array{Float64}(undef,length(t))
    x3,y3 = Array{Float64}(undef,length(t)), Array{Float64}(undef,length(t))
    x4,y4 = Array{Float64}(undef,length(t)), Array{Float64}(undef,length(t))
    for i ∈ eachindex(t)
        x1[i],y1[i] = perimeter(t[i],q,1)
        x2[i],y2[i] = perimeter(t[i],q,2)
        x3[i],y3[i] = perimeter(t[i],q,3)
        x4[i],y4[i] = perimeter(t[i],q,4)
    end
    alpha --> 0
    seriestype --> :path
    fillcolor --> SHAPE_COLOR
    fillrange --> q.origin.y
    fillalpha --> SHAPE_FILL_ALPHA
    aspect_ratio --> 1
    legend --> false
    vcat(x1,x4,reverse(x3)),vcat(y1,y4,reverse(y3))
end

#
# """
#     struct Ellipse <: AbstractDisk <: AbstractShape{1}
#
#     Ellipse((a,b), x0, y0, ϕ=0; reference=:center) -> ellipse
#     Ellipse(a, b, x0, y0, ϕ=0; reference=:center) -> ellipse
#
#
# `a, b` axis lengths
# `x0, y0` is location of `reference`
# `ϕ` angle of `a` axis
#
# ----------------
#
#     (::Ellipse)(x,y) -> is_in_ellipse::Bool
# """
# struct Ellipse <: AbstractDisk
#     a::Float64
#     b::Float64
#     x0::Float64
#     y0::Float64
#     ϕ::Float64
#     sinϕ::Float64
#     cosϕ::Float64
#     corners::Tuple{}
#     models::Array{Opt,1}
#
#     Ellipse((a,b),args...;kwargs...) = Ellipse(a,b,args...;kwargs...)
#     function Ellipse(a::Number,b::Number,x0::Number,y0::Number,ϕ::Number=0; reference::Symbol=:center)
#         if reference ∈ [:bottom, :Bottom, :b, :B]
#             x0,y0 = (x0,y0) .+ rotate( 0, b,cos(ϕ),sin(ϕ))
#         elseif reference ∈ [:top, :Top, :T, :t]
#             x0,y0 = (x0,y0) .+ rotate( 0,-b,cos(ϕ),sin(ϕ))
#         elseif reference ∈ [:left, :Left, :L, :l]
#             x0,y0 = (x0,y0) .+ rotate( a, 0,cos(ϕ),sin(ϕ))
#         elseif reference ∈ [:right, :Right, :R, :r]
#             x0,y0 = (x0,y0) .+ rotate(-a, 0,cos(ϕ),sin(ϕ))
#         end
#          e = new(a,b,x0,y0,ϕ,sin(ϕ),cos(ϕ))
#          models = build_models(e)
#          return new(a,b,x0,y0,ϕ,sin(ϕ),cos(ϕ),(),models)
#     end
#
#     Ellipse(c::Circle) = Ellipse(c.R,c.R,c.x0,c.y0,c.ϕ)
#
#     function (e::Ellipse)(x,y)
#         xurot, yurot = unrotate(x, y, e)
#         return hypot((xurot-e.x0)/e.a, (yurot-e.y0)/e.b) < 1-eps()
#     end
#
#     function Base.show(io::IO, ellipse::Ellipse)
#         print(io, "Ellipse(a=", fmt("2.2f",ellipse.a), ", b=", fmt("2.2f",ellipse.b), ", x0=", fmt("2.2f",ellipse.x0), ", y0=", fmt("2.2f",ellipse.y0), ", ϕ=∠", fmt("3.2f",(mod2pi(ellipse.ϕ))*180/π), "°)")
#     end
# end
# function perimeter(t::Number,e::Ellipse,side::Int)
#     sn,cs = sincos(t)
#     xur,yur =    e.a*cs, e.b*sn
#     xpur,ypur = -e.a*sn, e.b*cs
#     x,y = (e.x0,e.y0) .+ rotate(xur,yur,e.cosϕ,e.sinϕ)
#     ẋ,ẏ = rotate(xpur,ypur,e.cosϕ,e.sinϕ)
#     return x,y,ẋ,ẏ
# end
#
#
# """
#     struct DeformedDisk{N} <: AbstractDisk <: AbstractShape{1}
#
#     DeformedDisk{N}(R, x0, y0, M, a, φ, ϕ=0) -> deformeddisk
#     DeformedDisk{N}((R,M,a,φ), x0, y0, ϕ=0) -> deformeddisk
#
# `R` is radius,
# `x0` and `y0` is center of circle of radius `R`,
# `M` is array of length `N` of multipole integers,
# `a` is array of length `N` amplitudes,
# `φ` is array of length `N` of angles
# `ϕ` is overall rotation angle
#
# -----------------
#
#     (::DeformedDisk)(x,y) -> is_in_disk::Bool
# """
# struct DeformedDisk{N} <: AbstractDisk
#     R::Float64
#     x0::Float64
#     y0::Float64
#     ϕ::Float64
#     M::SArray{Tuple{N},Int,1,N}
#     a::SArray{Tuple{N},Float64,1,N}
#     φ::SArray{Tuple{N},Float64,1,N}
#     sinϕ::Float64
#     cosϕ::Float64
#     corners::Tuple{}
#     models::Array{Opt,1}
#
#     DeformedDisk{n}((R,M,a,φ),x0,y0,args...;kwargs...) where n= DeformedDisk{n}(R,x0,y0,M,a,φ,args...;kwargs...)
#     function DeformedDisk{n}(R::Number,x0::Number,y0::Number,M,a,φ,ϕ=0) where n
#         @assert n==length(M)==length(a)==length(φ) "parameter N in DeformedDisk{N}(...) must be equal to length(M)"
#         d = new{n}(R,x0,y0,ϕ,SVector{n}(M),SVector{n}(a),SVector{n}(φ.+ϕ),sin(ϕ),cos(ϕ))
#         models = build_models(d)
#         return new{n}(R,x0,y0,d.ϕ,d.M,d.a,d.φ,d.sinϕ,d.cosϕ,(),models)
#     end
#
#     DeformedDisk(c::Circle) = DeformedDisk{0}(c.R,c.x0,c.y0,[],[],[],c.ϕ)
#
#     function (d::DeformedDisk)(x,y)
#         ϕ = atan(y-d.y0,x-d.x0)
#         r = hypot(x-d.x0,y-d.y0) < d.R + sum(d.a.*(map((m,φ)->cos(m*(ϕ-φ)),d.M,d.φ))) - eps(d.R)
#         return r
#     end
#
#     function Base.show(io::IO, d::DeformedDisk{N}) where N
#          print(io,"DeformedDisk with $N multipoles")
#     end
# end
# function perimeter(t::Number,d::DeformedDisk,side::Int)
#     r = d.R + sum(d.a.*(map((m,φ)->cos(m*(t-φ)),d.M,d.φ)))
#     ṙ = sum(d.a.*(map((m,φ)->-m*sin(m*(t-φ)),d.M,d.φ)))
#     sn,cs = sincos(t)
#     x,y = d.x0 + r*cs, d.y0 + r*sn
#     ẋ = ṙ*cs - r*sn
#     ẏ = ṙ*sn + r*cs
#     return x,y,ẋ,ẏ
# end

# """
#     struct Annulus <: AbstractShape{2}
#
#     Annulus((R1,R2),x0,y0,ϕ=0) -> annulus
#     Annulus(R1, R2, x0, y0, ϕ=0) -> annulus
#
# `R1` is inner radius, `R2` outer
# `x0, y0` is center
#
# -------------
#     (::Annulus)(x,y) -> is_in_annulus::Bool
# """
# struct Annulus <: AbstractShape{2}
#     R1::Float64
#     R2::Float64
#     x0::Float64
#     y0::Float64
#     ϕ::Float64
#     sinϕ::Float64
#     cosϕ::Float64
#     corners::Tuple{}
#     models::Array{Opt,1}
#
#     Annulus((R1,R2),args...;kwargs...) = Annulus(R1,R2,args...;kwargs...)
#     function Annulus(R1::Number,R2::Number,x0::Number,y0::Number,ϕ::Number=0)
#         @assert R1≤R2 "R1=$R1, R2=$R2 must satisfy R1≤R2"
#         a = new(R1,R2,x0,y0,0,0,1)
#         models = build_models(a)
#         return new(R1,R2,x0,y0,0,sin(ϕ),cos(ϕ),(),models)
#     end
#
#     (c::Annulus)(x,y) = c.R1+eps(c.R1) < hypot( (x-c.x0),(y-c.y0) ) < c.R2-eps(c.R2)
#
#     function Base.show(io::IO, a::Annulus)
#         print(io, "Annulus(R1=$(fmt("2.2f",a.R1)), R2=$(fmt("2.2f",a.R2)), x0=$(fmt("2.2f",a.x0)), y0=$(fmt("2.2f",a.y0)))")
#     end
# end
# function perimeter(t::Number,a::Annulus,side::Int)
#     sn,cs = sincos(t)
#     R = side==1 ? a.R1 : a.R2
#     x = a.x0 + R*cs
#     y = a.y0 + R*sn
#     ẋ = -R*sn
#     ẏ =  R*cs
#     return x,y,ẋ,ẏ
# end
# @recipe function f(cr::Annulus)
#     @series Circle(cr.R1,cr.x0,cr.y0)
#     @series Circle(cr.R2,cr.x0,cr.y0)
# end
#
#
#

#
#
# # """
# #     Parallelogram(a, b, α, x0, y0, ϕ)
# # """
# # struct Parallelogram <: AbstractParallelogram
# #     a::Float64
# #     b::Float64
# #     α::Float64
# #     x0::Float64
# #     y0::Float64
# #     ϕ::Float64
# #     sinϕ::Float64
# #     cosϕ::Float64
# #     tanα::Float64
# #     cosα::Float64
# #     models::Array{Opt,1}
# #
# #     function Parallelogram(a::Number, b::Number, α::Number, x0::Number, y0::Number, ϕ::Number)
# #         new(a,b,α,x0,y0,ϕ,cos(ϕ),sin(ϕ),tan(α),cos(α))
# #     end
# #
# #     function (p::Parallelogram)(x,y)
# #         xrot, yrot = rotate(x-p.x0, y-p.y0, p.cosϕ, p.sinϕ)
# #         return p.tanα*(xrot-p.a) < yrot < p.tanα*xrot  &&  0 < yrot < p.cosα*p.b
# #     end
# #
# #     Base.show(io::IO, par::Parallelogram) = print(io, "Parallelogram(a=$(fmt("2.2f",par.a)), b=$(fmt("2.2f",par.b)), α=$(fmt("2.2f",par.α)), x0=$(fmt("2.2f",par.x0)), y0=$(fmt("2.2f",par.y0)), ϕ=$(fmt("2.2f",par.ϕ)))")
# #
# #     @recipe function f(pg::Parallelogram)
# #         x = cumsum([0, pg.a, +pg.b*cos(pg.α), -pg.a, -pg.b*cos(pg.α)])
# #         y = cumsum([0, 0, +pg.b*sin(pg.α), 0, -pg.b*sin(pg.α)])
# #         alpha --> 0
# #         seriestype --> :path
# #         fillcolor --> SHAPE_COLOR
# #         fillrange --> 0
# #         fillalpha --> SHAPE_FILL_ALPHA
# #         aspect_ratio --> 1
# #         legend --> false
# #         xrot, yrot = rotate(x,y,pg.cosϕ,-pg.sinϕ)
# #         pg.x0 .+ xrot, pg.y0 .+ yrot
# #     end
# # end
#


text:
using MIRT
using MATLAB
using Plots

# setup matlab to use MIRT
if !@isdefined(irtdir)
	ENV["MATLAB_ROOT"] = "/Applications/freeware/matlab"

	irtdir = "/Users/fessler/src/matlab/alg"
	tmp = "addpath('$irtdir')"
	eval_string(tmp)
	mat"setup"
end

fov = 500.
ig = image_geom( ; nx=512, ny=504, fov=fov)
down = 2 # small size for debug
ig = ig.down(down)

mat"ig = image_geom('nx', $(ig.nx), 'ny', $(ig.ny), 'fov', $fov, 'dy', 'dx');"

sg = sino_geom(:par, nb=888, na=984, offset=0.25, d=541/959)
sg = sg.down(down)

mat"sg = sino_geom('par', 'nb', $(sg.nb), 'na', $(sg.na), 'dr', $(sg.d), 'strip_width', 'd');"

ell = [30. 0 120 100 0 1]
ell = [0. 0 200 100 0 1]
ell = [0. 40 200 200 0 1;
60. 0 120 100 0 1]

im_julia = ellipse_im(ig, ell)

@mput ell
mat"im_mat = ellipse_im(ig, ell)"

@mget im_mat

plot(jim(im_julia, "julia"), jim(im_mat, "matlab"), jim(im_mat - im_julia))

sino_julia = ellipse_sino(sg, ell)

mat"[sino_mat, pos, ang] = ellipse_sino(sg, ell)"
mat"sizeof(sino_mat)"
#@mget sino_mat

#plot(jim(sino_julia, "julia"), jim(sino_mat, "matlab"), jim(sino_mat - sino_julia))


text:
@testset "NALUX" begin
    nalu = NALUX(2, 3, Dense(2, 3))
    batch = 30
    generate() = randn(Float32, 2, batch) .* 4
    f(x) = x[1,:] .* x[2,:]
    x = generate()
    z = nalu(x)

    @test size(z) == (3,batch)

    ps = params(nalu)
    @test length(ps) == 5

    # TODO: fix NALUX tests
    # mse(x) = Flux.mse(nalu(x), reshape(repeat(f(x), inner=(3,1)),3,batch))
    # reg(;β=0.01f0) = β*sum(norm, [nalu.A,nalu.rM,nalu.iM])
    # function bce(x;β=0.001f0)
    #     g = gate(nalu,x)
    #     β*sum(Flux.binarycrossentropy.(g,g))
    # end

    # loss(x) = mse(x) + reg() + bce(x)
    # #loss(x) = mse(x)
    # opt = ADAM(0.01)
    # train_data = [(generate(),) for _ in 1:10000]
    # Flux.train!(loss, ps, train_data, opt, cb=()->(@info loss(x)))

    # display(nalu.rM)
    # display(nalu.iM)
    # display(nalu.A)
    # display(gate(nalu,x))
    # @test all(isapprox.(nalu.rM, [1.0 1.0; 1.0 1.0; 1.0 1.0], atol=1e-1))
    # @test all(isapprox.(gate(nalu, x[:,1]), [0.0, 0.0, 0.0], atol=1e-1))
end


text:
function H5ForwardModeling(;
  model::JUDI.Model,
  opt::JUDI.JUDIOptions,
  q::JUDI.judiVector,
  recGeometry)

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

  # Setup operators
  Pr = judiProjection(recGeometry)
  F = judiModeling(model; options=opt)
  Ps = judiProjection(q.geometry)

  dobs = Pr*F*adjoint(Ps)*q
  @info "Modeling finished\n"
end


text:
shapeconregr_n_ds = [
    [
    (2, 2), # compile run
    (2, 5),
    (2, 10),
    (2, 15),
    (2, 20),
    (2, 25),
    (2, 30),
    ],
    [
    (3, 2), # compile run
    (3, 4),
    (3, 6),
    (3, 8),
    (3, 10),
    (3, 12),
    ],
    [
    (3, 2), # compile run
    (4, 3),
    (4, 4),
    (4, 5),
    (4, 6),
    (4, 7),
    ],
    [
    (3, 2), # compile run
    (6, 2),
    (6, 3),
    (6, 4),
    ],
    [
    (3, 2), # compile run
    (8, 2),
    (8, 3),
    ],
    [
    (3, 2), # compile run
    (10, 2),
    ],
    [
    (3, 2), # compile run
    (12, 2),
    ],
    [
    (3, 2), # compile run
    (14, 2),
    ],
    ]
shapeconregr_insts(use_nat::Bool) = [
    [(n, ceil(Int, 1.1 * binomial(n + 2d, n)), :func4, 100.0, 2d, use_nat,
    false, false, true, false) for (n, d) in nds]
    for nds in shapeconregr_n_ds
    ]

insts = OrderedDict()
insts["nat"] = (nothing, shapeconregr_insts(true))
insts["ext"] = (:SOCExpPSD, shapeconregr_insts(false))
return (ShapeConRegrJuMP, insts)


text:
using MAT

@testset "simulateHBF" begin
    # MATLAB output
    test = matread("resources/simulateHBF.mat")

    x = [1; zeros(2^11-1)]
    for mode = 0:3
        for i = 1:4
            f1, f2 = exampleHBF(i)
            y = simulateHBF(x, f1, f2, mode)
            @test isapprox(test["y_$mode"][:,i], y)
        end
    end
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
module SizeNonInterlacedTest
using SimplePNGs
using Test

include("common.jl")
using .TestCommon: load_json

pl(name) = SimplePNGs.load(joinpath("PngSuite", name*".png"))

@testset "Size test files" begin
    @testset "1x1 paletted file, no interlacing" begin
        img1 = load_json("s01n3p01")
        img2 = pl("s01n3p01")
        @test img1 == img2
    end

    @testset "2x2 paletted file, no interlacing" begin
        img1 = load_json("s02n3p01")
        img2 = pl("s02n3p01")
        @test img1 == img2
    end

    @testset "3x3 paletted file, no interlacing" begin
        img1 = load_json("s03n3p01")
        img2 = pl("s03n3p01")
        @test img1 == img2
    end

    @testset "4x4 paletted file, no interlacing" begin
        img1 = load_json("s04n3p01")
        img2 = pl("s04n3p01")
        @test img1 == img2
    end

    @testset "5x5 paletted file, no interlacing" begin
        img1 = load_json("s05n3p02")
        img2 = pl("s05n3p02")
        @test img1 == img2
    end

    @testset "6x6 paletted file, no interlacing" begin
        img1 = load_json("s06n3p02")
        img2 = pl("s06n3p02")
        @test img1 == img2
    end

    @testset "7x7 paletted file, no interlacing" begin
        img1 = load_json("s07n3p02")
        img2 = pl("s07n3p02")
        @test img1 == img2
    end

    @testset "8x8 paletted file, no interlacing" begin
        img1 = load_json("s08n3p02")
        img2 = pl("s08n3p02")
        @test img1 == img2
    end

    @testset "9x9 paletted file, no interlacing" begin
        img1 = load_json("s09n3p02")
        img2 = pl("s09n3p02")
        @test img1 == img2
    end

    @testset "32x32 paletted file, no interlacing" begin
        img1 = load_json("s32n3p04")
        img2 = pl("s32n3p04")
        @test img1 == img2
    end

    @testset "33x33 paletted file, no interlacing" begin
        img1 = load_json("s33n3p04")
        img2 = pl("s33n3p04")
        @test img1 == img2
    end

    @testset "34x34 paletted file, no interlacing" begin
        img1 = load_json("s34n3p04")
        img2 = pl("s34n3p04")
        @test img1 == img2
    end

    @testset "35x35 paletted file, no interlacing" begin
        img1 = load_json("s35n3p04")
        img2 = pl("s35n3p04")
        @test img1 == img2
    end

    @testset "36x36 paletted file, no interlacing" begin
        img1 = load_json("s36n3p04")
        img2 = pl("s36n3p04")
        @test img1 == img2
    end

    @testset "37x37 paletted file, no interlacing" begin
        img1 = load_json("s37n3p04")
        img2 = pl("s37n3p04")
        @test img1 == img2
    end

    @testset "38x38 paletted file, no interlacing" begin
        img1 = load_json("s38n3p04")
        img2 = pl("s38n3p04")
        @test img1 == img2
    end

    @testset "39x39 paletted file, no interlacing" begin
        img1 = load_json("s39n3p04")
        img2 = pl("s39n3p04")
        @test img1 == img2
    end

    @testset "40x40 paletted file, no interlacing" begin
        img1 = load_json("s40n3p04")
        img2 = pl("s40n3p04")
        @test img1 == img2
    end
end

end # module


text:
function checkfield(field, quotes::Char, escape::Char, trimwhitespace::Bool)
    isquoted = false
    inquotes = false
    escaped = false
    toskip = falses(length(field))
    anyskip = false
    for i in eachindex(field)
        c = field[i]
        if escaped
            escaped = false
        elseif !inquotes
            if trimwhitespace && occursin(r"\s", string(c))
                anyskip = true
                toskip[i] = true
            elseif c == quotes
                inquotes = true
                isquoted = true
                anyskip = true
                toskip[i] = true
            elseif c == escape
                escaped = true
                anyskip = true
                toskip[i] = true
            else
               continue
            end
        else    # inquotes
            if c == quotes
                if quotes != escape || count(c -> c == quotes, field[i:end]) == 1
                    inquotes = false
                    anyskip = true
                    toskip[i] = true
                else
                    escaped = true
                    anyskip = true
                    toskip[i] = true
                end
            elseif c == escape
                escaped = true
                anyskip = true
                toskip[i] = true
            else
               continue
            end
        end
    end
    if (inquotes || escaped)
        return field, isquoted, true
    else
        if anyskip
            f = field[Int[nextind(field, 0, i) for i in findall(!, toskip)]]
            field = SubString(f, firstindex(f), lastindex(f))
        end
        return field, isquoted, false
    end
end

function checkfield(field, quotes::Char, escape::Missing, trimwhitespace::Bool)
    isquoted = false
    badbreak = false
    quoteindices = findall(c -> c == quotes, field)
    if length(quoteindices) == 0
        if trimwhitespace
            field = strip(field)
        end
    elseif length(quoteindices) == 2
        field = SubString(field, nextind(field, last(quoteindices[1])), prevind(field, quoteindices[2]))
        if trimwhitespace
            field = strip(field)
        end
        isquoted = true
    else
        badbreak = true
    end
    return field, isquoted, badbreak
end

function checkfield(field, quotes::Missing, escape::Char, trimwhitespace::Bool)
    isquoted = false
    inquotes = false
    escaped = false
    toskip = falses(length(field))
    anyskip = false
    for i in eachindex(field)
        c = field[i]
        if escaped
            escaped = false
        elseif c == escape
            escaped = true
            anyskip = true
            toskip[i] = true
        else
            continue
        end
    end
    if (inquotes || escaped)
        return field, false, true
    end
    if anyskip
        f = field[Int[nextind(field, 0, i) for i in findall(!, toskip)]]

        field = SubString(f, firstindex(f), lastindex(f))
    end
    if trimwhitespace
        return strip(field), false, false
    else
        return field, false, false
    end
end


text:
using Flux
println("Integro Differential equation tests")
using DiffEqFlux
println("Starting Soon!")
using ModelingToolkit
using DiffEqBase
using Test, NeuralPDE
println("Starting Soon!")
using GalacticOptim
using Optim
using Quadrature,Cubature, Cuba
using QuasiMonteCarlo
using SciMLBase
using OrdinaryDiffEq
import ModelingToolkit: Interval, infimum, supremum
using DomainSets

using Random
Random.seed!(100)

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

#Integration Tests
println("Integral Tests")
@parameters t
@variables i(..)
Di = Differential(t)
Ii = Integral(t in DomainSets.ClosedInterval(0, t))
eq = Di(i(t)) + 2*i(t) + 5*Ii(i(t)) ~ 1
bcs = [i(0.) ~ 0.0]
domains = [t ∈ Interval(0.0,2.0)]
chain = Chain(Dense(1,15,Flux.σ),Dense(15,1))
initθ = Float64.(DiffEqFlux.initial_params(chain))
strategy_ = NeuralPDE.GridTraining(0.1)
discretization = NeuralPDE.PhysicsInformedNN(chain,
                                             strategy_;
                                             init_params = nothing,
                                             phi = nothing,
                                             derivative = nothing,
                                             )
@named pde_system = PDESystem(eq,bcs,domains,[t],[i(t)])
sym_prob = NeuralPDE.symbolic_discretize(pde_system, discretization)
prob = NeuralPDE.discretize(pde_system,discretization)
res = GalacticOptim.solve(prob, BFGS(); cb = cb, maxiters=100)
ts = [infimum(d.domain):0.01:supremum(d.domain) for d in domains][1]
phi = discretization.phi

analytic_sol_func(t) = 1/2*(exp(-t))*(sin(2*t))
u_real  = [analytic_sol_func(t) for t in ts]
u_predict  = [first(phi([t],res.minimizer)) for t in ts]
@test Flux.mse(u_real, u_predict) < 0.001
# plot(ts,u_real)
# plot!(ts,u_predict)

## Simple Integral Test
println("Simple Integral Test")

@parameters x
@variables u(..)
Ix = Integral(x in DomainSets.ClosedInterval(0, x))
eq = Ix(u(x)) ~ (x^3)/3
bcs = [u(0.) ~ 0.0]
domains = [x ∈ Interval(0.0,1.00)]
chain = Chain(Dense(1,15,Flux.σ),Dense(15,1))
initθ = Float64.(DiffEqFlux.initial_params(chain))
strategy_ = NeuralPDE.GridTraining(0.1)
discretization = NeuralPDE.PhysicsInformedNN(chain,
                                             strategy_;
                                             init_params = nothing,
                                             phi = nothing,
                                             derivative = nothing,
                                             )
@named pde_system = PDESystem(eq,bcs,domains,[x],[u(x)])
prob = NeuralPDE.discretize(pde_system,discretization)
res = GalacticOptim.solve(prob, BFGS(); cb = cb, maxiters=100)
xs = [infimum(d.domain):0.01:supremum(d.domain) for d in domains][1]
phi = discretization.phi
u_predict  = [first(phi([x],res.minimizer)) for x in xs]
u_real  = [x^2 for x in xs]
@test Flux.mse(u_real, u_predict) < 0.001

# plot(xs,u_real)
# plot!(xs,u_predict)

#simple multidimensitonal integral test
println("simple multidimensitonal integral test")

@parameters x,y
@variables u(..)
Dx = Differential(x)
Dy = Differential(y)
Ix = Integral((x,y) in DomainSets.UnitSquare())
eq = Ix(u(x,y)) ~ 1/3
bcs = [u(0., 0.) ~ 1, Dx(u(x,y)) ~ -2*x , Dy(u(x ,y)) ~ -2*y ]
domains = [x ∈ Interval(0.0,1.00), y ∈ Interval(0.0,1.00)]
chain = Chain(Dense(2,15,Flux.σ),Dense(15,1))
initθ = Float64.(DiffEqFlux.initial_params(chain))
strategy_ = NeuralPDE.GridTraining(0.1)
discretization = NeuralPDE.PhysicsInformedNN(chain,
                                             strategy_;
                                             init_params = nothing,
                                             phi = nothing,
                                             derivative = nothing,
                                             )
@named pde_system = PDESystem(eq,bcs,domains,[x,y],[u(x, y)])
prob = NeuralPDE.discretize(pde_system,discretization)
res = GalacticOptim.solve(prob, BFGS(); cb = cb, maxiters=100)
xs = 0.00:0.01:1.00
ys = 0.00:0.01:1.00
phi = discretization.phi

u_real = collect(1 - x^2 - y^2 for y in ys, x in xs);
u_predict = collect(Array(phi([x,y], res.minimizer))[1] for y in ys, x in xs);
@test Flux.mse(u_real, u_predict) < 0.001

# error_ = u_predict .- u_real
# p1 = plot(xs,ys,u_real,linetype=:contourf,label = "analytic")
# p2 = plot(xs,ys,u_predict,linetype=:contourf,label = "predict")
# p3 = plot(xs,ys,error_,linetype=:contourf,label = "error")
# plot(p1,p2,p3)

@parameters x,y
@variables u(..)
Dx = Differential(x)
Dy = Differential(y)
Ix = Integral((x,y) in DomainSets.ProductDomain(UnitInterval(),ClosedInterval(0 ,x)))
eq = Ix(u(x,y)) ~ 5/12
bcs = [u(0., 0.) ~ 0, Dy(u(x,y)) ~ 2*y , u(x, 0) ~ x ]
domains = [x ∈ Interval(0.0,1.00), y ∈ Interval(0.0,1.00)]
chain = Chain(Dense(2,15,Flux.σ),Dense(15,1))
initθ = Float64.(DiffEqFlux.initial_params(chain))
strategy_ = NeuralPDE.GridTraining(0.1)
discretization = NeuralPDE.PhysicsInformedNN(chain,
                                             strategy_;
                                             init_params = nothing,
                                             phi = nothing,
                                             derivative = nothing,
                                             )
@named pde_system = PDESystem(eq,bcs,domains,[x,y],[u(x, y)])
prob = NeuralPDE.discretize(pde_system,discretization)
res = GalacticOptim.solve(prob, BFGS(); cb = cb, maxiters=100)
xs = 0.00:0.01:1.00
ys = 0.00:0.01:1.00
phi = discretization.phi

u_real = collect( x + y^2 for y in ys, x in xs);
u_predict = collect(Array(phi([x,y], res.minimizer))[1] for y in ys, x in xs);
@test Flux.mse(u_real, u_predict) < 0.01

# error_ = u_predict .- u_real
# p1 = plot(xs,ys,u_real,linetype=:contourf,label = "analytic")
# p2 = plot(xs,ys,u_predict,linetype=:contourf,label = "predict")
# p3 = plot(xs,ys,error_,linetype=:contourf,label = "error")
# plot(p1,p2,p3)


text:
# include(Pkg.dir("KUnet/test/mnist.jl"))
using MNIST: xtrn, ytrn, xtst, ytst

reload("src/NNN.jl")

const batch_size  = 100
const input_size  = 784
const output_size = 10
const hidden_size = 100

function mnist_train(xtrn, ytrn, xtst, ytst)
    drop = NNN.Dropout(0.5)

    net = NNN.Network(batch_size)
    NNN.push!(net,
             NNN.HiddenLayer(batch_size, input_size,  hidden_size;
                            activation=NNN.relu,
                            post_filters=NNN.Filter[ drop ]),
             NNN.AdaGradRDA(input_size, hidden_size;
                           μ=0.1,
                           λ=0.00001))
    NNN.push!(net,
             NNN.OutputLayer(batch_size, hidden_size, output_size;
                            activation=NNN.softmax,
                            loss=NNN.cross_entropy_cost),
             NNN.SGD(;λ=0.0000001,
                    ε=0.2))

    drop.enabled = false
    NNN.train(xtrn, ytrn, net; iter=3, batch_size=batch_size)
    NNN.test(xtst,  ytst, net; batch_size=batch_size)
end

mnist_train(xtrn, ytrn, xtst, ytst)


text:
# example assuming * operation is always binary

# ENV["JULIA_DEBUG"] = Metatheory

using Metatheory
using Metatheory.Util

include("numberfold.jl")

comm_monoid = @theory begin
    a * b => b * a
    a * 1 => a
    a * (b * c) => (a * b) * c
end

G = EGraph(:(3 * 4))
analyze!(G, NumberFold)

# exit(0)

@testset "Basic Constant Folding Example - Commutative Monoid" begin
    @test (true == @areequalg G comm_monoid 3 * 4 12)

    @test (true == @areequalg G comm_monoid 3 * 4 12 4*3  6*2)
end

# Metatheory.EGraphs.PRINTIT[] = true
# Metatheory.options.verbose = true
# Metatheory.options.printiter = true

@testset "Basic Constant Folding Example 2 - Commutative Monoid" begin
    ex = :(a * 3 * b * 4)
    G = EGraph(cleanast(ex))
    analyze!(G, NumberFold)
    addexpr!(G, :(12*a))
    println(saturate!(G, comm_monoid))
    display(G.classes); println()
    @test (true == @areequalg G comm_monoid (12*a)*b ((6*2)*b)*a)
    @test (true == @areequalg G comm_monoid (3 * a) * (4 * b) (12*a)*b ((6*2)*b)*a)
end

@testset "Basic Constant Folding Example - Adding analysis after saturation" begin
    G = EGraph(:(3 * 4))
    # addexpr!(G, 12)
    saturate!(G, comm_monoid)
    addexpr!(G, :(a * 2))
    analyze!(G, NumberFold)
    saturate!(G, comm_monoid)

    # display(G.classes); println()
    # println(G.root)
    # display(G.analyses[1].data); println()

    @test (true == areequal(G, comm_monoid, :(3 * 4), 12, :(4*3), :(6*2)))

    ex = :(a * 3 * b * 4)
    G = EGraph(cleanast(ex))
    analyze!(G, NumberFold)
    params=SaturationParams(timeout=15)
    @test areequal(G, comm_monoid, :((3 * a) * (4 * b)), :((12*a)*b),
        :(((6*2)*b)*a); params=params)
end

@testset "Infinite Loops analysis" begin
    boson = @theory begin
        1 * x => x
    end

    G = EGraph(Util.cleanast( :(1 * x) ))
    params=SaturationParams(timeout=100)
    saturate!(G,boson, params)
    ex = extract!(G, ExtractionAnalysis{astsize})

    # println(ex)

    using Metatheory.EGraphs
    boson = @theory begin
        (:c * :cdag) => :cdag * :c + 1
        a * (b + c) => (a * b) + (a * c)
        (b + c) * a => (b * a) + (c * a)
        # 1 * x => x
        (a * b) * c => a * (b * c)
        a * (b * c) => (a * b) * c
    end

    G = EGraph(Util.cleanast( :(c * c * cdag * cdag) ))
    saturate!(G,boson)
    ex = extract!(G, astsize_inv)

    # println(ex)
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct StaticReserveGroup{T <: ReserveDirection} <: Service
        name::String
        available::Bool
        requirement::Float64
        ext::Dict{String, Any}
        contributing_services::Vector{Service}
        internal::InfrastructureSystemsInternal
    end

Data Structure for a group reserve product for system simulations.

# Arguments
- `name::String`
- `available::Bool`
- `requirement::Float64`: the static value of required reserves in system p.u., validation range: `(0, nothing)`, action if invalid: `error`
- `ext::Dict{String, Any}`
- `contributing_services::Vector{Service}`: Services that contribute for this requirement constraint
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct StaticReserveGroup{T <: ReserveDirection} <: Service
    name::String
    available::Bool
    "the static value of required reserves in system p.u."
    requirement::Float64
    ext::Dict{String, Any}
    "Services that contribute for this requirement constraint"
    contributing_services::Vector{Service}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function StaticReserveGroup{T}(name, available, requirement, ext=Dict{String, Any}(), contributing_services=Vector{Service}(), ) where T <: ReserveDirection
    StaticReserveGroup{T}(name, available, requirement, ext, contributing_services, InfrastructureSystemsInternal(), )
end

function StaticReserveGroup{T}(; name, available, requirement, ext=Dict{String, Any}(), contributing_services=Vector{Service}(), internal=InfrastructureSystemsInternal(), ) where T <: ReserveDirection
    StaticReserveGroup{T}(name, available, requirement, ext, contributing_services, internal, )
end

# Constructor for demo purposes; non-functional.
function StaticReserveGroup{T}(::Nothing) where T <: ReserveDirection
    StaticReserveGroup{T}(;
        name="init",
        available=false,
        requirement=0.0,
        ext=Dict{String, Any}(),
        contributing_services=Vector{Service}(),
    )
end

"""Get [`StaticReserveGroup`](@ref) `name`."""
get_name(value::StaticReserveGroup) = value.name
"""Get [`StaticReserveGroup`](@ref) `available`."""
get_available(value::StaticReserveGroup) = value.available
"""Get [`StaticReserveGroup`](@ref) `requirement`."""
get_requirement(value::StaticReserveGroup) = value.requirement
"""Get [`StaticReserveGroup`](@ref) `ext`."""
get_ext(value::StaticReserveGroup) = value.ext
"""Get [`StaticReserveGroup`](@ref) `contributing_services`."""
get_contributing_services(value::StaticReserveGroup) = value.contributing_services
"""Get [`StaticReserveGroup`](@ref) `internal`."""
get_internal(value::StaticReserveGroup) = value.internal

"""Set [`StaticReserveGroup`](@ref) `available`."""
set_available!(value::StaticReserveGroup, val) = value.available = val
"""Set [`StaticReserveGroup`](@ref) `requirement`."""
set_requirement!(value::StaticReserveGroup, val) = value.requirement = val
"""Set [`StaticReserveGroup`](@ref) `ext`."""
set_ext!(value::StaticReserveGroup, val) = value.ext = val


text:
# used for single-line strings
struct StringDiff{T1, T2} <: DeepDiff
    before::T1
    after::T2
    diff::VectorDiff
end

# used for multi-line strings
struct StringLineDiff{T1, T2} <: DeepDiff
    before::T1
    after::T2
    diff::VectorDiff
end

function deepdiff(X::AbstractString, Y::AbstractString)
    if occursin("\n", X) || occursin("\n", Y)
        # we'll compare hashes of each line rather than the text itself, because
        # these comparisons are done many times
        xhashes = map(hash, split(X, '\n'))
        yhashes = map(hash, split(Y, '\n'))

        StringLineDiff(X, Y, deepdiff(xhashes, yhashes))
    else
        StringDiff(X, Y, deepdiff(collect(X), collect(Y)))
    end
end

const AllStringDiffs = Union{StringDiff, StringLineDiff}

before(diff::AllStringDiffs) = diff.before
after(diff::AllStringDiffs) = diff.after
added(diff::AllStringDiffs) = added(diff.diff)
removed(diff::AllStringDiffs) = removed(diff.diff)
changed(diff::AllStringDiffs) = []

Base.:(==)(d1::T, d2::T) where {T<:AllStringDiffs} = fieldequal(d1, d2)

function Base.show(io::IO, diff::StringLineDiff)
    xlines = split(diff.before, '\n')
    ylines = split(diff.after, '\n')
    println(io, "\"\"\"")
    visitall(diff.diff) do idx, state, last
        if state == :removed
            printstyled(io, "- ", escape_string(xlines[idx]), color=:red)
        elseif state == :added
            printstyled(io, "+ ", escape_string(ylines[idx]), color=:green)
        else
            print(io, "  ", escape_string(xlines[idx]))
        end
        if last
            print(io, "\"\"\"")
        else
            println(io)
        end
    end
end

function Base.show(io::IO, diff::StringDiff)
    xchars = before(diff.diff)
    ychars = after(diff.diff)
    laststate = :init

    print(io, "\"")
    visitall(diff.diff) do idx, state, last
        if !hascolor(io)
            # check to see if we need to close a block
            if laststate == :removed && state != :removed
                print(io, "-}")
            elseif laststate == :added && state != :added
                print(io, "+}")
            end
            # check to see if we need to open a block
            if laststate != :removed && state == :removed
                print(io, "{-")
            elseif laststate != :added && state == :added
                print(io, "{+")
            end
        end
        if state == :removed
            printstyled(io, string(xchars[idx]), color=:red)
        elseif state == :added
            printstyled(io, string(ychars[idx]), color=:green)
        else
            print(io, xchars[idx])
        end
        laststate = state
    end
    if !hascolor(io)
        if laststate == :removed
            print(io, "-}")
        elseif laststate == :added
            print(io, "+}")
        end
    end
    print(io, "\"")
end


text:
# # Header
# This is markdown: here's an equation $\frac{1}{2}$

2 + 2

# --
#
# A fragment
#
# ---
#
# # Here's a way to add plots

using Plots; gr()
plot(rand(10))
savefig("statplot.svg");

# ![](statplot.svg)


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

import BenchmarkTools

# Benchmark variables:
repeats = 3;
samples = 1e4;
count = 0;

"""
	print_version()

Prints the TAP version.

# Examples

``` julia
julia> print_version()
```
"""
function print_version()
	@printf( "TAP version 13\n" );
end

"""
	print_summary( total, passing )

Print the benchmark summary.

# Arguments

* `total`: total number of tests
* `passing`: number of passing tests

# Examples

``` julia
julia> print_summary( 3, 3 )
```
"""
function print_summary( total, passing )
	@printf( "#\n" );
	@printf( "1..%d\n", total ); # TAP plan
	@printf( "# total %d\n", total );
	@printf( "# pass  %d\n", passing );
	@printf( "#\n" );
	@printf( "# ok\n" );
end

"""
	print_results( iterations, elapsed )

Print benchmark results.

# Arguments

* `iterations`: number of iterations
* `elapsed`: elapsed time (in seconds)

# Examples

``` julia
julia> print_results( 1000000, 0.131009101868 )
```
"""
function print_results( iterations, elapsed )
	rate = iterations / elapsed

	@printf( "  ---\n" );
	@printf( "  iterations: %d\n", iterations );
	@printf( "  elapsed: %0.9f\n", elapsed );
	@printf( "  rate: %0.9f\n", rate );
	@printf( "  ...\n" );
end

"""
	benchmark( setup, expr )

Run a benchmark.

# Notes

* Benchmark results are returned as a two-element array: [ iterations, elapsed ].
* The number of iterations is not the true number of iterations. Instead, an 'iteration' is defined as a 'sample', which is a computed estimate for a single evaluation.
* The elapsed time is in seconds.

# Arguments

* `setup`: setup expression
* `expr`: expression to benchmark

# Examples

``` julia
julia> out = benchmark( :(), :( sin( 3.14 ) ) );
```
"""
function benchmark( setup, expr )
	t = eval( :( BenchmarkTools.@benchmark $expr samples=$samples setup=($setup) ) )

	# Compute the total "elapsed" time and convert from nanoseconds to seconds:
	s = sum( t.times ) / 1.0e9;

	# Determine the number of "iterations":
	iter = length( t.times );

	# Return the results:
	[ iter, s ];
end

"""
	bench( name, setup, expr )

Run a named benchmark.

# Arguments

* `name`: benchmark name (suffix)
* `setup`: setup expression
* `expr`: expression to benchmark

# Examples

``` julia
julia> bench( "sin", :(), :( sin( 3.14 ) ) );
```
"""
function bench( name, setup, expr )
	for i in 1:repeats
		@printf( "# julia::%s\n", name );
		global count += 1;
		results = benchmark( setup, expr );
		print_results( results[ 1 ], results[ 2 ] );
		@printf( "ok %d benchmark finished\n", count );
	end
end

"""
	main()

Run benchmarks.

# Examples

``` julia
julia> main();
```
"""
function main()
	print_version();

	name = "ndarray::instantiation";
	setup = :();
	stmt = :( A = Array( Float64, 3, 2 ) );
	bench( name, setup, stmt );

	name = "ndarray::get::dtype";
	setup = :( A = Array( Float64, 3, 2 ) );
	stmt = quote
		d = eltype( A );
		if d != Float64
			throw( error( "unexpected value" ) );
		end
	end
	bench( name, setup, stmt );

	name = "ndarray::get::length";
	setup = :( A = Array( Float64, 3, 2 ) );
	stmt = quote
		N = length( A );
		if N != 6
			throw( error( "unexpected value" ) );
		end
	end
	bench( name, setup, stmt );

	name = "ndarray::get::ndims";
	setup = :( A = Array( Float64, 3, 2 ) );
	stmt = quote
		M = ndims( A );
		if M != 2
			throw( error( "unexpected value" ) );
		end
	end
	bench( name, setup, stmt );

	name = "ndarray::get::shape";
	setup = :( A = Array( Float64, 3, 2 ) );
	stmt = quote
		shape = size( A );
		if shape[ 1 ] != 3
			throw( error( "unexpected value" ) );
		end
	end
	bench( name, setup, stmt );

	name = "ndarray::get::strides";
	setup = :( A = Array( Float64, 3, 2 ) );
	stmt = quote
		s = strides( A );
		if s[ 1 ] != 1
			throw( error( "unexpected value" ) );
		end
	end
	bench( name, setup, stmt );

	name = "ndarray::1d:get";
	setup = :( A = rand( Float64, 6 ) );
	stmt = quote
		j = Int( floor( rand()*6.0 ) ) + 1;
		v = A[ j ];
		if v != v
			throw( error( "unexpected value" ) );
		end
	end
	bench( name, setup, stmt );

	name = "ndarray::2d:get";
	setup = :( A = rand( Float64, 3, 2 ) );
	stmt = quote
		j = Int( floor( rand()*3.0 ) ) + 1;
		v = A[ j, 2 ];
		if v != v
			throw( error( "unexpected value" ) );
		end
	end
	bench( name, setup, stmt );

	name = "ndarray::3d:get";
	setup = :( A = rand( Float64, 1, 3, 2 ) );
	stmt = quote
		j = Int( floor( rand()*3.0 ) ) + 1;
		v = A[ 1, j, 2 ];
		if v != v
			throw( error( "unexpected value" ) );
		end
	end
	bench( name, setup, stmt );

	name = "ndarray::4d:get";
	setup = :( A = rand( Float64, 1, 1, 3, 2 ) );
	stmt = quote
		j = Int( floor( rand()*3.0 ) ) + 1;
		v = A[ 1, 1, j, 2 ];
		if v != v
			throw( error( "unexpected value" ) );
		end
	end
	bench( name, setup, stmt );

	name = "ndarray::5d:get";
	setup = :( A = rand( Float64, 1, 1, 1, 3, 2 ) );
	stmt = quote
		j = Int( floor( rand()*3.0 ) ) + 1;
		v = A[ 1, 1, 1, j, 2 ];
		if v != v
			throw( error( "unexpected value" ) );
		end
	end
	bench( name, setup, stmt );

	name = "ndarray::1d:set";
	setup = :( A = rand( Float64, 6 ) );
	stmt = quote
		j = Int( floor( rand()*6.0 ) ) + 1;
		v = rand();
		A[ j ] = v;
		if A[ j ] != v
			throw( error( "unexpected value" ) );
		end
	end
	bench( name, setup, stmt );

	name = "ndarray::2d:set";
	setup = :( A = rand( Float64, 3, 2 ) );
	stmt = quote
		j = Int( floor( rand()*3.0 ) ) + 1;
		v = rand();
		A[ j, 2 ] = v;
		if A[ j, 2 ] != v
			throw( error( "unexpected value" ) );
		end
	end
	bench( name, setup, stmt );

	name = "ndarray::3d:set";
	setup = :( A = rand( Float64, 1, 3, 2 ) );
	stmt = quote
		j = Int( floor( rand()*3.0 ) ) + 1;
		v = rand();
		A[ 1, j, 2 ] = v;
		if A[ 1, j, 2 ] != v
			throw( error( "unexpected value" ) );
		end
	end
	bench( name, setup, stmt );

	name = "ndarray::4d:set";
	setup = :( A = rand( Float64, 1, 1, 3, 2 ) );
	stmt = quote
		j = Int( floor( rand()*3.0 ) ) + 1;
		v = rand();
		A[ 1, 1, j, 2 ] = v;
		if A[ 1, 1, j, 2 ] != v
			throw( error( "unexpected value" ) );
		end
	end
	bench( name, setup, stmt );

	name = "ndarray::5d:set";
	setup = :( A = rand( Float64, 1, 1, 1, 3, 2 ) );
	stmt = quote
		j = Int( floor( rand()*3.0 ) ) + 1;
		v = rand();
		A[ 1, 1, 1, j, 2 ] = v;
		if A[ 1, 1, 1, j, 2 ] != v
			throw( error( "unexpected value" ) );
		end
	end
	bench( name, setup, stmt );

	print_summary( count, count );
end

main();

max_stars_repo_name: andbmme/stdlib

max_issues_repo_name: andbmme/stdlib

max_forks_repo_name: andbmme/stdlib


text: include("FQE.jl")