text:
macro import_scene_huge()
    :(
        using ComplexVisual:
            CV_SceneSetupChain, CV_2DScene,
            CV_MinimalSetupChain, CV_2DMinimalScene,
            cv_get_actionpixel_update, cv_get_statepixel_update,
            cv_get_redraw_func,
            CV_Response,
            cv_setup_2dminimal_scene,
            cv_setup_cycle_state,
            CV_CreateData
    )
end

import Base: show

"""
Scene construction with a chain idea.

## What is a Scene?

A scene is layout together with all the callbacks/functions to update/draw all
relevant parts. Often painters depends on "degrees of freedom" that change.
[The change of such degrees of freedom may be triggered by a mouse or other
"events".]

Scenes are constructed step-by-step. If an new painter is "added" to the
scene, there must be a possibility to ensure, that also this new painter
is called if the degrees of freedom changed.

Here is, where a `CV_SceneSetupChain` helps.

## How does this work?

All callback functions are gathered in vectors. This is the point where
the type information about the callback functions are completely lost.
For the return value (which is always `nothing`) this is not a problem.
For the call this is a tradeoff: If elements of such a vector are
called then this are purely "runtime"-calls (i.e. what to call, with
what types, etc. is determined at runtime). That's the negative part.
The positive part: Julia doesn't need to do the type-book-keeping (which
is rather tough if e.g. types informations should not be lost and all
these methods are nested: meaning every new method calls the old one first;
because every call is typically are closure with a lot of types;
preserving the type-info with this way would result in an unbearable
compile-time).

If a new part (e.g. a painter) changes the layout then a new
`CV_SceneSetupChain` is built (with the new layout). This is very important,
because for the (nested) layouts all type informations are preserved (also
in the `CV_SceneSetupChain`). So the painter calls (which use parts of
the layout) can make use of this type informations.

All the given callback functions are appended to the callback vectors.

Required fields:
* `layout`
* `draw_once_func`:      callbacks that are called after the layout is fixed
                         and a Layout-Canvas was constructed.
                         Argument: last layout
* `actionpixel_update`:  callbacks that are called after "the main action"
                         occured (e.g. mouse dragged).
                         Arguments: pixel_x, pixel_y (both w.r.t. Layout Canvas)
* `statepixel_update`:   callbacks that are called after the
                         "state-change action" occured.
                         Arguments: pixel_x, pixel_y (both w.r.t. Layout Canvas)
* `redraw_func`:         callbacks that are called after a redraw
                         is triggered (e.g. because some trafo-parameters have
                         changed; typically all painters react and clear their
                         caches and recompute their paintings)
"""
abstract type CV_SceneSetupChain end

"""
Possible response for some (scene-)callback functions
"""
struct CV_Response
    redraw_flag  :: Bool
end

CV_Response(; redraw_flag=false) = CV_Response(redraw_flag)

"""
A minimalistic `CV_SceneSetupChain`.
"""
struct CV_MinimalSetupChain{layoutT} <: CV_SceneSetupChain # {{{
    layout               :: layoutT
    draw_once_func       :: Vector{Any}
    actionpixel_update   :: Vector{Any}
    statepixel_update    :: Vector{Any}
    redraw_func          :: Vector{Any}
end

function CV_MinimalSetupChain(layout)
    return CV_MinimalSetupChain(layout, Vector(), Vector(), Vector(), Vector())
end

"""
"update" a scene-setup by new layout and/or by adding callback-functions.
"""
function cv_combine(old::CV_MinimalSetupChain;
        layout=missing, draw_once_func=missing,
        actionpixel_update=missing, statepixel_update=missing,
        redraw_func=missing)
    new = ismissing(layout) ? old : CV_MinimalSetupChain(layout,
        old.draw_once_func, old.actionpixel_update, old.statepixel_update,
        old.redraw_func)
    !ismissing(draw_once_func) && push!(new.draw_once_func, draw_once_func)
    !ismissing(actionpixel_update) && push!(new.actionpixel_update, actionpixel_update)
    !ismissing(statepixel_update) && push!(new.statepixel_update, statepixel_update)
    !ismissing(redraw_func) && push!(new.redraw_func, redraw_func)
    return new
end

# }}}

"""
A scene is a framed layout which has everything to be shown and to be used
interactively.

A `CV_2DScene` must support:
`cv_get_actionpixel_update`, `cv_get_statepixel_update`,
`cv`get_redraw_func`
"""
abstract type CV_2DScene <: CV_Framed2DLayout
    # actionpixel_update
    # statepixel_update
    # redraw_func
end

"""
A minimalistic `CV_2DScene`.
"""
struct CV_2DMinimalScene{parentT, apuT, spuT, rfT} <: CV_2DScene    # {{{
    parent_layout            :: parentT
    actionpixel_update       :: apuT
    statepixel_update        :: spuT
    redraw_func              :: rfT
end

@layout_composition_getter(actionpixel_update,      CV_2DMinimalScene)
@layout_composition_getter(statepixel_update,       CV_2DMinimalScene)
@layout_composition_getter(redraw_func,             CV_2DMinimalScene)

function cv_setup_2dminimal_scene(setup::CV_SceneSetupChain)
    can_layout = cv_canvas_for_layout(setup.layout)
    cc_can_layout = cv_create_context(can_layout)
    layout = CV_MinimalFramed2DLayout(setup.layout, can_layout, cc_can_layout)

    scene_redraw_func = () -> begin
        for func in setup.redraw_func
            func(layout)
        end
        return nothing
    end

    scene_actionpixel_update = (px, py) -> begin
        redraw_flag = false
        for func in setup.actionpixel_update
            resp = func(px, py, layout)
            if resp isa CV_Response  && resp.redraw_flag
                redraw_flag = true
            end
        end
        if redraw_flag
            scene_redraw_func()
        end
        return nothing
    end

    scene_statepixel_update = (px, py) -> begin
        redraw_flag = false
        for func in setup.statepixel_update
            resp = func(px, py, layout)
            if resp isa CV_Response && resp.redraw_flag
                redraw_flag = true
            end
        end
        if redraw_flag
            scene_redraw_func()
        end
        return nothing
    end

    scene = CV_2DMinimalScene(layout,
        scene_actionpixel_update, scene_statepixel_update,
        scene_redraw_func)
    for func in setup.draw_once_func
        func(scene)
    end

    return cv_combine(setup; layout=scene)
end
# }}}

"""
If some advanced/complicated graphic objects are constructed (e.g. a slider
with rulers, etc.) then often such `cv_create_...` methods may want to
return further informations about the constructed objects. This is done
in concrete subtype of this (abstract) type.
"""
abstract type CV_CreateData end

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
# The original Java version is Copyright 2014 Google Inc.
# This Julia version is Copyright 2022 Ted Dunning
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

