text:
@testset "selecttransform" begin
    rng = MersenneTwister(123456)

    select = [1, 3, 5]
    t = SelectTransform(select)

    x_vecs = [randn(rng, maximum(select)) for _ in 1:3]
    x_cols = ColVecs(randn(rng, maximum(select), 6))
    x_rows = RowVecs(randn(rng, 4, maximum(select)))

    Xs = [x_vecs, x_cols, x_rows]

    @testset "$(typeof(x))" for x in Xs
        x′ = map(t, x)
        @test all([t(x[n]) == x[n][select] for n in eachindex(x)])
        @test all([t(x[n]) == x′[n] for n in eachindex(x)])
    end

    symbols = [:a, :b, :c, :d, :e]
    select_symbols = [:a, :c, :e]

    ts = SelectTransform(select_symbols)

    a_vecs = map(x -> AxisArray(x; col=symbols), x_vecs)
    a_cols = ColVecs(AxisArray(x_cols.X; col=symbols, index=(1:6)))
    a_rows = RowVecs(AxisArray(x_rows.X; index=(1:4), col=symbols))

    As = [a_vecs, a_cols, a_rows]

    @testset "$(typeof(a))" for (a, x) in zip(As, Xs)
        a′ = map(ts, a)
        x′ = map(t, x)
        @test a′ == x′
    end

    select2 = [2, 3, 5]
    KernelFunctions.set!(t, select2)
    @test t.select == select2

    select_symbols2 = [:b, :c, :e]
    KernelFunctions.set!(ts, select_symbols2)
    @test ts.select == select_symbols2

    @test repr(t) == "Select Transform (dims: $(select2))"
    @test repr(ts) == "Select Transform (dims: $(select_symbols2))"

    test_ADs(() -> SEKernel() ∘ SelectTransform([1, 2]))

    X = randn(rng, (4, 3))
    A = AxisArray(X; row=[:a, :b, :c, :d], col=[:x, :y, :z])
    Y = randn(rng, (4, 2))
    B = AxisArray(Y; row=[:a, :b, :c, :d], col=[:v, :w])
    Z = randn(rng, (2, 3))
    C = AxisArray(Z; row=[:e, :f], col=[:x, :y, :z])

    tx_row = SEKernel() ∘ SelectTransform([1, 2, 4])
    ta_row = SEKernel() ∘ SelectTransform([:a, :b, :d])
    tx_col = SEKernel() ∘ SelectTransform([1, 3])
    ta_col = SEKernel() ∘ SelectTransform([:x, :z])

    @test kernelmatrix(tx_row, X; obsdim=2) ≈ kernelmatrix(ta_row, A; obsdim=2)
    @test kernelmatrix(tx_col, X; obsdim=1) ≈ kernelmatrix(ta_col, A; obsdim=1)

    @test kernelmatrix(tx_row, X, Y; obsdim=2) ≈ kernelmatrix(ta_row, A, B; obsdim=2)
    @test kernelmatrix(tx_col, X, Z; obsdim=1) ≈ kernelmatrix(ta_col, A, C; obsdim=1)

    @testset "$(AD)" for AD in [:Zygote, :ForwardDiff]
        gx = gradient(AD, X) do x
            testfunction(tx_row, x, 2)
        end
        ga = gradient(AD, A) do a
            testfunction(ta_row, a, 2)
        end
        @test gx ≈ ga
        gx = gradient(AD, X) do x
            testfunction(tx_col, x, 1)
        end
        ga = gradient(AD, A) do a
            testfunction(ta_col, a, 1)
        end
        @test gx ≈ ga
        gx = gradient(AD, X) do x
            testfunction(tx_row, x, Y, 2)
        end
        ga = gradient(AD, A) do a
            testfunction(ta_row, a, B, 2)
        end
        @test gx ≈ ga
        gx = gradient(AD, X) do x
            testfunction(tx_col, x, Z, 1)
        end
        ga = gradient(AD, A) do a
            testfunction(ta_col, a, C, 1)
        end
        @test gx ≈ ga
    end

    @testset "$(AD)" for AD in [:ReverseDiff]
        @test_broken ga = gradient(AD, A) do a
            testfunction(ta_row, a, 2)
        end
        @test_broken ga = gradient(AD, A) do a
            testfunction(ta_col, a, 1)
        end
        @test_broken ga = gradient(AD, A) do a
            testfunction(ta_row, a, B, 2)
        end
        @test_broken ga = gradient(AD, A) do a
            testfunction(ta_col, a, C, 1)
        end
    end

    @testset "single-index" begin
        t = SelectTransform(4)
        @testset "$(name)" for (name, x) in [
            ("Vector{<:Vector}", [randn(6) for _ in 1:3]),
            ("ColVecs", ColVecs(randn(5, 10))),
            ("RowVecs", RowVecs(randn(11, 4))),
        ]
            @test KernelFunctions._map(t, x) isa AbstractVector{Float64}
        end
    end
end


text:
"""
Provides the [`mdflatten`](@ref) function that can "flatten" Markdown objects into
a string, with formatting etc. stripped.

Note that the tests in `test/mdflatten.jl` should be considered to be the spec
for the output (number of newlines, indents, formatting, etc.).
"""
module MDFlatten

export mdflatten

import ..Utilities

using Compat

import Base.Markdown:
    MD, BlockQuote, Bold, Code, Header, HorizontalRule,
    Image, Italic, LaTeX, LineBreak, Link, List, Paragraph, Table,
    Footnote, Admonition

"""
Convert a Markdown object to a `String` of only text (i.e. not formatting info).

It drop most of the extra information (e.g. language of a code block, URLs)
and formatting (e.g. emphasis, headers). This "flattened" representation can
then be used as input for search engines.
"""
function mdflatten(md)
    io = IOBuffer()
    mdflatten(io, md)
    Utilities.takebuf_str(io)
end

mdflatten(io, md) = mdflatten(io, md, md)
mdflatten(io, md::MD, parent) = mdflatten(io, md.content, md)
mdflatten(io, vec::Vector, parent) = map(x -> mdflatten(io, x, parent), vec)
function mdflatten(io, vec::Vector, parent::MD)
    # this special case separates top level blocks with newlines
    for md in vec
        mdflatten(io, md, parent)
        print(io, "\n\n")
    end
end

# Block level MD nodes
mdflatten(io, h::Header{N}, parent) where {N} = mdflatten(io, h.text, h)
mdflatten(io, p::Paragraph, parent) = mdflatten(io, p.content, p)
mdflatten(io, bq::BlockQuote, parent) = mdflatten(io, bq.content, bq)
mdflatten(io, ::HorizontalRule, parent) = nothing
function mdflatten(io, list::List, parent)
    for (idx, li) in enumerate(list.items)
        for (jdx, x) in enumerate(li)
            mdflatten(io, x, list)
            jdx == length(li) || print(io, '\n')
        end
        idx == length(list.items) || print(io, '\n')
    end
end
function mdflatten(io, t::Table, parent)
    for (idx, row) = enumerate(t.rows)
        for (jdx, x) in enumerate(row)
            mdflatten(io, x, t)
            jdx == length(row) || print(io, ' ')
        end
        idx == length(t.rows) || print(io, '\n')
    end
end

# Inline nodes
mdflatten(io, text::AbstractString, parent) = print(io, text)
mdflatten(io, link::Link, parent) = mdflatten(io, link.text, link)
mdflatten(io, b::Bold, parent) = mdflatten(io, b.text, b)
mdflatten(io, i::Italic, parent) = mdflatten(io, i.text, i)
mdflatten(io, i::Image, parent) = print(io, "(Image: $(i.alt))")
mdflatten(io, m::LaTeX, parent) = print(io, replace(m.formula, r"[^()+\-*^=\w\s]" => ""))
mdflatten(io, ::LineBreak, parent) = print(io, '\n')

# Is both inline and block
mdflatten(io, c::Code, parent) = print(io, c.code)

# Special (inline) "node" -- due to JuliaMark's interpolations
mdflatten(io, expr::Union{Symbol,Expr}, parent) = print(io, expr)

mdflatten(io, f::Footnote, parent) = footnote(io, f.id, f.text, f)
footnote(io, id, text::Nothing, parent) = print(io, "[$id]")
function footnote(io, id, text, parent)
    print(io, "[$id]: ")
    mdflatten(io, text, parent)
end

function mdflatten(io, a::Admonition, parent)
    println(io, "$(a.category): $(a.title)")
    mdflatten(io, a.content, a)
end

end

ext: jl


text:
module VortexHelperBubbleWrapBlock
using ..Ahorn, Maple

@mapdef Entity "VortexHelper/BubbleWrapBlock" BubbleWrapBlock(x::Integer, y::Integer, width::Integer=8, height::Integer=8, canDash::Bool=true, respawnTime::Number=3.0)

const placements = Ahorn.PlacementDict(
    "Respawning Dash Block (Vortex Helper)" => Ahorn.EntityPlacement(
        BubbleWrapBlock,
		"rectangle"
    )
)

Ahorn.minimumSize(entity::BubbleWrapBlock) = 16, 16
Ahorn.resizable(entity::BubbleWrapBlock) = true, true

function Ahorn.selection(entity::BubbleWrapBlock)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

	return Ahorn.Rectangle(x, y, width, height)
end

function renderBubbleBlock(ctx::Ahorn.Cairo.CairoContext, x::Number, y::Number, width::Number, height::Number)
    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    frame = "objects/VortexHelper/bubbleWrapBlock/bubbleBlock"

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y, 8, 0, 8, 8)
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y + height - 8, 8, 16, 8, 8)
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, x, y + (i - 1) * 8, 0, 8, 8, 8)
        Ahorn.drawImage(ctx, frame, x + width - 8, y + (i - 1) * 8, 16, 8, 8, 8)
    end

    for i in 2:tilesWidth - 1, j in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y + (j - 1) * 8, 8, 8, 8, 8)
    end

    Ahorn.drawImage(ctx, frame, x, y, 0, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, x + width - 8, y, 16, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, x, y + height - 8, 0, 16, 8, 8)
    Ahorn.drawImage(ctx, frame, x + width - 8, y + height - 8, 16, 16, 8, 8)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::BubbleWrapBlock, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

	renderBubbleBlock(ctx, x, y, width, height)
end

end

max_stars_repo_name: catapillie/VortexHelper

max_issues_repo_name: catapillie/VortexHelper

max_forks_repo_name: catapillie/VortexHelper


text:
export
    AVDeviceRect,
    AVAppToDevMessageType,
    AV_APP_TO_DEV_NONE,
    AV_APP_TO_DEV_WINDOW_SIZE,
    AV_APP_TO_DEV_WINDOW_REPAINT,
    AV_APP_TO_DEV_PAUSE,
    AV_APP_TO_DEV_PLAY,
    AV_APP_TO_DEV_TOGGLE_PAUSE,
    AV_APP_TO_DEV_SET_VOLUME,
    AV_APP_TO_DEV_MUTE,
    AV_APP_TO_DEV_UNMUTE,
    AV_APP_TO_DEV_TOGGLE_MUTE,
    AV_APP_TO_DEV_GET_VOLUME,
    AV_APP_TO_DEV_GET_MUTE,
    AVDevToAppMessageType,
    AV_DEV_TO_APP_NONE,
    AV_DEV_TO_APP_CREATE_WINDOW_BUFFER,
    AV_DEV_TO_APP_PREPARE_WINDOW_BUFFER,
    AV_DEV_TO_APP_DISPLAY_WINDOW_BUFFER,
    AV_DEV_TO_APP_DESTROY_WINDOW_BUFFER,
    AV_DEV_TO_APP_BUFFER_OVERFLOW,
    AV_DEV_TO_APP_BUFFER_UNDERFLOW,
    AV_DEV_TO_APP_BUFFER_READABLE,
    AV_DEV_TO_APP_BUFFER_WRITABLE,
    AV_DEV_TO_APP_MUTE_STATE_CHANGED,
    AV_DEV_TO_APP_VOLUME_LEVEL_CHANGED,
    AVDeviceCapabilitiesQuery,
    AVDeviceInfo,
    AVDeviceInfoList

immutable AVDeviceRect
    x::Cint
    y::Cint
    width::Cint
    height::Cint
end

