text:
module InformationGeometry

# using Reexport
# @reexport
using LinearAlgebra, Random, Distributions, DataFrames

using Distributed, StaticArrays, SparseArrays
using OrdinaryDiffEq, DiffEqCallbacks, BoundaryValueDiffEq
using ModelingToolkit, Symbolics, DataInterpolations
using ForwardDiff, DerivableFunctions
using BenchmarkTools, LsqFit, Optim, Measurements, HCubature
using SpecialFunctions, Tullio, DataFrames, Roots, Combinatorics
using LibGEOS, Sobol, ProgressMeter, Suppressor
using RecipesBase, Plots
using TreeViews

import DataInterpolations: AbstractInterpolation
import SciMLBase: AbstractODESolution, AbstractODEFunction

import SciMLBase: remake
export remake

import DerivableFunctions: suff, MaximalNumberOfArguments, KillAfter, GetArgLength, Builder
import DerivableFunctions: DerivableFunction, DFunction, EvalF, EvaldF, EvalddF
import DerivableFunctions: GetDeriv, GetGrad, GetJac, GetHess, GetMatrixJac
import DerivableFunctions: GetGrad!, GetJac!, GetHess!, GetMatrixJac!, diff_backends

# include("NumericalTools/Differentiation.jl")
export diff_backends
export GetDeriv, GetGrad, GetJac, GetMatrixJac, GetHess, GetDoubleJac
export GetGrad!, GetHess!, GetJac!, GetMatrixJac!

############## TODOs ##############

# How can kwargs be consistently added to FindConfBoundary(), FindMLE(), ... such that tols are not messed up?

#### Functionality

# Employ ShadowTheatre to project simultaneous x-θ-confidence regions to M.
# Geodesic coordinates - use geodesic radius and angles to specify parameter configurations.
# ADD TESTS: General Parameter Space / Model Transformations: return new model where `log` has been applied to some / all θ compoents
# Allow for inplace models and dmodels -- potentially very large performance benefits
# ADD TESTS: Add CompositeDataSet to master
# Test F-Boundaries in more detail
# Integrate along structural unidentifiability -- Switching ODE functions mid-integration: https://diffeq.sciml.ai/stable/basics/faq/
# Use InDomain function specified for ModelMaps for constrained optimization.
# Symbolic Fisher? Christoffel (from 2nd derivatives)?
# Use vector-callback with domain hypercube for continuous boundary detection

# Compute EigenFlow of Fisher Metric -> Should be transformed to straight lines (in coordinates) by Decorrelation Transformation
# Employ MultistartOptimization for obtaining the MLE.
# Try to recognize which components of θ are inside an exponential such that `log` leads to linear parametrization -> i.e. compare EmbeddingMatrices

#### Cosmetic / Convenience

# Custom type for confidence boundaries?
# Show distribution types for DataSetExact
# Infer variable names from DataFrames
# IO methods for ModelMaps in general
# Add measures of practical unidentifiability (eigenvalues of Fisher?)

#### Performance / Optimization

# Make OrthVF in-place everywhere
# Optimize GetModel() for ODEs for performance by only executing ODEFunction once
# Use IntegrateND() instead of MonteCarloArea()
# Remove Plots.jl dependency and instead write recipes for VisualizeSols and others
# Change EvaluateAlongGeodesicLength using affine reparametrization to save computations
# Use Cuba.jl for Monte Carlo integration
# Generalize FisherMetric to other error distributions

#### Tests to add

# GenerateInterruptedBoundary
# F-Boundary
# Constructors of ModelMap, DataSet, DataModel, etc
# ..... make full list

#### Miscellaneous

# Add exporting methods and detailed examples to docs
# Improve ConfidenceBands
# Fix FindConfBoundary() and GenerateBoundary() for BigFloat    -- does it work now? Add a test
# Make user-facing keywords (even) more uniform: tol, meth, Big, ...

abstract type AbstractDataSet end
abstract type AbstractDataModel end
abstract type Cuboid end
export AbstractDataSet, AbstractDataModel, Cuboid

import Base: length, rand, BigFloat, in, union, intersect, join, ==
import DataFrames.DataFrame
import Distributions.loglikelihood

BoolArray = Union{BitArray,AbstractArray{<:Bool}}
BoolVector = Union{BitVector,AbstractVector{<:Bool}}

include("Subspaces.jl")
export Plane
export BasisVector, PlaneCoordinates, Shift, IsOnPlane, TranslatePlane, RotatePlane, DecomposeWRTPlane
export DistanceToPlane, ProjectOntoPlane, IsNormalToPlane, MinimizeOnPlane, ParallelPlanes
export HyperPlane, HyperCube
export Inside, in, ConstructCube, CubeWidths, CubeVol, Center, TranslateCube, ResizeCube
export DropCubeDims, FaceCenters, Corners
export PositiveDomain, NegativeDomain, FullDomain, rand, EmbeddedODESolution

include("ModelMaps.jl")
export ModelMap, ModelOrFunction, InformNames, IsInDomain
export Transform, LogTransform, Log10Transform, ExpTransform, Power10Transform, ScaleTransform
export TranslationTransform, LinearTransform, AffineTransform, LinearDecorrelation
export EmbedModelVia, Embedding
export LinearModel, QuadraticModel, ExponentialModel, SumExponentialsModel, PolynomialModel

include("DataStructures/DistributionTypes.jl")
export Dirac, GeneralProduct

include("GeneralDataStructures.jl")
# export HealthyData, HealthyCovariance, CheckModelHealth
export xdata, ydata, xsigma, ysigma, xInvCov, yInvCov, Npoints, xdim, ydim, pdim, DataspaceDim, Data, MLE, LogLikeMLE, WoundX
export xdist, ydist, dist
export Predictor, dPredictor, LogPrior, ConsistentElDims
export MeasureAutoDiffPerformance
export DataDist, SortDataSet, SortDataModel, SubDataSet, SubDataModel, DataFrame, join, length
export MLEinPlane, PlanarDataModel, DetermineDmodel

include("DataStructures/DataSet.jl")
export DataSet

include("DataStructures/DataSetExact.jl")
export DataSetExact, LogLike

include("DataStructures/CompositeDataSet.jl")
export CompositeDataSet

include("DataStructures/GeneralizedDataSet.jl")
export GeneralizedDataSet

include("DataStructures/DataModel.jl")
export DataModel, Prior

include("NumericalTools/Miscellaneous.jl")

include("NumericalTools/Integration.jl")
include("NumericalTools/Optimization.jl")
export GetH, suff, Unpack, Unwind, Windup, ToCols, PromoteStatic, SplitAfter
export ConfAlpha, ConfVol, InvConfVol, ChisqCDF, InvChisqCDF
export KillAfter
export Integrate1D, IntegrateND, IntegrateOverConfidenceRegion, IntegrateOverApproxConfidenceRegion
export LineSearch, MonteCarloArea
export curve_fit, RobustFit, TotalLeastSquares, BlockMatrix

include("ConfidenceRegions.jl")
export likelihood, loglikelihood, Score, WilksCriterion, WilksTest, OrthVF, FindMLE
export AutoScore, AutoMetric
export FindConfBoundary, FCriterion, FTest, FindFBoundary
export GenerateBoundary, ConfidenceRegion, ConfidenceRegions
export GenerateInterruptedBoundary, InterruptedConfidenceRegion
export IsStructurallyIdentifiable, StructurallyIdentifiable
export FisherMetric, GeometricDensity
export ConfidenceRegionVolume, CoordinateVolume
export ExpectedInvariantVolume, GeodesicRadius, CoordinateDistortion, Sensitivity

export EmbeddingMap, EmbeddingMatrix, EmbeddingMap!, EmbeddingMatrix!, Pullback, Pushforward
export AIC, AICc, BIC, ModelComparison, IsLinearParameter, IsLinear, LeastInformativeDirection

export FindConfBoundaryOnPlane, LinearCuboid, IntersectCube, IntersectRegion, MincedBoundaries, ConfidenceBoundary

export ApproxInRegion, ShadowTheatre, CastShadow

include("ProfileLikelihood.jl")
export ProfileLikelihood, PlotProfileTrajectories, InterpolatedProfiles, ProfileBox, PracticallyIdentifiable
export ValInserter, InsertIntoFirst, InsertIntoLast

include("Divergences.jl")
export KullbackLeibler

include("Geodesics.jl")
export ComputeGeodesic, GeodesicLength, GeodesicCrossing, DistanceAlongGeodesic, Endpoints, EvaluateEach
export RadialGeodesics, BoundaryViaGeodesic, GeodesicBetween, GeodesicDistance, GeodesicEnergy, ExponentialMap, LogarithmicMap, KarcherMean, MBAM
# Needs redo:
# export ConstLengthGeodesics, ConstParamGeodesics

include("Curvature.jl")
export ChristoffelSymbol, Riemann, Ricci, RicciScalar
# export Weyl
# Also add Kretschmann, Schouten?

include("DiffEqModels.jl")
export GetModel, ModifyODEmodel

include("SymbolicComputations.jl")
export Optimize, OptimizedDM, InplaceDM, SymbolicModel, SymbolicdModel

include("Plotting.jl")
export FittedPlot, ResidualPlot, PlotScalar, PlotLoglikelihood, Plot2DVF, ResidualSquaredError, PlotEllipses
export Deplanarize, VisualizeSols, VisualizeGeos, VisualizeSolPoints, ConstructAmbientSolution
export ConfidenceBands, ApproxConfidenceBands, PlotConfidenceBands, ConfidenceBandWidth, PredictionEnsemble, PlotMatrix
export EvaluateAlongGeodesic, PlotAlongGeodesic, EvaluateAlongCurve, PlotAlongCurve, PhaseSpacePlot
# export RectangularFacetIndices, RectToTriangFacets, CreateMesh, ToObj, WriteObj

include("Exporting.jl")
export SaveAdaptive, SaveConfidence, SaveGeodesics, SaveDataSet
# export Homogenize, Dehomogenize

include("CustomIO.jl")
# export GeneratedFromAutoDiff, GeneratedFromSymbolic

end # module


text:
"""
The [`Themes`](@ref) module contains APIs to more easily compile the Documenter Sass themes
into CSS files.

To compile an Sass file into a Documenter theme, you can use the [`Themes.compile`](@ref)
function:

```julia-repl
julia> using DocumenterTools: Themes

julia> Themes.compile("mytheme.scss")
```

When working on the Documenter built-in themes, the [`Themes.themewatcher`](@ref) function
can be used to automatically update all the built-in themes when any of the Sass files are
modified. To enable it, just run

```julia-repl
julia> using DocumenterTools: Themes

julia> Themes.themewatcher()
```

Note that it will read and overwrite the Sass and CSS files of the Documenter of the
environment DocumenterTools is loaded in — make sure that you have Documenter added as a
development dependency to that environment.
"""
module Themes
using FileWatching
using Documenter.Writers: HTMLWriter
using Sass

"""
    compile(src[, dst])

Compile an input Sass/SCSS file `src` into a CSS file. The standard Documenter Sass/SCSS
files are available in the include path.

The optional `dst` argument can be used to specify the output file. Otherwise, the file
extension of the `src` file is simply replaced by `.css`.
"""
function compile(src, dst=nothing)
    isfile(src) || error("SCSS source file missing: $src")
    if dst === nothing
        s = (endswith(src, ".scss") || endswith(src, ".sass")) ? first(splitext(src)) : src
        dst = "$(s).css"
    end
    # Create the directory for destination file if it doesn't exist
    isdir(dirname(dst)) || mkpath(dirname(dst))
    Sass.compile_file(src, dst; include_path=HTMLWriter.ASSETS_SASS)
end

"""
    compile_native_theme(name; dst=nothing)

Compiles a native Documenter theme and places it into Documenter's assets directory.

Optionally, the `dst` keyword argument can be used to specify the output file.
"""
function compile_native_theme(name; dst = nothing)
    name in HTMLWriter.THEMES || error("Bad theme name. Valid themes: $(HTMLWriter.THEMES)")
    src = joinpath(HTMLWriter.ASSETS_SASS, "$(name).scss")
    dst = (dst === nothing) ? joinpath(HTMLWriter.ASSETS_THEMES, "$(name).css") : dst
    compile(src, dst)
    return dst
end

function compile_native_themes()
    for theme in HTMLWriter.THEMES
        compile_native_theme(theme)
    end
end

# themewatcher() watches the Documenter theme source SCSS/Sass files and recompiles
# the native themes if there are any changes
const compiler_channel = Channel{Any}(1024)
const tasks = Vector{Tuple{Task,Channel{Any}}}()

"""
    themewatcher()

Starts an asynchronous background task that checks for changes in the Documenter Sass files
and recompiles all native themes whenever changes are detected.
"""
themewatcher(; kwargs...) = isempty(tasks) ?
    themewatcher_start(; kwargs...) : themewatcher_stop(; kwargs...)

function themewatcher_start(; dst=nothing, kwargs...)
    # Start the compiler task
    t = @async compiler(compiler_channel; dst=dst)
    push!(tasks, (t, compiler_channel))
    # Start watchers on all theme subdirectories
    for directory in all_subdirectories(HTMLWriter.ASSETS_SASS)
        start_watcher(directory)
    end
end

function themewatcher_stop(; kwargs...)
    # Send stop messages to all tasks
    for (_, c) in tasks
        put!(c, false)
    end
    # Wait for all of the tasks to finish
    for (t, _) in tasks
        try
            wait(t)
        catch e
            @warn "Task $t failed with exception" e
        end
    end
    # clear out the task list
    empty!(tasks)
    # clear out the compiler_channel channel
    while isready(compiler_channel)
        take!(compiler_channel)
    end
end

function all_subdirectories(d)
    vcat(d, mapreduce(vcat, walkdir(d)) do (root, dirs, _)
        [joinpath(root, d) for d in dirs]
    end)
end

function start_watcher(directory)
    channel = Channel{Any}(32)
    task = @async watcher(channel, directory, compiler_channel)
    push!(tasks, (task, channel))
end

function watcher(channel, directory, compiler_channel)
    @info "Starting watcher for: $(directory)"
    while true
        (f, reason) = watch_folder(directory, 2)
        if reason.changed
            @info "Change detected: $(f)" directory
            if endswith(f, ".scss") || endswith(f, ".sass")
                put!(compiler_channel, true)
            else
                @info "Ignoring $(f) -- not Sass"
            end
        end
        if isready(channel)
            @info "Received message, exiting watcher." directory take!(channel)
            return
        end
    end
