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

"""
    Docs

The `Docs` module provides the `@doc` macro which can be used to set and retrieve
documentation metadata for Julia objects.

Please see the manual section on documentation for more
information.
"""
module Docs

"""
# Documentation

Functions, methods and types can be documented by placing a string before the definition:

    \"""
    # The Foo Function
    `foo(x)`: Foo the living hell out of `x`.
    \"""
    foo(x) = ...

The `@doc` macro can be used directly to both set and retrieve documentation / metadata. By
default, documentation is written as Markdown, but any object can be placed before the
arrow. For example:

    @doc "blah" ->
    function foo() ...

The `->` is not required if the object is on the same line, e.g.

    @doc "foo" foo

## Documenting objects after they are defined
You can document an object after its definition by

    @doc "foo" function_to_doc
    @doc "bar" TypeToDoc

For macros, the syntax is `@doc "macro doc" :(@Module.macro)` or `@doc "macro doc"
:(string_macro"")` for string macros. Without the quote `:()` the expansion of the macro
will be documented.

## Retrieving Documentation
You can retrieve docs for functions, macros and other objects as follows:

    @doc foo
    @doc @time
    @doc md""

## Functions & Methods
Placing documentation before a method definition (e.g. `function foo() ...` or `foo() = ...`)
will cause that specific method to be documented, as opposed to the whole function. Method
docs are concatenated together in the order they were defined to provide docs for the
function.
"""
:(Core.@doc)

include("bindings.jl")

import Base.Markdown: @doc_str, MD
import Base.Meta: quot, isexpr
import Base: Callable
import ..CoreDocs: lazy_iterpolate

export doc

# Basic API / Storage

const modules = Module[]
const META    = gensym(:meta)

meta(m::Module) = isdefined(m, META) ? getfield(m, META) : ObjectIdDict()

function initmeta(m::Module)
    if !isdefined(m, META)
        eval(m, :(const $META = $(ObjectIdDict())))
        push!(modules, m)
    end
    nothing
end

function signature!(tv, expr::Expr)
    is_macrocall = isexpr(expr, :macrocall)
    if is_macrocall || isexpr(expr, :call)
        sig = :(Union{Tuple{}})
        first_arg = is_macrocall ? 3 : 2 # skip function arguments
        for arg in expr.args[first_arg:end]
            isexpr(arg, :parameters) && continue
            if isexpr(arg, :kw) # optional arg
                push!(sig.args, :(Tuple{$(sig.args[end].args[2:end]...)}))
            end
            push!(sig.args[end].args, argtype(arg))
        end
        if isexpr(expr.args[1], :curly) && isempty(tv)
            append!(tv, tvar.(expr.args[1].args[2:end]))
        end
        for i = length(tv):-1:1
            push!(sig.args, :(Tuple{$(tv[i].args[1])}))
        end
        for i = length(tv):-1:1
            sig = Expr(:where, sig, tv[i])
        end
        sig
    elseif isexpr(expr, :where)
        append!(tv, tvar.(expr.args[2:end]))
        signature!(tv, expr.args[1])
    else
        signature!(tv, expr.args[1])
    end
end
signature!(tv, other) = :(Union{})
signature(expr::Expr) = signature!([], expr)
signature(other) = signature!([], other)

function argtype(expr::Expr)
    isexpr(expr, :(::))  && return expr.args[end]
    isexpr(expr, :(...)) && return :(Vararg{$(argtype(expr.args[1]))})
    argtype(expr.args[1])
end
argtype(other) = :Any

tvar(x::Expr)   = x
tvar(s::Symbol) = :($s <: Any)

# Docsystem types.
# ================

"""
    Docs.DocStr

Stores the contents of a single docstring as well as related metadata.

Both the raw text, `.text`, and the parsed markdown, `.object`, are tracked by this type.
Parsing of the raw text is done lazily when a request is made to render the docstring,
which helps to reduce total precompiled image size.

The `.data` fields stores several values related to the docstring, such as: path,
linenumber, source code, and fielddocs.
"""
mutable struct DocStr
    text   :: Core.SimpleVector
    object :: Nullable
    data   :: Dict{Symbol, Any}
end

function docstr(binding::Binding, @nospecialize typesig = Union{})
    for m in modules
        dict = meta(m)
        if haskey(dict, binding)
            docs = dict[binding].docs
            if haskey(docs, typesig)
                return docs[typesig]
            end
        end
    end
    error("could not find matching docstring for '$binding :: $typesig'.")
end
docstr(object, data = Dict()) = _docstr(object, data)

_docstr(vec::Core.SimpleVector, data) = DocStr(vec,            Nullable(),       data)
_docstr(str::AbstractString,    data) = DocStr(Core.svec(str), Nullable(),       data)
_docstr(object,                 data) = DocStr(Core.svec(),    Nullable(object), data)

_docstr(doc::DocStr, data) = (doc.data = merge(data, doc.data); doc)

macro ref(x)
    binding = bindingexpr(namify(x))
    typesig = signature(x)
    esc(docexpr(__source__, __module__, binding, typesig))
end

docexpr(__source__, __module__, args...) = Expr(:call, docstr, args...)

function formatdoc(d::DocStr)
    buffer = IOBuffer()
    for part in d.text
        formatdoc(buffer, d, part)
    end
    Markdown.parse(seekstart(buffer))
end
@noinline formatdoc(buffer, d, part) = print(buffer, part)

function parsedoc(d::DocStr)
    if isnull(d.object)
        md = formatdoc(d)
        md.meta[:module] = d.data[:module]
        md.meta[:path]   = d.data[:path]
        d.object = Nullable(md)
    end
    get(d.object)
end

"""
    MultiDoc

Stores a collection of docstrings for related objects, ie. a `Function`/`DataType` and
associated `Method` objects.

Each documented object in a `MultiDoc` is referred to by it's signature which is represented
by a `Union` of `Tuple` types. For example the following `Method` definition

    f(x, y) = ...

is stored as `Tuple{Any, Any}` in the `MultiDoc` while

    f(x::T, y = ?) where {T} = ...

is stored as `Union{Tuple{T, Any}, Tuple{T}} where T`.

Note: The `Function`/`DataType` object's signature is always `Union{}`.
"""
mutable struct MultiDoc
    "Ordered (via definition order) vector of object signatures."
    order::Vector{Type}
    "Documentation for each object. Keys are signatures."
    docs::ObjectIdDict

    MultiDoc() = new(Type[], ObjectIdDict())
end

# Docstring registration.
# =======================

"""
    Docs.doc!(__module__, binding, str, sig)

Adds a new docstring `str` to the docsystem of `__module__` for `binding` and signature `sig`.
"""
function doc!(__module__::Module, b::Binding, str::DocStr, @nospecialize sig = Union{})
    initmeta(__module__)
    m = get!(meta(__module__), b, MultiDoc())
    if haskey(m.docs, sig)
        # We allow for docstrings to be updated, but print a warning since it is possible
        # that over-writing a docstring *may* have been accidental.  The warning
        # is suppressed for symbols in Main, for interactive use (#23011).
        __module__ == Main || warn("replacing docs for '$b :: $sig' in module '$(__module__)'.")
    else
        # The ordering of docstrings for each Binding is defined by the order in which they
        # are initially added. Replacing a specific docstring does not change it's ordering.
        push!(m.order, sig)
    end
    m.docs[sig] = str
    str.data[:binding] = b
    str.data[:typesig] = sig
    return b
end

# Docstring lookup.
# =================

"""
    getdoc(obj)
    getdoc(obj, sig)

Return a custom docstring object associated with the object `obj` and, optionally, the tuple
type signature `sig`. See `MultiDoc` docs for an explanation of the possible values of `sig`.

The returned object can either be a markdown object generated by `Markdown.parse` or some
other custom type used to display non-markdown formatted documentation.

A return value of `nothing` can be used to signify to the docsystem that no documentation
was found for `obj`, in which case the docsystem will fall back to searching for the
`Binding` associated with `obj` instead.
"""
function getdoc end

getdoc(x, sig) = getdoc(x)
getdoc(x) = nothing

"""
    Docs.doc(binding, sig)

Returns all documentation that matches both `binding` and `sig`.

If `getdoc` returns a non-`nothing` result on the value of the binding, then a
dynamic docstring is returned instead of one based on the binding itself.
"""
function doc(binding::Binding, sig::Type = Union{})
    if defined(binding)
        result = getdoc(resolve(binding), sig)
        result === nothing || return result
    end
    results, groups = DocStr[], MultiDoc[]
    # Lookup `binding` and `sig` for matches in all modules of the docsystem.
    for mod in modules
        dict = meta(mod)
        if haskey(dict, binding)
            multidoc = dict[binding]
            push!(groups, multidoc)
            for msig in multidoc.order
                sig <: msig && push!(results, multidoc.docs[msig])
            end
        end
    end
    if isempty(groups)
        # When no `MultiDoc`s are found that match `binding` then we check whether `binding`
        # is an alias of some other `Binding`. When it is we then re-run `doc` with that
        # `Binding`, otherwise if it's not an alias then we generate a summary for the
        # `binding` and display that to the user instead.
        alias = aliasof(binding)
        alias == binding ? summarize(alias, sig) : doc(alias, sig)
    else
        # There was at least one match for `binding` while searching. If there weren't any
        # matches for `sig` then we concatenate *all* the docs from the matching `Binding`s.
        if isempty(results)
            for group in groups, each in group.order
                push!(results, group.docs[each])
            end
        end
        # Get parsed docs and concatenate them.
        md = catdoc(map(parsedoc, results)...)
        # Save metadata in the generated markdown.
        if isa(md, Markdown.MD)
            md.meta[:results] = results
            md.meta[:binding] = binding
            md.meta[:typesig] = sig
        end
        return md
    end
end

# Some additional convenience `doc` methods that take objects rather than `Binding`s.
doc(obj::UnionAll) = doc(Base.unwrap_unionall(obj))
doc(object, sig::Type = Union{}) = doc(aliasof(object, typeof(object)), sig)
doc(object, sig...)              = doc(object, Tuple{sig...})

"""
    Docs.fielddoc(binding, field)

Returns documentation for a particular `field` of a type if it exists.
"""
function fielddoc(binding::Binding, field::Symbol)
    for mod in modules
        dict = meta(mod)
        if haskey(dict, binding)
            multidoc = dict[binding]
            if haskey(multidoc.docs, Union{})
                fields = multidoc.docs[Union{}].data[:fields]
                if haskey(fields, field)
                    doc = fields[field]
                    return isa(doc, Markdown.MD) ? doc : Markdown.parse(doc)
                end
            end
        end
    end
    fields = join(["`$f`" for f in fieldnames(resolve(binding))], ", ", ", and ")
    fields = isempty(fields) ? "no fields" : "fields $fields"
    Markdown.parse("`$(resolve(binding))` has $fields.")
end

# As with the additional `doc` methods, this converts an object to a `Binding` first.
fielddoc(object, field::Symbol) = fielddoc(aliasof(object, typeof(object)), field)

# Object Summaries.
# =================

function summarize(binding::Binding, sig)
    io = IOBuffer()
    println(io, "No documentation found.\n")
    if defined(binding)
        summarize(io, resolve(binding), binding)
    else
        println(io, "Binding `", binding, "` does not exist.")
    end
    md = Markdown.parse(seekstart(io))
    # Save metadata in the generated markdown.
    md.meta[:results] = DocStr[]
    md.meta[:binding] = binding
    md.meta[:typesig] = sig
    return md
end

function summarize(io::IO, λ::Function, binding)
    kind = startswith(string(binding.var), '@') ? "macro" : "`Function`"
    println(io, "`", binding, "` is a ", kind, ".")
    println(io, "```\n", methods(λ), "\n```")
end

function summarize(io::IO, T::DataType, binding)
    println(io, "# Summary")
    println(io, "```")
    println(io,
            T.abstract ? "abstract type" :
            T.mutable  ? "mutable struct" :
            Base.isstructtype(T) ? "struct" : "primitive type",
            " ", T, " <: ", supertype(T)
            )
    println(io, "```")
    if !T.abstract && T.name !== Tuple.name && !isempty(fieldnames(T))
        println(io, "# Fields")
        println(io, "```")
        pad = maximum(length(string(f)) for f in fieldnames(T))
        for (f, t) in zip(fieldnames(T), T.types)
            println(io, rpad(f, pad), " :: ", t)
        end
        println(io, "```")
    end
    if !isempty(subtypes(T))
        println(io, "# Subtypes")
        println(io, "```")
        for t in subtypes(T)
            println(io, t)
        end
        println(io, "```")
    end
    if supertype(T) != Any
        println(io, "# Supertype Hierarchy")
        println(io, "```")
        Base.show_supertypes(io, T)
        println(io)
        println(io, "```")
    end
end

function summarize(io::IO, m::Module, binding)
    readme = Pkg.dir(string(m), "README.md")
    if isfile(readme)
        println(io, "Displaying the `README.md` for the module instead.\n")
        println(io, "---\n")
        println(io, read(readme, String))
    else
        println(io, "No docstring or `README.md` found for module `", m, "`.\n")
    end
end

function summarize(io::IO, ::T, binding) where T
    println(io, "`", binding, "` is of type `", T, "`.\n")
    summarize(io, T, binding)
end

# Utilities.
# ==========

"""
`catdoc(xs...)`: Combine the documentation metadata `xs` into a single meta object.
"""
catdoc() = nothing
catdoc(xs...) = vcat(xs...)