md"""
An `S1ChordAngle` expresses an angle on a sphere but expresses this angle in terms
of the square of the length of the chord between the two points that define the angle.

This is equivalent of the `haversine` which was used before computers to simplify
spherical trigonometry computations and allows many common operations to be done
without trigonometry. For instance, the sine of a `S1ChordAngle` can be computed
using a single square root operation. An `S1ChordAngle` is, however, limited to
angles in the range from `[0, π]` which corresponds to the range `[0, 4]`.

For large angles, there can be some loss of accuracy. Near the maximum value, the
representation of an `S1ChordAngle` can be expected to have an error of about
`(10^{-15} / x)`.

This file is an idiomatic translation of the Java implementation of S2
"""
struct S1ChordAngle
    distance2::Float64

    S1ChordAngle(distance2::Number) = S1ChordAngle(distance2, false)

    function S1ChordAngle(distance2::Number, force::Bool)
        if force
            new(distance2)
        else
            new(min(S1ChordAngle_MAX_LENGTH2, distance2))
        end
    end

    function S1ChordAngle(x::Vector, y::Vector)
        # assert (norm2(x) == 1 && norm2(y) == 1)
        # The distance may slightly exceed 4.0 due to roundoff errors
        S1ChordAngle(distance(x, y)^2)
    end

    # Returns a new chord angle approximated from {@code angle} (see
    # getS1AngleConstructorMaxError() for the max magnitude of the error).
    #
    # Angles outside the range [0, Pi] are handled as follows:
    #
    #

    #
{@link S1Angle#INFINITY} is mapped to {@link #INFINITY}
    #
negative angles are mapped to {@link #NEGATIVE}
    #
finite angles larger than Pi are mapped to {@link #STRAIGHT}
    #

    #
    # Note that this operation is relatively expensive and should be avoided. To use {@link
    # S1ChordAngle} effectively, you should structure your code so that input arguments are converted
    # to S1ChordAngles at the beginning of your algorithm, and results are converted back to `S1Angle`s
    # only at the end.

    function S1ChordAngle(angle::S1Angle)
        if (radians(angle) < 0)
            return S1ChordAngle_NEGATIVE;
        elseif (angle == S1Angle_INFINITY)
            return S1ChordAngle_INFINITY
        else
            # The chord length is 2 * sin(angle / 2)
            length = 2 * sin(0.5 * min(π, radians(angle)))
            return S1ChordAngle(length * length);
        end
    end
end

const S1ChordAngle_MAX_LENGTH2 = 4.0
const S1ChordAngle_ZERO = S1ChordAngle(0)
const S1ChordAngle_RIGHT = S1ChordAngle(2)
const S1ChordAngle_STRAIGHT = S1ChordAngle(4)
const S1ChordAngle_INFINITY = S1ChordAngle(Inf, true)
const S1ChordAngle_NEGATIVE = S1ChordAngle(-1)

isNegative(a::S1ChordAngle) = a.distance2 < 0
isInfinity(a::S1ChordAngle) = isinf(a.distance2)
isSpecial(a::S1ChordAngle) = isNegative(a) || isInfinity(a)
isZero(a::S1ChordAngle) = a.distance2 == 0

# Returns true if getLength2() is within the normal range of 0 to 4 (inclusive) or the angle is
# special.

function isValid(a::S1ChordAngle)::Bool
    return (a.distance2 >= 0 && a.distance2 <= MAX_LENGTH2) || isNegative(a) || isInfinity(a)
end

# Convert the chord angle to an {@link S1Angle}. {@link #INFINITY} is converted to {@link
# S1Angle#INFINITY}, and {@link #NEGATIVE} is converted to a negative {@link S1Angle}. This
# operation is relatively expensive.

function toAngle(a::S1ChordAngle)
    if (isNegative(a))
        return radians(-1);
    elseif (isInfinity(a))
        return S1Angle_INFINITY;
    else
        return radians(2 * asin(0.5 * sqrt(a.distance2)))
    end
end

# Returns a new S1ChordAngle whose chord distance represents the sum of the angular distances
# represented by the 'a' and 'b' chord angles.
#
# Note that this method is much more efficient than converting the chord angles to S1Angles
# and adding those. It requires only one square root plus a few additions and multiplications.

function (+)(a::S1ChordAngle, b::S1ChordAngle)
    if isSpecial(a) || isSpecial(b)
        throw(ArgumentError("Can't add special chord angles"))
    end

    let a2 = a.distance2
        b2 = b.distance2
        begin
            # Optimization for the common case where "b" is an error tolerance parameter that happens to be
            # set to zero.
            if (b2 == 0)
                return a
            end

            # Clamp the angle sum to at most 180 degrees.
            if (a2 + b2 >= S1ChordAngle_MAX_LENGTH2)
                return S1ChordAngle_STRAIGHT
            end

            # Let "a" and "b" be the (non-squared) chord lengths, and let c = a+b.
            # Let A, B, and C be the corresponding half-angles (a = 2*sin(A), etc).
            # Then the formula below can be derived from c = 2 * sin(A+B) and the relationships
            #   sin(A+B) = sin(A)*cos(B) + sin(B)*cos(A)
            #   cos(X) = sqrt(1 - sin^2(X)) .
            x = a2 * (1 - 0.25 * b2); # isValid() => non-negative
            y = b2 * (1 - 0.25 * a2); # isValid() => non-negative
            return S1ChordAngle(min(S1ChordAngle_MAX_LENGTH2, x + y + 2 * sqrt(x * y)))
        end
    end
end

md"""
Subtract one S1ChordAngle from another.

Note that this method is much more efficient than converting the chord angles to S1Angles
and adding those. It requires only one square root plus a few additions and multiplications.
"""
function (-)(a::S1ChordAngle, b::S1ChordAngle) :: S1ChordAngle
    # See comments in add(S1ChordAngle, S1ChordAngle).
    if isSpecial(a) || isSpecial(b)
        throw(ArgumentError("Can't subtract special chord angles"))
    end
    let a2 = a.distance2
        b2 = b.distance2
        begin
            if (b2 == 0)
                return a
            end
            if (a2 <= b2)
                return S1ChordAngle_ZERO;
            end
            let x = a2 * (1 - 0.25 * b2)
                y = b2 * (1 - 0.25 * a2)
                return S1ChordAngle(max(0.0, x + y - 2 * sqrt(x * y)))
            end
        end
    end
end

md"""
Returns the smaller of the given instances
"""
Base.min(a::S1ChordAngle, b::S1ChordAngle) = a.distance2 <= b.distance2 ? a : b

md"""
Returns the larger of the given instances
"""
Base.max(a::S1ChordAngle, b::S1ChordAngle) = a.distance2 > b.distance2 ? a : b

md"""
Returns the square of Math.sin(toAngle().radians()), but computed more efficiently
"""
function sin2(a::S1ChordAngle)
    if isSpecial(a)
        throw(ArgumentError("Can't take the sin of a special angle"))
    end
    # Let "a" be the (non-squared) chord length, and let A be the corresponding half-angle
    # (a = 2*sin(A)). The formula below can be derived from:
    #   sin(2*A) = 2 * sin(A) * cos(A)
    #   cos^2(A) = 1 - sin^2(A)
    # This is much faster than converting to an angle and computing its sine.
    return a.distance2 * (1 - 0.25 * a.distance2)
end

md"""
Returns Math.sin(toAngle().radians()), but computed more efficiently
"""
Base.sin(a::S1ChordAngle) = sqrt(sin2(a))

md"""
Returns Math.cos(toAngle().radians()), but computed more efficiently
"""
function Base.cos(a::S1ChordAngle)
    # cos(2*A) = cos^2(A) - sin^2(A) = 1 - 2*sin^2(A)
    if isSpecial(a)
        throw(ArgumentError("Can't take the cosine of a special angle"))
    end

    return 1 - 0.5 * a.distance2
end

md"""
Returns Math.tan(toAngle().radians()), but computed more efficiently.
"""
Base.tan(a::S1ChordAngle) = sin(a) / cos(a);

md"""
Returns a new S1ChordAngle that has been adjusted by the given error bound (which can be
positive or negative). {@code error} should be the value returned by one of the error bound
methods below. For example:

```
   {@code S1ChordAngle a = new S1ChordAngle(x, y);}
   {@code S1ChordAngle a1 = a.plusError(a.getS2PointConstructorMaxError());}
```

If this {@link #isSpecial}, we return {@code this}.
"""
function plusError(error::Float64)
    return isSpecial() ? this : fromLength2(max(0.0, min(S1ChordAngle_MAX_LENGTH2, length2 + error)))
end

md"""
Returns the error in {@link #fromS1Angle}
"""
getS1AngleConstructorMaxError() = S2.DBL_EPSILON * length2

md"""
There is a relative error of {@code 2.5 * DBL_EPSILON} when computing the squared distance,
plus a relative error of {@code 2 * DBL_EPSILON} and an absolute error of {@code 16 *
DBL_EPSILON^2} because the lengths of the input points may differ from 1 by up to {@code 2 *
DBL_EPSILON} each. (This is the maximum length error in {@link S2Point#normalize}).
"""
getS2PointConstructorMaxError() = (4.5 * S2.DBL_EPSILON * length2) + (16 * S2.DBL_EPSILON * S2.DBL_EPSILON)

ext: jl


text:
using Distributed
addprocs(1);
using Revise
using LinearAlgebra
@everywhere  using Pkg; Pkg.activate(".")
@everywhere  push!(LOAD_PATH, "/Users/mtg79/Documents/FastTransforms.jl/src/")
@everywhere  using Revise, FastTransforms, Distributed, DistributedArrays

##
x = rand(10^2)
c = complex(x)

## serial
@time A = plan_nufft1(x)
@time b = A*c

## parallel
@time plan =parallel_plan_nufft1(x)
@time f = mulpar3( plan, c)


text:
using LightGraphsMatching

struct Matching
    g::DiGraph
    w::Dict{Edge,Float64}
    solver_instance::Any
    function Matching(graph::DiGraph, weights::Array{Float64,2}, ndds::Array{Int64, 1}, solver_instance::Any)
        n = nv(graph)

        new(g, w, solver_instance)
    end
end

function solve(problem::Matching)
    match = maximum_weight_matching(g,with_optimizer(Cbc.Optimizer, logLevel=0),w)
    return problem.m, problem.x, problem.in_, nothing, nothing
end


text:
#! /usr/bin/julia

# Rosetta Code, Ordered words

function isordered{T<:String}(w::T)
    p = '\0'
    for c in w
        p <= c || return false
        p = c
    end
    return true
end

maxlen = 0
wlst = String[]
WL = open("ordered_words.txt", "r")

for w in eachline(WL)
    w = chomp(w)
    wlen = length(w)
    wlen>=maxlen && isordered(w) || continue
    if wlen > maxlen
        maxlen = wlen
        wlst = [w]
    else
        push!(wlst, w)
    end
end
close(WL)

for w in wlst
    println("    ", w)
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
ENV["GKSwstype"] = "100"  # set 'GR environment' to 'no output' (for Travis CI)
using Documenter, Reachability

DocMeta.setdocmeta!(Reachability, :DocTestSetup, :(using Reachability); recursive=true)

makedocs(
    sitename = "Reachability.jl",
    modules = [Reachability],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets = ["assets/juliareach.css"]),
    pages = [
        "Home" => "index.md",
        "Library" => Any[
        "User interface" => "lib/interface.md",
        "Systems" => "lib/systems.md",
        "Algorithms" => "lib/algorithms.md",
        "Transformations" => "lib/transformations.md",
        "Discretization" => "lib/discretize.md",
        "Distributed computations" => "lib/distributed.md"],
        "Publications" => "publications.md",
        "Citations" => "citations.md",
        "About" => "about.md"
    ]
)

deploydocs(
    repo = "github.com/JuliaReach/Reachability.jl.git"
)


text:
include("read_matpower.jl")
include("read_psse.jl")
include("read_go.jl")

const get_data_functions = Dict(
    "MATPOWERM" => get_data_matpower_m,
    "MATPOWERMAT" => get_data_matpower_mat,
    "RAW" => get_data_raw,
    "RAWX" => get_data_rawx,
    "RAWGO" => get_data_rawgo
)
const nbus_functions = Dict(
    "MATPOWERM" => nbus_matpower_m,
    "MATPOWERMAT" => nbus_matpower_mat,
    "RAW" => nbus_raw,
    "RAWX" => nbus_rawx,
    "RAWGO" => nbus_rawgo
)
const nbranch_functions = Dict(
    "MATPOWERM" => nbranch_matpower_m,
    "MATPOWERMAT" => nbranch_matpower_mat,
    "RAW" => nbranch_raw,
    "RAWX" => nbranch_rawx,
    "RAWGO" => nbranch_rawgo
)
const ngen_functions = Dict(
    "MATPOWERM" => ngen_matpower_m,
    "MATPOWERMAT" => ngen_matpower_mat,
    "RAW" => ngen_raw,
    "RAWX" => ngen_rawx,
    "RAWGO" => ngen_rawgo
)
const ntransformer_functions = Dict(
    "MATPOWERM" => ntransformer_matpower_m,
    "MATPOWERMAT" => ntransformer_matpower_mat,
    "RAW" => ntransformer_raw,
    "RAWX" => ntransformer_rawx,
    "RAWGO" => ntransformer_rawgo
)

function get_network_name(path::AbstractString, format::AbstractString)
    if format == "MATPOWERM"
        return basename(path)[begin:end - 2]
    elseif format == "RAWGO"
        path = _resolve_rawgo_path(path, "dir")
        paths = splitpath(path)
        name = paths[end - 1]
        scenario = parse(Int, paths[end][begin + 9: end])
        return "$(name)_$(scenario)"
    else
        error("Not Implemented")
    end
end

function convert_cols_to_int!(df::DataFrame, cols)
    for c in cols
        df[!, c] = floor.(Int, df[!, c])
    end
end

function read_network(path::AbstractString, format::AbstractString)
    bus, gen, branch, gencost, baseMVA = get_data_functions[format](path)
    colnames = get_matpower_cols()

    bus = DataFrame(bus, colnames["bus"])
    convert_cols_to_int!(bus, [:ID, :TYPE])

    gen = DataFrame(gen, colnames["gen"])
    convert_cols_to_int!(gen, [:ID, :GEN_STATUS])

    branch = DataFrame(branch, colnames["branch"])
    convert_cols_to_int!(branch, [:SRC, :DST, :BR_STATUS])

    gencost = DataFrame(gencost, get_matpower_gencost_cols(gencost))
    convert_cols_to_int!(gencost, [:MODEL, :NCOST])

    name = get_network_name(path, format)
    return PowerFlowNetwork(name, bus, gen, branch, gencost, baseMVA)
end

nbus(path::AbstractString, format::AbstractString) = nbus_functions[format](path)
nbranch(path::AbstractString, format::AbstractString; distinct_pair=false) = nbranch_functions[format](path; distinct_pair) +
                                                                             ntransformer_functions[format](path; distinct_pair)
ngen(path::AbstractString, format::AbstractString; distinct_pair=false) = ngen_functions[format](path; distinct_pair)

function _get_source_type(path::AbstractString)
    if splitext(path)[end] == ".m"
        return "MATPOWERM"
    else
        # TODO: better distinction
        return "RAWGO"
    end
end

function nbus(path::AbstractString)
    source_type = _get_source_type(path)
    return nbus(path, source_type)
end
function nbranch(path::AbstractString; distinct_pair=false)
    source_type = _get_source_type(path)
    return nbranch(path, source_type; distinct_pair=distinct_pair)
end
function ngen(path::AbstractString; distinct_pair=false)
    source_type = _get_source_type(path)
    return ngen(path, source_type; distinct_pair=distinct_pair)
end

nbranch_unique(path::AbstractString) = nbranch(path::AbstractString; distinct_pair=true)
ngen_unique(path::AbstractString) = ngen(path::AbstractString; distinct_pair=true)

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
import SymbolicUtils: quick_cancel; export quick_cancel
import SymbolicUtils: Postwalk

# change SymbolicUtils' quick_cancel to simplify powers of fractions correctly
function quick_cancel(x::Term, y::Term)
	if x.f == exp && y.f == exp
		return exp(x.arguments[1] - y.arguments[1]), 1
	else
		return x,y
	end
end

quick_cancel(x::Term, y::Pow) = y.base isa Term && y.base.f == exp ? quick_cancel(x, expand_exp_power(y)) : x,y

"Returns true if expr is an exponential"
is_exp(expr) = expr isa Term && expr.f == exp

"Expand powers of exponential such that exp(x)^n => exp(x*n) "
expand_exp_power(expr) = expr isa Pow && is_exp(expr.base) ? exp(expr.base.arguments[1] * expr.exp) : expr
expand_exp_power(expr::Add) = sum([expand_exp_power(arg) for arg in arguments(expr)])
expand_exp_power(expr::Mul) = prod([expand_exp_power(arg) for arg in arguments(expr)])
expand_exp_power(expr::Num) = expand_exp_power(expr.val)

"Expands using SymbolicUtils.expand and expand_exp_power (changes exp(x)^n to exp(x*n)"
expand_all(x) = Postwalk(expand_exp_power)(SymbolicUtils.expand(x))
expand_all(x::Complex{Num}) = expand_all(x.re) + im* expand_all(x.im)
expand_all(x::Num) = Num(expand_all(x.val))

"Apply a function f on every member of a sum or a product"
_apply_termwise(f, x::Add) = sum([f(arg) for arg in arguments(x)])
_apply_termwise(f, x::Mul) = prod([f(arg) for arg in arguments(x)])
_apply_termwise(f, x::Div) = _apply_termwise(f, x.num) / _apply_termwise(f, x.den)
_apply_termwise(f,x) = f(x)

function simplify_complex(x::Complex)
    isequal(x.im, 0) ? x.re : x.re + im*x.im
end

simplify_complex(x::Add) = _apply_termwise(simplify_complex, x)
simplify_complex(x::Union{Mul, Div}) = _apply_termwise(simplify_complex, x)
simplify_complex(x) = x

"Simplify products of exponentials such that exp(a)*exp(b) => exp(a+b)
This is included in SymbolicUtils as of 17.0 but the method here avoid other simplify calls"
function simplify_exp_products(expr::Mul)
	ind = findall(x -> is_exp(x), arguments(expr))
	rest_ind = setdiff(1:length(arguments(expr)), ind)
	rest = isempty(rest_ind) ? 1 : prod(arguments(expr)[rest_ind])
	total = isempty(ind) ? 0 : sum(getindex.(getfield.(arguments(expr)[ind], Symbol("arguments")), 1))
	SymbolicUtils.is_literal_number(total) ? (total == 0 && return rest) : return rest * exp(total)

end

#simplify_exp_products(expr::Add) = sum([simplify_exp_products(arg) for arg in arguments(expr)])
simplify_exp_products(x::Complex{Num}) = Complex{Num}(simplify_exp_products(x.re.val), simplify_exp_products(x.im.val))
simplify_exp_products(x::Union{Add, Div}) = _apply_termwise(simplify_exp_products, x)
simplify_exp_products(x::Num) = simplify_exp_products(x.val)
simplify_exp_products(x) = x

function exp_to_trig(x)
	if x isa Term && x.f == exp
		arg = first(x.arguments)
		trigarg = -im*arg # the argument of the to-be trig function
		trigarg = simplify_complex(trigarg)

		# put arguments of trigs into a standard form such that sin(x) = -sin(-x), cos(x) = cos(-x) are recognized

		if trigarg isa Add

			first_symbol = minimum(cat(string.(arguments(trigarg)), string.(arguments(-trigarg)), dims=1))

			# put trigarg => -trigarg the lowest alphabetic argument of trigarg is lower than that of -trigarg
			# this is a meaningless key but gives unique signs to all sums
			is_first = minimum(string.(arguments(trigarg))) == first_symbol
			return is_first ? cos(-trigarg) -im*sin(-trigarg) : cos(trigarg)+im* sin(trigarg)
		end

		trigarg isa Mul && trigarg.coeff < 0 ? cos(-trigarg) -im*sin(-trigarg) : cos(trigarg)+im* sin(trigarg)

	else
		x
	end
end

exp_to_trig(x::Union{Add, Mul, Div}) = _apply_termwise(exp_to_trig, x)
exp_to_trig(x::Num) = exp_to_trig(x.val)
exp_to_trig(x::Complex{Num}) = exp_to_trig(x.re)+im* exp_to_trig(x.im)

# sometimes, expressions get stored as Complex{Num} with no way to decode what real(x) and imag(x)
# this overloads the Num constructor to return a Num if x.re and x.im have similar arguments
Num(x::Complex{Num}) = isequal(x.re.val.arguments, x.im.val.arguments) ? Num(first(x.re.val.arguments)) : error("Cannot convert Complex{Num} " * string(x) * " to Num")

#=
chop(x) = x
chop(x::Complex{Int64})= Int64(x)
chop(x::Complex{Float64}) = Float64(x)
chop(x::Add) = _apply_termwise(chop, x)
chop(x::Mul) = _apply_termwise(chop, x)
chop(x::Num) = chop(x.val)
chop(x::Complex{Num}) = Complex{Num}(x.re, x.im)

#expand_fraction(expr::Div) = expr.num isa Add ? sum([arg / expr.den for arg in arguments(expr.num)]) : expr
#expand_fraction(expr) = expr

#simplify_parts(expr::Num) = simplify_parts(expr.val)
#simplify_parts(expr::Add) = sum([simplify_fractions(arg) for arg in arguments(expr)])
#simplify_parts(expr) = simplify_fractions(expr)
=#


text:
#
println("Ae) Test of the CI part for an internally generated neon multiplet without Breit interaction.")
#
settings = AsfSettings(AsfSettings(),selectLevelsCI = true, selectedLevelsCI = [1,2, 4,5, 7,8], breitCI=false,
                                     selectSymmetriesCI = true, selectedSymmetriesCI = [ LevelSymmetry(1//2,Basics.plus),
                                                          LevelSymmetry(1//2,Basics.minus),  LevelSymmetry(5//2,Basics.plus)], jjLS = LSjjSettings(false) )

wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=JAC.Radial.Grid(true), nuclearModel=Nuclear.Model(26.),
                        configs=[Configuration("[Ne] 3s^2 3p^5"), Configuration("[Ne] 3s 3p^6"), Configuration("[Ne] 3s^2 3p^4 3d")],  #,
                                 # Configuration("[Ne] 3s^2 3p^2 3d^3"), Configuration("[Ne] 3s 3p^2 3d^4"), Configuration("[Ne] 3p^2 3d^5")],
                        asfSettings=settings  )

@time wb = perform(wa)


text:
# script to run all SLICOT models in several benchmark settings

# Plots is required here, unless never used, to avoid an issue with "world age"
using Reachability, Plots

# model list
models = String[
    "motor",
    "building",
    "pde",
    "heat",
    "iss",
    "beam",
    "mna1",
    "fom",
    "mna5"
   ]

# load models and benchmark scripts
model_library = include(@relpath "model_library.jl")
for model in models
    path = (@relpath "../") * model_library[model]
    include(path)
end
include(@relpath "reach_1D_single.jl")
include(@relpath "reach_1D_all.jl")
include(@relpath "reach_2D_box_two.jl")
include(@relpath "reach_2D_eps_two.jl")
include(@relpath "reach_kD_all.jl")
include(@relpath "check_1D_discrete.jl")
include(@relpath "check_kD_discrete.jl")
include(@relpath "check_1D_dense.jl")
include(@relpath "check_kD_dense.jl")

# create plots (option)?
create_plots = true
if create_plots
    include(@relpath "plot_reach.jl")
end

# run benchmarks

# -- reach_1D_single --
# no modification, optionally create plots
reach_1D_single(models, create_plots)

# -- reach_1D_all --
# remove mna5
models_filtered = filter(e -> e != "mna5", models)
reach_1D_all(models_filtered)

# -- reach_2D_box_two --
# no modification, optionally create plots
reach_2D_box_two(models, create_plots)

# -- reach_2D_eps_two --
# no modification, optionally create plots
reach_2D_eps_two(models, create_plots)

# -- reach_kD_all --
# only use pde, optionally create plots
models_filtered = filter(e -> e ∈ ["fom"], models)
reach_kD_all(models_filtered, create_plots)

# -- check_1D_discrete --
# no modification
check_1D_discrete(models)

# -- check_kD_discrete --
# only use iss and fom; use default partition
models_filtered = filter(e -> e ∈ ["iss", "fom"], models)
check_kD_discrete(models_filtered)

# -- check_1D_dense --
# remove iss and fom; define deltas
models_filtered = filter(e -> e ∉ ["iss", "fom"], models)
delta_map = Dict(
    "motor" => 1e-3,
    "building" => 2e-3,
    "pde" => 3e-4,
    "heat" => 1e-3,
    "beam" => 5e-5,
    "mna1" => 4e-4,
    "mna5" => 3e-1
   )
deltas = Vector{Float64}(length(models_filtered))
for (i, model) in enumerate(models_filtered)
    deltas[i] = delta_map[model]
end
check_1D_dense(models_filtered, deltas)

# -- check_kD_dense --
# only use iss; use default partition
models_filtered = filter(e -> e ∈ ["iss"], models)
delta_map = Dict(
    "iss" => 6e-4,
   )
deltas = Vector{Float64}(length(models_filtered))
for (i, model) in enumerate(models_filtered)
    deltas[i] = delta_map[model]
end
check_kD_dense(models_filtered, deltas)

return nothing

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
#=
MIT License

Optimal Testing and Containment Strategies for Universities in Mexico amid COVID-19

Copyright © 2021 Test and Contain. Luis Benavides-Vázquez, Héctor Alonso Guzmán-Gutiérrez, Jakob Jonnerby, Philip Lazos, Edwin Lock, Francisco J. Marmolejo-Cossío, Ninad Rajgopal,and José Roberto Tello-Ayala. https://www.testandcontain.com/

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
=#
using DataFrames
using CSV

include("frontier.jl")

function ranges(pop, config)
    k = length(pop.p)
    prealloc = preallocate(k)
    budget, step, guaranteed = config.budget, config.stepsize, config.guaranteed
    gsize = config.groupsizes
    mins = fill(Inf, k+1)
    maxes = fill(-Inf, k+1)
    do_nothing = critical(pop, prealloc, ones(Int, k), zeros(Int, k))
    @showprogress for groups in GroupIterator(gsize, k)
        for tests in TestIterator(budget, step, guaranteed, pop.n .÷ groups, k)
            objectives!(pop, prealloc, groups, tests, do_nothing)
            for i in 1:k+1
                mins[i] = min(mins[i], prealloc.obj[i])
                maxes[i] = max(maxes[i], prealloc.obj[i])
            end
        end
    end
    return mins, maxes
end

"""Compute the frontier, loading parameters from `filename` file."""
function compute(filename, nobuckets=false, buckets=nothing)
    pop, config = load_from_json(filename)
    if buckets != nothing
        config.bucketsizes = buckets
    end
    solutions = frontier(pop, config, nobuckets)
    return solutions
end

function run()
    # Process commandline arguments
    if length(ARGS) === 0
        println("Please specify the input data file. Optionally, specify an output file as a second argument.")
        println("\nUsage: run `julia path/to/data.json [path/to/output.csv]` in the terminal.")
        return
    end
    input = ARGS[1]
    n = parse(Int32, ARGS[2])
    output = length(ARGS)==3 ? ARGS[3] : nothing
    pop, config = load_from_json(input)

    k = length(pop.p) # number of categories
    total_sol = 3k+1

    # Creates initial bucket with the variance of each objective
    mins, maxes = ranges(pop,config)
    buckets_var = maxes - mins

    # Sets lower and uper bounds for binary search (interval [0,1])
    ubound = 1.0
    lbound = 0.0

    solutions_no = Solutions()
    matrix_no = hcat(solutions_no...)'
    buckets = buckets_var
    # Iterates ten times by checking number of solution in each run of compute
    for i in 1:10
        α = (ubound + lbound)/2
        buckets = α * buckets_var
        solutions_no = compute(input, false, buckets)
        matrix_no = hcat(solutions_no...)'
        total_sol = size(matrix_no, 1)
        if total_sol > n
            lbound = α
        else
            ubound = α
        end
    end
    print("Buckets")
    print(*)
    print(round.(buckets,digits=2))
    print(*)
    # Final solutions to print
    solutions = solutions_no
    matrix = matrix_no
    #

    header = [["g$i" for i in 1:k];
             ["t$i" for i in 1:k];
             ["h"];
             ["c$i" for i in 1:k]]
    df = DataFrame(matrix, header)  # TODO: make pretty table headings

    if output === nothing
        # Print solutions to terminal
        println(df)
        # Print some stats
        println("Solution stats per objective")
        println(describe(df[!, 2k+1:3k+1]))
    else
        # Save solutions to output filename
        CSV.write(output, df, writeheader=true)
    end

    return
end

run()

max_stars_repo_name: hguzmang/eaamo

max_issues_repo_name: hguzmang/eaamo

max_forks_repo_name: hguzmang/eaamo


text:
## ---------------  VarInfo

"""
	$(SIGNATURES)

Holds meta information about one variable. `options` holds predefined options and user defined options.
Predefined options are:
- size (required)
- lb, ub
- lbs, ubs
- dimSums
- grandSum

These can be used to check a data matrix.

# Example
```julia
vi = VarInfo{Vector{Int}}(:x, "x stores a Vector of Int"; size = (3,), lb = 0);
eltype(vi) == Int;
```
"""
struct VarInfo{T}
    varName :: Symbol
    description :: String
    options :: Dict{Symbol, Any}
end

"""
	$(SIGNATURES)

Constructor that provides options as keyword arguments.
"""
function VarInfo(vName :: Symbol, vDescr :: AbstractString, T :: DataType; kwargs...)
    d = dict_from_kwargs(; kwargs...);
    if T <: Number
        d[:size] = ();
    end
    return VarInfo{T}(vName, vDescr, d);
end

function dict_from_kwargs(; kwargs...)
    d = Dict{Symbol, Any}();
    for kw in kwargs
        d[kw[1]] = kw[2];
    end
    return d
end

## ----------  Extended Base functions

Base.size(v :: VarInfo{T}) where {T} = get_option(v, :size);
Base.eltype(v :: VarInfo{T}) where {T} = eltype(T);

function Base.isequal(v1 :: VarInfo{T}, v2 :: VarInfo{T}) where {T}
    isEqual = isequal_common(v1, v2)  &&  isequal_options(v1, v2);
    return isEqual
end

Base.show(io :: IO, vi :: VarInfo{T}) where {T} =
    print(io, "VarInfo for $(vi.varName) of size $(size(vi))");

## -----------  Access

"""
	$(SIGNATURES)

Returns `T` for `VarInfo{T}`.
"""
value_type(vi :: VarInfo{T}) where T = T;

"""
	$(SIGNATURES)

Retrieve an option. Return `nothing` if option not found.
"""
function get_option(vi :: VarInfo{T}, oName :: Symbol) where T
    if has_option(vi, oName)
        o = vi.options[oName];
    else
        o = nothing;
    end
    return o
end

"""
	$(SIGNATURES)

Get multiple options. `nothing` if not found.
"""
function get_options(vi :: VarInfo{T}, oNames) where T
    return [get_option(vi, oName)  for oName in oNames];
end

"""
	$(SIGNATURES)

Does the `VarInfo` have an option?
"""
has_option(vi :: VarInfo{T}, oName :: Symbol) where T =
    haskey(vi.options, oName);

"""
	$(SIGNATURES)

Returns name of the variable.
"""
var_name(v :: VarInfo{T}) where T = v.varName;

"""
	$(SIGNATURES)

Returns description of a variable.
"""
description(v :: VarInfo{T}) where T = v.description;

"""
	$(SIGNATURES)

Returns the bounds of the variable. These are `nothing` if not set.
"""
bounds(v :: VarInfo{T}) where T = get_options(v, (:lb, :ub));

"""
	$(SIGNATURES)

Returns `true` if `size(vi) == ()`. For example, this is the case for `Real`, but also for user defined types.
"""
isscalar(vi :: VarInfo{T}) where T =
    size(vi) == ();

isequal_common(v1 :: VarInfo{T1}, v2 :: VarInfo{T2}) where {T1, T2} = false;

isequal_common(v1 :: VarInfo{T}, v2 :: VarInfo{T}) where T =
    (var_name(v1) == var_name(v2))  &&
    (size(v1) == size(v2));

isequal_options(v1 :: VarInfo{T}, v2 :: VarInfo{T}) where T =
    true; # stub +++++

    """
	$(SIGNATURES)

Initialize an array (or scalar) with zeros. If `VarInfo` does not contain size, return empty array.
"""
function init_values(vi :: VarInfo{T}) where T
    if isscalar(vi)
        return zero(T);
    elseif has_option(vi, :size)
        return zeros(eltype(T), size(vi));
    else
        error("Size is required");
        # return Array{T, 1}();
    end
end

## ----------  Checking VarInfo

function check_var_info(vi :: VarInfo{T}; silent :: Bool = true) where T
    isValid = true;
    isValid = isValid  &&  check_bounds(vi; silent = silent);
    isValid = isValid  &&  check_size(vi; silent = silent);
    return isValid
    # more checks +++++
end

function check_bounds(vi :: VarInfo{T}; silent :: Bool = true) where T
    isValid = true;
    for bnd in (:lb, :ub)
        isValid = isValid  &&  check_one_bound(vi, bnd; silent = silent);
    end
    return isValid
end

function check_one_bound(vi :: VarInfo{T}, bnd :: Symbol; silent :: Bool = true) where T
    isValid = true;
    if has_option(vi, bnd)
        b = get_option(vi, bnd);
        if !isa(b, eltype(T))
            isValid = false;
            !silent  &&  @warn "Invalid eltype of $bnd:  $(typeof(b))"
        end
    end
    return isValid
end

function check_size(vi :: VarInfo{T}; silent :: Bool = true) where T
    isValid = true;
    if has_option(vi, :size)
        sz = get_option(vi, :size);
        isValid = isValid && isa(sz, Tuple);
        if T <: AbstractArray
            isValid = isValid && (length(sz) == ndims(T));
        end
    else
        error("Size is required");
    end
    return isValid
end

## --------------  Checking data

"""
	$(SIGNATURES)

Check an array against the constraints implied by a `ArrayInfo`.

# Arguments
- `silent`: governs whether `@warn` messages are shown.
"""
function check_var(vi :: VarInfo{T},  x :: T2;
    silent :: Bool = true) where {T, T2}

    if T != T2
        !silent  &&  @warn "Type mismatch for $vi: $T vs $T2";
        return false;
    end
    isValid = check_bounds(vi, x; silent = silent)  &&
        check_size(vi, x; silent = silent) &&
        check_dim_sums(vi, x; silent = silent)  &&
        check_grand_sum(vi, x; silent = silent);
    return isValid
end

function check_bounds(vi :: VarInfo{T},
    x :: T;
    silent :: Bool = true) where T

    isValid = true;
    lb, ub = bounds(vi);
    if !isnothing(lb)
        if !all_at_least(x, lb)
            isValid = false;
            silent  ||  @warn "$vi below lower bound";
        end
    end
    if !isnothing(ub)
        if !all_at_most(x, ub)
            isValid = false;
            silent  ||  @warn "$vi above upper bound";
        end
    end
    return isValid
end

function check_size(vi :: VarInfo{T}, x :: T;
    silent :: Bool = true) where T

    isValid = true;
    if has_option(vi, :size)
        sz = get_option(vi, :size);
        if size(x) != sz
            isValid = false;
            @warn "Wrong size $vi:  $(size(x))";
        end
    else
        error("Size is required");
    end
    return isValid
end

function check_dim_sums(vi :: VarInfo{T}, x :: T;
    silent :: Bool = true, atol = 1e-6) where T

    has_option(vi, :dimSums)  ||  return true;

    isValid = true;
    for iDim = 1 : ndims(T)
        isValid = isValid  &&
            check_dim_sum(vi, x, iDim; silent = silent, atol = atol);
    end
    return isValid
end

function check_dim_sum(vi :: VarInfo{T}, x :: T, iDim :: Integer;
    silent :: Bool = true, atol = 1e-6) where T

    dimSums = get_option(vi, :dimSums);
    if isnothing(dimSums[iDim])
        isValid = true;
    elseif eltype(T) <: Real
        isValid = all(isapprox.(sum(x; dims = iDim), dimSums[iDim], atol = atol));
    elseif eltype(T) <: Integer
        isValid = all(isequal.(sum(x; dims = iDim), dimSums[iDim]));
    end
    return isValid
end

function check_grand_sum(vi :: VarInfo{T}, x :: T;
    silent :: Bool = true, atol = 1e-6) where T

    has_option(vi, :grandSum)  ||  return true;
    if eltype(T) <: Integer
        isValid = isequal(sum(x), get_option(vi, :grandSum));
    elseif eltype(T) <: Real
        isValid = isapprox(sum(x), get_option(vi, :grandSum), atol = atol);
    else
        isValid = true;
    end
    return isValid
end

## ---------  Testing

make_test_var_infos(rng :: AbstractRNG) =
    make_test_var_infos(6, rng);

function make_test_var_infos(nVars, rng :: AbstractRNG)
    eltypeV = [Float64, Int, Float32, UInt8, Int16, Float16];
    szV = [(), (3,), (2,3), (4,3,2), (), (2,)];

    varMeta = Dict{Symbol, VarInfo}();
    for iVar = 1 : nVars
        varName = Symbol("var$iVar");
        lb = convert(eltypeV[iVar], iVar);
        ub = convert(eltypeV[iVar], iVar + 5);
        if szV[iVar] == ()
            T = eltypeV[iVar];
        else
            T = Array{eltypeV[iVar], length(szV[iVar])};
        end
        varMeta[varName] = make_test_var_info(varName, T, szV[iVar], rng);
    end
    return varMeta
end

"""
	$(SIGNATURES)

Make `VarInfo` for testing.
"""
function make_test_var_info(
        vName :: Symbol,
        T :: DataType,
        sz,
        rng :: AbstractRNG)

    eType = eltype(T);
    iVar = rand(rng, 1 : 30);
    lb = convert(eType, iVar);
    ub = convert(eType, lb + 5);
    if T <: Real
        vi = VarInfo(vName, "Descr $vName", T; lb = lb, ub = ub);
    elseif T <: Array
        N = ndims(T);
        if isnothing(sz)
            sz = Tuple(rand(rng, 1:5, N));
        end
        vi = VarInfo(vName, "Descr $vName", T; lb = lb, ub = ub, size = sz);
    else
        error("Invalid $T")
    end
    @assert check_var_info(vi; silent = false);
    return vi
end

function make_test_values(d :: Dict{Symbol, VarInfo}, rng :: AbstractRNG)
    values = Dict{Symbol, Any}();
    for (vName, vInfo) in d
        values[vName] = make_test_values(vInfo, rng);
    end
    return values
end

"""
	$(SIGNATURES)

Make values for a `VarInfo` for testing.
"""
function make_test_values(vi :: VarInfo{T},
    rng :: AbstractRNG) where T

    eType = eltype(T);
    lb = get_option(vi, :lb);
    if isnothing(lb)
        lb = eType(1);
    end
    ub = get_option(vi, :ub);
    if isnothing(ub)
        ub = lb + eType(5);
    end
    if T <: Real
        sz = ();
    else
        sz = get_option(vi, :size);
        if isnothing(sz)
            sz = rand(rng, 1:5, ndims(T)...);
        end
    end

    if eType <: Integer
        # This works for scalars when `sz == ()`
        valueV = rand(rng, lb : ub, sz...);
    elseif eType <: AbstractFloat
        valueV = lb .+ (ub - lb) .* rand(rng, eType, sz...);
    else
        error("Invalid type");
    end
    # Make sure scalars are scalar
    if T <: Real
        valueV = valueV[1];
    end
    # Ensure dim sums
    # Ensure grand sum
    @assert check_var(vi, valueV; silent = false)
    return valueV
end

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


text:
function removeAgent_cpu!(pos,keepList_,removeList_,remV_)
    keepList_[pos] = 0
    removeList_Pos_ = Threads.atomic_add!(remV_,1) + 1
    removeList_[removeList_Pos_] = pos
end

function removeAgent_gpu!(pos,keepList_,removeList_,remV_)
    keepList_[pos] = 0
    removeList_Pos_ = CUDA.atomic_add!(CUDA.pointer(remV_,Int32(1)),Int32(1)) + 1
    removeList_[removeList_Pos_] = pos
end

function removeDeadAgents_cpu!(remV_,keepList_,removeList_,localV,nLocal,identityV,nIdentity,localVCopy,nLocalC,identityVCopy,nIdentityC)

    Threads.@threads for ic1_ in 1:remV_[]
        oldPos_ = keepList_[ic1_]
        newPos_ = removeList_[ic1_]
        if oldPos_ > 0
            for ic2_ in 1:nLocal
                localV[newPos_, ic2_] = localV[oldPos_, ic2_]
            end
            for ic2_ in 1:nIdentity
                identityV[newPos_, ic2_] = identityV[oldPos_, ic2_]
            end
            for ic2_ in 1:nLocalC
                localVCopy[newPos_, ic2_] = localVCopy[oldPos_, ic2_]
            end
            for ic2_ in 1:nIdentityC
                identityVCopy[newPos_, ic2_] = identityVCopy[oldPos_, ic2_]
            end
        end
    end

    return
end

function removeDeadAgents_cpu!(remV_,keepList_,removeList_,localV,nLocal,identityV,nIdentity,localVCopy,nLocalC)

    Threads.@threads for ic1_ in 1:remV_[]
        oldPos_ = keepList_[ic1_]
        newPos_ = removeList_[ic1_]
        if oldPos_ > 0
            for ic2_ in 1:nLocal
                localV[newPos_, ic2_] = localV[oldPos_, ic2_]
            end
            for ic2_ in 1:nIdentity
                identityV[newPos_, ic2_] = identityV[oldPos_, ic2_]
            end
            for ic2_ in 1:nLocalC
                localVCopy[newPos_, ic2_] = localVCopy[oldPos_, ic2_]
            end
        end
    end

    return
end

function removeDeadAgents_cpu!(remV_,keepList_,removeList_,localV,nLocal,identityV,nIdentity,identityVCopy::Array{<:Int,2},nIdentityC)

    Threads.@threads for ic1_ in 1:remV_[]
        oldPos_ = keepList_[ic1_]
        newPos_ = removeList_[ic1_]
        if oldPos_ > 0
            for ic2_ in 1:nLocal
                localV[newPos_, ic2_] = localV[oldPos_, ic2_]
            end
            for ic2_ in 1:nIdentity
                identityV[newPos_, ic2_] = identityV[oldPos_, ic2_]
            end
            for ic2_ in 1:nIdentityC
                identityVCopy[newPos_, ic2_] = identityVCopy[oldPos_, ic2_]
            end
        end
    end

    return
end

function removeDeadAgents_cpu!(remV_,keepList_,removeList_,localV,nLocal,identityV,nIdentity)

    Threads.@threads for ic1_ in 1:remV_[]
        oldPos_ = keepList_[ic1_]
        newPos_ = removeList_[ic1_]
        if oldPos_ > 0
            for ic2_ in 1:nLocal
                localV[newPos_, ic2_] = localV[oldPos_, ic2_]
            end
            for ic2_ in 1:nIdentity
                identityV[newPos_, ic2_] = identityV[oldPos_, ic2_]
            end
        end
    end

    return
end

function removeDeadAgents_gpu!(remV_,keepList_,removeList_,localV,nLocal,identityV,nIdentity,localVCopy,nLocalC,identityVCopy,nIdentityC)

    index_ = (threadIdx().x) + (blockIdx().x - 1) * blockDim().x
    stride_ = blockDim().x * gridDim().x

    for ic1_ in index_:stride_:remV_[1]
        oldPos_ = keepList_[ic1_]
        newPos_ = removeList_[ic1_]
        if oldPos_ > 0
            for ic2_ in 1:nLocal
                localV[newPos_, ic2_] = localV[oldPos_, ic2_]
            end
            for ic2_ in 1:nIdentity
                identityV[newPos_, ic2_] = identityV[oldPos_, ic2_]
            end
            for ic2_ in 1:nLocalC
                localVCopy[newPos_, ic2_] = localVCopy[oldPos_, ic2_]
            end
            for ic2_ in 1:nIdentityC
                identityVCopy[newPos_, ic2_] = identityVCopy[oldPos_, ic2_]
            end
        end
    end

    return nothing
end

function removeDeadAgents_gpu!(remV_,keepList_,removeList_,localV,nLocal,identityV,nIdentity,localVCopy::CuDeviceMatrix{<:AbstractFloat, 1},nLocalC)

    index_ = (threadIdx().x) + (blockIdx().x - 1) * blockDim().x
    stride_ = blockDim().x * gridDim().x

    for ic1_ in index_:stride_:remV_[1]
        oldPos_ = keepList_[ic1_]
        newPos_ = removeList_[ic1_]
        if oldPos_ > 0
            for ic2_ in 1:nLocal
                localV[newPos_, ic2_] = localV[oldPos_, ic2_]
            end
            for ic2_ in 1:nIdentity
                identityV[newPos_, ic2_] = identityV[oldPos_, ic2_]
            end
            for ic2_ in 1:nLocalC
                localVCopy[newPos_, ic2_] = localVCopy[oldPos_, ic2_]
            end
        end
    end

    return
end

function removeDeadAgents_gpu!(remV_,keepList_,removeList_,localV,nLocal,identityV,nIdentity,identityVCopy::CuDeviceMatrix{<:Int, 1},nIdentityC)

    index_ = (threadIdx().x) + (blockIdx().x - 1) * blockDim().x
    stride_ = blockDim().x * gridDim().x

    for ic1_ in index_:stride_:remV_[1]
        oldPos_ = keepList_[ic1_]
        newPos_ = removeList_[ic1_]
        if oldPos_ > 0
            for ic2_ in 1:nLocal
                localV[newPos_, ic2_] = localV[oldPos_, ic2_]
            end
            for ic2_ in 1:nIdentity
                identityV[newPos_, ic2_] = identityV[oldPos_, ic2_]
            end
            for ic2_ in 1:nIdentityC
                identityVCopy[newPos_, ic2_] = identityVCopy[oldPos_, ic2_]
            end
        end
    end

    return
end

function removeDeadAgents_gpu!(remV_,keepList_,removeList_,localV,nLocal,identityV,nIdentity)

    index_ = (threadIdx().x) + (blockIdx().x - 1) * blockDim().x
    stride_ = blockDim().x * gridDim().x

    for ic1_ in index_:stride_:remV_[1]
        oldPos_ = keepList_[ic1_]
        newPos_ = removeList_[ic1_]
        if oldPos_ > 0
            for ic2_ in 1:nLocal
                localV[newPos_, ic2_] = localV[oldPos_, ic2_]
            end
            for ic2_ in 1:nIdentity
                identityV[newPos_, ic2_] = identityV[oldPos_, ic2_]
            end
        end
    end

    return
end

"""
function addEventDeath_(code::Expr,p::Program_,platform::String)

Generate the functions for division events.
"""
function addEventRemoveAgent_(code::Expr,p::Program_,platform::String)

    #Add remove event if declared
    if inexpr(code,:removeAgent)

        name = Meta.parse(string("removeAgent_",platform,"!"))
        code = postwalk(x -> @capture(x,removeAgent()) ? :(AgentBasedModels.$name(ic1_,keepList_,removeList_,remV_)) : x, code)

        if inexpr(code,:removeAgent) #Check correct declaration
            error("removeAgent declared with wrong number of parameters. It has to be declared as removeAgent())")
        end

        #Add the keep list initialization at the beginning
        codeinit = quote keepList_[ic1_] = ic1_ end
        push!(codeinit.args,code.args...)
        code = codeinit

        #Add to code the removing algorithm
        a = [:remV_,:keepList_,
                :removeList_,
                :localV,length(p.agent.declaredSymbols["Local"]),
                :identityV,length(p.agent.declaredSymbols["Identity"])]
        if !isempty(p.update["Local"])
            push!(a, :localVCopy,length(p.update["Local"]))
        end
        if !isempty(p.update["Identity"])
            push!(a, :identityVCopy, length(p.update["Identity"]))
        end
        if platform == "cpu"
            update = quote
                if N > 0
                    sort!(keepList_,rev=true)
                    @platformAdapt AgentBasedModels.removeDeadAgents_cpu!($(a...))
                    N -= remV_[]
                    remV_[] = 0
                    keepList_ .= 0
                end
            end
        elseif platform == "gpu"
            update = quote
                if N > 0
                    sort!(keepList_,rev=true)
                    @platformAdapt AgentBasedModels.removeDeadAgents_gpu!($(a...))
                    N -= Core.Array(remV_)[1]
                    remV_ .= 0
                    keepList_ .= 0
                end
            end
        end
        push!(p.execInloop.args,update)

        # Add to program
        if platform == "cpu"
            push!(p.declareVar.args,
                :(begin
                    remV_ = Threads.Atomic{$INT}()
                    keepList_ = zeros($INT,nMax)
                    removeList_ = zeros($INT,nMax)
                end)
            )
        elseif platform == "gpu"
            push!(p.declareVar.args,
                :(begin
                    remV_ = CUDA.zeros($INTCUDA,1)
                    keepList_ = zeros($INTCUDA,nMax)
                    removeList_ = zeros($INTCUDA,nMax)
                end)
            )
        end

        push!(p.args,:remV_,:keepList_,:removeList_)

    end

    return code
end


text:
module Forcing

export ModelForcing, SimpleForcing, ParameterizedForcing, Relaxation, GaussianMask, LinearTarget

using Oceananigans.Fields

include("simple_forcing.jl")
include("model_forcing.jl")
include("parameterized_forcing.jl")
include("relaxation.jl")

end


text:
using .Gadfly

using Colors
using Pkg.Artifacts
using CSV
using DataFrames
using Statistics

const NeXLPalette =
    convert.(
        RGB{Colors.N0f8},
        distinguishable_colors(
            66,
            [RGB(253 / 255, 253 / 255, 241 / 255), RGB(0, 0, 0), colorant"DodgerBlue4"],
            transform = deuteranopic,
        )[3:end],
    )

const NeXLColorblind = NeXLPalette

"""
    Gadfly.plot(transitions::AbstractVector{Transition}; mode=:Energy|:Weight, palette=NeXLPalette)

Plot either the :Energies or :Weights associated with the specified transitions over the range of supported elements.
"""
function Gadfly.plot(
    transitions::AbstractVector{Transition};
    mode = :Energy,
    palette = NeXLPalette,
)
    if mode == :Energy
        plotXrayEnergies(transitions)
    elseif mode == :Weight
        plotXrayWeights(transitions)
    end
end

function plotXrayEnergies(transitions::AbstractVector{Transition}; palette = NeXLPalette)
    layers, names = [], String[]
    colors = distinguishable_colors(
        length(transitions) + 2,
        Color[RGB(253 / 255, 253 / 255, 241 / 255), RGB(0, 0, 0)],
    )
    for (i, tr) in enumerate(transitions)
        x, y = [], []
        for elm in element.(eachelement())
            if has(elm, tr)
                push!(x, z(elm))
                push!(y, energy(characteristic(elm, tr)))
            end
        end
        if !isempty(x)
            push!(names, repr(tr))
            append!(
                layers,
                Gadfly.layer(
                    x = x,
                    y = y,
                    Geom.point,
                    Gadfly.Theme(default_color = colors[i+2]),
                ),
            )
        end
    end
    Gadfly.plot(
        layers...,
        Gadfly.Guide.title("Characteristic X-ray Energies"),
        Gadfly.Guide.manual_color_key("Type", names, color = colors[3:end]),
        Gadfly.Guide.xlabel("Atomic Number"),
        Guide.ylabel("Energy (eV)"),
        Gadfly.Coord.cartesian(xmin = eachelement()[1], xmax = eachelement()[end]),
    )
end

function plotXrayWeights(transitions::AbstractVector{Transition}, schoonjan::Bool = false)
    layers, names = [], String[]
    colors = distinguishable_colors(
        length(transitions) + 2,
        Color[RGB(253 / 255, 253 / 255, 241 / 255), RGB(0, 0, 0)],
    )
    for (i, tr) in enumerate(transitions)
        x, y = [], []
        for elm in element.(eachelement())
            if has(elm, tr)
                push!(x, z(elm))
                push!(
                    y,
                    schoonjan ? normweight(characteristic(elm, tr)) :
                    strength(characteristic(elm, tr)),
                )
            end
        end
        if !isempty(x)
            push!(names, repr(tr))
            append!(
                layers,
                Gadfly.layer(
                    x = x,
                    y = y,
                    Geom.point,
                    Gadfly.Theme(default_color = colors[i+2]),
                ),
            )
        end
    end
    if schoonjan  # Compare to Shoonjan's xraylib on GitHub
        function parse2(tr)
            ss1, ss2 = missing, missing
            try
                if tr[1] == 'K'
                    ss1 = n"K1"
                    ss2 =
                        (tr[2] ≠ 'O' && tr[2] ≠ 'P') ? parse(SubShell, tr[2:end]) :
                        (tr[2] == 'O' ? n"O3" : n"P3")
                else
                    ss1 = parse(SubShell, tr[1:2])
                    ss2 =
                        (tr[3] ≠ 'O' && tr[3] ≠ 'P') ? parse(SubShell, tr[3:end]) :
                        (tr[3] == 'O' ? n"O3" : n"P3")
                end
                return exists(ss1, ss2) ? Transition(ss1, ss2) : missing
            catch
                return missing
            end
        end
        artpath = downloadschoonjan()
        csv = CSV.read(
            joinpath(artpath, "radrate.dat"),
            DataFrame,
            delim = ' ',
            ignorerepeated = true,
            header = 0,
        )
        insertcols!(csv, 3, Column2p = parse2.(csv[!, :Column2]))
        filter!(r -> (!ismissing(r[:Column2p])) && (r[:Column2p] in transitions), csv)
        insertcols!(
            csv,
            3,
            Column2pp = CategoricalArray(map(n -> "Schoonj $n", csv[!, :Column2p])),
        )
        append!(layers, layer(csv, x = :Column1, y = :Column3, color = :Column2pp))
    end
    Gadfly.plot(
        layers...,
        Gadfly.Guide.title("Characteristic X-ray Weights"),
        Gadfly.Guide.manual_color_key("Type", names, colors[3:end]),
        Gadfly.Guide.xlabel("Atomic Number"),
        Guide.ylabel("Weight"),
        Gadfly.Coord.cartesian(xmin = eachelement()[1], xmax = eachelement()[end]),
    )
end

"""
    Gadfly.plot(sss::AbstractVector{SubShell}, mode=:EdgeEnergy|:FluorescenceYield; palette=NeXLPalette)

Plot the edge energies/fluorescence yields associated with the specified vector of SubShell objects.
"""
function Gadfly.plot(sss::AbstractVector{SubShell}, mode = :EdgeEnergy)
    if mode == :FluorescenceYield
        plotFluorescenceYield(sss::AbstractVector{SubShell})
    else
        plotEdgeEnergies(sss)
    end
end

# Download the data once per installation...
function downloadschoonjan()
    artifacts_toml = joinpath(@__DIR__, "Artifacts.toml")
    hash = artifact_hash("schoonjan", artifacts_toml)
    if hash === nothing || !artifact_exists(hash)
        hash = create_artifact() do artifact_dir
            # We create the artifact by simply downloading a few files into the new artifact directory
            download(
                "https://github.com/tschoonj/xraylib/blob/master/data/fluor_yield.dat?raw=true",
                joinpath(artifact_dir, "fluor_yield.dat"),
            )
            download(
                "https://github.com/tschoonj/xraylib/blob/master/data/radrate.dat?raw=true",
                joinpath(artifact_dir, "radrate.dat"),
            )
        end
        bind_artifact!(artifacts_toml, "schoonjan", hash)
        @info "Downloaded Schoonjan data into Archive."
    end
    return artifact_path(hash)
end

function plotFluorescenceYield(sss::AbstractVector{SubShell}, schoonjan::Bool = false)
    layers, names = [], String[]
    colors = distinguishable_colors(
        length(sss) + 2,
        Color[RGB(253 / 255, 253 / 255, 241 / 255), RGB(0, 0, 0)],
    )
    for (i, sh) in enumerate(sss)
        x, y = [], []
        for elm in element.(eachelement())
            if has(elm, sh)
                push!(x, z(elm))
                push!(y, fluorescenceyield(atomicsubshell(elm, sh), CullenEADL))
            end
        end
        if !isempty(x)
            push!(names, repr(sh))
            append!(
                layers,
                Gadfly.layer(
                    x = x,
                    y = y,
                    Geom.point,
                    Gadfly.Theme(default_color = colors[i+2]),
                ),
            )
        end
    end
    if schoonjan
        artpath = downloadschoonjan()
        csv = CSV.read(
            joinpath(artpath, "fluor_yield.dat"),
            DataFrame,
            delim = ' ',
            ignorerepeated = true,
            header = 0,
        )
        sssname = [repr(ss) for ss in sss]
        filter!(r -> r[:Column2] in sssname, csv)
        insertcols!(
            csv,
            3,
            Column2p = CategoricalArray(map(n -> "Schoonj $n", csv[!, :Column2])),
        )
        append!(layers, layer(csv, x = :Column1, y = :Column3, color = :Column2p))
    end
    Gadfly.plot(
        layers...,
        Gadfly.Guide.title("Fluourescence Yield"),
        Gadfly.Guide.manual_color_key("Type", names, colors[3:end]),
        Gadfly.Guide.xlabel("Atomic Number"),
        Guide.ylabel("Yield (Fractional)"),
        Scale.y_log10(maxvalue = 1.0),
        Gadfly.Coord.cartesian(xmin = eachelement()[1], xmax = eachelement()[end]),
    )
end

function plotEdgeEnergies(sss::AbstractVector{SubShell})
    layers, names = [], String[]
    colors = distinguishable_colors(
        length(sss) + 2,
        Color[RGB(253 / 255, 253 / 255, 241 / 255), RGB(0, 0, 0)],
    )
    for (i, sh) in enumerate(sss)
        x, y = [], []
        for elm in element.(eachelement())
            if has(elm, sh)
                push!(x, z(elm))
                push!(y, energy(atomicsubshell(elm, sh)))
            end
        end
        if !isempty(x)
            push!(names, repr(sh))
            append!(
                layers,
                Gadfly.layer(
                    x = x,
                    y = y,
                    Geom.point,
                    Gadfly.Theme(default_color = colors[i+2]),
                ),
            )
        end
    end
    Gadfly.plot(
        layers...,
        Gadfly.Guide.title("Atomic Sub-Shell Energies"),
        Gadfly.Guide.manual_color_key("Type", names, colors[3:end]),
        Gadfly.Guide.xlabel("Atomic Number"),
        Guide.ylabel("Edge Energy (eV)"),
        Gadfly.Coord.cartesian(xmin = eachelement()[1], xmax = eachelement()[end]),
    )
end

"""
    compareMACs(elm::Element; palette=NeXLPalette)

Plot a comparison of the FFAST and Heinrich MAC tabulations for the specified Element.
"""
function compareMACs(elm::Element; palette = NeXLPalette)
    l1 = layer(
        ev -> log10(mac(elm, ev, FFASTDB)),
        100.0,
        20.0e3,
        Geom.line,
        Gadfly.Theme(default_color = palette[1]),
    )
    l2 = layer(
        ev -> log10(mac(elm, ev, DTSA)),
        100.0,
        20.0e3,
        Geom.line,
        Gadfly.Theme(default_color = palette[2]),
    )
    Gadfly.plot(
        l1,
        l2,
        Gadfly.Guide.title("MAC - $elm"),
        Gadfly.Guide.manual_color_key("Type", ["Default/FFAST", "Heinrich"], palette[1:2]),
        Gadfly.Guide.xlabel("Energy (eV)"),
        Guide.ylabel("log₁₀(MAC (cm²/g))"),
        Gadfly.Coord.cartesian(xmin = 0.0, xmax = 20.0e3),
    )
end

"""
    plot(alg::Type{<:NeXLAlgorithm}, elm::Union{Element,Material}; palette = NeXLPalette, xmax=20.0e3)

Plot a MAC tabulations for the specified Element or Material.
"""
function Gadfly.plot(
    alg::Type{<:NeXLAlgorithm},
    elm::Union{Element,Material};
    palette = NeXLPalette,
    xmax = 20.0e3,
)
    l1 = layer(
        ev -> log10(mac(elm, ev, alg)),
        100.0,
        xmax,
        Geom.line,
        Gadfly.Theme(default_color = palette[1]),
    )
    Gadfly.plot(
        l1,
        Gadfly.Guide.title("MAC - $(name(elm))"),
        Gadfly.Guide.xlabel("Energy (eV)"),
        Guide.ylabel("log₁₀(MAC (cm²/g))"),
        Gadfly.Coord.cartesian(xmin = 0.0, xmax = xmax),
    )
end

"""
    compareMACs(elm::Element; palette=NeXLPalette)

Plot a comparison of the FFAST and Heinrich MAC tabulations for the specified Element or Material.
"""
function Gadfly.plot(
    alg::Type{<:NeXLAlgorithm},
    elms::AbstractVector;
    palette = NeXLPalette,
    xmin = 100.0,
    xmax = 20.0e3,
)
    layers, colors, names = Layer[], Color[], String[]
    for (i, elm) in enumerate(elms)
        append!(
            layers,
            layer(
                ev -> log10(mac(elm, ev, alg)),
                xmin,
                xmax,
                Geom.line,
                Gadfly.Theme(default_color = palette[i]),
            ),
        )
        push!(colors, palette[i])
        push!(names, name(elm))
    end
    Gadfly.plot(
        layers...,
        Gadfly.Guide.xlabel("Energy (eV)"),
        Guide.ylabel("log₁₀(MAC (cm²/g))"),
        Gadfly.Guide.manual_color_key("Material", names, colors),
        Gadfly.Coord.cartesian(xmin = max(0.0, xmin - 100.0), xmax = xmax),
    )
end

function Gadfly.plot(
    mats::AbstractVector{<:Material};
    known::Union{Material,Missing} = missing,
    delta::Bool = false,
    label::AbstractString = "Material",
    palette = NeXLPalette,
)
    allelms = collect(union(map(keys, mats)...))
    xs = [name(mat) for mat in mats]
    layers, names, colors = Layer[], String[], RGB{Float64}[]
    if ismissing(known)
        known = material(
            "the Mean",
            Dict{Element,Float64}(
                elm => mean([value(mat[elm]) for mat in mats]) for elm in allelms
            ),
        )
    end
    for (i, elm) in enumerate(allelms)
        if delta
            append!(
                layers,
                layer(
                    x = xs,
                    y = [value(mat[elm]) - known[elm] for mat in mats],
                    ymin = [value(mat[elm]) - σ(mat[elm]) - known[elm] for mat in mats],
                    ymax = [value(mat[elm]) + σ(mat[elm]) - known[elm] for mat in mats],
                    Gadfly.Theme(default_color = palette[i]),
                    Geom.errorbar,
                    Geom.point,
                ),
            )
        else
            append!(
                layers,
                layer(
                    x = xs,
                    y = [value(mat[elm]) for mat in mats],
                    ymin = [value(mat[elm]) - σ(mat[elm]) for mat in mats],
                    ymax = [value(mat[elm]) + σ(mat[elm]) for mat in mats],
                    Gadfly.Theme(default_color = palette[i]),
                    Geom.errorbar,
                    Geom.point,
                ),
            )
        end
        push!(names, name(elm))
        push!(colors, palette[i])
    end
    lighten(col) = weighted_color_mean(0.2, RGB(col), colorant"white")
    if delta
        plot(
            layers...,
            Guide.ylabel("Δ(Mass Fraction)"),
            Guide.xlabel(label),
            Guide.manual_color_key("Element", names, Colorant[colors...]),
            Guide.title("Difference from $(known)"),
            Geom.hline(color = "black"),
            yintercept = [0.0],
        )
    else
        plot(
            layers...,
            Guide.ylabel("Mass Fraction"),
            Guide.xlabel(label), #
            Guide.manual_color_key("Element", names, Colorant[colors...]),
            yintercept = [known[elm] for elm in allelms],
            Geom.hline(color = [lighten(col) for col in colors], style = :dash),
        )
    end
end

function plot2(
    mats::AbstractVector{<:Material};
    label::AbstractString = "Material",
    palette = NeXLPalette,
)
    allelms = sort(convert(Vector{Element}, collect(union(map(keys, mats)...))))
    elmcol = Dict(elm => palette[i] for (i, elm) in enumerate(allelms))
    xs, ymin, ymax, ygroups, colors = String[], Float64[], Float64[], Element[], Color[]
    for mat in mats
        append!(xs, [name(mat) for elm in keys(mat)])
        append!(ymin, [value(mat[elm]) - σ(mat[elm]) for elm in keys(mat)])
        append!(ymax, [value(mat[elm]) + σ(mat[elm]) for elm in keys(mat)])
        append!(colors, [elmcol[elm] for elm in keys(mat)])
        append!(ygroups, collect(keys(mat)))
    end
    plot(
        x = xs,
        ymin = ymin,
        ymax = ymax,
        color = colors,
        ygroup = ygroups,
        Geom.subplot_grid(Geom.errorbar, free_y_axis = true),
        Scale.ygroup(labels = elm -> symbol(elm), levels = allelms),
        Guide.xlabel(label),
        Guide.ylabel("Mass Fraction by Element"),
    )
end


text:
#
# Correctness Tests
#

using Calculus
using Compat
using Compat.Test
using Compat.LinearAlgebra

tests = ["finite_difference",
         "derivative",
         "check_derivative",
         "symbolic",
         "deparse"]

println("Running tests:")

for t in tests
    println(" * $(t)")
    include("$(t).jl")
end


text:
"""
    TransferEntropyCausalityTest{N}

The supertype of all abstract and composite types representing a transfer
entropy causality test applied to `N` different prediction lags.

Concrete subtypes are

- [`TransferOperatorGridTest`](@ref)
- [`VisitationFrequencyTest`](@ref)
- [`ApproximateSimplexIntersectionTest`](@ref)
- [`ExactSimplexIntersectionTest`](@ref)
"""
abstract type TransferEntropyCausalityTest{N} <: EntropyBasedCausalityTest end

export TransferEntropyCausalityTest


text:
module TopOptProblems

using Ferrite, StaticArrays, LinearAlgebra
using SparseArrays, Setfield, Requires
using ..TopOpt.Utilities
using ..TopOpt: PENALTY_BEFORE_INTERPOLATION
using ..Utilities: @forward_property
import Distributions

using VTKDataTypes

import Ferrite: assemble!

abstract type AbstractTopOptProblem end

include("utils.jl")
include("grids.jl")
include("metadata.jl")
include("problem_types.jl")
include("multiload.jl")
include("elementmatrix.jl")
include("matrices_and_vectors.jl")
include("elementinfo.jl")
include("assemble.jl")
include("buckling.jl")

include(joinpath("IO", "IO.jl"))
using .InputOutput

include("Visualization/Visualization.jl")
using .Visualization

export RayProblem, PointLoadCantilever, HalfMBB, LBeam, TieBeam, InpStiffness, StiffnessTopOptProblem, AbstractTopOptProblem, GlobalFEAInfo, ElementFEAInfo, YoungsModulus, assemble, assemble_f!, RaggedArray, ElementMatrix, rawmatrix, bcmatrix, save_mesh, RandomMagnitude, MultiLoad

end # module