text:
# Load Qt framework
# adjust qtlibdir path and version info if necessary
# BUG: error: incomplete type 'QWidget' named in nested name specifier, etc
# This goes back to Cxx.jl which has several bugs in Cxx/test/qttest.jl

@static if is_apple()
    const qtlibdir = joinpath(homedir(),"Qt/5.7/clang_64/lib/")
    const QtCore = joinpath(qtlibdir,"QtCore.framework/")
    const QtWidgets = joinpath(qtlibdir,"QtWidgets.framework/")

    addHeaderDir(qtlibdir; isFramework = true, kind = C_System)

    Libdl.dlopen(joinpath(QtCore,"QtCore_debug"))
    addHeaderDir(joinpath(QtCore,"Headers"), kind = C_System)
    addHeaderDir(joinpath(QtCore,"Headers/5.7.0/QtCore"))
    cxxinclude(joinpath(QtCore,"Headers/5.7.0/QtCore/private/qcoreapplication_p.h"))

    Libdl.dlopen(joinpath(QtWidgets,"QtWidgets"))
    addHeaderDir(joinpath(QtWidgets,"Headers"), kind = C_System)
end

@static if is_linux()
    const qtincdir = isdir("/usr/include/qt5") ? "/usr/include/qt5" : "/usr/include/x86_64-linux-gnu/qt5"
    const qtlibdir = "/usr/lib/x86_64-linux-gnu/"
    const QtWidgets = joinpath(qtincdir,"QtWidgets")

    addHeaderDir(qtincdir, kind = C_System)
    addHeaderDir(QtWidgets, kind = C_System)

    Libdl.dlopen(joinpath(qtlibdir,"libQt5Core.so"), Libdl.RTLD_GLOBAL)
    Libdl.dlopen(joinpath(qtlibdir,"libQt5Gui.so"), Libdl.RTLD_GLOBAL)
    Libdl.dlopen(joinpath(qtlibdir,"libQt5Widgets.so"), Libdl.RTLD_GLOBAL)
end

cxxinclude("QApplication", isAngled=true)
cxxinclude("QMessageBox", isAngled=true)
cxxinclude("QFileDialog", isAngled=true)

# convertQString
cxx"""
    const char* convertQString(QString myString) {
        std::string s = myString.toLocal8Bit().constData();
        // alternatives: toAscii(), toLatin1(), toUtf8()
        // std::cout << s << std::endl;
        const char* c = s.c_str();
        return(c);
    }
"""

# Open file with QFileDialog, returns Julia String
function getOpenFileName()
    const a = "FileOpen"
    argv = Ptr{Uint8}[pointer(a),C_NULL]  # const char** argv
    argc = [int32(1)]                     # int agrc
    app = @cxx QApplication(*(pointer(argc)),pointer(argv))
    Qname = @cxx QFileDialog::getOpenFileName(cast(C_NULL, pcpp"QWidget"),
       pointer("Open File..."), pointer(homedir()), pointer("Image Files (*.png *.jpg *.jpeg *.bmp *.tif *tiff)"));
    bytestring(@cxx convertQString(Qname));  # somewhow necessary to remove the first string
    filename = bytestring(@cxx convertQString(Qname))
end

# Open file with QFileDialog, returns Julia String
function getSaveFileName()
    const a = "FileSave"
    argv = Ptr{Uint8}[pointer(a),C_NULL]  # const char** argv
    argc = [int32(1)]                     # int agrc
    app = @cxx QApplication(*(pointer(argc)),pointer(argv))
    Qname = @cxx QFileDialog::getSaveFileName(cast(C_NULL, pcpp"QWidget"),
       pointer("Save File..."), pointer(homedir()), pointer("Image Files (*.png *.jpg *.jpeg *.bmp *.tif *tiff)"));
    bytestring(@cxx convertQString(Qname));  # somewhow necessary to remove the first string
    filename = bytestring(@cxx convertQString(Qname))
end


text:
##### NEW KEYS #####

struct ShapeKey <: Gadfly.GuideElement
    title::AbstractString
    labels::Vector{String}
    pos::Vector
    visible::Bool
end
ShapeKey(;title="Shape", labels=String[], pos=[], visible=true) = ShapeKey(title, labels, pos, visible)
ShapeKey(v::Nothing) = ShapeKey(visible=false)
ShapeKey(title::AbstractString, labels::Vector{String}, pos::Vector) = ShapeKey(title, labels, pos, true)

"""
    Guide.shapekey[(; title="Shape", labels=String[], pos=Float64[])]
    Guide.shapekey(title, labels, pos)

Enable control of the auto-generated shapekey.  Set the key `title` and the item `labels`.
`pos` overrides [Theme(key_position=)](@ref Gadfly) and can be in either
relative (e.g. [0.7w, 0.2h] is the lower right quadrant), absolute (e.g. [0mm,
0mm]), or plot scale (e.g. [0,0]) coordinates. `Guide.shapekey(nothing)` will hide the key.
"""
const shapekey = ShapeKey

function render(guide::ShapeKey, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)

    (theme.key_position == :none || !guide.visible || isa(aes.shape[1], Function)) && return PositionedGuide[]
    gpos = guide.pos
    (theme.key_position == :inside) && isempty(gpos) &&  (gpos = [0.7w, 0.25h])

    # Aesthetics for keys: shape_key_title, shape_label (Function), shape_key_shapes (AbstractDict)
    nshapes = length(unique(aes.shape))
    guide_title = (guide.title≠"Shape" || aes.shape_key_title==nothing) ? guide.title : aes.shape_key_title
    shape_key_labels = isempty(guide.labels) ? aes.shape_label(1:nshapes) : guide.labels

    colors = Colorant[]
    if (aes.shape_key_title !=nothing)  && (aes.color_key_title==aes.shape_key_title)
        colors = collect(keys(aes.color_key_colors))
    end
    shapes = theme.point_shapes[1:nshapes]

    title_context, title_width = render_key_title2(guide_title, theme)
    ctxs = render_discrete_key(shape_key_labels, title_context, title_width, theme, shapes=shapes, colors=colors)

    position, stackable = right_guide_position, true
    if !isempty(gpos)
        position, stackable = over_guide_position, false
        ctxs = [compose(context(), (context(gpos[1],gpos[2]), ctxs[1]))]
    elseif theme.key_position == :left
        position = left_guide_position
    elseif theme.key_position == :top
        position = top_guide_position
    elseif theme.key_position == :bottom
        position = bottom_guide_position
    end

    return [PositionedGuide(ctxs, 0, position, stackable)]
end

function render_discrete_key(labels::Vector{String}, title_ctx::Context, title_width::Measure, theme::Gadfly.Theme;
    colors=Colorant[], aes_color_label=nothing, shapes=Function[], sizes=Measure[])

    n = max(length(colors), length(shapes), length(sizes))
    n==0 && (n = length(labels))

    isempty(shapes) && (shapes = fill((theme.key_swatch_shape==nothing) ? theme.point_shapes[1] :  theme.key_swatch_shape, n))
    isempty(colors) && (colors = fill((theme.key_swatch_color==nothing) ? theme.default_color : theme.key_swatch_color, n))
    isempty(sizes) && (sizes = fill((theme.key_swatch_size==nothing) ? theme.point_size : theme.key_swatch_size, n))

    # only consider layouts with a reasonable number of columns
    maxcols = theme.key_max_columns < 1 ? 1 : theme.key_max_columns
    maxcols = min(n, maxcols)

    extents = text_extents(theme.key_label_font, theme.key_label_font_size, values(labels)...)
    text_widths, text_heights = first.(extents), last.(extents)

    ypad = 1.0mm
    title_height = title_ctx.box.a[2]
    swatch_size = 2*maximum(sizes)
    entry_height = max(swatch_size, maximum(text_heights)) + ypad

    # return a context with a lyout of numcols columns
    function make_layout(numcols)
        colrows = Array{Int}(undef, numcols)
        m = n
        for i in 1:numcols
            colrows[i] = ceil(Int, m/(1+numcols-i))
            m -= colrows[i]
        end

        xpad = 1mm
        colwidths = Array{Measure}(undef, numcols)
        m = 0
        for (i, nrows) in enumerate(colrows)
            if m == n
                colwidths[i] = 0mm
            else
                colwidth = maximum(text_widths[m+1:m+nrows])
                colwidth += entry_height + xpad
                colwidths[i] = colwidth
                m += nrows
            end
        end

        ctxwidth = sum(colwidths)
        ctxheight = entry_height * colrows[1] + title_height

        ctxp = ctxpromise() do draw_context
            yoff = 0.5h - ctxheight/2
            outerctx = context()

            compose!(outerctx, (context(xpad, yoff), title_ctx))

            ctx = context(0, yoff+title_height, ctxwidth, ctxheight-title_height, units=UnitBox(0, 0, 1, colrows[1]))

            m = 0
            xpos = 0w
            for (colwidth, nrows) in zip(colwidths, colrows)

                x = [0.5cy]
                clrs = colors[m+1:m+nrows]
                shps = shapes[m+1:m+nrows]
                szs = sizes[m+1:m+nrows]

                swatches_shapes = [f(x, [y-0.5].*cy, [s]) for (y,(f,s)) in enumerate(zip(shps, szs))]
                sw1 = [(context(), s, fill(c), fillopacity(theme.alphas[1]), stroke(theme.discrete_highlight_color(c)))
                     for (s,c) in zip(swatches_shapes, clrs)]
                swatches = compose!(context(), linewidth(theme.highlight_width), sw1...)

                swatch_labels = compose!(context(),
                text([entry_height+xpad], [0.5:nrows;]*cy, labels[m+1:m+nrows], [hleft], [vcenter]),
                    font(theme.key_label_font),
                    fontsize(theme.key_label_font_size),
                    fill(theme.key_label_color))

                col = compose!(context(xpos, 0, colwidth), swatches, swatch_labels)

                if aes_color_label != nothing
                    classes = [svg_color_class_from_label(aes_color_label([c])[1]) for c in clrs]
                    #class_jscalls = ["data(\"color_class\", \"$(c)\")" for c in classes]
                    compose!(col,
                        svgclass(classes),
                        jscall(["""
                            data(\"color_class\", \"$(c)\")
                            .click(Gadfly.colorkey_swatch_click)
                            """ for c in classes]))
                end
                compose!(ctx, col)

                m += nrows
                xpos += colwidth
            end

            return compose!(outerctx, ctx,
                            # defeat webkit's asinine default drag behavior
                            jscall("drag(function() {}, function() {}, function() {})"),
                            svgclass("guide colorkey"))
        end

        return compose!(context(minwidth=max(title_width, ctxwidth),
                    minheight=ctxheight, units=UnitBox()),
                    ctxp)
    end

    return map(make_layout, 1:maxcols)
end

function render_key_title2(title::AbstractString, theme::Gadfly.Theme)
    title_width, title_height = max_text_extents(theme.key_title_font, theme.key_title_font_size, title)

    if theme.guide_title_position == :left
        title_form = text(0.0w, title_height, title, hleft, vbottom)
    elseif theme.guide_title_position == :center
        title_form = text(0.5w, title_height, title, hcenter, vbottom)
    elseif theme.guide_title_position == :right
        title_form = text(1.0w, title_height, title, hright, vbottom)
    else
        error("$(theme.guide_title_position) is not a valid guide title position")
    end

    title_padding = 1.5mm
    title_context = compose!(
        context(0w, 0h, 1w, title_height + title_padding),
        title_form,
        stroke(nothing),
        font(theme.key_title_font),
        fontsize(theme.key_title_font_size),
        fill(theme.key_title_color))

    return title_context, title_width
end

struct SizeKey <: Gadfly.GuideElement
    title::AbstractString
    labels::Vector{String}
    pos::Vector{Compose.MeasureOrNumber}
    visible::Bool
end
SizeKey(;title="Size", labels=String[], pos=[], visible=true) = SizeKey(title, labels, pos, visible)
SizeKey(v::Nothing) = SizeKey(visible=false)
SizeKey(title::AbstractString, labels::Vector{String}, pos::Vector) = SizeKey(title, labels, pos, true)

"""
    Guide.sizekey[(; title="size", labels=String[], pos=[])]
    Guide.sizekey(title, labels, pos)

Enable control of the sizekey.  Set the key `title` and the item `labels`.
`pos` overrides [Theme(key_position=)](@ref Gadfly) and can be in either
relative (e.g. [0.7w, 0.2h] is the lower right quadrant), absolute (e.g. [0mm,
0mm]), or plot scale (e.g. [0,0]) coordinates. Currently `Guide.sizekey` will only work by adding
`Scale.size_discrete2` to the plot. `Guide.sizekey(nothing)` will hide the key.
"""
const sizekey = SizeKey