const keywords = Dict{Symbol, DocStr}()

isdoc(s::AbstractString) = true

isdoc(x) = isexpr(x, :string) ||
    (isexpr(x, :macrocall) && x.args[1] === Symbol("@doc_str")) ||
    (isexpr(x, :call) && x.args[1] === Base.Markdown.doc_str)

function unblock(ex)
    isexpr(ex, :block) || return ex
    exs = filter(ex -> !(isa(ex, LineNumberNode) || isexpr(ex, :line)), ex.args)
    length(exs) == 1 || return ex
    return unblock(exs[1])
end

uncurly(ex) = isexpr(ex, :curly) ? ex.args[1] : ex

namify(x) = nameof(x, isexpr(x, :macro))

function nameof(x::Expr, ismacro)
    if isexpr(x, :.)
        ismacro ? macroname(x) : x
    else
        n = isexpr(x, (:module, :struct)) ? 2 : 1
        nameof(x.args[n], ismacro)
    end
end
nameof(q::QuoteNode, ismacro) = nameof(q.value, ismacro)
nameof(s::Symbol, ismacro)    = ismacro ? macroname(s) : s
nameof(other, ismacro)        = other

macroname(s::Symbol) = Symbol('@', s)
macroname(x::Expr)   = Expr(x.head, x.args[1], macroname(x.args[end].value))

isfield(x) = isexpr(x, :.) &&
    (isa(x.args[1], Symbol) || isfield(x.args[1])) &&
    (isa(x.args[2], QuoteNode) || isexpr(x.args[2], :quote))

# @doc expression builders.
# =========================

"""
    Docs.metadata(source, module, expr, ismodule)

Build a `Dict` expression containing metadata captured from the expression `expr`.

Fields that may be included in the returned `Dict`:

- `:path`:       Symbol representing the file where `expr` is defined.
- `:linenumber`: Linenumber where `expr` is defined.
- `:module`:     Module where the docstring is defined.
- `:fields`:     `Dict` of all field docs found in `expr`. Only for concrete types.
"""
function metadata(__source__, __module__, expr, ismodule)
    args = []
    # Filename and linenumber of the docstring.
    __file__ = isa(__source__.file, Symbol) ? String(__source__.file) : ""
    push!(args, Pair(:path, __file__))
    push!(args, Pair(:linenumber, __source__.line))
    # Module in which the docstring is defined.
    if ismodule # Module docs go inside the module with name `expr`
        push!(args, :($Pair(:module, $expr)))
    else
        push!(args, Pair(:module, __module__))
    end
    if isexpr(expr, :struct)
        # Field docs for concrete types.
        fields = []
        tmp = nothing
        for each in expr.args[3].args
            if isdoc(each)
                tmp = each
            elseif tmp !== nothing && (isa(each, Symbol) || isexpr(each, :(::)))
                push!(fields, (namify(each), tmp))
                tmp = nothing
            end
        end
        dict = :($(Dict)($([:($(Pair)($(quot(f)), $d)) for (f, d) in fields]...)))
        push!(args, :($(Pair)(:fields, $dict)))
    end
    :($(Dict)($(args...)))
end

function keyworddoc(__source__, __module__, str, def)
    docstr = esc(docexpr(__source__, __module__, lazy_iterpolate(str), metadata(__source__, __module__, def, false)))
    return :($(keywords)[$(esc(quot(def.name)))] = $docstr)
end

function objectdoc(__source__, __module__, str, def, expr, sig = :(Union{}))
    binding = esc(bindingexpr(namify(expr)))
    docstr  = esc(docexpr(__source__, __module__, lazy_iterpolate(str), metadata(__source__, __module__, expr, false)))
    quote
        $(esc(def))
        $(doc!)($__module__, $binding, $docstr, $(esc(sig)))
    end
end

function calldoc(__source__, __module__, str, def)
    args = def.args[2:end]
    if isempty(args) || all(validcall, args)
        objectdoc(__source__, __module__, str, nothing, def, signature(def))
    else
        docerror(def)
    end
end
validcall(x) = isa(x, Symbol) || isexpr(x, (:(::), :..., :kw, :parameters))

function moduledoc(__source__, __module__, meta, def, def′)
    name  = namify(def′)
    docex = Expr(:call, doc!, name, bindingexpr(name),
        docexpr(__source__, name, lazy_iterpolate(meta), metadata(__source__, __module__, name, true)))
    if def === nothing
        esc(:($eval($name, $(quot(docex)))))
    else
        def = unblock(def)
        block = def.args[3].args
        if !def.args[1]
            isempty(block) && error("empty baremodules are not documentable.")
            insert!(block, 2, :(import Base: @doc))
        end
        push!(block, docex)
        esc(Expr(:toplevel, def))
    end
end

# Shares a single doc, `meta`, between several expressions from the tuple expression `ex`.
function multidoc(__source__, __module__, meta, ex, define)
    out = Expr(:toplevel)
    str = docexpr(__source__, __module__, lazy_iterpolate(meta), metadata(__source__, __module__, ex, false))
    ref = Ref{DocStr}()
    for (n, arg) in enumerate(ex.args)
        # The first `arg` to be documented needs to also create the docstring for the group.
        # Subsequent `arg`s just need `ref` to be able to find the docstring without having
        # to create an entirely new one each.
        docstr = n === 1 ? :($(ref)[] = $str) : :($(ref)[])
        push!(out.args, docm(__source__, __module__, docstr, arg, define))
    end
    return out
end

"""
    @__doc__(ex)

Low-level macro used to mark expressions returned by a macro that should be documented. If
more than one expression is marked then the same docstring is applied to each expression.

    macro example(f)
        quote
            \$(f)() = 0
            @__doc__ \$(f)(x) = 1
            \$(f)(x, y) = 2
        end |> esc
    end

`@__doc__` has no effect when a macro that uses it is not documented.
"""
:(Core.@__doc__)

function __doc__!(meta, def, define)
    # Two cases must be handled here to avoid redefining all definitions contained in `def`:
    if define
        # `def` has not been defined yet (this is the common case, i.e. when not generating
        # the Base image). We just need to convert each `@__doc__` marker to an `@doc`.
        finddoc(def) do each
            each.head = :macrocall
            each.args = [Symbol("@doc"), nothing, meta, each.args[end], define] # TODO: forward line number info
        end
    else
        # `def` has already been defined during Base image gen so we just need to find and
        # document any subexpressions marked with `@__doc__`.
        docs  = []
        found = finddoc(def) do each
            push!(docs, :(@doc($meta, $(each.args[end]), $define)))
        end
        # If any subexpressions have been documented then replace the entire expression with
        # just those documented subexpressions to avoid redefining any definitions.
        if found
            def.head = :toplevel
            def.args = docs
        end
        found
    end
end
# Walk expression tree `def` and call `λ` when any `@__doc__` markers are found. Returns
# `true` to signify that at least one `@__doc__` has been found, and `false` otherwise.
function finddoc(λ, def::Expr)
    if isexpr(def, :block, 2) && isexpr(def.args[1], :meta, 1) && def.args[1].args[1] === :doc
        # Found the macroexpansion of an `@__doc__` expression.
        λ(def)
        true
    else
        found = false
        for each in def.args
            found |= finddoc(λ, each)
        end
        found
    end
end
finddoc(λ, def) = false

# Predicates and helpers for `docm` expression selection:

const FUNC_HEADS    = [:function, :stagedfunction, :macro, :(=)]
const BINDING_HEADS = [:typealias, :const, :global, :(=)]  # deprecation: remove `typealias` post-0.6
# For the special `:@mac` / `:(Base.@mac)` syntax for documenting a macro after definition.
isquotedmacrocall(x) =
    isexpr(x, :copyast, 1) &&
    isa(x.args[1], QuoteNode) &&
    isexpr(x.args[1].value, :macrocall, 2)
# Simple expressions / atoms the may be documented.
isbasicdoc(x) = isexpr(x, :.) || isa(x, Union{QuoteNode, Symbol})
is_signature(x) = isexpr(x, :call) || (isexpr(x, :(::), 2) && isexpr(x.args[1], :call)) || isexpr(x, :where)

function docm(source::LineNumberNode, mod::Module, meta, ex, define = true)
    # Some documented expressions may be decorated with macro calls which obscure the actual
    # expression. Expand the macro calls and remove extra blocks.
    x = unblock(macroexpand(mod, ex))
    # Don't try to redefine expressions. This is only needed for `Base` img gen since
    # otherwise calling `loaddocs` would redefine all documented functions and types.
    def = define ? x : nothing
    if isa(x, GlobalRef) && (x::GlobalRef).mod == mod
        x = (x::GlobalRef).name
    end

    # Keywords using the `@kw_str` macro in `base/docs/basedocs.jl`.
    #
    #   "..."
    #   kw"if", kw"else"
    #
    isa(x, Base.BaseDocs.Keyword) ? keyworddoc(source, mod, meta, x) :

    # Method / macro definitions and "call" syntax.
    #
    #   function f(...) ... end
    #   f(...) = ...
    #   macro m(...) end
    #   function f end
    #   f(...)
    #
    isexpr(x, FUNC_HEADS) && is_signature(x.args[1])   ? objectdoc(source, mod, meta, def, x, signature(x)) :
    isexpr(x, :function)  && !isexpr(x.args[1], :call) ? objectdoc(source, mod, meta, def, x) :
    isexpr(x, :call)                                   ? calldoc(source, mod, meta, x) :

    # Type definitions.
    #
    #   struct T ... end
    #   abstract type T end
    #   primitive type T N end
    #
    isexpr(x, [:struct, :abstract, :primitive]) ? objectdoc(source, mod, meta, def, x) :

    # "Bindings". Names that resolve to objects with different names, ie.
    #
    #   const T = S
    #   T = S
    #   global T = S
    #
    isexpr(x, BINDING_HEADS) && !isexpr(x.args[1], :call) ? objectdoc(source, mod, meta, def, x) :

    # Quoted macrocall syntax. `:@time` / `:(Base.@time)`.
    isquotedmacrocall(x) ? objectdoc(source, mod, meta, def, x) :
    # Modules and baremodules.
    isexpr(x, :module) ? moduledoc(source, mod, meta, def, x) :
    # Document several expressions with the same docstring. `a, b, c`.
    isexpr(x, :tuple) ? multidoc(source, mod, meta, x, define) :
    # Errors generated by calling `macroexpand` are passed back to the call site.
    isexpr(x, :error) ? esc(x) :
    # When documenting macro-generated code we look for embedded `@__doc__` calls.
    __doc__!(meta, x, define) ? esc(x) :
    # Any "basic" expression such as a bare function or module name or numeric literal.
    isbasicdoc(x) ? objectdoc(source, mod, meta, nothing, x) :

    # All other expressions are undocumentable and should be handled on a case-by-case basis
    # with `@__doc__`. Unbound string literals are also undocumentable since they cannot be
    # retrieved from the module's metadata `ObjectIdDict` without a reference to the string.
    docerror(ex)
end

function docerror(ex)
    txt = """
    cannot document the following expression:

    $(isa(ex, AbstractString) ? repr(ex) : ex)"""
    if isexpr(ex, :macrocall)
        txt *= "\n\n'$(ex.args[1])' not documentable. See 'Base.@__doc__' docs for details."
    end
    :($(error)($txt, "\n"))
end

function docm(source::LineNumberNode, mod::Module, ex)
    if isexpr(ex, :->)
        docm(source, mod, ex.args...)
    elseif haskey(keywords, ex)
        parsedoc(keywords[ex])
    elseif isa(ex, Union{Expr, Symbol})
        binding = esc(bindingexpr(namify(ex)))
        if isexpr(ex, :call) || isexpr(ex, :macrocall)
            sig = esc(signature(ex))
            :($(doc)($binding, $sig))
        else
            :($(doc)($binding))
        end
    else
        :($(doc)($(typeof)($(esc(ex)))))
    end
end

# MD support
catdoc(md::MD...) = MD(md...)

include("utils.jl")

# Swap out the bootstrap macro with the real one.
Core.atdoc!(docm)

macro local_hygiene(expr)
    # removes `esc` Exprs relative to the module argument to expand
    # and resolves everything else relative to this (Doc) module
    # this allows us to get good errors and backtraces
    # from calling docm (by not using macros),
    # while also getting macro-expansion correct (by using the macro-expander)
    return expr
end
function loaddocs(docs)
    unescape = GlobalRef(Docs, Symbol("@local_hygiene"))
    for (mod, ex, str, file, line) in docs
        data = Dict(:path => string(file), :linenumber => line)
        doc = docstr(str, data)
        docstring = docm(LineNumberNode(line, file), mod, doc, ex, false) # expand the real @doc macro now
        eval(mod, Expr(:macrocall, unescape, nothing, docstring))
    end
    empty!(docs)
end

end

max_stars_repo_name: Mikewl/julia

max_issues_repo_name: Mikewl/julia

max_forks_repo_name: Mikewl/julia


text:
abstract AbstractPerm{T} <: AbstractMatrix{T}

export AbstractPerm
export plength, isid, topiv, numcycles, iscyclic

import Base: rank, sign, det, logdet, trace, ishermitian, issym,
iseven, istriu, istril, isposdef, null, getindex, transpose,
ctranspose, inv, pmap, isperm, one, zero, full, sparse, size, eltype