# begin enum AVAppToDevMessageType
typealias AVAppToDevMessageType Uint32
const AV_APP_TO_DEV_NONE = @compat UInt32(1313820229)
const AV_APP_TO_DEV_WINDOW_SIZE = @compat UInt32(1195724621)
const AV_APP_TO_DEV_WINDOW_REPAINT = @compat UInt32(1380274241)
const AV_APP_TO_DEV_PAUSE = @compat UInt32(1346458912)
const AV_APP_TO_DEV_PLAY = @compat UInt32(1347174745)
const AV_APP_TO_DEV_TOGGLE_PAUSE = @compat UInt32(1346458964)
const AV_APP_TO_DEV_SET_VOLUME = @compat UInt32(1398165324)
const AV_APP_TO_DEV_MUTE = @compat UInt32(541939028)
const AV_APP_TO_DEV_UNMUTE = @compat UInt32(1431131476)
const AV_APP_TO_DEV_TOGGLE_MUTE = @compat UInt32(1414354260)
const AV_APP_TO_DEV_GET_VOLUME = @compat UInt32(1196838732)
const AV_APP_TO_DEV_GET_MUTE = @compat UInt32(1196250452)
# end enum AVAppToDevMessageType

# begin enum AVDevToAppMessageType
typealias AVDevToAppMessageType Uint32
const AV_DEV_TO_APP_NONE = @compat UInt32(1313820229)
const AV_DEV_TO_APP_CREATE_WINDOW_BUFFER = @compat UInt32(1111708229)
const AV_DEV_TO_APP_PREPARE_WINDOW_BUFFER = @compat UInt32(1112560197)
const AV_DEV_TO_APP_DISPLAY_WINDOW_BUFFER = @compat UInt32(1111771475)
const AV_DEV_TO_APP_DESTROY_WINDOW_BUFFER = @compat UInt32(1111770451)
const AV_DEV_TO_APP_BUFFER_OVERFLOW = @compat UInt32(1112491596)
const AV_DEV_TO_APP_BUFFER_UNDERFLOW = @compat UInt32(1112884812)
const AV_DEV_TO_APP_BUFFER_READABLE = @compat UInt32(1112687648)
const AV_DEV_TO_APP_BUFFER_WRITABLE = @compat UInt32(1113018912)
const AV_DEV_TO_APP_MUTE_STATE_CHANGED = @compat UInt32(1129141588)
const AV_DEV_TO_APP_VOLUME_LEVEL_CHANGED = @compat UInt32(1129729868)
# end enum AVDevToAppMessageType

immutable AVDeviceCapabilitiesQuery
    av_class::Ptr{AVClass}
    device_context::Ptr{AVFormatContext}
    codec::AVCodecID
    sample_format::AVSampleFormat
    pixel_format::AVPixelFormat
    sample_rate::Cint
    channels::Cint
    channel_layout::Int64
    window_width::Cint
    window_height::Cint
    frame_width::Cint
    frame_height::Cint
    fps::AVRational
end

immutable AVDeviceInfo
    device_name::Ptr{Uint8}
    device_description::Ptr{Uint8}
end

immutable AVDeviceInfoList
    devices::Ptr{Ptr{AVDeviceInfo}}
    nb_devices::Cint
    default_device::Cint
end


text:
using DelayEmbeddings, Random
export ShuffleDimensions

"""
    ShuffleDimensions()

Multidimensional surrogates of input *datasets* (`DelayEmbeddings.Dataset`, which are
also multidimensional) that have shuffled dimensions in each point.

These surrogates destroy the state space structure of the dataset and are thus
suited to distinguish deterministic datasets from high dimensional noise.
"""
struct ShuffleDimensions <: Surrogate end

function surrogenerator(x, sd::ShuffleDimensions, rng = Random.default_rng())
    s = copy(x.data)
    @assert x isa Dataset "input `x` must be `DelayEmbeddings.Dataset` for `ShuffleDimensions`"
    return SurrogateGenerator(sd, x, s, nothing, rng)
end

function (sg::SurrogateGenerator{<:ShuffleDimensions})()
    s = sg.s
    for i in 1:length(s)
        @inbounds s[i] = shuffle(sg.rng, s[i])
    end
    return Dataset(s)
end


text:
function φmin(u::Point{d},v::Point{d}) where d
    sq = sqrt.(u.^2+v.^2)
    # p = (u+v-sq)/2
    # J_u = (1 .- u./sq)/2
    # J_v = (1 .- v./sq)/2
    p =    SVector{d}( (v[i]<Inf ? (u[i]+v[i]-sq[i])/2 : u[i]) for i=1:d )
    J_u =  SVector{d}( (v[i]<Inf ? (1.0 - u[i]./sq[i])/2 : 1) for i=1:d )
    J_v =  SVector{d}( (v[i]<Inf ? (1.0 - v[i]./sq[i])/2 : 0) for i=1:d )

    return p, SDiagonal(J_u), SDiagonal(J_v)
end

function φmax(u::Point{d},v::Point{d}) where d
    sq = sqrt.(u.^2+v.^2)
    # p = (u+v+sq)/2
    # J_u = (1 .+ u./sq)/2
    # J_v = (1 .+ v./sq)/2

    p =    SVector{d}( (v[i]>-Inf ? (u[i]+v[i]+sq[i])/2 : u[i]) for i=1:d )
    J_u =  SVector{d}( (v[i]>-Inf ? (1.0 + u[i]./sq[i])/2 : 1) for i=1:d )
    J_v =  SVector{d}( (v[i]>-Inf ? (1.0 + v[i]./sq[i])/2 : 0) for i=1:d )

    return p, SDiagonal(J_u), SDiagonal(J_v)
end

function PhiPhi(f::Point{d}, x::Point{d}, a::Point{d}, b::Point{d}) where d
    y, y_f, y_x = φmin(f,x-a)

    z, z_y, z_x = φmax(y,x-b)
    return z, (z_y*y_f), (z_y*y_x + z_x)
end

function PhiPhi0(f::Point{d}, x::Point{d}, a::Point{d}, b::Point{d}) where d
    y = φmin(f,x-a)[1]
    z = φmax(y,x-b)[1]
    return z
end

function PhiPhi(f::Point{d}, D::SMatrix{d,d,Float64,q}, x::Point{d}, a::Point{d}, b::Point{d}) where d where q
    z, z_f, z_x = PhiPhi(f,x,a,b)
    z, z_f*D + z_x
end

function PhiPhi!(F::Vector{Vector{Point{d}}},
                 X::Vector{Vector{Point{d}}},
                 A::Vector{Vector{Point{d}}},
                 B::Vector{Vector{Point{d}}},
                 D::Vector{v},
                 J::Matrix{w}) where v<:AbstractVector{SMatrix{d,d,Float64,q}} where w<:AbstractVector{SMatrix{d,d,Float64,q}} where d where q

    n_m, n_M = size(J)
    N = length(F[1])
    for i=1:n_m
        for n=1:N
            f = F[i][n]
            x = X[i][n]
            a = A[i][n]
            b = B[i][n]
            z, z_f, z_x = PhiPhi(f,x,a,b)
            F[i][n] = z
            D[i][n] = z_f*D[i][n] + z_x
            for j=1:n_M
                J[i,j][n] = z_f*J[i,j][n]
            end
        end
    end
end

function PhiPhi!(F::Vector{<:AbstractVector{Point{d}}},
                 X::Vector{<:AbstractVector{Point{d}}},
                 A::Vector{<:AbstractVector{Point{d}}},
                 B::Vector{<:AbstractVector{Point{d}}},
                 J::Matrix{w}) where w<:AbstractVector{SMatrix{d,d,Float64,q}} where d where q

    n_m, n_M = size(J)
    N = length(F[1])
    for i=1:n_m
        for n=1:N
            f = F[i][n]
            x = X[i][n]
            a = A[i][n]
            b = B[i][n]
            z, z_f, z_x = PhiPhi(f,x,a,b)
            F[i][n] = z
            for j=1:n_M
                J[i,j][n] = z_f*J[i,j][n]
            end
        end
    end
end

function PhiPhi!(F::Vector{Vector{Point{d}}},
                X::Vector{Vector{Point{d}}},
                A::Vector{Vector{Point{d}}},
                B::Vector{Vector{Point{d}}},
                D::Vector{v}) where v<:AbstractVector{SMatrix{d,d,Float64,q}} where d where q

    n_m = size(D,1)
    N = length(F[1])
    for i=1:n_m
        for n=1:N
            f = F[i][n]
            x = X[i][n]
            a = A[i][n]
            b = B[i][n]
            z, z_f, z_x = PhiPhi(f,x,a,b)
            F[i][n] = z
            D[i][n] = z_f*D[i][n] + z_x
        end
    end
end

function PhiPhi(F::Vector{Vector{Point{d}}},X::Vector{Vector{Point{d}}},A::Vector{Vector{Point{d}}},B::Vector{Vector{Point{d}}},D::Vector{Vector{SMatrix{d,d,Float64,q}}}, J::Matrix{Vector{SMatrix{d,d,Float64,q}}}) where d where q
    FF = deepcopy(F)
    DD = deepcopy(D)
    JJ = deepcopy(J)
    PhiPhi!(FF,X,A,B,DD,JJ)
    return FF,DD,JJ
end


text:
"""
    MOInput(x::AbstractVector, out_dim::Integer)

A data type to accomodate modelling multi-dimensional output data.
"""
struct MOInput{T<:AbstractVector} <: AbstractVector{Tuple{Any,Int}}
    x::T
    out_dim::Integer
end

Base.length(inp::MOInput) = inp.out_dim * length(inp.x)

Base.size(inp::MOInput, d) = d::Integer == 1 ? inp.out_dim * size(inp.x, 1) : 1
Base.size(inp::MOInput) = (inp.out_dim * size(inp.x, 1),)

Base.lastindex(inp::MOInput) = length(inp)
Base.firstindex(inp::MOInput) = 1

function Base.getindex(inp::MOInput, ind::Integer)
    if ind > 0
        out_dim = ind ÷ length(inp.x) + 1
        ind = ind % length(inp.x)
        if ind == 0
            ind = length(inp.x)
            out_dim -= 1
        end
        return (inp.x[ind], out_dim::Int)
    else
        throw(BoundsError(string("Trying to access at ", ind)))
    end
end

Base.iterate(inp::MOInput) = (inp[1], 1)
function Base.iterate(inp::MOInput, state)
    return (state < length(inp)) ? (inp[state + 1], state + 1) : nothing
end


text:
module Matte

import Genie, Genie.Router
import Base64: base64encode
import DataFrames, Dates, JSON, Logging, Pkg, Plots, Tables

import Base: getproperty, setproperty!

# layouts
export sidebar_layout, side_panel, main_panel
export tabs_layout, tab_panel
export footer_control_layout, content_panel, control_panel
export custom_grid_layout, custom_grid_row, custom_grid_column, custom_card
export expansion_panel_list, expansion_panel
export header, footer

# style elements
export p, h1, h2, h3, br, dialog, snackbar, div, span, show_if, hide_if, circular_loader
export divider, icon

# input elements
export slider, text_input, number_input, button, floating_action_button, tooltip, date_picker
export time_picker, selector, checkbox, switch, radio, list, list_item, button_group

# output elements
export text_output, plots_output, datatable_output, update_output

# functions and helpers
export run_app, stop_app
export new_matte_app, matte_example

include("reflection.jl")
include("ui-misc.jl")
include("ui-design-elements.jl")
include("ui-input-elements.jl")
include("ui-layouts.jl")
include("ui-output-elements.jl")
include("ui-update.jl")
include("generate-js.jl")
include("sanity.jl")
include("session.jl")
include("static.jl")
include("template.jl")
include("run.jl")
include("backend.jl")
include("usethis.jl")

end # module


text:
# WIP
# K/Ca/Leak channels. synapses with presyn NT release + d[NT] + postsyn ionotropic receptor