function render(guide::SizeKey, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)

    (theme.key_position==:none || !guide.visible || aes.size_key_vals===nothing) && return PositionedGuide[]
    gpos = guide.pos
    (theme.key_position==:inside) && isempty(gpos) &&  (gpos = [0.7w, 0.25h])

    # Aesthetics for keys: size_key_title, size_label (Function), size_key_vals (AbstractDict)
    nsizes = length(unique(aes.size))
    guide_title = (guide.title≠"Size" || aes.size_key_title===nothing) ? guide.title : aes.size_key_title
    sizes = collect(keys(aes.size_key_vals))
    size_key_labels = isempty(guide.labels) ? aes.size_label(sizes) : guide.labels

    colors = Colorant[]
    if (aes.size_key_title≠nothing)  && (aes.color_key_title==aes.size_key_title)
        colors = collect(keys(aes.color_key_colors))
    end

    title_context, title_width = render_key_title2(guide_title, theme)
    ctxs = render_discrete_key(size_key_labels, title_context, title_width, theme, sizes=sizes, colors=colors)

    position, stackable = right_guide_position, true
    if !isempty(gpos)
        position, stackable = over_guide_position, false
        ctxs = [compose(context(), (context(gpos[1], gpos[2]), ctxs[1]))]
    elseif theme.key_position == :left
        position = left_guide_position
    elseif theme.key_position == :top
        position = top_guide_position
    elseif theme.key_position == :bottom
        position = bottom_guide_position
    end

    return [PositionedGuide(ctxs, 0, position, stackable)]
end

struct ManualDiscreteKey <: Gadfly.GuideElement
    title::AbstractString
    labels::Vector{String}
    pos::Vector
    colors::Vector{Colorant}
    shapes::Vector{Function}
    sizes::Vector{Measure}
    visible::Bool
end

function ManualDiscreteKey(;title="", labels=String[], pos=[], color=Colorant[], shape=Function[], size=Measure[])
    CT, ST, SZT = eltype(color), eltype(shape), eltype(size)
    ncolors, nshapes, nsizes =  length(color), length(shape), length(size)
    n = max(ncolors, nshapes, nsizes)
    ncolors==1 && (color = repeat(color, n))
    nshapes==1 && (shape = repeat(shape, n))
    nsizes==1 && (size = repeat(size, n))

    theme = Gadfly.current_theme()
    clrs = CT<:Int ? theme.discrete_colormap(maximum(color))[color] : Gadfly.parse_colorant(color)
    shps = ST<:Int ? theme.point_shapes[shape] : shape
    szs = SZT<:Int ? theme.discrete_sizemap(maximum(size))[size] : size

    cataes = [clrs, shps, szs]
    notempty = .!isempty.(cataes)
    if any(notempty)
        swatches = collect(Tuple, zip(cataes[notempty]...))
        !allunique(swatches) && error("Swatches should not be repeated in a manual key")
    end
    return ManualDiscreteKey(title, labels, pos, clrs, shps, szs,  true)
end

"""
    Guide.manual_discrete_key(title::String, labels::Vector{String};
        pos=[], color=Colorant[], shape=Function[], size=Measure[])

Manually define a discrete key with the legend `title` and `labels`, and swatch `color`, `shape` and `size`.
The swatch aesthetics can be Vectors of specific types (as above), or integer ranges.
Integer ranges refer to the order of items in the discrete Theme palettes [Discrete Scales](@ref).
Set the key position inside using `pos` (see [`Guide.sizekey`](@ref), [`Guide.shapekey`](@ref)).
"""
manual_discrete_key(title::AbstractString, labels::Vector{String}; kwargs...) =
    ManualDiscreteKey(title=title, labels=labels; kwargs...)

"""
    Guide.manual_color_key(title::String, labels::Vector{String}, color::AbstractVector=
        1:length(labels); pos=[], shape=Function[], size=Measure[])

Similar to [`Guide.manual_discrete_key`](@ref). `color` can be Ints, or Colors and/or Strings
(default is Ints).
"""
function manual_color_key(title::AbstractString, labels::Vector{String},
        color::AbstractVector=1:length(labels); kwargs...)
    ManualDiscreteKey(title=title, labels=labels, color=color; kwargs...)
end