import Base: eig, eigfact, eigmax, eigmin, eigs, eigvals, eigvecs
import DataStructures: list

size(m::AbstractPerm) = (s = plength(m); (s,s))
eltype{T}(c::AbstractPerm{T}) = T

# this one is not called for some reason
#getindex{T}(m::AbstractPerm{T}, i::Real, j::Real) =  pmap(m,j) == i ? one(T) : zero(T)
# this one is called
getindex(m::AbstractPerm, i::Real, j::Real) =  pmap(m,i) == j ? one(eltype(m)) : zero(eltype(m))
\(p::AbstractPerm, k::Integer) = k / p

function mkerrf()
    for sym in (:plength, :isid, :isperm, :pmap, :sign, :list)
        @eval begin
            ($sym)(p::AbstractPerm) = error("AbstractPerm: `", $sym , "' not defined for ", typeof(p))
        end
    end
end

mkerrf()

# dense matrix
full(p::AbstractPerm) = (T=eltype(p); T[p[i,j] for i=1:plength(p), j=1:plength(p)])

# sparse matrix
function sparse{T}(m::AbstractPerm{T})
    n = plength(m)
    I = [1:n]
    J = T[ pmap(m,i) for i in 1:n] # should call a routine for this, so that PermCyc is efficient
#    J = (m.data)[I]
    V = ones(T,n)
    sparse(I,J,V)
end

zero(p::AbstractPerm) = error("zero not defined for type ", typeof(p))
ctranspose(m::AbstractPerm) = inv(m)
transpose(m::AbstractPerm) = inv(m)
inv(m::AbstractPerm) = m^-1
topiv(p::AbstractPerm) = PermPlain.perm2ipiv(list(p).data)
numcycles(p::AbstractPerm) = length(cyclelengths(p))
# We allow fixed points in this definition. See wikipedia
iscyclic(p::AbstractPerm) = numcycles(p) == 1

rank(p::AbstractPerm) = plength(p)
det(p::AbstractPerm) = sign(p)
logdet(p::AbstractPerm) = sign(p) > 0 ? 0 : error("AbstractPerm: DomainError: determinant is -1")
trace(p::AbstractPerm) = plength(p) - supportsize(p)
ishermitian(p::AbstractPerm) = isid(p)
issym(p::AbstractPerm) = isid(p)
iseven(p::AbstractPerm) = sign(p) == one(eltype(p))
istriu(p::AbstractPerm) = isid(p)
istril(p::AbstractPerm) = isid(p)
isposdef(p::AbstractPerm) = isid(p)
null(p::AbstractPerm) = zeros(Float64,plength(p),0) # for consistency

# These can be computed more efficiently without using full matrices
# lufact can be done very easily, but we have to break into the
# interface
for f in ( :eig, :eigfact, :eigmax, :eigmin, :eigs, :eigvals, :eigvecs)
    @eval begin
        ($f)(p::AbstractPerm) = ($f)(full(p))
    end
end

ext: jl


text:
todisplay = VERSION <= VersionNumber("1.2") ? true : false

using CAGD
using SparseArrays
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

if todisplay
	using ViewerGL
	GL = ViewerGL
	include("../examples/3d/views.jl")
end

V1, (VV,EV,FV,CV) = Lar.simplex(3, true)
tetra1 = V1, EV, FV, CV
V2 = [
    0.0 0.5 0.1 0.1
    0.0 0.1 0.5 0.1
    0.0 0.1 0.1 0.5
]
tetra2 = V2, EV, FV, CV
twotetra = Lar.Struct([ tetra1, tetra2 ])
V,EV,FV,CV = Lar.struct2lar(twotetra)
model = CAGD.Model(V);
CAGD.addModelCells!(model, 1, convert(Lar.ChainOp, Lar.coboundary_0(EV::Lar.Cells)))
CAGD.addModelCells!(model, 2, Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells))

# Spatial Arrangement
split_model = CAGD.facesplitting(model)

if todisplay  displayModel(split_model)  end

congr_model = CAGD.mergeModelVertices(split_model, signed_merge=true)

if todisplay  displayModel(congr_model)  end