end

function compiler(compiler_channel; dst=nothing)
    @info "Starting Sass compiler task." dst
    while true
        sleep(1)
        if isready(compiler_channel)
            @info "Received messages to compiler."
            while isready(compiler_channel)
                msg = take!(compiler_channel)
                if !msg
                    @info "Exiting compiler."
                    return
                end
            end
            for theme in HTMLWriter.THEMES
                @info "Compiling: $(theme)"
                _dst = nothing
                if dst !== nothing
                    _dst = joinpath(dst, "$(theme).css")
                end
                try
                    compile_native_theme(theme; dst=_dst)
                catch e
                    @warn "Compilation of $(theme) failed" e
                end
            end
        end
    end
end

end

ext: jl


text:
module chapters_sharing_data
using Dash, DashHtmlComponents, DashCoreComponents

export examples

include("../../utils.jl")
include("../../styles.jl")

include("../../reusable_components/Header.jl")
include("../../reusable_components/Example.jl")
include("../../reusable_components/Syntax.jl")

scoping_wrong = LoadExampleCode("./dash_docs/chapters/sharing_data/examples/scoping_wrong.jl")
scoping = LoadExampleCode("./dash_docs/chapters/sharing_data/examples/scoping.jl")

examples = [scoping_wrong, scoping]

app =  dash()
scoping_wrong.callback!(app)
scoping.callback!(app)