function gen_sys_stateful(wiring; T=Float32, name)
  @variables t
  D = Differential(t)

  n_neurons = wiring.n_total
  n_in = s_in = wiring.s_in
  w_sens = wiring.matrices[:w_sens]
  w_syns = wiring.matrices[:w_syns]

  sts = @variables begin
    (v[1:n_neurons](t) = rand_uniform(T, -80, -40, n_neurons)), [lower=T(-90), upper=T(30)]
    (a_K[1:n_neurons](t) = rand_uniform(T, 0.02, 0.1, n_neurons)), [lower=T(0), upper=T(1)]
    (a_Ca[1:n_neurons](t) = rand_uniform(T, 0.02, 0.1, n_neurons)), [lower=T(0), upper=T(1)]
    (i_Ca[1:n_neurons](t) = rand_uniform(T, 0.5, 0.9, n_neurons)), [lower=T(0), upper=T(1)]

    I_L[1:n_neurons](t)

    a_K_α[1:n_neurons](t)
    a_K_β[1:n_neurons](t)
    I_K[1:n_neurons](t)

    a_Ca_α[1:n_neurons](t)
    a_Ca_β[1:n_neurons](t)
    i_Ca_α[1:n_neurons](t)
    i_Ca_β[1:n_neurons](t)
    I_Ca[1:n_neurons](t)

    I_sens[1:n_neurons](t)
    I_syns[1:n_neurons](t)
  end

  ps = @parameters begin
    stim[1:n_neurons] = fill(T(1), n_neurons)
    C[1:n_neurons] = rand_uniform(T, 1, 1.01, n_neurons), [lower=T(0.5), upper=T(2.1)]
    G_L[1:n_neurons] = rand_uniform(T, 0.001, 0.2, n_neurons), [lower=T(0.001), upper=T(1)]
    E_L[1:n_neurons] = rand_uniform(T, -80, -50, n_neurons), [lower=T(-90), upper=T(-20)]

    a_K_α_A[1:n_neurons] = rand_uniform(T, 0.01, 0.1, n_neurons), [lower=T(0.001), upper=T(1)]
    a_K_α_vh[1:n_neurons] = rand_uniform(T, -40, -30, n_neurons), [lower=T(-80), upper=T(-20)]
    a_K_α_k[1:n_neurons] = rand_uniform(T, 10, 30, n_neurons), [lower=T(1), upper=T(100)]
    a_K_β_A[1:n_neurons] = rand_uniform(T, 2, 9, n_neurons), [lower=T(1), upper=T(10)]
    a_K_β_vh[1:n_neurons] = rand_uniform(T, -50, -40, n_neurons), [lower=T(-80), upper=T(-30)]
    a_K_β_k[1:n_neurons] = rand_uniform(T, 20, 40, n_neurons), [lower=T(10), upper=T(100)]

    a_Ca_α_A[1:n_neurons] = rand_uniform(T, 0.01, 0.1, n_neurons), [lower=T(0.001), upper=T(1)]
    a_Ca_α_vh[1:n_neurons] = rand_uniform(T, -70, -60, n_neurons), [lower=T(-80), upper=T(80)]
    a_Ca_α_k[1:n_neurons] = rand_uniform(T, 2, 20, n_neurons), [lower=T(1), upper=T(40)]
    a_Ca_β_A[1:n_neurons] = rand_uniform(T, 1, 10, n_neurons), [lower=T(0.5), upper=T(20)]
    a_Ca_β_vh[1:n_neurons] = rand_uniform(T, -75, -50, n_neurons), [lower=T(-80), upper=T(80)]
    a_Ca_β_k[1:n_neurons] = rand_uniform(T, 10, 30, n_neurons), [lower=T(1), upper=T(100)]

    i_Ca_α_A[1:n_neurons] = rand_uniform(T, 0.01, 0.1, n_neurons), [lower=T(0.001), upper=T(1)]
    i_Ca_α_vh[1:n_neurons] = rand_uniform(T, -40, -30, n_neurons), [lower=T(-80), upper=T(30)]
    i_Ca_α_k[1:n_neurons] = rand_uniform(T, 10, 30, n_neurons), [lower=T(1), upper=T(100)]
    i_Ca_β_A[1:n_neurons] = rand_uniform(T, 2, 9, n_neurons), [lower=T(1), upper=T(10)]
    i_Ca_β_vh[1:n_neurons] = rand_uniform(T, -60, -50, n_neurons), [lower=T(-70), upper=T(-40)]
    i_Ca_β_k[1:n_neurons] = rand_uniform(T, 20, 40, n_neurons), [lower=T(10), upper=T(100)]

    G_K[1:n_neurons] = rand_uniform(T, 0.001, 0.2, n_neurons), [lower=T(0.001), upper=T(1)]
    E_K[1:n_neurons] = rand_uniform(T, -80, -50, n_neurons), [lower=T(-90), upper=T(-20)]
    G_Ca[1:n_neurons] = rand_uniform(T, 0.001, 0.2, n_neurons), [lower=T(0.001), upper=T(1)]
    E_Ca[1:n_neurons] = rand_uniform(T, 20, 60, n_neurons), [lower=T(10), upper=T(80)]
  end

  sts = reduce(vcat, collect.(sts))
  ps = reduce(vcat, collect.(ps))

  a_K_αf(i) = a_K_α_A[i] * ((v[i] - a_K_α_vh[i]) / (1f0 - exp(-(v[i] - a_K_α_vh[i]) / a_K_α_k[i])))
  a_K_βf(i) = a_K_β_A[i] * exp(-(v[i] - a_K_β_vh[i]) / a_K_β_k[i])

  a_Ca_αf(i) = a_Ca_α_A[i] * ((v[i] - a_Ca_α_vh[i]) / (1f0 - exp(-(v[i] - a_Ca_α_vh[i]) / a_Ca_α_k[i])))
  a_Ca_βf(i) = a_Ca_β_A[i] * exp(-(v[i] - a_Ca_β_vh[i]) / a_Ca_β_k[i])
  i_Ca_αf(i) = i_Ca_α_A[i] * exp(-(v[i] - i_Ca_α_vh[i]) / i_Ca_α_k[i])
  i_Ca_βf(i) = i_Ca_β_A[i] * ((v[i] - i_Ca_β_vh[i]) / (1f0 - exp(-(v[i] - i_Ca_β_vh[i]) / i_Ca_β_k[i])))

  eqs = Equation[]
  for i in 1:n_neurons
    push!(eqs, D(v[i]) ~ -(I_L[i] + I_K[i] + I_Ca[i] + I_sens[i] + I_syns[i]) / C[i])
  end
  for i in 1:n_neurons
    push!(eqs, I_L[i] ~ G_L[i] * (v[i] - E_L[i]))

    push!(eqs, I_K[i] ~ G_K[i] * a_K[i]^4 * (v[i] - E_K[i]))
    push!(eqs, a_K_α[i] ~ a_K_αf(i))
    push!(eqs, a_K_β[i] ~ a_K_βf(i))
    push!(eqs, D(a_K[i]) ~ a_K_α[i] * (1f0 - a_K[i]) - a_K_β[i]*a_K[i])

    push!(eqs, I_Ca[i] ~ G_Ca[i] * a_Ca[i]^3 * i_Ca[i] * (v[i] - E_Ca[i]))
    push!(eqs, a_Ca_α[i] ~ a_Ca_αf(i))
    push!(eqs, a_Ca_β[i] ~ a_Ca_βf(i))
    push!(eqs, D(a_Ca[i]) ~ a_Ca_α[i] * (1f0 - a_Ca[i]) - a_Ca_β[i]*a_Ca[i])

    push!(eqs, i_Ca_α[i] ~ i_Ca_αf(i))
    push!(eqs, i_Ca_β[i] ~ i_Ca_βf(i))
    push!(eqs, D(i_Ca[i]) ~ i_Ca_α[i] * (1f0 - i_Ca[i]) - i_Ca_β[i]*i_Ca[i])
  end

  systems = ODESystem[]

  I_senss = reshape(Num[0 for _ in 1:s_in*n_neurons], s_in, n_neurons)
  I_synss = reshape(Num[0 for _ in 1:n_neurons*n_neurons], n_neurons, n_neurons)

  for d in 1:n_neurons
    for s in 1:n_in
      w_sens[s,d] == 0 && continue

      ss = Symbol("a_nt_release_sens$(s)$(d)")
      a_nt_release_sens = @variables $ss(t)
      ss = Symbol("nt_sens$(s)$(d)")
      nt_sens = @variables ($ss(t) = rand_uniform(T, 0.001, 0.1)), [lower=T(0), upper=T(10)]
      ss = Symbol("a_r_sens$(s)$(d)")
      a_r_sens = @variables $ss(t)

      ss = Symbol("a_nt_release_sens_A$(s)$(d)")
      a_nt_release_sens_A = @parameters $ss = rand_uniform(T, 0.001, 0.9), [lower=T(0.0001), upper=T(2)]
      ss = Symbol("a_nt_release_sens_vh$(s)$(d)")
      a_nt_release_sens_vh = @parameters $ss = rand_uniform(T, -40, 20), [lower=T(-80), upper=T(80)]
      ss = Symbol("a_nt_release_sens_k$(s)$(d)")
      a_nt_release_sens_k = @parameters $ss = rand_uniform(T, 0.1, 1), [lower=T(0.05), upper=T(10)]

      ss = Symbol("nt_sens_t_accum$(s)$(d)")
      nt_sens_t_accum = @parameters $ss = rand_uniform(T, 0.01, 0.1), [lower=T(0.0001), upper=T(10)]
      ss = Symbol("nt_sens_t_elim$(s)$(d)")
      nt_sens_t_elim = @parameters $ss = rand_uniform(T, 0.1, 1), [lower=T(0.0001), upper=T(10)]

      ss = Symbol("rit_sens$(s)$(d)")
      rit_sens = @parameters $ss = rand_uniform(T, 0.99, 0.999), [lower=T(0.1), upper=T(1)]
      ss = Symbol("cit_sens$(s)$(d)")
      cit_sens = @parameters $ss = rand_uniform(T, 0.5, 0.9), [lower=T(0.1), upper=T(1)]

      ss = Symbol("G_sens$(s)$(d)")
      G_sens = @parameters $ss = rand_uniform(T, 0.001, 0.1), [lower=T(0), upper=T(1)]
      ss = Symbol("E_sens$(s)$(d)")
      E_sens = @parameters $ss = wiring.matrices[:E_sens][s,d], [lower=T(-80), upper=T(80)]

      push!(eqs, a_nt_release_sens[1] ~ a_nt_release_sens_A[1] / (1 + exp(-a_nt_release_sens_k[1] * (fff(stim[s]) - a_nt_release_sens_vh[1]))))
      push!(eqs, D(nt_sens[1]) ~ a_nt_release_sens[1] / nt_sens_t_accum[1] - nt_sens[1] / nt_sens_t_elim[1])
      push!(eqs, a_r_sens[1] ~ rit_sens[1] * nt_sens[1] / (cit_sens[1] + nt_sens[1]))

      I_senss[s,d] = G_sens[1] * a_r_sens[1] * (v[d] - E_sens[1])

      push!(sts, a_nt_release_sens[1], nt_sens[1], a_r_sens[1])
      push!(ps, a_nt_release_sens_A[1], a_nt_release_sens_vh[1], a_nt_release_sens_k[1], nt_sens_t_accum[1], nt_sens_t_elim[1], rit_sens[1], cit_sens[1], G_sens[1], E_sens[1])
    end
    for s in 1:n_neurons
      w_syns[s,d] == 0 && continue

      ss = Symbol("a_nt_release_syns$(s)$(d)")
      a_nt_release_syns = @variables $ss(t)
      ss = Symbol("nt_syns$(s)$(d)")
      nt_syns = @variables ($ss(t) = rand_uniform(T, 0.001, 0.1)), [lower=T(0), upper=T(10)]
      ss = Symbol("a_r_syns$(s)$(d)")
      a_r_syns = @variables $ss(t)

      ss = Symbol("a_nt_release_syns_A$(s)$(d)")
      a_nt_release_syns_A = @parameters $ss = rand_uniform(T, 1.01, 10), [lower=T(1), upper=T(1000)]
      ss = Symbol("a_nt_release_syns_vh$(s)$(d)")
      a_nt_release_syns_vh = @parameters $ss = rand_uniform(T, -50, 50), [lower=T(-80), upper=T(80)]
      ss = Symbol("a_nt_release_syns_k$(s)$(d)")
      a_nt_release_syns_k = @parameters $ss = rand_uniform(T, 0.1, 0.8), [lower=T(0.05), upper=T(1)]

      ss = Symbol("nt_syns_t_accum$(s)$(d)")
      nt_syns_t_accum = @parameters $ss = rand_uniform(T, 1, 10), [lower=T(0.0001), upper=T(10)]
      ss = Symbol("nt_syns_t_elim$(s)$(d)")
      nt_syns_t_elim = @parameters $ss = rand_uniform(T, 1, 10), [lower=T(0.0001), upper=T(10)]

      ss = Symbol("rit_syns$(s)$(d)")
      rit_syns = @parameters $ss = rand_uniform(T, 0.99, 0.999), [lower=T(0.1), upper=T(1)]
      ss = Symbol("cit_syns$(s)$(d)")
      cit_syns = @parameters $ss = rand_uniform(T, 0.5, 0.9), [lower=T(0.1), upper=T(1)]

      ss = Symbol("G_syns$(s)$(d)")
      G_syns = @parameters $ss = rand_uniform(T, 0.001, 0.1), [lower=T(0), upper=T(1)]
      ss = Symbol("E_syns$(s)$(d)")
      E_syns = @parameters $ss = wiring.matrices[:E_syns][s,d], [lower=T(-80), upper=T(80)]

      push!(eqs, a_nt_release_syns[1] ~ a_nt_release_syns_A[1] / (1 + exp(-a_nt_release_syns_k[1] * (v[s] - a_nt_release_syns_vh[1]))))
      push!(eqs, D(nt_syns[1]) ~ a_nt_release_syns[1] / nt_syns_t_accum[1] - nt_syns[1] / nt_syns_t_elim[1])
      push!(eqs, a_r_syns[1] ~ rit_syns[1] * nt_syns[1] / (cit_syns[1] + nt_syns[1]))

      I_synss[s,d] = G_syns[1] * a_r_syns[1] * (v[d] - E_syns[1])

      push!(sts, a_nt_release_syns[1], nt_syns[1], a_r_syns[1])
      push!(ps, a_nt_release_syns_A[1], a_nt_release_syns_vh[1], a_nt_release_syns_k[1], nt_syns_t_accum[1], nt_syns_t_elim[1], rit_syns[1], cit_syns[1], G_syns[1], E_syns[1])
    end
  end
  for i in 1:n_neurons
    push!(eqs, I_sens[i] ~ sum(I_senss[:,i]))
    push!(eqs, I_syns[i] ~ sum(I_synss[:,i]))
  end

  ODESystem(eqs, t, sts, ps; systems, name)