gift_model = deepcopy(congr_model)
FCcycles, bicon_comps = CAGD.tgw(congr_model, 3)
CAGD.addModelCells!(gift_model, 3, convert(Lar.ChainOp, FCcycles'))

if todisplay  viewExplode(gift_model)  end

externals = [bc[CAGD.getExternalBoundary(congr_model, 3, FCcycles, bc)] for bc in bicon_comps]
internals = [setdiff(bicon_comps[i], externals[i]) for i = 1 : length(externals)]
assembled_model = deepcopy(congr_model)
## To Fix
FC = FCcycles[:, internals...]
JK = ones(Int8, length(externals))
FC = [FCcycles * SparseArrays.sparse(externals, JK, JK, size(FCcycles, 2), 1) FC]
CAGD.addModelCells!(assembled_model, 3, convert(Lar.ChainOp, FC'))

if todisplay  viewExplode(assembled_model)  end


text:
module CommunalHelperCassetteFallingBlock

using ..Ahorn, Maple

@mapdef Entity "CommunalHelper/CassetteFallingBlock" CassetteFallingBlock(x::Integer,
                                                                  y::Integer,
                                                                  width::Integer=Maple.defaultBlockWidth,
                                                                  height::Integer=Maple.defaultBlockHeight,
                                                                  index::Integer=0,
                                                                  tempo::Number=1.0)

const colorNames = Dict{String, Int}(
    "Blue" => 0,
    "Rose" => 1,
    "Bright Sun" => 2,
    "Malachite" => 3
)

const colors = Dict{Int, Ahorn.colorTupleType}(
    1 => (240, 73, 190, 255) ./ 255,
	2 => (252, 220, 58, 255) ./ 255,
	3 => (56, 224, 78, 255) ./ 255
)
const defaultColor = (73, 170, 240, 255) ./ 255

const placements = Ahorn.PlacementDict(
    "Cassette Falling Block ($index - $color) (Communal Helper)" => Ahorn.EntityPlacement(
        CassetteFallingBlock,
        "rectangle",
        Dict{String, Any}(
            "index" => index,
        )
    ) for (color, index) in colorNames
)

Ahorn.editingOptions(entity::CassetteFallingBlock) = Dict{String, Any}(
    "index" => colorNames
)

Ahorn.nodeLimits(entity::CassetteFallingBlock) = 1, 1

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

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

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

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

const block = "objects/cassetteblock/solid"

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::CassetteFallingBlock, room::Maple.Room)
    x, y = Ahorn.position(entity)

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

    tileWidth = ceil(Int, width / 8)
    tileHeight = ceil(Int, height / 8)

    index = Int(get(entity.data, "index", 0))
    color = get(colors, index, defaultColor)

    for i in 1:tileWidth, j in 1:tileHeight
        tx = (i == 1) ? 0 : ((i == tileWidth) ? 16 : 8)
        ty = (j == 1) ? 0 : ((j == tileHeight) ? 16 : 8)

        Ahorn.drawImage(ctx, block, x + (i - 1) * 8, y + (j - 1) * 8, tx, ty, 8, 8, tint=color)
    end
end

end

max_stars_repo_name: FlynnSC/CommunalHelper

max_issues_repo_name: FlynnSC/CommunalHelper

max_forks_repo_name: FlynnSC/CommunalHelper


text:
using LinearAlgebra
using Arpack
using LinearMaps

mutable struct Data{T}
    """
    Data structure for the solution of
        min         ½x'Px + q'x
        subject to  Ax ≤ b
                    ‖x‖ = r
    with an active set algorithm.

    The most important element is F, which holds
    - F.QR:  an updatable QR factorization of the working constraints
    - F.Z:   a view on an orthonormal matrix spanning the nullspace of the working constraints
    - F.P:   the hessian of the problem
    - F.ZPZ: equal to F.Z'*F.P*F.P, i.e. the reduced hessian

    The Data structure also keeps matrices of the constraints not in the working set (A_ignored and b_ignored)
    stored in a continuous manner. This allows efficient use of BLAS.
    """
    x::Vector{T}  # Variable of the minimization problem
    Xmin::Matrix{T}
    y::Vector{T}  # Reduced variable on the nullspace of the working constraints
    x0::Vector{T} # x - F.Z'*y
    f0::T # Constant term in reduced TRS

    n::Int  # length(x)
    m::Int  # Number of constraints

    F::NullspaceHessian{T}  # See first comments on the definition of Data
    q::Vector{T}
    Zq::Vector{T} # Reduced linear cost
    A::Matrix{T}
    b::Vector{T}
    r::T
    working_set::Vector{Int}
    ignored_set::Vector{Int}

    done::Bool
    removal_idx::Int

    A_ignored::SubArray{T, 2, Matrix{T}, Tuple{UnitRange{Int}, Base.Slice{Base.OneTo{Int}}}, false}
    b_ignored::SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int}}, true}
    A_shuffled::Matrix{T}  # That's where A_ignored "views into"
    b_shuffled::Vector{T}  # That's where b_ignored "views into"

    # Options
    tolerance::T
    verbosity::Int
    printing_interval::Int

    # Logging
    iteration::Int
    λ::Vector{T}
    μ::T
    residual::T
    trs_choice::Char
    grad_steps::Int
    eigen_steps::Int
    timings::DataFrame

    function Data(P::Matrix{T}, q::Vector{T}, A::Matrix{T}, b::Vector{T},
        r::T, x::Vector{T}; kwargs...) where T

        m, n = size(A)
        working_set = findall((A*x - b)[:] .>= 1)
        if length(working_set) >= n
            working_set = working_set[1:n-1]
        end
        ignored_set = setdiff(1:m, working_set)

        QR = UpdatableQR(A[working_set, :]')
        A_shuffled = zeros(T, m, n)
        l = length(ignored_set)
        A_shuffled[end-l+1:end, :] .= view(A, ignored_set, :)
        b_shuffled = zeros(T, m)
        b_shuffled[end-l+1:end] .= view(b, ignored_set)

        Data(P, q, A, b, r, x,
            QR, working_set, ignored_set,
            A_shuffled, b_shuffled; kwargs...)
    end

    function Data(P::Matrix{T}, q::Vector{T},  A::Matrix{T}, b::Vector{T}, r::T, x::Vector{T},
        QR::UpdatableQR{T}, working_set::Vector{Int}, ignored_set::Vector{Int},
        A_shuffled::Matrix{T}, b_shuffled::Vector{T};
        verbosity=1, printing_interval=50, tolerance=1e-11) where T

        F = NullspaceHessian{T}(P, QR)
        l = length(ignored_set)

        m, n = size(A)
        λ = zeros(T, m);

        new{T}(x, zeros(T, 0, 0), zeros(T, 0), zeros(T, 0), 0.0,
            n, m, F, q, zeros(T, 0), A, b, r, working_set, ignored_set, false, 0,
            view(A_shuffled, m-l+1:m, :),
            view(b_shuffled, m-l+1:m),
            A_shuffled, b_shuffled,
            T(tolerance), verbosity, printing_interval,  # Options
            0, λ, 0, NaN, '-', 0, 0,  # Logging
            # Timings
            DataFrame(projection=zeros(T, 0), gradient_steps=zeros(T, 0),
            trs=zeros(T, 0), trs_local=zeros(T, 0),
            move_to_optimizer=zeros(T, 0), move_to_local_optimizer=zeros(T, 0),
            curvature_step=zeros(T, 0), kkt=zeros(T, 0),
            add_constraint=zeros(T, 0), remove_constraint=zeros(T, 0))
        )
    end
end

function solve_boundary(P::Matrix{Float64}, q::Vector{Float64}, A::Matrix{Float64}, b::Vector{Float64}, r::Float64,
    x::Vector{Float64}; max_iter=Inf, kwargs...)
    data = Data(P, q, A, b, r, x; kwargs...)

    if data.verbosity > 0
        print_header(data)
        print_info(data)
    end

    while !data.done && data.iteration <= max_iter
        iterate!(data)

        if data.verbosity > 0
            mod(data.iteration, 10*data.printing_interval) == 0 && print_header(data)
            (mod(data.iteration, data.printing_interval) == 0 || data.done) && print_info(data)
        end
    end
    return data.x, [data.λ; data.μ]
end

function iterate!(data::Data{T}) where{T}
    t_remove, t_add, t_proj, t_grad, t_kkt = zero(T), zero(T), zero(T), zero(T), zero(T)
    t_trs, t_move, t_trs_l, t_move_l, t_curv = zero(T), zero(T), zero(T), zero(T), zero(T)

    if data.removal_idx > 0
        t_remove = @elapsed remove_constraint!(data, data.removal_idx)
        data.removal_idx = 0
    end

    t_proj = @elapsed data.y = data.F.Z'*data.x  # Reduced free variable
    if norm(data.y) <= 1e-10
        @warn "LICQ failed. Terminating"
        data.done = true
        return
    end
    data.x0 = data.x - data.F.Z*data.y
    data.Zq = data.F.Z'*(data.q + data.F.P*data.x0) # Reduced q
    data.f0 = dot(data.x0, (data.F.P*data.x0))/2 + dot(data.q, data.x0)

    t_grad = @elapsed new_constraint = gradient_steps(data, 5)
    if isnan(new_constraint) && norm(projected_gradient(data, data.y)) < 1e-9 && length(data.Xmin) > 0
        Ymin = data.F.Z'*data.Xmin # Minimizer of the previous working set
        t_curv += @elapsed new_constraint = curvature_step(data, Ymin[:, 1], false)
    end

    if isnan(new_constraint)
        t_trs = @elapsed Ymin, info = trs_robust(data.F.ZPZ, data.Zq, norm(data.y), tol=1e-11)
        data.Xmin = data.x0 .+ data.F.Z*Ymin;
        t_move = @elapsed new_constraint, is_minimizer = move_towards_optimizers!(data, Ymin)
        if isnan(new_constraint) && !is_minimizer
            t_trs_l = @elapsed Ymin, info = trs_robust(data.F.ZPZ, data.Zq, norm(data.y), tol=1e-11, compute_local=true)
            t_move_l = @elapsed new_constraint, is_minimizer = move_towards_optimizers!(data, Ymin)
            if isnan(new_constraint) && !is_minimizer
                t_grad += @elapsed new_constraint = gradient_steps(data)
                if isnan(new_constraint)
                    t_curv += @elapsed new_constraint = curvature_step(data, Ymin[:, 1])
                end
            end
        end
    end
    data.x = data.x0 + data.F.Z*data.y;
    if !isnan(new_constraint)
        t_add = @elapsed add_constraint!(data, new_constraint)
    end
    if isnan(new_constraint) || data.F.m <= 1
        t_kkt = @elapsed data.removal_idx = check_kkt!(data)
    end

    if data.done && data.verbosity >= 2
        push!(data.timings, [t_proj, t_grad, t_trs, t_trs_l, t_move, t_move_l, t_curv, t_kkt, t_add, t_remove])
        data.timings |> CSV.write(string("timings.csv"))
        sums =  [sum(data.timings[i]) for i in 1 : size(data.timings, 2)]
        labels =  names(data.timings)
        show(stdout, "text/plain", [labels sums]); println()
        @printf "Total time (excluding factorizations): %.4e seconds.\n" sum(sums[1:end-2])
    end
    data.iteration += 1
end

function check_kkt!(data)
    F = qr(data.F.Z'*data.x)
    R = view(data.F.QR.R, 1:data.F.QR.m+1, 1:data.F.QR.m+1)
    R[1:end-1, end] = data.F.QR.Q1'*data.x
    R[end, end] = F.factors[1, 1]

    g = data.F.P*data.x + data.q
    Qg = -[data.F.QR.Q1'*g; (F.Q'*(data.F.QR.Q2'*g))[1]]
    multipliers = UpperTriangular(R)\Qg

    λ = multipliers[1:end-1]
    μ = multipliers[end]
    data.residual = norm(projected_gradient(data, data.y))
    # data.residual = norm(g + [data.A[data.working_set, :]' data.x]*[λ; μ])

    data.λ .= 0
    data.λ[data.working_set] .= λ
    data.μ = μ
    if all(data.λ .>= -max(1e-8, data.residual))
        data.done = true
        data.removal_idx = 0
    else
        data.removal_idx = argmin(λ)
    end

    return data.removal_idx
end

function move_towards_optimizers!(data, Y, flag=true)
    how_many = size(Y, 2)
    if isfeasible(data, Y[:, 1])
        data.y = Y[:, 1]
        data.trs_choice = 'g'
        return NaN, true
    end

    if how_many == 1
        data.trs_choice = 'm'
        d1 = data.y
        d2 = data.y - Y[:, 1];
        minimizer_value = f(data, Y[:, 1])
    else
        data.trs_choice = 'm'
        d1 = data.y - Y[:, 1]
        d2 = data.y - Y[:, 2];
        minimizer_value = max(f(data, Y[:, 1]), f(data, Y[:, 2]))
    end
    D = qr([d1 d2]).Q
    new_constraint = minimize_2d(data, D[:, 1], D[:, 2])
    if isnan(new_constraint)
        is_minimizer = ((norm(projected_gradient(data, data.y)) <= 1e-6)
                        && (f(data, data.y) <= minimizer_value + max(abs(minimizer_value)/100, 1e-7)))
    else
        is_minimizer = false
    end

    return new_constraint, is_minimizer
end

function gradient_steps(data, max_iter=1e4)
    k = 0;
    new_constraint = NaN
    project_grad = projected_gradient(data, data.y)
    while isnan(new_constraint) && k < max_iter && norm(project_grad)/max(abs(f(data, data.y)), 10) >= 1e-8
        d1 = data.y/norm(data.y)
        d2 = project_grad/norm(project_grad)
        new_constraint = minimize_2d(data, d1, d2)
        project_grad = projected_gradient(data, data.y)
        k += 1
        if mod(k, 1000) == 0
            @warn string(k, "th gradient step with f: ", f(data, data.y), " proj grad norm: ", norm(project_grad))
        end
    end
    data.grad_steps += k

    return new_constraint
end

function curvature_step(data, y_global, warning=true)
    # Find tangent eigenvector v_min which should correspond to a negative eigenvalue
    project!(x) = axpy!(-dot(x, data.y)/dot(data.y, data.y), data.y, x)
    l = -(data.y'*(data.F.ZPZ*data.y)+data.Zq'*data.y)/dot(data.y, data.y) - 100 # Approximate Lagrange multiplier
    function custom_mul!(y::AbstractVector, x::AbstractVector)
        # print(dot(x, data.y), " ")
        project!(x)
        mul!(y, data.F.ZPZ, x)
        axpy!(l, x, y)
        project!(y)
    end
    L = LinearMap{Float64}(custom_mul!, data.F.m; ismutating=true, issymmetric=true)
    λ_min, v_min, _ = eigs(L, nev=1, which=:SR, v0=project!(randn(data.F.m)), tol=1e-5)
    λ_min = λ_min[1]; v_min = v_min[:, 1]
    λ_min += 100

    if abs(dot(v_min, data.y)) >= 1e-6
        warning && @warn "Search for negative curvature failed; we assume that we are in a local minimum"
        return NaN
    end
    if λ_min >= -1e-8
        warning && @warn "We ended up in an unexpected TRS local minimum. Perhaps the problem is badly scaled." λ_min
        return NaN
    end
    d1 = project!(v_min)
    d2 = y_global - data.y
    D = qr([d1 d2]).Q

    y_copy = copy(data.y)
    new_constraint = minimize_2d(data, D[:, 1], D[:, 2], false)
    if isnan(new_constraint)
        data.y = y_copy
        new_constraint = minimize_2d(data, D[:, 1], -D[:, 2], false)
        warning && isnan(new_constraint) && @warn "Could not activate a new constraint; checking KKT conditions on a saddle point."
    end
    return new_constraint
end

function find_violations!(violating_constraints::Vector{Bool}, a1::Vector{T}, a2::Vector{T}, b::Vector{T}, x::T, y::T, tol::T) where{T}
    n = length(a1)
    @inbounds for i = 1:n
        violating_constraints[i] = x*a1[i] + y*a2[i] - b[i] >= tol
    end
    return violating_constraints
end

function _minimize_2d(P::Matrix{T}, q::Vector{T}, r::T,
    a1::Vector{T}, a2::Vector{T}, b::Vector{T}, x0::T, y0::T,
    allow_flip=true) where {T}

    grad = P*[x0; y0] + q
    flip = false
    if angle(x0, y0, grad[1], grad[2]) < pi && allow_flip
        flip = true
        # Flip y axis, thus obtaining dot(grad, [x0, y0]) > 0
        y0 = -y0
        a2 = -a2
        P[1, 2] = -P[1, 2]; P[2, 1] = -P[2, 1]
        q[2] = -q[2]
    end

    infeasibility(x, y) = maximum(a1*x + a2*y - b)
    tol = max(1e-12, 1.05*infeasibility(x0, y0))

    new_constraint = NaN
    X, info = trs_boundary_small(P, q, r; compute_local=true)

    if infeasibility(X[1, 1], X[2, 1]) <= tol
        x, y = X[1, 1], X[2, 1]
    else
        θ0 = angle(one(T), zero(T), x0, y0) # atan(y0, x0)
        δθ_global = angle(x0, y0, X[1, 1], X[2, 1])
        # δθ_global = 2*pi - 2e-7 < δθ_global ? zero(T) : δθ_global
        if size(X, 2) == 1
            δθ = δθ_global
        else
            δθ_local = angle(x0, y0, X[1, 2], X[2, 2])
            # δθ_local = 2*pi - 2e-7 < δθ_local ? zero(T) : δθ_local
            δθ = min(δθ_global, δθ_local)
        end
        θ = θ0 + δθ

        n = length(b)
        violating_constraints = Vector{Bool}(undef, n)

        if !any(find_violations!(violating_constraints, a1, a2, b, r*cos(θ), r*sin(θ), tol))
            x, y = r*cos(θ), r*sin(θ)
        else
            f_2d(x, y) = dot([x; y], P*[x; y])/2 + dot([x; y], q)
            θ_low = θ0; θ_high = θ
            # Binary Search
            for i = 1:80
                find_violations!(violating_constraints, a1, a2, b, r*cos(θ), r*sin(θ), tol)
                if any(violating_constraints) || f_2d(r*cos(θ), r*sin(θ)) > f_2d(x0, y0)
                    θ_high = θ
                    new_constraint = findlast(violating_constraints)
                else
                    θ_low = θ
                end
                θ = θ_low + (θ_high - θ_low)/2
            end

            if new_constraint == nothing || isnan(new_constraint)
                x, y = r*cos(θ), r*sin(θ)
                new_constraint = NaN
            else
                tol = 1.2*tol # Lossen up the tolerance more now that we identified the new constraint
                x1, y1, x2, y2 = circle_line_intersections(a1[new_constraint], a2[new_constraint], b[new_constraint], r)
                if infeasibility(x1, y1) <= tol
                    if infeasibility(x2, y2) <= tol && f_2d(x2, y2) <= f_2d(x1, y1)
                        x, y = x2, y2
                    else
                        x, y = x1, y1
                    end
                elseif infeasibility(x2, y2) <= tol
                    x, y = x2, y2
                else
                    x, y = r*cos(θ_low), r*sin(θ_low) # Just in case the circle_line_intersections fails
                    if abs(a1[new_constraint]*x + a2[new_constraint]*y - b[new_constraint]) >= tol
                        new_constraint = NaN
                    end
                end
            end
        end
    end
    if flip; y = -y; end

    return x, y, new_constraint
end

function minimize_2d(data, d1, d2, allow_flip=true)
    x0 = dot(data.y, d1); y0 = dot(data.y, d2)
    z = data.y - d1*x0 - d2*y0
    r = sqrt(x0^2 + y0^2)

    # Calculate 2d cost matrix [P11 P12
    #                           P12 P22]
    # and vector [q1; q2]
    Pd1 = Symmetric(data.F.ZPZ)*d1
    P11 = dot(d1, Pd1)
    q1 = dot(d1, data.Zq) + dot(Pd1, z)

    Pd2 = Symmetric(data.F.ZPZ)*d2
    P22 = dot(d2, Pd2)
    P12 = dot(d1, Pd2)
    q2 = dot(d2, data.Zq) + dot(Pd2, z)
    P = [P11 P12; P12 P22]; q = [q1; q2]

    b = data.b_ignored - data.A_ignored*(data.F.Z*z + data.x0)
    a1 = data.A_ignored*(data.F.Z*d1)
    a2 = data.A_ignored*(data.F.Z*d2)

    # Discard perfectly correlated constraints
    idx = a1.^2 + a2.^2 .<= 1e-18
    a1[idx] .= 1; a2[idx] .= 0; b[idx] .= 2*r

    x, y, new_constraint = _minimize_2d(P, q, r, a1, a2, b, x0, y0, allow_flip)
    @. data.y = z + x*d1 + y*d2

    return new_constraint
end

function find_closest_minimizer(y1, y2, y_g, y_l)
    theta_g = angle(y1, y2, y_g[1], y_g[2])
    if isempty(y_l)
        return theta_g
    else
        theta_l = angle(y1, y2, y_l[1], y_l[2])
        min(theta_l, theta_g)
    end
end

function angle(x1, y1, x2, y2)
    dot = x1*x2 + y1*y2      # dot product between [x1, y1] and [x2, y2]
    det = x1*y2 - y1*x2      # determinant
    angle = atan(det, dot)  # atan2(y, x) or atan2(sin, cos))
    if angle < 0
        angle += 2*pi
    end
    return angle
end

function isfeasible(data::Data{T}, x) where{T}
    tol = max(1.5*maximum(data.A*data.x - data.b), data.tolerance)
    return all(data.A_ignored*(data.F.Z*x + data.x0) - data.b_ignored .<= data.tolerance)
end

function f(data::Data{T}, x) where{T}
    return dot(x, data.F.ZPZ*x)/2 + dot(x, data.Zq) + data.f0
end

function grad(data::Data{T}, x) where{T}
    return data.F.ZPZ*x + data.Zq
end

function projected_gradient(data::Data{T}, x) where{T}
    if length(data.F.ZPZ) == 1 # This is because data.y is not always updated after adding a constraint
        return 0*x
    end

    gradient = grad(data, x)
    F = qr([x gradient])
    if minimum(size(F)) > 1
        v = F.Q[:, 2]
        return dot(v, gradient)*v
    else
        return 0*x
    end
end


text:
__precompile__()

module NeuralNetDiffEq
#dependencies

using Reexport
@reexport using DiffEqBase
using Knet, Compat, ForwardDiff
import DiffEqBase: solve

# Abstract Types
@compat abstract type NeuralNetDiffEqAlgorithm <: AbstractODEAlgorithm end

# DAE Algorithms
immutable nnode <: NeuralNetDiffEqAlgorithm
    hl_width::Int
end

nnode(;hl_width=10) = nnode(hl_width)

export nnode

#constant functions
#z(P,i,x) = P[1][i]*x + P[2][i]
sig_der(x) = sigm(x)*(1-sigm(x))

function solve(
    prob::AbstractODEProblem,
    alg::NeuralNetDiffEqAlgorithm;
    dt = nothing,
    timeseries_errors = true,
    iterations = 50)

    u0 = prob.u0
    tspan = prob.tspan
    f = prob.f
    t0 = tspan[1]

    if dt == nothing
        error("dt must be set.")
    end

    #types and dimensions
    uElType = eltype(u0)
    tType = typeof(tspan[1])
    outdim = length(u0)

    #hidden layer
    hl_width = alg.hl_width

    #The phi trial solution
    phi(P,t) = u0 + (t-t0)*predict(P,t)

    #train points generation
    dtrn = generate_data(tspan[1],tspan[2],dt,atype=tType)

    #iterations
    _maxiters = iterations

    #initialization of weights and bias
    w = init_params(uElType,hl_width)
    # println(w[1])
    # println(w[2])
    # println(w[3])

    #initialization of optimization parameters (Adam by default for now)
    lr_ = 0.1
    beta1_ = 0.9
    beta2_ = 0.95
    eps_ = 1e-6
    prms = Any[]

    for i=1:length(w)
    prm = Adam(lr=lr_, beta1=beta1_, beta2=beta2_, eps=eps_)
    push!(prms, prm)
    end

    #iters = 1000

    #reporting the accuracy
    #report(epoch)=println((:epoch,epoch,:trn,accuracy(w,dtrn,predict),:tst,accuracy(w,dtst,predict)))
    #report(epoch)=println((:epoch,epoch,:trn,loss_trial(w,dtrn)))
    #report(0)
    #P_tuned = train(w,prms,dtrn,regFlag; epochs=100, iters=1000)
    #@time for epoch=1:epochs
    #@time train(w, prms, dtrn, f, phi, hl_width; maxiters=_maxiters)
        #report(epoch)
    #end

    @time for iters=1:_maxiters
            train(w, prms, dtrn, f, phi, hl_width; maxiters=1)
            loss = test(w,dtrn,f,phi,hl_width)
            if mod(iters,100) == 0
                println((:iteration,iters,:loss,loss))
            end
            #gradcheck(loss_trial, w, dtrn, f, phi, hl_width...; gcheck=10, verbose=true)
            #check_Phi_gradient(w,dtrn,hl_width)
            if loss < 10^(-8.0)
                break
            end
        end

    # for t in log; println(t); end
    # return w

    #solutions at timepoints
    u = [phi(w,x) for x in dtrn]

    build_solution(prob,alg,dtrn,u,
               timeseries_errors = timeseries_errors,
               retcode = :Success)

end #solve

include("training_utils.jl")
#include("interface.jl")
end # module


text:
template = Dict([
    ("seed", 0),
    ("d_fitness", 1),
    ("n_population", 100),
    ("n_elite", 10),
    ("n_gen", 50000),
    ("log_gen", 100),
    ("save_gen", 1000),
    ("m_rate", 0.5),
    ("rows", 1),
    ("columns", 25),
    ("recur", 0.9),
    ("n_out", 1),
    ("n_in", 12),
    ("out_m_rate", 0.9),
    ("functions", [
        "f_subtract",
        "f_add",
        "f_div",
        "f_mult",
        "f_sqrt",
        "f_abs",
        "f_exp",
        "f_pow",
        "f_sin",
        "f_cos",
        "f_tanh",
        "f_sqrt_xy",
        "f_lt",
        "f_gt"
        ]
    )
])

max_stars_repo_name: mahmoud-al-najar/cgp-morphcast

max_issues_repo_name: mahmoud-al-najar/cgp-morphcast

max_forks_repo_name: mahmoud-al-najar/cgp-morphcast


text:
@testset "Evalcode" begin
    # see `converter/md_blocks:convert_code_block`
    # see `converter/lx/resolve_lx_input_*`
    # --------------------------------------------
    h = raw"""
        Simple code:
        ```julia:./code/exca1
        a = 5
        print(a^2)
        ```
        then:
        \input{output}{./code/exca1}
        done.
        """ * J.EOS |> seval

    h2 = raw"""
        Simple code:
        ```julia:excb1
        a = 5
        print(a^2)
        ```
        then:
        \input{output}{excb1}
        done.
        """ * J.EOS |> seval

    @test h == h2

    spatha = joinpath(J.PATHS[:assets], "index", "code", "exca1.jl")
    spathb = joinpath(J.PATHS[:assets], "index", "code", "excb1.jl")
    @test isfile(spatha)
    @test isfile(spathb)
    @test isapproxstr(read(spatha, String), """
        $(J.MESSAGE_FILE_GEN_JMD)
        a = 5\nprint(a^2)""")

    opath = joinpath(J.PATHS[:assets], "index", "code", "output", "exca1.out")
    @test isfile(opath)
    @test read(opath, String) == "25"

    @test isapproxstr(h, raw"""

Simple code:

a = 5
            print(a^2)

            then:

25

            done.
""")
end

@testset "Eval (errs)" begin
    # see `converter/md_blocks:convert_code_block`
    # --------------------------------------------
    h = raw"""
        Simple code:
        ```python:./scripts/testpy
        a = 5
        print(a**2)
        ```
        done.
        """ * J.EOS |> seval

    @test isapproxstr(h, raw"""

Simple code:

a = 5
            print(a**2)

            done.
""")
end

@testset "Eval (rinput)" begin
    h = raw"""
        Simple code:
        ```julia:/assets/scripts/test2
        a = 5
        print(a^2)
        ```
        then:
        \input{output}{/assets/scripts/test2}
        done.
        """ * J.EOS |> seval

    spath = joinpath(J.PATHS[:assets], "scripts", "test2.jl")
    @test isfile(spath)
    @test occursin("a = 5\nprint(a^2)", read(spath, String))

    opath = joinpath(J.PATHS[:assets], "scripts", "output", "test2.out")
    @test isfile(opath)
    @test read(opath, String) == "25"

    @test isapproxstr(h, """

Simple code:

a = 5
            print(a^2)

            then:

25

            done.
""")

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

    J.CUR_PATH[] = "pages/pg1.md"

    h = raw"""
        Simple code:
        ```julia:./code/abc2
        a = 5
        print(a^2)
        ```
        then:
        \input{output}{./code/abc2}
        done.
        """ * J.EOS |> seval

    spath = joinpath(J.PATHS[:assets], "pages", "pg1", "code", "abc2.jl")
    @test isfile(spath)
    @test occursin("a = 5\nprint(a^2)", read(spath, String))

    opath = joinpath(J.PATHS[:assets], "pages", "pg1", "code", "output" ,"abc2.out")
    @test isfile(opath)
    @test read(opath, String) == "25"

    @test isapproxstr(h, """

Simple code:

a = 5
            print(a^2)

            then:

25
  done.
""")
end

@testset "Eval (module)" begin
    h = raw"""
        Simple code:
        ```julia:scripts/test1
        using LinearAlgebra
        a = [5, 2, 3, 4]
        print(dot(a, a))
        ```
        then:
        \input{output}{scripts/test1}
        done.
        """ * J.EOS |> seval
    # dot(a, a) == 54
    @test occursin("""then:
54
 done.""", h)
end

@testset "Eval (img)" begin
    J.CUR_PATH[] = "index.html"
    h = raw"""
        Simple code:
        ```julia:tv2
        write(joinpath(@__DIR__, "output", "tv2.png"), "blah")
        ```
        then:
        \input{plot}{tv2}
        done.
        """ * J.EOS |> seval
    @test occursin("then:  done.", h)
end

@testset "Eval (throw)" begin
    h = raw"""
        Simple code:
        ```julia:scripts/test1
        sqrt(-1)
        ```
        then:
        \input{output}{scripts/test1}
        done.
        """ * J.EOS |> seval
    # errors silently
    @test occursin("then:
There was an error running the code:\nDomainError", h)
end

@testset "Eval (nojl)" begin
    h = raw"""
        Simple code:
        ```python:scripts/test1
        sqrt(-1)
        ```
        done.
        """ * J.EOS

    @test (@test_logs (:warn, "Eval of non-julia code blocks is not yet supported.") h |> seval) == "
Simple code:
sqrt(-1)
 done.
\n"
end

# temporary fix for 186: make error appear and also use `abspath` in internal include
@testset "eval #186" begin
    h = raw"""
        Simple code:
        ```julia:scripts/test186
        fn = "tempf.jl"
        write(fn, "a = 1+1")
        println("Is this a file? $(isfile(fn))")
        include(abspath(fn))
        println("Now: $a")
        rm(fn)
        ```
        done.
        \output{scripts/test186}
        """ * J.EOS |> seval
    @test isapproxstr(h, raw"""

Simple code:
fn = "tempf.jl"
            write(fn, "a = 1+1")
            println("Is this a file? $(isfile(fn))")
            include(abspath(fn))
            println("Now: $a")
            rm(fn)
 done.
Is this a file? true
            Now: 2

            """)
end

@testset "show" begin
    h = raw"""
        @def hascode = true
        @def reeval = true
        ```julia:ex
        a = 5
        a *= 2
        ```
        \show{ex}
        """ |> jd2html_td
    @test isapproxstr(h, """

a = 5
        a *= 2

10

        """)

    # Show with stdout
    h = raw"""
        @def hascode = true
        @def reeval = true
        ```julia:ex
        a = 5
        println("hello")
        a *= 2
        ```
        \show{ex}
        """ |> jd2html_td
    @test isapproxstr(h, """

a = 5
        println("hello")
        a *= 2

hello
        10

        """)
end


text:
#SQLite

@info "SQLite.jl is available - including functionality to read / write SQLite databases"

using .SQLite

function read_table(::TableIOInterface.SQLiteFormat, filename:: AbstractString, tablename:: AbstractString; kwargs...)
    db = SQLite.DB(filename)
    try
        result = read_table(db, tablename; kwargs...)
        return result
    finally
        close(db)
    end
end

function read_table(::TableIOInterface.SQLiteFormat, filename:: AbstractString; kwargs...)
    db = SQLite.DB(filename)
    try
        result = read_table(db; kwargs...)
        return result
    finally
        close(db)
    end
end

"""
    read_table(db:: SQLite.DB, tablename:: AbstractString; kwargs...)

Reads the content of a whole table and returns a Tables.jl compatible object.
This method takes an instance of an SQLite database connection as input.
"""
function read_table(db:: SQLite.DB, tablename:: AbstractString; kwargs...)
    _checktablename(tablename)
    result = DBInterface.execute(db, "select * from $tablename") # SQL parameters cannot be used for table names
    return result
end

function read_table(db:: SQLite.DB; kwargs...)
    table_list = list_tables(db)
    length(table_list) > 1 && @warn "File contains more than one table, the alphabetically first one is taken"
    return read_table(db, first(table_list); kwargs...)
end

read_sql(db:: SQLite.DB, sql:: AbstractString) = DBInterface.execute(db, sql)

function write_table!(::TableIOInterface.SQLiteFormat, filename:: AbstractString, tablename:: AbstractString, table; kwargs...)
    _checktable(table)
    _checktablename(tablename)
    db = SQLite.DB(filename)
    try
        write_table!(db, tablename, table; kwargs...)
    finally
        close(db)
    end
    nothing
end

"""
    write_table(db:: SQLite.DB, tablename:: AbstractString, table; kwargs...)

Writes the data into a new table in the SQLite database.
As default, this operation will fail if the table already exists. Different behavior can be specified using
the `kwargs` passed to `SQLite.load!`.

Note that this behavior is different to the one for PostgreSQL, where the table must already exist.
"""
function write_table!(db:: SQLite.DB, tablename:: AbstractString, table; kwargs...)
    _checktable(table)
    _checktablename(tablename)
    table |> SQLite.load!(db, tablename; kwargs...)
    nothing
end

function list_tables(::TableIOInterface.SQLiteFormat, filename:: AbstractString)
    db = SQLite.DB(filename)
    try
        return list_tables(db)
    finally
        close(db)
    end
end

function list_tables(db:: SQLite.DB)
    files = SQLite.tables(db).name
    return files |> sort
end


text:
@recipe function f(mach::MLJBase.Machine{<:EitherTunedModel})
    r = report(mach).plotting
    z = r.measurements
    x = r.parameter_values[:,1]
    y = r.parameter_values[:,2]

    r.parameter_scales[1] == :none &&
        (x = string.(x))

    r.parameter_scales[2] == :none &&
        (y = string.(y))

    xsc, ysc = r.parameter_scales

    xguide --> r.parameter_names[1]
    yguide --> r.parameter_names[2]
    xscale --> (xsc in [:custom, :linear] ? :identity : xsc)
    yscale --> (ysc in [:custom, :linear] ? :identity : ysc)

    st = get(plotattributes, :seriestype, :scatter)

    if st ∈ (:surface, :heatmap, :contour, :contourf, :wireframe)
        ux = unique(x)
        uy = unique(y)
        m = reshape(z, (length(ux), length(uy)))'
        ux, uy, m
    else
        label --> ""
        seriestype := st
        ms = get(plotattributes, :markersize, 4)
        markersize := _getmarkersize(ms, z)
        marker_z --> z
        x, y
    end
end

function _getmarkersize(ms, z)
    ret = sqrt.(z)
    minz, maxz = extrema(x for x in ret if isfinite(x))
    ret .-= minz
    ret ./= maxz - minz
    4ms .* ret .+ 1
end


text:
@doc raw"""
Constructs allowed rate-of-change constraints from variables, initial condtions, and rate data.

# Constraints
If t = 1:

``` variable[name, 1] - initial_conditions[ix].value <= rate_data[1][ix].up ```

``` initial_conditions[ix].value - variable[name, 1] <= rate_data[1][ix].down ```

If t > 1:

``` variable[name, t] - variable[name, t-1] <= rate_data[1][ix].up ```

``` variable[name, t-1] - variable[name, t] <= rate_data[1][ix].down ```

# LaTeX

`` r^{down} \leq x_1 - x_{init} \leq r^{up}, \text{ for } t = 1 ``

`` r^{down} \leq x_t - x_{t-1} \leq r^{up}, \forall t \geq 2 ``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* rate_data::Tuple{Vector{String}, Vector{UpDown}} : gives name (1) and max ramp up/down rates (2)
* initial_conditions::Vector{InitialCondition} : for time zero 'variable'
* cons_name::Symbol : name of the constraint
* var_name::Tuple{Symbol, Symbol, Symbol} : the name of the variable
"""
function device_linear_rateofchange!(
    psi_container::PSIContainer,
    rate_data::Vector{DeviceRampConstraintInfo},
    cons_name::Symbol,
    var_name::Symbol,
)
    parameters = model_has_parameters(psi_container)
    time_steps = model_time_steps(psi_container)
    up_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "up")
    down_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "dn")

    variable = get_variable(psi_container, var_name)

    set_name = [get_component_name(r) for r in rate_data]
    con_up = add_cons_container!(psi_container, up_name, set_name, time_steps)
    con_down = add_cons_container!(psi_container, down_name, set_name, time_steps)

    for r in rate_data
        name = get_component_name(r)
        ic_power = get_value(get_ic_power(r))
        @debug "add rate_of_change_constraint" name ic_power
        @assert (parameters && isa(ic_power, PJ.ParameterRef)) || !parameters
        expression_ub = JuMP.AffExpr(0.0, variable[name, 1] => 1.0)
        for val in r.additional_terms_ub
            JuMP.add_to_expression!(
                expression_ub,
                get_variable(psi_container, val)[name, 1],
            )
        end
        con_up[name, 1] = JuMP.@constraint(
            psi_container.JuMPmodel,
            expression_ub - ic_power <= r.ramp_limits.up
        )
        expression_lb = JuMP.AffExpr(0.0, variable[name, 1] => 1.0)
        for val in r.additional_terms_lb
            JuMP.add_to_expression!(
                expression_lb,
                get_variable(psi_container, val)[name, 1],
                1.0,
            )
        end
        con_down[name, 1] = JuMP.@constraint(
            psi_container.JuMPmodel,
            ic_power - expression_lb <= r.ramp_limits.down
        )
    end

    for t in time_steps[2:end], r in rate_data
        name = get_component_name(r)
        expression_ub = JuMP.AffExpr(0.0, variable[name, t] => 1.0)
        for val in r.additional_terms_ub
            JuMP.add_to_expression!(
                expression_ub,
                get_variable(psi_container, val)[name, t],
            )
        end
        con_up[name, t] = JuMP.@constraint(
            psi_container.JuMPmodel,
            expression_ub - variable[name, t - 1] <= r.ramp_limits.up
        )
        expression_lb = JuMP.AffExpr(0.0, variable[name, t] => 1.0)
        for val in r.additional_terms_lb
            JuMP.add_to_expression!(
                expression_lb,
                get_variable(psi_container, val)[name, t],
                1.0,
            )
        end
        con_down[name, t] = JuMP.@constraint(
            psi_container.JuMPmodel,
            variable[name, t - 1] - expression_lb <= r.ramp_limits.down
        )
    end

    return
end

@doc raw"""
Constructs allowed rate-of-change constraints from variables, initial condtions, start/stop status, and rate data

# Equations
If t = 1:

``` variable[name, 1] - initial_conditions[ix].value <= rate_data[1][ix].up + rate_data[2][ix].max*varstart[name, 1] ```

``` initial_conditions[ix].value - variable[name, 1] <= rate_data[1][ix].down + rate_data[2][ix].min*varstop[name, 1] ```

If t > 1:

``` variable[name, t] - variable[name, t-1] <= rate_data[1][ix].up + rate_data[2][ix].max*varstart[name, t] ```

``` variable[name, t-1] - variable[name, t] <= rate_data[1][ix].down + rate_data[2][ix].min*varstop[name, t] ```

# LaTeX

`` r^{down} + r^{min} x^{stop}_1 \leq x_1 - x_{init} \leq r^{up} + r^{max} x^{start}_1, \text{ for } t = 1 ``

`` r^{down} + r^{min} x^{stop}_t \leq x_t - x_{t-1} \leq r^{up} + r^{max} x^{start}_t, \forall t \geq 2 ``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* rate_data::Tuple{Vector{String}, Vector{UpDown}, Vector{MinMax}} : (1) gives name
                                                                     (2) gives min/max ramp rates
                                                                     (3) gives min/max for 'variable'
* initial_conditions::Vector{InitialCondition} : for time zero 'variable'
* cons_name::Symbol : name of the constraint
* var_names::Tuple{Symbol, Symbol, Symbol} : the names of the variables
- : var_names[1] : 'variable'
- : var_names[2] : 'varstart'
- : var_names[3] : 'varstop'
"""
function device_mixedinteger_rateofchange!(
    psi_container::PSIContainer,
    rate_data::Vector{DeviceRampConstraintInfo},
    cons_name::Symbol,
    var_names::Tuple{Symbol, Symbol, Symbol},
)
    parameters = model_has_parameters(psi_container)
    time_steps = model_time_steps(psi_container)
    up_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "up")
    down_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "dn")

    variable = get_variable(psi_container, var_names[1])
    varstart = get_variable(psi_container, var_names[2])
    varstop = get_variable(psi_container, var_names[3])

    set_name = [get_component_name(r) for r in rate_data]
    con_up = add_cons_container!(psi_container, up_name, set_name, time_steps)
    con_down = add_cons_container!(psi_container, down_name, set_name, time_steps)

    for r in rate_data
        name = get_component_name(r)
        ic_power = get_value(get_ic_power(r))
        @debug "add rate_of_change_constraint" name ic_power
        @assert (parameters && isa(ic_power, PJ.ParameterRef)) || !parameters
        expression_ub = JuMP.AffExpr(0.0, variable[name, 1] => 1.0)
        for val in r.additional_terms_ub
            JuMP.add_to_expression!(
                expression_ub,
                get_variable(psi_container, val)[name, 1],
            )
        end
        con_up[name, 1] = JuMP.@constraint(
            psi_container.JuMPmodel,
            expression_ub - (ic_power) <=
            r.ramp_limits.up + r.limits.max * varstart[name, 1]
        )
        expression_lb = JuMP.AffExpr(0.0, variable[name, 1] => 1.0)
        for val in r.additional_terms_lb
            JuMP.add_to_expression!(
                expression_lb,
                get_variable(psi_container, val)[name, 1],
                1.0,
            )
        end
        con_down[name, 1] = JuMP.@constraint(
            psi_container.JuMPmodel,
            (ic_power) - expression_lb <=
            r.ramp_limits.down + r.limits.min * varstop[name, 1]
        )
    end

    for t in time_steps[2:end], r in rate_data
        name = get_component_name(r)
        expression_ub = JuMP.AffExpr(0.0, variable[name, t] => 1.0)
        for val in r.additional_terms_ub
            JuMP.add_to_expression!(
                expression_ub,
                get_variable(psi_container, val)[name, t],
            )
        end
        con_up[name, t] = JuMP.@constraint(
            psi_container.JuMPmodel,
            expression_ub - variable[name, t - 1] <=
            r.ramp_limits.up + r.limits.max * varstart[name, 1]
        )
        expression_lb = JuMP.AffExpr(0.0, variable[name, t] => 1.0)
        for val in r.additional_terms_lb
            JuMP.add_to_expression!(
                expression_lb,
                get_variable(psi_container, val)[name, t],
                1.0,
            )
        end
        con_down[name, t] = JuMP.@constraint(
            psi_container.JuMPmodel,
            variable[name, t - 1] - expression_lb <=
            r.ramp_limits.down + r.limits.min * varstop[name, t]
        )
    end

    return
end

@doc raw"""
Constructs allowed rate-of-change constraints from variables, initial condtions, start/stop status, and rate data

# Equations
If t = 1:

``` variable[name, 1] - initial_conditions[ix].value <= rate_data[1][ix].up ```

``` initial_conditions[ix].value - variable[name, 1] <= rate_data[1][ix].down ```

If t > 1:

``` variable[name, t] - variable[name, t-1] <= rate_data[1][ix].up  ```

``` variable[name, t-1] - variable[name, t] <= rate_data[1][ix].down ```

# LaTeX

`` r^{down}  \leq x_1 - x_{init} \leq r^{up}  \text{ for } t = 1 ``

`` r^{down} \leq x_t - x_{t-1} \leq r^{up}  \forall t \geq 2 ``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* rate_data::Tuple{Vector{String}, Vector{UpDown}, Vector{MinMax}} : (1) gives name
                                                                     (2) gives min/max ramp rates
                                                                     (3) gives min/max for 'variable'
* initial_conditions::Vector{InitialCondition} : for time zero 'variable'
* cons_name::Symbol : name of the constraint
* var_names::Tuple{Symbol, Symbol, Symbol} : the names of the variables
- : var_name : 'variable'
"""
function device_multistart_rateofchange!(
    psi_container::PSIContainer,
    rate_data::Vector{DeviceRampConstraintInfo},
    cons_name::Symbol,
    var_name::Symbol,
)
    time_steps = model_time_steps(psi_container)
    up_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "up")
    down_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "dn")

    variable = get_variable(psi_container, var_name)

    set_name = [get_component_name(r) for r in rate_data]
    con_up = add_cons_container!(psi_container, up_name, set_name, time_steps)
    con_down = add_cons_container!(psi_container, down_name, set_name, time_steps)

    for r in rate_data
        name = get_component_name(r)
        ic_power = get_value(get_ic_power(r))
        expression_ub = JuMP.AffExpr(0.0, variable[name, 1] => 1.0)
        for val in r.additional_terms_ub
            JuMP.add_to_expression!(
                expression_ub,
                get_variable(psi_container, val)[name, 1],
            )
        end
        con_up[name, 1] = JuMP.@constraint(
            psi_container.JuMPmodel,
            expression_ub - (ic_power) <= r.ramp_limits.up
        )
        expression_lb = JuMP.AffExpr(0.0, variable[name, 1] => 1.0)
        for val in r.additional_terms_lb
            JuMP.add_to_expression!(
                expression_lb,
                get_variable(psi_container, val)[name, 1],
                1.0,
            )
        end
        con_down[name, 1] = JuMP.@constraint(
            psi_container.JuMPmodel,
            (ic_power) - expression_lb <= r.ramp_limits.down
        )
    end

    for t in time_steps[2:end], r in rate_data
        name = get_component_name(r)
        expression_ub = JuMP.AffExpr(0.0, variable[name, t] => 1.0)
        for val in r.additional_terms_ub
            JuMP.add_to_expression!(
                expression_ub,
                get_variable(psi_container, val)[name, t],
            )
        end
        con_up[name, t] = JuMP.@constraint(
            psi_container.JuMPmodel,
            expression_ub - variable[name, t - 1] <= r.ramp_limits.up
        )
        expression_lb = JuMP.AffExpr(0.0, variable[name, t] => 1.0)
        for val in r.additional_terms_lb
            JuMP.add_to_expression!(
                expression_lb,
                get_variable(psi_container, val)[name, t],
                1.0,
            )
        end
        con_down[name, t] = JuMP.@constraint(
            psi_container.JuMPmodel,
            variable[name, t - 1] - expression_lb <= r.ramp_limits.down
        )
    end

    return
