text:
using CImGui
using CImGui.CSyntax
using CImGui.CSyntax.CStatic
using Printf

"""
    ShowDemoWindowMisc()
"""
function ShowDemoWindowMisc()
    if CImGui.CollapsingHeader("Filtering")
        # static ImGuiTextFilter filter;
        CImGui.Text("Filter usage:\n"*
                    "  \"\"         display all lines\n"*
                    "  \"xxx\"      display lines containing \"xxx\"\n"*
                    "  \"xxx,yyy\"  display lines containing \"xxx\" or \"yyy\"\n"*
                    "  \"-xxx\"     hide lines containing \"xxx\"")
        # filter.Draw()
        lines = ["aaa1.c", "bbb1.c", "ccc1.c", "aaa2.cpp", "bbb2.cpp", "ccc2.cpp", "abc.h", "hello, world"]
        for i = 1:length(lines)
            # filter.PassFilter(lines[i]) && CImGui.BulletText(lines[i])
        end
    end

    if CImGui.CollapsingHeader("Inputs, Navigation & Focus")
        io = CImGui.GetIO()

        CImGui.Text("WantCaptureMouse: $(io.WantCaptureMouse)")
        CImGui.Text("WantCaptureKeyboard: $(io.WantCaptureKeyboard)")
        CImGui.Text("WantTextInput: $(io.WantTextInput)")
        CImGui.Text("WantSetMousePos: $(io.WantSetMousePos)")
        CImGui.Text("NavActive: $(io.NavActive), NavVisible: $(io.NavVisible)")

        if CImGui.TreeNode("Keyboard, Mouse & Navigation State")
            if CImGui.IsMousePosValid()
                CImGui.Text("Mouse pos: ($(io.MousePos.x), $(io.MousePos.y))")
            else
                CImGui.Text("Mouse pos: <INVALID>")
            end
            CImGui.Text("Mouse delta: ($(io.MouseDelta.x), $(io.MouseDelta.y))")
            CImGui.Text("Mouse down:")
            for i = 0:4
                dur = CImGui.Get_MouseDownDuration(io, i)
                dur ≥ 0 || continue
                CImGui.SameLine()
                CImGui.Text(@sprintf("b%d (%.02f secs)", i, dur))
            end
            CImGui.Text("Mouse clicked:")
            for i = 0:4
                CImGui.IsMouseClicked(i) || continue
                CImGui.SameLine()
                CImGui.Text("b$i")
            end
            CImGui.Text("Mouse dbl-clicked:")
            for i = 0:4
                CImGui.IsMouseDoubleClicked(i) || continue
                CImGui.SameLine()
                CImGui.Text("b$i")
            end
            CImGui.Text("Mouse released:")
            for i = 0:4
                CImGui.IsMouseReleased(i) || continue
                CImGui.SameLine()
                CImGui.Text("b$i")
            end
            CImGui.Text(@sprintf("Mouse wheel: %.1f", io.MouseWheel))

            CImGui.Text("Keys down:")
            for i = 0:511
                dur = CImGui.Get_KeysDownDuration(io, i)
                dur ≥ 0 || continue
                CImGui.SameLine()
                txt = @sprintf "%d (%.02f secs)" i dur
                CImGui.Text(txt)
            end
            CImGui.Text("Keys pressed:")
            for i = 0:511
                CImGui.IsKeyPressed(i) || continue
                CImGui.SameLine()
                CImGui.Text("$i")
            end
            CImGui.Text("Keys release:")
            for i = 0:511
                CImGui.IsKeyReleased(i) || continue
                CImGui.SameLine()
                CImGui.Text("$i")
            end
            CImGui.Text(@sprintf("Keys mods: %s%s%s%s", io.KeyCtrl ? "CTRL " : "", io.KeyShift ? "SHIFT " : "", io.KeyAlt ? "ALT " : "", io.KeySuper ? "SUPER " : ""))

            CImGui.Text("NavInputs down:")
            for i = 0:20
                nav = CImGui.Get_NavInputs(io, i)
                nav > 0 || continue
                CImGui.SameLine()
                CImGui.Text(@sprintf("[%d] %.2f", i, nav))
            end
            CImGui.Text("NavInputs pressed:")
            for i = 0:20
                dur = CImGui.Get_NavInputsDownDuration(io, i)
                dur == 0.0 || continue
                CImGui.SameLine()
                CImGui.Text("[$i]")
            end
            CImGui.Text("NavInputs duration:")
            for i = 0:20
                dur = CImGui.Get_NavInputsDownDuration(io, i)
                dur ≥ 0.0 || continue
                CImGui.SameLine()
                CImGui.Text(@sprintf("[%d] %.2f", i, dur))
            end

            CImGui.Button("Hovering me sets the\nkeyboard capture flag")
            CImGui.IsItemHovered() && CImGui.CaptureKeyboardFromApp(true)
            CImGui.SameLine()
            CImGui.Button("Holding me clears the\nthe keyboard capture flag")
            CImGui.IsItemActive() && CImGui.CaptureKeyboardFromApp(false)

            CImGui.TreePop()
        end

        if CImGui.TreeNode("Tabbing")
            CImGui.Text("Use TAB/SHIFT+TAB to cycle through keyboard editable fields.")
            @cstatic buf="dummy"*"\0"^(27) begin
                CImGui.InputText("1", buf, length(buf))
                CImGui.InputText("2", buf, length(buf))
                CImGui.InputText("3", buf, length(buf))
                CImGui.PushAllowKeyboardFocus(false)
                CImGui.InputText("4 (tab skip)", buf, length(buf))
                # CImGui.SameLine(); ShowHelperMarker("Use CImGui.PushAllowKeyboardFocus(bool)\nto disable tabbing through certain widgets.")
                CImGui.PopAllowKeyboardFocus()
                CImGui.InputText("5", buf, length(buf))
            end
            CImGui.TreePop()
        end

        if CImGui.TreeNode("Focus from code")
            focus_1 = CImGui.Button("Focus on 1"); CImGui.SameLine();
            focus_2 = CImGui.Button("Focus on 2"); CImGui.SameLine();
            focus_3 = CImGui.Button("Focus on 3");
            has_focus = Cint(0)

            @cstatic buf="click on a button to set focus"*"\0"^98 begin
                focus_1 && CImGui.SetKeyboardFocusHere()
                CImGui.InputText("1", buf, length(buf))
                CImGui.IsItemActive() && (has_focus = 1;)

                focus_2 && CImGui.SetKeyboardFocusHere()
                CImGui.InputText("2", buf, length(buf))
                CImGui.IsItemActive() && (has_focus = 2;)

                CImGui.PushAllowKeyboardFocus(false)
                focus_3 && CImGui.SetKeyboardFocusHere()
                CImGui.InputText("3 (tab skip)", buf, length(buf))
                CImGui.IsItemActive() && (has_focus = 3;)
                CImGui.PopAllowKeyboardFocus()
            end
            if has_focus != 0
                CImGui.Text("Item with focus: $has_focus")
            else
                CImGui.Text("Item with focus: <none>")
            end

            # use >= 0 parameter to SetKeyboardFocusHere() to focus an upcoming item
            @cstatic f3 = Cfloat[0.0, 0.0, 0.0] begin
                focus_ahead = Cint(-1)
                CImGui.Button("Focus on X") && (focus_ahead = 0;)
                CImGui.SameLine()
                CImGui.Button("Focus on Y") && (focus_ahead = 1;)
                CImGui.SameLine()
                CImGui.Button("Focus on Z") && (focus_ahead = 2;)
                focus_ahead != -1 && CImGui.SetKeyboardFocusHere(focus_ahead)
                CImGui.SliderFloat3("Float3", f3, 0.0, 1.0)
            end
            CImGui.TextWrapped("NB: Cursor & selection are preserved when refocusing last used item in code.")
            CImGui.TreePop()
        end

        if CImGui.TreeNode("Dragging")
            CImGui.TextWrapped("You can use CImGui.GetMouseDragDelta(0) to query for the dragged amount on any widget.")
            for button = 0:2
                CImGui.Text(@sprintf("IsMouseDragging(%d):\n  w/ default threshold: %d,\n  w/ zero threshold: %d\n  w/ large threshold: %d",
                    button, CImGui.IsMouseDragging(button), CImGui.IsMouseDragging(button, 0.0), CImGui.IsMouseDragging(button, 20.0)))
            end
            CImGui.Button("Drag Me")
            if CImGui.IsItemActive()
                # draw a line between the button and the mouse cursor
                draw_list = CImGui.GetWindowDrawList()
                CImGui.PushClipRectFullScreen(draw_list)
                click_pos = CImGui.Get_MouseClickedPos(io, 0)
                CImGui.AddLine(draw_list, click_pos, io.MousePos, CImGui.GetColorU32(CImGui.ImGuiCol_Button), 4.0)
                CImGui.PopClipRect(draw_list)

                # drag operations gets "unlocked" when the mouse has moved past a certain threshold (the default threshold is stored in io.MouseDragThreshold)
                # you can request a lower or higher threshold using the second parameter of IsMouseDragging() and GetMouseDragDelta()
                value_raw = CImGui.GetMouseDragDelta(0, 0.0)
                value_with_lock_threshold = CImGui.GetMouseDragDelta(0)
                mouse_delta = io.MouseDelta
                CImGui.SameLine()
                txt = @sprintf "Raw (%.1f, %.1f), WithLockThresold (%.1f, %.1f), MouseDelta (%.1f, %.1f)" value_raw.x value_raw.y value_with_lock_threshold.x value_with_lock_threshold.y mouse_delta.x mouse_delta.y
                CImGui.Text(txt)
            end
            CImGui.TreePop()
        end

        if CImGui.TreeNode("Mouse cursors")
            mouse_cursors_names = ["Arrow", "TextInput", "ResizeAll", "ResizeNS", "ResizeEW", "ResizeNESW", "ResizeNWSE", "Hand", "NotAllowed"]
            @assert length(mouse_cursors_names) == CImGui.ImGuiMouseCursor_COUNT
            cur = CImGui.GetMouseCursor()
            CImGui.Text("Current mouse cursor = $cur: $(mouse_cursors_names[cur+1])")
            CImGui.Text("Hover to see mouse cursors:")
            CImGui.SameLine()
            ShowHelpMarker("Your application can render a different mouse cursor based on what CImGui.GetMouseCursor() returns. If software cursor rendering (io.MouseDrawCursor) is set ImGui will draw the right cursor for you, otherwise your backend needs to handle it.")
            for i = 0:CImGui.ImGuiMouseCursor_COUNT-1
                label = @sprintf("Mouse cursor %d: %s", i, mouse_cursors_names[i+1])
                CImGui.Bullet()
                CImGui.Selectable(label, false)
                (CImGui.IsItemHovered() || CImGui.IsItemFocused()) && CImGui.SetMouseCursor(i)
            end
            CImGui.TreePop()
        end
    end
end


text:
module Sheets

export Sheet

using MappedArrays

using ..Blobs
using ..Elements
using Base: copy!
import ..Motions: position, mutually_induce_velocity!, self_induce_velocity!, advect!

const MappedPositions{T,R} = MappedArrays.ReadonlyMappedArray{ComplexF64,1,Array{Blob{T,R},1},typeof(Elements.position)} where {T,R}

mutable struct Sheet{T,R} <: Element
    blobs::Vector{Blob{T,R}}
    Ss::Vector{T}
    δ::Float64
    zs::MappedPositions{T,R}
end

Elements.unwrap(s::Sheet) = s.blobs

function Sheet(zs::AbstractVector{Complex{T}}, Ss::AbstractVector{S}, δ::Float64) where {T <: Number, S <: Number}
    dSs = compute_trapezoidal_weights(Ss)
    blobs = Blob{S}.(zs, dSs, δ)

    zs = mappedarray(Elements.position, blobs)

    Sheet{S,T}(blobs, copy(Ss), δ, zs)
end

function Sheet(blobs::Vector{Blob{T,R}}, Ss::AbstractVector, δ::Float64) where {T,R}
    newblobs = copy(blobs)
    zs = mappedarray(Elements.position, blobs)
    Sheet{T,R}(newblobs, copy(Ss), δ, zs)
end

Base.length(s::Sheet) = length(s.blobs)

function mutually_induce_velocity!(ws₁, ws₂, sheet₁::Sheet, sheet₂::Sheet, t)
    mutually_induce_velocity!(ws₁, ws₂, sheet₁.blobs, sheet₂.blobs, t)
    nothing
end

function self_induce_velocity!(ws, sheet::Sheet, t)
    self_induce_velocity!(ws, sheet.blobs, t)
end

function compute_trapezoidal_weights(Ss)
    N = length(Ss)

    dSs = similar(Ss)
    dSs[1] = 0.5*(Ss[2] - Ss[1])
    for i in 2:N-1
        dSs[i] = 0.5*(Ss[i+1] - Ss[i-1])
    end
    dSs[N] = 0.5*(Ss[N] - Ss[N-1])

    return dSs
end