end

#Zygote.@nograd gen_sys_stateful

function StatefulMTK(wiring, solver, sensealg;
  T=Float32, tspan=T.((0, 1)),
  #wiring=nothing,
  mtkize=false, gen_jac=false, kwargs...)

  @named sys = gen_sys_stateful(wiring)
  ssys = structural_simplify(sys)
  defs = ModelingToolkit.get_defaults(ssys)
  prob = ODEProblem(ssys, defs, tspan, tgrad=true, jac=true)

  lb, ub = BCTRNN.get_bounds(ssys, wiring.s_in; T)

  rnncell = BCTRNNCell(wiring, solver, sensealg, prob, lb, ub; kwargs...)
  MyRecur(rnncell)
end


text:
using SpecialFunctions; # need to evaluate the Hankel function (Green's function)
using Luxor;  # this package is needed to indentify the points inside/outside the curve
using Plots;
#---------------------------------------------------------------------------------------
struct Curve
    name::String
    N::Int64               # number of nodes
    x::Array{Float64,2}    # nodes on the curve generated by a parametrization
    dx::Array{Float64,2}   # 1st derivative of the curve parametrization at the nodes
    d2x::Array{Float64,2}  # 2nd derivative of the curve parametrization at the nodes
    d3x::Array{Float64,2}  # 3rd derivative of the curve parametrization at the nodes
    τ::Array{Float64,1}    # norm of the derivative of the curve parametrization at the nodes
    ν::Array{Float64,2}    # unit normal at the nodes
    Δs::Float64            # discretization size in parameter space (~2*pi/N)
    s::Array{Float64,1}    # discrete values of the parameter
    ds::Array{Float64,1}   # derivative of the sigmoid change of variable (for curves with corners)
    R::Array{Float64,2}    # Nystrom weights
    param::Function        # Curve parametrization
    PD:: Int64             # Polynomial degree in the mesh refinement
end
#---------------------------------------------------------------------------------------
function SkeletonCctor(ver,dir,M::Array{Int64,1},p=4)
    # ver: coordinates of the vertices
    # dir: orientation of the segment:
    #       dir[1]:index of starting vertex
    #       dir[2]:index of ending vertex
    # M: roughly half of the number of nodes in the segment
    # p: polinomial degree of the sigmoid change of variables

    Np    = size(dir,1);
    parts = Vector{Curve}(undef,Np)

    Mvec = Array{Int64,1}(undef,Np)
    if length(M)==1
        Mvec .= M
    else
        Mvec = M
    end

    Npts = 0
    for n=1:Np
        st = dir[n,1];ed = dir[n,2];
        parts[n] = lineSegment(Mvec[n],ver[st,:],ver[ed,:],p);
        Npts = Npts+parts[n].N
    end

    plt = plot(parts[1].x[:,1],parts[1].x[:,2]);
    for n =2:Np
        plt = plot!(parts[n].x[:,1],parts[n].x[:,2]);
    end
    plt = plot!(plt,legend=false,aspect_ratio=:equal,frame=:box)
    # return
    return s = (parts=parts,vertices=ver,orientation=dir,Nparts=Np,plt=plt,Npts = Npts)

end
#---------------------------------------------------------------------------------------
function NystromWeights(M)
    R = zeros(2*M,2*M);
        for p=1:2*M
            tp = pi/M*(p-1);
            for j=1:2*M
                tj = pi/M*(j-1);
                R[p,j]=-2*pi/M*sum((ones(M-1)./collect(1:M-1)).*cos.(collect(1:M-1)*(tp-tj)))-pi/(M^2)*cos(M*(tp-tj));
            end
        end
    return R
end

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

function lineSegment(M::Int64,xs,xe,p=4)
    Δs = pi/M
    t = Δs*collect(1:2*M-1)
    N = 2*M-1
    (s,ds) = polySigmoid(t,p)

    param(s) = ([xs[1].+(xe[1]-xs[1])/2/pi*s xs[2].+(xe[2]-xs[2])/2/pi*s],1/2/pi*[(xe[1]-xs[1])*size(s,1) (xe[2]-xs[2])*size(s,1)])

    # Segment
    x   = [xs[1].+(xe[1]-xs[1])/2/pi*s xs[2].+(xe[2]-xs[2])/2/pi*s]
    dx  = 1/2/pi*[(xe[1]-xs[1])*ones(N,1) (xe[2]-xs[2])*ones(N,1)]
    d2x = zeros(N,2)
    d3x = zeros(N,2)

    τ  = sqrt.(dx[:,1].^2+dx[:,2].^2);
    ν = [dx[:,2]./τ -dx[:,1]./τ]

    c = Curve("line",N,x,dx,d2x,d3x,τ,ν,Δs,s,ds,NystromWeights(M),param,p)

    return c

end

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

function kite(M::Integer,rad=1,loc=[0 0])
    Δs = pi/M
    s = Δs*collect(0:2*M-1);
    N = 2*M
    param(s) =(loc.+rad*[cos.(s)+0.65*cos.(2*s).-0.65 1.5*sin.(s)],rad*[-sin.(s)-1.3*sin.(2*s) 1.5*cos.(s)])
    # kite
    x  =loc.+rad*[cos.(s)+0.65*cos.(2*s).-0.65 1.5*sin.(s)];
    dx = rad*[-sin.(s)-1.3*sin.(2*s) 1.5*cos.(s)];
    d2x = rad*[-cos.(s)-2.6*cos.(2*s) -1.5*sin.(s)];
    d3x = rad*[sin.(s)+5.2*sin.(2*s) -1.5*cos.(s)];

    τ  = sqrt.(dx[:,1].^2+dx[:,2].^2);
    ν = [dx[:,2]./τ -dx[:,1]./τ]
    ds = ones(N);

    c = Curve("kite",N,x,dx,d2x,d3x,τ,ν,Δs,s,ds,NystromWeights(M),param,0);

    return c
end
# ---------------------------------------------------------------------------------------
function matricesSkeleton(prtId::Int,sk,kVec,coeff=[1.0 -1;1 -1;1 -1;1 -1])

    prt = sk.parts[prtId]

    if length(kVec)==1

        k=kVec

        if 1==prtId
            matAux = matricesNear(k,prt)
        else
            matAux = matricesFar(k,prt,sk.parts[1])
        end

        S  = matAux.S;
        K  = matAux.K;
        Kp = matAux.Kp;
        T  = matAux.T;

        Np = sk.Nparts;

        for n=2:Np
            if n==prtId
                matAux = matricesNear(k,prt)
            else
                matAux = matricesFar(k,prt,sk.parts[n])
            end
            #println(size(matAux.S))

            S  =[S matAux.S];
            K  =[K matAux.K];
            Kp =[Kp matAux.Kp];
            T  =[T matAux.T];

        end

        return mat = (S=S,K=K,Kp=Kp,T=T)

    elseif length(kVec)==2

        matAux1 = matricesSkeleton(prtId,sk,kVec[1])
        matAux2 = matricesSkeleton(prtId,sk,kVec[2])

        S  = coeff[1,1]*matAux1.S  + coeff[1,2]*matAux2.S
        K  = coeff[2,1]*matAux1.K  + coeff[2,2]*matAux2.K
        Kp = coeff[3,1]*matAux1.Kp + coeff[3,2]*matAux2.Kp
        T  = coeff[4,1]*matAux1.T  + coeff[4,2]*matAux2.T

        return mat = (S=S,K=K,Kp=Kp,T=T)

    end

end
# ---------------------------------------------------------------------------------------
function matricesNear(k,prt)

    τ = prt.τ
    ν = prt.ν
    s = prt.s
    ds= prt.ds
    R = prt.R

    # R = NystromWeights(prt.N)

    N = prt.N
    Δs = prt.Δs

    x1 = prt.x[:,1];     x2 = prt.x[:,2]
    dx1 = prt.dx[:,1];   dx2 = prt.dx[:,2]
    d2x1 = prt.d2x[:,1]; d2x2 = prt.d2x[:,2]
    d3x1 = prt.d3x[:,1]; d3x2 = prt.d3x[:,2]
    n1 = ν[:,1];n2=ν[:,2]

    S  = zeros(N,N)+im*zeros(N,N)
    K  = zeros(N,N)+im*zeros(N,N)
    Kp = zeros(N,N)+im*zeros(N,N)
    T  = zeros(N,N)+im*zeros(N,N)

    psi = -0.577215664901532;

    ii = floor(N/2)
    for i=1:N
        for j=1:N

            r = sqrt((x1[i]-x1[j])^2+(x2[i]-x2[j])^2);

            J0 = besselj0(k*r); J1 = besselj1(k*r)
            H0 = hankelh1(0,k*r); H1 = hankelh1(1,k*r)

            lSin = log(4*(sin(Δs*(i-j)/2))^2)

            if i!=j

                S1 = (-1/(4*pi))*J0*τ[j];
                S2 = im/4*H0*τ[j]-S1*lSin

                K1 = -k/(4*pi)*J1/r*(dx2[j]*(x1[i]-x1[j])-dx1[j]*(x2[i]-x2[j]));
                K2 = im*k/4*H1/r*((x1[i]-x1[j])*dx2[j]-(x2[i]-x2[j])*dx1[j])-K1*lSin

                Kp1 = k/(4*pi)*J1/r*(dx2[i]*(x1[i]-x1[j])-dx1[i]*(x2[i]-x2[j]))*τ[j]/τ[i]
                Kp2 = im*k/4*H1/r*((x1[j]-x1[i])*dx2[i]-(x2[j]-x2[i])*dx1[i])*τ[j]/τ[i]-Kp1*lSin

                Nt = ((x1[i]-x1[j])*dx1[j]+(x2[i]-x2[j])*dx2[j])*((x1[i]-x1[j])*dx1[i]+(x2[i]-x2[j])*dx2[i])/r^2

                NS1 = -1/2/pi*Nt*(k^2*J0-2*k*J1/r)-k/2/pi*(dx1[i]*dx1[j]+dx2[i]*dx2[j])/r*J1

                NS = im/2*Nt*(k^2*H0-2*k*H1/r)+im*k/2*(dx1[i]*dx1[j]+dx2[i]*dx2[j])/r*H1+1/4/pi/(sin((s[i]-s[j])/2))^2

                NS2 = NS-NS1*lSin

            else

                S1 = (-1/(4*pi))*J0*τ[j];
                S2  =(im/4+psi/2/pi-1/4/pi*log(k^2/4*τ[j]^2))*τ[j]+2*log(ds[j])*S1

                K1 = 0.0
                K2 = -1/(4*pi)*(dx1[j]*d2x2[j]-dx2[j]*d2x1[j])/τ[j]^2+2*log(ds[j])*K1

                Kp1 = K1; Kp2 = K2

                NS1 = -k^2/4/pi*τ[i]^2

                NS2 = k^2/4/pi*(pi*im-1+2*psi-2*log(k*τ[i]/2))*τ[i]^2+1/12/pi+
                    +1/2/pi*(dx1[i]*d2x1[i]+dx2[i]*d2x2[i])^2/τ[i]^4 +
                    -1/4/pi*(d2x1[i]^2+d2x2[i]^2)/τ[i]^2 +
                    -1/6/pi*(dx1[i]*d3x1[i]+dx2[i]*d3x2[i])/τ[i]^2

                NS2 = NS2 + 2*log(ds[j])*NS1

            end

            S[i,j] = (R[i,j]*S1 + Δs*S2)*ds[j];
            K[i,j] = (R[i,j]*K1 + Δs*K2)*ds[j];
            Kp[i,j] = (R[i,j]*Kp1 + Δs*Kp2)*ds[j];

            T[i,j] = -(R[i,j]*NS1 + Δs*NS2)*ds[j]/2/τ[i];
            T[i,j] = k^2*S[i,j]*(ν[i,1]*ν[j,1]+ν[i,2]*ν[j,2]) + T[i,j];

        end

    end

    mat = (S=S,K=K,Kp=Kp,T=T);

    plot(plt)
    return mat