end


text:
module TestIgnoreCase
end # module


text:
using Test
using Plots
using Profile
using FlameGraphs
using Polynomials
using Sirion
#include("../src/Sirion.jl")

function calcRelativeDiff(x::Float64, y::Float64)::Float64
    return abs(x-y)/(abs(x)+abs(y))
end

function runRiemannProblem1()

    # Initialise simulation
    grid = range(0.0,1.0,length=100)
    d_func = x->1
    p_func = x->10-9*(sign(x-0.5)+1)/2
    v_func = x->0
    init_cond = Sirion.initHydroSnapshot(grid, d_func, p_func, v_func)
    eos = Sirion.IdealGas(5.0/3.0)
    tsf = Sirion.CFL(0.3)
    vm = Sirion.Eulerian()
    rs = Sirion.HLLC()
    bc = Sirion.RigidWall(rs)
    fc = Sirion.SimpleFluxCalculator(rs, bc)
    pg = Sirion.Planar()
    st = Sirion.ZeroSource()
    cu = Sirion.SimpleCellUpdater()
    sim = Sirion.initHDSim(init_cond,
                          eos,
                          tsf,
                          vm,
                          fc,
                          pg,
                          st,
                          cu)

    # Run simulation
    while sim.time<0.067
        Sirion.timeAdvance(sim)
    end

    # Perform diagnostics
    pressure_probe = sim.state.cells[length(sim.state.cells)÷2].pressure
    velocity_probe = sim.state.cells[length(sim.state.cells)÷2].velocity
    reference_pressure = 5.1123
    reference_velocity = 1.53795

    cond1 = calcRelativeDiff(reference_pressure, pressure_probe)<0.01
    cond2 = calcRelativeDiff(reference_velocity, velocity_probe)<0.01

    return cond1 && cond2