function render(guide::ManualDiscreteKey, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    gpos = guide.pos
    (theme.key_position == :inside) && isempty(gpos) &&  (gpos = [0.7w, 0.25h])

    title_context, title_width = render_key_title2(guide.title, theme)
    ctxs = render_discrete_key(guide.labels, title_context, title_width, theme,
        shapes=guide.shapes, colors=guide.colors, sizes=guide.sizes)

    position, stackable = right_guide_position, true
    if !isempty(gpos)
        position, stackable = over_guide_position, false
        ctxs = [compose(context(), (context(gpos[1],gpos[2]), ctxs[1]))]
    elseif theme.key_position == :left
        position = left_guide_position
    elseif theme.key_position == :top
        position = top_guide_position
    elseif theme.key_position == :bottom
        position = bottom_guide_position
    end

    return [PositionedGuide(ctxs, 0, position, stackable)]
end

ext: jl


text:
function HHG(f::TVF) where {TVF<:AbstractVector{<:AbstractFloat}}
    n = length(f)
    spec = fft(f)
    return abs.(spec) .^ 2 / n^2
end
function HHGfreq(f::TVF, ncyc) where {TVF<:AbstractVector{<:AbstractFloat}}
    x = collect(range(1, length(f), length = length(f)))
    return x ./ ncyc
end


text:
using StaticArrays, Simplices, SimplexSplitting, Distributions, TimerOutputs

include("../example_embeddings.jl")
function potentially_intersecting_simplices2(t::Triangulation, image_i::Int)
    inds_potential_simplices = Int[]

    n_simplices = length(t.radii)

    @inbounds for i = 1:n_simplices
        dist_difference = ((t.centroids_im[image_i] - t.centroids[i]).' *
                            (t.centroids_im[image_i] - t.centroids[i]) - (t.radii_im[image_i] + t.radii[i])^2)[1]
        if dist_difference < 0
            push!(inds_potential_simplices, i)
        end
    end
    return inds_potential_simplices
end

"""
Returns a matrix of convex coefficients to construct points contained within
a `dim`-dimensional simplex. Coefficients can either be random (generated
according to a uniform distribution) or be constructed such that the resulting
points are uniformly distributed within the simplex (generated by shape-preserving
splitting a generic simplex). The `sample_randomly` argument controls this
behaviour.

The default is to not sample randomly (`sample_randomly = false`). In this case,
you won't get the exact number of convex coefficient sets you want, but the
smallest number of set of coefficients such that you have *at least* `n_pts`
convex coefficient sets.

If `sample_randomly = true`, you get as many convex coefficient combinations as
`n_pts`.
"""
function get_coeffs(dim::Int, n_randpts::Int, sample_randomly::Bool)
    # Create a set of convex coefficients to be used for all simplices
    if sample_randomly
        convex_coeffs = rand(Distributions.Uniform(0, 1), n_randpts, dim + 1)
        convex_coeffs .= convex_coeffs ./ sum(convex_coeffs, 2)
    else
        minimum_split_factor = ceil(Int, n_randpts^(1 / dim))[1]
        coeffs = Simplices.even_sampling_rules(dim, minimum_split_factor).'
        convex_coeffs = Size(size(coeffs, 1), dim + 1)(coeffs)
    end

    convex_coeffs.'
end

"""
Prepares the triangulation for a discrete approximation of the Markov matrix.
Creates alternative representations of the simplices that allow for efficient
(mostly non-allocating) checking if a point lies inside a simplex.
"""
function prepare_mm_dd(t::Triangulation)

    n_simplices = size(t.simplex_inds, 1)
    n_vertices = size(t.simplex_inds, 2)
    dim = n_vertices - 1

    # The [:, j, i] th entry of these two arrays holds the jth vertex of the
    # ith simplex, but instead of having just `dim` vertices, we append a `1`
    # to the end of the vectors. This allows for efficent (non-allocating)
    # computation within the `contains_point_lessalloc!` function. If we instead
    # would have appended the 1's inside that function, we would be performing
    # memory-allocating operations, which are very expensive. Doing this instead
    # gives orders of magnitude speed-ups for sufficiently large triangulations.
    S1 = Array{Float64}(dim + 1, dim + 1, n_simplices)
    IS1 = Array{Float64}(dim + 1, dim + 1, n_simplices)

    # Collect simplices in the form of (dim+1)^2-length column vectors. This
    # also helps with the
    simplices = Size((dim+1)^2, n_simplices)(zeros((dim+1)^2, n_simplices))
    imsimplices = Size((dim+1)^2, n_simplices)(zeros((dim+1)^2, n_simplices))

    for i in 1:n_simplices
        for j in 1:n_vertices
            S1[:, j, i] = vcat(t.points[t.simplex_inds[i, j], :], 1.0)
            IS1[:, j, i] = vcat(t.impoints[t.simplex_inds[i, j], :], 1.0)
        end

        simplices[:, i] = reshape(S1[:, :, i], (dim+1)^2)
        imsimplices[:, i] = reshape(IS1[:, :, i], (dim+1)^2)
    end

    return simplices, imsimplices
end

"""
Checks if a point lies inside a simplex `sx`. The simplex is represented as a
column vector [V1; V2; V3; V4], where V1 is the vertex represented as a
(dim+1)-length column vector, where the first 1:dim entries is the vertex
coordinates and the last entry is a 1. The reason for appending the 1 is for
efficiency when computing determinants (appending row vectors to matrices is
very expensive, because it allocates memory, so we rather reformat our simplices
beforehand).

`s_arr` is a temporary array of the same size as `sx`. `signs` is a
(dim + 1)-length vector containing the signs of the determinants of the
(at most) (dim + 1) matrices we need to check to determine whether the point
lies inside the simplex or not. `dim` is just the dimension of the space.
"""
function contains_point_lessalloc!(signs, s_arr, sx, point, dim::Int)
    # Redefine the temporary simplex. This is in-place, so we don't allocate
    # memory. We could also have re-initialised `signs`, but since we're never
    # comparing more than two consecutive signs, this is not necessary.
    s_arr .= 0.0
    s_arr .= s_arr .+ sx
    signs .= 0.0

    # Replace first vertex with the point
    s_arr[1:dim] = point

    # Compute the first sign
    signs[1] = sign(det(reshape(s_arr, dim + 1, dim + 1)))
    s_arr .= 0.0
    s_arr .= s_arr .+ sx

    for κ = 2:(dim) # Check remaining signs and stop if sign changes
        # Replace the ith vertex with the point we're cheking (leaving the
        # 1 appended to Vi intact.)
        s_arr[(dim + 1)*(κ - 1)+1:((dim + 1)*(κ - 1)+ 1 + dim - 1)] = point
        signs[κ] = sign(det(reshape(s_arr, dim + 1, dim + 1)))

        if !(signs[κ] == signs[κ - 1])
           return false
        end
        s_arr .= 0.0
        s_arr .= s_arr .+ sx
    end

    # Last the last vertex with the point in question
    s_arr[((dim + 1)*(dim)+1):((dim+1)^2-1)] = point
    signs[end] = sign(det(reshape(s_arr, dim + 1, dim + 1)))

    if !(signs[end] == signs[end - 1])
       return false
    end

    return true
end

max_stars_repo_name: kahaaga/InvariantDistribution

max_issues_repo_name: kahaaga/InvariantDistribution

max_forks_repo_name: kahaaga/InvariantDistribution


text:
module superoperators

export SuperOperator, DenseSuperOperator, SparseSuperOperator,
        spre, spost, liouvillian, expm

import Base: ==, *, /, +, -
import ..bases

using Compat
using ..bases, ..operators, ..operators_dense, ..operators_sparse

"""
Base class for all super operator classes.

Super operators are bijective mappings from operators given in one specific
basis to operators, possibly given in respect to another, different basis.
To embed super operators in an algebraic framework they are defined with a
left hand basis `basis_l` and a right hand basis `basis_r` where each of
them again consists of a left and right hand basis.
```math
A_{bl_1,bl_2} = S_{(bl_1,bl_2) ↔ (br_1,br_2)} B_{br_1,br_2}
\\\\
A_{br_1,br_2} = B_{bl_1,bl_2} S_{(bl_1,bl_2) ↔ (br_1,br_2)}
```
"""
abstract type SuperOperator end

"""
    DenseSuperOperator(b1[, b2, data])

SuperOperator stored as dense matrix.
"""
mutable struct DenseSuperOperator <: SuperOperator
    basis_l::Tuple{Basis, Basis}
    basis_r::Tuple{Basis, Basis}
    data::Matrix{Complex128}
    function DenseSuperOperator(basis_l::Tuple{Basis, Basis}, basis_r::Tuple{Basis, Basis}, data::Matrix{Complex128})
        if length(basis_l[1])*length(basis_l[2]) != size(data, 1) ||
           length(basis_r[1])*length(basis_r[2]) != size(data, 2)
            throw(DimensionMismatch())
        end
        new(basis_l, basis_r, data)
    end
end

function DenseSuperOperator(basis_l::Tuple{Basis, Basis}, basis_r::Tuple{Basis, Basis})
    Nl = length(basis_l[1])*length(basis_l[2])
    Nr = length(basis_r[1])*length(basis_r[2])
    data = zeros(Complex128, Nl, Nr)
    DenseSuperOperator(basis_l, basis_r, data)
end

"""
    SparseSuperOperator(b1[, b2, data])

SuperOperator stored as sparse matrix.
"""
mutable struct SparseSuperOperator <: SuperOperator
    basis_l::Tuple{Basis, Basis}
    basis_r::Tuple{Basis, Basis}
    data::SparseMatrixCSC{Complex128, Int}
    function SparseSuperOperator(basis_l::Tuple{Basis, Basis}, basis_r::Tuple{Basis, Basis}, data::SparseMatrixCSC{Complex128, Int})
        if length(basis_l[1])*length(basis_l[2]) != size(data, 1) ||
           length(basis_r[1])*length(basis_r[2]) != size(data, 2)
            throw(DimensionMismatch())
        end
        new(basis_l, basis_r, data)
    end
end

function SparseSuperOperator(basis_l::Tuple{Basis, Basis}, basis_r::Tuple{Basis, Basis})
    Nl = length(basis_l[1])*length(basis_l[2])
    Nr = length(basis_r[1])*length(basis_r[2])
    data = spzeros(Complex128, Nl, Nr)
    SparseSuperOperator(basis_l, basis_r, data)
end

SuperOperator(basis_l, basis_r, data::SparseMatrixCSC{Complex128, Int}) = SparseSuperOperator(basis_l, basis_r, data)
SuperOperator(basis_l, basis_r, data::Matrix{Complex128}) = DenseSuperOperator(basis_l, basis_r, data)

Base.copy(a::T) where {T<:SuperOperator} = T(a.basis_l, a.basis_r, copy(a.data))

Base.full(a::SparseSuperOperator) = DenseSuperOperator(a.basis_l, a.basis_r, full(a.data))
Base.full(a::DenseSuperOperator) = copy(a)

Base.sparse(a::DenseSuperOperator) = SparseSuperOperator(a.basis_l, a.basis_r, sparse(a.data))
Base.sparse(a::SparseSuperOperator) = copy(a)

==(a::T, b::T) where {T<:SuperOperator} = samebases(a, b) && (a.data == b.data)

Base.length(a::SuperOperator) = length(a.basis_l[1])*length(a.basis_l[2])*length(a.basis_r[1])*length(a.basis_r[2])
bases.samebases(a::SuperOperator, b::SuperOperator) = samebases(a.basis_l[1], b.basis_l[1]) && samebases(a.basis_l[2], b.basis_l[2]) &&
                                                      samebases(a.basis_r[1], b.basis_r[1]) && samebases(a.basis_r[2], b.basis_r[2])
bases.multiplicable(a::SuperOperator, b::SuperOperator) = multiplicable(a.basis_r[1], b.basis_l[1]) && multiplicable(a.basis_r[2], b.basis_l[2])
bases.multiplicable(a::SuperOperator, b::Operator) = multiplicable(a.basis_r[1], b.basis_l) && multiplicable(a.basis_r[2], b.basis_r)

# Arithmetic operations
function *(a::SuperOperator, b::DenseOperator)
    check_multiplicable(a, b)
    data = a.data*reshape(b.data, length(b.data))
    return DenseOperator(a.basis_l[1], a.basis_l[2], reshape(data, length(a.basis_l[1]), length(a.basis_l[2])))
end

function *(a::SuperOperator, b::SuperOperator)
    check_multiplicable(a, b)
    return SuperOperator(a.basis_l, b.basis_r, a.data*b.data)
end

*(a::SuperOperator, b::Number) = SuperOperator(a.basis_l, a.basis_r, a.data*b)
*(a::Number, b::SuperOperator) = SuperOperator(b.basis_l, b.basis_r, a*b.data)

/(a::SuperOperator, b::Number) = SuperOperator(a.basis_l, a.basis_r, a.data/b)

+(a::SuperOperator, b::SuperOperator) = (check_samebases(a, b); SuperOperator(a.basis_l, a.basis_r, a.data+b.data))

-(a::SuperOperator) = SuperOperator(a.basis_l, a.basis_r, -a.data)
-(a::SuperOperator, b::SuperOperator) = (check_samebases(a, b); SuperOperator(a.basis_l, a.basis_r, a.data-b.data))

"""
    spre(op)

Create a super-operator equivalent for right side operator multiplication.

For operators ``A``, ``B`` the relation

```math
    \\mathrm{spre}(A) B = A B
```

holds. `op` can be a dense or a sparse operator.
"""
spre(op::Operator) = SuperOperator((op.basis_l, op.basis_l), (op.basis_r, op.basis_r), tensor(op, identityoperator(op)).data)

"""
Create a super-operator equivalent for left side operator multiplication.

For operators ``A``, ``B`` the relation

```math
    \\mathrm{spost}(A) B = B A
```

holds. `op` can be a dense or a sparse operator.
"""
spost(op::Operator) = SuperOperator((op.basis_r, op.basis_r), (op.basis_l, op.basis_l), kron(transpose(op.data), identityoperator(op).data))

function _check_input(H::Operator, J::Vector, Jdagger::Vector, rates::Union{Vector{Float64}, Matrix{Float64}})
    for j=J
        @assert typeof(j) <: Operator
        check_samebases(H, j)
    end
    for j=Jdagger
        @assert typeof(j) <: Operator
        check_samebases(H, j)
    end
    @assert length(J)==length(Jdagger)
    if typeof(rates) == Matrix{Float64}
        @assert size(rates, 1) == size(rates, 2) == length(J)
    elseif typeof(rates) == Vector{Float64}
        @assert length(rates) == length(J)
    end
end

"""
    liouvillian(H, J; rates, Jdagger)

Create a super-operator equivalent to the master equation so that ``\\dot ρ = S ρ``

The super-operator ``S`` is defined by

```math
\\dot ρ = S ρ = -\\frac{i}{ħ} [H, ρ] + 2 J ρ J^† - J^† J ρ - ρ J^† J
```

# Arguments
* `rates`: Vector or matrix specifying the coefficients for the jump operators.
* `Jdagger`: Vector containing the hermitian conjugates of the jump operators. If they
             are not given they are calculated automatically.
"""
function liouvillian(H::T, J::Vector{T};
            rates::Union{Vector{Float64}, Matrix{Float64}}=ones(Float64, length(J)),
            Jdagger::Vector{T}=dagger.(J)) where T<:Operator
    _check_input(H, J, Jdagger, rates)
    L = spre(-1im*H) + spost(1im*H)
    if typeof(rates) == Matrix{Float64}
        for i=1:length(J), j=1:length(J)
            jdagger_j = rates[i,j]/2*Jdagger[j]*J[i]
            L -= spre(jdagger_j) + spost(jdagger_j)
            L += spre(rates[i,j]*J[i]) * spost(Jdagger[j])
        end
    elseif typeof(rates) == Vector{Float64}
        for i=1:length(J)
            jdagger_j = rates[i]/2*Jdagger[i]*J[i]
            L -= spre(jdagger_j) + spost(jdagger_j)
            L += spre(rates[i]*J[i]) * spost(Jdagger[i])
        end
    end
    return L
end

"""
    expm(op::DenseSuperOperator)

Operator exponential which can for example used to calculate time evolutions.
"""
Base.expm(op::DenseSuperOperator) = DenseSuperOperator(op.basis_l, op.basis_r, expm(op.data))

end # module


text:
"Represents an animation object"
struct Animation
    dir::String
    frames::Vector{String}
end

function Animation()
    tmpdir = convert(String, mktempdir())
    Animation(tmpdir, String[])
end

"""
    frame(animation[, plot])

Add a plot (the current plot if not specified) to an existing animation
"""
function frame(anim::Animation, plt::P=current()) where P<:AbstractPlot
    i = length(anim.frames) + 1
    filename = @sprintf("%06d.png", i)
    png(plt, joinpath(anim.dir, filename))
    push!(anim.frames, filename)
end

giffn() = (isijulia() ? "tmp.gif" : tempname()*".gif")
movfn() = (isijulia() ? "tmp.mov" : tempname()*".mov")
mp4fn() = (isijulia() ? "tmp.mp4" : tempname()*".mp4")

mutable struct FrameIterator
    itr
    every::Int
    kw
end
FrameIterator(itr; every=1, kw...) = FrameIterator(itr, every, kw)

"""
Animate from an iterator which returns the plot args each iteration.
"""
function animate(fitr::FrameIterator, fn = giffn(); kw...)
    anim = Animation()
    for (i, plotargs) in enumerate(fitr.itr)
        if mod1(i, fitr.every) == 1
            plot(wraptuple(plotargs)...; fitr.kw...)
            frame(anim)
        end
    end
    gif(anim, fn; kw...)
end

# most things will implement this
function animate(obj, fn = giffn(); every=1, fps=20, loop=0, kw...)
    animate(FrameIterator(obj, every, kw), fn; fps=fps, loop=loop)
end

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

"Wraps the location of an animated gif so that it can be displayed"
struct AnimatedGif
    filename::String
end

file_extension(fn) = Base.Filesystem.splitext(fn)[2][2:end]

gif(anim::Animation, fn = giffn(); kw...) = buildanimation(anim.dir, fn; kw...)
mov(anim::Animation, fn = movfn(); kw...) = buildanimation(anim.dir, fn, false; kw...)
mp4(anim::Animation, fn = mp4fn(); kw...) = buildanimation(anim.dir, fn, false; kw...)

function buildanimation(animdir::AbstractString, fn::AbstractString,
                        is_animated_gif::Bool=true;
                        fps::Integer = 20, loop::Integer = 0,
                        variable_palette::Bool=false,
                        show_msg::Bool=true)
    fn = abspath(fn)

    if is_animated_gif
        if variable_palette
            # generate a colorpalette for each frame for highest quality, but larger filesize
            palette="palettegen=stats_mode=single[pal],[0:v][pal]paletteuse=new=1"
            run(`ffmpeg -v 0 -framerate $fps -loop $loop -i $(animdir)/%06d.png -lavfi "$palette" -y $fn`)
        else
            # generate a colorpalette first so ffmpeg does not have to guess it
            run(`ffmpeg -v 0 -i $(animdir)/%06d.png -vf "palettegen=stats_mode=diff" -y "$(animdir)/palette.bmp"`)
            # then apply the palette to get better results
            run(`ffmpeg -v 0 -framerate $fps -loop $loop -i $(animdir)/%06d.png -i "$(animdir)/palette.bmp" -lavfi "paletteuse=dither=sierra2_4a" -y $fn`)
        end
    else
        run(`ffmpeg -v 0 -framerate $fps -loop $loop -i $(animdir)/%06d.png -pix_fmt yuv420p -y $fn`)
    end

    show_msg && @info("Saved animation to ", fn)
    AnimatedGif(fn)
end

# write out html to view the gif... note the rand call which is a hack so the image doesn't get cached
function Base.show(io::IO, ::MIME"text/html", agif::AnimatedGif)
    ext = file_extension(agif.filename)
    write(io, if ext == "gif"
        "<img src=\"$(relpath(agif.filename))?$(rand())>\" />"
    elseif ext in ("mov", "mp4")
        "<video controls><source src=\"$(relpath(agif.filename))?$(rand())>\" type=\"video/$ext\"></video>"
    else
        error("Cannot show animation with extension $ext: $agif")
    end)
end

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

function _animate(forloop::Expr, args...; callgif = false)
  if forloop.head != :for
    error("@animate macro expects a for-block. got: $(forloop.head)")
  end

  # add the call to frame to the end of each iteration
  animsym = gensym("anim")
  countersym = gensym("counter")
  freqassert = :()
  block = forloop.args[2]

  # create filter
  n = length(args)
  filterexpr = if n == 0
    # no filter... every iteration gets a frame
    true

  elseif args[1] == :every
    # filter every `freq` frames (starting with the first frame)
    @assert n == 2
    freq = args[2]
    freqassert = :(@assert isa($freq, Integer) && $freq > 0)
    :(mod1($countersym, $freq) == 1)

  elseif args[1] == :when
    # filter on custom expression
    @assert n == 2
    args[2]

  else
    error("Unsupported animate filter: $args")
  end

  push!(block.args, :(if $filterexpr; frame($animsym); end))
  push!(block.args, :(global $countersym += 1))

  # add a final call to `gif(anim)`?
  retval = callgif ? :(gif($animsym)) : animsym

  # full expression:
  esc(quote
    $freqassert             # if filtering, check frequency is an Integer > 0
    $animsym = Animation()  # init animation object
    global $countersym = 1         # init iteration counter
    $forloop                # for loop, saving a frame after each iteration
    $retval                 # return the animation object, or the gif
  end)
end

"""
Builds an `Animation` using one frame per loop iteration, then create an animated GIF.

Example:

```
  p = plot(1)
  @gif for x=0:0.1:5
    push!(p, 1, sin(x))
  end
```
"""
macro gif(forloop::Expr, args...)
    _animate(forloop, args...; callgif = true)
end

"""
Collect one frame per for-block iteration and return an `Animation` object.

Example:

```
p = plot(1)
anim = @animate for x=0:0.1:5
    push!(p, 1, sin(x))
end
gif(anim)
```
"""
macro animate(forloop::Expr, args...)
    _animate(forloop, args...)
end


text:
using Test
using SeisModels

@testset "Basic properties" begin
    let m = SteppedLayeredModel(1, 2, [0.5, 1], [1, 2], [3, 4], [5, 6], false,
            [], [], [], [], [], true, [10, 20], [100, 200])
        @test !isanisotropic(m)
        @test surface_radius(m) == 1.0
        @test radius(m, 0.25) == 0.75
        @test depth(m, 0.25) == 0.75
        @test hasattenuation(m)
        @test hasdensity(m)
    end
    let a = 1, n = 2, arr = rand(2), m = SteppedLayeredModel(a, n, [0.005 + 0.99rand(), 1],
            arr, arr, [], true, arr, arr, arr, arr, arr, false, [], [])
        @test isanisotropic(m)
        @test !hasattenuation(m)
        @test !hasdensity(m)
    end
    @test hasreffrequency(PREM)
    @test reffrequency(PREM) == 1.0
    let m = PREMPolyModel([getfield(PREM, f)
                           for f in fieldnames(PREMPolyModel)][1:end-1]..., NaN)
        @test !hasreffrequency(m)
        @test_throws ArgumentError reffrequency(m)
    end
end

@testset "Model comparison" begin
    @test IASP91 == deepcopy(IASP91)
    @test PREM != AK135
    let prem′ = deepcopy(PREM)
        @test prem′ ≈ PREM
        prem′.vp .= 1.00000001 .*PREM.vp
        @test prem′ ≈ PREM
        prem′.vp .= 1.000001 .* PREM.vp
        @test !(prem′ ≈ PREM)
    end
end


text:
# Definitions borrow from https://github.com/FluxML/NNlib.jl (names used to
# standardize package). TODO: Decide how/if to incorporate NNlib depedency into
# McCormick.jl/EAGO.jl.
oftf(x, y) = oftype(float(x), y)
leakyrelu(x, a=oftf(x, 0.01)) = max(a * x, x)
swish1(x) = x * sigmoid(x)

"""
    AtomType
"""
@enum(AtomType, VAR_ATOM, PARAM_ATOM, CONST_ATOM, SELECT_ATOM,
                ABS, ABS2, INV, RAD2DEG, DEG2RAD, ONE, ZERO,
                MULT, DIV, PLUS, MINUS, POW, MIN, MAX, STEP, SIGN,
                LOG, LOG2, LOG10, LOG1P, EXP, EXP2, EXP10, EXPM1,
                SIN, COS, TAN, CSC, SEC, COT,
                ASIN, ACOS, ATAN, ACSC, ASEC, ACOT,
                SINH, COSH, TANH, CSCH, SECH, COTH,
                ASINH, ACOSH, ATANH, ACSCH, ASECH, ACOTH,
                ERF, ERFC, ERFINV, ERFCINV, SQRT, CBRT,
                # ML specific functions
                RELU, LEAKYRELU, SOFTPLUS, SOFTSIGN, GELU, SILU, SIGMOID,
                # EAGO modeling functions
                XLOGX, ARH,
                # EAGO bound enforcing functions
                POS, NEG, LOWER_BND, UPPER_BND, BND,
                # Atoms for nonexpression types
                USER, USERN, SUBEXPR)

# Define univariate functions that correspond to the AtomType
const UNIVARIATE_ATOM_DICT = Dict{AtomType, Symbol}()

UNIVARIATE_ATOM_DICT[PLUS]      = :+
UNIVARIATE_ATOM_DICT[MINUS]     = :-
UNIVARIATE_ATOM_DICT[INV]       = :inv
UNIVARIATE_ATOM_DICT[ONE]       = :one
UNIVARIATE_ATOM_DICT[ZERO]      = :zero
UNIVARIATE_ATOM_DICT[ABS]       = :abs
UNIVARIATE_ATOM_DICT[ABS2]      = :abs2
UNIVARIATE_ATOM_DICT[RAD2DEG]   = :rad2deg
UNIVARIATE_ATOM_DICT[DEG2RAD]   = :deg2rad
UNIVARIATE_ATOM_DICT[STEP]      = :step
UNIVARIATE_ATOM_DICT[SIGN]      = :sign

UNIVARIATE_ATOM_DICT[LOG]       = :log
UNIVARIATE_ATOM_DICT[LOG2]      = :log2
UNIVARIATE_ATOM_DICT[LOG10]     = :log10
UNIVARIATE_ATOM_DICT[LOG1P]     = :log1p
UNIVARIATE_ATOM_DICT[EXP]       = :exp
UNIVARIATE_ATOM_DICT[EXP2]      = :exp2
UNIVARIATE_ATOM_DICT[EXP10]     = :exp10
UNIVARIATE_ATOM_DICT[EXPM1]     = :expm1

UNIVARIATE_ATOM_DICT[SIN]       = :sin
UNIVARIATE_ATOM_DICT[COS]       = :cos
UNIVARIATE_ATOM_DICT[TAN]       = :tan
UNIVARIATE_ATOM_DICT[CSC]       = :csc
UNIVARIATE_ATOM_DICT[SEC]       = :sec
UNIVARIATE_ATOM_DICT[COT]       = :cot
UNIVARIATE_ATOM_DICT[ASIN]      = :asin
UNIVARIATE_ATOM_DICT[ACOS]      = :acos
UNIVARIATE_ATOM_DICT[ATAN]      = :atan
UNIVARIATE_ATOM_DICT[ACSC]      = :acsc
UNIVARIATE_ATOM_DICT[ASEC]      = :asec
UNIVARIATE_ATOM_DICT[ACOT]      = :acot
UNIVARIATE_ATOM_DICT[SINH]      = :sinh
UNIVARIATE_ATOM_DICT[COSH]      = :cosh
UNIVARIATE_ATOM_DICT[TANH]      = :tanh
UNIVARIATE_ATOM_DICT[CSCH]      = :csch
UNIVARIATE_ATOM_DICT[SECH]      = :sech
UNIVARIATE_ATOM_DICT[COTH]      = :coth
UNIVARIATE_ATOM_DICT[ASINH]     = :asinh
UNIVARIATE_ATOM_DICT[ACOSH]     = :acosh
UNIVARIATE_ATOM_DICT[ATANH]     = :atanh
UNIVARIATE_ATOM_DICT[ACSCH]     = :acsch
UNIVARIATE_ATOM_DICT[ASECH]     = :asech
UNIVARIATE_ATOM_DICT[ACOTH]     = :acoth

UNIVARIATE_ATOM_DICT[ERF]       = :erf
UNIVARIATE_ATOM_DICT[ERFC]      = :erfc
UNIVARIATE_ATOM_DICT[ERFINV]    = :erfinv
UNIVARIATE_ATOM_DICT[ERFCINV]   = :erfcinv

UNIVARIATE_ATOM_DICT[SQRT]      = :sqrt
UNIVARIATE_ATOM_DICT[CBRT]      = :cbrt

UNIVARIATE_ATOM_DICT[POS]       = :positive
UNIVARIATE_ATOM_DICT[NEG]       = :negative
UNIVARIATE_ATOM_DICT[USER]      = :user

UNIVARIATE_ATOM_DICT[RELU]       = :relu
UNIVARIATE_ATOM_DICT[LEAKYRELU]  = :leakyrelu
UNIVARIATE_ATOM_DICT[SIGMOID]    = :sigmoid
UNIVARIATE_ATOM_DICT[GELU]       = :gelu
UNIVARIATE_ATOM_DICT[SOFTPLUS]   = :softplus
UNIVARIATE_ATOM_DICT[SOFTSIGN]   = :softsign
UNIVARIATE_ATOM_DICT[SILU]       = :swish1

UNIVARIATE_ATOM_DICT[XLOGX]      = :xlogx

# Define bivariate functions that correspond to the AtomType
const BIVARIATE_ATOM_DICT = Dict{AtomType, Symbol}()
BIVARIATE_ATOM_DICT[DIV]       = :/
BIVARIATE_ATOM_DICT[POW]       = :^
BIVARIATE_ATOM_DICT[LOWER_BND] = :lower_bnd
BIVARIATE_ATOM_DICT[UPPER_BND] = :upper_bnd
BIVARIATE_ATOM_DICT[ARH]       = :arh
BIVARIATE_ATOM_DICT[MINUS]     = :-
BIVARIATE_ATOM_DICT[MIN]       = :min
BIVARIATE_ATOM_DICT[MAX]       = :max
BIVARIATE_ATOM_DICT[PLUS]      = :+
BIVARIATE_ATOM_DICT[MULT]      = :*
#BIVARIATE_ATOM_DICT[EXPAX]     = :lower_bnd
#BIVARIATE_ATOM_DICT[EXPXY]     = :lower_bnd

# Define n-arity functions that correspond to the AtomType
const NARITY_ATOM_DICT = Dict{AtomType, Symbol}()
NARITY_ATOM_DICT[MIN]  = :min
NARITY_ATOM_DICT[MAX]  = :max
NARITY_ATOM_DICT[PLUS] = :+
NARITY_ATOM_DICT[MULT] = :*
NARITY_ATOM_DICT[USERN]   = :usern

const ALL_ATOM_DICT = Dict{AtomType, Symbol}()
foreach(x -> setindex!(ALL_ATOM_DICT, x[2], x[1]), UNIVARIATE_ATOM_DICT)
foreach(x -> setindex!(ALL_ATOM_DICT, x[2], x[1]), BIVARIATE_ATOM_DICT)
foreach(x -> setindex!(ALL_ATOM_DICT, x[2], x[1]), NARITY_ATOM_DICT)

# A functions that may be 1 to n-arity functions that correspond to the AtomType
ALL_ATOM_DICT[BND]       = :bnd
#ATM_EVAL[QUAD1]     = :quad1
#ATM_EVAL[QUAD2]     = :quad2

# List of keys only
const UNIVARIATE_ATOM_TYPES = AtomType[k for k in keys(UNIVARIATE_ATOM_DICT)]
const BIVARIATE_ATOM_TYPES = AtomType[k for k in keys(BIVARIATE_ATOM_DICT)]
const NARITY_ATOM_TYPES = AtomType[k for k in keys(NARITY_ATOM_DICT)]
const ALL_ATOM_TYPES = AtomType[k for k in keys(ALL_ATOM_DICT)]

# Reverse lookup dicts
const REV_UNIVARIATE_ATOM_DICT = Dict(UNIVARIATE_ATOM_DICT[k] => k for k in keys(UNIVARIATE_ATOM_DICT))
const REV_BIVARIATE_ATOM_DICT = Dict(BIVARIATE_ATOM_DICT[k] => k for k in keys(BIVARIATE_ATOM_DICT))
const REV_NARITY_ATOM_DICT = Dict(NARITY_ATOM_DICT[k] => k for k in keys(NARITY_ATOM_DICT))

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
struct TriangleMesh_2D <: LinearUnstructuredMesh_2D
    name::String
    points::Vector{Point_2D}
    edges::Vector{SVector{2, UInt32}}
    materialized_edges::Vector{LineSegment_2D}
    faces::Vector{SVector{3, UInt32}}
    materialized_faces::Vector{Triangle_2D}
    edge_face_connectivity::Vector{SVector{2, UInt32}}
    face_edge_connectivity::Vector{SVector{3, UInt32}}
    boundary_edges::Vector{Vector{UInt32}}
    face_sets::Dict{String, Set{UInt32}}
end

function TriangleMesh_2D(;
        name::String = "DefaultMeshName",
        points::Vector{Point_2D} = Point_2D[],
        edges::Vector{SVector{2, UInt32}} = SVector{2, UInt32}[],
        materialized_edges::Vector{LineSegment_2D} = LineSegment_2D[],
        faces::Vector{SVector{3, UInt32}} = SVector{3, U}[],
        materialized_faces::Vector{Triangle_2D} = Triangle_2D[],
        edge_face_connectivity::Vector{SVector{2, UInt32}} = SVector{2, UInt32}[],
        face_edge_connectivity ::Vector{SVector{3, UInt32}} = SVector{3, UInt32}[],
        boundary_edges::Vector{Vector{UInt32}} = Vector{UInt32}[],
        face_sets::Dict{String, Set{UInt32}} = Dict{String, Set{UInt32}}()
    )
        return TriangleMesh_2D(name,
                               points,
                               edges,
                               materialized_edges,
                               faces,
                               materialized_faces,
                               edge_face_connectivity,
                               face_edge_connectivity,
                               boundary_edges,
                               face_sets,
                              )
end

Base.broadcastable(mesh::TriangleMesh_2D) = Ref(mesh)

struct QuadrilateralMesh_2D <: LinearUnstructuredMesh_2D
    name::String
    points::Vector{Point_2D}
    edges::Vector{SVector{2, UInt32}}
    materialized_edges::Vector{LineSegment_2D}
    faces::Vector{SVector{4, UInt32}}
    materialized_faces::Vector{Quadrilateral_2D}
    edge_face_connectivity::Vector{SVector{2, UInt32}}
    face_edge_connectivity::Vector{SVector{4, UInt32}}
    boundary_edges::Vector{Vector{UInt32}}
    face_sets::Dict{String, Set{UInt32}}
end

function QuadrilateralMesh_2D(;
        name::String = "DefaultMeshName",
        points::Vector{Point_2D} = Point_2D[],
        edges::Vector{SVector{2, UInt32}} = SVector{2, UInt32}[],
        materialized_edges::Vector{LineSegment_2D} = LineSegment_2D[],
        faces::Vector{SVector{4, UInt32}} = SVector{4, U}[],
        materialized_faces::Vector{Quadrilateral_2D} = Quadrilateral_2D[],
        edge_face_connectivity::Vector{SVector{2, UInt32}} = SVector{2, UInt32}[],
        face_edge_connectivity ::Vector{SVector{4, UInt32}} = SVector{4, UInt32}[],
        boundary_edges::Vector{Vector{UInt32}} = Vector{UInt32}[],
        face_sets::Dict{String, Set{UInt32}} = Dict{String, Set{UInt32}}()
    )
        return QuadrilateralMesh_2D(name,
                                    points,
                                    edges,
                                    materialized_edges,
                                    faces,
                                    materialized_faces,
                                    edge_face_connectivity,
                                    face_edge_connectivity,
                                    boundary_edges,
                                    face_sets,
                                   )
end

Base.broadcastable(mesh::QuadrilateralMesh_2D) = Ref(mesh)

# Axis-aligned bounding box, in 2d a rectangle.
function boundingbox(mesh::LinearUnstructuredMesh_2D; boundary_shape::String="Unknown")
    # If the mesh does not have any quadratic faces, the bounding_box may be determined
    # entirely from the points.
    nsides = length(mesh.boundary_edges)
    if nsides !== 0
        boundary_edge_IDs = reduce(vcat, mesh.boundary_edges)
        point_IDs = reduce(vcat, mesh.edges[boundary_edge_IDs])
        return boundingbox(mesh.points[point_IDs])
    else
        return boundingbox(mesh.points)
    end
end

# A vector of SVectors, denoting the edge ID each face is connected to.
function face_edge_connectivity(mesh::QuadrilateralMesh_2D)
    if length(mesh.edges) === 0
        @error "Mesh does not have edges!"
    end
    # A vector of MVectors of zeros for each face
    # Each MVector is the length of the number of edges
    face_edge = [MVector{4, UInt32}(0, 0, 0, 0) for _ in eachindex(mesh.faces)]
    # For each face in the mesh, generate the edges.
    # Search for the index of the edge in the mesh.edges vector
    # Insert the index of the edge into the face_edge connectivity vector
    for i in eachindex(mesh.faces)
        for (j, edge) in enumerate(edges(mesh.faces[i]))
            face_edge[i][j] = searchsortedfirst(mesh.edges, SVector(edge.data))
        end
        if any(x->x === 0x00000000, face_edge[i])
            @error "Could not determine the face/edge connectivity of face $i"
        end
    end
    return [SVector(sort(conn).data) for conn in face_edge]
end

# A vector of SVectors, denoting the edge ID each face is connected to.
function face_edge_connectivity(mesh::TriangleMesh_2D)
    if length(mesh.edges) === 0
        @error "Mesh does not have edges!"
    end
    # A vector of MVectors of zeros for each face
    # Each MVector is the length of the number of edges
    face_edge = [MVector{3, UInt32}(0, 0, 0) for _ in eachindex(mesh.faces)]
    # For each face in the mesh, generate the edges.
    # Search for the index of the edge in the mesh.edges vector
    # Insert the index of the edge into the face_edge connectivity vector
    for i in eachindex(mesh.faces)
        for (j, edge) in enumerate(edges(mesh.faces[i]))
            face_edge[i][j] = searchsortedfirst(mesh.edges, SVector(edge.data))
        end
        if any(x->x === 0x00000000, face_edge[i])
            @error "Could not determine the face/edge connectivity of face $i"
        end
    end
    return [SVector(sort(conn).data) for conn in face_edge]
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
# Here we will define iteration (start,done,next) and indexing
# (getindex,setindex!,endof) interfaces for generic Value types.

# Julia handles access to AbstractArray, Associative, and Tuple
# subtypes using getindex:

# getindex(obj, key...) => value
# setindex!(obj, val, key...) => obj

# We do not allow overwriting, so setindex! for Values not allowed:

setindex!(x::Value,i...)=error("Overwriting operations currently not supported.")

# We handle the containers by overloading getindex:

@primitive  getindex(x,i...),dxi,xi  ungetindex(dxi,x,i)
getindex{T<:Grad}(::Type{T},o...)=nothing # Only the first arg has gradient
addtest(getindex, rand(2), 1)
addtest(getindex, rand(3), 2:3)
addtest(getindex, rand(2,2), 1, 2)
addtest(getindex, rand(3,3), 1:2, 3)

# Gradient of getindex: If xi=getindex(x,i...) and we receive dxi,
# ungetindex creates dx representing zeros similar to x, with only
# dx[i] set to dxi.  We use the sparse container OneHot for
# efficiency.
# x -> getindex -> xi -> grad -> dxi -> ungetindex -> dx -> grad -> ddx -> getindex -> ddxi

ungetindex(dxi,x,i)=OneHot(dxi,x,i)

# For higher order derivatives, the operation of ungetindex might be
# recorded and differentiated, so it must be a primitive.  It is only
# differentiable wrt its first arg.  The following methods cover
# (a,a), (v,a), (g,a), (g1,a) argtypes.

ungetindex_r = recorder(ungetindex)
ungetindex(dxi::Value,x,i)=ungetindex_r(dxi,x,i)
ungetindex(::Type{Grad{1}},ddx,dx,dxi,x,i)=getindex(ddx,i...)
ungetindex(::Type,o...)=nothing

# It should unbox its arguments, but it only needs to record if the
# first argument is boxed.  We'll have to define this manually.  To
# unbox the second arg and resolve ambiguity the following methods
# cover (a,v), (v,v), (g,v), (g1,v).

ungetindex(dxi::Value,x::Value,i)=ungetindex(dxi,getval(x),getval(i))
ungetindex(dxi,x::Value,i)=ungetindex(dxi,getval(x),getval(i))
ungetindex(::Type{Grad{1}},ddx::Value,dx,dxi,x,i)=getindex(ddx,getval(i)...)
ungetindex(::Type,ddx::Value,o...)=nothing

addtest(ungetindex, rand(), rand(2), (2,))
addtest(ungetindex, rand(2), rand(3), (2:3,))
addtest(ungetindex, rand(), rand(2,2), (1,2))
addtest(ungetindex, rand(2), rand(3,3), (1:2,3))

# For efficiency we use the following sparse container

if !isdefined(:OneHot)
immutable OneHot; value; container; index; end
end

Base.sum(b::OneHot)=sum(b.value)
Base.getindex(b::OneHot,i...)=getindex(full(b),i...)
Base.zeros(b::OneHot)=zeros(b.container)
Base.ones(b::OneHot)=ones(b.container)
Base.length(b::OneHot)=length(b.container)

function Base.full(b::OneHot)
    if isa(b.container,Tuple)
        c=zeroslike(collect(b.container))
        if isa(b.value,Tuple)
            setindex!(c,[b.value...],b.index...)
        else
            setindex!(c,b.value,b.index...)
        end
        return tuple(c...)
    else
        c=zeroslike(b.container)
        setindex!(c,b.value,b.index...)
        return c
    end
end

zeroslike{T<:Number}(a::AbstractArray{T})=zeros(a)
zeroslike(a::AbstractArray)=fill!(Array(Any,size(a)),nothing)
zeroslike(a::Associative)=similar(a)
zeroslike(o::OneHot)=zeros(o)

_dbg(x::OneHot)="OH$(id2(x))_$(_dbg(x.container))_$((x.index...))_$(x.value)"
Base.show(io::IO, n::OneHot)= print(io, _dbg(n))

# sum_outgrads needs to handle OneHot values:
sum_outgrads(a::OneHot,b::OneHot)=(if a.index==b.index; OneHot(a.container,sum_outgrads(a.value,b.value),a.index); else; sum_outgrads(full(a),b); end)
sum_outgrads(a::Value,b::OneHot)=error((:sum,a,b))
sum_outgrads(a::OneHot,b::Value)=error((:sum,a,b))
sum_outgrads(a::Void,b::OneHot)=full(b)
sum_outgrads(a::OneHot,b::Void)=error((:sum,a,b))
sum_outgrads(a::OneHot,b)=error((:sum,a,Any))
sum_outgrads(a::Tuple,b::OneHot)=(b=full(b);ntuple(length(a)) do i; sum_outgrads(a[i],b[i]); end)
sum_outgrads(a::AbstractArray,b::OneHot)=setindex!(a,sum_outgrads(getindex(a,b.index...),b.value),b.index...)
sum_outgrads(a::Associative,b::OneHot)=setindex!(a,sum_outgrads(get(a,b.index...,nothing),b.value),b.index...)

# Iteration is used in `for x in a` loops and for `(x,y)=a` multiple
# assignments.

# start(a) => initial state
# done(a,state) => whether state s is final
# next(a,state) => (element, nextState)

# start and done return state and bool, not differentiable.
start(a::Value)=start(a.value)
done(a::Value,i)=done(a.value,i)

# next returns a tuple with an element and needs to be defined for each iterable.
# Specific types need to define their own next methods for Values:

next(a::Value,i)=throw(MethodError(next,(a,i)))
next{T<:Array}(a::Value{T},i) = (a[i],i+1)
next{T<:Tuple}(a::Value{T},i) = (a[i],i+1)
next{T<:Number}(a::Value{T},i) = (a,true)
# This needs more work:
# next{T<:Base.ValueIterator}(a::Value{T},i) = (d=a.value.dict; (d.vals[i], skip_deleted(d,i+1)))

# Finally here are some common functions that do not return floats
# (e.g. length) or return constant outputs (e.g. zero).

interfaces1arg = [
:eltype,
:endof,
:isempty,
:length,
:ndims,
:one,
:ones,
:strides,
:zero,
:zeros,
]

for _f in interfaces1arg
    @eval @zerograd $_f(x)
end

interfacesNarg = [
:checkbounds,
:eachindex,
:isassigned,
:pointer,
:similar,
:size,
:stride,
]

for _f in interfacesNarg
    @eval @zerograd $_f(x,i...)
end

interfaces2arg = [
:(==),
:isequal,
:isless,
]

==(a::WeakRef,b::Value)=(a==b.value) # prevents clash with base.jl:68
==(a::Value,b::WeakRef)=(a.value==b) # prevents clash with base.jl:69

for _f in interfaces2arg
    @eval @zerograd $_f(a,b)
end

@primitive copy(x),dy dy
addtest(copy, rand(2))


text:
export AdjDxGeo, AdjDyGeo, AdjDxExp, AdjDξExp, AdjDxLog, AdjDyLog
export AdjDforwardLogs
@doc doc"""
    AdjDxGeo(M,x,y,t,η)
computes the adjoint of $D_xg(t;x,y)[\eta]$.

# See also
 [`DxGeo`](@ref), [`adjointJacobiField`](@ref)
"""
AdjDxGeo(M::mT,x::P,y::P,t::Float64,η::T) where {mT <: Manifold, P <: MPoint, T<: TVector}= adjointJacobiField(M,x,y,t,η,βDgx)
@doc doc"""
    AdjDyGeo(M,x,y,t,η)
computes the adjoint of $D_yg(t;x,y)[\eta]$.

# See also
 [`DyGeo`](@ref), [`adjointJacobiField`](@ref)
"""
AdjDyGeo(M::mT,x::P,y::P,t::Float64,η::T) where {mT <: Manifold, P <: MPoint, T<: TVector} = adjointJacobiField(M,y,x,1-t,η,βDgx)
@doc doc"""
    AdjDxExp(M,x,ξ,η)
computes the adjoint of $D_x\exp_x\xi[\eta]$.

# See also
 [`DxExp`](@ref), [`adjointJacobiField`](@ref)
"""
AdjDxExp(M::mT,x::P,ξ::T,η::T) where {mT <: Manifold, P <: MPoint, T<: TVector} = adjointJacobiField(M,x,exp(M,x,ξ),1.,η,βDexpx)
@doc doc"""
    AdjDξExp(M,x,ξ,η)
computes the adjoint of $D_\xi\exp_x\xi[\eta]$.
Note that $\xi\in T_\xi(T_x\mathcal M) = T_x\mathcal M$ is still a tangent vector.

# See also
 [`DξExp`](@ref), [`adjointJacobiField`](@ref)
"""
AdjDξExp(M::mT,x::P,ξ::T,η::T) where {mT <: Manifold, P <: MPoint, T<: TVector} = adjointJacobiField(M,x,exp(M,x,ξ),1.0,η,βDexpξ)
@doc doc"""
    AdjDxLog(M,x,y,η)
computes the adjoint of $D_xlog_xy[\eta]$.

# See also
 [`DxLog`](@ref), [`adjointJacobiField`](@ref)
"""
AdjDxLog(M::mT,x::P,y::P,η::T) where {mT <: Manifold, P <: MPoint, T<: TVector} = adjointJacobiField(M,x,y,0.,η,βDlogx)
@doc doc"""
    AdjDyLog(M,x,y,η)
computes the adjoint of $D_ylog_xy[\eta]$.

# See also
 [`DyLog`](@ref), [`adjointJacobiField`](@ref)
"""
AdjDyLog(M::mT,x::P,y::P,η::T) where {mT <: Manifold, P <: MPoint, T<: TVector} = adjointJacobiField(M,x,y,1.,η,βDlogy)
@doc doc"""
    ξ = AdjDforwardLogs(M,x,ν)

compute the adjoibnt differential of [`forwardLogs`](@ref) $F$ orrucirng,
in the power manifold array, the differential of the function

$F_i(x) = \sum_{j\in\mathcal I_i} \log_{x_i} x_j$

where $i$ runs over all indices of the [`Power`](@ref) manifold `M` and $\mathcal I_i$
denotes the forward neighbors of $i$
Let $n$ be the number dimensions of the [`Power`](@ref) manifold (i.e. `length(size(x)`)).
Then the input tangent vector lies on the manifold $\mathcal M' = \mathcal M^n$.

# Input
* `M`     – a [`Power`](@ref) manifold
* `x`     – a [`PowPoint`](@ref).
* `ν`     – a [`PowTVector`](@ref) from $T_X\mathcal M'$, where
  $X = (x,...,x)\in\mathcal M'$ is an $n$-fold copy of $x$ where $\mathcal N (x,...,x)N.

# Ouput
* ξ – resulting tangent vector in $T_x\mathcal M$ representing the adjoint
  differentials of the logs.
"""
function AdjDforwardLogs(M::Power,x::PowPoint,ν::PowTVector)::PowTVector
  sX = size(x)
  R = CartesianIndices(sX)
  d = length(sX)
  maxInd = [last(R).I...] # maxInd as Array
  N = Power(M.manifold,(sX...,d))
  ξ = zeroTVector(M,x)
  for i in R # iterate over all pixel
    for k in 1:d # for all direction combinations
        I = [i.I...] # array of index
        J = I .+ 1 .* (1:d .== k) #i + e_k is j
        if all( J .<= maxInd ) # is this neighbor in range?
            j = CartesianIndex{d}(J...) # neigbbor index as Cartesian Index
            ξ[i] += AdjDxLog(M.manifold,x[i],x[j],ν[i,k])
            ξ[j] += AdjDyLog(M.manifold,x[i],x[j],ν[i,k])
        end
    end # directions
  end # i in R
  return ξ
end


text:
using GPUCompiler

if !@isdefined(TestRuntime)
    include("../util.jl")
end

# create a SPIRV-based test compiler, and generate reflection methods for it

for method in (:code_typed, :code_warntype, :code_llvm, :code_native)
    # only code_typed doesn't take a io argument
    args = method == :code_typed ? (:job,) : (:io, :job)
    spirv_method = Symbol("spirv_$(method)")

    @eval begin
        function $spirv_method(io::IO, @nospecialize(func), @nospecialize(types);
                               kernel::Bool=false, kwargs...)
            source = FunctionSpec(func, Base.to_tuple_type(types), kernel)
            target = SPIRVCompilerTarget()
            params = TestCompilerParams()
            job = CompilerJob(target, source, params)
            GPUCompiler.$method($(args...); kwargs...)
        end
        $spirv_method(@nospecialize(func), @nospecialize(types); kwargs...) =
            $spirv_method(stdout, func, types; kwargs...)
    end
end


text:
typealias Verb UInt8

const GET    = 1 << 0 |> Verb
const PUT    = 1 << 1 |> Verb
const POST   = 1 << 2 |> Verb
const DELETE = 1 << 3 |> Verb
const PATCH  = 1 << 4 |> Verb

const STI = Dict(
  "GET"    =>  GET,
  "PUT"    =>  PUT,
  "POST"   =>  POST,
  "DELETE" =>  DELETE,
  "PATCH"  =>  PATCH,
  "get"    =>  GET,
  "put"    =>  PUT,
  "post"   =>  POST,
  "delete" =>  DELETE,
  "patch"  =>  PATCH
)

"""
# proper_method

for cases when an endpoint can many different verbs
first check if there is an exact match simpler
"""
function proper_method(verb, map)
  in(verb, map |> keys) && return map[verb]

  for v in keys(map)
    verb & v > 0 && return map[v]
  end

  throw("No key $verb in map.")
end


text:
########################################
## File Name: type_test.jl
## Author: Haruki Nishimura (hnishimura@stanford.edu)
## Date Created: 2020/05/12
## Description: Test code for src/state_types.jl, src/belief_types.jl,
## src/control_types.jl
########################################

using LinearAlgebra
using Distributions

@testset "Type Definition Test" begin

# Belief Types
b1 = BelMvNormal(1., zeros(2), Matrix(1.0I, 2, 2));
@test b1.t == 1.0
@test b1.params == [0.0, 0.0, 1.0, 0.0, 0.0, 1.0]
@test b1.dim == 2

d1 = MvNormal(b1)
@test d1.μ == zeros(2)
@test Matrix(d1.Σ) == Matrix(1.0I, 2, 2)

b2 = VecBelMvNormal([b1 for ii = 1:10])
@test length(b2) == 10

d2 = MvNormal.(b2);
@test length(d2) == 10

# State Types
p1 = PhysPos(1., zeros(2), 2)
p2 = PhysPos(1., zeros(2))
@test p2.t == p1.t
@test p2.pos == p1.pos
@test p2.dim == p1.dim
@test vec(p2) == zeros(2)

p3 = PhysManipulate2D(1., zeros(2), pi/2, zeros(2), 0., 1., 1., ones(2), 0.3)
@test vec(p3) == [0.0, 0.0, pi/2, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.3]

p4 = PhysManipulate2D(1., vec(p3))
@test p4.t == 1.0
@test p4.pos == [0.0, 0.0]
@test p4.θ == pi/2
@test p4.vel == [0.0, 0.0]
@test p4.ω == 0.
@test p4.m == 1.
@test p4.J == 1.
@test p4.r == [1.0, 1.0]
@test p4.μ == 0.3
@test p4.dim == 2

s1 = AugState(p1, b2)
@test s1.t == p1.t
@test s1.p == p1
@test s1.b == b2

# Control Types
c1 = PosControl(0., zeros(2))
@test c1.t == 0.
@test c1.vel == [0.0, 0.0]
@test c1.dim == 2

c2 = PosControl(0., 0.)
@test c2.t == 0.
@test c2.vel == [0.0]
@test c2.dim == 1

c3 = MControl2D(0., 1., 2., 3.)
@test c3 == MControl2D(0., [1., 2., 3.])
@test vec(c3) == [1., 2., 3.]

end


text:
module CPython

using Libdl, Dates, Compat, Requires
import ..Python:
    CONFIG,
    isnull,
    ism1,
    PYERR,
    NOTIMPLEMENTED,
    _typeintersect,
    tryconvert,
    ispyreftype,
    ispyref,
    pyptr,
    putresult,
    takeresult,
    CACHE,
    Python
using Base: @kwdef
using UnsafePointers: UnsafePtr

pyglobal(name) = dlsym(CONFIG.libptr, name)
pyglobal(r::Ref{Ptr{T}}, name) where {T} = begin
    p = r[]
    if isnull(p)
        p = r[] = Ptr{T}(pyglobal(name))
    end
    p
end
pyloadglobal(r::Ref{Ptr{T}}, name) where {T} = begin
    p = r[]
    if isnull(p)
        p = r[] = unsafe_load(Ptr{Ptr{T}}(pyglobal(name)))
    end
    p
end

macro cdef(name, rettype, argtypes)
    name isa QuoteNode && name.value isa Symbol || error("name must be a symbol, got $name")
    jname = esc(name.value)
    refname = esc(Symbol(name.value, :__ref))
    name = esc(name)
    rettype = esc(rettype)
    argtypes isa Expr && argtypes.head == :tuple ||
        error("argtypes must be a tuple, got $argtypes")
    nargs = length(argtypes.args)
    argtypes = esc(argtypes)
    args = [gensym() for i = 1:nargs]
    quote
        const $refname = Ref(C_NULL)
        $jname($(args...)) =
            ccall(pyglobal($refname, $name), $rettype, $argtypes, $(args...))
    end
end

include("consts.jl")
include("fundamentals.jl")
include("none.jl")
include("object.jl")
include("sequence.jl")
include("mapping.jl")
include("method.jl")
include("str.jl")
include("bytes.jl")
include("tuple.jl")
include("type.jl")
include("number.jl")
include("iter.jl")
include("bool.jl")
include("int.jl")
include("float.jl")
include("complex.jl")
include("list.jl")
include("dict.jl")
include("set.jl")
include("buffer.jl")
include("stdlib.jl")
include("collections.jl")
include("io.jl")
include("range.jl")
include("ctypes.jl")
include("numpy.jl")
include("slice.jl")
include("fraction.jl")
include("datetime.jl")
include("newtype.jl")
include("juliaerror.jl")
include("juliabase.jl")
include("juliaraw.jl")
include("juliaany.jl")
include("juliaiterator.jl")
include("juliatype.jl")
include("juliadict.jl")
include("juliaset.jl")
include("juliaarray.jl")
include("juliavector.jl")
include("juliamodule.jl")
include("julianumber.jl")
include("juliaio.jl")
include("as.jl")
include("arg.jl")

@init begin
    PyObject_TryConvert_AddRules(
        "builtins.NoneType",
        [
            (Nothing, PyNone_TryConvertRule_nothing, 100),
            (Missing, PyNone_TryConvertRule_missing),
        ],
    )
    PyObject_TryConvert_AddRules("builtins.bool", [(Bool, PyBool_TryConvertRule_bool, 100)])
    PyObject_TryConvert_AddRules(
        "numbers.Integral",
        [
            (Integer, PyLongable_TryConvertRule_integer, 100),
            (Rational, PyLongable_TryConvertRule_tryconvert),
            (Real, PyLongable_TryConvertRule_tryconvert),
            (Number, PyLongable_TryConvertRule_tryconvert),
            (Any, PyLongable_TryConvertRule_tryconvert),
        ],
    )
    PyObject_TryConvert_AddRules(
        "builtins.float",
        [(Float64, PyFloatable_TryConvertRule_convert, 100)],
    )
    PyObject_TryConvert_AddRules(
        "numbers.Real",
        [
            (Float64, PyFloatable_TryConvertRule_convert),
            (BigFloat, PyFloatable_TryConvertRule_convert),
            (Float32, PyFloatable_TryConvertRule_convert),
            (Float16, PyFloatable_TryConvertRule_convert),
            (AbstractFloat, PyFloatable_TryConvertRule_tryconvert),
            (Real, PyFloatable_TryConvertRule_tryconvert),
            (Number, PyFloatable_TryConvertRule_tryconvert),
        ],
    )
    PyObject_TryConvert_AddRules(
        "builtins.complex",
        [(Complex{Float64}, PyComplexable_TryConvertRule_convert, 100)],
    )
    PyObject_TryConvert_AddRules(
        "numbers.Complex",
        [
            (Complex{Float64}, PyComplexable_TryConvertRule_convert),
            (Complex{BigFloat}, PyComplexable_TryConvertRule_convert),
            (Complex{Float32}, PyComplexable_TryConvertRule_convert),
            (Complex{Float16}, PyComplexable_TryConvertRule_convert),
            (Complex{T} where {T<:AbstractFloat}, PyComplexable_TryConvertRule_tryconvert),
            (Complex{T} where {T<:Real}, PyComplexable_TryConvertRule_tryconvert),
            (Number, PyComplexable_TryConvertRule_tryconvert),
        ],
    )
    PyObject_TryConvert_AddRules(
        "builtins.bytes",
        [
            (Vector{UInt8}, PyBytes_TryConvertRule_vector),
            (Vector{Int8}, PyBytes_TryConvertRule_vector),
            (String, PyBytes_TryConvertRule_string),
        ],
    )
    PyObject_TryConvert_AddRules(
        "builtins.str",
        [
            (String, PyUnicode_TryConvertRule_string, 100),
            (Symbol, PyUnicode_TryConvertRule_symbol),
            (Char, PyUnicode_TryConvertRule_char),
            (Vector{UInt8}, PyUnicode_TryConvertRule_vector),
            (Vector{Int8}, PyUnicode_TryConvertRule_vector),
        ],
    )
    PyObject_TryConvert_AddRules(
        "builtins.tuple",
        [(Tuple, PyIterable_ConvertRule_tuple, 100)],
    )
    PyObject_TryConvert_AddRules(
        "builtins.range",
        [
            (
                StepRange{T,S} where {T<:Integer,S<:Integer},
                PyRange_TryConvertRule_steprange,
                100,
            ),
            (UnitRange{T} where {T<:Integer}, PyRange_TryConvertRule_unitrange),
        ],
    )
    PyObject_TryConvert_AddRules(
        "collections.abc.Iterable",
        [
            (Vector, PyIterable_ConvertRule_vecorset),
            (Set, PyIterable_ConvertRule_vecorset),
            (Tuple, PyIterable_ConvertRule_tuple),
            (NamedTuple, PyIterable_ConvertRule_namedtuple),
            (Pair, PyIterable_ConvertRule_pair),
        ],
    )
    PyObject_TryConvert_AddRules(
        "collections.abc.Sequence",
        [(Vector, PyIterable_ConvertRule_vecorset)],
    )
    PyObject_TryConvert_AddRules(
        "collections.abc.Set",
        [(Set, PyIterable_ConvertRule_vecorset)],
    )
    PyObject_TryConvert_AddRules(
        "collections.abc.Mapping",
        [(Dict, PyMapping_ConvertRule_dict)],
    )
    PyObject_TryConvert_AddRules("datetime.time", [(Time, PyTime_TryConvertRule_time, 100)])
    PyObject_TryConvert_AddRules("datetime.date", [(Date, PyDate_TryConvertRule_date, 100)])
    PyObject_TryConvert_AddRules(
        "datetime.datetime",
        [(DateTime, PyDateTime_TryConvertRule_datetime, 100)],
    )
    PyObject_TryConvert_AddRules(
        "datetime.timedelta",
        [
            (Dates.Period, PyTimeDelta_TryConvertRule_period, 100),
            (Dates.CompoundPeriod, PyTimeDelta_TryConvertRule_compoundperiod),
        ],
    )
    PyObject_TryConvert_AddRules(
        "juliaaa.ValueBase",
        [(Any, PyJuliaValue_TryConvert_any, 1000)],
    )
    PyObject_TryConvert_AddExtraTypes([
        PyIterableABC_Type,
        PyCallableABC_Type,
        PySequenceABC_Type,
        PyMappingABC_Type,
        PySetABC_Type,
        PyNumberABC_Type,
        PyComplexABC_Type,
        PyRealABC_Type,
        PyRationalABC_Type,
        PyIntegralABC_Type,
    ])

    ### ctypes
    for (p, T) in [
        ("char", Cchar),
        ("wchar", Cwchar_t),
        ("byte", Cchar),
        ("ubyte", Cuchar),
        ("short", Cshort),
        ("ushort", Cushort),
        ("int", Cint),
        ("uint", Cuint),
        ("long", Clong),
        ("ulong", Culong),
        ("longlong", Clonglong),
        ("ulonglong", Culonglong),
        ("size_t", Csize_t),
        ("ssize_t", Cssize_t),
        ("float", Cfloat),
        ("double", Cdouble), #=("longdouble", ???),=#
        ("char_p", Ptr{Cchar}),
        ("wchar_p", Ptr{Cwchar_t}),
        ("void_p", Ptr{Cvoid}),
    ]
        isptr = occursin("_p", p)
        isfloat = occursin("float", p) || occursin("double", p)
        isint = !(isfloat || isptr)
        isreal = isint || isfloat
        PyObject_TryConvert_AddRules(
            "ctypes.c_$p",
            [
                (
                    p == "char_p" ? Cstring : p == "wchar_p" ? Cwstring : Union{},
                    PySimpleCData_TryConvert_value{T,false}(),
                ),
                (T, PySimpleCData_TryConvert_value{T,false}()),
                (isint ? Integer : Union{}, PySimpleCData_TryConvert_value{T,true}()),
                (isint ? Rational : Union{}, PySimpleCData_TryConvert_value{T,true}()),
                (isreal ? Float64 : Union{}, PySimpleCData_TryConvert_value{T,false}()),
                (isreal ? BigFloat : Union{}, PySimpleCData_TryConvert_value{T,false}()),
                (isreal ? Float32 : Union{}, PySimpleCData_TryConvert_value{T,false}()),
                (isreal ? Float16 : Union{}, PySimpleCData_TryConvert_value{T,false}()),
                (
                    isreal ? AbstractFloat : Union{},
                    PySimpleCData_TryConvert_value{T,true}(),
                ),
                (isreal ? Real : Union{}, PySimpleCData_TryConvert_value{T,true}()),
                (isreal ? Number : Union{}, PySimpleCData_TryConvert_value{T,true}()),
                (isptr ? Ptr : Union{}, PySimpleCData_TryConvert_value{T,false}()),
                (Any, PySimpleCData_TryConvert_value{T,true}()),
            ],
        )
    end

    ### numpy
    # TODO: Compound types
    # TODO: datetime64, timedelta64
    for (p, T) in [
        ("int8", Int8),
        ("int16", Int16),
        ("int32", Int32),
        ("int64", Int64),
        ("int128", Int128),
        ("uint8", UInt8),
        ("uint16", UInt16),
        ("uint32", UInt32),
        ("uint64", UInt64),
        ("uint128", UInt128),
        ("float16", Float16),
        ("float32", Float32),
        ("float64", Float64),
        ("complex32", Complex{Float16}),
        ("complex64", Complex{Float32}),
        ("complex128", Complex{Float64}),
    ]
        isint = occursin("int", p)
        isfloat = occursin("float", p)
        iscomplex = occursin("complex", p)
        isreal = isint || isfloat
        PyObject_TryConvert_AddRules(
            "numpy.$p",
            [
                (T, PyNumpySimpleData_TryConvert_value{T,false}(), 100),
                (isint ? Integer : Union{}, PyNumpySimpleData_TryConvert_value{T,true}()),
                (isint ? Rational : Union{}, PyNumpySimpleData_TryConvert_value{T,true}()),
                (isreal ? Float64 : Union{}, PyNumpySimpleData_TryConvert_value{T,false}()),
                (
                    isreal ? BigFloat : Union{},
                    PyNumpySimpleData_TryConvert_value{T,false}(),
                ),
                (isreal ? Float32 : Union{}, PyNumpySimpleData_TryConvert_value{T,false}()),
                (isreal ? Float16 : Union{}, PyNumpySimpleData_TryConvert_value{T,false}()),
                (
                    isreal ? AbstractFloat : Union{},
                    PyNumpySimpleData_TryConvert_value{T,true}(),
                ),
                (isreal ? Real : Union{}, PyNumpySimpleData_TryConvert_value{T,true}()),
                (
                    iscomplex ? Complex{Float64} : Union{},
                    PyNumpySimpleData_TryConvert_value{T,false}(),
                ),
                (
                    iscomplex ? Complex{BigFloat} : Union{},
                    PyNumpySimpleData_TryConvert_value{T,false}(),
                ),
                (
                    iscomplex ? Complex{Float32} : Union{},
                    PyNumpySimpleData_TryConvert_value{T,false}(),
                ),
                (
                    iscomplex ? Complex{Float16} : Union{},
                    PyNumpySimpleData_TryConvert_value{T,false}(),
                ),
                (
                    iscomplex ? (Complex{T} where {T<:AbstractFloat}) : Union{},
                    PyNumpySimpleData_TryConvert_value{T,true}(),
                ),
                (
                    iscomplex ? (Complex{T} where {T<:Real}) : Union{},
                    PyNumpySimpleData_TryConvert_value{T,true}(),
                ),
                (Number, PyNumpySimpleData_TryConvert_value{T,true}()),
                (Any, PyNumpySimpleData_TryConvert_value{T,true}()),
            ],
        )
    end
end

end


text:
"""
    ODESolvers

Ordinary differential equation solvers
"""
module ODESolvers

using KernelAbstractions
using KernelAbstractions.Extras: @unroll
using StaticArrays
using Requires
@init @require CUDAnative = "be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
    using .CUDAnative
end

using ..SpaceMethods
using ..LinearSolvers
using ..MPIStateArrays: device, realview

export solve!, updatedt!, gettime

abstract type AbstractODESolver end
"""
    gettime(solver::AbstractODESolver)

Returns the current simulation time of the ODE solver `solver`
"""
gettime(solver::AbstractODESolver) = solver.t

"""
    getdt(solver::AbstractODESolver)

Returns the current simulation time step of the ODE solver `solver`
"""
getdt(solver::AbstractODESolver) = solver.dt

"""
    ODESolvers.general_dostep!(Q, solver::AbstractODESolver, p,
                               timeend::Real, adjustfinalstep::Bool)

Use the solver to step `Q` forward in time from the current time, to the time
`timeend`. If `adjustfinalstep == true` then `dt` is adjusted so that the step
does not take the solution beyond the `timeend`.
"""
function general_dostep!(
    Q,
    solver::AbstractODESolver,
    p,
    timeend::Real;
    adjustfinalstep::Bool,
)
    time, dt = solver.t, solver.dt
    final_step = false
    if adjustfinalstep && time + dt > timeend
        dt = timeend - time
        updatedt!(solver, dt)
        final_step = true
    end
    @assert dt > 0

    dostep!(Q, solver, p, time)

    if !final_step
        solver.t += dt
    else
        solver.t = timeend
    end
end

"""
    updatedt!(solver::AbstractODESolver, dt)

Change the time step size to `dt` for the ODE solver `solver`.
"""
updatedt!(solver::AbstractODESolver, dt) = (solver.dt = dt)

"""
    updatetime!(solver::AbstractODESolver, time)

Change the current time to `time` for the ODE solver `solver`.
"""
updatetime!(solver::AbstractODESolver, time) = (solver.t = time)

isadjustable(solver::AbstractODESolver) = true

# {{{ run!
"""
    solve!(Q, solver::AbstractODESolver; timeend,
           stopaftertimeend=true, numberofsteps, callbacks)

Solves an ODE using the `solver` starting from a state `Q`. The state `Q` is
updated inplace. The final time `timeend` or `numberofsteps` must be specified.

A series of optional callback functions can be specified using the tuple
`callbacks`; see [`GenericCallbacks`](@ref).
"""
function solve!(
    Q,
    solver::AbstractODESolver,
    p = nothing;
    timeend::Real = Inf,
    adjustfinalstep = true,
    numberofsteps::Integer = 0,
    callbacks = (),
)

    @assert isfinite(timeend) || numberofsteps > 0
    if adjustfinalstep && !isadjustable(solver)
        error("$solver does not support time step adjustments. Can only be used with `adjustfinalstep=false`.")
    end
    t0 = gettime(solver)

    # Loop through an initialize callbacks (if they need it)
    foreach(callbacks) do cb
        try
            cb(true)
        catch
        end
    end

    step = 0
    time = t0
    while time < timeend
        step += 1

        time = general_dostep!(
            Q,
            solver,
            p,
            timeend;
            adjustfinalstep = adjustfinalstep,
        )

        # FIXME: Determine better way to handle postcallback behavior
        # Current behavior:
        #   retval = 1 exit after all callbacks
        #   retval = 2 exit immediately
        retval = 0
        for (i, cb) in enumerate(callbacks)
            # FIXME: Consider whether callbacks need anything, or if function closure
            #        can be used for everything
            thisretval = cb()
            thisretval = (thisretval == nothing) ? 0 : thisretval
            !(thisretval in (0, 1, 2)) &&
            error("callback #$(i) returned invalid value. It should return either:
                  `nothing` (continue time stepping)
                  `0`       (continue time stepping)
                  `1`       (stop time stepping after all callbacks)
                  `2`       (stop time stepping immediately)")
            retval = max(thisretval, retval)
            retval == 2 && return gettime(solver)
        end
        retval == 1 && return gettime(solver)

        # Figure out if we should stop
        if numberofsteps == step
            return gettime(solver)
        end
    end
    gettime(solver)
end
# }}}

include("LowStorageRungeKuttaMethod.jl")
include("StrongStabilityPreservingRungeKuttaMethod.jl")
include("AdditiveRungeKuttaMethod.jl")
include("MultirateInfinitesimalStepMethod.jl")
include("MultirateRungeKuttaMethod.jl")

end # module

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
using Test, ArrowTypes, UUIDs

include("tests.jl")


text:
########## Fake type-system

# Used to label all objects
struct VarRef
    parent::Union{VarRef,Nothing}
    name::Symbol
end
VarRef(m::Module) = VarRef((parentmodule(m) == Main || parentmodule(m) == m) ? nothing : VarRef(parentmodule(m)), nameof(m))

# These mirror Julia types (w/o the Fake prefix)
struct FakeTypeName
    name::VarRef
    parameters::Vector{Any}
end

function FakeTypeName(@nospecialize(x); justname=false)
    if x isa DataType
        xname = x.name
        xnamename = xname.name
        if justname
            FakeTypeName(VarRef(VarRef(x.name.module), x.name.name), [])
        else
            # FakeTypeName(VarRef(VarRef(x.name.module), x.name.name), _parameter.(x.parameters))
            ft = FakeTypeName(VarRef(VarRef(x.name.module), x.name.name), [])
            for p in x.parameters
                push!(ft.parameters, _parameter(p))
            end
            ft
        end
    elseif x isa Union
        FakeUnion(x)
    elseif x isa UnionAll
        FakeUnionAll(x)
    elseif x isa TypeVar
        FakeTypeVar(x)
    elseif x isa Core.TypeofBottom
        FakeTypeofBottom()
    else
        error((x, typeof(x)))
    end
end

struct FakeTypeofBottom end
struct FakeUnion
    a
    b
end
FakeUnion(u::Union) = FakeUnion(FakeTypeName(u.a, justname=true), FakeTypeName(u.b, justname=true))
struct FakeTypeVar
    name::Symbol
    lb
    ub
end
FakeTypeVar(tv::TypeVar) = FakeTypeVar(tv.name, FakeTypeName(tv.lb, justname=true), FakeTypeName(tv.ub, justname=true))
struct FakeUnionAll
    var::FakeTypeVar
    body::Any
end
FakeUnionAll(ua::UnionAll) = FakeUnionAll(FakeTypeVar(ua.var), FakeTypeName(ua.body, justname=true))

function _parameter(@nospecialize(p))
    if p isa Union{Int,Symbol,Bool,Char}
        p
    elseif !(p isa Type) && isbitstype(typeof(p))
        0
    elseif p isa Tuple
        _parameter.(p)
    else
        FakeTypeName(p, justname=true)
    end
end

Base.print(io::IO, vr::VarRef) = vr.parent === nothing ? print(io, vr.name) : print(io, vr.parent, ".", vr.name)
function Base.print(io::IO, tn::FakeTypeName)
    print(io, tn.name)
    if !isempty(tn.parameters)
        print(io, "{")
        for i = 1:length(tn.parameters)
            print(io, tn.parameters[i])
            i != length(tn.parameters) && print(io, ",")
        end
        print(io, "}")
    end
end
Base.print(io::IO, x::FakeUnionAll) = print(io, x.body, " where ", x.var)
function Base.print(io::IO, x::FakeUnion; inunion=false)
    !inunion && print(io,  "Union{")
    print(io, x.a, ",")
    if x.b isa FakeUnion
        print(io, x.b, inunion=true)
    else
        print(io, x.b, "}")
    end
end
function Base.print(io::IO, x::FakeTypeVar)
    if isfakebottom(x.lb)
        if isfakeany(x.ub)
            print(io, x.name)
        else
            print(io, x.name, "<:", x.ub)
        end
    elseif isfakeany(x.ub)
        print(io, x.lb, "<:", x.name)
    else
        print(io, x.lb, "<:", x.name, "<:", x.ub)
    end
end

isfakeany(t) = false
isfakeany(t::FakeTypeName) = isfakeany(t.name)
isfakeany(vr::VarRef) = vr.name === :Any && vr.parent isa VarRef && vr.parent.name === :Core && vr.parent.parent === nothing

isfakebottom(t) = false
isfakebottom(t::FakeTypeofBottom) = true

Base.:(==)(a::FakeTypeName, b::FakeTypeName) = a.name == b.name && a.parameters == b.parameters
Base.:(==)(a::VarRef, b::VarRef) = a.parent == b.parent && a.name == b.name
Base.:(==)(a::FakeTypeVar, b::FakeTypeVar) = a.lb == b.lb && a.name == b.name && a.ub == b.ub
Base.:(==)(a::FakeUnionAll, b::FakeUnionAll) = a.var == b.var && a.body == b.body
Base.:(==)(a::FakeUnion, b::FakeUnion) = a.a == b.a && a.b == b.b
Base.:(==)(a::FakeTypeofBottom, b::FakeTypeofBottom) = true


text:
# 225x225 complex dense
# Compare to S&P, Fig. 60.2

using Pseudospectra, Test

@testset "Complex dense" begin
    opts = Dict{Symbol,Any}(:ax => [-1.1,1.2,-1.1,1.1], :npts => 60,
                            :levels => -5:0.5:-1)
    A = Pseudospectra.landau_fox_li(225,16)

    ps_data = new_matrix(A,opts)
    driver!(ps_data,opts,gs)
    @test iscomputed(ps_data)
    modeplot(ps_data,0,0.0+1.0im)
    modeplot(ps_data,1,0.0+1.0im)

end


text:
import Base: inv, ∘

import Random: AbstractRNG
import Distributions: logpdf, rand, rand!, _rand!, _logpdf, params

#######################################
# AD stuff "extracted" from Turing.jl #
#######################################

abstract type ADBackend end
struct ForwardDiffAD <: ADBackend end
struct TrackerAD <: ADBackend end

const ADBACKEND = Ref(:forward_diff)
function setadbackend(backend_sym)
    @assert backend_sym == :forward_diff || backend_sym == :reverse_diff
    backend_sym == :forward_diff
    ADBACKEND[] = backend_sym
end

ADBackend() = ADBackend(ADBACKEND[])
ADBackend(T::Symbol) = ADBackend(Val(T))
function ADBackend(::Val{T}) where {T}
    if T === :forward_diff
        return ForwardDiffAD
    else
        return TrackerAD
    end
end

######################
# Bijector interface #
######################
"Abstract type for a bijector."
abstract type AbstractBijector end

"Abstract type of bijectors with fixed dimensionality."
abstract type Bijector{N} <:AbstractBijector end

dimension(b::Bijector{N}) where {N} = N
dimension(b::Type{<:Bijector{N}}) where {N} = N

Broadcast.broadcastable(b::Bijector) = Ref(b)

"""
    isclosedform(b::Bijector)::bool
    isclosedform(b⁻¹::Inversed{<:Bijector})::bool

Returns `true` or `false` depending on whether or not evaluation of `b`
has a closed-form implementation.

Most bijectors have closed-form evaluations, but there are cases where
this is not the case. For example the *inverse* evaluation of `PlanarLayer`
requires an iterative procedure to evaluate and thus is not differentiable.
"""
isclosedform(b::Bijector) = true

"""
    inv(b::Bijector)
    Inversed(b::Bijector)

A `Bijector` representing the inverse transform of `b`.
"""
struct Inversed{B <: Bijector, N} <: Bijector{N}
    orig::B

    Inversed(b::B) where {N, B<:Bijector{N}} = new{B, N}(b)
end

inv(b::Bijector) = Inversed(b)
inv(ib::Inversed{<:Bijector}) = ib.orig

"""
    logabsdetjac(b::Bijector, x)
    logabsdetjac(ib::Inversed{<:Bijector}, y)

Computes the log(abs(det(J(b(x))))) where J is the jacobian of the transform.
Similarily for the inverse-transform.

Default implementation for `Inversed{<:Bijector}` is implemented as
`- logabsdetjac` of original `Bijector`.
"""
logabsdetjac(ib::Inversed{<:Bijector}, y) = - logabsdetjac(ib.orig, ib(y))

"""
    forward(b::Bijector, x)
    forward(ib::Inversed{<:Bijector}, y)

Computes both `transform` and `logabsdetjac` in one forward pass, and
returns a named tuple `(rv=b(x), logabsdetjac=logabsdetjac(b, x))`.

This defaults to the call above, but often one can re-use computation
in the computation of the forward pass and the computation of the
`logabsdetjac`. `forward` allows the user to take advantange of such
efficiencies, if they exist.
"""
forward(b::Bijector, x) = (rv=b(x), logabsdetjac=logabsdetjac(b, x))
forward(ib::Inversed{<:Bijector}, y) = (
    rv=ib(y),
    logabsdetjac=logabsdetjac(ib, y)
)

"""
    logabsdetjacinv(b::Bijector, y)

Just an alias for `logabsdetjac(inv(b), y)`.
"""
logabsdetjacinv(b::Bijector, y) = logabsdetjac(inv(b), y)

##############################
# Example bijector: Identity #
##############################

struct Identity{N} <: Bijector{N} end
(::Identity)(x) = x
inv(b::Identity) = b

logabsdetjac(::Identity, x::Real) = zero(eltype(x))
@generated function logabsdetjac(
    b::Identity{N1},
    x::AbstractArray{T2, N2}
) where {N1, T2, N2}
    if N1 == N2
        return :(zero(eltype(x)))
    elseif N1 + 1 == N2
        return :(zeros(eltype(x), size(x, $N2)))
    else
        return :(throw(MethodError(logabsdetjac, (b, x))))
    end
end

########################
# Convenient constants #
########################
const ZeroOrOneDimBijector = Union{Bijector{0}, Bijector{1}}

######################
# Bijectors includes #
######################
# General
include("bijectors/adbijector.jl")
include("bijectors/composed.jl")
include("bijectors/stacked.jl")

# Specific
include("bijectors/exp_log.jl")
include("bijectors/logit.jl")
include("bijectors/scale.jl")
include("bijectors/shift.jl")
include("bijectors/simplex.jl")
include("bijectors/truncated.jl")
include("bijectors/distribution_bijector.jl")

# Normalizing flow related
include("bijectors/planar_layer.jl")
include("bijectors/radial_layer.jl")

##################
# Other includes #
##################
include("transformed_distribution.jl")


text:
#!/usr/bin/env julia
# Useful Functions and Methods
# https://teaching.healthtech.dtu.dk/22110/index.php/Useful_Functions_and_Methods

# Exercise 1
read_data() = try
    if length(ARGS) == 1
        filename = ARGS[1]
    else
        filename = "data/appendix1.txt"
    end
    open(_read_data, filename)
catch e
    msg = sprint(showerror, e)
    println(msg)
    exit(1)
end

function _read_data(f)
    SEQ_flag = false
    seq = ""
    muts = Array{Array{String, 1}, 1}()

    for line in eachline(f)
        FT_reg = match(r"^FT\s+(VARIANT|MUTAGEN)\s+([0-9]+)\s+\S+\s+(\S)\s*->\s*(\S)", line)

        if FT_reg !== nothing
            push!(muts, [FT_reg[x] for x in 1:4])
        end
        if startswith(line, "//")
            SEQ_flag = false
        end
        if SEQ_flag
            seq *= replace(strip(line), " " => "")
        end
        if startswith(line, "SQ")
            SEQ_flag = true
        end
    end

    println(muts)
    println(seq)

    out = open("mutations.fsa", "w")

    for mut in muts
        res_string = mut[2]
        res = parse(Int64, res_string)
        if mut[3] == string(seq[res])
            seq_mut = seq[begin:(res-1)] * mut[4] * seq[(res+1):end]
            println(out, ">", mut[1], mut[2], " ", mut[3], "->", mut[4])
            for i in 1:60:length(seq_mut)
                if i + 59 > length(seq_mut)
                    println(out, seq_mut[i:end])
                else
                    println(out, seq_mut[i:(i + 59)])
                end
            end
        else
            println(mut[3], " and ", seq[res], " did not match.")
        end
    end

    close(out)
end

read_data()

# Exercise 2
using DelimitedFiles

function readpoints(f)
    datapoints = Matrix{Float64}(undef, 6, 0)
    ids = Vector{String}()

    for ln in eachline(f)
        ln_split = split(ln, "\t")
        push!(ids, ln_split[1])
        datapoints = hcat(datapoints, [parse(Float64, x) for x in ln_split[2:end]])
    end

    ids, datapoints
end

function average(a)
    sum(a)/length(a)
end

function highest_lowest()
    translate_arr = readdlm("data/appendix5.txt", String)

    translate_dict = Dict(translate_arr[i, 1] => translate_arr[i, 3]
        for i in 1:size(translate_arr,1))

    exclude_arr = readlines("data/appendix4.txt")

    ids, data = readpoints("data/appendix3.txt")
    translated_exclude = Set([get(translate_dict, id, "") for id in exclude_arr])

    println("Excluded (translated): ", translated_exclude)

    avgs = Vector{Float64}(undef, length(ids))
    avgs_included = Int[]

    for i in eachindex(avgs)
        if ids[i] ∉ translated_exclude
            avgs[i] = average(data[:,i])
            push!(avgs_included, i)
        end
    end

    filtered_avgs = avgs[avgs_included]
    filtered_ids = ids[avgs_included]

    lowest_i = argmin(filtered_avgs)
    highest_i = argmax(filtered_avgs)

    println("Highest average: ", filtered_ids[highest_i], " : ", filtered_avgs[highest_i])
    println("Lowest average: ", filtered_ids[lowest_i], " : ", filtered_avgs[lowest_i])
end

highest_lowest()


text:
using ODEOpt
using ODE
using Base.Test

function testgrav()
    # create time samples
    tmin = 0.0
    tmax = 10
    nsamples = 500
    trange = linspace(tmin, tmax, nsamples)
    #initial position
    ystart = [0.2; 1.0]::Vector{Float64}
    time, y = ode45(grav, ystart, trange)
    x,v = extract(y)
    tol = 1e-2
    @test_approx_eq_eps  time[0.0 .<= abs(x).<=tol][end] 2.9458917835671343 1e-8
end
testgrav()

function noise(σ, n::Integer)
    r = randn(n)
    return σ*r
end

using Optim
using PyPlot
let
# create time samples
srand(0)
tmin = 0.0
tmax = 10
ntimes = 400
trange = linspace(tmin, tmax, ntimes)
#initial position
ystart = [0.0; 1.0]::Vector{Float64}
β = 1.0
βmin = 0.0
βmax = 4.0

σ = 0.01
PLOT = true
samplefactor = 2

function simulate(fgen::Function, β::Float64, sol::Solution, samples, PLOT::Bool)
    t, x,v = extract(sol)
    r = noise(σ, length(samples))
    data = x[samples] .+  r
    tdata = t[samples]
    #= @show tdata[end] t[end] length(t) length(tdata) =#
    # TODO: why are the samples missing the end of the domain ?
    #       the samples IntRange is constructed from a different solution from sol.
    # TODO: sample uniformly rather than at Runge Kutta points.
    if PLOT
        solutionplot(t, x,v)
        plot(tdata, data, "o", label="data")
        legend()
    end

    @show norm(σ*r)
    z = soserr(t, x, tdata, data,)
    #= @test_approx_eq_eps norm(σ*r) z 1e-8 =#

    @show sqrt(soserr(t, x, tdata, data))
    tsolve, xsolve, vsolve, d² = obj(f, ystart, tdata, data)
    @show sqrt(d²)

    γ(x) = obj(fgen(x), ystart, tdata, data)[end]
    @show res = optimize(γ, βmin, βmax)

    β⁺ = res.minimum
    t⁺, y⁺ = ode45(fgen(β⁺), ystart, trange)
    x⁺, v⁺ = extract(y⁺)
    if PLOT
        solutionplot(t⁺, x⁺, v⁺)
    end

    rec =  Record(β, ystart, length(samples),  res)
    return rec
end

recs = Record[]
f = fgen(β)
sol = Solution(ode45(f, ystart, trange)...)
function test_smallerr(recs)
    @test_approx_eq_eps recs[end].β recs[end].res.minimum 1e-1
end

for samplefactor in [1,2,3,4,5,6,7,8,8,9,10]
    samples = 1:samplefactor:length(trange)
    #= @show samples length(trange) length(samples) length(trange[samples]) =#
    r = simulate(fgen, β, sol, samples, false)
    push!(recs, r)
    test_smallerr(recs)
end

println(STDOUT, "n\tβ\tβ⁺\terr\tfval")
for (i,r) in enumerate(recs)
    tprint(STDOUT, r)
end

samplefactor=1
for β in float([1,2,3,4,5,6,7,8,8,9,10])
    βmax = 3β
    f = fgen(β)
    sol = Solution(ode45(f, ystart, trange)...)
    samples = 1:samplefactor:length(trange)
    r = simulate(fgen, β, sol, samples, β==10?true:false)
    push!(recs, r)
    test_smallerr(recs)
end
println(STDOUT, "n\tβ\tβ⁺\terr\tfval")
for (i,r) in enumerate(recs)
    tprint(STDOUT, r)
end
end
include("secondorder.jl")

max_stars_repo_name: jpfairbanks/ODEOpt

max_issues_repo_name: jpfairbanks/ODEOpt

max_forks_repo_name: jpfairbanks/ODEOpt