function advect!(sheet₊::Sheet{S,R}, sheet₋::Sheet{S,R}, ws, Δt) where {S,R}
    if sheet₊ != sheet₋
        N₊ = length(sheet₊)
        N₋ = length(sheet₋)
        if N₊ != N₋
            resize!(sheet₊.blobs, N₋)
            resize!(sheet₊.Ss, N₋)
        end
        copy!(sheet₊.Ss, sheet₋.Ss)
        sheet₊.δ = sheet₋.δ
    end

    advect!(sheet₊.blobs, sheet₋.blobs, ws, Δt)
    sheet₊
end

include("sheets/surgery.jl")

end

ext: jl


text:
using InspectDR
using Interact
using Blink
include("../notebook/BodePlots.jl") #Defines BodePlots.* utilities.

#When MIME"image/svg+xml" is disabled, Jupyter eventually requests PNG inline graphics:
#(Not certain if PNG is faster than SVG)
InspectDR.defaults.rendersvg = false

#Create base Observable-s that control plot output:
fmult = Observable(Float64(1e9)) #For all user-input
fmax = Observable(Float64(10)) #Maximum plot frequency
G = Observable(Float64(20)) #Gain (dB)
fz = Observable(Float64(fmax[])) #Zero frequency
fp1 = Observable(Float64(.05)) #Pole 1 frequency
fp2 = Observable(Float64(.9)) #Pole 2 frequency
annot = Observable(Bool(true)) #Annotate plot?
dispinline = Observable(Bool(true)) #Display inline?

#Define desired frequency ranges:
fmult_opt = OrderedDict("kHz" => 1.0e3, "MHz" => 1.0e6, "GHz" => 1.0e9) #Make sure to use Float64
fmult_lookup = OrderedDict([v=>k for (k,v) in fmult_opt]) #For labels

#-------------------------------------------------------------------------------
#(Workaround for updates being supressed on mutable objects)
#-------------------------------------------------------------------------------
#Immutable reference to object for throttle():
#(Workaround for updates being supressed on mutable objects)
struct ThrottleRef
	v
end
Base.:(==)(::ThrottleRef, ::ThrottleRef) = false #Always different (always update)
Base.show(io::IO, m::MIME, r::ThrottleRef) = show(io, m, r.v) #Redirect show to referenced object
Base.showable(m::MIME, r::ThrottleRef) = showable(m, r.v) #Redirect showable to referenced object
#-------------------------------------------------------------------------------

function frequency_slider(sig::Observable, fmax::Observable, fmult::Observable; label::String="")
	_fmin = 1e-3
	rng(_fmax) = 10 .^ range(log10(_fmin), stop=log10(max(_fmin, _fmax)), length=100)
	#NOTE: Limit max value to 2_fmin to avoid issues:
	return map(fmax, fmult) do _fmax, _fmult #Max freq & units on label depend on others
		units = fmult_lookup[_fmult]
		return slider(rng(_fmax), value=sig, label="$label ($units)")
	end
end

#Construct widgets used to control plots:
tgl_fmult = togglebuttons(fmult_opt, value=fmult, label="freq multiplier")
txt_fmax = spinbox(value=fmax, label="max freq")
sld_G = slider(-10:1.0:100, value=G, label="Gain (dB)")
sld_fz = frequency_slider(fz, fmax, fmult, label="zero")
sld_fp1 = frequency_slider(fp1, fmax, fmult, label="pole 1")
sld_fp2 = frequency_slider(fp2, fmax, fmult, label="pole 2")
chk_annot = checkbox(annot, label="Display annotation")
chk_inline = checkbox(dispinline, label="Display inline (throttled for speed)")
btn_creategtkplot = button("Create GTK plot")

#Build widget list to control plot for future code blocks
widgetlist = [tgl_fmult, txt_fmax, sld_G, sld_fz, sld_fp1, sld_fp2, chk_annot, chk_inline, btn_creategtkplot]

#Widget-controlled plot:
pobj = BodePlots.new() #Create initial plot object
updater_plot = map(fmult, fmax, G, fz, fp1, fp2, annot) do s, _fmax, _G, _fz, _fp1, _fp2, a
	BodePlots.update(pobj, s*1e-3, _fmax*s, _G, _fz*s, _fp1*s, _fp2*s, a)
end

#Conditionally display inline plot with throttling to reduce traffic:
throttled_plot = throttle(0.1, map(updater_plot, dispinline) do _plot, _inline
	if _inline
		return ThrottleRef(_plot)
	else
		return nothing
	end
end)

#Provide means to spawn new InspectDR-GTK plot windows:
creator_gtkplot = map(observe(btn_creategtkplot)) do _s
	#Display plot in Gtk GUI:
	gtkgui = display(InspectDR.GtkDisplay(), pobj) #Use same plot object

	#Refresh GUI when plot changes:
	updater_gtkgui = map(updater_plot) do _up
		InspectDR.refresh(gtkgui)
	end
end

#Display main window:
wnd = Window()
body!(wnd, tgl_fmult) #first call to body! does not always work.
body!(wnd, hbox(vbox(widgetlist...), vbox(throttled_plot)))

:DEMO_READY


text:
using Junior.Model

@record immutable Comment
    @int id
    @string content
    @reference author
end

# concurrency primitives

c = Chan(1)

@timeout 80

@go begin
    @get c
    @set c value
    @alt! begin
    end
end

@go begin
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

import LinearAlgebra: AbstractTriangular

"""
    SparseMatrixCSCSymmHerm

`Symmetric` or `Hermitian` of a `SparseMatrixCSC` or `SparseMatrixCSCView`.
"""
const SparseMatrixCSCSymmHerm{Tv,Ti} = Union{Symmetric{Tv,<:SparseMatrixCSCUnion{Tv,Ti}},
                                            Hermitian{Tv,<:SparseMatrixCSCUnion{Tv,Ti}}}

const AbstractTriangularSparse{Tv,Ti} = AbstractTriangular{Tv,<:SparseMatrixCSCUnion{Tv,Ti}}

# converting Symmetric/Hermitian/AbstractTriangular/SubArray of SparseMatrixCSC
# and Transpose/Adjoint of AbstractTriangular of SparseMatrixCSC to SparseMatrixCSC
for wr in (Symmetric, Hermitian, Transpose, Adjoint,
           UpperTriangular, LowerTriangular, UnitUpperTriangular, UnitLowerTriangular,
           SubArray)

    @eval SparseMatrixCSC(A::$wr) = _sparsem(A)
    @eval SparseMatrixCSC{Tv}(A::$wr{Tv}) where Tv = _sparsem(A)
    @eval SparseMatrixCSC{Tv}(A::$wr) where Tv = SparseMatrixCSC{Tv}(_sparsem(A))
    @eval SparseMatrixCSC{Tv,Ti}(A::$wr) where {Tv,Ti} = SparseMatrixCSC{Tv,Ti}(_sparsem(A))
end

"""
    iswrsparse(::S)
    iswrsparse(::Type{S})

Returns `true` if type `S` is backed by a sparse array, and `false` otherwise.
"""
iswrsparse(::T) where T<:AbstractArray = iswrsparse(T)
iswrsparse(::Type) = false
iswrsparse(::Type{T}) where T<:AbstractSparseArray = true

"""
    depth(::Type{S})

Returns 0 for unwrapped S, and nesting depth for wrapped (nested) abstract arrays.
"""
depth(::T) where T = depth(T)
depth(::Type{T}) where T<:AbstractArray = 0

for wr in (Symmetric, Hermitian,
           LowerTriangular, UnitLowerTriangular, UpperTriangular, UnitUpperTriangular,
           Transpose, Adjoint, SubArray,
           Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal)

    pl = wr === SubArray ? :($wr{<:Any,<:Any,T}) : :($wr{<:Any,T})
    @eval iswrsparse(::Type{<:$pl}) where T = iswrsparse(T)
    @eval depth(::Type{<:$pl}) where T = depth(T) + 1
end

# convert parent and re-wrap in same wrapper
_sparsewrap(A::Symmetric) = Symmetric(_sparsem(parent(A)), A.uplo == 'U' ? :U : :L)
_sparsewrap(A::Hermitian) = Hermitian(_sparsem(parent(A)), A.uplo == 'U' ? :U : :L)
_sparsewrap(A::SubArray) = SubArray(_sparsem(parent(A)), A.indices)
for ty in ( LowerTriangular, UnitLowerTriangular,
            UpperTriangular, UnitUpperTriangular,
            Transpose, Adjoint)

    @eval _sparsewrap(A::$ty) = $ty(_sparsem(parent(A)))
end
function _sparsewrap(A::Union{Diagonal,Bidiagonal,Tridiagonal,SymTridiagonal})
    dropzeros!(sparse(A))
end

"""
    unwrap(A::AbstractMatrix)

In case A is a wrapper type (`SubArray, Symmetric, Adjoint, SubArray, Triangular, Tridiagonal`, etc.)
convert to `Matrix` or `SparseMatrixCSC`, depending on final storage type of A.
For other types return A itself.
"""
unwrap(A::Any) = A
unwrap(A::AbstractMatrix) = iswrsparse(A) ? convert(SparseMatrixCSC, A) : convert(Array, A)

# For pure sparse matrices and vectors return A.
# For wrapped sparse matrices or vectors convert to SparseMatrixCSC.
# Handle nested wrappers properly.
# Use abstract matrix fallback if A is not sparse.
function _sparsem(@nospecialize A::AbstractArray{Tv}) where Tv
    if iswrsparse(A)
        if depth(A) >= 1
            _sparsem(_sparsewrap(A))
        else
            A
        end
    else
        # explicitly call abstract matrix fallback using getindex(A,...)
        invoke(SparseMatrixCSC{Tv,Int}, Tuple{AbstractMatrix}, A)
    end
end

_sparsem(A::AbstractSparseMatrix) = A
_sparsem(A::AbstractSparseVector) = A

# Transpose/Adjoint of sparse vector (returning sparse matrix)
function _sparsem(A::Union{Transpose{<:Any,<:AbstractSparseVector},Adjoint{<:Any,<:AbstractSparseVector}})
    B = parent(A)
    n = length(B)
    Ti = eltype(B.nzind)
    fadj = A isa Transpose ? transpose : adjoint
    colptr = fill!(Vector{Ti}(undef, n + 1), 0)
    colptr[1] = 1
    colptr[B.nzind .+ 1] .= 1
    cumsum!(colptr, colptr)
    rowval = fill!(similar(B.nzind), 1)
    nzval = fadj.(nonzeros(B))
    SparseMatrixCSC(1, n, colptr, rowval, nzval)
end

function _sparsem(A::Union{Transpose{<:Any,<:SparseMatrixCSC},Adjoint{<:Any,<:SparseMatrixCSC}})
    ftranspose(parent(A), A isa Transpose ? transpose : adjoint)
end

# Symmetric/Hermitian of sparse matrix
_sparsem(A::SparseMatrixCSCSymmHerm) = _sparsem(A.uplo == 'U' ? nzrangeup : nzrangelo, A)
# Triangular of sparse matrix
_sparsem(A::UpperTriangular{T,<:AbstractSparseMatrix}) where T = triu(A.data)
_sparsem(A::LowerTriangular{T,<:AbstractSparseMatrix}) where T = tril(A.data)
# view of sparse matrix
_sparsem(S::SubArray{<:Any,2,<:SparseMatrixCSC}) = getindex(S.parent,S.indices...)

# 4 cases: (Symmetric|Hermitian) variants (:U|:L)
function _sparsem(fnzrange::Function, sA::SparseMatrixCSCSymmHerm{Tv}) where {Tv}
    A = sA.data
    rowval = rowvals(A)
    nzval = nonzeros(A)
    m, n = size(A)
    Ti = eltype(rowval)
    fadj = sA isa Symmetric ? transpose : adjoint
    newcolptr = Vector{Ti}(undef, n+1)
    diagmap = fadj == transpose ? identity : real

    newcolptr[1] = 1
    colrange = fnzrange === nzrangeup ? (1:n) : (n:-1:1)
    @inbounds for j = colrange
        r = fnzrange(A, j); r1 = r.start; r2 = r.stop
        newcolptr[j+1] = r2 - r1 + 1
        for k = r1:r2
            row = rowval[k]
            if row != j
                newcolptr[row+1] += 1
            end
        end
    end
    cumsum!(newcolptr, newcolptr)
    nz = newcolptr[n+1] - 1
    newrowval = Vector{Ti}(undef, nz)
    newnzval = Vector{Tv}(undef, nz)
    @inbounds for j = 1:n
        newk = newcolptr[j]
        for k = fnzrange(A, j)
            i = rowval[k]
            nzv = nzval[k]
            if i != j
                newrowval[newk] = i
                newnzval[newk] = nzv
                newk += 1
                ni = newcolptr[i]
                newrowval[ni] = j
                newnzval[ni] = fadj(nzv)
                newcolptr[i] = ni + 1
            else
                newrowval[newk] = i
                newnzval[newk] = diagmap(nzv)
                newk += 1
            end
        end
        newcolptr[j] = newk
    end
    _sparse_gen(m, n, newcolptr, newrowval, newnzval)
end