end

function runRiemannProblem2()

    # Initialise simulation
    grid = range(0.0,1.0,length=100)
    d_func = x->1
    p_func = x->1
    v_func = x->-5*sign(x-0.5)
    init_cond = Sirion.initHydroSnapshot(grid, d_func, p_func, v_func)
    eos = Sirion.IdealGas(5.0/3.0)
    tsf = Sirion.CFL(0.3)
    vm = Sirion.Eulerian()
    rs = Sirion.HLLC()
    bc = Sirion.RigidWall(rs)
    fc = Sirion.SimpleFluxCalculator(rs, bc)
    pg = Sirion.Planar()
    st = Sirion.ZeroSource()
    cu = Sirion.SimpleCellUpdater()
    sim = Sirion.initHDSim(init_cond,
                          eos,
                          tsf,
                          vm,
                          fc,
                          pg,
                          st,
                          cu)

    # Run simulation
    while sim.time<0.08
        Sirion.timeAdvance(sim)
    end

    # Perform diagnostics
    pressure_probe = sim.state.cells[1+length(sim.state.cells)÷2].pressure
    velocity_probe = sim.state.cells[1+length(sim.state.cells)÷2].velocity
    reference_pressure = 35.5397
    reference_velocity = 0.0

    cond1 = calcRelativeDiff(reference_pressure, pressure_probe)<0.1
    cond2 = abs(velocity_probe)<0.01

    return cond1 && cond2