end
# ---------------------------------------------------------------------------------------
function matricesFar(k,prtEvl,prtInt)

    τI = prtInt.τ
    τE = prtEvl.τ

    νI = prtInt.ν
    νE = prtEvl.ν

    ds = prtInt.ds
    Δs = prtInt.Δs

    NI = prtInt.N
    NE = prtEvl.N

    x1I = prtInt.x[:,1]
    x2I = prtInt.x[:,2]

    x1E = prtEvl.x[:,1]
    x2E = prtEvl.x[:,2]

    dx1I = prtInt.dx[:,1]
    dx2I = prtInt.dx[:,2]

    dx1E = prtEvl.dx[:,1]
    dx2E = prtEvl.dx[:,2]

    S = zeros(NE,NI)+im*zeros(NE,NI)
    K = zeros(NE,NI)+im*zeros(NE,NI)
    Kp = zeros(NE,NI)+im*zeros(NE,NI)
    T = zeros(NE,NI)+im*zeros(NE,NI)

    for i=1:NE,j=1:NI

        r = sqrt((x1E[i]-x1I[j])^2+(x2E[i]-x2I[j])^2)

        H0 = hankelh1(0,k*r)
        H1 = hankelh1(1,k*r)

        S2 = im/4*H0*τI[j];
        K2 = im*k/4*H1/r*((x1E[i]-x1I[j])*νI[j,1]+(x2E[i]-x2I[j])*νI[j,2])*τI[j]

        Kp2 = im*k/4*H1/r*((x1I[j]-x1E[i])*νE[i,1]+(x2I[j]-x2E[i])*νE[i,2])*τI[j]

        # single layer
        S[i,j] = Δs*S2*ds[j];

        # double layer
        K[i,j] = Δs*K2*ds[j];

        # adjoint double layer operator
        Kp[i,j] = Δs*Kp2*ds[j];

        # hipersingular
        # A11 = (im/4)*k^2*(x1E[i]-x1I[j])^2*H0/r^2 +
        #        (im/4)*k*(r^2-2*(x1E[i]-x1I[j])^2)*H1/r^3

        # A12 = -(im/2)*k*(x1E[i]-x1I[j])*(x2E[i]-x2I[j])*H1/r^3 +
        #        (1i/4)*k^2*(x1E[i]-x1I[j])*(x2E[i]-x2I[j])*H0/r^2

        # A21 = -(1i/2)*k*(x1E[i]-x1I[j]).*(x2E[i]-x2I[j])*H1/r^3 +
        #        (1i/4)*k^2*(x1E[i]-x1I[j])*(x2E[i]-x2I[j])*H0/r^2

        # A22 = (im/4)*k^2*(x2E[i]-x2I[j])^2*H0/r^2 +
        #       (im/4)*k*(r^2-2*(x2E[i]-x2I[j])^2)*H1/r^3

        # T[i,j] = (νE[i,1]*(A11*νI[j,1]+ A12*νI[j,2]) + νE[i,2]*(A21*νI[j,1] + A22*νI[j,2]))*Δs*τI[j]*ds[j]

        R1 = x1E[i]-x1I[j];R2 = x2E[i]-x2I[j];

        A11 = (im/4)*k^2*R1^2*H0/r^2+(im/4)*k*(r^2-2*R1^2)*H1/r^3

        A12 = -(im/2)*k*R1*R2*H1/r^3+(im/4)*k^2*R1*R2*H0/r^2

        A21 = -(im/2)*k*R1*R2*H1/r^3+(im/4)*k^2*R1*R2*H0/r^2

        A22 = (im/4)*k^2*R2^2*H0/r^2+(im/4)*k*(r^2-2*R2^2)*H1/r^3

        n1I =νI[j,1];n2I =νI[j,2];
        n1E =νE[i,1];n2E =νE[i,2];

        T[i,j] = (n1E*(A11*n1I + A12*n2I) + n2E*(A21*n1I + A22*n2I))*Δs*ds[j]*τI[j]

    end

    mat = (S=S,K=K,Kp=Kp,T=T);

    return mat

end
# ---------------------------------------------------------------------------------------
function potentials(k,x,prt)

    Npe = size(x,1)
    Npi = prt.N

    SL = Array{Complex{Float64},2}(undef,Npe,Npi)
    DL = Array{Complex{Float64},2}(undef,Npe,Npi)

    x1 = x[:,1];x2 = x[:,2];

    xp1 = prt.x[:,1];xp2 = prt.x[:,2];

    # dxp1 = prt.dx[:,1];dxp2 = prt.dx[:,2];

    ν1 = prt.ν[:,1];ν2 = prt.ν[:,2]

    ds = prt.ds;Δs = prt.Δs

    tol  = 1e-4;

    for i = 1:Npe
        for j=1:Npi
            r = sqrt((x1[i]-xp1[j])^2+(x2[i]-xp2[j])^2)
            if r<tol
                r=tol
            end
            SL[i,j] = im/4*hankelh1(0,k*r)*prt.τ[j]*prt.ds[j]*Δs
            DL[i,j] = im*k/4*hankelh1(1,k*r)/r*((x1[i]-xp1[j])*ν1[j]+(x2[i]-xp2[j])*ν2[j])*prt.τ[j]*prt.ds[j]*Δs
        end
    end

    pot = (SL=SL,DL=DL);
    return pot
end