# 2 cases: Unit(Upper|Lower)Triangular{Tv,SparseMatrixCSC}
function _sparsem(A::AbstractTriangularSparse{Tv}) where Tv
    S = A.data
    rowval = rowvals(S)
    nzval = nonzeros(S)
    m, n = size(S)
    Ti = eltype(rowval)
    fnzrange = A isa Union{UpperTriangular,UnitUpperTriangular} ? nzrangeup : nzrangelo
    unit = A isa Union{UnitUpperTriangular,UnitLowerTriangular}
    newcolptr = Vector{Ti}(undef, n+1)
    newrowval = Vector{Ti}(undef, nnz(S))
    newnzval = Vector{Tv}(undef, nnz(S))
    newcolptr[1] = 1
    uplo = fnzrange == nzrangeup
    newk = 1
    @inbounds for j = 1:n
        newkk = newk
        if unit
            newk += !uplo
        end
        r = fnzrange(S, j); r1 = r.start; r2 = r.stop
        for k = r1:r2
            i = rowval[k]
            if i != j || i == j && !unit
                newrowval[newk] = i
                newnzval[newk] = nzval[k]
                newk += 1
            end
        end
        if unit
            uplo && (newkk = newk)
            newrowval[newkk] = j
            newnzval[newkk] = one(Tv)
            newk += uplo
        end
        newcolptr[j+1] = newk
    end
    nz = newcolptr[n+1] - 1
    resize!(newrowval, nz)
    resize!(newnzval, nz)
    SparseMatrixCSC(m, n, newcolptr, newrowval, newnzval)
end

# 8 cases: (Transpose|Adjoint){Tv,[Unit](Upper|Lower)Triangular}
function _sparsem(taA::Union{Transpose{Tv,<:AbstractTriangularSparse},
                             Adjoint{Tv,<:AbstractTriangularSparse}}) where {Tv}

    sA = taA.parent
    A = sA.data
    rowval = rowvals(A)
    nzval = nonzeros(A)
    m, n = size(A)
    Ti = eltype(rowval)
    fnzrange = sA isa Union{UpperTriangular,UnitUpperTriangular} ? nzrangeup : nzrangelo
    fadj = taA isa Transpose ? transpose : adjoint
    unit = sA isa Union{UnitUpperTriangular,UnitLowerTriangular}
    uplo = A isa Union{UpperTriangular,UnitUpperTriangular}

    newcolptr = Vector{Ti}(undef, n+1)
    fill!(newcolptr, 1unit)
    newcolptr[1] = 1
    @inbounds for j = 1:n
        for k = fnzrange(A, j)
            i = rowval[k]
            if i != j || i == j && !unit
                newcolptr[i+1] += 1
            end
        end
    end
    cumsum!(newcolptr, newcolptr)
    nz = newcolptr[n+1] - 1
    newrowval = Vector{Ti}(undef, nz)
    newnzval = Vector{Tv}(undef, nz)

    @inbounds for j = 1:n
        if !uplo && unit
            ni = newcolptr[j]
            newrowval[ni] = j
            newnzval[ni] = fadj(one(Tv))
            newcolptr[j] = ni + 1
        end
        for k = fnzrange(A, j)
            i = rowval[k]
            nzv = nzval[k]
            if i != j || i == j && !unit
                ni = newcolptr[i]
                newrowval[ni] = j
                newnzval[ni] = fadj(nzv)
                newcolptr[i] = ni + 1
            end
        end
        if uplo && unit
            ni = newcolptr[j]
            newrowval[ni] = j
            newnzval[ni] = fadj(one(Tv))
            newcolptr[j] = ni + 1
        end
    end
    _sparse_gen(n, m, newcolptr, newrowval, newnzval)
end

function _sparse_gen(m, n, newcolptr, newrowval, newnzval)
    @inbounds for j = n:-1:1
        newcolptr[j+1] = newcolptr[j]
    end
    newcolptr[1] = 1
    SparseMatrixCSC(m, n, newcolptr, newrowval, newnzval)
end

max_stars_repo_name: ajozefiak/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# $$
# R_{\mu\nu} - \frac{1}{2}g_{\mu\nu}R = \frac{8\pi G}{c^4}T_{\mu\nu}
# $$

# %%
# https://nbviewer.jupyter.org/urls/gist.githubusercontent.com/terasakisatoshi/e4af4f618161d87cdc288285585527c5/raw/3a50de35748be12ab6461156440685052aadfcd2/hexagon_something.ipynb

using HTTP, JSON
function display_tweet(link)
    api = "https://publish.twitter.com/oembed?url=$link"
    r = response = HTTP.request("GET", api);
    j = JSON.parse(String(r.body))
    HTML(j["html"])
end
display_tweet("https://twitter.com/fujitapiroc1964/status/1436641682234114055")

# %%

max_stars_repo_path: 0020/mathjax example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0020/mathjax example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0020/mathjax example.jl

max_forks_repo_name: genkuroki/public


text:
"""This module defines classes of the package."""
# Author: Henri Gérard <hgerard.proy@gmail.com>
# License: MIT

abstract type RegressionModel end

# Define an object to
struct LinearRegression <: RegressionModel
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------

    """
    function LinearRegression()
        """
        Aggregator of projection

        Parameters
        ----------
        i::Int

        robustModel::RobustModel

        Returns
        -------

        """
        return new()
    end
end

struct LogisticRegression <: RegressionModel
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------

    """
    function LogisticRegression()
        """
        Aggregator of projection

        Parameters
        ----------
        i::Int

        robustModel::RobustModel

        Returns
        -------

        """
        return new()
    end
end

abstract type ParallelParam end

struct Sequential <: ParallelParam
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------

    """
    function Sequential()
        """
        Aggregator of projection

        Parameters
        ----------
        i::Int

        robustModel::RobustModel

        Returns
        -------

        """
        return new()
    end
end

struct Parallel <: ParallelParam
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------

    """
    function Parallel()
        """
        Aggregator of projection

        Parameters
        ----------
        i::Int

        robustModel::RobustModel

        Returns
        -------

        """
        return new()
    end
end

abstract type AlgoParams end

struct OptParams <: AlgoParams
    itmax::Int
    stability::Float64
    learning_rate::Float64
    verbosity::Int
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------

    """
    function OptParams(itmax, stability, learning_rate; verbosity = 0)
        """
        Aggregator of projection

        Parameters
        ----------
        i::Int

        robustModel::RobustModel

        Returns
        -------

        """
        return new(itmax, stability, learning_rate, verbosity)
    end
end

struct ProjParams <: AlgoParams
    ITER_MAX::Int
    precision::Float64
    sample::Int
    para_proj::ParallelParam
    para_inter::ParallelParam
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------

    """
    function ProjParams(ITER_MAX, precision, sample; para_proj=Sequential(), para_inter=Sequential())
        """
        Aggregator of projection

        Parameters
        ----------
        i::Int

        robustModel::RobustModel

        Returns
        -------

        """
        return new(ITER_MAX, precision, sample, para_proj, para_inter)
    end
end

abstract type AmbiguitySet end

struct DivergenceSet <: AmbiguitySet
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------

    """
    function DivergenceSet()
        """
        Aggregator of projection

        Parameters
        ----------
        i::Int

        robustModel::RobustModel

        Returns
        -------

        """
        return new()
    end
end

struct WassersteinSet <: AmbiguitySet
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------

    """
    function WassersteinSet()
        """
        Aggregator of projection

        Parameters
        ----------
        i::Int

        robustModel::RobustModel

        Returns
        -------

        """
        return new()
    end
end

abstract type GeneralConstraint end

struct PositiveConstraint <: GeneralConstraint
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------

    """
    function PositiveConstraint()
        """
        Aggregator of projection

        Parameters
        ----------
        i::Int

        robustModel::RobustModel

        Returns
        -------

        """
        return new()
    end
end

struct DROConstraint <: GeneralConstraint
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------

    """
    function DROConstraint()
        """
        Aggregator of projection

        Parameters
        ----------
        i::Int

        robustModel::RobustModel

        Returns
        -------

        """
        return new()
    end
end

struct EntropicConstraint <: GeneralConstraint
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------

    """
    function EntropicConstraint()
        """
        Aggregator of projection

        Parameters
        ----------
        i::Int

        robustModel::RobustModel

        Returns
        -------

        """
        return new()
    end
end

abstract type DivergenceConstraint<:GeneralConstraint end

struct KLConstraint <: DivergenceConstraint
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------

    """
    function KLConstraint()
        """
        Aggregator of projection

        Parameters
        ----------
        i::Int

        robustModel::RobustModel

        Returns
        -------

        """
        return new()
    end
end

struct RobustModel
    descent_direction::Array{Float64,1}
    I0::UnitRange{Int64}
    name::GeneralConstraint
    regressionModel::RegressionModel
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------

    """
    function RobustModel(N, nb_features, ϵ, ambiguity, regressionModel)
        """
        Aggregator of projection

        Parameters
        ----------
        i::Int

        robustModel::RobustModel

        Returns
        -------

        """
        if ambiguity == "KLdivergence"
            descent_direction = [ϵ; 1; zeros(nb_features); ones(N)/N]
            I0 = 1:N
            dim = N+1
            name = KLConstraint()
        elseif ambiguity == "wasserstein"
            descent_direction = [ϵ; zeros(nb_features); ones(N)/N]
            I0 = 1:N^2
            dim = N^2+1
            name = DROConstraint()
        elseif ambiguity == "entropic"
            descent_direction = [1; zeros(nb_features); ones(N)/N]
            I0 = 1:N
            dim = N+1
            name = EntropicConstraint()
        end
        return new(descent_direction, I0, name, regressionModel)
    end
end

max_stars_repo_name: Henri-Gerard/robox

max_issues_repo_name: Henri-Gerard/robox

max_forks_repo_name: Henri-Gerard/robox


text:
using HTTP
using Test
using Pluto
using Pluto: ServerSession, ClientSession, SessionActions
using Pluto.Configuration
using Pluto.Configuration: notebook_path_suggestion, from_flat_kwargs, _convert_to_flags
using Pluto.WorkspaceManager: poll

@testset "Configurations" begin

cd(Pluto.project_relative_path("test")) do
    @test notebook_path_suggestion() == joinpath(pwd(), "")
end

@testset "from_flat_kwargs" begin
    opt = from_flat_kwargs(; compile="min", launch_browser=false)
    @test opt.compiler.compile == "min"
    @test opt.server.launch_browser == false

    @test_throws MethodError from_flat_kwargs(; asdfasdf="test")

    structs_kwargs = let
        structs = [
            Pluto.Configuration.ServerOptions,
            Pluto.Configuration.SecurityOptions,
            Pluto.Configuration.EvaluationOptions,
            Pluto.Configuration.CompilerOptions
        ]
        sets = [collect(fieldnames(s)) for s in structs]
        vcat(sets...)::Vector{Symbol}
    end

    from_flat_kwargs_kwargs = let
        method = only(methods(Pluto.Configuration.from_flat_kwargs))
        syms = method.slot_syms
        names = split(syms, "\0")[2:end-1]
        Symbol.(names)::Vector{Symbol}
    end

    # Verify that all struct fields can be set via `from_flat_kwargs`.
    # Also verifies ordering to improve code readability.
    @test structs_kwargs == from_flat_kwargs_kwargs
end

@testset "flag conversion" begin
    @test _convert_to_flags(Configuration.CompilerOptions(threads="123")) ==
        ["--startup-file=no", "--history-file=no", "--threads=123"]

    @test _convert_to_flags(Configuration.CompilerOptions(threads=123)) ==
        ["--startup-file=no", "--history-file=no", "--threads=123"]

    @test _convert_to_flags(Configuration.CompilerOptions()) ⊇
        ["--startup-file=no", "--history-file=no"]

    @test _convert_to_flags(Configuration.CompilerOptions(compile="min")) ⊇
    ["--compile=min", "--startup-file=no", "--history-file=no"]
end

@testset "Authentication" begin
    port = 1238
    options = Pluto.Configuration.from_flat_kwargs(; port, launch_browser=false, workspace_use_distributed=false)
    🍭 = Pluto.ServerSession(; options)
    fakeclient = ClientSession(:fake, nothing)
    🍭.connected_clients[fakeclient.id] = fakeclient
    host = 🍭.options.server.host
    secret = 🍭.secret
    println("Launching test server...")
    server_task = @async Pluto.run(🍭)
    sleep(2)

    local_url(suffix) = "http://$host:$port/$suffix"
    withsecret(url) = occursin('?', url) ? "$url&secret=$secret" : "$url?secret=$secret"
    @test HTTP.get(local_url("favicon.ico")).status == 200

    function requeststatus(url, method)
        r = HTTP.request(method, url; status_exception=false, redirect=false)
        r.status
    end

    nb = SessionActions.open(🍭, Pluto.project_relative_path("sample", "Basic.jl"); as_sample=true)

    simple_routes = [
        ("", "GET"),
        ("edit?id=$(nb.notebook_id)", "GET"),
        ("notebookfile?id=$(nb.notebook_id)", "GET"),
        ("notebookexport?id=$(nb.notebook_id)", "GET"),
        ("statefile?id=$(nb.notebook_id)", "GET"),
    ]

    function tempcopy(x)
        p = tempname()
        Pluto.readwrite(x, p)
        p
    end
    @assert isfile(Pluto.project_relative_path("sample", "Basic.jl"))

    effect_routes = [
        ("new", "GET"),
        ("new", "POST"),
        ("open?url=$(HTTP.URIs.escapeuri("https://raw.githubusercontent.com/fonsp/Pluto.jl/v0.14.5/sample/Basic.jl"))", "GET"),
        ("open?url=$(HTTP.URIs.escapeuri("https://raw.githubusercontent.com/fonsp/Pluto.jl/v0.14.5/sample/Basic.jl"))", "POST"),
        ("open?path=$(HTTP.URIs.escapeuri(Pluto.project_relative_path("sample", "Basic.jl") |> tempcopy))", "GET"),
        ("open?path=$(HTTP.URIs.escapeuri(Pluto.project_relative_path("sample", "Basic.jl") |> tempcopy))", "POST"),
        ("sample/Basic.jl", "GET"),
        ("sample/Basic.jl", "POST"),
        ("notebookupload", "POST"),
    ]

    for (suffix, method) in simple_routes ∪ effect_routes
        url = local_url(suffix)
        @test requeststatus(url, method) == 403
    end

    # no notebooks were opened
    @test length(🍭.notebooks) == 1

    for (suffix, method) in simple_routes
        url = local_url(suffix) |> withsecret
        @test requeststatus(url, method)  ∈ 200:299
    end

    for (suffix, method) in setdiff(effect_routes, [("notebookupload", "POST")])
        url = local_url(suffix) |> withsecret
        @test requeststatus(url, method) ∈ 200:399 # 3xx are redirects
    end

    @async schedule(server_task, InterruptException(); error=true)
end

@testset "Open Notebooks at Startup" begin
    port = 1338
    host = "localhost"
    local_url(suffix) = "http://$host:$port/$suffix"

    urls = [
    "https://raw.githubusercontent.com/fonsp/Pluto.jl/v0.12.16/sample/Basic.jl",
    "https://gist.githubusercontent.com/fonsp/4e164a262a60fc4bdd638e124e629d64/raw/8ffe93c680e539056068456a62dea7bf6b8eb622/basic_pkg_notebook.jl",
    ]
    nbnames = download.(urls)

    server_running() = HTTP.get(local_url("favicon.ico")).status == 200 && HTTP.get(local_url("edit")).status == 200

    # without notebook at startup
    server_task = @async Pluto.run(port=port, launch_browser=false, workspace_use_distributed=false, require_secret_for_access=false, require_secret_for_open_links=false)
    @test poll(5) do
        server_running()
    end
    @async schedule(server_task, InterruptException(); error=true)

    # with a single notebook at startup
    server_task = @async Pluto.run(notebook=first(nbnames), port=port, launch_browser=false, workspace_use_distributed=false, require_secret_for_access=false, require_secret_for_open_links=false)
    @test poll(5) do
        server_running()
    end
    @async schedule(server_task, InterruptException(); error=true)

    # with multiple notebooks at startup
    server_task = @async Pluto.run(notebook=nbnames, port=port, launch_browser=false, workspace_use_distributed=false, require_secret_for_access=false, require_secret_for_open_links=false)
    @test poll(5) do
        server_running()
    end
    @async schedule(server_task, InterruptException(); error=true)

end

# TODO are the processes closed properly?
# TODO we reuse the same port without awaiting the shutdown of the previous server

end # testset


text:
export Material, UniformMaterial, TexturedMaterial, AATexturedMaterial
export getcolor

abstract type Material end

"""
    struct UniformMaterial <: Material