end

function runRiemannProblem3()

    # Initialise simulation
    grid = range(0.0,1.0,length=100)
    d_func = x->1
    p_func = x->1
    v_func = x->1.0*sign(x-0.5)
    init_cond = Sirion.initHydroSnapshot(grid, d_func, p_func, v_func)
    eos = Sirion.IdealGas(5.0/3.0)
    tsf = Sirion.CFL(0.3)
    vm = Sirion.Eulerian()
    rs = Sirion.HLLC()
    bc = Sirion.RigidWall(rs)
    fc = Sirion.SimpleFluxCalculator(rs, bc)
    pg = Sirion.Planar()
    st = Sirion.ZeroSource()
    cu = Sirion.SimpleCellUpdater()
    sim = Sirion.initHDSim(init_cond,
                          eos,
                          tsf,
                          vm,
                          fc,
                          pg,
                          st,
                          cu)

    # Run simulation
    while sim.time<0.22
        Sirion.timeAdvance(sim)
    end

    # Perform diagnostics
    pressure_probe = sim.state.cells[1+length(sim.state.cells)÷2].pressure
    velocity_probe = sim.state.cells[1+length(sim.state.cells)÷2].velocity
    reference_pressure = 0.224614
    reference_velocity = 0.0

    cond1 = calcRelativeDiff(reference_pressure, pressure_probe)<0.01
    cond2 = abs(velocity_probe)<0.01

    return cond1 && cond2
end

function l1_fit(a1, a2)
    abs_dif = [abs(x-y) for (x,y) in zip(a1, a2)]
    return sum(abs_dif)/length(abs_dif)
end

function runPlanarNoh()

    # Initialise simulation
    grid = range(0.0,1.0,length=100)
    d_func = x->1
    p_func = x->1e-6
    v_func = x->-1.0
    init_cond = Sirion.initHydroSnapshot(grid, d_func, p_func, v_func)
    g = 5.0/3.0
    eos = Sirion.IdealGas(g)
    tsf = Sirion.CFL(0.3)
    vm = Sirion.Eulerian()
    rs = Sirion.HLLC()
    bc = Sirion.DifferentBC(Sirion.RigidWall(rs),
                            Sirion.FreeFlow(rs))
    fc = Sirion.SimpleFluxCalculator(rs, bc)
    pg = Sirion.Planar()
    st = Sirion.ZeroSource()
    cu = Sirion.SimpleCellUpdater()
    sim = Sirion.initHDSim(init_cond,
                          eos,
                          tsf,
                          vm,
                          fc,
                          pg,
                          st,
                          cu)

    # Run simulation
    while sim.time<2.0
        Sirion.timeAdvance(sim)
    end

    # Perform diagnostics
    upstream_speed = 1.0
    shock_speed = (g-1)*upstream_speed/2
    shock_position = shock_speed*sim.time
    upstream_density = 1
    downstream_density = (g+1)*upstream_density/(g-1)
    upstream_velocity = -1
    downstream_velocity = 0
    upstream_pressure = 1e-6
    downstream_pressure = ((g+1)/2)*upstream_density*
        (downstream_velocity-upstream_velocity)^2
    d_func2(x) = x<shock_position ? downstream_density : upstream_density
    p_func2(x) = x<shock_position ? downstream_pressure : upstream_pressure
    v_func2(x) = x<shock_position ? downstream_velocity : upstream_velocity
    analytic_d = [d_func2(x) for x in Sirion.midGrid(sim.state.grid)]
    analytic_p = [p_func2(x) for x in Sirion.midGrid(sim.state.grid)]
    analytic_v = [v_func2(x) for x in Sirion.midGrid(sim.state.grid)]
    cond1 = l1_fit(analytic_d, [itm.density for itm in sim.state.cells])<0.03
    cond2 = l1_fit(analytic_d, [itm.density for itm in sim.state.cells])<0.03
    cond3 = l1_fit(analytic_v, [itm.velocity for itm in sim.state.cells])<0.01
    return (cond1 && cond2 && cond3)

end

function runSedovTaylor()

    # Initialise simulation
    #grid = range(1e-6,1.0,length=200)
    grid = 10.0.^(range(-1.5, 0.0, length = 100))
    #d_func = x->x<1e-2 ? 1e4 : x^-2
    d_func = x->1
    p_func = x->x<1e-1 ? 1e6 : 1e-9
    v_func = x->0
    init_cond = Sirion.initHydroSnapshot(grid, d_func, p_func, v_func)
    g = 5.0/3.0
    eos = Sirion.IdealGas(g)
    tsf = Sirion.CFL(0.3)
    vm = Sirion.Eulerian()
    rs = Sirion.HLLC()
    bc = Sirion.RigidWall(rs)
    fc = Sirion.SimpleFluxCalculator(rs, bc)
    pg = Sirion.Spherical()
    #pg = Sirion.Planar()
    st = Sirion.SphericalComplamentary()
    #st = Sirion.ZeroSource()
    cu = Sirion.SimpleCellUpdater()
    sim = Sirion.initHDSim(init_cond,
                          eos,
                          tsf,
                          vm,
                          fc,
                          pg,
                          st,
                          cu)

    # Run simulation
    #while sim.time<3e-3
    time_list = Array{Float64}(undef,0)
    radius_list = Array{Float64}(undef,0)
    while sim.time<2e-3
        Sirion.timeAdvance(sim)
        mxval, mxind = findmax([c.velocity for c in sim.state.cells])
        push!(time_list, sim.time)
        push!(radius_list, sim.state.grid[mxind])
    end

    r_thres = 0.15
    fitd = fit(map(log,time_list[radius_list .> r_thres]),
               map(log,radius_list[radius_list .> r_thres]),
               1)
    #println(fitd)
    #plot(time_list[radius_list .> r_thres],
    #     radius_list[radius_list .> r_thres],
    #     xaxis=:log,
    #     yaxis=:log)
    #plot!(time_list[radius_list .> r_thres],
    #      [exp(fitd(log(x)))
    #       for x in time_list[radius_list .> r_thres]],
    #      xaxis=:log,
    #      yaxis=:log)

    return abs(fitd[1]-0.4)<2e-3
end

function runSedovTaylorLagrangian()

    # Initialise simulation
    #grid = range(1e-6,1.0,length=200)
    grid = 10.0.^(range(-1.5, 0.0, length = 100))
    #d_func = x->x<1e-2 ? 1e4 : x^-2
    d_func = x->1
    p_func = x->x<1e-1 ? 1e6 : 1e-9
    v_func = x->0
    init_cond = Sirion.initHydroSnapshot(grid, d_func, p_func, v_func)
    g = 5.0/3.0
    eos = Sirion.IdealGas(g)
    tsf = Sirion.CFL(0.3)
    vm = Sirion.Lagrangian(Sirion.KeepEdgesFixed())
    #vm = Sirion.Eulerian()
    rs = Sirion.HLLC()
    bc = Sirion.RigidWall(rs)
    fc = Sirion.SimpleFluxCalculator(rs, bc)
    pg = Sirion.Spherical()
    #pg = Sirion.Planar()
    st = Sirion.SphericalComplamentary()
    #st = Sirion.ZeroSource()
    cu = Sirion.SimpleCellUpdater()
    sim = Sirion.initHDSim(init_cond,
                          eos,
                          tsf,
                          vm,
                          fc,
                          pg,
                          st,
                          cu)

    # Run simulation
    #while sim.time<3e-3
    time_list = Array{Float64}(undef,0)
    radius_list = Array{Float64}(undef,0)
    while sim.time<5e-3
        Sirion.timeAdvance(sim)
        mxval, mxind = findmax([c.velocity for c in sim.state.cells])
        push!(time_list, sim.time)
        push!(radius_list, sim.state.grid[mxind])
    end

    r_thres = 0.2
    fitd = fit(map(log,time_list[radius_list .> r_thres]),
               map(log,radius_list[radius_list .> r_thres]),
               1)
    #plot(time_list[radius_list .> r_thres],
    #     radius_list[radius_list .> r_thres],
    #     xaxis=:log,
    #     yaxis=:log)
    #plot!(time_list[radius_list .> r_thres],
    #      [exp(fitd(log(x)))
    #       for x in time_list[radius_list .> r_thres]],
    #      xaxis=:log,
    #      yaxis=:log)

    return abs(fitd[1]-0.4)<2e-3
end

function runPrimakoff()

    # Initialise simulation
    #grid = range(1e-6,1.0,length=200)
    grid = 10.0.^(range(-2.5, 0.0, length = 100))
    d_func = x->x<1e-2 ? 1e4 : x^-2
    #d_func = x->1
    p_func = x->x<1e-2 ? 1e6 : 1e-9
    v_func = x->0
    init_cond = Sirion.initHydroSnapshot(grid, d_func, p_func, v_func)
    g = 5.0/3.0
    eos = Sirion.IdealGas(g)
    tsf = Sirion.CFL(0.3)
    vm = Sirion.Eulerian()
    rs = Sirion.HLLC()
    bc = Sirion.RigidWall(rs)
    fc = Sirion.SimpleFluxCalculator(rs, bc)
    pg = Sirion.Spherical()
    #pg = Sirion.Planar()
    st = Sirion.SphericalComplamentary()
    #st = Sirion.ZeroSource()
    cu = Sirion.SimpleCellUpdater()
    sim = Sirion.initHDSim(init_cond,
                          eos,
                          tsf,
                          vm,
                          fc,
                          pg,
                          st,
                          cu)

    # Run simulation
    #while sim.time<3e-3
    while sim.time<3e-1
        Sirion.timeAdvance(sim)
        if sim.cycle%1000==0
            println(sim.time)
        end
    end

    # Perform diagnostics
    p1 = plot(Sirion.midGrid(sim.state.grid),
              [c.density for c in sim.state.cells],
              legend=false)
    p2 = plot(Sirion.midGrid(sim.state.grid),
              [c.pressure for c in sim.state.cells],
              legend=false)
    p3 = plot(Sirion.midGrid(sim.state.grid),
              [c.velocity for c in sim.state.cells],
              legend=false)
    p = plot(p1,p2,p3,layout=(3,1))
    display(p)

    #return false
end

@test runRiemannProblem1()
@test runRiemannProblem2()
@test runRiemannProblem3()
@test runPlanarNoh()
@test runSedovTaylor()
@test runSedovTaylorLagrangian()
#@profile runPrimakoff()
#Juno.profiler()
#Profile.clear()
#open("./prof.txt", "w") do s
#    Profile.print(IOContext(s, :displaysize => (24, 500)))
#end

max_stars_repo_name: bolverk/Sirion

max_issues_repo_name: bolverk/Sirion

max_forks_repo_name: bolverk/Sirion