app.layout = html_div() do
    Header("Sharing Data Between Callbacks"),
    html_blockquote(dcc_markdown("This is the 5th chapter of the [Dash Tutorial](/).
    The [previous chapter](/interactive-visualization) covered how to use callbacks with
    the `dcc_graph` component. The [next and final chapter](/faqs) covers
    frequently asked questions and gotchas. Just getting started? Make sure to [install the necessary
    dependencies](/installation)")),
    dcc_markdown("
    One of the core Dash principles explained in the [Getting Started Guide on Callbacks] is that
    **Dash callbacks must never modify variables outside of their scope**. It is not safe to
    modify any `global` variables. This chapter explains why and provides some alternative
    patterns for sharing state between callbacks."),

    Header("Why Share State?"),

    dcc_markdown("
    In some apps, you may have multiple callbacks that depend on expensive data processing tasks like
    making SQL queries, running simulations, or downloading data.

    Rather than have each callback run the same expensive task, you can have one callback run the task and
    then share the results with the rest of the callbacks.

    This need has been somewhat ameliorated now that you can have [multiple outputs](/basic-callbacks) for one callback.
    This way, that expensive task can be done once and immediately used in all the outputs. But in some cases this
    still isn't ideal, for example if there are simple follow-on tasks that modify the results, like unit conversions.
    We shouldn't need to repeat a large database query just to change the results from Fahrenheit to Celsius!"),

    Header("Why Global Variables Will Break Your App"),

    dcc_markdown("
    Dash is designed to work in multi-user environments where multiple people may view the application at the
    same time and will have **independent sessions**.

    If your app uses `global` variables, then one's user's session could set the variable to one value which
    would affect the next user's session.

    Dash is also designed to be able to run with **multiple workers** so that callbacks can be executed in parallel.

    When Dash apps run across multiple workers, their memory is *not shared*. This means that if you modify
    a global variable in one callback, that modification will not be applied to the rest of the workers.

    Here's a sketch of an app with a callback that modifies data out of its scope. This type of pattern *will not
    work reliably* for the reasons outlined above.

    "),
    scoping_wrong.source_code,
    Example(scoping_wrong.layout),

    dcc_markdown("""
    To fix this example, simply re-assign the filtered dataset to a new variable
    inside the callback, or follow one of the strategies outlined in the next
    part of this guide.

    """),

    scoping.source_code,
    Example(scoping.layout),

    Header("Sharing Data Between Callbacks"),

    dcc_markdown("""
    In order to share data safely across multiple Julia processes, we need to store the data
    somewhere that is accessible to each of the processes.

    There are two main places to store this data:

    1. In the user's browser session.
    2. On the disk (e.g. on a file or on a new database).

    The following examples illustrate these approaches.

    """),
    Header("Example 1 - Storing Data in the Browser With A Hidden Div")
    dcc_markdown("""
    To save data a user's browser session:
    * Implemented by saving the data as part of Dash's front-end store through
    the methonds explained in https://community.plotly.com/t/sharing-a-dataframe-between-plots/6173.

    * Data has to be converted to a string like JSON for storage and transport.

    * Data that is cached in this way will *only be available in the user's current sessiion*.
        * If you open up a new browser, the app's callbacks will always compute the data. The
        data is only cached and transported between callbacks within the session.

        * As such, unlike with caching, this method doesn't increase the memory footprint
        of the app.

        * There could be a cost in network transport. If you're sharing 10MB of data between callbacks,
        then that data will be transported over the network between each callback.

        * If the network cost is too high, then compute the aggregations upfront and transport those. Your
        app likely won't be displaying 10MB of data, it will just be displaying a subset or an aggregation
        of it.

    This example outlines how you can perform an expensive data processing step in one callback, serialize the
    output as JSON, and provide it as an input to the other callbacks. This example uses standard Dash callbacks
    and stores the JSON-ifed data inside a hidden div in the app.

    """),
    dcc_markdown("""

    ```
    global_df = DataFrame(CSV.File("..."))

    app.layout = html_div() do
        dcc_graph(id="graph"),
        html_table(id="table"),
        dcc_dropdown(id="dropdown"),
        html_div(id="intermediate-value", style = (display = "None"))

    end

    callback!(
        app,
        Output("intermediate-value", "children"),
        Input("dropdown", "value"),
    ) do value

        dff = your_expensive_clean_or_compute_step(value)
        return JSON.json(dff)
    end

    callback!(
        app,
        Output("graph", "figure"),
        Input("intermediate-value", "children"),
    ) do dff

        figure = create_figure(dff)
        return figure
    end

    callback!(
        app,
        Output("table", "children"),
        Input("intermediate-value", "children"),
    ) do dff

        table = create_table(dff)
        return table
    end
    ```

    """),
    Header("Example 2 - Computing Aggregations Upfront"),
    dcc_markdown("""
    Sending the computed data over the network can be expensive if the data is large. In some cases, serializing
    this data and JSON can also be expensive.

    In many cases, your app will only display a subset or an aggregation of the computed or filtered data.
    In these cases, you could precompute your aggregations in your data processing callback and transport these
    aggregations to the remaining callbacks.

    Here's a simple example of how you might transport filtered or aggregated data to multiple callbacks.

    """),
    dcc_markdown("""

    ```
    callback!(
        app,
        Output("intermediate-value", "children"),
        Input("dropdown", "value"),
    ) do value

        df = your_expensive_clearn_or_compute_step(value)

        # a few filter steps that compute the data
        # as it's needed in the future callback

        df1 = filter(row -> row.Fruit == "apples")
        df2 = filter(row -> row.Fruit == "oranges")
        df3 = filter(row -> row.Fruit == "figs")

        datasets = Dict("df1"=>JSON.json(df1), "df2"=>JSON.json(df2), "df3"=>JSON.json(df3))
        return JSON.json(datasets)
    end

    callback!(
        app,
        Output("graph", figure"),
        Input("intermediate-value", "children"),
    ) do dff

        figure = create_figure_1(dff)
        return figure

    end

    callback!(
        app,
        Output("graph", figure"),
        Input("intermediate-value", "children"),
    ) do dff

        figure = create_figure_2(dff)
        return figure

    end

    callback!(
        app,
        Output("graph", figure"),
        Input("intermediate-value", "children"),
    ) do dff

        figure = create_figure_3(dff)
        return figure

    end
    ```
    """),
    dcc_markdown("
    The next chapter of the tutorial covers frequently asked questions about Dash and
    some common gotchas.

    [Dash Tutorial Part 6: FAQs and Gotchas](/faqs)
    "),
    html_hr(),
    dcc_markdown("[Back to the table of contents](/)")

end

end

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
# SKCE

function CalibrationErrors.unsafe_skce_eval_targets(
    kernel::SqExponentialKernel{Euclidean},
    p::MvNormal,
    y::AbstractVector{<:Real},
    p̃::MvNormal,
    ỹ::AbstractVector{<:Real},
)
    # extract the parameters
    μ = p.μ
    Σ = p.Σ
    μ̃ = p̃.μ
    Σ̃ = p̃.Σ

    # compute inverse scaling matrices
    invA = LinearAlgebra.I + Σ
    invB = LinearAlgebra.I + Σ̃
    invC = invA + Σ̃

    res =
        kernel(y, ỹ) - mvnormal_gaussian_kernel(invA, μ, ỹ) -
        mvnormal_gaussian_kernel(invB, μ̃, y) + mvnormal_gaussian_kernel(invC, μ, μ̃)

    return res
end

function mvnormal_gaussian_kernel(invA, y, ỹ)
    return exp(-(LinearAlgebra.logdet(invA) + invquad_diff(invA, y, ỹ)) / 2)
end

# UCME

# predicted normal distributions with squared exponential kernel for the targets
function CalibrationErrors.unsafe_ucme_eval_targets(
    kernel::SqExponentialKernel{Euclidean},
    p::MvNormal,
    y::AbstractVector{<:Real},
    testp::MvNormal,
    testy::AbstractVector{<:Real},
)
    # compute inverse scaling matrix
    invA = LinearAlgebra.I + p.Σ

    return kernel(y, testy) - mvnormal_gaussian_kernel(invA, p.μ, testy)
end

# kernels with input transformations
# TODO: scale upfront?

function CalibrationErrors.unsafe_skce_eval_targets(
    kernel::TransformedKernel{
        SqExponentialKernel{Euclidean},<:Union{ScaleTransform,ARDTransform,LinearTransform}
    },
    p::MvNormal,
    y::AbstractVector{<:Real},
    p̃::MvNormal,
    ỹ::AbstractVector{<:Real},
)
    # obtain the transform
    t = kernel.transform

    return CalibrationErrors.unsafe_skce_eval_targets(
        SqExponentialKernel(), apply(t, p), t(y), apply(t, p̃), t(ỹ)
    )
end

function CalibrationErrors.unsafe_ucme_eval_targets(
    kernel::TransformedKernel{
        SqExponentialKernel{Euclidean},<:Union{ScaleTransform,ARDTransform,LinearTransform}
    },
    p::MvNormal,
    y::AbstractVector{<:Real},
    testp::MvNormal,
    testy::AbstractVector{<:Real},
)
    # obtain the transform
    t = kernel.transform

    # `testp` is irrelevant for the evaluation and therefore not transformed
    return CalibrationErrors.unsafe_ucme_eval_targets(
        SqExponentialKernel(), apply(t, p), t(y), testp, t(testy)
    )
end

## utilities

# internal `apply` avoids type piracy of transforms
function apply(t::ScaleTransform, d::MvNormal)
    s = first(t.s)
    return MvNormal(t(d.μ), s^2 * d.Σ)
end
function apply(t::ARDTransform, d::MvNormal)
    return MvNormal(t(d.μ), scale_cov(d.Σ, t.v))
end
# `X_A_Xt` only works with StridedMatrix: https://github.com/JuliaStats/PDMats.jl/issues/96
apply(t::LinearTransform, d::MvNormal) = Matrix(t.A) * d

# scale covariance matrix
# TODO: improve efficiency
scale_cov(Σ::PDMats.ScalMat, v) = PDMats.PDiagMat(v .^ 2 .* Σ.value)
scale_cov(Σ::PDMats.PDiagMat, v) = PDMats.PDiagMat(v .^ 2 .* Σ.diag)
scale_cov(Σ::PDMats.AbstractPDMat, v) = PDMats.X_A_Xt(Σ, LinearAlgebra.diagm(v))

invquad_diff(A::PDMats.ScalMat, x, y) = Distances.sqeuclidean(x, y) / A.value
function invquad_diff(
    A::PDMats.PDiagMat, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}
)
    n = length(x)
    length(y) == n || throw(DimensionMismatch("x and y must be of the same length"))
    size(A) == (n, n) ||
        throw(DimensionMismatch("size of A is not consistent with x and y"))
    return sum(abs2(xi - yi) / wi for (xi, yi, wi) in zip(x, y, A.diag))
end
invquad_diff(A::PDMats.AbstractPDMat, x, y) = PDMats.invquad(A, x .- y)


text:
using RobertoMD

using Rotations
using JSON

control=Dict("dt"=>0.01,
            "steps"=>1000,
            "velocity verlet"=>true,
            "LAMMPSTrj"=>Dict("file"=>"copolymer.lammpstrj","freq"=>1000),
            "restart"=>Dict("JSONRestart"=>true,"file"=>"copolymer.json","freq"=>1000),
            "BerendsenNVT"=>Dict("tau"=>1.0),
            "harmonic bond"=>Dict(
            "k"=>Dict("1"=>1000.0),
            "l0"=>Dict("1"=>1.0)
            ),
            "Canonical field"=>Dict(
            "χ"=>Dict(  "1"=>[0.0,5.0],
                        "2"=>[5.0,0.0] ),
            "κ"=>0.2,
            "uniform mesh"=>true,
            "update"=>1,
            "Lcell"=>1.0),
            "thermo information"=>Dict( "freq"=>1,
                                        "energy"=>true,
                                       "momentum"=>true,
                                       "write"=>true,
                                        "file"=>"copolymer.log"),
            "density"=>0.85,
            "bondlength"=>1.0,
            "chain_length"=>20,
            "num_polymers"=>500,
            "zero velocity"=>true,
            )

boxsize=(control["num_polymers"]*control["chain_length"]/control["density"])^(1.0/3)
bondlength=control["bondlength"]
chain_length=control["chain_length"]

polymer=Dict()
for polymer_i in 1:control["num_polymers"]
    polymer[string(polymer_i)]=Dict()
    polymer[string(polymer_i)]["atoms"]=Dict()
    polymer[string(polymer_i)]["bonds"]=Dict()
    pos_init=rand(3)*boxsize
    for monomer_i in 1:chain_length
        r = rand(RotMatrix{3})
        q = QuatRotation(r)
        pos_init=pos_init+q*[bondlength,0.0,0.0]
        type_=1
        if monomer_i>div(control["chain_length"],2)
            type_=2
        end
        polymer[string(polymer_i)]["atoms"][string(monomer_i)]=Dict("type"=>type_,"mass"=>1,"coords"=>pos_init)
    end
    for bond_i in 1:chain_length-1
        polymer[string(polymer_i)]["bonds"][string(bond_i)]=[1,bond_i,bond_i+1]
    end

end

configs=Dict(
    "box"=>[boxsize,boxsize,boxsize],
    "molecules"=>polymer
    )

#control_file=open("control.json","w")
#JSON.print(control_file,control,4)

#=
=#
#configuration_file=open("configuration.json","w")
#JSON.print(configuration_file,configs,4)
#=
configs=JSON.parse(open("copolymer.json","r"))
=#

Simulate(control,configs)


text:
module GigaSOM
export hello

"""
    hello(who::String)
Return "Hello, `who`".
"""
hello(who::String) = "Hello, $who"

end # module


text:
# --------------------------------------------------
# WCS functions

function wcs2d_from_header(header::FITSHeader)
    cdelt = [header["CDELT" * string(i)] for i in [1, 2]]
    ctype = [header["CTYPE" * string(i)] for i in [1, 2]]
    crpix = convert(Array{Float64, 1}, [header["CRPIX" * string(i)] for i in [1, 2]])
    crval = [header["CRVAL" * string(i)] for i in [1, 2]];
    wcs = WCSTransform(2;
        cdelt = cdelt,
        ctype = ctype,
        crpix = crpix,
        crval = crval
    );
    return wcs
end

function wcs3d_from_header(header::FITSHeader)
    cdelt = [header["CDELT" * string(i)] for i in [1, 2, 3]]
    ctype = [header["CTYPE" * string(i)] for i in [1, 2, 3]]
    crpix = convert(Array{Float64, 1}, [header["CRPIX" * string(i)] for i in [1, 2, 3]])
    crval = [header["CRVAL" * string(i)] for i in [1, 2, 3]];
    wcs = WCSTransform(3;
        cdelt = cdelt,
        ctype = ctype,
        crpix = crpix,
        crval = crval
    );
    return wcs
end

function parse_fitsheader(header::String)
    n = length(header)
    lines = [header[(i - 1)*80 + 1: i * 80] for i in 1:(length(header) ÷ 80)]
    keys = []
    values = []
    comments = []
    for i in lines
        if (startswith(i, "COMMENT")) || length(strip(i)) == 0
            continue
        end
        k, v = split(i, "=")
        push!(keys, string(strip(k)))
        val, comment = split(v, "/")
        val = string(strip(val))
        val_types = [tryparse(T, val) for T in [Int, Float64]]
        if val_types[1] != nothing
            val = val_types[1]
        elseif val_types[2] != nothing
            val = val_types[2]
        else
            val = replace(val, "'" => "")
        end

        push!(values, val)
        push!(comments, string(strip(comment)))
    end
    return convert(Array{String, 1}, keys), values, convert(Array{String, 1}, comments)
end

function header2d_of(cube::AbstractSpectralCube, line::AbstractString, content_description::AbstractString)
    hkeys, hvalues, comments = wcs2d_from_header(cube.header) |> WCS.to_header |> parse_fitsheader
    hkeys = ["SIMPLE" ; hkeys; "DATA"; "LINE"]
    hvalues = [true; hvalues; content_description; line];
    comments = [cube.header.comments[1]; comments; ""; ""]
    header = FITSHeader(hkeys, hvalues, comments);
    return header
end


text:
using Test
using KoreanJosa

@testset "Do nothing" begin
    @test josa("Korean josa system is hard to learn") == "Korean josa system is hard to learn"
    @test josa("1986 + 2021") == "1986 + 2021"
end

@testset "(이), (가)" begin
    @test josa("이것(이)") == josa("이것(가)") == "이것이"
    @test josa("구글(이)") == josa("구글(가)") == "구글이"
    @test josa("니(이)") == josa("니(가)") == "니가"
    @test josa("누(이)") == josa("누(가)") == "누가"
end

@testset "(과), (와)" begin
    @test josa("우리(과)") == josa("우리(와)") == "우리와"
    @test josa("아이폰(과)") == josa("아이폰(와)") == "아이폰과"
    @test josa("갤럭시(과)") == josa("갤럭시(와)") == "갤럭시와"
    @test josa("학교(과)") == josa("학교(와)") == "학교와"
end

@testset "(을), (를)" begin
    @test josa("이것(을)") == josa("이것(를)") == "이것을"
    @test josa("장난(을)") == josa("장난(를)") == "장난을"
    @test josa("고기(을)") == josa("고기(를)") == "고기를"
    @test josa("장소(을)") == josa("장소(를)") == "장소를"
end

@testset "(은), (는)" begin
    @test josa("왕(은)") == josa("왕(는)") == "왕은"
    @test josa("목숨(은)") == josa("목숨(는)") == "목숨은"
    @test josa("신하(은)") == josa("신하(는)") == "신하는"
    @test josa("자네(은)") == josa("자네(는)") == "자네는"
end

@testset "(으로), (로)" begin
    @test josa("바다(으로)") == josa("바다(로)") == "바다로"
    @test josa("서울(으로)") == josa("서울(로)") == "서울로"
    @test josa("부산(으로)") == josa("부산(로)") == "부산으로"
    @test josa("대학(으로)") == josa("대학(로)") == "대학으로"
end

@testset "(아), (야)" begin
    @test josa("수진(아)") == josa("수진(야)") == "수진아"
    @test josa("재석(아)") == josa("재석(야)") == "재석아"
    @test josa("혜미(아)") == josa("혜미(야)") == "혜미야"
    @test josa("용희(아)") == josa("용희(야)") == "용희야"
end

@testset "(이여), (여)" begin
    @test josa("장사꾼(이여)") == josa("장사꾼(여)") == "장사꾼이여"
    @test josa("학생(이여)") == josa("학생(여)") == "학생이여"
    @test josa("그대(이여)") == josa("그대(여)") == "그대여"
    @test josa("용사(이여)") == josa("용사(여)") == "용사여"
end

@testset "(이라), (라)" begin
    @test josa("오늘(이라)") == josa("오늘(라)") == "오늘이라"
    @test josa("고놈(이라)") == josa("고놈(라)") == "고놈이라"
    @test josa("모레(이라)") == josa("모레(라)") == "모레라"
    @test josa("제주도(이라)") == josa("제주도(라)") == "제주도라"
end

@testset "All" begin
    s = "오늘(는) 맑고푸르른 하늘(와) 청명한 바다(을) 보며, 한국(로) 돌아가고 싶(는) 마음에 말(를) 꺼냈다. 혜미(아), 정철(야) 너희(은) 어떻게 생각하니? 글쎄... 선생님(라)면 오랜만(여)서 반가워 하실 것 같아."
    @test josa(s) == "오늘은 맑고푸르른 하늘과 청명한 바다를 보며, 한국으로 돌아가고 싶은 마음에 말을 꺼냈다. 혜미야, 정철아 너희는 어떻게 생각하니? 글쎄... 선생님이라면 오랜만이여서 반가워 하실 것 같아."
end

@testset "Exception" begin
    @test_throws ArgumentError KoreanJosa.JOSAMarker("(읍)")
    @test_throws ArgumentError KoreanJosa.JOSA_이_가[3]
end


text:
# This file was generated, do not modify it. # hide
#hideall
using PyCall
lines = replace("""import numpy as np
np.random.seed(2)
x = np.random.randn(5)
r = np.linalg.norm(x) / len(x)
np.round(r, 2)""", r"(^|\n)([^\n]+)\n?$" => s"\1res = \2")
py"""
$$lines
"""
println(py"res")

max_stars_repo_name: mcognetta/franklindocs

max_issues_repo_name: mcognetta/franklindocs

max_forks_repo_name: mcognetta/franklindocs


text:
using MAT
import Nettle

images_url = "http://www.rctn.org/bruno/sparsenet/IMAGES.mat"
images_file = "$(@__DIR__)/IMAGES.mat"
images_md5 = "b20cb33935e95b556cb83bef5e9afcf1"

file_ok = false

if isfile(images_file)
    h = Nettle.Hasher("md5")
    Nettle.update!(h, read(images_file))
    if Nettle.hexdigest!(h) == images_md5
        file_ok = true
    end
end

if !file_ok
    println("Downloading $images_url")
    download(images_url, images_file)
end

function nat_images()
    f = matopen(images_file)
    read(f, "IMAGES")
end


text:
# # Fitting SNAP using GaN data
#
# This example:
#    - Uses the atomic positions in the folder `example_GaN_data`
#    - Generates surrogate DFT data based on the GaN model presented in `10.1088/1361-648x/ab6cbe`
#    - Uses `snap.jl` and 80% of the GaN data to create the matrix A. The matrix is generated only with the energy block.
#    - Uses 80% of the GaN data to create the vector b. The reference energy (``E_{\rm ref}``) is assumed to be zero.
#    - Uses backslash to fit the parameters ``\boldsymbol{\beta}``, thus, solves ``\mathbf{A \cdot \boldsymbol{\beta}=y}``
#    - Uses 20% of the GaN data and the bispectrum components provided by `snap.jl`and the fitted parameters ``\boldsymbol{\beta}`` to validate the fitting.
#    - The error computed during the validation is < ~5%.
# TODO:
#    - Check if the GaN model is correct.
#    - Calculate reference energy (``E_{\rm ref}``).
#    - Check ill-conditioned system.

#
# ## Install and import dependencies
#
# First let's make sure we have all required packages installed

# ```julia
# using Pkg
# pkg"add LAMMPS"
# ```
#
# We start off by importing the necessary packages

using LAMMPS
using LinearAlgebra: norm, pinv
using Printf

# ##  Estimation of the SNAP coefficients
#
# The choice of the coefficients ``\boldsymbol{\beta}=(\beta_0^1, \tilde{\beta}^1, \dots, \beta_0^l, \tilde{\beta}^l)``
# is based on a system of linear equations which considers a large number of atomic
# configurations and ``l`` atom types. The matrix formulation for this system ``\mathbf{A \cdot \boldsymbol{\beta}=y}``
# is defined in the following equations (see 10.1016/j.jcp.2014.12.018):

# ```math
# \begin{equation*}
#    \mathbf{A}=
#    \begin{pmatrix}
#        \vdots &  &  & & & & & & \\
#        N_{s_1} & \sum_{i=1}^{N_{s_1}} B_1^i & \dots & \sum_{i=1}^{N_{s_1}} B_k^i  & \dots & N_{s_L} & \sum_{i=1}^{N_{s_L}} B_1^i & \dots & \sum_{i=1}^{N_{s_L}} B_k^i\\
#        \vdots &  &  & & & & & & \\
#        0 & -\sum_{i=1}^{N_{s_1}} \frac{\partial B_1^i}{\partial r_j^{\alpha}} & \dots & -\sum_{i=1}^{N_{s_1}} \frac{\partial B_k^i}{\partial r_j^{\alpha}} & \dots &  0 & -\sum_{i=1}^{N_{s_l}} \frac{\partial B_1^i}{\partial r_j^{\alpha}} & \dots & -\sum_{i=1}^{N_{s_l}} \frac{\partial B_k^i}{\partial r_j^{\alpha}} \\
#        \vdots &  &  & & & & & & \\
#        0 & - \sum_{j=1}^{N_{s_1}} r^j_{\alpha} \sum_{i=1}^{N_{s_1}} \frac{\partial B_1^i}{\partial r_j^{\beta}} & \dots & - \sum_{j=1}^{N_{s_1}} r^j_{\alpha} \sum_{i=1}^{N_{s_1}} \frac{\partial B_k^i}{\partial r_j^{\beta}} & \dots & 0 & - \sum_{j=1}^{N_{s_l}} r^j_{\alpha} \sum_{i=1}^{N_{s_l}} \frac{\partial B_1^i}{\partial r_j^{\beta}} & \dots & - \sum_{j=1}^{N_{s_l}} r^j_{\alpha} \sum_{i=1}^{N_{s_l}} \frac{\partial B_k^i}{\partial r_j^{\beta}}\\
#        \vdots &  &  & & & & & & \\
#    \end{pmatrix}
# \end{equation*}
# ```

# The indexes ``\alpha, \beta = 1,2,3`` depict the ``x``, ``y`` and ``z``
# spatial component, ``j`` is an individual atom, and ``s`` a particular configuration.
# All atoms in each configuration are considered. The number of atoms of type $m$ in the configuration
# ``s`` is ``N_{s_m}``.

# The RHS of the linear system is computed as:

# ```math
# \begin{equation*}
#  \mathbf{y}=  \begin{pmatrix}
#    \vdots \\
#   E^s_{\rm qm} -  E^s_{\rm ref} \\
#   \vdots \\\\
#   F^{s,j,\alpha}_{\rm qm} - F^{s,j,\alpha}_{\rm ref} \\
#   \vdots \\
#   W_{\rm qm}^{s,\alpha,\beta} - W_{\rm ref}^{s,\alpha,\beta} \\
#   \vdots \\
#    \end{pmatrix}
# \end{equation*}
# ```

# # Calculate ``A`` using `snap.jl`

# This example only fits the energy, thus, the first block of the matrix A.

include(joinpath(dirname(pathof(LAMMPS)), "..", "examples", "snap.jl"))
A

# # Calculate ``y``

# A molecular mechanics model for the interaction of gallium and nitride is used
# to generate surrogate DFT data (see 10.1088/1361-648x/ab6cbe).
# The total potential energy is computed for each configuration

# ```math
# E = \sum_{i \lt j;  r_{i,j} \leq rcut } E_{GaN}(r_{i,j}) \  \text{ where } \\
#
# E_{GaN}(r_{i,j}) =
# \left\{
#    \begin{array}{ll}
#        E_{C}(r_{i,j}) + E_{LJ}(r_{i,j}, \epsilon_{Ga,Ga}, \sigma_{Ga,Ga})  & \mbox{if both atoms are Ga} \\
#        E_{C}(r_{i,j}) + E_{LJ}(r_{i,j}, \epsilon_{N,N}, \sigma_{N,N})  & \mbox{if both atoms are N} \\
#        E_{C}(r_{i,j}) + E_{BM}(r_{i,j}, A_{Ga,N}, \rho_{Ga,N})  & \mbox{if one atom is Ga and the other is N}\\
#    \end{array}
# \right.
# ```

# Note: this is a work in progress, this model should be checked.

const N = N1 + N2
const M1 = M
const M2 = 61

const ε_Ga_Ga = 0.643
const σ_Ga_Ga = 2.390
const ε_N_N = 1.474
const σ_N_N = 1.981
const A_Ga_N = 608.54
const ρ_Ga_N = 0.435
const q_Ga = 3.0
const q_N = -3.0
const ε0 = 55.26349406 # e2⋅GeV−1⋅fm−1 ToDo: check this
const E_ref = zeros(M1)
E_LJ(r, ε = 1.0, σ = 1.0) = 4.0 * ε * ((σ / norm(r))^12 - (σ / norm(r))^6)
E_BM(r, A = 1.0, ρ = 1.0) = A * exp(-norm(r) / ρ)
E_C(r) = q_Ga * q_N / (4.0 * π * ε0 * norm(r))

function read_atomic_conf(m, N)
    rs = []
    open(joinpath(DATA, string(m), "DATA")) do f
        for i = 1:23
            readline(f)
        end
        for i = 1:N
            s = split(readline(f))
            r = [parse(Float64, s[3]),
                 parse(Float64, s[4]),
                 parse(Float64, s[5])]
            push!(rs, r)
        end
    end
    return rs
end

function calc_tot_energy(rcut, rs, N, ε_Ga_Ga, σ_Ga_Ga, ε_N_N, σ_N_N, A_Ga_N, ρ_Ga_N)
    E_tot_acc = 0.0
    for i = 1:N
        for j = i:N
            r_diff = rs[i] - rs[j]
            if norm(r_diff) <= rcut && norm(r_diff) > 0.0
                if i <= N1 && j <= N1
                    E_tot_acc += E_C(r_diff) + E_LJ(r_diff, ε_Ga_Ga, σ_Ga_Ga)
                elseif i > N1 && j > N1
                    E_tot_acc += E_C(r_diff) + E_LJ(r_diff, ε_N_N, σ_N_N)
                else
                    E_tot_acc += E_C(r_diff) + E_BM(r_diff, A_Ga_N, ρ_Ga_N)
                end
            end
        end
    end
    return E_tot_acc
end

# Vector ``y`` is calculated using the DFT data and ``E_{\rm ref}``

function calc_y(rcut, M1, N, ε_Ga_Ga, σ_Ga_Ga, ε_N_N, σ_N_N, A_Ga_N, ρ_Ga_N, E_ref)
    y = zeros(M1)
    for m = 1:M1
        rs = read_atomic_conf(m, N)
        y[m] = calc_tot_energy(rcut, rs, N, ε_Ga_Ga, σ_Ga_Ga,
                               ε_N_N, σ_N_N, A_Ga_N, ρ_Ga_N) - E_ref[m]
    end
    return y
end

y = calc_y(rcut, M1, N, ε_Ga_Ga, σ_Ga_Ga, ε_N_N, σ_N_N, A_Ga_N, ρ_Ga_N, E_ref)

# ## Calculate the optimal solution

# The optimal solution $\widehat{\mathbf{\beta}}$ is, thus:
# ```math
# \begin{equation*}
#     \widehat{\mathbf{\beta}} = \mathrm{argmin}_{\mathbf{\beta}} ||\mathbf{A \cdot \boldsymbol{\beta} - y}||^2 = \mathbf{A^{-1} \cdot y}
# \end{equation*}
# ```

β = A \ y

# ## Check results

# The local energy can be decomposed into separate contributions for each atom.
# SNAP energy can be written in terms of the bispectrum components of the atoms
# and a set of coefficients. ``K`` components of the bispectrum are considered
# so that ``\mathbf{B}^{i}=\{ B^i_1, \dots, B_K^i\}`` for each atom ``i``, whose
# SNAP energy is computed as follows:
# ```math
#    E^i_{\rm SNAP}(\mathbf{B}^i) = \beta_0^{\alpha_i} + \sum_{k=1}^K \beta_k^{\alpha_i} B_k^i =  \beta_0^{\alpha_i} + \boldsymbol{\tilde{\beta}}^{\alpha_i} \cdot \mathbf{B}^i
# ```
# where $\alpha_i$ depends on the atom type.

function calc_fitted_tot_energy(path, β, ncoeff, N1, N)
    ## Calculate y
    lmp = LMP(["-screen","none"])
    bs = run_snap(lmp, path, rcut, twojmax)

    E_tot_acc = 0.0
    for n in 1:N1
        E_atom_acc = β[1]
        for k in 2:ncoeff
            k2 = k - 1
            E_atom_acc += β[k] * bs[k2, n]
        end
        E_tot_acc += E_atom_acc
    end
    for n in N1+1:N
        E_atom_acc = β[ncoeff+1]
        for k in ncoeff+2:2*ncoeff
            k2 = k - ncoeff - 1
            E_atom_acc += β[k] * bs[k2, n]
        end
        E_tot_acc += E_atom_acc
    end
    command(lmp, "clear")
    return E_tot_acc
end

# Comparison between the potential energy calculated based on the DFT data, and
# the SNAP potential energy calculated based on the bispectrum components
# provided by ``snap.jl`` and the fitted coefficients ``\mathbf{\beta}``.

@printf("Potential Energy, Fitted Potential Energy, Error (%%)\n")
for m = M1+1:M2
    path = joinpath(DATA, string(m), "DATA")
    rs = read_atomic_conf(m, N)
    E_tot = calc_tot_energy(rcut, rs, N, ε_Ga_Ga, σ_Ga_Ga, ε_N_N, σ_N_N, A_Ga_N, ρ_Ga_N)
    E_tot_fit = calc_fitted_tot_energy(path, β, ncoeff, N1, N)
    @printf("%0.2f, %0.2f, %0.2f\n", E_tot, E_tot_fit,
            abs(E_tot - E_tot_fit) / E_tot * 100.)
end


text:
struct PosixPath <: AbstractPath
    parts::Tuple{Vararg{String}}
end

PosixPath() = PosixPath(tuple())

function PosixPath(str::AbstractString)
    str = string(str)

    if isempty(str)
        return PosixPath(tuple("."))
    end

    tokenized = split(str, POSIX_PATH_SEPARATOR)
    if isempty(tokenized[1])
        tokenized[1] = POSIX_PATH_SEPARATOR
    end
    return PosixPath(tuple(map(String, tokenized)...))
end

# The following should be implemented in the concrete types
==(a::PosixPath, b::PosixPath) = parts(a) == parts(b)
parts(path::PosixPath) = path.parts
ispathtype(::Type{PosixPath}, str::AbstractString) = Sys.isunix()

function isabs(path::PosixPath)
    if parts(path)[1] == POSIX_PATH_SEPARATOR
        return true
    else
        return false
    end
end

drive(path::PosixPath) = ""

function root(path::PosixPath)
    if parts(path)[1] == POSIX_PATH_SEPARATOR
        return POSIX_PATH_SEPARATOR
    else
        return ""
    end
end

function expanduser(path::PosixPath)
    p = parts(path)

    if p[1] == "~"
        if length(p) > 1
            return PosixPath(tuple(homedir(), p[2:end]...))
        else
            return PosixPath(tuple(homedir()))
        end
    end

    return path
end


text:
"""
An upper arrowhead matrix:

A = [d_1           e_1]
    [    d_2       e_2]
    [       d_3    e_3]
    [           ⋱   ⋮ ]
    [               f ]
"""
struct UpperArrow{T} <: AbstractMatrix{T}
    d::Vector{T}
    e::Vector{T}
    f::T
end

size(A::UpperArrow) = length(A.d) + 1, length(A.d) + 1

function getindex(A::UpperArrow{T}, i::Integer, j::Integer) where T
    n = size(A, 1)
    if i == j
        if i == n
            return A.f
        else
            return A.d[i]
        end
    elseif j == n
        return A.e[i]
    else
        return zero(T)
    end
end

issymmetric(A::UpperArrow) = false

+(A::UpperArrow, B::UpperArrow) = UpperArrow(A.d+B.d, A.e+B.e, A.f+B.f)
-(A::UpperArrow, B::UpperArrow) = UpperArrow(A.d-B.d, A.e-B.e, A.f-B.f)

function chol(A::SymArrow{T}) where T
    a, b, c = A.a, A.b, A.c
    d = sqrt.(a)
    e = b./d
    f = sqrt(c - norm(e)^2)
    UpperArrow(d, e, f)
end

function inv(A::UpperArrow{T}) where T
    d, e, f = A.d, A.e, A.f
    d1 = inv.(d)
    f1 = inv(f)
    e1 = scale!(-f1, d1.*e)
    UpperArrow(d1, e1, f1)
end

function similarity(A::SymArrow{T}, B::SymArrow{T}) where T
    n = size(A, 1)
    R = chol(B)
    Ri = inv(R)
    a, b, c = A.a, A.b, A.c
    d, e, f = Ri.d, Ri.e, Ri.f
    a1 = Vector{T}(n - 1)
    b1 = Vector{T}(n - 1)
    c1 = c*f^2
    twof = 2f
    @inbounds @simd for i = 1:n-1
        a1[i] = a[i]*d[i]^2
        ae = a[i]*e[i]
        b1[i] = d[i]*(ae + b[i]*f)
        c1 += e[i]*(ae + b[i]*twof)
    end
    SymArrow(a1, b1, c1)
end

eigvals(A::SymArrow{T}, B::SymArrow{T}) where T = eigvals(similarity(A, B))


text:
using FileIO
using FixedPointNumbers
using ColorTypes
function work(name::String;range::Union{Nothing,Vector}=nothing,dir::String=homedir()*"/Desktop",wid::Union{Nothing,Int}=nothing,type::String="jpg",result::String=name*"0."*type)
	num=0
	opened=Vector{Matrix{RGB{N0f8}}}()
	if range===nothing
		v=readdir(dir;sort=false)
		cd(dir)
		for i in 1:32
			fname="$name$i.$type"
			if in(fname,v)
				push!(opened,load(fname))
			else
				num=i-1
				break
			end
		end
	else
		num=length(range::Vector)
		cd(dir)
		for i in range
			push!(opened,load("$name$i.$type"))
		end
	end
	if wid===nothing
		wid=num&1==1 ? 1 : 2
	end
	k1=div(num,wid)
	k2=num%wid
	grw=0
	grh=0
	maxgrw=0
	grhs=Vector{Int}()
	count=1
	for img in opened
		tup=size(img)
		h=tup[1]
		w=tup[2]
		grw+=w
		grh=max(h,grh)
		if count%wid==0
			maxgrw=max(maxgrw,grw)
			push!(grhs,grh)
			grw=0
			grh=0
		end
		count+=1
	end
	if k2!=0
		maxgrw=max(maxgrw,grw)
		push!(grhs,grh)
	end
	res=fill(RGB(0,0,0),sum(grhs),maxgrw)
	curline=0
	curcol=0
	for i in 1:k1
		for j in 1:wid
			img=opened[(i-1)*wid+j]
			tup=size(img)
			for x in 1:tup[1]
				for y in 1:tup[2]
					res[x+curline,y+curcol]=img[x,y]
				end
			end
			curcol+=tup[2]
		end
		curline+=grhs[i]
		curcol=0
	end
	for j in 1:k2
		img=opened[k1*wid+j]
		tup=size(img)
		for x in 1:tup[1]
			for y in 1:tup[2]
				res[x+curline,y+curcol]=img[x,y]
			end
		end
		curcol+=tup[2]
	end
	save(result,res)
end

max_stars_repo_name: Rratic/misc

max_issues_repo_name: Rratic/misc

max_forks_repo_name: Rratic/misc


text:
using Gadfly, DataArrays, RDatasets

set_default_plot_size(6inch, 3inch)

plot(dataset("car", "SLID"), x="Wages", color="Language", Stat.histogram, Geom.line)


text:
__precompile__()

module SRKGenerator

using NLopt, CUDArt

const END_C = true
const INEQ_C = false
const EQ_C = true

include("getCoefs.jl")
include("main.jl")
include("optim_functions.jl")
include("printing.jl")
include("translate.jl")

export srk_optimize

end # module


text:
using Test
import SWI: Options, saveconfiguration

σ = [2,2,1]
unwrapping = :laplacian
mode = :linear
level = 4
tmp_folder = "tmp"
o = Options(;σ=σ, unwrapping=unwrapping, mode=mode, level=level, writedir=tmp_folder)
@test o.σ == σ
@test o.unwrapping == unwrapping
@test o.mode == mode
@test o.level == level
@test o.combination_type == :SNR
@test o.sensitivity == nothing
@test o.writedir == "tmp"
@test o.writesteps == false

mkpath(tmp_folder)
saveconfiguration(o)
text = """unwrapping: laplacian
mode: linear
level: 4
combination_type: SNR
sensitivity: nothing
writedir: tmp
writesteps: false
"""
file = "$tmp_folder/settings_swi.txt"
@test read(file, String) == text
rm(file)
try rm(tmp_folder) catch end

printstyled("Tests Passed!"; color=:green)


text:
@testset "Percolation" begin

    @test isapprox(percolation(gridnetwork(5, 5)), 0.574468, rtol = 1e-5)
    @test isapprox(percolation(gridnetwork(5, 6)), 0.584746, rtol = 1e-5)
    @test isapprox(percolation(gridnetwork(6, 6)), 0.594595, rtol = 1e-5)
end


text:
"""
Abstract type for time series storage implementations.

All subtypes must implement:
- add_time_series!
- add_time_series_reference!
- remove_time_series!
- get_time_series
- clear_time_series!
- get_num_time_series
- check_read_only
"""
abstract type TimeSeriesStorage end

function make_time_series_storage(;
    in_memory = false,
    filename = nothing,
    directory = nothing,
)
    if in_memory
        storage = InMemoryTimeSeriesStorage()
    elseif !isnothing(filename)
        storage = Hdf5TimeSeriesStorage(; filename = filename)
    else
        storage = Hdf5TimeSeriesStorage(true; directory = directory)
    end

    return storage
end

function make_component_label(component_uuid::UUIDs.UUID, label::AbstractString)
    return string(component_uuid) * "_" * label
end

function deserialize_component_label(component_label::AbstractString)
    data = split(component_label, "_")
    component = UUIDs.UUID(data[1])
    label = data[2]
    return component, label
end

function serialize(storage::TimeSeriesStorage, file_path::AbstractString)
    if storage isa Hdf5TimeSeriesStorage
        # The data is currently in a temp file, so we can just make a copy.
        cp(get_file_path(storage), file_path; force = true)
    elseif storage isa InMemoryTimeSeriesStorage
        convert_to_hdf5(storage, file_path)
    else
        error("unsupported type $(typeof(storage))")
    end

    @info "Serialized time series data to $file_path."
end


text:
module CpuDispatcher

using CpuId, Hwloc

# cat /proc/info
proccpuinfo() = read(pipeline(`cat /proc/cpuinfo`, `grep flags`, `uniq`, `cut -c 10-`, `tr ' ' '\n'`, `tr a-z A-Z`, `sort`),String) |> chomp |> x -> split(x, '\n')

fs_proc = proccpuinfo()

# jl_dump_host_cpu

so = IOBuffer()
x = run(pipeline(`julia -e "ccall(:jl_dump_host_cpu,Cvoid,())"`, stdout=so, stderr=so))
String(take!(so))

function jldumphost()
	read(`julia -e "ccall(:jl_dump_host_cpu,Cvoid,())"`, String) # FIX PLZ IZ BROKEN
	jldump = "sse3, pclmul, ssse3, fma, cx16, sse4.1, sse4.2, movbe, popcnt, aes, xsave, avx, f16c, rdrnd, fsgsbase, bmi, avx2, bmi2, sahf, lzcnt, xsaveopt" |> x -> split(x, ", ") |> sort .|> uppercase;
end

fs_jldump = jldumphost()

# CpuId
fs_cpuid = String.(cpuidfeatures());

export fs_cpuid, fs_jldump, fs_proc

end # module


text:
"Useful constants for atmospheric science on Earth. Uses Unitful constants."
module constants_unitful

using Unitful

#export UniversalGasConst, AvogadroNumber, BoltzmannConstant, SpeedOfLight
#export PermittivityOfVacuum, PlanckConstant
#export AirDryMolecularWeight, AirDrySpecificHeatConstPressure, AirDrySpecificHeatConstVolume
#export AirDryGasConst, AirDryThermalConductivity
#export Cp, Cpv, Cw
#export H2OMolecularWeight, H2OGasConst, LatentHeatVapor, LatentHeatFusion
#export EarthGravity, EarthAngularVelocity, EarthSunSurfaceDistance
#export EarthRadius, SolarFluxTOA
#export p0, KelvinCelsius

# physical and chemical constants
UniversalGasConst=Unitful.R # 8314.3 # J/K/kmol
AvogadroNumber=Unitful.Na # 6.022e26 # Avogadro's number
BoltzmannConstant=Unitful.k # 1.381e-23 # J/K/molecule Boltzmann's constant
SpeedOfLight=Unitful.c0 #2.998e8 # m/s
PermittivityOfVacuum=Unitful.ϵ0 # 8.85e-12  # C^-2 N^-1 m^-2 permittivity of vacuum
PlanckConstant=Unitful.h # 6.6262e-34 # J*s

# properties of dry air
AirDryMolecularWeight=28.97u"kg/kmol"
AirDrySpecificHeatConstPressure=1004.0u"J/K/kg" # Wallace and Hobbs
AirDrySpecificHeatConstVolume=717.0u"J/K/kg"
#AirDryGasConst=287.0u"J/K/kg"
AirDryThermalConductivity=2.40e-2u"J/m/s/K" # at 0C
Rd=287.04u"J/K/kg" # J/K/kg Bolton
Cp=1005.7u"J/K/kg" # J/K/kg Bolton

# properties of water
H2OMolecularWeight=18.016u"g/mol"
H2OGasConst=461u"J/K/kg"
LatentHeatVapor=2.5e6u"J/kg" # at 0 C
LatentHeatFusion=3.34e5u"K/kg"
Rv=461.5u"J/kg/K" # Bolton
# specific heats of vapor and liquid
Cpv=1870u"J/kg/K" # Bolton
Cw=4190u"J/kg/K"  # Bolton

# properties of Earth
EarthGravity=9.8u"m/s/s"            # m/s/s on gravity on Earth surface
EarthAngularVelocity=7.292e-5u"s^-1"
EarthSunSurfaceDistance=1.50e11u"m"
EarthRadius=6.37e6u"m"
SolarFluxTOA=1.38e3u"W/m^2"

# reference values
#hPa=100u"Pa"
pref=1.e5u"Pa"
#KelvinCelsius=273.15 # K
end # module constants_unitful

module constants_unitless

#export UniversalGasConst, AvogadroNumber, BoltzmannConstant
#export PermittivityOfVacuum, PlanckConstant, SpeedOfLight
#export AirDryMolecularWeight, AirDrySpecificHeatConstPressure, AirDrySpecificHeatConstVolume
#export AirDryThermalConductivity
#export Rd, Cp, Cpv, Cw
#export H2OMolecularWeight, H2OGasConst, LatentHeatVapor, LatentHeatFusion
#export EarthGravity, EarthAngularVelocity, EarthRadius
#export EarthSunSurfaceDistance, SolarFluxTOA

# physical and chemical constants
UniversalGasConst=8314.31 #J/K/kmol
AvogadroNumber=6.022e26 # Avogadro's number
BoltzmannConstant=1.381e-23 # J/K/molecule Boltzmann's constant
SpeedOfLight=2.998e8 # m/s
PermittivityOfVacuum=8.85e-12  # C^-2 N^-1 m^-2 permittivity of vacuum
PlanckConstant=6.6262e-34 # J/s

# properties of dry air
AirDryMolecularWeight=28.97 # kg/kmol
AirDrySpecificHeatConstPressure=1004.  # J/K/kg, Wallace and Hobbs
AirDrySpecificHeatConstVolume=717.   # J/K/kg
AirDryGasConst=287.   # J/K/kg
AirDryThermalConductivity=2.40e-2 # J/m/s/K at 0C
Rd=287.04 # J/K/kg Bolton
Cp=1005.7 # J/K/kg Bolton

# properties of water
H2OMolecularWeight=18.016 # kg/kmol
H2OGasConst=461        # J/K/kg
LatentHeatVapor=2.5e6  # J kg-1 at 0 C
LatentHeatFusion=3.34e5 # K kg-1
Rv=461.5 # Bolton
# specific heats of vapor and liquid
Cpv=1870 # Bolton
Cw=4190  # Bolton

# properties of Earth
EarthGravity=9.8                    # m/s/s on gravity on Earth surface
EarthAngularVelocity=7.292e-5       # s-1
EarthSunSurfaceDistance=1.50e11     # m
EarthRadius=6.37e6                  # m
SolarFluxTOA=1.38e3                 # W/m2

# reference values
pref=1.e5   # hPa
KelvinCelsius=273.15 # K
end # constants unitless


text:
# This file was derived from "docs/src/assets/figures/generate-named-color-charts.jl" (v0.9.6).

module NamedColorCharts

using Colors

struct ColorChartSVG <: Main.SVG
    buf::IOBuffer
end

function ColorChartSVG(colorcategory)
    io = IOBuffer()

    colornames = colordictionary[colorcategory]
    numbercells = length(colornames)
    numbercols = 10
    colwidth = 14.8 #mm
    rowheight = 15 # mm
    swatchheight = 9 # mm
    pagewidth = 150 # mm
    margin = round(pagewidth - colwidth * numbercols) / 2
    numberrows = convert(Int, ceil(numbercells/numbercols))
    headerheight = 15 # mm
    pageheight = headerheight + numberrows * rowheight + margin
    write(io,
        """
        <svg xmlns="http://www.w3.org/2000/svg" version="1.1"
             viewBox="0 0 $pagewidth $pageheight"
             width="$(pagewidth)mm" height="$(pageheight)mm"
             style="width:100%; height:auto;"
             shape-rendering="crispEdges" stroke="none">
        <defs>
            <style><![CDATA[
                @import url('https://fonts.googleapis.com/css?family=IBM+Plex+Sans+Condensed&display=swap');
                g text:last-child {
                    font-size: 2px;
                    fill: #000;
                }
                g text {
                    text-anchor: middle;
                    font-size: 2.5px;
                    font-weight: 500;
                    font-family: 'IBM Plex Sans Condensed', 'Futura Condensed Medium', 'Segoe UI', 'Roboto Condensed', sans-serif;
                    fill: currentColor;
                }
                .header {
                    font-size: 5px;
                    text-anchor: start;
                    font-family: 'Lato', 'Helvetica Neue', 'Arial', sans-serif;
                    font-weight: bold;
                    fill: currentColor;
                }
                .w {
                    fill: #fff !important;
                }
            ]]></style>
        </defs>
        """)
    write(io,
        """
        <text x="$margin" y="10" class="header">$(titlecase(colorcategory))</text>
        <path d="M0,12 h$(pagewidth)" stroke="#888" stroke-width="0.25" />
        """)

    for n = 1:numbercells
        name = colornames[n]
        col = parse(RGB, name)
        x = round(margin + colwidth * ((n - 1) % numbercols), digits=1)
        y = round(headerheight + rowheight * ((n - 1) ÷ numbercols), digits=1)
        cx = round(colwidth / 2, digits=1)

        write(io, """<g transform="translate($x,$y)">""")
        write(io, """<rect y="5" width="$colwidth" height="$swatchheight" fill="#$(hex(col))" />""")
        write(io, "<text>")

        len = length(name)
        if len <= 12
            write(io, """<tspan x="$cx" y="4.2">$name</tspan>""")
        else
            m = match(r"^(dark|pale|light|lemon|medium|blanched|(?:.+(?=white|blue|purple|blush)))(.+)$", name)
            if m !== nothing
                write(io, """<tspan x="$cx" y="2.3">$(m.captures[1])</tspan>""")
                write(io, """<tspan x="$cx" y="4.5">$(m.captures[2])</tspan>""")
            else
                write(io, """<tspan x="$cx" y="2.3">$(name[1:(len÷2)])</tspan>""")
                write(io, """<tspan x="$cx" y="4.5">$(name[(len÷2+1):end])</tspan>""")
            end
        end
        write(io, "</text>")

        vec = string(round(red(col), digits=2), ", ",
                     round(green(col), digits=2), ", ",
                     round(blue(col), digits=2))
        veccol = convert(Lab, col).l > 60 ? "" : " class=\"w\""
        write(io, """<text x="$cx" y="13"$veccol>$vec</text>""")
        write(io, "</g>\n")
    end

    write(io, "</svg>")
    ColorChartSVG(io)
end

function compare_colors(colorname_a, colorname_b)
    # compare two colors, looking just at their LUV luminance values
    luv1 = parse(Luv, colorname_a)
    luv2 = parse(Luv, colorname_b)
    luv1.l == luv2.l ? colorname_a < colorname_b : luv1.l > luv2.l
end

# sort colors into categories
# these are very bikesheddable
function classifycolornames(colornames)
    colorcategories = Dict{String, Array{String, 1}}(
        "whites" => [],
        "pinks" => [],
        "reds" => [],
        "oranges" => [],
        "yellows" => [],
        "greens" => [],
        "cyans" => [],
        "blues" => [],
        "purples" => [],
        "browns" => [],
        "grays" => [])
    for colorname in colornames
        if occursin(r"(^grey)|(^gray)", colorname)
            push!(colorcategories["grays"], colorname)
        elseif occursin(r"turquoise|^aqua$|teal|cyan", colorname)
            push!(colorcategories["cyans"], colorname)
        elseif occursin(r"lemon|khaki|wheat", colorname)
            push!(colorcategories["yellows"], colorname)
        elseif occursin(r"orange", colorname)
            push!(colorcategories["oranges"], colorname)
        elseif occursin(r"chartreuse|olive|marine|lime|green", colorname)
            push!(colorcategories["greens"], colorname)
        elseif occursin(r"maroon|violetred|pink", colorname)
            push!(colorcategories["pinks"], colorname)
        elseif occursin(r"mistyrose|plum|thistle|lavender|violet|orchid|magenta|fuchsia|purple", colorname)
            push!(colorcategories["purples"], colorname)
        elseif occursin(r"azure|gains|slate|navy|indigo|blue", colorname)
            push!(colorcategories["blues"], colorname)
        elseif occursin(r"gold|yellow", colorname)
            push!(colorcategories["yellows"], colorname)
        elseif occursin(r"chocolate|sienna|wood|moccasin|bisque|peru|peach|papaya|almond|tan|brown", colorname)
            push!(colorcategories["browns"], colorname)
        elseif occursin(r"crimson|firebrick|salmon|coral|tomato|red", colorname)
            push!(colorcategories["reds"], colorname)
        elseif occursin(r"honey|mint|beige|linen|corn|sea|silver|lace|ivory|snow|white", colorname)
            push!(colorcategories["whites"], colorname)
        elseif occursin(r"grey|gray|black", colorname)
            push!(colorcategories["grays"], colorname)
        else
            error(colorname)
        end
    end
    # sort the RGB values in Luv.l format
    # TODO perhaps do this conversion just once instead of 1000s of times :)
    for k in keys(colorcategories)
        vals = colorcategories[k]
        sort!(vals, lt = compare_colors)
        colorcategories[k] = vals
    end
    return colorcategories
end

function makecolordictionary()
    colornames = collect(keys(Colors.color_names))
    colordictionary = classifycolornames(colornames)
    checkclassified = String[]
    for category in keys(colordictionary)
        for c in colordictionary[category]
            push!(checkclassified, c)
        end
    end
    stilltodo = setdiff(colornames, checkclassified)
    length(stilltodo) > 0 && error("gotta catch them all!")
    return colordictionary
end

const colordictionary = makecolordictionary()

end


text:
include("../src/QualitativeAudioFeatures.jl")

using Test, PortAudio, WAV, Glob
using .QualitativeAudioFeatures

@testset "System Test" begin
  @testset "Live Extraction" begin
    audio_device = "system"
    number_of_tracks = 1
    # stream = PortAudioStream(audio_device, number_of_tracks)

    # for i in 1:1000
    #   audio_data = read(stream, 1s)
    # end
  end

  @testset "WAV Files" begin
    wav_dir = "../data/reaper-files"
    @testset "Location" begin
      pan_files = glob("*.wav", wav_dir)
    end
  end

end


text:
# Declare a new `isfinite` function to avoid type piracy. This new function works with
# `Char` and `Period` as well as other types.
isfinite(x) = iszero(x - x)
isfinite(x::Real) = Base.isfinite(x)
isfinite(x::Union{Type{T}, T}) where T<:TimeType = Base.isfinite(x)


text:
n = 200
m = 100
p = 25

info("Testing ", MOD.kernel)
for kernelobj in (additive_kernels..., composition_kernels...)
    for T in FloatingPointTypes
        x = rand(T,p)
        y = rand(T,p)
        k = convert(Kernel{T}, (kernelobj)())

        if isa(k, PairwiseKernel)
            @test MOD.kernel(k, x[1], y[1]) == MOD.phi(k, x[1], y[1])
            @test MOD.kernel(k, x, y)       == sum(map((x,y) -> MOD.phi(k,x,y), x, y))
        else
            @test MOD.kernel(k, x[1], y[1]) == MOD.phi(k.phi, MOD.phi(k.kappa, x[1], y[1]))
            @test MOD.kernel(k, x, y)       == MOD.phi(k.phi, sum(map((x,y) -> MOD.phi(k.kappa,x,y), x, y)))
        end
    end
end

info("Testing ", MOD.kernelmatrix)
steps = length(additive_kernels) + length(composition_kernels) + 1
counter = 0
info("    Progress:   0.0%")
for kernelobj in (additive_kernels..., composition_kernels..., MLKTest.PairwiseTestKernel)
    for T in FloatingPointTypes
        Set_X = [rand(T,p) for i = 1:n]
        Set_Y = [rand(T,p) for i = 1:m]

        X = transpose(hcat(Set_X...))
        Y = transpose(hcat(Set_Y...))
        k = convert(Kernel{T}, (kernelobj)())
        #println("κ:", k, ", X:", typeof(X), ", Y:", typeof(Y))

        P = [MOD.kernel(k,x,y) for x in Set_X, y in Set_X]
        @test_approx_eq MOD.kernelmatrix!(Val{:row}, Array(T,n,n), k, X)  P
        @test_approx_eq MOD.kernelmatrix!(Val{:col}, Array(T,n,n), k, X') P
        @test_approx_eq MOD.kernelmatrix(Val{:row}, k, X)  P
        @test_approx_eq MOD.kernelmatrix(Val{:col}, k, X') P
        @test_approx_eq MOD.kernelmatrix(k, X) P

        P = [MOD.kernel(k,x,y) for x in Set_X, y in Set_Y]
        @test_approx_eq MOD.kernelmatrix!(Val{:row}, Array(T,n,m), k, X,  Y)  P
        @test_approx_eq MOD.kernelmatrix!(Val{:col}, Array(T,n,m), k, X', Y') P
        @test_approx_eq MOD.kernelmatrix(Val{:row}, k, X,  Y)  P
        @test_approx_eq MOD.kernelmatrix(Val{:col}, k, X', Y') P
        @test_approx_eq MOD.kernelmatrix(k, X,  Y) P
    end
    counter += 1
    info("    Progress: ", @sprintf("%5.1f", counter/steps*100), "%")
end

#=
info("Testing ", centerkernelmatrix)
for T in (Float32, Float64)
    A = T[1 2 3;
          2 3 4;
          3 4 5]

    a = mean(A,1)

    @test_approx_eq centerkernelmatrix(A) ((A .- a) .- a') .+ mean(A)
end
=#


text:
abstract Kernel
type Uniform<:Kernel; end
(::Uniform){T}(x::T) = abs(x) < one(T) ? 1/2 : 0
type Triangular<:Kernel; end
(::Triangular){T}(x::T) = abs(x) < one(T) ? one(T)-abs(x) : zero(T)
type Gaussian<:Kernel; end
(::Gaussian)(x) = normpdf(x)
type Epanechnikov<:Kernel; end
(::Epanechnikov){T}(x::T) = abs(x) < one(T) ? T(3/4)*(one(T)-x^2) : zero(T)
type Biweight<:Kernel; end
(::Biweight){T}(x::T) = abs(x) < one(T) ? T(15/16)*(one(T)-x^2)^2 : zero(T)
type Triweight<:Kernel; end
(::Triweight){T}(x::T) = abs(x) < one(T) ? T(35/32)*(one(T)-x^2)^3 : zero(T)
type Tricube<:Kernel; end
(::Tricube){T}(x::T) = abs(x) < one(T) ? T(70/81)*(one(T)-abs(x)^3)^3 : zero(T)
type Cosine<:Kernel; end
(::Cosine){T}(x::T) = abs(x) < one(T) ? T(π/4)*cos(x*π/2) : zero(T)
type Logistic<:Kernel; end
(::Logistic){T}(x::T) = 1/(2+exp(x)+exp(-x))
type Sigmoid<:Kernel; end
(::Sigmoid){T}(x::T) = 2/(π(exp(x)+exp(-x)))
type Silverman<:Kernel; end
(::Silverman){T}(x::T) = 1/2*exp(-abs(x)/sqrt(2))*sin(abs(x)/sqrt(2)+π/4)


text:
"""
    posterior = bpca_ais(X, Σ; [maxiter=1000], [S=300])

`X` is a D×N matrix of flux observations.

`σ` is a D×N matrix of error measurements.

Arguments in brackets are optional.
Default number of iterations is 1200.
Default number of samples for Adaptive Importance Sampling is 300.

Returns posterior distribution of PCA parameters returned as a ```Distributions.MvNormal``` type.
"""
function bpcaAIS(X, σ; maxiter=1000, S=300)

    Σ = σ.^2

    #---------------------------------------------
    # Set constants and report message
    #---------------------------------------------

    Q = 1

    JITTER = 1e-8

    D, N = size(X)

    @printf("BPCAAIS: %d number of data items of dimension %d, projected to %d components\n", N, D, Q)

    #---------------------------------------------
    function logPrior(W, μ)
    #---------------------------------------------

        -1e-10*(sum(W.^2) + sum(μ.^2))

    end

    #---------------------------------------------
    function marginalLogLikelihood(W, μ)
    #---------------------------------------------

       local ℓ = zero(eltype(W))

       local WWT = W*W'

       @inbounds for n in 1:N

           # check eqs. (12.36), (12.43) in Bishop
           local C = WWT + Diagonal(Σ[:,n]) + JITTER*I

           ℓ +=  - 0.5*logdet(C) - 0.5*sum((X[:,n].-μ)'*(C\(X[:,n].-μ)))

       end

       return ℓ

    end

    #---------------------------------------------
    function unpack(param)
    #---------------------------------------------

        local MARK = 0

        local W = reshape(param[MARK+1:MARK+D*Q], D, Q)

        MARK += D*Q

        local μ = reshape(param[MARK+1:MARK+D], D)

        MARK += D

        @assert(MARK == length(param))

        return W, μ

    end

    #---------------------------------------------
    function objective(param)
    #---------------------------------------------

        @assert(length(param) == D*Q + D)

        local W, μ = unpack(param)

        return marginalLogLikelihood(W, μ) + logPrior(W, μ)

    end

    #---------------------------------------------
    # Initialise with ppca estimates
    #---------------------------------------------

    W, μ, =  ppca(X)

    #---------------------------------------------
    # Run AIS
    #---------------------------------------------

    return ais(objective, [vec(W);μ], S, maxiter)

end


text:
using ApproxFun

#The following samples a 2D Cauchy Distribution

f = Fun((x,y)->1/(2π*(x^2+y^2+1)^(3/2)),JacobiWeight(2.,2.,Line())^2)
x = ApproxFun.sample(f,100)

max_stars_repo_path: examples/2D Cauchy Distribution.jl

max_issues_repo_path: examples/2D Cauchy Distribution.jl

max_forks_repo_path: examples/2D Cauchy Distribution.jl


text:
function ADScalarObjective(f,x0::AbstractArray)
    Hres = DiffResults.HessianResult(x0)
    function _g(df,x,Hresult)
        ForwardDiff.gradient!(Hresult,f,x)
        df .= DiffResults.gradient(Hresult)
        df
    end
    function _fg(df,x,Hresult)
        ForwardDiff.gradient!(Hresult,f,x)
        df .= DiffResults.gradient(Hresult)
        fx = DiffResults.value(Hresult)
        return fx,df
    end

    function _fgh(df,d2f,x,Hresult)
        ForwardDiff.hessian!(Hresult,f,x)
        d2f .= DiffResults.hessian(Hresult)
        df .= DiffResults.gradient(Hresult)
        fx = DiffResults.value(Hresult)
        return fx,df,d2f
    end

    function h(d2f,x)
        ForwardDiff.hessian!(d2f,f,x)
        d2f
    end
    g(df,x) = _g(df,x,Hres)
    fg(df,x) = _fg(df,x,Hres)
    fgh(df,d2f,x) =  _fgh(df,d2f,x,Hres)
    return ScalarObjective(f=f,
    g=g,
    fg=fg,
    fgh=fgh,
    h=h)
end

function ADScalarObjective(f,x0::Number)
    function g(x)
        return ForwardDiff.derivative(f,x)
    end

    function fg(∂fx,x)
        ∂fx,x = f∂f(f,x)
        return ∂fx,x
    end

    function fgh(∂fx,∂2fx,x)
        fx,∂fx,∂2fx = f∂f∂2f(f,x)
        return fx,∂fx,∂2fx
    end

    function h(∂2fx,x)
        ∂2fx = ForwardDiff.derivative(g,x)
        return ∂2fx
    end

    return ScalarObjective(f=f,
    g=g,
    fg=fg,
    fgh=fgh,
    h=h)
end
"""
    function optimize(f,x0,method=LineSearch(Newton()), options=OptimizationOptions())
"""
function optimize(f,x0::Number,method=LineSearch(Newton()),options=OptimizationOptions())
    scalarobj = ADScalarObjective(f,x0)
    optprob = OptimizationProblem(scalarobj; inplace=false)
    return NLSolvers.solve(optprob, x0, method,options)
end

#this should allow to use another schemes to update
struct ModifiedActiveBox{M,T}
    scheme::M
    ϵ::T
end

Base.summary(ab::ModifiedActiveBox) = "Active Box, using $(Base.summary(ab.scheme))"

NLSolvers.modelscheme(ab::ModifiedActiveBox) = ab.scheme
ModifiedActiveBox(scheme=Newton(); epsilon=1e-12) =  ModifiedActiveBox(scheme,epsilon)

function box_optimize(f,x0,lb,ub,options=OptimizationOptions();tolbounds = 1e-12)
    scalarobj = ADScalarObjective(f,x0)
    prob = OptimizationProblem(obj=scalarobj,
    bounds=(lb,ub); inplace=true)
    NLSolvers.solve(prob, x0, ModifiedActiveBox(Newton(linsolve = cholesky_linsolve);epsilon=tolbounds), options)
end

#see https://github.com/JuliaNLSolvers/NLSolvers.jl/issues/21
function NLSolvers.solve(prob::OptimizationProblem, x0, scheme::ModifiedActiveBox, options::OptimizationOptions)
    t0 = time()

    x0, B0 = x0, (false*x0*x0' +I)
    lower, upper = NLSolvers.bounds(prob)
    ϵbounds = mapreduce(b->(b[2] - b[1])/2, min, zip(lower, upper)) # [1, pp. 100: 5.41]

    !any(clamp.(x0, lower, upper) .!= x0) || error("Initial guess not in the feasible region")

    linesearch = NLSolvers.ArmijoBertsekas()
    mstyle = NLSolvers.OutOfPlace()
    direction_scheme = NLSolvers.modelscheme(scheme)
    objvars = NLSolvers.prepare_variables(prob, scheme, x0, copy(x0), B0)
    f0, ∇f0 = objvars.fz, norm(objvars.∇fz, Inf) # use user norm
    fz, ∇fz = objvars.fz, objvars.∇fz # use user norm
    fx, ∇fx = fz, copy(∇fz)
    B = B0
    x, z = copy(x0), copy(x0)
    Tf = typeof(fz)
    is_first=false
    Ix = Diagonal(z.*0 .+ 1)

    unconstrained_d = copy(∇fz)
    d = copy(∇fz)
    dmin = copy(∇fz)
    y = nothing
    s = copy(∇fz)
    xnorm = copy(∇fz)
    activeset = similar(∇fz,Bool)
    Ĥ = copy(B)
    Ĥchol = cholesky(Ĥ)
    for iter = 1:options.maxiter
        x .= z
        fx = fz
        ∇fx .= ∇fz

        ϵ = min(norm(clamp.(x.-∇fx, lower, upper).-x), ϵbounds) # Kelley 5.41 and just after (83) in [1]
        activeset .= NLSolvers.is_ϵ_active.(x, lower, upper, ∇fx, ϵ)
        Ĥ .= NLSolvers.diagrestrict.(B, activeset, activeset', Ix)
        #d = clamp.(x.-HhatChol\∇fx, lower, upper).-x
        # Update current gradient and calculate the search direction
        NLSolvers.find_direction!(unconstrained_d,Ĥ,objvars.Pg,∇fx,direction_scheme)
        d .= clamp.(x .+ unconstrained_d, lower, upper).-x #clamp unconstrained direction
        dmin .= d./unconstrained_d
        d_scale = minimum(dmin)
        #scale the direction, the original direction is the same.
        #d .=  @. unconstrained_d*d_scale + (one(d_scale)-d_scale)*d
        !iszero(d_scale) && (d .= unconstrained_d .* d_scale)
        #if the scale is zero, we are actually in a boundary, and we can only grow parallel along the box.
        φ = NLSolvers._lineobjective(mstyle, prob, ∇fz, z, x, d, fz, dot(∇fz, d))

        # Perform line search along d
        # Also returns final step vector and update the state
        α, f_α, ls_success = NLSolvers.find_steplength(mstyle, linesearch, φ, Tf(1), ∇fz, activeset, lower, upper, x, d, ∇fx, activeset)
        # # Calculate final step vector and update the state
        if !ls_success
            return NLSolvers.ConvergenceInfo(scheme, (prob=prob, B=B, ρs=norm(x.-z), ρx=norm(x), minimizer=z, fx=fx, minimum=fz, ∇fz=∇fz, f0=f0, ∇f0=∇f0, iter=iter, time=time()-t0), options)
        end

        s .= @. α * d
        z .= @. x + s
        s .= clamp.(z, lower, upper) .- x
        z .= x .+ s
        # Update approximation. it only works with newton at the moment
        fz, ∇fz, B, s, y = NLSolvers.update_obj!(prob.objective,s,y, ∇fx, z, ∇fz, B, direction_scheme)
        #fz, ∇fz, B, s, y = NLSolvers.update_obj(prob.objective, s, ∇fx, z, ∇fz, B, direction_scheme,is_first)
        #@show fz, ∇fz, B, s, y
        #@show Ĥ
        xnorm .= x.-clamp.(x.-∇fz, lower, upper)
        normx = norm(xnorm, Inf)
        if normx < options.g_abstol
            return NLSolvers.ConvergenceInfo(scheme, (prob=prob, B=B, ρs=norm(x.-z), ρx=norm(x), minimizer=z, fx=fx, minimum=fz, ∇fz=∇fz, f0=f0, ∇f0=∇f0, iter=iter, time=time()-t0), options)
        end
    end
  iter =  options.maxiter
  return NLSolvers.ConvergenceInfo(scheme, (prob=prob, B=B, ρs=norm(x.-z), ρx=norm(x), minimizer=z, fx=fx, minimum=fz, ∇fz=∇fz, f0=f0, ∇f0=∇f0, iter=iter, time=time()-t0), options)
end

#=
julia> Optim.minimizer(results)
2-element Vector{Float64}:
 1.2500000000000002
 1.5626548683420372

julia> Optim.minimum(results)
0.06250239842033663
=#


text:
# Common stuff
## convenient types and type conversion functions

const GChars = Union{Cchar, Char}
const IVec  = Vector{Cint}
const FVec = Vector{Float64}
const CVec = Vector{Cchar}

const GCharOrVec = Union{Cchar, Char, Vector{Cchar}, Vector{Char}}

const CoeffMat = Union{Matrix{Float64}, SparseMatrixCSC{Float64}}

cchar(c::Cchar) = c
cchar(c::Char) = convert(Cchar, c)

ivec(v::IVec) = v
fvec(v::FVec) = v
cvec(v::CVec) = v

ivec(v::Vector) = convert(IVec, v)
fvec(v::Vector) = convert(FVec, v)
cvec(v::Vector) = convert(CVec, v)

ivec(v::Int) = Cint[v]
fvec(v::Float64) = Float64[v]
fvec(v::Integer) = Float64[v]
#cvec(v::Vector) = convert(CVec, v)

# cvecx(v, n) and fvecx(v, n)
# converts v into a vector of Cchar or Float64 of length n,
# where v can be either a scalar or a vector of length n.

_chklen(n::Integer, v::Vector) = _chklen(v, n::Integer)
_chklen(v, n::Integer) = (length(v) == n || error("Inconsistent argument dimensions."))
_cmplen(v1::Vector, v2::Vector) = (length(v1) == length(v2) || error("Inconsistent argument dimensions."))

cvecx(c::GChars, n::Integer) = fill(cchar(c), n)
cvecx(c::Vector{Cchar}, n::Integer) = (_chklen(c, n); c)
cvecx(c::Vector{Char}, n::Integer) = (_chklen(c, n); convert(Vector{Cchar}, c))

fvecx(v::Real, n::Integer) = fill(Float64(v), n)
fvecx(v::Vector{Float64}, n::Integer) = (_chklen(v, n); v)
fvecx(v::Vector{T}, n::Integer) where {T<:Real} = (_chklen(v, n); convert(Vector{Float64}, v))

inds32(n::Integer) = collect(Cint(1):Cint(n))

# empty vector & matrix (for the purpose of supplying default arguments)

const emptyfvec = Array{Float64}(undef, 0)
const emptyfmat = Array{Float64}(undef, 0, 0)

# macro to call a Xpress C functions
macro xprs_ccall(func, args...)
    f = "XPRS$(func)"
    args = map(esc,args)

    Compat.Sys.isunix() && return quote
        ccall(($f,xprs), $(args...))
    end
    Compat.Sys.iswindows() && VERSION < v"0.6-" && return quote
        ccall(($f,xprs), stdcall, $(args...))
    end
    Compat.Sys.iswindows() && VERSION >= v"0.6-" && return quote
        ccall(($f,xprs), $(esc(:stdcall)), $(args...))
    end
end

"""
    getlibversion()

Get Xpress optimizer version info
"""
function getlibversion()
    out = Array{Cchar}(undef, 16)                                     # "                "
    ret = @xprs_ccall(getversion, Cint, ( Ptr{Cchar},), out)   # ( Cstring,), out)

    numbers = split(unsafe_string(pointer(out)) ,".")

    _major = parse(Int,numbers[1])
    _minor = parse(Int,numbers[2])
    _tech  = parse(Int,numbers[3])

    return  VersionNumber(_major, _minor, _tech)
end
#const version = getlibversion()

function fixinfinity(val::Float64)
    if val == Inf
        return XPRS_PLUSINFINITY
    elseif val == -Inf
        return XPRS_MINUSINFINITY
    else
        return val
    end
end
function fixinfinity!(vals::Vector{Float64})
    map!(fixinfinity, vals, vals)
end


text:
""" Data structures for graphs, based on C-sets.

Provides the category theorist's four basic kinds of graphs: graphs (aka
directed multigraphs), symmetric graphs, reflexive graphs, and symmetric
reflexive graphs. Also defines half-edge graphs. The API generally follows that
of [LightGraphs.jl](https://github.com/JuliaGraphs/LightGraphs.jl), with some
departures due to differences between the data structures.
"""
module BasicGraphs
export HasVertices, HasGraph,
  AbstractGraph, Graph, nv, ne, src, tgt, edges, vertices,
  has_edge, has_vertex, add_edge!, add_edges!, add_vertex!, add_vertices!,
  rem_edge!, rem_edges!, rem_vertex!, rem_vertices!,
  neighbors, inneighbors, outneighbors, all_neighbors, induced_subgraph,
  AbstractSymmetricGraph, SymmetricGraph, inv,
  AbstractReflexiveGraph, ReflexiveGraph, refl,
  AbstractSymmetricReflexiveGraph, SymmetricReflexiveGraph,
  AbstractHalfEdgeGraph, HalfEdgeGraph, vertex, half_edges,
  add_dangling_edge!, add_dangling_edges!,
  AbstractWeightedGraph, WeightedGraph, weight,
  AbstractSymmetricWeightedGraph, SymmetricWeightedGraph

import Base: inv
using Requires

using ...Present, ...CSetDataStructures, ...Theories

# Base types
############

""" Abstract type for C-sets that contain vertices.

This type encompasses C-sets where the schema C contains an object `V`
interpreted as vertices. This includes, for example, graphs and half-edge
graphs, but not bipartite graphs or wiring diagrams.
"""
@abstract_acset_type HasVertices

""" Abstract type for C-sets that contain a graph.

This type encompasses C-sets where the schema for graphs is a subcategory of C.
This includes, for example, graphs, symmetric graphs, and reflexive graphs, but
not half-edge graphs.
"""
@abstract_acset_type HasGraph <: HasVertices

# Graphs
########

@present TheoryGraph(FreeSchema) begin
  V::Ob
  E::Ob
  src::Hom(E,V)
  tgt::Hom(E,V)
end

""" Abstract type for graphs, aka directed multigraphs.
"""
@abstract_acset_type AbstractGraph <: HasGraph

""" A graph, also known as a directed multigraph.
"""
@acset_type Graph(TheoryGraph, index=[:src,:tgt]) <: AbstractGraph

function (::Type{T})(nv::Int; kw...) where T <: HasVertices
  g = T()
  add_vertices!(g, nv; kw...)
  return g
end

""" Number of vertices in a graph.
"""
nv(g::HasVertices) = nparts(g, :V)

""" Number of edges in a graph, or between two vertices in a graph.

In a symmetric graph, this function counts both edges in each edge pair, so that
the number of edges in a symmetric graph is twice the number of edges in the
corresponding undirected graph (at least when the edge involution has no fixed
points).
"""
ne(g::HasGraph) = nparts(g, :E)
ne(g::HasGraph, src::Int, tgt::Int) =
  count(subpart(g, e, :tgt) == tgt for e in incident(g, src, :src))

""" Source vertex (vertices) of edges(s) in a graph.
"""
src(g::HasGraph, args...) = subpart(g, args..., :src)

""" Target vertex (vertices) of edges(s) in a graph.
"""
tgt(g::HasGraph, args...) = subpart(g, args..., :tgt)

""" Vertices in a graph.
"""
vertices(g::HasVertices) = parts(g, :V)

""" Edges in a graph, or between two vertices in a graph.
"""
edges(g::HasGraph) = parts(g, :E)
edges(g::HasGraph, src::Int, tgt::Int) =
  (e for e in incident(g, src, :src) if subpart(g, e, :tgt) == tgt)

""" Whether the graph has the given vertex.
"""
has_vertex(g::HasVertices, v) = has_part(g, :V, v)

""" Whether the graph has the given edge, or an edge between two vertices.
"""
has_edge(g::HasGraph, e) = has_part(g, :E, e)
has_edge(g::HasGraph, src::Int, tgt::Int) =
  has_vertex(g, src) && tgt ∈ outneighbors(g, src)

""" Add a vertex to a graph.
"""
add_vertex!(g::HasVertices; kw...) = add_part!(g, :V; kw...)

""" Add multiple vertices to a graph.
"""
add_vertices!(g::HasVertices, n::Int; kw...) = add_parts!(g, :V, n; kw...)

""" Add an edge to a graph.
"""
add_edge!(g::HasGraph, src::Int, tgt::Int; kw...) =
  add_part!(g, :E; src=src, tgt=tgt, kw...)

""" Add multiple edges to a graph.
"""
function add_edges!(g::HasGraph, srcs::AbstractVector{Int},
                    tgts::AbstractVector{Int}; kw...)
  @assert (n = length(srcs)) == length(tgts)
  add_parts!(g, :E, n; src=srcs, tgt=tgts, kw...)
end

""" Remove a vertex from a graph.

When `keep_edges` is false (the default), all edges incident to the vertex are
also deleted. When `keep_edges` is true, incident edges are preserved but their
source/target vertices become undefined.
"""
rem_vertex!(g::HasVertices, v::Int; kw...) = rem_vertices!(g, v:v; kw...)

""" Remove multiple vertices from a graph.

Edges incident to any of the vertices are treated as in [`rem_vertex!`](@ref).
"""
function rem_vertices!(g::AbstractGraph, vs; keep_edges::Bool=false)
  if !keep_edges
    es = flatten([incident(g, vs, :src); incident(g, vs, :tgt)])
    rem_parts!(g, :E, unique!(sort!(es)))
  end
  rem_parts!(g, :V, vs)
end
flatten(x::AbstractVector{<:AbstractVector{T}}) where T =
  reduce(vcat, x, init=T[])

""" Remove an edge from a graph.
"""
rem_edge!(g::HasGraph, e::Int) = rem_part!(g, :E, e)
rem_edge!(g::HasGraph, src::Int, tgt::Int) =
  rem_edge!(g, first(edges(g, src, tgt)))

""" Remove multiple edges from a graph.
"""
rem_edges!(g::HasGraph, es) = rem_parts!(g, :E, es)

""" Neighbors of vertex in a graph.

In a graph, this function is an alias for [`outneighbors`](@ref); in a symmetric
graph, a vertex has the same out-neighbors and as in-neighbors, so the
distinction is moot.

In the presence of multiple edges, neighboring vertices are given *with
multiplicity*. To get the unique neighbors, call `unique(neighbors(g))`.
"""
neighbors(g::AbstractGraph, v::Int) = outneighbors(g, v)

""" In-neighbors of vertex in a graph.
"""
inneighbors(g::AbstractGraph, v::Int) = subpart(g, incident(g, v, :tgt), :src)

""" Out-neighbors of vertex in a graph.
"""
outneighbors(g::AbstractGraph, v::Int) = subpart(g, incident(g, v, :src), :tgt)

""" Union of in-neighbors and out-neighbors in a graph.
"""
all_neighbors(g::AbstractGraph, v::Int) =
  Iterators.flatten((inneighbors(g, v), outneighbors(g, v)))

""" Subgraph induced by a set of a vertices.

The [induced subgraph](https://en.wikipedia.org/wiki/Induced_subgraph) consists
of the given vertices and all edges between vertices in this set.
"""
function induced_subgraph(g::G, vs::AbstractVector{Int}) where G <: HasGraph
  vset = Set(vs)
  length(vs) == length(vset) || error("Duplicate vertices in: $vs")
  es = Iterators.filter(Iterators.flatten(incident(g, vs, :src))) do e
    tgt(g, e) ∈ vset
  end
  sub = G()
  copy_parts!(sub, g, V=vs, E=collect(Int, es))
  sub
end

# Symmetric graphs
##################

@present TheorySymmetricGraph <: TheoryGraph begin
  inv::Hom(E,E)

  compose(inv,inv) == id(E)
  compose(inv,src) == tgt
  compose(inv,tgt) == src
end

""" Abstract type for symmetric graph, possibly with data attributes.
"""
@abstract_acset_type AbstractSymmetricGraph <: HasGraph

""" A symmetric graph, or graph with an orientation-reversing edge involution.

Symmetric graphs are closely related, but not identical, to undirected graphs.
"""
@acset_type SymmetricGraph(TheorySymmetricGraph, index=[:src]) <: AbstractSymmetricGraph
# Don't index `inv` because it is self-inverse and don't index `tgt`
# because `src` contains the same information due to symmetry of graph.

""" Involution on edge(s) in a symmetric graph.
"""
inv(g::HasGraph, args...) = subpart(g, args..., :inv)

add_edge!(g::AbstractSymmetricGraph, src::Int, tgt::Int; kw...) =
  add_edges!(g, src:src, tgt:tgt; kw...)

function add_edges!(g::AbstractSymmetricGraph, srcs::AbstractVector{Int},
                    tgts::AbstractVector{Int}; kw...)
  @assert (n = length(srcs)) == length(tgts)
  k = nparts(g, :E)
  add_parts!(g, :E, 2n; src=vcat(srcs,tgts), tgt=vcat(tgts,srcs),
             inv=vcat((k+n+1):(k+2n),(k+1):(k+n)), kw...)
end

function rem_vertices!(g::AbstractSymmetricGraph, vs; keep_edges::Bool=false)
  if !keep_edges
    es = flatten(incident(g, vs, :src))
    rem_parts!(g, :E, unique!(sort!([es; inv(g, es)])))
  end
  # FIXME: Vertex removal is inefficient because `rem_parts!` still searches for
  # edges with given targets but `tgt` is not indexed.
  rem_parts!(g, :V, vs)
end

rem_edge!(g::AbstractSymmetricGraph, e::Int) = rem_edges!(g, e:e)
rem_edges!(g::AbstractSymmetricGraph, es) =
  rem_parts!(g, :E, unique!(sort!([es; inv(g, es)])))

neighbors(g::AbstractSymmetricGraph, v::Int) =
  subpart(g, incident(g, v, :src), :tgt)
inneighbors(g::AbstractSymmetricGraph, v::Int) = neighbors(g, v)
outneighbors(g::AbstractSymmetricGraph, v::Int) = neighbors(g, v)
all_neighbors(g::AbstractSymmetricGraph, v::Int) = neighbors(g, v)

# Reflexive graphs
##################

@present TheoryReflexiveGraph <: TheoryGraph begin
  refl::Hom(V,E)

  compose(refl, src) == id(V)
  compose(refl, tgt) == id(V)
end

""" Abstract type for reflexive graphs, possibly with data attributes.
"""
@abstract_acset_type AbstractReflexiveGraph <: HasGraph

""" A reflexive graph.

[Reflexive graphs](https://ncatlab.org/nlab/show/reflexive+graph) are graphs in
which every vertex has a distinguished self-loop.
"""
@acset_type ReflexiveGraph(TheoryReflexiveGraph, index=[:src,:tgt]) <: AbstractReflexiveGraph

""" Reflexive loop(s) of vertex (vertices) in a reflexive graph.
"""
refl(g::HasGraph, args...) = subpart(g, args..., :refl)

add_vertex!(g::AbstractReflexiveGraph; kw...) =
  only(add_vertices!(g, 1; kw...))

function add_vertices!(g::AbstractReflexiveGraph, n::Int; kw...)
  vs = add_parts!(g, :V, n; kw...)
  es = add_parts!(g, :E, n, src=vs, tgt=vs)
  set_subpart!(g, vs, :refl, es)
  vs
end

function rem_vertices!(g::AbstractReflexiveGraph, vs; keep_edges::Bool=false)
  es = if keep_edges
    sort(refl(g, vs)) # Always delete reflexive edges.
  else
    unique!(sort!(flatten([incident(g, vs, :src); incident(g, vs, :tgt)])))
  end
  rem_parts!(g, :E, es)
  rem_parts!(g, :V, vs)
end

# Symmetric reflexive graphs
############################

@present TheorySymmetricReflexiveGraph <: TheorySymmetricGraph begin
  refl::Hom(V,E)

  compose(refl, src) == id(V)
  compose(refl, tgt) == id(V)
  compose(refl, inv) == refl # Reflexive loop fixed by involution.
end

""" Abstract type for symmetric reflexive graphs, possibly with data attributes.
"""
@abstract_acset_type AbstractSymmetricReflexiveGraph <: HasGraph

""" A symmetric reflexive graph.

Symmetric reflexive graphs are both symmetric graphs ([`SymmetricGraph`](@ref))
and reflexive graphs ([`ReflexiveGraph`](@ref)) such that the reflexive loops
are fixed by the edge involution.
"""
@acset_type SymmetricReflexiveGraph(TheorySymmetricReflexiveGraph, index=[:src]) <:
  AbstractSymmetricReflexiveGraph

add_vertex!(g::AbstractSymmetricReflexiveGraph; kw...) =
  only(add_vertices!(g, 1; kw...))

function add_vertices!(g::AbstractSymmetricReflexiveGraph, n::Int; kw...)
  vs = add_parts!(g, :V, n; kw...)
  es = add_parts!(g, :E, n, src=vs, tgt=vs)
  set_subpart!(g, vs, :refl, es)
  set_subpart!(g, es, :inv, es)
  vs
end

add_edge!(g::AbstractSymmetricReflexiveGraph, src::Int, tgt::Int; kw...) =
  add_edges!(g, src:src, tgt:tgt; kw...)

function add_edges!(g::AbstractSymmetricReflexiveGraph,
                    srcs::AbstractVector{Int}, tgts::AbstractVector{Int}; kw...)
  @assert (n = length(srcs)) == length(tgts)
  k = nparts(g, :E)
  add_parts!(g, :E, 2n; src=vcat(srcs,tgts), tgt=vcat(tgts,srcs),
             inv=vcat((k+n+1):(k+2n),(k+1):(k+n)), kw...)
end

function rem_vertices!(g::AbstractSymmetricReflexiveGraph, vs;
                       keep_edges::Bool=false)
  es = if keep_edges
    sort(refl(g, vs)) # Always delete reflexive edges.
  else
    es = flatten(incident(g, vs, :src))
    unique!(sort!([es; inv(g, es)]))
  end
  rem_parts!(g, :E, es)
  # FIXME: Vertex removal is inefficient for same reason as `SymmetricGraph`.
  rem_parts!(g, :V, vs)
end

rem_edge!(g::AbstractSymmetricReflexiveGraph, e::Int) = rem_edges!(g, e:e)
rem_edges!(g::AbstractSymmetricReflexiveGraph, es) =
  rem_parts!(g, :E, unique!(sort!([es; inv(g, es)])))

# Half-edge graphs
##################

@present TheoryHalfEdgeGraph(FreeSchema) begin
  V::Ob
  H::Ob
  vertex::Hom(H,V)
  inv::Hom(H,H)

  compose(inv, inv) == id(H)
end

""" Abstract type for half-edge graphs, possibly with data attributes.
"""
@abstract_acset_type AbstractHalfEdgeGraph <: HasVertices

""" A half-edge graph.

[Half-edge
graphs](https://www.algebraicjulia.org/blog/post/2020/09/cset-graphs-2/) are a
variant of undirected graphs whose edges are pairs of "half-edges" or "darts".
Half-edge graphs are isomorphic to symmetric graphs but have a different data
model.
"""
@acset_type HalfEdgeGraph(TheoryHalfEdgeGraph, index=[:vertex]) <: AbstractHalfEdgeGraph

""" Incident vertex (vertices) of half-edge(s) in a half-edge graph.
"""
vertex(g::AbstractHalfEdgeGraph, args...) = subpart(g, args..., :vertex)

""" Involution on half-edge(s) in a half-edge graph.
"""
inv(g::AbstractHalfEdgeGraph, args...) = subpart(g, args..., :inv)

""" Half-edges in a half-edge graph, or incident to a vertex.
"""
half_edges(g::AbstractHalfEdgeGraph) = parts(g, :H)
half_edges(g::AbstractHalfEdgeGraph, v) = incident(g, v, :vertex)

function half_edge_pairs(g::AbstractHalfEdgeGraph, src::Int, tgt::Int)
  hs = half_edges(g, src)
  hs′ = inv(g, hs)
  has_tgt = vertex(g, hs′) .== tgt
  (hs[has_tgt], hs′[has_tgt])
end

@inline add_edge!(g::AbstractHalfEdgeGraph, src::Int, tgt::Int; kw...) =
  add_half_edge_pair!(g, src, tgt; kw...)

@inline add_edges!(g::AbstractHalfEdgeGraph, srcs::AbstractVector{Int},
                   tgts::AbstractVector{Int}; kw...) =
  add_half_edge_pairs!(g, srcs, tgts; kw...)

function add_half_edge_pair!(g::AbstractHalfEdgeGraph, src::Int, tgt::Int; kw...)
  k = nparts(g, :H)
  add_parts!(g, :H, 2; vertex=[src,tgt], inv=[k+2,k+1], kw...)
end

function add_half_edge_pairs!(g::AbstractHalfEdgeGraph, srcs::AbstractVector{Int},
                              tgts::AbstractVector{Int}; kw...)
  @assert (n = length(srcs)) == length(tgts)
  k = nparts(g, :H)
  add_parts!(g, :H, 2n; vertex=vcat(srcs,tgts),
             inv=vcat((k+n+1):(k+2n),(k+1):(k+n)), kw...)
end

""" Add a dangling edge to a half-edge graph.

A "dangling edge" is a half-edge that is paired with itself under the half-edge
involution. They are usually interpreted differently than "self-loops", i.e., a
pair of distinct half-edges incident to the same vertex.
"""
add_dangling_edge!(g::AbstractHalfEdgeGraph, v::Int; kw...) =
  add_part!(g, :H; vertex=v, inv=nparts(g,:H)+1)

""" Add multiple dangling edges to a half-edge graph.
"""
function add_dangling_edges!(g::AbstractHalfEdgeGraph, vs::AbstractVector{Int}; kw...)
  n, k = length(vs), nparts(g, :H)
  add_parts!(g, :H, n; vertex=vs, inv=(k+1):(k+n), kw...)
end

function rem_vertices!(g::AbstractHalfEdgeGraph, vs; keep_edges::Bool=false)
  if !keep_edges
    hs = flatten(incident(g, vs, :vertex))
    rem_parts!(g, :H, unique!(sort!([hs; inv(g, hs)])))
  end
  rem_parts!(g, :V, vs)
end

rem_edge!(g::AbstractHalfEdgeGraph, src::Int, tgt::Int) =
  rem_parts!(g, :H, sort!(unique(first.(half_edge_pairs(g, src, tgt)))))

rem_edge!(g::AbstractHalfEdgeGraph, h::Int) = rem_edges!(g, h:h)
rem_edges!(g::AbstractHalfEdgeGraph, hs) =
  rem_parts!(g, :H, unique!(sort!([hs; inv(g, hs)])))

# Weighted graphs
#################

@present TheoryWeightedGraph <: TheoryGraph begin
  Weight::AttrType
  weight::Attr(E,Weight)
end

""" Abstract type for weighted graphs.
"""
@abstract_acset_type AbstractWeightedGraph <: AbstractGraph

""" A weighted graph.

A graph in which every edge has a numerical weight.
"""
@acset_type WeightedGraph(TheoryWeightedGraph, index=[:src,:tgt]) <: AbstractWeightedGraph

""" Weight(s) of edge(s) in a weighted graph.
"""
weight(g::HasGraph, args...) = subpart(g, args..., :weight)

@present TheorySymmetricWeightedGraph <: TheorySymmetricGraph begin
  Weight::AttrType
  weight::Attr(E,Weight)

  compose(inv, weight) == weight
end

""" Abstract type for symmetric weighted graphs.
"""
@abstract_acset_type AbstractSymmetricWeightedGraph <: AbstractSymmetricGraph

""" A symmetric weighted graph.

A symmetric graph in which every edge has a numerical weight, preserved by the
edge involution.
"""
@acset_type SymmetricWeightedGraph(TheorySymmetricWeightedGraph, index=[:src]) <:
  AbstractSymmetricWeightedGraph

# JuliaGraphs constructors
##########################

function __init__()
  @require LightGraphs="093fc24a-ae57-5d10-9952-331d41423f4d" begin
    import .LightGraphs
    import .LightGraphs: SimpleGraph, SimpleDiGraph

    function (::Type{LG})(g::HasGraph) where LG <: Union{SimpleGraph,SimpleDiGraph}
      lg = LG(nv(g))
      for (s, t) in zip(src(g), tgt(g))
        LightGraphs.add_edge!(lg, s, t)
      end
      lg
    end

    function SimpleGraph(g::AbstractHalfEdgeGraph)
      lg = SimpleGraph(nv(g))
      for e in half_edges(g)
        e′ = inv(g,e)
        if e <= e′
          LightGraphs.add_edge!(lg, vertex(g,e), vertex(g,e′))
        end
      end
      lg
    end
  end

  @require MetaGraphs="626554b9-1ddb-594c-aa3c-2596fe9399a5" begin
    import .MetaGraphs
    import .MetaGraphs: MetaGraph, MetaDiGraph

    MetaDiGraph(g::AbstractWeightedGraph{S, Tuple{U}}) where {S,U} =
      to_weighted_metagraph(MetaDiGraph{Int,U}, g)
    MetaGraph(g::AbstractSymmetricWeightedGraph{S, Tuple{U}}) where {S,U} =
      to_weighted_metagraph(MetaGraph{Int,U}, g)

    function to_weighted_metagraph(MG::Type{<:MetaGraphs.AbstractMetaGraph}, g)
      mg = MG(nv(g))
      for (s, t, w) in zip(src(g), tgt(g), weight(g))
        MetaGraphs.add_edge!(mg, s, t, :weight, w)
      end
      mg
    end
  end
end

end


text:
# getindex STaylor1
getindex(a::STaylor1, n::Int) = a.coeffs[n+1]
getindex(a::STaylor1, u::UnitRange{Int}) = a.coeffs[u .+ 1]
getindex(a::STaylor1, c::Colon) = a.coeffs[c]
getindex(a::STaylor1, u::StepRange{Int,Int}) = a.coeffs[u .+ 1]

@inline iterate(a::STaylor1{N,T}, state=0) where {N, T<:Number} = state > N-1 ? nothing : (a.coeffs[state+1], state+1)
@inline firstindex(a::STaylor1) = 0
@inline lastindex(a::STaylor1{N,T}) where {N, T<:Number} = N-1
@inline eachindex(s::STaylor1{N,T}) where {N, T<:Number} = UnitRange(0, N-1)
@inline size(s::STaylor1{N,T}) where {N, T<:Number} = N
@inline length(s::STaylor1{N,T}) where {N, T<:Number} = N
@inline get_order(s::STaylor1{N,T}) where {N, T<:Number} = N - 1
@inline eltype(s::STaylor1{N,T}) where {N, T<:Number} = T
@inline axes(a::STaylor1) = ()

function Base.findfirst(a::STaylor1{N,T}) where {N, T<:Number}
    first = findfirst(x->!iszero(x), a.coeffs)
    isa(first, Nothing) && return -1
    return first-1
end
# Finds the last non-zero entry; extended to Taylor1
function Base.findlast(a::STaylor1{N,T}) where {N, T<:Number}
    last = findlast(x->!iszero(x), a.coeffs)
    isa(last, Nothing) && return -1
    return last-1
end
constant_term(a::STaylor1) = a[0]