#--------------------------------------------------------
function potEval(φ,k,lims,h,crv,coeff)

    xLim = lims[1:2]
    yLim = lims[3:4]

    xa = xLim[1]:h:xLim[2]
    ya = yLim[1]:h:yLim[2]

    Nx = size(xa,1)
    Ny = size(ya,1)

    X = repeat(reshape(xa, 1, :), Ny, 1)
    Y = repeat(ya, 1, Nx)

    pts = [reshape(X,Nx*Ny,1) reshape(Y,Nx*Ny,1)];

    pot = potentials(k,pts,crv);

    Pφ = coeff[1]*pot.DL*φ+coeff[2]*pot.SL*φ;

    Pφ_in  = im*zeros(Nx*Ny,1)
    Pφ_out = im*zeros(Nx*Ny,1)

    # Find points inside the curve crv
    polygon = SVector.(crv.x[:,1],crv.x[:,2])
    points = vec(SVector.(xa',ya))
    inside = [isinside(p, polygon; allowonedge=true) for p in points]

    for j=1:Nx*Ny
        if inside[j]==true
            Pφ_in[j] = Pφ[j]
            Pφ_out[j] = NaN
        else
            Pφ_out[j] = Pφ[j]
            Pφ_in[j] = NaN
        end
    end

    potPlot=(In=reshape(Pφ_in,Ny,Nx),Out=reshape(Pφ_out,Ny,Nx),x=xa,y=ya,X=X,Y=Y,Nx=Nx,Ny=Ny);

    return potPlot
end

#--------------------------------------------------------
function solPlot(sol,crv,opt,reg,cscale=(-1,1));
    if opt=="real"
        if reg=="interior"
            plt = heatmap(sol.x, sol.y,real.(sol.In),color=:RdBu,clim=cscale)
        elseif reg=="exterior"
            plt = heatmap(sol.x, sol.y,real.(sol.Out),color=:RdBu,clim=cscale)
        end
    elseif opt=="abs"
        if reg=="interior"
            plt = heatmap(sol.x, sol.y,abs.(sol.In),clim=cscale)
        elseif reg=="exterior"
            plt = heatmap(sol.x, sol.y,abs.(sol.Out),clim=cscale)
        end
    end

    pX = [crv.x[:,1];crv.x[1,1]];pY = [crv.x[:,2];crv.x[1,2]]
    plt = plot!(pX,pY,aspect_ratio=:equal,lw=1.5,linecolor=:black,legend=false,frame=:box)

    return plt

end

#--------------------------------------------------------
function polySigmoid(s,p);
    #  Change of variable that cancels the singularity at the end points
    v(s,p) = (1/p-1/2)*((pi.-s)/pi).^3 .+1/p*(s.-pi)/pi .+ 1/2;
    dv(s,p)= -(3*(1/p-1/2))*(pi.-s).^2/pi^3 .+1/(p*pi)

    v0 = v(s,p)
    v1 = v(2*pi.-s,p)

    w = 2*pi*v0.^p./(v0.^p.+v1.^p)

    dv0 = dv(s,p);
    dv1 = dv(2*pi.-s,p);

    dw = 2*p*pi*(v0.^(p-1).*dv0.*v1.^p+v0.^p.*v1.^(p-1).*dv1)./(v0.^p+v1.^p).^2;

    w0  = (w=w,dw=dw)
    return w0;
end
#--------------------------------------------------------

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
#!/usr/bin/env julia

using Compose

rawimg = read(joinpath(@__DIR__,"smiley.png"));

compose(context(),
    rectangle(0.1,0.1,0.1,0.1),
    circle(0.3,0.15,0.05),
    polygon([(0.45,0.1),(0.4,0.2),(0.5,0.2)]),
    ellipse(0.65,0.15,0.1,0.05),
    (context(), stroke("black"),
        line([(0.1,0.3),(0.2,0.3)]),
        curve((0.25,0.35),(0.25,0.25),(0.35,0.25),(0.35,0.35))),
    bitmap("image/png",rawimg,0.4,0.25,0.1,0.1),
    text(0.6,0.3,"hello")) |> SVG("primitives.svg")


text:
using DrWatson
@quickactivate "PECUZAL_Julia"
using DelayEmbeddings

# These tests are only for DelayEmbeddings.jl are not strictly needed for this
# scientific project

ti = time()

include("test_pecuzal_embedding.jl")

ti = time() - ti
println("\nTest took total time of:")
println(round(ti, sigdigits=3), " seconds or ", round(ti/60, sigdigits=3), " minutes")


text:
export ismonoid

ismonoid(T::Type, f::Function) =
    closure(T, f) === Closed &&
    identification(T, f) === Identified &&
    associativity(T, f) === Associative


text:
include("../imports.jl")

@testset "ImagePreprocessing" begin
    encfns = [
        () -> ImagePreprocessing(),
        () -> ImagePreprocessing(buffered=false),
        () -> ImagePreprocessing(T=Float64),
        () -> ImagePreprocessing(augmentations=FastAI.augs_lighting()),
        () -> ImagePreprocessing(C=Gray{N0f8}, means=SVector(0.), stds=SVector(1.)),
    ]
    for encfn in encfns
        enc = encfn()
        block = Image{2}()
        img = rand(RGB{N0f8}, 10, 10)
        testencoding(enc, block, img)

        ctx = Validation()
        outblock = encodedblock(enc, block)
        a = encode(enc, ctx, block, img)
        rimg = decode(enc, ctx, outblock, a)
        if eltype(rimg) <: RGB
            @test img ≈ rimg
        end
    end

    @testset "3D" begin
        enc = ImagePreprocessing()
        block = Image{3}()
        img = rand(RGB{N0f8}, 10, 10, 10)
        testencoding(enc, block, img)

        enc = ImagePreprocessing(buffered=false)
        block = Image{3}()
        img = rand(RGB{N0f8}, 10, 10, 10)
        testencoding(enc, block, img)
    end

    @testset ExtendedTestSet "imagedatasetstats" begin

        @testset ExtendedTestSet "RGB" begin
            data = [zeros(RGB{Float32}, 10, 10), ones(RGB{Float32}, 10, 10)]
            means, stds = imagedatasetstats(data, RGB{N0f8}; progress=false)
            @test means ≈ [0.5, 0.5, 0.5]
            @test stds ≈ [0., 0., 0.]
        end

        @testset ExtendedTestSet "Gray" begin
            data = [zeros(Gray{Float32}, 10, 10), ones(Gray{Float32}, 10, 10)]
            means, stds = imagedatasetstats(data, Gray{N0f8}; progress=false)
            @test means ≈ [0.5]
            @test stds ≈ [0.]
        end

    end

    @testset ExtendedTestSet "setup" begin
        data = [
            zeros(10, 10),
            ones(10, 10),
        ]
        enc = setup(ImagePreprocessing, Image{2}(), data, C = Gray{N0f8})
        @test enc.stats[1] ≈ [0.5]
        @test enc.stats[2] ≈ [0.]
    end
end


text:
function sample!(m::AbstractModel; mcmcsize = 4000, burnin = mcmcsize ÷ 2, thin = 1)
    (; d, f, g, β) = skeleton(m)
    chainf = [zeros(length(f)) for _ in 1:(mcmcsize - burnin) ÷ thin]
    chaing = [zeros(length(g)) for _ in 1:(mcmcsize - burnin) ÷ thin]
    chainβ = [zeros(length(β)) for _ in 1:(mcmcsize - burnin) ÷ thin]
    chainnclus = [0 for _ in 1:(mcmcsize - burnin) ÷ thin]
    for iter in 1:mcmcsize
        step!(m)
        if iter > burnin && iszero((iter - burnin) % thin)
            row = (iter - burnin) ÷ thin
            chainf[row] .= f
            chaing[row] .= g
            chainβ[row] .= β
            chainnclus[row] = length(unique(d))
        end
    end
    return chainf, chainβ, chaing, chainnclus
end

function step!(m::AbstractModel)
    update_f!(m)
    update_d!(m)
    update_r!(m)
    update_n!(m)
    update_β!(m)
    update_atoms!(m) # see interface.jl
end

function update_f!(m::AbstractModel)
    (; f, N1, rmax, s, X1, X1vec, y1, β, ϕ1) = skeleton(m)
    mul!(ϕ1, X1, β)
    @. ϕ1 = 1.0 / (1.0 + exp(ϕ1))
    for i in 1:N1
        f[i] = 0.0
        for j in 1:rmax[]
            wj = get_w(ϕ1[i], s[], j)
            f[i] += wj * kernel_pdf(m, y1[i], X1vec[i], j)
        end
    end
end

function get_w(θ0, s0, j)
    # exp(
    #     - log(j) +
    #     logabsbinomial(j + s0 - 2, j - 1)[1] +
    #     s0 * log(θ0) +
    #     (j - 1) * log(1 - θ0) +
    #     log(_₂F₁(BigFloat(j + s0 - 1), 1, j + 1, 1 - θ0))
    # )
    return θ0 * (1 - θ0)^(j - 1)
end

function update_d!(m::AbstractModel)
    (; d, N0, r, rmax, X0vec, y0) = skeleton(m)
    p = zeros(rmax[])
    for i in 1:N0
        resize!(p, r[i])
        for j in 1:r[i]
            p[j] = eps() + kernel_pdf(m, y0[i], X0vec[i], j)
        end
        p ./= sum(p)
        d[i] = rand(Categorical(p))
    end
    return nothing
end

function update_r!(m::AbstractModel)
    (; d, N0, r, rmax, s, X0, β, ϕ0) = skeleton(m)
    mul!(ϕ0, X0, β)
    @. ϕ0 = 1.0 / (1.0 + exp(ϕ0))
    for i in 1:N0
        zi = rand(Beta(r[i] - d[i] + 1, d[i]))
        vi = rand(Gamma(r[i] + s[] - 1, 1))
        mi = rand(Poisson((1 - ϕ0[i]) * zi * vi))
        r[i] = mi + d[i]
    end
    rmax[] = maximum(r)
    return nothing
end

function update_n!(m::AbstractModel)
    (; d, n, N0, rmax) = skeleton(m)
    while length(n) < rmax[]
        push!(n, 0)
    end
    n .= 0
    for i in 1:N0
        n[d[i]] += 1
    end
    return nothing
end

function update_β!(m::AbstractModel)
    (; r, rmodel) = skeleton(m)
    @. rmodel.y = r - 1
    BNB.step!(Random.GLOBAL_RNG, rmodel)
    return nothing
end

struct Womack <: DiscreteMultivariateDistribution
    D::Int
    ζ::Float64
    p::Vector{Float64}
    punnormalized::Vector{Float64}
    function Womack(D::Int, ζ::Float64)
        p = big.([zeros(D); 1.0])
        for d1 in (D - 1):-1:0
            for d2 in 1:(D - d1)
                p[1 + d1] += (
                    ζ * p[1 + d1 + d2] * binomial(big(d1 + d2), big(d1))
                )
            end
        end
        p /= sum(p)
        punnormalized = copy(p)
        for d1 in 1:D
            p[d1] /= binomial(big(D), big(d1 - 1))
        end
        return new(D, ζ, p, punnormalized)
    end
end

function pdf(d::Womack, g::Vector{Bool})
    return d.p[sum(g) + 1]
end

function logpdf(d::Womack, g::Vector{Bool})
    return log(pdf(d, g))
end
# ok: 2022-03-22


text:
const MetaUndirectedGraph = MetaGraph{<:Any,<:Any,<:SimpleGraph}

SimpleGraph(g::MetaUndirectedGraph) = g.graph

is_directed(::Type{<:MetaUndirectedGraph}) = false

function arrange(g::MetaUndirectedGraph, label_1, label_2)
    vprops = g.vprops
    arrange(g::MetaUndirectedGraph, label_1, label_2, vprops[label_1], vprops[label_2])
end

function arrange(g::MetaUndirectedGraph, label_1, label_2, u, v)
    if u > v
        (label_2, label_1)
    else
        (label_1, label_2)
    end
end


text:
using Documenter, ElectromagneticFields

makedocs(
    sitename = "ElectromagneticFields.jl",
    format = Documenter.HTML(),
    pages = ["Home" => "index.md",
             "Modules" => "modules.md",
            ]
)

deploydocs(
    repo   = "github.com/DDMGNI/ElectromagneticFields.jl.git"
)


text:
# adds information about executed line to standard @time macro
# The code by Tamas Papp from
#   https://discourse.julialang.org/t/modifying-the-time-macro/2790
macro vtime(ex)
    quote
        local stats = Base.gc_num()
        local elapsedtime = time_ns()
        local val = $(esc(ex))
        elapsedtime = time_ns() - elapsedtime
        local diff = Base.GC_Diff(Base.gc_num(), stats)
        Logging.# @debug $(Meta.quot(ex)) # print
        Base.time_print(elapsedtime, diff.allocd, diff.total_time,
                        Base.gc_alloc_count(diff))
        println()
        val
    end
end

# wraps nested loops into one for convenience
#   @uwu for (i, j) in (arr1,arr2)
# is equivalent to
#   for i in arr1
#       for j in arr2
macro uwu(forloop)
    head, body = forloop.args
    iters, iterables = head.args

    stri   = map(string, iters.args)

    strhead = reduce(*,
        map(x -> "for "*x[1]*" in "*string(:($(x[2])))*" ",
        zip(stri, iterables.args))
    )

    strbody = string(body)
    whole = strhead*"\n"*strbody*"\nend"^length(stri)

    Meta.parse(whole)
end


text:
using Weave, Compat
using Base.Test

#Test if running document with and without cache works
isdir("documents/cache") && rm("documents/cache", recursive = true)
weave("documents/chunk_options.noweb", plotlib=nothing, cache=:all)
result = @compat readstring(open("documents/chunk_options.md"))
rm("documents/chunk_options.md")
weave("documents/chunk_options.noweb", plotlib=nothing, cache=:all)
cached_result = @compat readstring(open("documents/chunk_options.md"))
@test result == cached_result

# cache = :user
isdir("documents/cache") && rm("documents/cache", recursive = true)
out = "documents/chunk_cache.md"
Weave.weave("documents/chunk_cache.noweb", plotlib=nothing, cache=:user);
result = @compat readstring(open(out))
rm(out)
Weave.weave("documents/chunk_cache.noweb", plotlib=nothing, cache=:user);
cached_result = @compat readstring(open(out))
@test result == cached_result

if VERSION.minor == 3
  using Gadfly
  isdir("documents/cache") && rm("documents/cache", recursive = true)
  #Caching with Gadfly
  weave("documents/gadfly_formats_test.txt", doctype="tex", plotlib="gadfly", cache=:all)
  result = @compat readstring(open("documents/gadfly_formats_test.tex"))
  rm("documents/gadfly_formats_test.tex")
  weave("documents/gadfly_formats_test.txt", doctype="tex", plotlib="gadfly", cache=:all)
  cached_result = @compat readstring(open("documents/gadfly_formats_test.tex"))
  @test result == cached_result
end


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/

# based on the Javascript program
# optimizations by Kristoffer Carlsson

using Printf

A(i,j) = 1.0 / ( (((i+j)*(i+j+1)) >> 1) + i+1)

function Au!(w, u)
    n = length(u)
    Threads.@threads for i = 1:n
        w[i] = 0
        z = 0.0
        @simd for j = 1:n
           @inbounds z += A(i-1, j-1) * u[j]
        end
        w[i] = z
    end
end

function Atu!(v, w)
    n = length(w)
    Threads.@threads for i = 1:n
        z = 0.0
        @simd for j = 1:n
           @inbounds z += A(j-1,i-1) * w[j]
        end
        v[i] = z
    end
end

function AtAu!(w, v, u)
    Au!(w, u)
    Atu!(v, w)
end

function perf_spectralnorm(n::Int=100)
    u = ones(Float64, n)
    v = zeros(Float64 ,n)
    w = zeros(Float64, n)
    vv = vBv = 0
    for i = 1:10
        AtAu!(w, v, u)
        AtAu!(w, u, v)
    end
    for i = 1:n
        vBv += u[i]*v[i]
        vv += v[i]*v[i]
    end
    return sqrt(vBv/vv)
end

n = parse(Int,ARGS[1])
@printf("%.9f\n", perf_spectralnorm(n))


text:
#=
-------------------------------------------------------------------------------------------------------------
This file runs the 1D Shallow Water Equations using a unified CG/DG method using the AGGP storage
introduced in F.X. Giraldo's Introduction to Element-based Galerkin Methods using
Tensor-Product Bases: Analysis, Algorithms, and Applications.
Note that it only requires the storage of a global mass matrix, inside of TI_LSRK,
the CREATE_RHS routine builds the RHS without using global matrices.

The strategy followed is very similar to that described in Alg. 7.6 except that we do not build the
Flux matrix explicitly. This is all done in CREATE_RHS.jl.

Written by F.X. Giraldo on July 6, 2021.
Department of Applied Mathematics
Naval Postgraduate School
Monterey, CA 93943

For time-integration, we use LSRK45 (4th order, 5-stage).
The linearized shallow water equations run indefinitely using either CG or DG. For the nonlinear
equations, we would need some form of stabilization.

The interpolation points used are the following:
ipoints=1: Lobatto
ipoints=2: Legendre
ipoints=3: Chebyshev
ipoints=4: Equi-spaced

The integration points used are the following:
qpoints=1: Lobatto
qpoints=2: Legendre
-------------------------------------------------------------------------------------------------------------
=#

using Plots
using LinearAlgebra

#Include all files used here
include("QuadraturePoints.jl")
include("element_matrices.jl")
include("create_grid.jl")
include("global_matrices.jl")
include("exact_solution.jl")
include("compute_norms.jl")
include("ti_LSRK.jl")

#Some Constants
DFloat = Float64
machine_zero=eps(DFloat)

#Main function
function main()

    #-----------------------------Only Change these Input parameters---------------------------------#
    Ne=20
    N=4
    Q=N
    ipoints=1
    qpoints=1
    #space_method="CG" #time_final=0.4
    space_method="DG" #time_final=0.4
    case=1 #1=exponential, 2=square wave
    Δ_NL=0
    time_final=DFloat(1.0)
    plot_movie=true
    plot_solution=false
    #-----------------------------Only Change these Input parameters---------------------------------#

    Np=N+1
    Nq=Q+1
    Npoin_dg=Ne*Np
    Npoin_cg=Ne*N + 1
    dt=DFloat(0.0001)
    Courant_max=0.5
    ntime = ceil(Int64, time_final / dt)
    @show(N,Q,Ne,Npoin_cg,Npoin_dg,dt,ntime,case)

    #Select Interpolation Points
    if ipoints == 1
        (ξ,ω) = QuadraturePoints.lobatto_gauss(Np,DFloat) #Lobatto
    elseif ipoints == 2
        (ξ,ω) = QuadraturePoints.legendre_gauss(Np,DFloat) #Legendre
    elseif ipoints == 3
        (ξ,ω) = QuadraturePoints.chebyshev_gauss(Np,DFloat) #Chebyshev
    elseif ipoints == 4
        (ξ,ω) = QuadraturePoints.equispaced_points(Np,DFloat) #Equi-spaced
    end

    #Select Quadrature Points
    if qpoints == 1
        (ξq,ωq) = QuadraturePoints.lobatto_gauss(Nq,DFloat) #Lobatto
    elseif qpoints == 2
        (ξq,ωq) = QuadraturePoints.legendre_gauss(Nq,DFloat) #Legendre
    end

    #Construct Lagrange Polynomials
    (ψ,dψ) = QuadraturePoints.lagrange_basis(Np,Nq,ξ,ξq,DFloat)

    #Construct Local Mass and Differentiation Matrices
    (Me,De,Dwe,Fe) = element_matrices(ψ,dψ,Np,Nq,ωq,DFloat)

    #Construct Grid and (i,e)->I map
    (coord_cg,coord_dg,intma_cg,intma_dg,periodicity_cg,periodicity_dg) = create_grid(Np,Npoin_cg,Npoin_dg,Ne,ξ,DFloat)

    #CG or DG?
    if (space_method == "CG")
        Npoin=Npoin_cg
        coord=copy(coord_cg)
        intma=copy(intma_cg)
        periodicity=copy(periodicity_cg)
    elseif (space_method == "DG")
        Npoin=Npoin_dg
        coord=copy(coord_dg)
        intma=copy(intma_dg)
        periodicity=copy(periodicity_dg)
    end
    @show(space_method,Npoin,Ne,N,Q)

    #Construct Global Mass and Differentiation Matrices (only need Mass)
    (M,~,~,~) = global_matrices(Me,De,Dwe,Fe,intma,coord,Npoin,Ne,Np,periodicity,DFloat)

    #Compute Initial Condition
    time=DFloat(0.0)
    (qe,qb,gravity) = exact_solution(coord,Npoin,time,case,DFloat)
    q0=copy(qe)
    dx_min=coord[2]-coord[1]

    #Recompute DT, NTIME
    umax=DFloat(0.0)
    for i=1:Npoin
        h=q0[i,1]+qb[i]
        u=q0[i,2]/h
        umax=max(umax, abs(u) + sqrt(0.5*h^2) )
    end
    dt=Courant_max*dx_min/umax
    courant=umax*dt/dx_min
    ntime = ceil(Int64, time_final / dt)
    println("Courant = ",courant," dt = ",dt," ntime = ",ntime)

    #Time Integration
    #-----------------------------Students Add their Functions inside ti_LSRK--------------------------#
    (q0,time) = ti_LSRK!(q0,qb,coord,M,Dwe,intma,periodicity,time,ntime,dt,gravity,Δ_NL,space_method,plot_movie,case,DFloat)
    #-----------------------------Students Add their Functions inside ti_LSRK--------------------------#

    #Print Time and Extrema
    println(" ntime = ",ntime," time = ",time," qmax = ",maximum(q0)," qmin = ",minimum(q0))

    #Compute L2 Norm
    (qe,qb,gravity) = exact_solution(coord,Npoin,time,case,DFloat)
    norms = compute_norms(q0,qe,DFloat)
    println(" norm.h = ",norms[1]," norm.U = ",norms[2])

    #Plot Final Solution
     if plot_solution
        p1 = plot(coord,qe[:,1],linecolor = :blue,label = "Analytic Solution",legend = :bottomright,legendfont=font(5))
        plot!(p1,coord,q0[:,1],linestyle = :dashdot,linecolor = :red,label = "$space_method Solution")
        title!(p1,"Height (hs)")
        xlabel!(p1,"x")
        ylabel!(p1,"hs(x,t)")
        plot!(p1,xlims = (first(coord),last(coord)+.1),ylims = (-1,0.5))

        p2 = plot(coord,qe[:,2],linecolor = :blue,label = "Analytic Solution",legend = :bottomright,legendfont=font(5))
        plot!(p2,coord,q0[:,2],linestyle = :dashdot,linecolor = :red,label = "$space_method Solution")
        title!(p2,"Momentum (U)")
        xlabel!(p2,"x")
        ylabel!(p2,"U(x,t)")
        plot!(p2,xlims = (first(coord),last(coord)+.1),ylims = (-1,0.5))

        figure_1 = plot(p1,p2,layout = (1,2))
        display(figure_1)
    end
    println("**Simulation Finished**") #output
end
#Main

#----------------------------------#
# Run the main function
#----------------------------------#
main()


text:
using LinearAlgebra
using Plots

include("cheb.jl")
include("parameters.jl")
include("geometry.jl")
include("profile.jl")
include("build_A.jl")
include("compute_spectrum.jl")
include("mesh.jl")
include("plot_fields.jl")

Omega = 2*π/(24*3600)
lat   = π/4

params=parameters(
     H = 2e3,
    Ly = 600e3,
    f₀ = 2*Omega*cos(lat),
     g = 10,
    Lj = 60e3,
    Uj = 1.0,
    Ny = 200,
    dk = 5e-2,
    kₘ = 2e0
    );

params.Fr = (params.f₀ * params.Lj)^2/(params.g * params.H);
params.Ro =  params.Uj/(params.f₀ * params.Lj)

print("\n")
print("Linear-Stability-Caluculator\n")
print("============================\n\n")
print("Repo:      Linear-Stability-Calculators\n")
print("Code:      ShallowWater/Julia/linear_stability_shallow_water.jl\n")
print("Model:     Rotating Shallow Water\n")
print("Geometry:  Cartesian\n")
print("Structure: Bickley Jet\n")
print("\n")
print("Nondimensional Parameters\n")
print("=========================\n")
print("Fr = ", params.Fr, "\n")
print("Ro = ", params.Ro, "\n")

# Wavenumbers
ks = collect(params.dk:params.dk:params.kₘ) / params.Lj;
Nk = length(ks);

# Grid and basic state
Dy, Dy2, y = geometry(params);
U,  E      = profile(         y, params);

plot_basic_state(y, U, E, params.Ly, "basic_state.png")

# initialize fields to store
Nmodes = 2
     σ = zeros(Nmodes, Nk);
     ω = zeros(Nmodes, Nk);
σmodes = zeros(ComplexF64, 3*params.Ny+1, Nmodes, Nk);

# Compute growth rates
for cnt in 1:Nk
    local k = ks[cnt]
    local A = build_A(k, U, E, Dy, y, params);
    local σ[:,cnt], ω[:,cnt], σmodes[:,:,cnt] = compute_spectrum(A, k, params, Nmodes);
end

plot_growth_rates(ks, σ, Nmodes, "growth_rates_Cartesian_Bickley_jet.png")

mode_number = 1;
      σ_max = maximum(σ[mode_number,:]);
    k_index = sortperm(σ[mode_number,:],rev=true)[1];
          k = ks[k_index]; # pick wavenumber

plot_1D_streamfunction(k_index, k, y, σmodes, mode_number, "modes_1D_streamfunction.png")
plot_2D_streamfunction(k_index, k, y, σmodes, mode_number, "modes_2D_streamfunction.png")
plot_2D_vorticity( Dy, k_index, k, y, σmodes, mode_number, "modes_2D_vorticity.png")


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2021 Varnerlab
# Robert Frederick Smith School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED AS IS, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #

# setup paths -
const _PATH_TO_ROOT = pwd()
const _PATH_TO_SRC = joinpath(_PATH_TO_ROOT, "src")
const _PATH_TO_NETWORK = joinpath(_PATH_TO_SRC, "network")

# use packages -
using DifferentialEquations
using VLModelParametersDB
using DelimitedFiles
using Plots

# include my codes -
include("./src/Parameters.jl")
include("./src/Balances.jl")
include("./src/Kinetics.jl")
include("./src/Control.jl")


text:
using ZipCode
using DataFrames, Unitful
using Base.Test

include("tst_cleanzips.jl")
include("tst_distance.jl")
include("tst_utils.jl")


text:
# This file is part of Kpax3. License is MIT.

function crossover!(R1::Vector{Int},
                    R2::Vector{Int},
                    support::GASupport)
  fill!(support.oi.v, 0)
  fill!(support.oj.v, 0)

  encodepartition!(support.oi.R, R1)
  encodepartition!(support.oj.R, R2)

  # sample element on the diagonal (excluding 1)
  g = StatsBase.sample(2:support.n)

  # convert to linear index
  cutpoint = (g - 1) + support.n * (g - 1)

  # unit 1 is always found in cluster 1
  support.oi.v[1] = 1
  support.oj.v[1] = 1

  # vector of weights to compute probabilities (if needed)
  w = zeros(Int, support.n)

  for a in 2:support.n
    if support.oi.R[a] > cutpoint
      if support.oj.R[a] > cutpoint
        # no problem with this move, swap elements
        (support.oi.R[a], support.oj.R[a]) = (support.oj.R[a], support.oi.R[a])

        support.oi.v[1 + div(support.oi.R[a] - a, support.n)] += 1
        support.oj.v[1 + div(support.oj.R[a] - a, support.n)] += 1
      else
        # oj is receiving a "good" cluster but oi got broken
        support.oj.R[a] = support.oi.R[a]

        # we need to decide where to put a: an existing cluster or a new one?
        copyto!(w, 1, support.oi.v, 1, a - 1)
        w[a] = 1
        g = StatsBase.sample(StatsBase.ProbabilityWeights(w[1:a], a))

        support.oi.R[a] = a + support.n * (g - 1)

        support.oi.v[g] += 1
        support.oj.v[1 + div(support.oj.R[a] - a, support.n)] += 1
      end
    elseif support.oj.R[a] > cutpoint
      # oi is receiving a "good" cluster but oj got broken
      support.oi.R[a] = support.oj.R[a]

      # we need to decide where to put a: an existing cluster or a new one?
      copyto!(w, 1, support.oj.v, 1, a - 1)
      w[a] = 1
      g = StatsBase.sample(StatsBase.ProbabilityWeights(w[1:a], a))

      support.oj.R[a] = a + support.n * (g - 1)

      support.oi.v[1 + div(support.oi.R[a] - a, support.n)] += 1
      support.oj.v[g] += 1
    else
      support.oi.v[1 + div(support.oi.R[a] - a, support.n)] += 1
      support.oj.v[1 + div(support.oj.R[a] - a, support.n)] += 1
    end
  end

  decodepartition!(support.oi.R)
  decodepartition!(support.oj.R)

  nothing
end


text:
include("./imports.jl")

@testset ExtendedTestSet "ellipsoid_level_set" begin
    @testset ExtendedTestSet "ellipsoid_level_set" begin
        image_shape = (10, 10)
        answer = [
            0 0 0 0 0 0 0 0 0 0
            0 0 0 1 1 1 1 1 0 0
            0 0 1 1 1 1 1 1 1 0
            0 1 1 1 1 1 1 1 1 1
            0 1 1 1 1 1 1 1 1 1
            0 1 1 1 1 1 1 1 1 1
            0 1 1 1 1 1 1 1 1 1
            0 1 1 1 1 1 1 1 1 1
            0 0 1 1 1 1 1 1 1 0
            0 0 0 1 1 1 1 1 0 0
        ]
        test = ellipsoid_level_set(image_shape)
        @test test == answer
    end

    @testset ExtendedTestSet "ellipsoid_level_set" begin
        image_shape = (10, 10, 2)
        a1 = [
            0 0 0 0 0 0 0 0 0 0
            0 0 0 0 0 0 0 0 0 0
            0 0 0 0 0 0 0 0 0 0
            0 0 0 0 0 0 0 0 0 0
            0 0 0 0 0 0 0 0 0 0
            0 0 0 0 0 0 0 0 0 0
            0 0 0 0 0 0 0 0 0 0
            0 0 0 0 0 0 0 0 0 0
            0 0 0 0 0 0 0 0 0 0
            0 0 0 0 0 0 0 0 0 0
        ]

        a2 = [
            0 0 0 0 0 0 0 0 0 0
            0 0 0 1 1 1 1 1 0 0
            0 0 1 1 1 1 1 1 1 0
            0 1 1 1 1 1 1 1 1 1
            0 1 1 1 1 1 1 1 1 1
            0 1 1 1 1 1 1 1 1 1
            0 1 1 1 1 1 1 1 1 1
            0 1 1 1 1 1 1 1 1 1
            0 0 1 1 1 1 1 1 1 0
            0 0 0 1 1 1 1 1 0 0
        ]
        answer = cat(a1, a2; dims=3)
        test = ellipsoid_level_set(image_shape)
        @test test == answer
    end
end

@testset ExtendedTestSet "lazy_ellipsoid_level_set" begin
    @testset ExtendedTestSet "lazy_ellipsoid_level_set" begin
        image_shape = (10, 10)
        answer = LazySets.Ellipsoid([5.0, 5.0], Diagonal([25.0, 25.0]))
        test = lazy_ellipsoid_level_set(image_shape)
        @test test == answer
    end

    @testset ExtendedTestSet "lazy_ellipsoid_level_set" begin
        image_shape = (10, 10, 2)
        # a1 = [
        #     0 0 0 0 0 0 0 0 0 0
        #     0 0 1 1 1 1 1 1 0 0
        #     0 1 1 1 1 1 1 1 1 0
        #     0 1 1 1 1 1 1 1 1 0
        #     0 1 1 1 1 1 1 1 1 0
        #     0 1 1 1 1 1 1 1 1 0
        #     0 1 1 1 1 1 1 1 1 0
        #     0 1 1 1 1 1 1 1 1 0
        #     0 0 1 1 1 1 1 1 0 0
        #     0 0 0 0 0 0 0 0 0 0
        # ]

        # a2 = [
        #     0 0 0 0 0 0 0 0 0 0
        #     0 0 1 1 1 1 1 1 0 0
        #     0 1 1 1 1 1 1 1 1 0
        #     0 1 1 1 1 1 1 1 1 0
        #     0 1 1 1 1 1 1 1 1 0
        #     0 1 1 1 1 1 1 1 1 0
        #     0 1 1 1 1 1 1 1 1 0
        #     0 1 1 1 1 1 1 1 1 0
        #     0 0 1 1 1 1 1 1 0 0
        #     0 0 0 0 0 0 0 0 0 0
        # ]
        answer = LazySets.Ellipsoid([5.0, 5.0, 1.0], Diagonal([25.0, 25.0, 1.0]))
        test = lazy_ellipsoid_level_set(image_shape)
        @test test == answer
    end
end

@testset ExtendedTestSet "to_array" begin
    @testset ExtendedTestSet "to_array" begin
        image_shape = (10, 10)
        answer = [
            0 0 0 1 1 1 1 0 0 0
            0 1 1 1 1 1 1 1 1 0
            0 1 1 1 1 1 1 1 1 0
            1 1 1 1 1 1 1 1 1 1
            1 1 1 1 1 1 1 1 1 1
            1 1 1 1 1 1 1 1 1 1
            1 1 1 1 1 1 1 1 1 1
            0 1 1 1 1 1 1 1 1 0
            0 1 1 1 1 1 1 1 1 0
            0 0 0 1 1 1 1 0 0 0
        ]
        ellipse = LazySets.Ellipsoid([5.0, 5.0], Diagonal([25.0, 25.0]))
        test = to_array(image_shape, ellipse)
        @test test == answer
    end

    @testset ExtendedTestSet "lazy_ellipsoid_level_set" begin
        image_shape = (10, 10, 2)
        a1 = [
            0 0 0 0 0 0 0 0 0 0
            0 0 1 1 1 1 1 1 0 0
            0 1 1 1 1 1 1 1 1 0
            0 1 1 1 1 1 1 1 1 0
            0 1 1 1 1 1 1 1 1 0
            0 1 1 1 1 1 1 1 1 0
            0 1 1 1 1 1 1 1 1 0
            0 1 1 1 1 1 1 1 1 0
            0 0 1 1 1 1 1 1 0 0
            0 0 0 0 0 0 0 0 0 0
        ]

        a2 = [
            0 0 0 0 0 0 0 0 0 0
            0 0 1 1 1 1 1 1 0 0
            0 1 1 1 1 1 1 1 1 0
            0 1 1 1 1 1 1 1 1 0
            0 1 1 1 1 1 1 1 1 0
            0 1 1 1 1 1 1 1 1 0
            0 1 1 1 1 1 1 1 1 0
            0 1 1 1 1 1 1 1 1 0
            0 0 1 1 1 1 1 1 0 0
            0 0 0 0 0 0 0 0 0 0
        ]
        answer = cat(a1, a2; dims=3)
        ellipse = LazySets.Ellipsoid([5.0, 5.0, 1.0], Diagonal([25.0, 25.0, 1.0]))
        test = to_array(image_shape, ellipse)
        @test test == answer
    end
end


text:
module HybridSystems

include("automata.jl")
include("switchings.jl")

using MathematicalSystems

export AbstractHybridSystem, HybridSystem

"""
    AbstractHybridSystem

Abstract supertype for a hybrid system.
"""
abstract type AbstractHybridSystem <: MathematicalSystems.AbstractSystem end

"""
    HybridSystem{A, S, R, W, SV<:AbstractVector{S}, RV<:AbstractVector{R}, RW<:AbstractVector{W}} <: AbstractHybridSystem

A hybrid system modelled as a hybrid automaton.

### Fields

- `automaton`  -- hybrid automaton of type `A`.
- `modes`      -- vector of modes of type `S` indexed by the discrete states,
                  both the domain and the dynamic are stored in this field.
                  See [`stateset`](@ref) to access the domain.
- `resetmaps`  -- vector of reset maps of type `R` indexed by the label of the
                  transition, the guard is stored as constraint of the map in
                  this field. See [`stateset`](@ref) to access the guard.
- `switchings` -- vector of switchings of type `W` indexed by the label of the
                  transition, see [`AbstractSwitching`](@ref).
- `ext`        -- dictionary that can be used by extensions.

### Notes

The automaton `automaton` of type `A` models the different discrete states and
the allowed transitions with corresponding labels.

The mode dynamic and domain are stored in a continuous dynamical system of type
`S` in the vector `modes`. They are indexed by the discrete states of the
automaton.

The reset maps and guards are given as a map or a discrete dynamical system of
type `R` in the vector `resetmaps`. They are indexed by the labels of the
corresponding transition.

The switching of type `W` is given in the `switchings` vector, indexed by the
label of the transition.

Additional data can be stored in the `ext` field.

### Examples

See [the Thermostat example](https://github.com/blegat/HybridSystems.jl/blob/master/examples/Thermostat.ipynb).
"""
struct HybridSystem{A, S, R, W, SV<:AbstractVector{S}, RV<:AbstractVector{R}, RW<:AbstractVector{W}} <: AbstractHybridSystem
    automaton::A
    modes::SV
    resetmaps::RV
    switchings::RW
    # Can be used by extensions
    ext::Dict{Symbol, Any}
end
function HybridSystem(a, m, r, s)
    HybridSystem(a, m, r, s, Dict{Symbol, Any}())
end

if VERSION >= v"0.7-"
    Base.broadcastable(hs::HybridSystem) = Ref(hs)
end

# The hybrid system both acts like an automaton and a system

# Automaton
for f in [:state_property, :transition_property, :states, :nstates, :transitiontype, :transitions, :ntransitions, :source, :event, :target, :has_transition, :rem_transition!, :in_transitions, :out_transitions]
    @eval begin
        $f(hs::HybridSystem, args...) = $f(hs.automaton, args...)
    end
end
for f in [:state_property_type, :transition_property_type]
    @eval begin
        $f(::Type{<:HybridSystem{A}}, args...) where {A} = $f(A, args...)
    end
end

# System
export statedim, stateset, inputdim, inputset, guard, assignment, mode, target_mode, resetmap
"""
    statedim(hs::HybridSystem, u::Int)

Returns the dimension of the state space of the system at discrete state `u`.
"""
MathematicalSystems.statedim(hs::HybridSystem, u::Int) = statedim(hs.modes[u])

"""
    stateset(s::AbstractSystem, u::Int)

Returns the set of allowed states of the system at discrete state `u`.
"""
MathematicalSystems.stateset(hs::HybridSystem, u::Int) = stateset(hs.modes[u])

"""
    mode(hs::HybridSystem, u::Int)

Returns the mode of the dynamical system at discrete state `u`.
"""
mode(hs::HybridSystem, u::Int) = hs.modes[u]

"""
    target_mode(hs::AbstractHybridSystem, t)

Returns the target mode for the transition `t`.
"""
target_mode(hs::AbstractHybridSystem, t) = mode(hs, target(hs, t))

"""
    resetmap(hs::HybridSystem, t)

Returns the reset map for the transition `t`.
"""
resetmap(hs::HybridSystem, t) = hs.resetmaps[symbol(hs, t)]

"""
    assignment(hs::HybridSystem, t)

Returns the assignment for the transition `t`.
"""
assignment(hs::HybridSystem{A,S,R,W}, t) where {A, S<:AbstractSystem, R<:AbstractMap, W} = (args...) -> apply(resetmap(hs, t), args...)

"""
    guard(hs::HybridSystem, t)

Returns the guard for the transition `t`.
"""
guard(hs::AbstractHybridSystem, t) = stateset(resetmap(hs, t))

# for completeness, extend the stateset function from MathematicalSystems
# because guards are given as the state constraints of the reset map
MathematicalSystems.stateset(hs::AbstractHybridSystem, t) = guard(hs, t)

"""
    inputdim(s::AbstractSystem, u::Int)

Returns the dimension of the input space of the system at mode `u`.
"""
MathematicalSystems.inputdim(hs::HybridSystem, u::Int) = inputdim(hs.modes[u])

"""
    inputset(s::AbstractSystem, u::Int)

Returns the set of allowed inputs of the system at mode `u`.
"""
MathematicalSystems.inputset(hs::HybridSystem, u::Int) = inputset(hs.modes[u])

"""
    inputset(s::AbstractSystem, t)

Returns the set of allowed inputs for the transition `t`.
"""
MathematicalSystems.inputset(hs::HybridSystem, t) = inputset(hs.resetmaps[symbol(hs, t)])

function Base.show(io::IO, hs::HybridSystem)
    print(io, "Hybrid System with automaton ")
    print(io, hs.automaton)
end

# Particular cases
include("switchedsystems.jl")
include("paths.jl")

end # module


text:
using LinearAlgebra: diagm
using SparseArrays: sparse

normalize_indicator(i::AbstractArray) = i ./ sum(i)
normalize_indicator(i::AbstractArray, dims) = i ./ sum(i, dims=dims)

to_indicator_matrix(i::AbstractVector) = diagm(i)

function union_diagonal!(A::AbstractMatrix{T}) where {T}
    for i in 1:min(size(A)...)
        A[i, i] = one(T)
    end
    return A
end

_generate_slices(n, k) = Base.ntuple(i->(i==n) ? k : Colon(), 2)

"""
    graph_filter(ind, order)

Generate a graph filter of `order`-th order moment depending on the indicator `ind`.

## Arguments
- `ind`: the indicator vector or matrix with values of 0s or 1s.
- `order`: the order of moment for a graph filter. If `order=1`, generates a first moment graph filter,
which is the mean operator. If `order=2`, generates a second moment graph filter, which is the covariance operator.
"""
function graph_filter(ind::AbstractVector, order::Integer)
    if order == 1
        return sparse(normalize_indicator(ind))
    elseif order == 2
        adj = to_indicator_matrix(ind)
        return sparse(normalize_indicator(adj))
    else
        throw(ArgumentError("order other than 1 and 2 is not supported while get $order."))
    end
end

function graph_filter(adj::AbstractMatrix, order::Integer; dims::Integer=1, self::Bool=true)
    if order == 1
        self && union_diagonal!(adj)
        return sparse(normalize_indicator(adj, dims))
    elseif order == 2
        n = dims%2 + 1
        blocks = [graph_filter(vec(adj[_generate_slices(n, k)...]), 2) for k in 1:size(adj, n)]
        return cat(blocks...,dims=(1,2))
    else
        throw(ArgumentError("order other than 1 and 2 is not supported while get $order."))
    end
end

function moment(A::AbstractMatrix, X::AbstractMatrix, order::Integer; dims::Integer=1)
    if order == 1
        return first_moment(A, X; dims=dims)
    elseif order == 2
        return second_moment(A, X; dims=dims)
    else
        throw(ArgumentError("order other than 1 and 2 is not supported while get $order."))
    end
end

function first_moment(A::AbstractMatrix, X::AbstractMatrix; dims::Integer=1)
    adj = A .> 0
    if dims == 1
        G = graph_filter(adj, 1; dims=1)
        return G*X
    elseif dims == 2
        G = graph_filter(adj, 1; dims=2)
        return X*G
    else
        throw(ArgumentError("dims other than 1 and 2 is not supported while get $dims."))
    end
end

# second_moment(A::AbstractMatrix, X::AbstractMatrix; dims::Integer=1) = second_moment(A, X, X; dims=dims)

# function second_moment(A::AbstractMatrix, X::AbstractMatrix, Y::AbstractMatrix; dims::Integer=1)
#     # x' * G_2 * y - (G_1' * x)*(G_1' * y)
#     adj = A .> 0
#     if dims == 1
#         G = graph_filter(adj, 2; dims=1)
#         return X'*G*Y - first_moment(A, X; dims=1)*first_moment(A, Y; dims=1)'
#     elseif dims == 2
#         G = graph_filter(adj, 2; dims=2)
#         return vec(X')'*G*vec(Y') - first_moment(A, X; dims=2)*first_moment(A, Y; dims=2)'
#     else
#         throw(ArgumentError("dims other than 1 and 2 is not supported while get $dims."))
#     end
# end