text:
@testset "float" begin
    @testset "eps" begin
        @test one(Arf) + eps(Arf) > one(Arf)
        @test one(Arf) + eps(Arf) / 2 == one(Arf)
        @test Arf(4) + eps(Arf(4)) > Arf(4)
        @test Arf(4) + eps(Arf(4)) / 2 == Arf(4)
        @test Arf(4, prec = 64) + eps(Arf(4, prec = 64)) > Arf(4)
        @test Arf(4, prec = 64) + eps(Arf(4, prec = 64)) / 2 == Arf(4)

        @test eps(Arf) == Arf(eps(BigFloat))
        @test eps(Arf(4)) == Arf(eps(BigFloat(4)))
        @test Float64(eps(Arf(1, prec = 53))) == eps()
        @test Float64(eps(Arf(4, prec = 53))) == eps(4.0)

        @test eps(Arf) == Arblib.eps!(zero(Arf), one(Arf))
        @test eps(Arf(4)) == Arblib.eps!(one(Arf), Arf(4))
        @test eps(Arb) == Arblib.eps!(zero(Arb), one(Arb))
        @test eps(Arb(4)) == Arblib.eps!(one(Arb), Arb(4))

        @test eps(Arb) == Arb(eps(Arf))

        @test eps(Arf) isa Arf
        @test eps(ArfRef) isa Arf
        @test eps(one(Arf)) isa Arf
        @test eps(Arb) isa Arb
        @test eps(ArbRef) isa Arb
        @test eps(one(Arb)) isa Arb

        @test precision(eps(Arf(1, prec = 64))) == 64
        @test precision(eps(Arb(1, prec = 64))) == 64

        # Test aliasing
        x = Arf(4)
        @test Arblib.eps!(x, x) == eps(Arf(4))
        x = Arb(4)
        @test Arblib.eps!(x, x) == eps(Arb(4))

        @test isnan(eps(zero(Arf)))
        @test isnan(eps(zero(Arb)))
        @test isnan(eps(Arf(Inf)))
        @test isnan(eps(Arb(Inf)))
        @test isnan(eps(Arf(NaN)))
        @test isnan(eps(Arb(NaN)))
    end
end


text:
module TimeFreqFFTs

using LinearAlgebra
using FFTW

using ..Lattices:   Lattice
using ..Utilities: reshaped

struct TimeFreqFFT{T<:AbstractFloat,Tfft<:AbstractFFTs.Plan,Tifft<:AbstractFFTs.Plan}

    "Number of site in lattice."
    N::Int

    "Length of imaginary time axis."
    L::Int

    "Vector to represent unitary transformation to restore translational invariance in imaginary time direction."
    Θ::Vector{Complex{T}}

    "FFT plan."
    fftplan::Tfft

    "Inverse FFT plan."
    ifftplan::Tifft

    "Temporary storage vector."
    vtemp::Array{Complex{T},2}

    "Temporary storage vector."
    utemp::Vector{Complex{T}}

    function TimeFreqFFT(T::DataType,N::Int,L::Int)

        @assert T <: AbstractFloat
        vtemp    = zeros(Complex{T},L,N)
        utemp    = zeros(Complex{T},N*L)
        Θ        = [exp(-π*im*(τ-1)/L) for τ = 1:L]
        fftplan  = plan_fft(vtemp, (1,), flags=FFTW.PATIENT)
        ifftplan = plan_ifft(vtemp, (1,), flags=FFTW.PATIENT)
        Tfft     = typeof(fftplan)
        Tifft    = typeof(ifftplan)

        return new{T,Tfft,Tifft}(N,L,Θ,fftplan,ifftplan,vtemp,utemp)
    end
end

function TimeFreqFFT(lattice::Lattice{T},L::Int) where {T<:AbstractFloat}

    return TimeFreqFFT(T,lattice.nsites,L)
end

"""
Apply the transformation ν=[F⋅Θ]⋅v from imaginary time τ to frequency ω space.
"""
function τ_to_ω!(vout::AbstractVector{Complex{T1}},op::TimeFreqFFT{T1},vin::AbstractVector{T2}) where {T1<:AbstractFloat,T2<:Number}

    L  = op.L::Int
    N  = op.N::Int
    Θ  = op.Θ::Vector{Complex{T1}}
    vtemp = op.vtemp::Array{Complex{T1},2}
    # apply unitary transformation to restore translation invariance in the
    # imaginary time direciton.
    uin  = reshaped(vin,L,N)
    @fastmath @inbounds for i in 1:N
        for τ in 1:L
            vtemp[τ,i] = Θ[τ] * uin[τ,i]
        end
    end
    # apply (τ ⟶ ω) FFT
    uout = reshaped(vout,L,N)
    mul!(uout,op.fftplan,vtemp)
    return nothing
end

function τ_to_ω!(vout::AbstractVector{T1},op::TimeFreqFFT{T1},vin::AbstractVector{T2}) where {T1<:AbstractFloat,T2<:Complex}

    τ_to_ω!(op.utemp,op,vin)
    @. vout = real(op.utemp)
    return nothing
end

function τ_to_ω!(v::AbstractVector{Complex{T}},op::TimeFreqFFT{T}) where {T<:AbstractFloat}

    τ_to_ω!(v,op,v)
    return nothing
end

"""
Apply the transformation v=[Θᵀ⋅Fᵀ]⋅ν from frequency ω to imaginary time τ space.
"""
function ω_to_τ!(vout::AbstractVector{Complex{T}},op::TimeFreqFFT{T},vin::AbstractVector{Complex{T}}) where {T<:AbstractFloat}

    L     = op.L::Int
    N     = op.N::Int
    Θ     = op.Θ::Vector{Complex{T}}
    vtemp = op.vtemp::Array{Complex{T},2}
    # apply (ω ⟶ τ) FFT
    uin  = reshaped(vin, L,N)
    mul!(vtemp,op.ifftplan,uin)
    # apply unitary transformation to restore anitperiodic boundary conditions
    # imaginary time direciton.
    uout = reshaped(vout,L,N)
    @fastmath @inbounds for i in 1:N
        for τ in 1:L
            uout[τ,i] = conj(Θ[τ]) * vtemp[τ,i]
        end
    end
    return nothing
end

function ω_to_τ!(vout::AbstractVector{T},op::TimeFreqFFT{T},vin::AbstractVector{Complex{T}}) where {T<:AbstractFloat}

    L     = op.L::Int
    N     = op.N::Int
    Θ     = op.Θ::Vector{Complex{T}}
    vtemp = op.vtemp::Array{Complex{T},2}
    # apply (ω ⟶ τ) FFT
    uin  = reshaped(vin, L,N)
    mul!(vtemp,op.ifftplan,uin)
    # apply unitary transformation to restore anitperiodic boundary conditions
    # imaginary time direciton.
    uout = reshaped(vout,L,N)
    @fastmath @inbounds for i in 1:N
        for τ in 1:L
            uout[τ,i] = real( conj(Θ[τ]) * vtemp[τ,i] )
        end
    end
    return nothing
end

function ω_to_τ!(vout::AbstractVector{T2},op::TimeFreqFFT{T1},vin::AbstractVector{T1};thresh::T1=1e-12) where {T1<:AbstractFloat,T2<:Complex}

    copyto!(op.vtemp,vin)
    ω_to_τ!(vout,op,op.vtemp)
    return nothing
end

function ω_to_τ!(v::AbstractVector{Complex{T}},op::TimeFreqFFT{T}) where {T<:AbstractFloat}

    copyto!(op.vtemp,v)
    ω_to_τ!(v,op,op.vtemp)
    return nothing
end

end

max_stars_repo_name: cohensbw/Langevin

max_issues_repo_name: cohensbw/ElPhDynamics

max_forks_repo_name: cohensbw/ElPhDynamics


text:
module Histograms

using LinearAlgebra

export Histogram, duplicate, reset_distribution, ndim, get_data, get_values,
  get_buffer, swap_buffer!, copy_filter!, get_range, generate_prior,
  test_histogram, weighted_average, size

import Base.size

mutable struct Histogram{RangeType <: Real}
  range::Tuple{Vector{RangeType},Vector{RangeType}}
  data::Array{Float64, 2}

  buffer::Array{Float64, 2}

  # The default constructor specializes for when we can infer the type of the
  # range from the signature. (Alternatively, see the outer constructor).
  function Histogram(range::Tuple{Vector{D}, Vector{D}}, data) where D <: Real
    new{D}(range, data, Array{Float64}(undef, size(data)))
  end
end

include("sparse_histograms.jl")

# Outer constructor that defers determination of the histogram type until after
# pulling the ranges
Histogram(range, data) = Histogram(map(collect, range), data)

Histogram(range) = Histogram(range, generate_prior(range))

function reset_distribution!(x::Histogram)
  x.data = generate_prior(x)
  Void
end

# Copy constructor. Note that this only duplicates the data.
# We assume that nobody is crazy enough to modify the range.
duplicate(x::Histogram) = Histogram(x)
function Histogram(x::Histogram)
  Histogram(get_range(x), Array(get_data(x)))
end

ndim(x::AnyHistogram) = length(x.range)

size(x::Histogram) = size(x.data)

# Returns the data matrix
get_data(x::Histogram) = x.data
# Returns the data values themselves (as in for sparse matrices)
get_values(x::AnyHistogram) = x.data
get_buffer(x::AnyHistogram) = x.buffer

function swap_buffer!(x::AnyHistogram)
  old_data = x.data

  x.data = x.buffer
  x.buffer = old_data
end

# Copies filter data (exclusively)
function copy_filter!(x::Histogram; out::Histogram)
  out.data .= x.data
end

get_range(x::AnyHistogram) = x.range
get_range(x::AnyHistogram, index) = x.range[index]

function generate_prior(range)
  lengths = map(length, range)
  num_cells = prod(lengths)

  ones(lengths) / num_cells
end
generate_prior(x::Histogram) = generate_prior(x.range)

weighted_average(x::Histogram, dim) = dot(sum_all_dims_but(x.data, dim), collect(x.range[dim]))
weighted_average(h::Histogram) = map(x->weighted_average(h, x), 1:length(h.range))

# Maps indices in the ranges to the space they represent
function from_indices(h::Histogram, inds)
  n = length(h.range)
  out = Array{Float64}(undef, n)
  for ii = 1:n
    out[ii] = h.range[ii][inds[ii]]
  end
  out
end

###########
# Test code
###########

test_total_probability(x::Histogram) = abs(sum(get_data(x)) - 1) < 1e-3

test_positive(x::Histogram) = all(get_data(x) .>= 0)

function test_histogram()
  h = Histogram((1:0.5:2, 1:0.5:2))

  assert(test_total_probability(h))
  assert(test_positive(h))
end

function sum_all_dims_but(data::Array{T,N}, dim) where {T, N}
  all_but = filter(x->x!=dim, 1:N)
  sum(data, all_but)[:]
end

end

max_stars_repo_name: mcorah/CenterOfMass

max_issues_repo_name: mcorah/CenterOfMass

max_forks_repo_name: mcorah/CenterOfMass


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ ad211b82-c497-4081-8ff4-4e82d43a3db9
using PlutoUI

# ╔═╡ 17677ae6-e264-11eb-3006-77b002930765
md"""
# Chapter 7
#### Nicetak
 $(import Dates; Dates.format(Dates.today(), Dates.DateFormat("U d, Y")))
"""

# ╔═╡ ad98b4fc-f7d9-4917-896d-a3fa783f0ebf
md"""
## 7.1
 $N$個の整数$a_0, a_1, \dots, a_{N-1}$と, $N$個の整数$b_0, b_1, \dots, b_{N-1}$が与えられます. $a_0, a_1, \dots, a_{N-1}$から何個かと$b_0, b_1, \dots, b_{N-1}$から何個かを選んでペアを作ります, ただし各ペア$(a_i, b_j)$は$a_i < b_j$を満たさなければなりません. 最大で何ペア作れるかを$O(N \log N)$で求めるアルゴリズムを設計してください.
"""

# ╔═╡ db0fd574-7e1f-466c-97e6-52217e95adbb
function f₁(as, bs)
	N = length(as)
	sort!(as)
	sort!(bs)

	ans = 0
	for b in bs
		ans += searchsortedfirst(as, b) - 1
	end

	return ans
end

# ╔═╡ 72ebfd16-144f-4034-a407-22e265fc8236
with_terminal() do
	N = 5
	as = rand(1:N, N)
	bs = rand(1:N, N)

	@show as
	@show bs
	@show f₁(as, bs)
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.9"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "81690084b6198a2e1da36fcfda16eeca9f9f24e4"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.1"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "c8abc88faa3f7a3950832ac5d6e690881590d6dc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.1.0"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╟─17677ae6-e264-11eb-3006-77b002930765
# ╟─ad211b82-c497-4081-8ff4-4e82d43a3db9
# ╟─ad98b4fc-f7d9-4917-896d-a3fa783f0ebf
# ╠═db0fd574-7e1f-466c-97e6-52217e95adbb
# ╟─72ebfd16-144f-4034-a407-22e265fc8236
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
module LearningAlgebraicVarieties

    import Plots
    import LinearAlgebra
    import Statistics
    import FixedPolynomials
    const FP = FixedPolynomials
    import MultivariatePolynomials
    const MP = MultivariatePolynomials
    import DynamicPolynomials
    import Combinatorics
    import RowEchelon
    import Clustering
    import Colors
    import Eirene

    include("auxiliary_functions.jl")
    include("estimate_equations.jl")
    include("estimate_dimension.jl")

end


text:
# ---
# title: 442. Find All Duplicates in an Array
# id: problem442
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/find-all-duplicates-in-an-array/description/>
# hidden: true
# ---
#
# Given an array of integers, 1 ≤ a[i] ≤ _n_ ( _n_ = size of array), some
# elements appear **twice** and others appear **once**.
#
# Find all the elements that appear **twice** in this array.
#
# Could you do it without extra space and in O( _n_ ) runtime?
#
# **Example:**
#
#
#
#     Input:
#     [4,3,2,7,8,2,3,1]
#
#     Output:
#     [2,3]
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end