A uniform-looking material.

# Initialization

    UniformMaterial(r, g, b)

Initialize the material with a color whose red, green, and blue components
are `r`, `g`, `b`. All the three values must be in the range ``[0, 1]``.

# Examples

```julia
texture = SolidState.UniformMaterial(0.5, 0.5, 0.1)
```
"""
struct UniformMaterial <: Material
    color::RGB{Float64}

    UniformMaterial(r, g, b) = new(RGB{Float64}(r, g, b))
end

getcolor(mat::UniformMaterial, u, v) = mat.color

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

"""
    struct TexturedMaterial <: Material

A textured material, mapping a bitmap to a ``(u, v)`` grid in the domain
``[0, 1] \times [0, 1]``.

If you are looking for an anti-aliased textured material, use
`AATexturedMaterial`: it produce good-looking results even for
low-resolution bitmaps, at the expense of more calculations.

# Initialization

    TexturedMaterial(img::Array{RGB{Float64}, 2})

Initialize the texture with a bitmap.

    TexturedMaterial(fname::AbstractString)

Initialize the texture with a bitmap loaded from a file with name
`fname`. All the image types implemented in the
[Images.jl](https://github.com/JuliaImages/Images.jl) package are
recognized.

# Examples

```julia
texture = SolidState.TexturedMaterial("picture.jpg")
```
"""
struct TexturedMaterial <: Material
    texture::Array{RGB{Float64}, 2}

    TexturedMaterial(img::Array{RGB{Float64}, 2}) = new(img)
    TexturedMaterial(fname::AbstractString) = new(convert.(RGB{Float64},
                                                           Images.load(fname)))
end

function getcolor(mat::TexturedMaterial, u, v)
    height, width = size(mat.texture)
    mat.texture[round(Int, height - v * (height - 1)),
                round(Int, u * (width - 1) + 1)]
end

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

interp_texture(img) = LinearInterpolation((1:size(img)[1], 1:size(img)[2]), img)

"""
    struct AATexturedMaterial <: Material

An anti-aliased textured material, mapping a bitmap to a ``(u, v)`` grid in
the domain ``[0, 1] \times [0, 1]``.

Use `TexturedMaterial` if your textures are high-resolution, as it is faster.
Anti-aliasing only helps when the resolution of textures are so low that the
rendered object looks pixelated.

# Initialization

    AATexturedMaterial(img::Array{RGB{Float64}, 2})

Initialize the texture with a bitmap.

    AATexturedMaterial(fname::AbstractString)

Initialize the texture with a bitmap loaded from a file with name
`fname`. All the image types implemented in the
[Images.jl](https://github.com/JuliaImages/Images.jl) package are
recognized.

# Examples

```julia
texture = SolidState.AATexturedMaterial("picture.jpg")
```
"""
struct AATexturedMaterial <: Material
    texture

    AATexturedMaterial(img::Array{RGB{Float64}, 2}) = new(interp_texture(img))
    AATexturedMaterial(fname::AbstractString) = new(interp_texture(convert.(RGB{Float64},
                                                                            Images.load(fname))))
end

function getcolor(mat::AATexturedMaterial, u, v)
    height, width = size(mat.texture)
    mat.texture(height - v * (height - 1), u * (width - 1) + 1)
end


text:
# int variable
x = 99
x, typeof(x)

# float variable
y = 10.0
y, typeof(y)

s = "Data Science"
s, typeof(s)

# string variable
t = true
t, typeof(t)

# assignment operator
p = 2
q = 3
r = 4
p, q, r

# division operator
r / p

# multiplication operator
p * q

# addition operator
p + q

# substration operator
r - q

q - r

# power operator
q ^ p

# equal comparison
p == p

p == q

# not equal comparison
p != q

# greater than or equal comparison
p >= q

# greater than comparison
p > q

# less than or equal comparison
p <= q

# less than comparison
p < q

# int 16
int_16 = Int16(16)
int_16, typeof(int_16), typemax(Int16)

# int 32
int_32 = Int32(32)
int_32, typeof(int_32), typemax(Int32)

# int 64
int_64 = Int64(64)
int_64, typeof(int_64), typemax(Int64)

# int 128
int_128 = Int64(128)
int_128, typeof(int_128), typemax(Int128)

# float 16
float_16 = Float16(16.0)
float_16, typeof(float_16), typemax(Float16)

# float 32
float_32 = Float32(32.0)
float_32, typeof(float_32), typemax(Float32)

# float 64
float_64 = Float64(64.0)
float_64, typeof(float_64), typemax(Float64)

# check if is +ve infinity
isinf(1/0)

# check if is -ve infinity
isinf(-1/0)

# check if is not a number
isnan(0/0)

s1 = "Data Science"

s2 = "I'm learning Data Science"

# get index 1
s2[1]

# get last index
s2[end]

# get 2 up to 6
s2[2:6]

# get last 6
s2[end-6:end]

# concat
s2 * ". It's fun!"

# string interpolation
"1 + 2 = $(1 + 2)"

"$x + $y = $(x + y)"

length(s2)

# find first occurance of a pattern
findfirst("Data", s2)

s3 = "pqrs"
s4 = "pqrt"
s3, s4

# check if string are equal
s3 == s4

s3 > s4

s3 < s4

# tuple data
t1 = (1, 2, 3, 4)
t1, typeof(t1)

t2 = ("Data", 1, 2.0)
t2, typeof(t2)

u, v, w = t2
u

# vector data
a1 = Array{Int64}(undef, 3) # vector of length 3
a1, typeof(a1)

a2 = Array{Int64}(undef, (3,2)) # 3 x 2 matrix
a2, typeof(a2)

a3 = [1,2,3]
a3

a4 = [1 2 3]

# matrix data
a5 = [1 2 3 ; 4 5 6]

a6 = [1 2 3 ; 4 5 6 ; 7 8 9]

zeros((2,2))

zeros(Float64, (2,2))

ones((2,2))

ones(Float64, (2,2))

rand(Int64, (2,2))

randn((2,2))

# dictionary data
d0 = Dict("id" => 100, "age" => 24)

d0["id"]

d0["age"]

d0["level"] = 1
d0

haskey(d0, "id")

delete!(d0, "level")
d0

# compaund expression
c0 = begin
    j = 10
    k = 2
    j * k
end
c0

c1 = (j = 10; k = 2; j * k)
c1

# if check
if k == 2
    "two"
end

# if..else check
if k == 2
    "two"
else
    "not two"
end

k == 2 ? "two" : "not two"

# for loop
s5 = "Julia"

for str in s5
    println(str)
end

# nested for loop
odd = [1, 3, 5]
even = [2, 4, 6]

for i in odd, j in even
    println("i * j = $(i * j)")
end

# while loop
first = 0
last = 10

while first < last
    global first
    println(first)
    first = first + 1
end

# catch exception
try
    log(-1)
catch ex
    println(ex)
end

# declare function
function add(x, y)
    return x+y
end

add(1, 2)

add2(x,y) = x+y

add2(1,2)

max_stars_repo_name: lykmapipo/data-science-learning

max_issues_repo_name: lykmapipo/data-science-learning

max_forks_repo_name: lykmapipo/data-science-learning


text:
using Chain
using Dates
using JSON3
using Base: @kwdef

using Stonks: JSONContent, ContentParserError
using Stonks.Models: AbstractStonksRecord

function snake_case(s::Union{String,Symbol})
  string = replace(String(s), r"[\-\.\s]" => "_")
  words = lowercase.(split(string, r"(?=[A-Z])"))
  res = join([i == 1 ? word : "_$word" for (i, word) in enumerate(words)])
  return isa(s, Symbol) ? Symbol(res) : res
end

function js_to_dict(js::JSON3.Object; key::Symbol=:raw, to_snake_case::Bool=false)
  vals = Dict{Symbol,Any}()
  for (k_raw, v) in js
    k = to_snake_case ? snake_case(k_raw) : k_raw
    if isa(v, JSON3.Object)
      vals[k] = get(v, Symbol(key), missing)
    else
      vals[k] = v
    end
  end
  return vals
end

function apply_filters(
  data::Vector{T},
  datecol::String;
  from::Union{Date,Missing}=missing,
  to::Union{Date,Missing}=missing,
) where {T<:AbstractStonksRecord}
  return @chain data begin
    isa(from, Date) ? filter(x -> getfield(x, Symbol(datecol)) >= from, _) : _
    isa(to, Date) ? filter(x -> getfield(x, Symbol(datecol)) <= to, _) : _
  end
end

function get_remap(remaps::Union{Dict{Symbol,Symbol},Missing}, key::Symbol)
  if !ismissing(remaps)
    haskey(remaps, key) ? remaps[key] : key
  else
    key
  end
end

function parse_jsvalue(js::AbstractDict, key::Symbol, T::Any)
  if haskey(js, key)
    val = js[key]
    value = (
      if typeof(val) === T
        val
      elseif T <: AbstractFloat && typeof(val) <: AbstractFloat
        T(val)
      elseif T === Date && isa(val, Date)
        val
      elseif T == String && isa(val, Date)
        val
      else
        tryparse(T, val)
      end
    )
    return isnothing(value) ? missing : value
  end
  return missing
end

union_types(x::Union) = (x.a, union_types(x.b)...)
union_types(x::Type) = (x,)

"""
    tryparse_js(T<:AbstractStonksRecord, js; [remaps], [fixed]) -> T<:AbstractStonksRecord

Tries to convert a flat JSON3 object to a data model type by builiding it using keyword arguments.

### Arguments
- `T<:AbstractStonksRecord`: a subtype of AbstractStonksRecord in which the js object will be parsed
- `js::JSON3.Object`: needs to be a flat object

### Keywords
- `[remaps]::Dict{Symbol,Symbol}`: entries T field => js field name. Use it only for keys which don't have the same name in the js as the type.
- `[fixed]::Dict{Symbol,String}`: entries like :key => value for fields with known values, like :symbol, :id, etc.
"""
function tryparse_js(
  ::Type{T},
  js::AbstractDict;
  remaps::Union{Dict{Symbol,Symbol},Missing}=missing,
  fixed::Union{Dict{Symbol,String},Missing}=missing,
) where {T<:AbstractStonksRecord}
  fields = @chain zip(fieldnames(T), T.types) begin
    [(field, collect(union_types(type))) for (field, type) in _]
    map(x -> (from=get_remap(remaps, x[1]), to=x[1], type=x[2]), _)
  end
  fields_req = @chain fields begin
    filter(ts -> all(map(t -> t !== Missing && t !== Nothing, ts.type)), _)
    map(first, _)
    !ismissing(remaps) ? map(from -> haskey(remaps, from) ? remaps[from] : from, _) : _
    !ismissing(fixed) ? vcat(_, keys(fixed)...) : _
  end
  js_fields = @chain [k for (k, v) in js] !ismissing(fixed) ? vcat(_, keys(fixed)...) : _
  diff = setdiff(fields_req, js_fields)
  if !isempty(diff)
    return ContentParserError("Required keys not found: $(join(diff, ", "))")
  end
  kwargs = @chain fields begin
    map(
      ts -> (
        from=ts.from,
        to=ts.to,
        type=first(filter(t -> t !== Missing && t !== Nothing, ts.type)),
      ),
      _,
    )
    Dict([k_to => parse_jsvalue(js, k_from, type) for (k_from, k_to, type) in _])
  end
  if !ismissing(fixed)
    [kwargs[k] = v for (k, v) in fixed]
  end
  return T(; kwargs...)
end


text:
abstract type AbstractDensityStratification{L,HV,TT,α,dρdz,g,Gdirec} end

    hasdensityles(a::Type{T}) where {L,HV,TT,α,dρdz,g,Gdirec,T<:AbstractDensityStratification{L,HV,TT,α,dρdz,g,Gdirec}} =
        L === NoLESScalar ? false : true
    @inline hasdensityles(a::AbstractDensityStratification{L,HV,TT,α,dρdz,g,Gdirec}) where {L,HV,TT,α,dρdz,g,Gdirec} = hasdensityles(typeof(a))

    has_les_density_vorticity_model(a::Type{T}) where {L,HV,TT,α,dρdz,g,Gdirec,T<:AbstractDensityStratification{L,HV,TT,α,dρdz,g,Gdirec}} =
        is_vorticity_model(L)
    @inline has_les_density_vorticity_model(a::AbstractDensityStratification{L,HV,TT,α,dρdz,g,Gdirec}) where {L,HV,TT,α,dρdz,g,Gdirec} = has_les_density_vorticity_model(typeof(a))
    @inline @par has_les_density_vorticity_model(s::Type{T}) where {T<:@par(AbstractSimulation)} = has_les_density_vorticity_model(DensityStratificationType)

    diffusivity(a::Type{T}) where {L,HV,TT,α,dρdz,g,Gdirec,T<:AbstractDensityStratification{L,HV,TT,α,dρdz,g,Gdirec}} =
        α
    @inline diffusivity(a::AbstractDensityStratification{L,HV,TT,α,dρdz,g,Gdirec}) where {L,HV,TT,α,dρdz,g,Gdirec} = diffusivity(typeof(a))
    @inline @par diffusivity(s::Type{T}) where {T<:@par(AbstractSimulation)} = diffusivity(DensityStratificationType)
    @inline diffusivity(s::AbstractSimulation) = diffusivity(typeof(s))

    meangradient(a::Type{T}) where {L,HV,TT,α,dρdz,g,Gdirec,T<:AbstractDensityStratification{L,HV,TT,α,dρdz,g,Gdirec}} =
        dρdz
    @inline meangradient(a::AbstractDensityStratification{L,HV,TT,α,dρdz,g,Gdirec}) where {L,HV,TT,α,dρdz,g,Gdirec} = meangradient(typeof(a))

    gravity(a::Type{T}) where {L,HV,TT,α,dρdz,g,Gdirec,T<:AbstractDensityStratification{L,HV,TT,α,dρdz,g,Gdirec}} =
        g
    @inline gravity(a::AbstractDensityStratification{L,HV,TT,α,dρdz,g,Gdirec}) where {L,HV,TT,α,dρdz,g,Gdirec} = gravity(typeof(a))
    @inline @par gravity(s::Type{T}) where {T<:@par(AbstractSimulation)} = gravity(DensityStratificationType)
    @inline gravity(s::AbstractSimulation) = gravity(typeof(s))

    graddir(a::Type{T}) where {L,HV,TT,α,dρdz,g,Gdirec,T<:AbstractDensityStratification{L,HV,TT,α,dρdz,g,Gdirec}} =
        Gdirec
    @inline graddir(a::AbstractDensityStratification{L,HV,TT,α,dρdz,g,Gdirec}) where {L,HV,TT,α,dρdz,g,Gdirec} = graddir(typeof(a))

    calculate_N(a::Type{T}) where {L,HV,TT,α,dρdz,g,Gdirec,T<:AbstractDensityStratification{L,HV,TT,α,dρdz,g,Gdirec}} =
        sqrt(g.z*dρdz)
    @inline calculate_N(a::AbstractDensityStratification{L,HV,TT,α,dρdz,g,Gdirec}) where {L,HV,TT,α,dρdz,g,Gdirec} = calculate_N(typeof(a))

    initialize!(a::AbstractDensityStratification,s::AbstractSimulation) = initialize!(a.timestep,parent(real(a.rhs)),diffusivity(a),s)

    function statsheader(a::AbstractDensityStratification{L,HV}) where {L,HV}
        mles = L !== NoLESScalar ? ",rhopr" : ""
        mhvis = HV !== NoHyperViscosity ? ",rhohdiss" : ""
        return "rho,rhop2,drhodxp2,drhodyp2,drhodzp2,rhodiss,trhodiss,bflux"*mles*mhvis
    end

    function stats(a::AbstractDensityStratification{L,HV},s::A) where {L,HV,A<:AbstractSimulation}
        sstats, sstats_fil = scalar_stats(a.ρ,a,s)
        trhodiss = sstats[end]
        trhodiss_fil = sstats_fil[end]
        bf, bf_fil = buoyancy_flux(a,s)

        if L !== NoLESScalar
            if !STAT_FIL
                lesstats = (scalar_les_stats(a.reduction,a.flux,a.ρ),)
                lesstats_fil = (0.,)
                trhodiss += lesstats[1]
            else
                lesstatsv, lesstats_filv = scalar_les_stats(a.reduction,s.reductionh,a.flux,a.ρ)
                lesstats = (lesstatsv,)
                lesstats_fil = (lesstats_filv,)
                trhodiss += lesstats[1]
                trhodiss_fil += lesstats_fil[1]
            end
        else
            lesstats = ()
            lesstats_fil = ()
        end

        if HV !== NoHyperViscosity
            if !STAT_FIL
                hvstats = (scalar_hvis_stats(a.reduction,a.ρ,a.hyperviscosity),)
                hvsstats_fil = (0.,)
                trhodiss += hvstats[1]
            else
                hvstatsv, hvstats_filv = scalar_hvis_stats(a.reduction,s.reductionh,a.ρ,a.hyperviscosity)
                hvstats = (hvstatsv,)
                hvsstats_fil = (hvstats_filv,)
                trhodiss += hvstats[1]
                trhodiss_fil += hvstats_fil[1]
            end
        else
            hvstats = ()
            hvstats_fil = ()
        end

        return (sstats...,trhodiss,bf,lesstats...,hvstats...),(sstats_fil...,trhodiss_fil,bf_fil,lesstats_fil...,hvstats_fil...)
    end

struct NoDensityStratification <: AbstractDensityStratification{NoLESScalar,NoHyperViscosity,nothing,nothing,nothing,nothing,nothing} end

    initialize!(a::NoDensityStratification,s::AbstractSimulation) = nothing

    statsheader(a::NoDensityStratification) = ""

    stats(a::NoDensityStratification,s::AbstractSimulation) = (),()

    msg(a::NoDensityStratification) = "\nDensity Stratification: No density stratification\n"

struct BoussinesqApproximation{L,HV,TTimeStep, α #=Difusitivity = ν/Pr =#,
                   dρdz #=Linear mean profile=#, g #=This is actually g/ρ₀ =#,
                   Gdirec#=Gravity direction =#} <: AbstractDensityStratification{L,HV,TTimeStep,α,dρdz,g,Gdirec}
    ρ::ScalarField{Float64,3,2,false}
    rhs::ScalarField{Float64,3,2,false}
    flux::VectorField{Float64,3,2,false}
    lesmodel::L
    hyperviscosity::HV
    timestep::TTimeStep
    reduction::Vector{Float64}

    function BoussinesqApproximation{TT,α,dρdz,g,Gdirec}(ρ,timestep,les,hv) where {TT,α,dρdz,g,Gdirec}
        ρrhs = similar(ρ)
        flux = VectorField(size(ρ.field.r),(ρ.kx.l,ρ.ky.l,ρ.kz.l))
        reduction = zeros(THR ? Threads.nthreads() : 1)
        return new{typeof(les),typeof(hv),TT,α,dρdz,g,Gdirec}(ρ,ρrhs,flux,les,hv,timestep,reduction)
    end
end

initialize!(a::BoussinesqApproximation,s::AbstractSimulation) = initialize!(a.timestep,parent(real(a.rhs)),diffusivity(a),s)

msg(a::BoussinesqApproximation{L,HV,TT,α,dρdz,g,Gdirec}) where {L,HV,TT,α,dρdz,g,Gdirec} = """

Density Stratification: Boussinesq Approximation
Density diffusivity: $(α)
g/ρ₀ : $(g)
Density mean gradient : $(dρdz)
Density mean gradient direction: $(Gdirec)
Density time-stepping method: $(TT)
Density LES model: $(msg(a.lesmodel))
Density HyperViscosity model: {$(msg(a.hyperviscosity))\n}

"""

function stratificationtype(d::AbstractDict,start,integrator,nx,ny,nz,lx,ly,lz,variableTimeStep,idt)

    s = NoDensityStratification()

    if haskey(d,:densityHyperViscosity)
        if d[:densityHyperViscosity] == "spectralBarrier"
            initkp = haskey(d,:densityInitk) ? parse(Float64,d[:densityInitk]) : Globals.cutoffr*2*KX[end] / 3
            endK = haskey(d,:densityEndk) ? parse(Float64,d[:densityEndk]) : sqrt(Globals.cutoffr*3)*KX[end]
            hyperviscositytype = SpectralBarrier(initkp,endK)
        else
            νh = parse(Float64,d[:densityHyperViscosityCoefficient])
            m = haskey(d,:densityHyperViscosityM) ? parse(Int,d[:densityHyperViscosityM]) : 2
            hyperviscositytype = HyperViscosity{νh,m}()
        end
    else
        hyperviscositytype = NoHyperViscosity()
    end

    if haskey(d,:densityStratification)

        haskey(d,:gravityDirection) ? (gdir = Symbol(d[:gravityDirection])) : (gdir = :z)
        α = ν/Float64(eval(Meta.parse(d[:Pr])))
        dρdz = Float64(eval(Meta.parse(d[:densityGradient])))
        gval = Float64(eval(Meta.parse(d[:zAcceleration])))/Float64(eval(Meta.parse(d[:referenceDensity])))
        if gdir === :z
            g = Vec(0.0,0.0,gval)
        elseif gdir === :y
            g = Vec(0.0,gval,0.0)
        elseif gdir === :x
            g = Vec(gval,0.0,0.0)
        end

        @info("Reading initial density field rho.$start")
        rho = isfile("rho.$start") ? ScalarField{Float64}("rho.$start") : ScalarField{Float64}((nx,ny,nz),(lx,ly,lz))
        gdir = haskey(d,:gravityDirection) ? Symbol(d[:gravityDirection]) : :z
        densitytimestep = if integrator === :Euller
            Euller{variableTimeStep,idt}(Ref(idt))
        elseif integrator === :Adams_Bashforth3rdO
            Adams_Bashforth3rdO{variableTimeStep,idt}()
        else
            ETD3rdO(variableTimeStep,idt,hyperviscositytype,α)
        end

        lestypedensity = if haskey(d,:densityLesModel)
            if d[:densityLesModel] == "vorticityDiffusion"
                c = haskey(d,:densityVorticityCoefficient) ? parse(Float64,d[:densityVorticityCoefficient]) : 0.1
                VorticityDiffusion(c)
            else
                EddyDiffusion()
            end
        else
            NoLESScalar()
        end
        s = BoussinesqApproximation{typeof(densitytimestep),α,dρdz,g,gdir}(rho,densitytimestep,lestypedensity,hyperviscositytype)
    end

    return s
end

function buoyancy_flux(a::AbstractDensityStratification,s::AbstractSimulation)
    gdir = graddir(a)
    g = gravity(s)
    if !STAT_FIL
        if gdir === :z
           return -g.z*proj_mean(a.reduction,s.u.c.z,a.ρ), 0.0
        elseif gdir === :y
           return -g.y*proj_mean(a.reduction,s.u.c.y,a.ρ), 0.0
        elseif gdir === :x
           return -g.x*proj_mean(a.reduction,s.u.c.x,a.ρ), 0.0
        end
    else
        if gdir === :z
           return -g.z .* proj_mean(a.reduction,s.reductionh,s.u.c.z,a.ρ)
        elseif gdir === :y
           return -g.y .* proj_mean(a.reduction,s.reductionh,s.u.c.y,a.ρ)
        elseif gdir === :x
           return -g.x .* proj_mean(a.reduction,s.reductionh,s.u.c.x,a.ρ)
        end
    end
end


text:
@testset "empty arrays" begin
    A = FixedValueArray(Nothing, 0)
    @test A isa FixedValueArray{DataType, 1}
    @test size(A) == (0,)
    @test length(A) == 0
    @test_throws BoundsError A[]
    @test_throws BoundsError A[1]

    A = FixedValueArray(Nothing, 0, 0)
    @test A isa FixedValueArray{DataType, 2}
    @test size(A) == (0, 0)
    @test length(A) == 0
    @test_throws BoundsError A[]
    @test_throws BoundsError A[1]

    A = FixedValueArray(Nothing, 0, 1)
    @test A isa FixedValueArray{DataType, 2}
    @test size(A) == (0, 1)
    @test length(A) == 0
    @test_throws BoundsError A[]
    @test_throws BoundsError A[1]

    A = FixedValueArray(Nothing, 1, 0)
    @test A isa FixedValueArray{DataType, 2}
    @test size(A) == (1, 0)
    @test length(A) == 0
    @test_throws BoundsError A[]
    @test_throws BoundsError A[1]

    A = FixedValueArray(Nothing, 0, 0, 0)
    @test A isa FixedValueArray{DataType, 3}
    @test size(A) == (0, 0, 0)
    @test length(A) == 0
    @test_throws BoundsError A[]
    @test_throws BoundsError A[1]
end

@testset "0-dimensional arrays" begin
    A = FixedValueArray(Nothing)
    @test A isa FixedValueArray{DataType, 0}
    @test size(A) == ()
    @test length(A) == 1
    @test_throws BoundsError A[0]
    @test A[] == Nothing
    @test_throws CanonicalIndexError A[] = Nothing
    @test A[1] == Nothing
    @test_throws CanonicalIndexError A[1] = Nothing
    @test_throws BoundsError A[2]
end

@testset "1-dimensional arrays" begin
    A = FixedValueArray(42, 3)
    @test A isa FixedValueArray{Int, 1}
    @test size(A) == (3,)
    @test length(A) == 3
    @test_throws BoundsError A[]
    @test_throws CanonicalIndexError A[] = 41
    @test_throws BoundsError A[0]
    @test_throws CanonicalIndexError A[0] = 41
    for e ∈ A
        @test e == 42
    end
    for i ∈ 1:length(A)
        @test A[i] == 42
        @test_throws CanonicalIndexError A[i] = 41
    end
    @test_throws BoundsError A[4]
    @test_throws CanonicalIndexError A[4] = 41
end

@testset "2-dimensional arrays" begin
    A = FixedValueArray(42, 2, 3)
    @test A isa FixedValueArray{Int, 2}
    @test size(A) == (2, 3)
    @test length(A) == 6
    @test_throws BoundsError A[]
    @test_throws CanonicalIndexError A[] = 41
    @test_throws BoundsError A[0]
    @test_throws CanonicalIndexError A[0] = 41
    @test_throws BoundsError A[0, 1]
    @test_throws CanonicalIndexError A[0, 1] = 41
    @test_throws BoundsError A[1, 0]
    @test_throws CanonicalIndexError A[1, 0] = 41
    for e ∈ A
        @test e == 42
    end
    for i ∈ 1:length(A)
        @test A[i] == 42
        @test_throws CanonicalIndexError A[i] = 41
    end
    for i ∈ 1:size(A)[1], j ∈ 1:size(A)[2]
        @test A[i, j] == 42
        @test_throws CanonicalIndexError A[i, j] = 42
    end
    @test_throws BoundsError A[7]
    @test_throws CanonicalIndexError A[7] = 41
    @test_throws BoundsError A[3, 1]
    @test_throws CanonicalIndexError A[0, 1] = 41
    @test_throws BoundsError A[1, 4]
    @test_throws CanonicalIndexError A[1, 0] = 41
end


text:
# Method type
@enum MethodType begin
  DIRK
  COEFFICIENT_MULTISTEP
  NORDSIECK_MULTISTEP
  GLM
end

# solver

abstract type AbstractNLSolver{algType,iip} end

mutable struct NLSolver{algType,iip,uType,gamType,tmpType,tType,C<:AbstractNLSolverCache} <: AbstractNLSolver{algType,iip}
  z::uType
  tmp::uType # DIRK and multistep methods only use tmp
  tmp2::tmpType # for GLM if neccssary
  ztmp::uType
  γ::gamType
  c::tType
  α::tType
  alg::algType
  κ::tType
  fast_convergence_cutoff::tType
  ηold::tType
  iter::Int
  maxiters::Int
  status::NLStatus
  cache::C
  method::MethodType
end

# default to DIRK
function NLSolver{iip,tType}(z, tmp, ztmp, γ, c, α, alg, κ, fast_convergence_cutoff, ηold, iter, maxiters, status, cache, method=DIRK, tmp2=nothing) where {iip,tType}
  NLSolver{typeof(alg), iip, typeof(z), typeof(γ), typeof(tmp2), tType, typeof(cache)}(
    z, tmp, tmp2, ztmp, γ,
    convert(tType, c), convert(tType, α), alg,
    convert(tType, κ), convert(tType, fast_convergence_cutoff),
    convert(tType, ηold), iter, maxiters, status, cache, method
  )
end

# caches

mutable struct NLNewtonCache{uType,tType,tType2,rateType,J,W,ufType,jcType,lsType} <: AbstractNLSolverCache
  ustep::uType
  tstep::tType
  k::rateType
  atmp::uType
  dz::uType
  J::J
  W::W
  new_W::Bool
  firststage::Bool
  firstcall::Bool
  W_γdt::tType
  du1::uType
  uf::ufType
  jac_config::jcType
  linsolve::lsType
  weight::uType
  invγdt::tType2
  new_W_γdt_cutoff::tType
  J_t::tType
end

mutable struct NLNewtonConstantCache{tType,tType2,J,W,ufType} <: AbstractNLSolverCache
  tstep::tType
  J::J
  W::W
  new_W::Bool
  firststage::Bool
  firstcall::Bool
  W_γdt::tType
  uf::ufType
  invγdt::tType2
  new_W_γdt_cutoff::tType
  J_t::tType
end

mutable struct NLFunctionalCache{uType,tType,rateType} <: AbstractNLSolverCache
  ustep::uType
  tstep::tType
  k::rateType
  atmp::uType
  dz::uType
end

mutable struct NLFunctionalConstantCache{tType} <: AbstractNLSolverCache
  tstep::tType
end

mutable struct NLAndersonCache{uType,tType,rateType,uEltypeNoUnits} <: AbstractNLSolverCache
  ustep::uType
  tstep::tType
  k::rateType
  atmp::uType
  dz::uType
  """residuals `g(zprev) - zprev` of previous fixed-point iteration"""
  dzold::uType
  """value `g(zprev)` of previous fixed-point iteration"""
  z₊old::uType
  Δz₊s::Vector{uType}
  Q::Matrix{uEltypeNoUnits}
  R::Matrix{uEltypeNoUnits}
  γs::Vector{uEltypeNoUnits}
  history::Int
  aa_start::Int
  droptol::Union{Nothing,tType}
end

mutable struct NLAndersonConstantCache{uType,tType,uEltypeNoUnits} <: AbstractNLSolverCache
  tstep::tType
  dz::uType
  """residuals `g(zprev) - zprev` of previous fixed-point iteration"""
  dzold::uType
  """value `g(zprev)` of previous fixed-point iteration"""
  z₊old::uType
  Δz₊s::Vector{uType}
  Q::Matrix{uEltypeNoUnits}
  R::Matrix{uEltypeNoUnits}
  γs::Vector{uEltypeNoUnits}
  history::Int
  aa_start::Int
  droptol::Union{Nothing,tType}
end


text:
function run(b::Buffer)
    b.state[:log] = ""
    r = read(STDIN, Char)
    handle_raw(b, r)
    render(b)
    h = bottom(b)-top(b)
    if mode(b)==command_mode
        write(STDOUT, string(":", b.state[:command]))
        move_sys_cursor(h+3, 1+parse(b.state[:command_insert]))
    elseif mode(b)==search_mode
        write(STDOUT, string("/", b.state[:search]))
        move_sys_cursor(h+3, 1+parse(b.state[:search_insert]))
    else
        move_sys_cursor(b.cursor.pos[1]-top(b)+1, b.cursor.pos[2]-left(b)+1)
    end
end

function init(text)
    cursor = Cursor([1,1])
    state = Dict(:actions => "")
    state[:running]  = "false"
    state[:console]  = ""
    state[:register] = "\""
    state[:log]      = ""
    state[:top]      = "1"
    state[:bottom]   = "19"
    state[:left]     = "0"
    state[:right]    = "80"
    state[:syntax]   = "true"
    state[:undo]     = ""
    init_call(state, :command)
    init_call(state, :search)
    registers = Dict{Char, String}()
    Buffer(text,
           [""],
           [""],
           cursor,
           [normal_mode, normal_mode],
           state,
           registers,
           Array{Char, 1}())
end

function handle_raw(b::Buffer, c::Char)
    b.state[:actions] = string(b.state[:actions], c)
    mode(b).actions[c](b)
    xp = (x(b)>0 ? xs(b) : "0-1")
    b.state[:console] = string(mode(b).signature,
                               " ",ys(b), ",", xp, " ",
                               join(b.args, ""), "\n", b.state[:log])
end

function onopen()
end
function onattach()
end


text:
"""
    MattssonAlmquistVanDerWeide2018Minimal()

Coefficients of the optimized SBP operators with nonuniform grid given in
- Mattsson, Almquist, van der Weide (2018)
  Boundary optimized diagonal-norm SBP operators.
  Journal of Computational Physics 374, pp. 1261-1266.
"""
struct MattssonAlmquistVanDerWeide2018Minimal <: SourceOfCoefficients end

function Base.show(io::IO, source::MattssonAlmquistVanDerWeide2018Minimal)
  if get(io, :compact, false)
    summary(io, source)
  else
    print(io,
        "Mattsson, Almquist, van der Weide (2018) \n",
        "  Boundary optimized diagonal-norm SBP operators ('Minimal'). \n",
        "  Journal of Computational Physics 374, pp. 1261-1266.")
  end
end

struct BoundaryAdaptedGrid{T,M,Grid} <: AbstractArray{T,1}
  xmin::T
  xmax::T
  xstart::SVector{M,T}
  uniform_grid::Grid
end

function BoundaryAdaptedGrid(xmin::T, xmax::T, _xstart::SVector{M,T}, N::Int) where {M,T}
  @argcheck xmin < xmax
  @argcheck N > 2M

  Δx = (xmax - xmin) / (2*_xstart[end] + N + 1 - 2M)
  xstart = Δx .* _xstart
  uniform_grid = range(xmin + xstart[end] + Δx, xmax - xstart[end] - Δx, length=N-2M)
  BoundaryAdaptedGrid{T,M,typeof(uniform_grid)}(xmin, xmax, xstart, uniform_grid)
end

function BoundaryAdaptedGrid(xmin, xmax, xstart, N)
  T = promote_type(typeof(xmin), typeof(xmax), eltype(xstart))
  BoundaryAdaptedGrid(T(xmin), T(xmax), T.(xstart), N)
end

function Base.length(grid::BoundaryAdaptedGrid)
  length(grid.uniform_grid) + 2 * length(grid.xstart)
end

function Base.getindex(grid::BoundaryAdaptedGrid, i::Int)
  N = length(grid)
  M = length(grid.xstart)
  @boundscheck begin
    @argcheck i > 0
    @argcheck i <= N
  end

  if i <= M
    grid.xmin + grid.xstart[i]
  elseif i >= N-M+1
    grid.xmax - grid.xstart[1+N-i]
  else
    @inbounds x = getindex(grid.uniform_grid, i-M)
    x
  end
end

Base.size(grid::BoundaryAdaptedGrid) = (length(grid),)
Base.step(grid::BoundaryAdaptedGrid) = step(grid.uniform_grid)

function construct_grid(::MattssonAlmquistVanDerWeide2018Minimal, accuracy_order, xmin, xmax, N)
  @argcheck N > 6
  T = promote_type(typeof(xmin), typeof(xmax))

  if accuracy_order == 4
    xstart = SVector(
      T(0.0000000000000e+00),
      T(7.7122987842562e-01),
      T(1.7712298784256e+00),
      T(2.7712298784256e+00),
    )
  elseif accuracy_order == 6
    xstart = SVector(
      T(0.0000000000000e+00),
      T(4.0842950991998e-01),
      T(1.1968523189207e+00),
      T(2.1968523189207e+00),
      T(3.1968523189207e+00),
      T(4.1968523189207e+00),
    )
  elseif accuracy_order == 8
    xstart = SVector(
      T(0.0000000000000e+00),
      T(4.9439570885261e-01),
      T(1.4051531374839e+00),
      T(2.4051531374839e+00),
      T(3.4051531374839e+00),
      T(4.4051531374839e+00),
      T(5.4051531374839e+00),
    )
  else
    throw(ArgumentError("Accuracy order $accuracy_order not implemented/derived."))
  end

  BoundaryAdaptedGrid(xmin, xmax, xstart, N)
end

function first_derivative_coefficients(source::MattssonAlmquistVanDerWeide2018Minimal,
                                       order::Int, T=Float64, mode=FastMode())
  if order == 4
    left_boundary = (
        # d1
        DerivativeCoefficientRow{T,1,3}(SVector(T(-1.8612097181862932),
                                                T(2.2966302628676907),
                                                T(-0.4354205446813973) )),
        # d2
        DerivativeCoefficientRow{T,1,4}(SVector(T(-0.611186527373936),
                                                T(0),
                                                T(0.6937383659498486),
                                                T(-0.08255183857591565) )),
        # d3
        DerivativeCoefficientRow{T,1,5}(SVector(T(0.11778272127063323),
                                                T(-0.7051567887047331),
                                                T(0),
                                                T(0.6712846484961179),
                                                T(-0.08391058106201474) )),
    )
    right_boundary = .- left_boundary
    upper_coef = SVector(T(2//3), T(-1//12))
    central_coef = zero(T)
    lower_coef = -upper_coef
    left_weights = SVector( T(2.6864248295847e-01),
                            T(1.0094667153500e+00),
                            T(9.9312068011715e-01) )
    right_weights = left_weights
    left_boundary_derivatives = Tuple{}()
    right_boundary_derivatives = left_boundary_derivatives

    DerivativeCoefficients(left_boundary, right_boundary,
                            left_boundary_derivatives, right_boundary_derivatives,
                            lower_coef, central_coef, upper_coef,
                            left_weights, right_weights, mode, 1, order, source)
  elseif order == 6
    left_boundary = (
        # d1
        DerivativeCoefficientRow{T,1,5}(SVector(T(-3.924566448353279),
                                                T(4.962014957077398),
                                                T(-1.2881968205658438),
                                                T(0.28645730739095715),
                                                T(-0.03570899554923223) )),
        # d2
        DerivativeCoefficientRow{T,1,5}(SVector(T(-1.022587534557978),
                                                T(0),
                                                T(1.3023996740158794),
                                                T(-0.3357380168806994),
                                                T(0.05592587742279807) )),
        # d3
        DerivativeCoefficientRow{T,1,6}(SVector(T(0.17402334848464412),
                                                T(-0.8537429915394141),
                                                T(0),
                                                T(0.8419418935035105),
                                                T(-0.179894658468737),
                                                T(0.017672408019993497) )),
        # d4
        DerivativeCoefficientRow{T,1,7}(SVector(T(-0.036159059808532096),
                                                T(0.20564365765539105),
                                                T(-0.7867088686169533),
                                                T(0),
                                                T(0.749328780370319),
                                                T(-0.1486175733002554),
                                                T(0.01651306370002838) )),
        # d5
        DerivativeCoefficientRow{T,1,8}(SVector(T(0.004554664962414826),
                                                T(-0.03461379228073425),
                                                T(0.16985243300141972),
                                                T(-0.7571710333787599),
                                                T(0),
                                                T(0.7508648039541874),
                                                T(-0.15017296079083747),
                                                T(0.016685884532315277) )),
    )
    right_boundary = .- left_boundary
    upper_coef = SVector(T(3//4), T(-3//20), T(1//60))
    central_coef = zero(T)
    lower_coef = -upper_coef
    left_weights = SVector( T(0.12740260779883),
                            T(0.61820981002054),
                            T(0.94308973897679),
                            T(1.0093019060199),
                            T(0.99884825610465) )
    right_weights = left_weights
    left_boundary_derivatives = Tuple{}()
    right_boundary_derivatives = left_boundary_derivatives

    DerivativeCoefficients(left_boundary, right_boundary,
                            left_boundary_derivatives, right_boundary_derivatives,
                            lower_coef, central_coef, upper_coef,
                            left_weights, right_weights, mode, 1, order, source)
  elseif order == 8
    left_boundary = (
        # d1
        DerivativeCoefficientRow{T,1,6}(SVector(T(-3.4425783018277825),
                                                T(4.592216499078388),
                                                T(-1.5780326084197085),
                                                T(0.5114508898159292),
                                                T(-0.08276413469602352),
                                                T(-0.00029234395080269453) )),
        # d2
        DerivativeCoefficientRow{T,1,6}(SVector(T(-0.8677229176949711),
                                                T(0),
                                                T(1.1480053882687404),
                                                T(-0.34675581554193646),
                                                T(0.07194779962833377),
                                                T(-0.005474454660166564) )),
        # d3
        DerivativeCoefficientRow{T,1,7}(SVector(T(0.2312364265180344),
                                                T(-0.8902776878396409),
                                                T(0),
                                                T(0.8364384829720679),
                                                T(-0.21345201833077546),
                                                T(0.039658060554022004),
                                                T(-0.0036032638737076134) )),
        # d4
        DerivativeCoefficientRow{T,1,8}(SVector(T(-0.07428870791937443),
                                                T(0.26655320361241747),
                                                T(-0.8291108452061701),
                                                T(0),
                                                T(0.8031154600461715),
                                                T(-0.200795518743896),
                                                T(0.038098105611977996),
                                                T(-0.0035716974011229742) )),
        # d5
        DerivativeCoefficientRow{T,1,9}(SVector(T(0.012018140909062053),
                                                T(-0.0552909326936426),
                                                T(0.21152178487239182),
                                                T(-0.8028866425212018),
                                                T(0),
                                                T(0.8000791459238572),
                                                T(-0.19995806770071226),
                                                T(0.03808725099061186),
                                                T(-0.0035706797803698614)  )),
        # d6
        DerivativeCoefficientRow{T,1,10}(SVector(T(4.246147652870495e-5),
                                                T(0.004208071691358718),
                                                T(-0.03930901628651875),
                                                T(0.2007872492592992),
                                                T(-0.800274203791225),
                                                T(0),
                                                T(0.800027268482503),
                                                T(-0.20000681712062576),
                                                T(0.03809653659440491),
                                                T(-0.00357155030572546) )),
    )
    right_boundary = .- left_boundary
    upper_coef = SVector(T(4//5), T(-1//5), T(4//105), T(-1//280))
    central_coef = zero(T)
    lower_coef = -upper_coef
    left_weights = SVector( T(0.14523997892351),
                            T(0.76864793350174),
                            T(0.99116487068535),
                            T(0.99992473335107),
                            T(1.0002097054636),
                            T(0.99996591555866) )
    right_weights = left_weights
    left_boundary_derivatives = Tuple{}()
    right_boundary_derivatives = left_boundary_derivatives

    DerivativeCoefficients(left_boundary, right_boundary,
                            left_boundary_derivatives, right_boundary_derivatives,
                            lower_coef, central_coef, upper_coef,
                            left_weights, right_weights, mode, 1, order, source)
  # elseif order == 8
  #   left_boundary = (
  #       # d1
  #       DerivativeCoefficientRow{T,1,6}(SVector(T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T() )),
  #       # d2
  #       DerivativeCoefficientRow{T,1,6}(SVector(T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T() )),
  #       # d3
  #       DerivativeCoefficientRow{T,1,7}(SVector(T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T() )),
  #       # d4
  #       DerivativeCoefficientRow{T,1,8}(SVector(T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T() )),
  #       # d5
  #       DerivativeCoefficientRow{T,1,9}(SVector(T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T() )),
  #       # d6
  #       DerivativeCoefficientRow{T,1,10}(SVector(T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T() )),
  #   )
  #   right_boundary = .- left_boundary
  #   upper_coef = SVector(T(4//5), T(-1//5), T(4//105), T(-1//280))
  #   central_coef = zero(T)
  #   lower_coef = -upper_coef
  #   left_weights = SVector( T(),
  #                           T(),
  #                           T(),
  #                           T(),
  #                           T(),
  #                           T() )
  #   right_weights = left_weights
  #   left_boundary_derivatives = Tuple{}()
  #   right_boundary_derivatives = left_boundary_derivatives

  #   DerivativeCoefficients(left_boundary, right_boundary,
  #                           left_boundary_derivatives, right_boundary_derivatives,
  #                           lower_coef, central_coef, upper_coef,
  #                           left_weights, right_weights, mode, 1, order, source)
  else
    throw(ArgumentError("Order $order not implemented/derived."))
  end
end


text:
using MAT
using Random
using Statistics
include("../src/HyperLocal.jl")

Add path to wherever matrix is store. Loading take a long time.
s = time()
@time M = matread("../data/AmazonReview5core_H.mat")
NodeLabels = vec(M["NodeLabels"])
NodeNames = M["NodeNames"]
LabelNames = M["LabelNames"]
H = M["H"]
order = round.(Int64,vec(sum(H,dims=2)))

## Extract hyperedges with only one node
e = findall(x->x>1,order)
H = H[e,:]
d = vec(sum(H,dims=1))
order = round.(Int64,vec(sum(H,dims=2)))
volA = sum(d)
m,n = size(H)
Ht = sparse(H')
toc = time()-s
println("Done loading things into memory in $toc seconds.")

labels = [1; 2; 3; 12; 18; 17; 25; 15; 24]  # Labels for small nodes
lnum = length(labels)

# See outer parameters
seednum = 10
ntimes = 1
epsilon = 1.0
grownum = 200       # How much to grow seed set by using BestNeighbors
deltas = 10 .^LinRange(0,3,10)
enum = length(deltas)

for lab = 1:9
    label = labels[lab]
    T = findall(x->x ==label,NodeLabels)
    nT = length(T)

    # Different seed set sizes depending on dataset
    if lab < 5
        grownum = 200
        seednum = 10
    elseif lab < 7
        grownum = 2000
        seednum = 50
    else
        grownum = 10000
        seednum = 200
    end

    # For each epsilon we store a different .mat file of outputs
    outputmat = "Output_VaryDelta/Label_$label"*"_$seednum"*"_$grownum.mat"
    println(outputmat)

    # Output from HyperLocal
    hl_pr = zeros(enum)
    hl_re = zeros(enum)
    hl_f1 = zeros(enum)
    hl_time = zeros(enum)
    hl_size = zeros(enum)
    newS = zeros(enum)
    hl_cond = zeros(enum)

    # Generate a new seed set
    p = randperm(nT)
    Rstart = T[p[1:seednum]]
    OneHop = get_immediate_neighbors(H,Ht,Rstart)
    Rmore = BestNeighbors(H,d,Rstart,OneHop,grownum)
    R = union(Rmore,Rstart)

    # Force seed nodes to be contained in output set
    Rs = findall(x->in(x,Rstart),R)

    for e = 1:length(deltas)             # Try each seed set with each delta

        delta = deltas[e]

        # Run HyperLocal
        s = time()
        S, lcond = HyperLocal(H,Ht,order,d,R,epsilon,delta,Rs,true)
        hl_time[e] = time()-s
        condS, volS, cutS = tl_cond(H,S,d,1.0,volA,order)
        pr, re, f1 = PRF(T,S)
        hl_pr[e] = pr
        hl_re[e] = re
        hl_f1[e] = f1
        hl_size[e] = length(S)
        hl_cond[e] = condS
        nS =  length(setdiff(S,R))
        newS[e] = nS

        println("$label ($nT): $f1 \t $nS \t $epsilon")
    end

    matwrite(outputmat, Dict("hl_size"=>hl_size, "newS"=>newS, "hl_time"=>hl_time,
    "hl_pr"=>hl_pr, "hl_re"=>hl_re, "hl_f1"=>hl_f1, "hl_cond"=>hl_cond))

end

max_stars_repo_name: nveldt/HypergraphFlowClustering

max_issues_repo_name: nveldt/HypergraphFlowClustering

max_forks_repo_name: nveldt/HypergraphFlowClustering


text:
# This file contains the function generator tools to drive other tools of DsSimulator.

import UUIDs: uuid4

"""
    @def_source ex

where `ex` is the expression to define to define a new AbstractSource component type. The usage is as follows:
```julia
@def_source struct MySource{T1,T2,T3,...,TN,OP, RO} <: AbstractSource
    param1::T1 = param1_default     # optional field
    param2::T2 = param2_default     # optional field
    param3::T3 = param3_default     # optional field
        ⋮
    paramN::TN = paramN_default     # optional field
    output::OP = output_default     # mandatory field
    readout::RO = readout_function  # mandatory field
end
```
Here, `MySource` has `N` parameters, an `output` port and a `readout` function.

!!! warning
    `output` and `readout` are mandatory fields to define a new source. The rest of the fields are the parameters of the source.

!!! warning
    `readout` must be a single-argument function, i.e. a function of time `t`.

!!! warning
    New source must be a subtype of `AbstractSource` to function properly.

# Example
```julia
julia> @def_source struct MySource{OP, RO} <: AbstractSource
       a::Int = 1
       b::Float64 = 2.
       output::OP = Outport()
       readout::RO = t -> (a + b) * sin(t)
       end

julia> gen = MySource();

julia> gen.a
1

julia> gen.output
1-element Outport{Outpin{Float64}}:
 Outpin(eltype:Float64, isbound:false)
```
"""
macro def_source(ex)
    ex.args[2].head == :(<:) && ex.args[2].args[2] == :AbstractSource ||
        error("Invalid usage. The type should be a subtype of AbstractSource.\n$ex")
    foreach(nex -> appendex!(ex, nex), [
        :( trigger::$TRIGGER_TYPE_SYMBOL = Inpin() ),
        :( handshake::$HANDSHAKE_TYPE_SYMBOL = Outpin{Bool}() ),
        :( callbacks::$CALLBACKS_TYPE_SYMBOL = nothing ),
        :( name::Symbol = Symbol() ),
        :( id::$ID_TYPE_SYMBOL = Causal.uuid4() )
        ])
    quote
        Base.@kwdef $ex
    end |> esc
end

##### Define Sources library
"""
    $TYPEDEF

Constructs a generic function generator with `readout` function and `output` port.

# Fields

    $TYPEDFIELDS

# Example
```jldoctest
julia> gen = FunctionGenerator(readout = t -> [t, 2t], output = Outport(2));

julia> gen.readout(1.)
2-element Array{Float64,1}:
 1.0
 2.0
```
"""
@def_source struct FunctionGenerator{RO, OP<:Outport} <: AbstractSource
    "Readout function"
    readout::RO
    "Output port"
    output::OP = Outport(1)
end

 """
    $SIGNATURES

Constructs a `SinewaveGenerator` with output of the form
```math
    x(t) = A sin(2 \\pi f  (t - \\tau) + \\phi) + B
```
where ``A`` is `amplitude`, ``f`` is `frequency`, ``\\tau`` is `delay` and ``\\phi`` is `phase` and ``B`` is `offset`.

# Fields

    $TYPEDFIELDS
"""
@def_source struct SinewaveGenerator{T1 <: Real,
                                     T2 <: Real,
                                     T3 <: Real,
                                     T4 <: Real,
                                     T5 <: Real,
                                     OP <: Outport,
                                     RO} <: AbstractSource
   "Amplitude"
    amplitude::T1 = 1.
    "Frequency"
    frequency::T2 = 1.
    "Phase"
    phase::T3 = 0.
    "Delay in seconds"
    delay::T4 = 0.
    "Offset"
    offset::T5 = 0.
    "Output port"
    output::OP = Outport()
    "Readout function"
    readout::RO = (t, amplitude=amplitude, frequency=frequency, delay=delay, offset=offset) ->
        amplitude * sin(2 * pi * frequency * (t - delay) + phase) + offset
end

"""
    $TYPEDEF

Constructs a `DampedSinewaveGenerator` which generates outputs of the form
```math
    x(t) = A e^{\\alpha t} sin(2 \\pi f (t - \\tau) + \\phi) + B
```
where ``A`` is `amplitude`, ``\\alpha`` is `decay`, ``f`` is `frequency`, ``\\phi`` is `phase`, ``\\tau`` is `delay`
and ``B`` is `offset`.

# Fields

    $TYPEDFIELDS
"""
@def_source struct DampedSinewaveGenerator{T1 <: Real,
                                           T2 <: Real,
                                           T3 <: Real,
                                           T4 <: Real,
                                           T5 <: Real,
                                           T6 <: Real,
                                           OP <: Outport,
                                           RO} <: AbstractSource
    "Amplitude"
    amplitude::T1 = 1.
    "Attenuation rate"
    decay::T2 = 0.5
    "Frequency"
    frequency::T3 = 1.
    "Phase"
    phase::T4 = 0.
    "Delay in seconds"
    delay::T5 = 0.
    "Offset"
    offset::T6 = 0.
    "Output port"
    output::OP = Outport()
    "Readout funtion"
    readout::RO = (t, amplitude=amplitude, decay=decay, frequency=frequency, phase=phase, delay=delay, offset=offset) ->
        amplitude * exp(decay * t) * sin(2 * pi * frequency * (t - delay)) + offset
end

"""
    $TYPEDEF

Constructs a `SquarewaveGenerator` with output of the form
```math
    x(t) = \\left\\{\\begin{array}{lr}
	A_1 + B, &  kT + \\tau \\leq t \\leq (k + \\alpha) T + \\tau \\\\
	A_2 + B,  &  (k + \\alpha) T + \\tau \\leq t \\leq (k + 1) T + \\tau
	\\end{array} \\right. \\quad k \\in Z
```
where ``A_1``, ``A_2`` is `level1` and `level2`, ``T`` is `period`, ``\\tau`` is `delay` ``\\alpha`` is `duty`.

# Fields

    $TYPEDFIELDS
"""
@def_source struct SquarewaveGenerator{T1 <: Real,
                                       T2 <: Real,
                                       T3 <: Real,
                                       T4 <: Real,
                                       T5 <: Real,
                                       OP <: Outport,
                                       RO} <: AbstractSource
    "High level"
    high::T1 = 1.
    "Low level"
    low::T2 = 0.
    "Period"
    period::T3 = 1.
    "Duty cycle given in range (0, 1)"
    duty::T4 = 0.5
    "Delay in seconds"
    delay::T5 = 0.
    "Output port"
    output::OP = Outport()
    "Readout function"
    readout::RO = (t, high=high, low=low, period=period, duty=duty, delay=delay) ->
        t <= delay ? low : ( ((t - delay) % period <= duty * period) ? high : low )
end

"""
    $TYPEDEF

Constructs a `TriangularwaveGenerator` with output of the form
```math
    x(t) = \\left\\{\\begin{array}{lr}
	\\dfrac{A t}{\\alpha T} + B, &  kT + \\tau \\leq t \\leq (k + \\alpha) T + \\tau \\\\[0.25cm]
	\\dfrac{A (T - t)}{T (1 - \\alpha)} + B,  &  (k + \\alpha) T + \\tau \\leq t \\leq (k + 1) T + \\tau
	\\end{array} \\right. \\quad k \\in Z
```
where ``A`` is `amplitude`, ``T`` is `period`, ``\\tau`` is `delay` ``\\alpha`` is `duty`.

# Fields

    $TYPEDFIELDS
"""
@def_source struct TriangularwaveGenerator{T1 <: Real,
                                           T2 <: Real,
                                           T3 <: Real,
                                           T4 <: Real,
                                           T5 <: Real,
                                           OP <: Outport,
                                           RO} <: AbstractSource
    "Amplitude"
    amplitude::T1 =  1.
    "Period"
    period::T2 = 1.
    "Duty cycle"
    duty::T3 = 0.5
    "Delay in seconds"
    delay::T4 = 0.
    "Offset"
    offset::T5 = 0.
    "Output port"
    output::OP = Outport()
    "Readout function"
    readout::RO = (t, amplitude=amplitude, period=period, duty=duty, delay=delay, offset=offset) -> begin
        if t <= delay
            return offset
        else
            t = (t - delay) % period
            if t <= duty * period
                amplitude / (duty * period) * t + offset
            else
                (amplitude * (period - t)) / (period * (1 - duty)) + offset
            end
        end
    end
end

"""
    $TYPEDEF

Constructs a `ConstantGenerator` with output of the form
```math
    x(t) = A
```
where ``A`` is `amplitude.

# Fields

    $TYPEDFIELDS
"""
@def_source struct ConstantGenerator{T1 <: Real, OP <: Outport, RO} <: AbstractSource
    "Amplitude"
    amplitude::T1 = 1.
    "Output port"
    output::OP = Outport()
    "Readout function"
    readout::RO = (t, amplitude=amplitude) -> amplitude
end

"""
    $TYPEDEF

Constructs a `RampGenerator` with output of the form
```math
    x(t) = \\alpha (t - \\tau)
```
where ``\\alpha`` is the `scale` and ``\\tau`` is `delay`.

# Fields

    $TYPEDFIELDS
"""
@def_source struct RampGenerator{T1 <: Real,
                                 T2 <: Real,
                                 T3 <: Real,
                                 OP <: Outport,
                                 RO} <: AbstractSource
    "Scale"
    scale::T1 = 1.
    "Delay in seconds"
    delay::T2 = 0.
    "Offset"
    offset::T3 = 0.
    "Output port"
    output::OP = Outport()
    "Readout function"
    readout::RO = (t, scale=scale, delay=delay, offset=offset) ->  scale * (t - delay) + offset
end

"""
    $TYPEDEF

Constructs a `StepGenerator` with output of the form
```math
    x(t) = \\left\\{\\begin{array}{lr}
	B, &  t \\leq \\tau  \\\\
	A + B,  &  t > \\tau
	\\end{array} \\right.
```
where ``A`` is `amplitude`, ``B`` is the `offset` and ``\\tau`` is the `delay`.

# Fields

    $TYPEDFIELDS
"""
@def_source struct StepGenerator{T1 <: Real,
                                 T2 <: Real,
                                 T3 <: Real,
                                 OP <: Outport,
                                 RO} <: AbstractSource
    "Amplitude"
    amplitude::T1 = 1.
    "Delay in seconds"
    delay::T2 = 0.
    "Offset"
    offset::T3 = 0.
    "Output port"
    output::OP = Outport()
    "Readout function"
    readout::RO = (t, amplitude=amplitude, delay=delay, offset=offset) ->
        t - delay >= 0 ? amplitude + offset : offset
end

"""
    $TYPEDEF

Constructs an `ExponentialGenerator` with output of the form
```math
    x(t) = A e^{\\alpha (t - \\tau)}
```
where ``A`` is `scale`, ``\\alpha`` is `decay` and ``\\tau`` is `delay`.

# Fields

    $TYPEDFIELDS
"""
@def_source struct ExponentialGenerator{T1 <: Real,
                                        T2 <: Real,
                                        T3 <: Real,
                                        T4 <: Real,
                                        OP <: Outport,
                                        RO} <: AbstractSource
    "Scale"
    scale::T1 = 1.
    "Attenuation decay"
    decay::T2 = -1.
    "Delay in seconds"
    delay::T3 = 0.
    "Offset"
    offset::T4 = 0.
    "Output port"
    output::OP = Outport()
    "Readout function"
    readout::RO = (t, scale=scale, decay=decay, delay=delay, offset=offset) -> scale * exp(decay * (t - delay)) + offset
end

"""
    $TYPEDEF

Constructs an `DampedExponentialGenerator` with outpsuts of the form
```math
    x(t) = A (t - \\tau) e^{\\alpha (t - \\tau)}
```
where ``A`` is `scale`, ``\\alpha`` is `decay`, ``\\tau`` is `delay`.

# Fields

    $TYPEDFIELDS
"""
@def_source struct DampedExponentialGenerator{T1 <: Real,
                                              T2 <: Real,
                                              T3 <: Real,
                                              T4 <: Real,
                                              OP <: Outport,
                                              RO} <: AbstractSource
    "Scale"
    scale::T1 = 1.
    "Attenuation rate"
    decay::T2 = -1.
    "Delay in seconds"
    delay::T3 = 0.
    "Offet"
    offset::T4 = 0.
    "Output port"
    output::OP = Outport()
    "Reaodout function"
    readout::RO = (t, scale=scale, decay=decay, delay=delay, offset=offset) ->
        scale * (t - delay) * exp(decay * (t - delay)) + offset
end

##### Pretty-Printing of generators.

show(io::IO, gen::FunctionGenerator) = print(io,
    "FunctionGenerator(readout:$(gen.readout),  output:$(gen.output))")
show(io::IO, gen::SinewaveGenerator) = print(io,
    "SinewaveGenerator(amp:$(gen.amplitude), freq:$(gen.frequency), phase:$(gen.phase), ",
    "offset:$(gen.offset), delay:$(gen.delay))")
show(io::IO, gen::DampedSinewaveGenerator) = print(io,
    "DampedSinewaveGenerator(amp:$(gen.amplitude), decay:$(gen.delay), freq:$(gen.frequency), ",
    "phase:$(gen.phase), delay:$(gen.delay), offset:$(gen.offset))")
show(io::IO, gen::SquarewaveGenerator) = print(io,
    "SquarewaveGenerator(high:$(gen.high), low:$(gen.low), period:$(gen.period), duty:$(gen.duty), ",
    "delay:$(gen.delay))")
show(io::IO, gen::TriangularwaveGenerator) = print(io,
    "TriangularwaveGenerator(amp:$(gen.amplitude), period:$(gen.period), duty:$(gen.duty), ",
    "delay:$(gen.delay), offset:$(gen.offset))")
show(io::IO, gen::ConstantGenerator) = print(io,
    "ConstantGenerator(amp:$(gen.amplitude))")
show(io::IO, gen::RampGenerator) = print(io, "RampGenerator(scale:$(gen.scale), delay:$(gen.delay))")
show(io::IO, gen::StepGenerator) = print(io,
    "StepGenerator(amp:$(gen.amplitude), delay:$(gen.delay), offset:$(gen.offset))")
show(io::IO, gen::ExponentialGenerator) = print(io,
    "ExponentialGenerator(scale:$(gen.scale), decay:$(gen.decay), delay:$(gen.delay))")
show(io::IO, gen::DampedExponentialGenerator) = print(io,
    "DampedExponentialGenerator(scale:$(gen.scale), decay:$(gen.decay), delay:$(gen.delay))")