text:
"""
A module for rendering `Document` objects to LaTeX and PDF.

# Keywords

[`LaTeXWriter`](@ref) uses the following additional keyword arguments that can be passed to
[`makedocs`](@ref Documenter.makedocs): `authors`, `sitename`.

**`sitename`** is the site's title displayed in the title bar and at the top of the
navigation menu. It goes into the `\\title` LaTeX command.

**`authors`** can be used to specify the authors of. It goes into the `\\author` LaTeX command.

"""
module LaTeXWriter
import ...Documenter: Documenter

"""
    LaTeXWriter.LaTeX(; kwargs...)

Output format specifier that results in LaTeX/PDF output.
Used together with [`makedocs`](@ref Documenter.makedocs), e.g.

```julia
makedocs(
    format = LaTeX()
)
```

The `makedocs` argument `sitename` will be used for the `\\title` field in the tex document,
and if the build is for a release tag (i.e. when the `"TRAVIS_TAG"` environment variable is set)
the version number will be appended to the title.
The `makedocs` argument `authors` should also be specified, it will be used for the
`\\authors` field in the tex document.

# Keyword arguments

**`platform`** sets the platform where the tex-file is compiled, either `"native"` (default), `"docker"`,
or "none" which doesn't compile the tex.

See [Other Output Formats](@ref) for more information.
"""
struct LaTeX <: Documenter.Writer
    platform::String
    function LaTeX(; platform = "native")
        platform ∈ ("native", "docker", "none") || throw(ArgumentError("unknown platform: $platform"))
        return new(platform)
    end
end

import ...Documenter:
    Anchors,
    Builder,
    Documents,
    Expanders,
    Documenter,
    Utilities,
    Writers

import Markdown

mutable struct Context{I <: IO} <: IO
    io::I
    in_header::Bool
    footnotes::Dict{String, Int}
    depth::Int
    filename::String # currently active source file
end
Context(io) = Context{typeof(io)}(io, false, Dict(), 1, "")

_print(c::Context, args...) = Base.print(c.io, args...)
_println(c::Context, args...) = Base.println(c.io, args...)

const STYLE = joinpath(dirname(@__FILE__), "..", "..", "assets", "latex", "documenter.sty")

hastex() = (try; success(`latexmk -version`); catch; false; end)

const DOCUMENT_STRUCTURE = (
    "part",
    "chapter",
    "section",
    "subsection",
    "subsubsection",
    "paragraph",
    "subparagraph",
)

# https://github.com/JuliaLang/julia/pull/32851
function mktempdir(args...; kwargs...)
    if VERSION < v"1.3.0-alpha.112"
        return Base.mktempdir(args...; kwargs...)
    else
        return Base.mktempdir(args...; cleanup=false, kwargs...)
    end
end

function render(doc::Documents.Document, settings::LaTeX=LaTeX())
    @info "LaTeXWriter: creating the LaTeX file."
    Base.mktempdir() do path
        cp(joinpath(doc.user.root, doc.user.build), joinpath(path, "build"))
        cd(joinpath(path, "build")) do
            name = doc.user.sitename
            let tag = get(ENV, "TRAVIS_TAG", "")
                if occursin(Base.VERSION_REGEX, tag)
                    v = VersionNumber(tag)
                    name *= "-$(v.major).$(v.minor).$(v.patch)"
                end
            end
            name = replace(name, " " => "")
            texfile = name * ".tex"
            pdffile = name * ".pdf"
            open(texfile, "w") do io
                context = Context(io)
                writeheader(context, doc)
                for (title, filename, depth) in files(doc.user.pages)
                    context.filename = filename
                    empty!(context.footnotes)
                    if 1 <= depth <= length(DOCUMENT_STRUCTURE)
                        header_type = DOCUMENT_STRUCTURE[depth]
                        header_text = "\n\\$(header_type){$(title)}\n"
                        if isempty(filename)
                            _println(context, header_text)
                        else
                            path = normpath(filename)
                            page = doc.blueprint.pages[path]
                            if get(page.globals.meta, :IgnorePage, :none) !== :latex
                                context.depth = depth + (isempty(title) ? 0 : 1)
                                context.depth > depth && _println(context, header_text)
                                latex(context, page, doc)
                            end
                        end
                    end
                end
                writefooter(context, doc)
            end
            cp(STYLE, "documenter.sty")

            # compile .tex
            status = compile_tex(doc, settings, texfile)

            # Debug: if DOCUMENTER_LATEX_DEBUG environment variable is set, copy the LaTeX
            # source files over to a directory under doc.user.root.
            if haskey(ENV, "DOCUMENTER_LATEX_DEBUG")
                dst = isempty(ENV["DOCUMENTER_LATEX_DEBUG"]) ? mktempdir(doc.user.root) :
                    joinpath(doc.user.root, ENV["DOCUMENTER_LATEX_DEBUG"])
                sources = cp(pwd(), dst, force=true)
                @info "LaTeX sources copied for debugging to $(sources)"
            end

            # If the build was successful, copy the PDF or the LaTeX source to the .build directory
            if status && (settings.platform != "none")
                cp(pdffile, joinpath(doc.user.root, doc.user.build, pdffile); force = true)
            elseif status && (settings.platform == "none")
                cp(pwd(), joinpath(doc.user.root, doc.user.build); force = true)
            else
                error("Compiling the .tex file failed. See logs for more information.")
            end
        end
    end
end

const DOCKER_IMAGE_TAG = "0.1"

function compile_tex(doc::Documents.Document, settings::LaTeX, texfile::String)
    if settings.platform == "native"
        Sys.which("latexmk") === nothing && (@error "LaTeXWriter: latexmk command not found."; return false)
        @info "LaTeXWriter: using latexmk to compile tex."
        try
            piperun(`latexmk -f -interaction=nonstopmode -view=none -lualatex -shell-escape $texfile`)
            return true
        catch err
            logs = cp(pwd(), mktempdir(); force=true)
            @error "LaTeXWriter: failed to compile tex with latexmk. " *
                   "Logs and partial output can be found in $(Utilities.locrepr(logs))." exception = err
            return false
        end
    elseif settings.platform == "docker"
        Sys.which("docker") === nothing && (@error "LaTeXWriter: docker command not found."; return false)
        @info "LaTeXWriter: using docker to compile tex."
        script = """
            mkdir /home/zeptodoctor/build
            cd /home/zeptodoctor/build
            cp -r /mnt/. .
            latexmk -f -interaction=nonstopmode -view=none -lualatex -shell-escape $texfile
            """
        try
            piperun(`docker run -itd -u zeptodoctor --name latex-container -v $(pwd()):/mnt/ --rm juliadocs/documenter-latex:$(DOCKER_IMAGE_TAG)`)
            piperun(`docker exec -u zeptodoctor latex-container bash -c $(script)`)
            piperun(`docker cp latex-container:/home/zeptodoctor/build/. .`)
            return true
        catch err
            logs = cp(pwd(), mktempdir(); force=true)
            @error "LaTeXWriter: failed to compile tex with docker. " *
                   "Logs and partial output can be found in $(Utilities.locrepr(logs))." exception = err
            return false
        finally
            try; piperun(`docker stop latex-container`); catch; end
        end
    elseif settings.platform == "none"
        @info "Skipping compiling tex file."
        return true
    end
end

function piperun(cmd)
    verbose = "--verbose" in ARGS || get(ENV, "DOCUMENTER_VERBOSE", "false") == "true"
    run(pipeline(cmd, stdout = verbose ? stdout : "LaTeXWriter.stdout",
                      stderr = verbose ? stderr : "LaTeXWriter.stderr"))
end

function writeheader(io::IO, doc::Documents.Document)
    custom = joinpath(doc.user.root, doc.user.source, "assets", "custom.sty")
    isfile(custom) ? cp(custom, "custom.sty"; force = true) : touch("custom.sty")
    preamble =
        """
        \\documentclass{memoir}

        \\usepackage{./documenter}
        \\usepackage{./custom}

        \\title{
            {\\HUGE $(doc.user.sitename)}\\\\
            {\\Large $(get(ENV, "TRAVIS_TAG", ""))}
        }
        \\author{$(doc.user.authors)}

        \\begin{document}

        \\frontmatter
        \\maketitle
        \\tableofcontents

        \\mainmatter

        """
    _println(io, preamble)
end

function writefooter(io::IO, doc::Documents.Document)
    _println(io, "\n\\end{document}")
end

function latex(io::IO, page::Documents.Page, doc::Documents.Document)
    for element in page.elements
        latex(io, page.mapping[element], page, doc)
    end
end

function latex(io::IO, vec::Vector, page, doc)
    for each in vec
        latex(io, each, page, doc)
    end
end

function latex(io::IO, anchor::Anchors.Anchor, page, doc)
    id = string(hash(string(anchor.id, "-", anchor.nth)))
    _println(io, "\n\\hypertarget{", id, "}{}\n")
    latex(io, anchor.object, page, doc)
end

## Documentation Nodes.

function latex(io::IO, node::Documents.DocsNodes, page, doc)
    for node in node.nodes
        latex(io, node, page, doc)
    end
end

function latex(io::IO, node::Documents.DocsNode, page, doc)
    id = string(hash(string(node.anchor.id)))
    # Docstring header based on the name of the binding and it's category.
    _println(io, "\\hypertarget{", id, "}{} ")
    _print(io, "\\hyperlink{", id, "}{\\texttt{")
    latexesc(io, string(node.object.binding))
    _print(io, "}} ")
    _println(io, " -- {", Utilities.doccat(node.object), ".}\n")
    # # Body. May contain several concatenated docstrings.
    _println(io, "\\begin{adjustwidth}{2em}{0pt}")
    latexdoc(io, node.docstr, page, doc)
    _println(io, "\n\\end{adjustwidth}")
end

function latexdoc(io::IO, md::Markdown.MD, page, doc)
    if haskey(md.meta, :results)
        # The `:results` field contains a vector of `Docs.DocStr` objects associated with
        # each markdown object. The `DocStr` contains data such as file and line info that
        # we need for generating correct scurce links.
        for (markdown, result) in zip(md.content, md.meta[:results])
            latex(io, Writers.MarkdownWriter.dropheaders(markdown), page, doc)
            # When a source link is available then print the link.
            url = Utilities.url(doc.internal.remote, doc.user.repo, result)
            if url !== nothing
                link = "\\href{$url}{\\texttt{source}}"
                _println(io, "\n", link, "\n")
            end
        end
    else
        # Docstrings with no `:results` metadata won't contain source locations so we don't
        # try to print them out. Just print the basic docstring.
        render(io, mime, dropheaders(md), page, doc)
    end
end

function latexdoc(io::IO, other, page, doc)
    # TODO: properly support non-markdown docstrings at some point.
    latex(io, other, page, doc)
end

## Index, Contents, and Eval Nodes.

function latex(io::IO, index::Documents.IndexNode, page, doc)
    # Having an empty itemize block in LaTeX throws an error, so we bail early
    # in that situation:
    isempty(index.elements) && (_println(io); return)

    _println(io, "\\begin{itemize}")
    for (object, _, page, mod, cat) in index.elements
        id = string(hash(string(Utilities.slugify(object))))
        text = string(object.binding)
        _print(io, "\\item \\hyperlink{")
        _print(io, id, "}{\\texttt{")
        latexesc(io, text)
        _println(io, "}}")
    end
    _println(io, "\\end{itemize}\n")
end

function latex(io::IO, contents::Documents.ContentsNode, page, doc)
    # Having an empty itemize block in LaTeX throws an error, so we bail early
    # in that situation:
    isempty(contents.elements) && (_println(io); return)

    depth = 1
    _println(io, "\\begin{itemize}")
    for (count, path, anchor) in contents.elements
        header = anchor.object
        level = Utilities.header_level(header)
        id = string(hash(string(anchor.id, "-", anchor.nth)))
        # If we're changing depth, we need to make sure we always print the
        # correct number of \begin{itemize} and \end{itemize} statements.
        if level > depth
            for k in 1:(level - depth)
                # if we jump by more than one level deeper we need to put empty
                # \items in -- otherwise LaTeX will complain
                (k >= 2) && _println(io, "\\item ~")
                _println(io, "\\begin{itemize}")
                depth += 1
            end
        elseif level < depth
            for _ in 1:(depth - level)
                _println(io, "\\end{itemize}")
                depth -= 1
            end
        end
        # Print the corresponding \item statement
        _print(io, "\\item \\hyperlink{", id, "}{")
        latexinline(io, header.text)
        _println(io, "}")
    end
    # print any remaining missing \end{itemize} statements
    for _ = 1:depth; _println(io, "\\end{itemize}"); end
    _println(io)
end

function latex(io::IO, node::Documents.EvalNode, page, doc)
    node.result === nothing ? nothing : latex(io, node.result, page, doc)
end

# Select the "best" representation for LaTeX output.
using Base64: base64decode
function latex(io::IO, mo::Documents.MultiOutput)
    foreach(x->Base.invokelatest(latex, io, x), mo.content)
end
function latex(io::IO, d::Dict{MIME,Any})
    filename = String(rand('a':'z', 7))
    if haskey(d, MIME"image/png"())
        write("$(filename).png", base64decode(d[MIME"image/png"()]))
        _println(io, """
        \\begin{figure}[H]
        \\centering
        \\includegraphics[max width=\\linewidth]{$(filename)}
        \\end{figure}
        """)
    elseif haskey(d, MIME"image/jpeg"())
        write("$(filename).jpeg", base64decode(d[MIME"image/jpeg"()]))
        _println(io, """
        \\begin{figure}[H]
        \\centering
        \\includegraphics[max width=\\linewidth]{$(filename)}
        \\end{figure}
        """)
    elseif haskey(d, MIME"text/latex"())
        latex(io, Utilities.mdparse(d[MIME"text/latex"()]; mode = :single))
    elseif haskey(d, MIME"text/markdown"())
        latex(io, Markdown.parse(d[MIME"text/markdown"()]))
    elseif haskey(d, MIME"text/plain"())
        latex(io, Markdown.Code(d[MIME"text/plain"()]))
    else
        error("this should never happen.")
    end
    return nothing
end

## Basic Nodes. AKA: any other content that hasn't been handled yet.

latex(io::IO, str::AbstractString, page, doc) = _print(io, str)

function latex(io::IO, other, page, doc)
    _println(io)
    latex(io, other)
    _println(io)
end

latex(io::IO, md::Markdown.MD) = latex(io, md.content)

function latex(io::IO, content::Vector)
    for c in content
        latex(io, c)
    end
end

function latex(io::IO, h::Markdown.Header{N}) where N
    tag = DOCUMENT_STRUCTURE[min(io.depth + N - 1, length(DOCUMENT_STRUCTURE))]
    _print(io, "\\", tag, "{")
    io.in_header = true
    latexinline(io, h.text)
    io.in_header = false
    _println(io, "}\n")
end

# Whitelisted lexers.
const LEXER = Set([
    "julia",
    "jlcon",
])

function latex(io::IO, code::Markdown.Code)
    language = isempty(code.language) ? "none" : code.language
    # the julia-repl is called "jlcon" in Pygments
    language = (language == "julia-repl") ? "jlcon" : language
    escape = '⊻' ∈ code.code
    if language in LEXER
        _print(io, "\n\\begin{minted}")
        if escape
            _print(io, "[escapeinside=\\#\\%]")
        end
        _println(io, "{", language, "}")
        if escape
            _print_code_escapes_minted(io, code.code)
        else
            _print(io, code.code)
        end
        _println(io, "\n\\end{minted}\n")
    else
        _print(io, "\n\\begin{lstlisting}")
        if escape
            _println(io, "[escapeinside=\\%\\%]")
            _print_code_escapes_lstlisting(io, code.code)
        else
            _println(io)
            _print(io, code.code)
        end
        _println(io, "\n\\end{lstlisting}\n")
    end
end

function _print_code_escapes_minted(io, s::AbstractString)
    for ch in s
        ch === '#' ? _print(io, "##%") :
        ch === '%' ? _print(io, "#%%") : # Note: "#\\%%" results in pygmentize error...
        ch === '⊻' ? _print(io, "#\\unicodeveebar%") :
                     _print(io, ch)
    end
end
function _print_code_escapes_lstlisting(io, s::AbstractString)
    for ch in s
        ch === '%' ? _print(io, "%\\%%") :
        ch === '⊻' ? _print(io, "%\\unicodeveebar%") :
                     _print(io, ch)
    end
end

function latexinline(io::IO, code::Markdown.Code)
    _print(io, "\\texttt{")
    _print_code_escapes_inline(io, code.code)
    _print(io, "}")
end

function _print_code_escapes_inline(io, s::AbstractString)
    for ch in s
        ch === '⊻' ? _print(io, "\\unicodeveebar{}") :
                     latexesc(io, ch)
    end
end

function latex(io::IO, md::Markdown.Paragraph)
    for md in md.content
        latexinline(io, md)
    end
    _println(io, "\n")
end

function latex(io::IO, md::Markdown.BlockQuote)
    wrapblock(io, "quote") do
        latex(io, md.content)
    end
end

function latex(io::IO, md::Markdown.Admonition)
    wrapblock(io, "quote") do
        wrapinline(io, "textbf") do
            latexinline(io, md.title)
        end
        _println(io, "\n")
        latex(io, md.content)
    end
end

function latex(io::IO, f::Markdown.Footnote)
    id = get(io.footnotes, f.id, 1)
    _print(io, "\\footnotetext[", id, "]{")
    latex(io, f.text)
    _println(io, "}")
end

function latex(io::IO, md::Markdown.List)
    # `\begin{itemize}` is used here for both ordered and unordered lists since providing
    # custom starting numbers for enumerated lists is simpler to do by manually assigning
    # each number to `\item` ourselves rather than using `\setcounter{enumi}{<start>}`.
    #
    # For an ordered list starting at 5 the following will be generated:
    #
    # \begin{itemize}
    #   \item[5. ] ...
    #   \item[6. ] ...
    #   ...
    # \end{itemize}
    #
    pad = ndigits(md.ordered + length(md.items)) + 2
    fmt = n -> (Markdown.isordered(md) ? "[$(rpad("$(n + md.ordered - 1).", pad))]" : "")
    wrapblock(io, "itemize") do
        for (n, item) in enumerate(md.items)
            _print(io, "\\item$(fmt(n)) ")
            latex(io, item)
            n < length(md.items) && _println(io)
        end
    end
end

function latex(io::IO, hr::Markdown.HorizontalRule)
    _println(io, "{\\rule{\\textwidth}{1pt}}")
end

# This (equation*, split) math env seems to be the only way to correctly render all the
# equations in the Julia manual. However, if the equation is already wrapped in
# align/align*, then there is no need to further wrap it (in fact, it will break).
function latex(io::IO, math::Markdown.LaTeX)
    if occursin(r"^\\begin\{align\*?\}", math.formula)
        _print(io, math.formula)
    else
        _print(io, "\\begin{equation*}\n\\begin{split}")
        _print(io, math.formula)
        _println(io, "\\end{split}\\end{equation*}")
    end
end

function latex(io::IO, md::Markdown.Table)
    _println(io, "\n\\begin{table}[h]")
    _print(io, "\n\\begin{tabulary}{\\linewidth}")
    _println(io, "{|", uppercase(join(md.align, '|')), "|}")
    for (i, row) in enumerate(md.rows)
        i === 1 && _println(io, "\\hline")
        for (j, cell) in enumerate(row)
            j === 1 || _print(io, " & ")
            latexinline(io, cell)
        end
        _println(io, " \\\\")
        _println(io, "\\hline")
    end
    _println(io, "\\end{tabulary}\n")
    _println(io, "\\end{table}\n")
end

function latex(io::IO, raw::Documents.RawNode)
    raw.name === :latex ? _println(io, "\n", raw.text, "\n") : nothing
end

# Inline Elements.

function latexinline(io::IO, md::Vector)
    for c in md
        latexinline(io, c)
    end
end

function latexinline(io::IO, md::AbstractString)
    latexesc(io, md)
end

function latexinline(io::IO, md::Markdown.Bold)
    wrapinline(io, "textbf") do
        latexinline(io, md.text)
    end
end

function latexinline(io::IO, md::Markdown.Italic)
    wrapinline(io, "emph") do
        latexinline(io, md.text)
    end
end

function latexinline(io::IO, md::Markdown.Image)
    wrapblock(io, "figure") do
        _println(io, "\\centering")
        url = if Utilities.isabsurl(md.url)
            @warn "images with absolute URLs not supported in LaTeX output in $(Utilities.locrepr(io.filename))" url = md.url
            # We nevertheless output an \includegraphics with the URL. The LaTeX build will
            # then give an error, indicating to the user that something wrong. Only the
            # warning would be drowned by all the output from LaTeX.
            md.url
        elseif startswith(md.url, '/')
            # URLs starting with a / are assumed to be relative to the document's root
            normpath(lstrip(md.url, '/'))
        else
            normpath(joinpath(dirname(io.filename), md.url))
        end
        url = replace(url, "\\" => "/") # use / on Windows too.
        wrapinline(io, "includegraphics[max width=\\linewidth]") do
            _print(io, url)
        end
        _println(io)
        wrapinline(io, "caption") do
            latexinline(io, md.alt)
        end
        _println(io)
    end
end

function latexinline(io::IO, f::Markdown.Footnote)
    id = get!(io.footnotes, f.id, length(io.footnotes) + 1)
    _print(io, "\\footnotemark[", id, "]")
end

function latexinline(io::IO, md::Markdown.Link)
    if io.in_header
        latexinline(io, md.text)
    else
        if occursin(".md#", md.url)
            file, target = split(md.url, ".md#"; limit = 2)
            id = string(hash(target))
            wrapinline(io, "hyperlink") do
                _print(io, id)
            end
            _print(io, "{")
            latexinline(io, md.text)
            _print(io, "}")
        else
            wrapinline(io, "href") do
                latexesc(io, md.url)
            end
            _print(io, "{")
            latexinline(io, md.text)
            _print(io, "}")
        end
    end
end

function latexinline(io, math::Markdown.LaTeX)
    # Handle MathJax and TeX inconsistency since the first wants `\LaTeX` wrapped
    # in math delims, whereas actual TeX fails when that is done.
    math.formula == "\\LaTeX" ? _print(io, math.formula) : _print(io, "\\(", math.formula, "\\)")
end

function latexinline(io, hr::Markdown.HorizontalRule)
    _println(io, "\\rule{\\textwidth}{1pt}}")
end

# Metadata Nodes get dropped from the final output for every format but are needed throughout
# rest of the build and so we just leave them in place and print a blank line in their place.
latex(io::IO, node::Documents.MetaNode, page, doc) = _println(io, "\n")

# Utilities.

const _latexescape_chars = Dict{Char, AbstractString}(
    '~' => "{\\textasciitilde}",
    '^' => "{\\textasciicircum}",
    '\\' => "{\\textbackslash}",
    '\'' => "{\\textquotesingle}",
    '"' => "{\\textquotedbl}",
)
for ch in "&%\$#_{}"
    _latexescape_chars[ch] = "\\$ch"
end

latexesc(io, ch::AbstractChar) = _print(io, get(_latexescape_chars, ch, ch))

function latexesc(io, s::AbstractString)
    for ch in s
        latexesc(io, ch)
    end
end

latexesc(s) = sprint(latexesc, s)

function wrapblock(f, io, env)
    _println(io, "\\begin{", env, "}")
    f()
    _println(io, "\\end{", env, "}")
end

function wrapinline(f, io, cmd)
    _print(io, "\\", cmd, "{")
    f()
    _print(io, "}")
end

function files!(out::Vector, v::Vector, depth)
    for each in v
        files!(out, each, depth + 1)
    end
    return out
end

# Tuples come from `hide(page)` with either
# (visible, nothing,    page,         children) or
# (visible, page.first, pages.second, children)
function files!(out::Vector, v::Tuple, depth)
    files!(out, v[2] == nothing ? v[3] : v[2] => v[3], depth)
    files!(out, v[4], depth)
end

files!(out, s::AbstractString, depth) = push!(out, ("", s, depth))

function files!(out, p::Pair{S, T}, depth) where {S <: AbstractString, T <: AbstractString}
    push!(out, (p.first, p.second, depth))
end

function files!(out, p::Pair{S, V}, depth) where {S <: AbstractString, V}
    push!(out, (p.first, "", depth))
    files!(out, p.second, depth)
end

files(v::Vector) = files!(Tuple{String, String, Int}[], v, 0)

end


text:
# https://julianonconvex.github.io/Nonconvex.jl/stable/
# concerning that Model and DictModel find different solutions

using Nonconvex
Nonconvex.@load Ipopt

f(x) = sqrt(x[2])
g(x, a, b) = (a*x[1] + b)^3 - x[2]

model = Model(f)

addvar!(model, 0.0, 10.0, init=1.0)
addvar!(model, 0.0, 10.0, init=1.0)

add_ineq_constraint!(model, x -> g(x, 2, 0))
add_ineq_constraint!(model, x -> g(x, -1, 1))

x0 = NonconvexCore.getinit(model)
r = optimize(model, IpoptAlg(), x0)

println(r.minimum)
println(r.minimizer)

f1(x) = sqrt(x["2"])
g1(x, a, b) = (a*x["1"] + b)^3 - x["2"]

model = DictModel(f1)

addvar!(model, "1", 0.0, 10.0, init=1.0)
addvar!(model, "2", 0.0, 10.0, init=1.0)

add_ineq_constraint!(model, x -> g1(x, 2, 0))
add_ineq_constraint!(model, x -> g1(x, -1, 1))

x0 = NonconvexCore.getinit(model)
r = optimize(model, IpoptAlg(), x0)

println(r.minimum)
println(r.minimizer)

ext: jl

max_stars_repo_name: lanl-ansi/rosetta-opf

max_issues_repo_name: lanl-ansi/rosetta-opf

max_forks_repo_name: lanl-ansi/rosetta-opf


text:
"""
    struct OLS
        X::Array{Float64,2}
        y::Array{Float64,1}
        betas::Array{Float64,1}
    end

    Immutable data structure that holds design matrix, response vector, and estimated regression parameters.

# Arguments
- `X::Array{Float64, 2}`: Design matrix.
- `y::Array{Float64, 1}`: Response vector.
- `betas::Array{Float64, 1}`: Regression coefficients.

See also: [`ols`](@ref), [`wls`](@ref)
"""
struct OLS
    X::Array{Float64,2}
    y::Array{Float64,1}
    betas::Array{Float64,1}
end

"""
    ols(X, y)

    Create OLS object with estimated regression coefficients.

# Arguments
- `X::Array{Float64, 2}`: Design matrix.
- `y::Array{Float64, 1}`: Response vector.

# Examples
```julia-repl
julia> X = hcat(ones(24), phones[:,"year"]);
julia> y = phones[:,"calls"];
julia> reg = ols(X, y)
julia> reg.betas
2-element Array{Float64,1}:
 -260.0592463768119
    5.04147826086957
```

See also: [`OLS`](@ref), [`wls`](@ref)

"""
ols(X::Array{Float64,2}, y::Array{Float64,1})::OLS = OLS(X, y, qr(X, Val(true)) \ y)

"""
    wls(X, y, wts)

    Estimate weighted least squares regression and create OLS object with estimated parameters.

# Arguments
- `X::Array{Float64, 2}`: Design matrix.
- `y::Array{Float64, 1}`: Response vector.
- `wts::Array{Float64, 1}`: Weights vector.

# Examples
```julia-repl
julia> X = hcat(ones(24), phones[:,"year"]);
julia> y = phones[:,"calls"];
julia> w = ones(24)
julia> w[15:20] .= 0.0
julia> reg = wls(X, y, w)
julia> reg.betas
2-element Array{Float64,1}:
 -63.481644325290425
   1.3040571939231453
```

See also: [`ols`](@ref), [`OLS`](@ref)

"""
function wls(X::Array{Float64,2}, y::Array{Float64,1}, wts::Array{Float64,1})
    W = Diagonal(sqrt.(wts))
    #  I commented this because passing weighted values of X and y to OLS
    #  causes wrong calculations of residuals.
    #  return ols(W * X, W * y)
    return OLS(X, y, qr(W * X, Val(true)) \ (W * y))
end

"""
    residuals(ols)

Estimate weighted least squares regression and create OLS object with estimated parameters.

# Arguments
- `ols::OLS`: OLS object, possible created using `ols` or `wls`.

See also: [`ols`](@ref), [`wls`](@ref)
"""
residuals(ols::OLS) = ols.y .- ols.X * ols.betas

"""
    coef(ols)

Extract regression coefficients from an `OLS` object.

# Arguments
- `ols::OLS`: OLS object, possible created using `ols` or `wls`.

See also: [`ols`](@ref), [`residuals`](@ref), [`predict`](@ref)
"""
coef(ols::OLS) = ols.betas

"""
    predict(ols)

Calculate estimated response using an `OLS` object.

# Arguments
- `ols::OLS`: OLS object, possible created using `ols` or `wls`.

See also: [`ols`](@ref), [`residuals`](@ref)
"""
predict(ols::OLS) = ols.X * ols.betas

predict(ols::OLS, X::Array{Float64,2}) = X * ols.betas

max_stars_repo_name: angeris/LinRegOutliers

max_issues_repo_name: angeris/LinRegOutliers

max_forks_repo_name: angeris/LinRegOutliers


text:
"""
    setColNames(str)
Call this function to update column names

remove blank spaces, dollar sign, or Hash chars

make all columns uppercase and replace hyphen with underscore
"""
function setColNames(str)
	return uppercase(reduce(replace, ["-" => "_"," " => "_", "\$" => "_","#" => "_","%" => "_","(" => "_",")" => "_"], init=str))
end

"""
    getDBConnection(credFilePath:: AbstractString)
Call this function to read database credentials

database credentials are stored in environment.txt file in following format.

`environment.txt`

    user=username
    pwd=password
    dsn=userdsnname
    hive=hdinsightstr
    port=portnumber

# Example
```julia-repl
julia> fl = getDBConnection("environment.txt")
```
"""
function getDBConnection(credFilePath::AbstractString)
	return "access db credentials from conn object. conn.user, conn.pwd... "
end

"""
    getSQLs(credFilePath:: AbstractString)
Call this function to read sqls for a given table from txt file.
SQLs are stored in txt file in this format.

`sqls.txt`

    createTable1=INSERT INTO table1 (column1) SELECT table2.column1 FROM table2 WHERE table2.column1 > 100;
    readTable1=SELECT * FROM table1;
    updateTable1=UPDATE table SET column1 = value1, column2 = value2, ... WHERE condition;
    upsertTable1=BEDIN tran IF EXISTS (SELECT * FROM table1 WITH (updlock,serializable) WHERE key = @key) BEGIN UDPATE table1 SET ... WHERE key = @key END ELSE BEGIN INSERT INTO table1 (key, ...) VALUES (@key, ...) END COMMIT TRAN
    softDeleteTable1=UPDATE table SET deleted=True, ... WHERE condition;
    hardDeleteTable1=delete * from table1 where table1.column1 > 100

# Example
```julia-repl
julia> fl = getSQLs("environment.txt")
```
"""
function getSQLs(credFilePath::AbstractString)
	return "access sqls from sq object. sq.createTable1, sq.upsertTable1 ... "
end

"""
    getDSNs()
Call this function to read computing machine available DSNs

# Example
```julia-repl
julia> dsns = getDSNs()
```
"""
function getDSNs()
    # ODBC.dsns()
	return "available DSNs on your machine."
end

"""
    getDrivers()
Call this function to read computing machine available Drivers

# Example
```julia-repl
julia> drivers = getDrivers()
```
"""
function getDrivers()
    # ODBC.drivers()
	return "available Drivers on your machine."
end

"""
    runSQL(conn:: AbstractString, sql:: AbstractString)
Call this function to run sql in database

# Example
```julia-repl
julia> res = runSQL(conn, SQLStatement)
```
"""
function runSQL(conn::AbstractString, sql::AbstractString)
    # DBInterface.execute(conn, "SQL statament") |> DataFrame
	return "returns results of SQL statement "
end

"""
    setCloseConnection()
Call this function to close database connection

# Example
```julia-repl
julia> cls = setCloseConnection()
```
"""
function setCloseConnection()
    # DBInterface.close!(conn)
	return "DB Connection is closed now."
end

"""
    getXLSinDirectory(str)
Call this function to read all xls files inside directory

# Example
```julia-repl
julia> files = getXLSinDirectory("dirname")
```
"""
function getXLSinDirectory(dirName)
	tmp = readdir(dirName)
	tmp_xls = [];
	for i in 1:length(tmp)
		push!(tmp_xls, XLSX.readxlsx(joinpath(dirName, tmp[i])))
	end
	return tmp_xls
end

"""
getArrangedWords(str)
Call this function to remove duplicates, unwanted symbols and return uppercase unique values

# Example
```julia-repl
julia> wd = getArrangedWords("Amit ; Shukla SHUKLA Shukle , . AmIT Amit # Shuklam Amit ,")
```
"""
function getArrangedWords(str)
    sd = StringDocument(str)
    prepare!(sd, strip_articles | strip_numbers | strip_html_tags | strip_pronouns)
    t = TextAnalysis.tokens(sd)
    # df_text!(uppercase.(df_text.t), :)
    # return uppercase(str), unique(sort!(uppercase.(t))), typeof(t),
    return filter!(x -> x ∉ [",","#", ";","."], unique(sort!(uppercase.(t))))
end

"""
getFuzzyWuzzy(str)
Call this function to find closest match for a given string in data frame column lookup

First parameter is the search string, second is the DataFrame followed by columnname in DataFrame which needs to be searched

# Example
```julia-repl
julia> df_dname = DataFrame(name=["John Doe", "Jen Doe","MICHAEL Doe", "Jacob Doe", "Julia Dpe", "Michael Jackson"],age=[35,26,35,10,5,45])
julia> wd = getFuzzyWuzzy("Mike Jackson", df_dname, "name")
```
"""
function getFuzzyWuzzy(str, lookupDF, colName)
    score = 0.0;
    lookupDF[!,:score] = zeros(size(lookupDF, 1))
    for i in 1:size(lookupDF, 1)
        # lookupDF[!,:score][i] = TokenMax(RatcliffObershelp())(uppercase(str),
        # 	uppercase(lookupDF[:,Symbol(colName)][i])
        lookupDF[!,:score][i] = TokenSet(RatcliffObershelp())(uppercase(str),
            uppercase(lookupDF[:,Symbol(colName)][i])
        )
    end
    # return lookupDF, minimum(lookupDF.score)
    return unique(DataFrame([lookupDF[findmin(lookupDF.score)[2],:] for g in groupby(lookupDF, Symbol(colName))]))
end

"""
setReplaceText(str)
Call this function to find and replace word in text string

First parameter is the term to be replaced, next is the term replaced with followed by text string where text is searched and replaced.

# Example
```julia-repl
julia> wd = setReplaceText("Shkla","Shukla ","Amit Shkla Los Angeles")
```
"""
function setReplaceText(term, replaceWith, srchInTxt)
    return repeat(replaceWith, length(collect(m for m in eachmatch(Regex("$term"), 	srchInTxt))))
end

"""
setRemoveText(str)
Call this function to find and remove word in text string

First parameter is the complete string text, next parameter is list of words to be removed.

# Example
```julia-repl
julia> wd = setRemoveText("Amit Shukla Shkla Los Angel Angeles", ["Shkla", "Angel"])
```
"""
function setRemoveText(str, rem_words)
    # sd = StringDocument(str)
    # remove_words!(sd, rem_words)
    # return String(strip(TextAnalysis.text(sd)))
    return "Amit Shukla Los Angeles"
end

# token functions
const regexp = Dict(
	:mention => r"@\w+",
	:hashtag => r"#\w+",
	:url => r"http[s]?://(?:[a-zA-Z]|[0-9]|[$-_@.&+]|[!*\(\),]|(?:%[0-9a-fA-F][0-9a-fA-F]))+"
);

"""
getTokens(s, token_type)
Call this function to extract tokens from string (example - extract urls)

First parameter is the complete string text, next parameter is list of words to be removed.

# Example
```julia-repl
julia> wd = getTokens("https://yahoo.com is the Yahoo website url", url)
```
"""
function getTokens(s, token_type)
	# return collect(x.match for x in eachmatch(regexp[token_type], s))
    return "https://yahoo.com"
end

"""
setRemoveTokens(str)
Call this function to remove tokens

# Example
```julia-repl
julia> wd = setRemoveTokens("Amit Shukla Shkla Los Angel Angeles")
```
"""
function setRemoveTokens(s)
	# for re in values(regexp)
	# 	s = replace(s, re => "")
	# end
	# return s
    return "Amit Shukla Los Angeles"
end

"""
getDuplicateRows(df_dname:: DataFrame, colNames:: Vector)
Call this function to find duplicates in a data frame column based on columnnames (key columns)

First parameter is the DataFrame, followed by all columnnames in DataFrame which are key columns in a dataframe.

Function returns dataframe row indexes which are duplicates based on key columns provided.

# Example
```julia-repl
julia> df_dname = DataFrame(name=["John Doe", "John Doe","MICHAEL Doe", "Jacob Doe", "Julia Dpe", "Michael Jackson"],age=[35,35,35,10,5,45], state=["CA","CA","CA","CA","CO","CA"])
julia> dup = getDuplicateRows(df_dname, ["name","age"])
```
"""
function getDuplicateRows(lookupDF, colNames)
    # score = 0.0;
    # lookupDF[!,:score] = zeros(size(lookupDF, 1))
    # for i in 1:size(lookupDF, 1)
    #     # lookupDF[!,:score][i] = TokenMax(RatcliffObershelp())(uppercase(str),
    #     # 	uppercase(lookupDF[:,Symbol(colName)][i])
    #     lookupDF[!,:score][i] = TokenSet(RatcliffObershelp())(uppercase(str),
    #         uppercase(lookupDF[:,Symbol(colName)][i])
    #     )
    # end
    # # return lookupDF, minimum(lookupDF.score)
    # return unique(DataFrame([lookupDF[findmin(lookupDF.score)[2],:] for g in groupby(lookupDF, Symbol(colName))]))
    return (1, 2)
end

"""
getKeyColumns(df_dname:: DataFrame)
Call this function to find key columns in a data frame

First parameter is the DataFrame

Function returns dataframe columns indexes which are key columns.

# Example
```julia-repl
julia> df_dname = DataFrame(name=["John Doe", "John Doe","MICHAEL Doe", "Jacob Doe", "Julia Dpe", "Michael Jackson"],age=[35,35,35,10,5,45], state=["CA","CA","CA","CA","CO","CA"])
julia> kcols = getKeyColumns(df_dname)
```
"""
function getKeyColumns(df_dname)
    # score = 0.0;
    # lookupDF[!,:score] = zeros(size(lookupDF, 1))
    # for i in 1:size(lookupDF, 1)
    #     # lookupDF[!,:score][i] = TokenMax(RatcliffObershelp())(uppercase(str),
    #     # 	uppercase(lookupDF[:,Symbol(colName)][i])
    #     lookupDF[!,:score][i] = TokenSet(RatcliffObershelp())(uppercase(str),
    #         uppercase(lookupDF[:,Symbol(colName)][i])
    #     )
    # end
    # # return lookupDF, minimum(lookupDF.score)
    # return unique(DataFrame([lookupDF[findmin(lookupDF.score)[2],:] for g in groupby(lookupDF, Symbol(colName))]))
    return (1, 2, 3)
end

"""
setRemDuplicateRows(df_dname:: DataFrame, colNames:: Vector)
Call this function to find & delete duplicates in a data frame column based on columnnames (key columns)

First parameter is the DataFrame, followed by all columnnames in DataFrame which are key columns in a dataframe.

Function returns dataframe after removing duplicates based on key columns provided.
in case of duplicates, it retains first row.

# Example
```julia-repl
julia> df_dname = DataFrame(name=["John Doe", "John Doe","MICHAEL Doe", "Jacob Doe", "Julia Dpe", "Michael Jackson"],age=[35,35,35,10,5,45], state=["CA","CA","CA","CA","CO","CA"])
julia> dup = setRemDuplicateRows(df_dname, ["name","age"])
```
"""
function setRemDuplicateRows(lookupDF, colNames)
    # score = 0.0;
    # lookupDF[!,:score] = zeros(size(lookupDF, 1))
    # for i in 1:size(lookupDF, 1)
    #     # lookupDF[!,:score][i] = TokenMax(RatcliffObershelp())(uppercase(str),
    #     # 	uppercase(lookupDF[:,Symbol(colName)][i])
    #     lookupDF[!,:score][i] = TokenSet(RatcliffObershelp())(uppercase(str),
    #         uppercase(lookupDF[:,Symbol(colName)][i])
    #     )
    # end
    # # return lookupDF, minimum(lookupDF.score)
    # return unique(DataFrame([lookupDF[findmin(lookupDF.score)[2],:] for g in groupby(lookupDF, Symbol(colName))]))
    return "dataframe after removing duplicates"
end

"""
getCategoryData(df_dname:: DataFrame, colName:: AbstractString, categoryRange)
Call this function to create a new column on DataFrame which provide a category based on ranges provided.

column must be contain only numerical values.

First parameter is the DataFrame, next is Column name for which categories are created, followed by Ranges.

Function returns dataframe with an extra columns with Categories like 1,2,3 and 0 for unmatched.

# Example
```julia-repl
julia> df_dname = DataFrame(name=["John Doe", "John Doe","MICHAEL Doe", "Jacob Doe", "Julia Dpe", "Michael Jackson"],age=[35,35,35,10,5,45], state=["CA","CA","CA","CA","CO","CA"])
julia> catData = getCategoryData(df_dname, "age", [0:10,10:20,20:30])
```
"""
function getCategoryData(lookupDF, colName)
    # score = 0.0;
    # lookupDF[!,:score] = zeros(size(lookupDF, 1))
    # for i in 1:size(lookupDF, 1)
    #     # lookupDF[!,:score][i] = TokenMax(RatcliffObershelp())(uppercase(str),
    #     # 	uppercase(lookupDF[:,Symbol(colName)][i])
    #     lookupDF[!,:score][i] = TokenSet(RatcliffObershelp())(uppercase(str),
    #         uppercase(lookupDF[:,Symbol(colName)][i])
    #     )
    # end
    # # return lookupDF, minimum(lookupDF.score)
    # return unique(DataFrame([lookupDF[findmin(lookupDF.score)[2],:] for g in groupby(lookupDF, Symbol(colName))]))
    return "original dataframe with a new category column added."
end

"""
getTreeData(df_dname:: DataFrame, colName:: AbstractString)
Call this function to create a flatten tree data structure.

First parameter is the DataFrame, next is Column name for which tree hierachies (flattened) are created.

Function returns a new dataframe with Levels.

# Example
```julia-repl
julia> df_dname = DataFrame(name=["John Doe", "John Doe","MICHAEL Doe", "Jacob Doe", "Julia Dpe", "Michael Jackson"],age=[35,35,35,10,5,45], state=["CA","CA","CA","CA","CO","CA"])
julia> catData = getTreeData(df_dname, "state")
```
"""
function getTreeData(lookupDF, colName)
    # score = 0.0;
    # lookupDF[!,:score] = zeros(size(lookupDF, 1))
    # for i in 1:size(lookupDF, 1)
    #     # lookupDF[!,:score][i] = TokenMax(RatcliffObershelp())(uppercase(str),
    #     # 	uppercase(lookupDF[:,Symbol(colName)][i])
    #     lookupDF[!,:score][i] = TokenSet(RatcliffObershelp())(uppercase(str),
    #         uppercase(lookupDF[:,Symbol(colName)][i])
    #     )
    # end
    # # return lookupDF, minimum(lookupDF.score)
    # return unique(DataFrame([lookupDF[findmin(lookupDF.score)[2],:] for g in groupby(lookupDF, Symbol(colName))]))
    return "original dataframe with levels, level 1 has two nodes CA, CO which contains name field as children."
end

"""
getMaskedData(df_dname:: DataFrame, colNames:: Vector)
Call this function to create a flatten tree data structure.

First parameter is the DataFrame, next is Column names which needs to be masked

Function returns a two dataframes, one with masked data and other with masked + original data.

# Example
```julia-repl
julia> df_dname = DataFrame(name=["John Doe", "John Doe","MICHAEL Doe", "Jacob Doe", "Julia Dpe", "Michael Jackson"],age=[35,35,35,10,5,45], state=["CA","CA","CA","CA","CO","CA"])
julia> catData = getMaskedData(df_dname, "name")
```
"""
function getMaskedData(lookupDF, colName)
    # score = 0.0;
    # lookupDF[!,:score] = zeros(size(lookupDF, 1))
    # for i in 1:size(lookupDF, 1)
    #     # lookupDF[!,:score][i] = TokenMax(RatcliffObershelp())(uppercase(str),
    #     # 	uppercase(lookupDF[:,Symbol(colName)][i])
    #     lookupDF[!,:score][i] = TokenSet(RatcliffObershelp())(uppercase(str),
    #         uppercase(lookupDF[:,Symbol(colName)][i])
    #     )
    # end
    # # return lookupDF, minimum(lookupDF.score)
    # return unique(DataFrame([lookupDF[findmin(lookupDF.score)[2],:] for g in groupby(lookupDF, Symbol(colName))]))
    return "one dataframe is with name field encrypted, other contains original and masked values both."
end


text:
# Analyze the blow-up time of numerical simulations of A 2D inviscid
# Kelvin-Helmholtz instability for different SBP discretizations
# implemented in Trixi.jl:
# - DGSEM with and without shock capturing
# - Gauss collocation
# - Polynomial and SBP approximations on triangles
#
# For optimal performance, run this script with
#   julia --threads=1 --check-bounds=no kelvin_helmholtz_instability.jl
# to let Julia use half of your `Sys.CPU_THREADS` cores. You can also specify
# the number of processes to use via
#   julia --threads=1 --procs=N --check-bounds=no kelvin_helmholtz_instability.jl
# where `N` is the number of processes you want to use.

# Package setup on one processor to avoid duplicating precompilation work etc.
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
Pkg.precompile()

using Distributed

# Add worker processes if none are present
if nprocs() == 1
  addprocs(Sys.CPU_THREADS ÷ 2)
end

# Package setup required on each processor
@everywhere begin
  using Pkg
  Pkg.activate(@__DIR__)
end
@everywhere using Trixi
@everywhere using OrdinaryDiffEq

# Setup on the driver process
using CSV, DataFrames

# Further setup required on each processor
@everywhere begin
  function set_initial_condition(atwood_number)
    A = atwood_number
    rho1 = 0.5 * one(A) # recover original with A = 3/7
    rho2 = rho1 * (1 + A) / (1 - A)
    function initial_condition_khi_nonsmooth(
        # domain size is [-1,+1]^2
        x, t, equations::CompressibleEulerEquations2D)

        # B is a discontinuous function with value 1 for -.5 <= x <= .5 and 0 elsewhere
        slope = 15
        B = 0.5 * (tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5))

        rho = rho1 + B * (rho2 - rho1)  # rho ∈ [rho_1, rho_2]
        v1 = B - 0.5                    # v1  ∈ [-.5, .5]
        v2 = 0.1 * sin(2 * pi * x[1])
        p = 1.0
        return prim2cons(SVector(rho, v1, v2, p), equations)
    end
    return initial_condition_khi_nonsmooth
  end

  function blow_up_time(atwood_number, args...; abstol=1.0e-7, reltol=1.0e-7)

    # 2D compressible Euler equations with ratio of specific heats 7/5
    equations = CompressibleEulerEquations2D(1.4)

    # dissipative setup using EC and ED fluxes
    surface_flux = flux_lax_friedrichs
    volume_flux  = flux_ranocha

    # use a try catch block since not all parameter combinations might work
    local mesh, solver
    try
      solver, mesh = setup_solver_mesh(args..., equations, surface_flux, volume_flux)
    catch e
      @info "solver setup faild" args e
      return NaN
    end

    semi = SemidiscretizationHyperbolic(
      mesh, equations, set_initial_condition(atwood_number), solver)

    tspan = (0.0, 10.0)
    ode = semidiscretize(semi, tspan)

    alive_callback = AliveCallback(alive_interval=500)

    @info "Initializing" atwood_number args

    alg = RDPK3SpFSAL49()
    integrator = init(ode, alg, abstol=abstol, reltol=reltol,
                      dt = 1.0e-4, # small initial dt
                      save_everystep=false, callback=alive_callback)

    try
      solve!(integrator)
    catch e
      @info "Blow-up" integrator.t args e
    end

    @info "Finished" atwood_number args

    return integrator.t
  end

  function setup_solver_mesh(::Type{DGSEM},
                             polydeg, refinement_level,
                             equations, surface_flux, volume_flux)

    solver = DGSEM(polydeg=polydeg, surface_flux=surface_flux,
                   volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

    coordinates_min = (-1.0, -1.0)
    coordinates_max = ( 1.0,  1.0)
    mesh = TreeMesh(coordinates_min, coordinates_max,
                    initial_refinement_level=refinement_level,
                    periodicity=(true, true),
                    n_cells_max=100_000)

    return solver, mesh
  end

  function setup_solver_mesh(::Type{VolumeIntegralShockCapturingHG},
                             polydeg, refinement_level,
                             equations, surface_flux, volume_flux)

    basis = LobattoLegendreBasis(polydeg)
    indicator_sc = IndicatorHennemannGassner(equations, basis,
                                             alpha_max=0.002,
                                             alpha_min=0.0001,
                                             alpha_smooth=true,
                                             variable=density_pressure)
    volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                     volume_flux_dg=volume_flux,
                                                     volume_flux_fv=surface_flux)
    solver = DGSEM(basis, surface_flux, volume_integral)

    coordinates_min = (-1.0, -1.0)
    coordinates_max = ( 1.0,  1.0)
    mesh = TreeMesh(coordinates_min, coordinates_max,
                    initial_refinement_level=refinement_level,
                    periodicity=(true, true),
                    n_cells_max=100_000)

    return solver, mesh
  end

  function setup_solver_mesh(::Type{GaussSBP},
                             polydeg, refinement_level,
                             equations, surface_flux, volume_flux)

    solver = DGMulti(polydeg=polydeg, element_type=Quad(),
                     approximation_type=GaussSBP(),
                     surface_integral=SurfaceIntegralWeakForm(surface_flux),
                     volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

    coordinates_min = (-1.0, -1.0)
    coordinates_max = ( 1.0,  1.0)
    cells_per_dimension = 2^refinement_level .* (1, 1)
    mesh = DGMultiMesh(solver;
                       coordinates_min, coordinates_max, cells_per_dimension,
                       periodicity=(true, true))

    return solver, mesh
  end

  function setup_solver_mesh(::Type{SBP},
                             polydeg, refinement_level,
                             equations, surface_flux, volume_flux)

    solver = DGMulti(polydeg=polydeg, element_type=Tri(),
                     approximation_type=SBP{StartUpDG.Kubatko{StartUpDG.LegendreFaceNodes}}(),
                     surface_integral=SurfaceIntegralWeakForm(surface_flux),
                     volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

    coordinates_min = (-1.0, -1.0)
    coordinates_max = ( 1.0,  1.0)
    cells_per_dimension = 2^refinement_level .* (1, 1)
    mesh = DGMultiMesh(solver;
                       coordinates_min, coordinates_max, cells_per_dimension,
                       periodicity=(true, true))

    return solver, mesh
  end

  function setup_solver_mesh(::Type{Polynomial},
                             polydeg, refinement_level,
                             equations, surface_flux, volume_flux)

    solver = DGMulti(polydeg=polydeg, element_type=Tri(),
                     approximation_type=Polynomial(),
                     surface_integral=SurfaceIntegralWeakForm(surface_flux),
                     volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

    coordinates_min = (-1.0, -1.0)
    coordinates_max = ( 1.0,  1.0)
    cells_per_dimension = 2^refinement_level .* (1, 1)
    mesh = DGMultiMesh(solver;
                       coordinates_min, coordinates_max, cells_per_dimension,
                       periodicity=(true, true))

    return solver, mesh
  end
end # @everywhere

# Prepare parameters
polydeg_and_refinement_level_quad = [(3, 5), (7, 4)]
polydeg_and_refinement_level_tri = [(3, 5), (6, 4)] # SBP() is not defined for polydeg=7
atwood_numbers = range(0.1, 0.9, length=20)

# Run simulations and save data
dgsem_parameters = [(atwood_number, DGSEM, polydeg, refinement_level) for (polydeg, refinement_level) in polydeg_and_refinement_level_quad for atwood_number in atwood_numbers]
sc_parameters = [(atwood_number, VolumeIntegralShockCapturingHG, polydeg, refinement_level) for (polydeg, refinement_level) in polydeg_and_refinement_level_quad for atwood_number in atwood_numbers]
gauss_parameters = [(atwood_number, GaussSBP, polydeg, refinement_level) for (polydeg, refinement_level) in polydeg_and_refinement_level_quad for atwood_number in atwood_numbers]
polynomial_parameters = [(atwood_number, Polynomial, polydeg, refinement_level) for (polydeg, refinement_level) in polydeg_and_refinement_level_tri for atwood_number in atwood_numbers]
sbp_parameters = [(atwood_number, SBP, polydeg, refinement_level) for (polydeg, refinement_level) in polydeg_and_refinement_level_tri for atwood_number in atwood_numbers]

parameters = vcat(dgsem_parameters, sc_parameters, gauss_parameters,
                  polynomial_parameters, sbp_parameters)
blow_up_times = pmap(args -> blow_up_time(args...), parameters)
data = DataFrame(variant = map(x -> x[2], parameters),
                 polydeg = map(x -> x[3], parameters),
                 refinement_level = map(x -> x[4], parameters),
                 atwood_number = map(x -> x[1], parameters),
                 blow_up_time = blow_up_times)
CSV.write(joinpath(@__DIR__, "kelvin_helmholtz_instability.txt"), data)


text:
mutable struct ParallelILU0Preconditioner{Tv, Ti} <: AbstractPreconditioner{Tv,Ti}
    A::ExtendableSparseMatrix{Tv,Ti}
    xdiag::Array{Tv,1}
    idiag::Array{Ti,1}
    phash::UInt64

    coloring::Array{Array{Ti,1},1}
    coloring_index::Array{Array{Ti,1},1}
    coloring_index_reverse::Array{Array{Ti,1},1}

    function ParallelILU0Preconditioner{Tv,Ti}() where {Tv,Ti}
        p=new()
        p.phash=0
        p
    end
end

"""
```
ParallelILU0Preconditioner(;valuetype=Float64,indextype=Int64)
ParallelILU0Preconditioner(matrix)
```

Parallel ILU preconditioner with zero fill-in.
"""
ParallelILU0Preconditioner(;valuetype::Type=Float64, indextype::Type=Int64)=ParallelILU0Preconditioner{valuetype,indextype}()

function update!(precon::ParallelILU0Preconditioner{Tv,Ti}) where {Tv,Ti}
    flush!(precon.A)

    # Get coloring and reorder matrix
    precon.coloring=graphcol(precon.A.cscmatrix)
    precon.coloring_index, precon.coloring_index_reverse=coloringindex(precon.coloring)
    precon.A=ExtendableSparseMatrix(reordermatrix(precon.A.cscmatrix, precon.coloring))

    cscmatrix=precon.A.cscmatrix
    colptr=cscmatrix.colptr
    rowval=cscmatrix.rowval
    nzval=cscmatrix.nzval
    n=cscmatrix.n

    if precon.phash==0
        n=size(precon.A,1)
        precon.xdiag=Array{Tv,1}(undef,n)
        precon.idiag=Array{Ti,1}(undef,n)
    end

    xdiag=precon.xdiag
    idiag=precon.idiag

    # Find main diagonal index and
    # copy main diagonal values
    if precon.phash != precon.A.phash
        @inbounds for j=1:n
            @inbounds for k=colptr[j]:colptr[j+1]-1
                i=rowval[k]
                if i==j
                    idiag[j]=k
                    break
                end
            end
        end
        precon.phash=precon.A.phash
    end

    @inbounds for j=1:n
        xdiag[j]=one(Tv)/nzval[idiag[j]]
        @inbounds for k=idiag[j]+1:colptr[j+1]-1
            i=rowval[k]
            for l=colptr[i]:colptr[i+1]-1
                if rowval[l]==j
                    xdiag[i]-=nzval[l]*xdiag[j]*nzval[k]
                    break
                end
            end
        end
    end
    precon
end

function  LinearAlgebra.ldiv!(u::AbstractArray{T,1}, precon::ParallelILU0Preconditioner, v::AbstractArray{T,1}) where T
    cscmatrix=precon.A.cscmatrix
    colptr=cscmatrix.colptr
    rowval=cscmatrix.rowval
    n=cscmatrix.n
    nzval=cscmatrix.nzval
    xdiag=precon.xdiag
    idiag=precon.idiag

    coloring = precon.coloring
    coloring_index = precon.coloring_index
    coloring_index_reverse = precon.coloring_index_reverse

    @inbounds for indset in coloring_index
	    @inbounds Threads.@threads for j in indset
	        x=zero(T)
	        @inbounds for k=colptr[j]:idiag[j]-1
	            x+=nzval[k]*u[rowval[k]]
	        end
	        u[j]=xdiag[j]*(v[j]-x)
	    end
	end

    @inbounds for indset in coloring_index_reverse
    	@inbounds Threads.@threads for j in indset
	        x=zero(T)
	        @inbounds for k=idiag[j]+1:colptr[j+1]-1
	            x+=u[rowval[k]]*nzval[k]
	        end
	        u[j]-=x*xdiag[j]
	    end
	end
end

function LinearAlgebra.ldiv!(precon::ParallelILU0Preconditioner, v::AbstractArray{T,1} where T)
    ldiv!(v, precon, v)
end

# Returns an independent set of the graph of a matrix
# Reference: https://research.nvidia.com/sites/default/files/pubs/2015-05_Parallel-Graph-Coloring/nvr-2015-001.pdf
function indset(A::SparseMatrixCSC{Tv,Ti}, W::StridedVector{Ti}) where {Tv,Ti}
	# Random numbers for all vertices
	lenW = length(W)
	# r = sample(1:lenW, lenW, replace = false)
	r = rand(lenW)
	@inbounds for i = 1:lenW
		if W[i] == 0
			r[i] = 0
		end
	end
	# Empty independent set
	S = zeros(Int, lenW)
	# Get independent set by comparing random number of vertex with the random
	# numbers of all neighbor vertices
	@inbounds Threads.@threads for i in 1:lenW
		if W[i] != 0
			j = A.rowval[A.colptr[i]:A.colptr[i+1]-1]
			if all(x->x==1, r[i] .>= r[j])
				S[i] = W[i]
			end
		end
	end
	# Remove zero entries and return independent set
	return filter!(x->x≠0, S)
end

# Returns coloring of the graph of a matrix
# Reference: https://research.nvidia.com/sites/default/files/pubs/2015-05_Parallel-Graph-Coloring/nvr-2015-001.pdf
function graphcol(A::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
	# Empty list for coloring
	C = []
	# Array of vertices
	W = [1:size(A)[1];]
	# Get all independent sets of the graph of the matrix
	while any(W .!= 0)
		# Get independent set
		S = indset(A + transpose(A), W)
		push!(C, S)
		# Remove entries in S from W
		@inbounds for s in S
			W[s] = 0
		end
	end
	# Return coloring
	return C
end

# Reorders a sparse matrix with provided coloring
function reordermatrix(A::SparseMatrixCSC{Tv,Ti}, coloring::Array{Array{Int64,1},1}) where {Tv,Ti}
	c = collect(Iterators.flatten(coloring))
	return A[c,:][:,c]
end

# Reorders a linear system with provided coloring
function reorderlinsys(A::SparseMatrixCSC{Tv,Ti}, b::StridedVector{Tv}, coloring::Array{Array{Int64,1},1}) where {Tv,Ti}
	c = collect(Iterators.flatten(coloring))
	return A[c,:][:,c], b[c]
end

# Returns an array with the same structure of the input coloring and ordered
# entries 1:length(coloring) and an array with the structure of
# reverse(coloring) and ordered entries length(coloring):-1:1
function coloringindex(coloring::Array{Array{Int64,1},1})
    # First array
	c = deepcopy(coloring)
	cnt = 1
	@inbounds for i in 1:length(c)
		@inbounds for j in 1:length(c[i])
			c[i][j] = cnt
			cnt += 1
		end
	end
    # Second array
	cc = deepcopy(reverse(coloring))
	@inbounds for i in 1:length(cc)
		@inbounds for j in 1:length(cc[i])
			cnt -= 1
			cc[i][j] = cnt
		end
	end
    # Return both
	return c, cc
end


text:
function parse_script(document_body)
    header_text, document_body, offset = separate_header_text_jl(document_body)
    header_text = replace(header_text, r"\n(#')|(#%%)|(# %%)" => "\n")
    return parse_header(header_text), parse_script_body(document_body)
end

# headers
# -------

const HEADER_REGEX_JL = r"^(?:#'|#%%|# %%)\s*---$(?<header>((?!---).)+)^(?:#'|#%%|# %%)\s*---$"ms

function separate_header_text_jl(text)
    m = match(HEADER_REGEX_JL, text)
    isnothing(m) && return "", text, 0
    header_text = m[:header]
    offset = @static if VERSION ≥ v"1.4"
        count("\n", header_text)
    else
        count(c->c==='\n', header_text)
    end
    return header_text, replace(text, HEADER_REGEX_JL => ""; count = 1), offset
end

# body
# ----

function parse_script_body(document_body)
    lines = split(document_body, '\n')

    doc_line = r"(^#'.*)|(^#%%.*)|(^# %%.*)"
    doc_start = r"(^#')|(^#%%)|(^# %%)"
    opt_line = r"(^#\+.*$)|(^#%%\+.*$)|(^# %%\+.*$)"
    opt_start = r"(^#\+)|(^#%%\+)|(^# %%\+)"

    content = ""
    state = :code
    doc_no = 0
    code_no = 0
    start_line = 1

    options = OptionDict()
    option_string = ""

    chunks = WeaveChunk[]
    for (line_no, line) in enumerate(lines)
        if (m = match(doc_line, line)) !== nothing && (m = match(opt_line, line)) === nothing
            line = replace(line, doc_start => "", count = 1)
            startswith(line, ' ') && (line = replace(line, ' ' => "", count = 1))
            if state === :code && !isempty(strip(content))
                push!(chunks, CodeChunk(string('\n', strip(content)), code_no += 1, start_line, option_string, options))
                content = ""
                start_line = line_no
            end
            state = :doc
        elseif (m = match(opt_line, line)) !== nothing
            start_line = line_no
            if state === :code && !isempty(strip(content))
                push!(chunks, CodeChunk(string('\n', strip(content)), code_no += 1, start_line, option_string, options))
                content = ""
            end
            if state === :doc && !isempty(strip(content))
                iszero(doc_no) || (content = string('\n', content)) # Add whitespace to doc chunk. Needed for markdown output
                push!(chunks, DocChunk(content, doc_no += 1, start_line))
                content = ""
            end

            option_string = replace(line, opt_start => "", count = 1)
            options = parse_options(option_string)
            haskey(options, :label) && (options[:name] = options[:label])
            haskey(options, :name) || (options[:name] = nothing)

            state = :code
            continue
        elseif state === :doc # && strip(line) != "" && strip(content) != ""
            state = :code
            iszero(doc_no) || (content = string('\n', content)) # Add whitespace to doc chunk. Needed for markdown output
            push!(chunks, DocChunk(content, doc_no += 1, start_line))
            content = ""

            options = Dict{Symbol,Any}()
            start_line = line_no
        end
        content *= string(line, '\n')
    end

    # Handle the last chunk
    chunk = state === :code ?
        CodeChunk(string('\n', strip(content)), code_no, start_line, option_string, options) :
        DocChunk(content, doc_no, start_line)
    push!(chunks, chunk)

    return chunks
end


text:
using DASSL, Test

@testset "Testing maxorder" begin

    F(t,y,dy)=(dy+y.^2)
    Fy(t,y,dy)=diagm(2y)
    Fdy(t,y,dy)=eye(length(y))

    sol(t)=1.0./(1 .+ t)
    tspan=[0.0, 1.0]

    atol = 1.0e-5
    rtol = 1.0e-3

    # test the maxorder option
    for order=1:6
        # scalar version
        (tn,yn,dyn)=DASSL.dasslSolve(F, sol(0.0), tspan, maxorder = order)
        aerror = maximum(abs.(yn-sol(tn)))
        rerror = maximum(abs.(yn-sol(tn))./abs.(sol(tn)))
        nsteps = length(tn)

        @test aerror < (2*nsteps*atol)
        @test rerror < (2*nsteps*rtol)

        # vector version
        (tnV,ynV,dynV)=DASSL.dasslSolve(F,[sol(0.0)], tspan, maxorder = order)

        @test  vcat(ynV...) == yn
        @test  vcat(dynV...) == dyn

        # analytical jacobian version (vector)
        (tna,yna,dyna)=dasslSolve(F, [sol(0.0)], tspan, maxorder = order, Fy = Fy, Fdy = Fdy)
        aerror = maximum(abs.(map(first,yn)-sol(tn)))
        rerror = maximum(abs.(map(first,yn)-sol(tn))./abs.(sol(tn)))
        nsteps = length(tn)

        @test aerror < (2*nsteps*atol)
        @test rerror < (2*nsteps*rtol)
    end
end

@testset "Testing common interface" begin
  include("common.jl")
end


text:
# Bayesian student-t multivariate regression with independent Gaussian priors

type SteinMultivariateStudentTRegressionGaussianPrior <: SteinMultivariateStudentTRegressionPosterior
    X::Array{Float64, 2}
    y::Array{Float64, 1}
    nu::Float64  # the number of degrees of freedom for the t-distribution
    sigma::Array{Float64, 2} # the covariance of the multivariate t-distribution
    priormu::Array{Float64, 1}
    priorvar::Array{Float64, 1}  # we assume components are independent
    c1::Float64
    c2::Float64
    c3::Float64
end

SteinMultivariateStudentTRegressionGaussianPrior(X::Array{Float64, 2},
                                                 y::Array{Float64, 1},
                                                 nu::Float64) =
SteinMultivariateStudentTRegressionGaussianPrior(X, y, nu, eye(size(X, 1)))

SteinMultivariateStudentTRegressionGaussianPrior(X::Array{Float64, 2},
                                                 y::Array{Float64, 1},
                                                 nu::Float64,
                                                 sigma::Array{Float64, 2}) =
SteinMultivariateStudentTRegressionGaussianPrior(X, y, nu, sigma,
                                                 repmat([0.0], size(X, 2)),
                                                 repmat([1.0], size(X, 2)))

# Constructor with canonical setting of (c1,c2,c3) = (1,1,1)
SteinMultivariateStudentTRegressionGaussianPrior(X::Array{Float64, 2},
                                                 y::Array{Float64, 1},
                                                 nu::Float64,
                                                 sigma::Array{Float64, 2},
                                                 priormu::Array{Float64, 1},
                                                 priorvar::Array{Float64, 1}) =
SteinMultivariateStudentTRegressionGaussianPrior(X, y, nu, sigma, priormu, priorvar,
                                                 1.0, 1.0, 1.0)

function logprior(d::SteinMultivariateStudentTRegressionGaussianPrior,
                  beta::Array{Float64,1})
    logpriorvalue =
        -0.5 * sum((beta .- d.priormu) .^ 2 ./ d.priorvar) -
        -0.5 * sum(log(2 * pi .* d.priorvar))
    logpriorvalue
end

function gradlogprior(d::SteinMultivariateStudentTRegressionGaussianPrior,
                      beta::Array{Float64,1})
    -(beta .- d.priormu) ./ d.priorvar
end

function numdimensions(d::SteinMultivariateStudentTRegressionGaussianPrior)
    length(d.priormu)
end


text:
using StatsBase
using Test, Random, Statistics, LinearAlgebra
using Missings

const ≅ = isequal

Random.seed!(1)

# to avoid using specialized method
arbitrary_fun(x, y) = cor(x, y)

@testset "pairwise and pairwise! with $f" for f in (arbitrary_fun, cor, cov)
    @testset "basic interface" begin
        x = [rand(10) for _ in 1:4]
        y = [rand(Float32, 10) for _ in 1:5]
        # to test case where inference of returned eltype fails
        z = [Vector{Any}(rand(Float32, 10)) for _ in 1:5]

        res = @inferred pairwise(f, x, y)
        @test res isa Matrix{Float64}
        res2 = zeros(Float64, size(res))
        @test pairwise!(f, res2, x, y) === res2
        @test res == res2 == [f(xi, yi) for xi in x, yi in y]

        res = pairwise(f, y, z)
        @test res isa Matrix{Float32}
        res2 = zeros(Float32, size(res))
        @test pairwise!(f, res2, y, z) === res2
        @test res == res2 == [f(yi, zi) for yi in y, zi in z]

        res = pairwise(f, Any[[1.0, 2.0, 3.0], [1.0f0, 3.0f0, 10.5f0]])
        @test res isa Matrix{Float64}
        res2 = zeros(AbstractFloat, size(res))
        @test pairwise!(f, res2, Any[[1.0, 2.0, 3.0], [1.0f0, 3.0f0, 10.5f0]]) === res2
        @test res == res2 ==
            [f(xi, yi) for xi in ([1.0, 2.0, 3.0], [1.0f0, 3.0f0, 10.5f0]),
                           yi in ([1.0, 2.0, 3.0], [1.0f0, 3.0f0, 10.5f0])]
        @test res isa Matrix{Float64}

        @inferred pairwise(f, x, y)

        @test_throws ArgumentError pairwise(f, [Int[]], [Int[]])
        @test_throws ArgumentError pairwise!(f, zeros(1, 1), [Int[]], [Int[]])

        res = pairwise(f, [], [])
        @test size(res) == (0, 0)
        @test res isa Matrix{Any}
        res2 = zeros(0, 0)
        @test pairwise!(f, res2, [], []) === res2

        res = pairwise(f, Vector{Int}[], Vector{Int}[])
        @test size(res) == (0, 0)
        @test res isa Matrix{Float64}
        res2 = zeros(0, 0)
        @test pairwise!(f, res2, Vector{Int}[], Vector{Int}[]) === res2

        res = pairwise(f, [[1, 2]], Vector{Int}[])
        @test size(res) == (1, 0)
        @test res isa Matrix{Float64}
        res2 = zeros(1, 0)
        @test pairwise!(f, res2, [[1, 2]], Vector{Int}[]) === res2

        res = pairwise(f, Vector{Int}[], [[1, 2], [2, 3]])
        @test size(res) == (0, 2)
        @test res isa Matrix{Float64}
        res2 = zeros(0, 2)
        @test pairwise!(f, res2, [], [[1, 2], [2, 3]]) === res2

        @test_throws DimensionMismatch pairwise!(f, zeros(1, 2), x, y)
        @test_throws DimensionMismatch pairwise!(f, zeros(1, 2), [], [])
        @test_throws DimensionMismatch pairwise!(f, zeros(0, 0),
                                                 [], [[1, 2], [2, 3]])
    end

    @testset "missing values handling interface" begin
        xm = [ifelse.(rand(100) .> 0.9, missing, rand(100)) for _ in 1:4]
        ym = [ifelse.(rand(100) .> 0.9, missing, rand(Float32, 100)) for _ in 1:4]
        zm = [ifelse.(rand(100) .> 0.9, missing, rand(Float32, 100)) for _ in 1:4]

        res = pairwise(f, xm, ym)
        @test res isa Matrix{Missing}
        res2 = zeros(Union{Float64, Missing}, size(res))
        @test pairwise!(f, res2, xm, ym) === res2
        @test res ≅ res2 ≅ [missing for xi in xm, yi in ym]

        res = pairwise(f, xm, ym, skipmissing=:pairwise)
        @test res isa Matrix{Float64}
        res2 = zeros(Union{Float64, Missing}, size(res))
        @test pairwise!(f, res2, xm, ym, skipmissing=:pairwise) === res2
        @test res ≅ res2
        @test isapprox(res, [f(collect.(skipmissings(xi, yi))...) for xi in xm, yi in ym],
                       rtol=1e-6)

        res = pairwise(f, ym, zm, skipmissing=:pairwise)
        @test res isa Matrix{Float32}
        res2 = zeros(Union{Float32, Missing}, size(res))
        @test pairwise!(f, res2, ym, zm, skipmissing=:pairwise) === res2
        @test res ≅ res2
        @test isapprox(res, [f(collect.(skipmissings(yi, zi))...) for yi in ym, zi in zm],
                       rtol=1e-6)

        nminds = mapreduce(x -> .!ismissing.(x),
                        (x, y) -> x .& y,
                        [xm; ym])
        res = pairwise(f, xm, ym, skipmissing=:listwise)
        @test res isa Matrix{Float64}
        res2 = zeros(Union{Float64, Missing}, size(res))
        @test pairwise!(f, res2, xm, ym, skipmissing=:listwise) === res2
        @test res ≅ res2
        @test isapprox(res, [f(view(xi, nminds), view(yi, nminds)) for xi in xm, yi in ym],
                       rtol=1e-6)

        if VERSION >= v"1.6.0-DEV"
            # inference of cor fails so use an inferrable function
            # to check that pairwise itself is inferrable
            for skipmissing in (:none, :pairwise, :listwise)
                g(x, y=x) = pairwise((x, y) -> x[1] * y[1], x, y, skipmissing=skipmissing)
                @test Core.Compiler.return_type(g, Tuple{Vector{Vector{Union{Float64, Missing}}}}) ==
                    Core.Compiler.return_type(g, Tuple{Vector{Vector{Union{Float64, Missing}}},
                                                    Vector{Vector{Union{Float64, Missing}}}}) ==
                    Matrix{<: Union{Float64, Missing}}
                if skipmissing in (:pairwise, :listwise)
                    @test_broken Core.Compiler.return_type(g, Tuple{Vector{Vector{Union{Float64, Missing}}}}) ==
                        Core.Compiler.return_type(g, Tuple{Vector{Vector{Union{Float64, Missing}}},
                                                        Vector{Vector{Union{Float64, Missing}}}}) ==
                        Matrix{Float64}
                end
            end
        end

        @test_throws ArgumentError pairwise(f, xm, ym, skipmissing=:something)
        @test_throws ArgumentError pairwise!(f, zeros(Union{Float64, Missing},
                                                      length(xm), length(ym)), xm, ym,
                                             skipmissing=:something)

        # variable with only missings
        xm = [fill(missing, 10), rand(10)]
        ym = [rand(10), rand(10)]

        res = pairwise(f, xm, ym)
        @test res isa Matrix{Union{Float64, Missing}}
        res2 = zeros(Union{Float64, Missing}, size(res))
        @test pairwise!(f, res2, xm, ym) === res2
        @test res ≅ res2 ≅ [f(xi, yi) for xi in xm, yi in ym]

        if VERSION >= v"1.5" # Fails with UndefVarError on Julia 1.0
            @test_throws ArgumentError pairwise(f, xm, ym, skipmissing=:pairwise)
            @test_throws ArgumentError pairwise(f, xm, ym, skipmissing=:listwise)

            res = zeros(Union{Float64, Missing}, length(xm), length(ym))
            @test_throws ArgumentError pairwise!(f, res, xm, ym, skipmissing=:pairwise)
            @test_throws ArgumentError pairwise!(f, res, xm, ym, skipmissing=:listwise)
        end

        for sm in (:pairwise, :listwise)
            @test_throws ArgumentError pairwise(f, [[1, 2]], [1], skipmissing=sm)
            @test_throws ArgumentError pairwise(f, [1], [[1, 2]], skipmissing=sm)
            @test_throws ArgumentError pairwise(f, [1], [1], skipmissing=sm)
        end
    end

    @testset "iterators" begin
        x = (v for v in [rand(10) for _ in 1:4])
        y = (v for v in [rand(10) for _ in 1:4])

        res = @inferred pairwise(f, x, y)
        res2 = zeros(size(res))
        @test pairwise!(f, res2, x, y) === res2
        @test res == res2 == pairwise(f, collect(x), collect(y))

        res = @inferred(pairwise(f, x))
        res2 = zeros(size(res))
        @test pairwise!(f, res2, x) === res2
        @test res == res2 == pairwise(f, collect(x))
    end

    @testset "non-vector entries" begin
        x = (Iterators.drop(v, 1) for v in [rand(10) for _ in 1:4])
        y = (Iterators.drop(v, 1) for v in [rand(10) for _ in 1:4])

        @test pairwise((x, y) -> f(collect(x), collect(y)), x, y) ==
            [f(collect(xi), collect(yi)) for xi in x, yi in y]
        @test pairwise((x, y) -> f(collect(x), collect(y)), x) ==
            [f(collect(xi1), collect(xi2)) for xi1 in x, xi2 in x]
        @test_throws ArgumentError pairwise((x, y) -> f(collect(x), collect(y)), x, y,
                                            skipmissing=:pairwise)
        @test_throws ArgumentError pairwise((x, y) -> f(collect(x), collect(y)), x, y,
                                            skipmissing=:listwise)
    end

    @testset "two-argument method" begin
        x = [rand(10) for _ in 1:4]
        res = pairwise(f, x)
        res2 = zeros(size(res))
        @test pairwise!(f, res2, x) === res2
        @test res == res2 == pairwise(f, x, x)
    end

    @testset "symmetric" begin
        x = [rand(10) for _ in 1:4]
        y = [rand(10) for _ in 1:4]

        @test pairwise(f, x, x, symmetric=true) ==
            pairwise(f, x, symmetric=true) ==
            Symmetric(pairwise(f, x, x), :U)

        res = zeros(4, 4)
        res2 = zeros(4, 4)
        @test pairwise!(f, res, x, x, symmetric=true) === res
        @test pairwise!(f, res2, x, symmetric=true) === res2
        @test res == res2 == Symmetric(pairwise(f, x, x), :U)

        @test_throws ArgumentError pairwise(f, x, y, symmetric=true)
        @test_throws ArgumentError pairwise!(f, res, x, y, symmetric=true)
    end

    @testset "cor corner cases" begin
        # Integer inputs must give a Float64 output
        res = pairwise(cor, [[1, 2, 3], [1, 5, 2]])
        @test res isa Matrix{Float64}
        @test res == [cor(xi, yi) for xi in ([1, 2, 3], [1, 5, 2]),
                                      yi in ([1, 2, 3], [1, 5, 2])]

        # NaNs are ignored for the diagonal
        res = pairwise(cor, [[1, 2, NaN], [1, 5, 2]])
        @test res isa Matrix{Float64}
        @test res ≅ [1.0 NaN
                     NaN 1.0]

        # missings are ignored for the diagonal
        res = pairwise(cor, [[1, 2, 7], [1, 5, missing]])
        @test res isa Matrix{Union{Float64, Missing}}
        @test res ≅ [1.0 missing
                     missing 1.0]
        res = pairwise(cor, Vector{Union{Int, Missing}}[[missing, missing, missing],
                                                        [missing, missing, missing]])
        @test res isa Matrix{Union{Float64, Missing}}
        @test res ≅ [1.0 missing
                     missing 1.0]
        if VERSION >= v"1.5"
            # except when eltype is Missing
            res = pairwise(cor, [[missing, missing, missing],
                                 [missing, missing, missing]])
            @test res isa Matrix{Missing}
            @test res ≅ [missing missing
                         missing missing]
        end

        for sm in (:pairwise, :listwise)
            res = pairwise(cor, [[1, 2, NaN, 4], [1, 5, 5, missing]], skipmissing=sm)
            @test res isa Matrix{Float64}
            @test res ≅ [1.0 NaN
                         NaN 1.0]
            if VERSION >= v"1.5"
                @test_throws ArgumentError pairwise(cor, [[missing, missing, missing],
                                                          [missing, missing, missing]],
                                                    skipmissing=sm)
            end
        end
    end
end


text:
#!/usr/bin/env julia
module Ownway
using Printf: @printf

function create_movedict(lydiastr)
    current = [1, 1]
    movedict = Dict{Tuple{Int,Int},Char}()

    for c in lydiastr
        movedict[current...] = c
        current[c === 'E' ? 1 : 2] += 1
    end

    movedict
end

function solve_task(bsize, lydiastr)
    lydiamoves = create_movedict(lydiastr)
    is_solved(s) = length(s) == 2bsize - 2

    function solve_maze(x, y, moves)
        is_solved(moves) && x == y && return moves
        (x > bsize || y > bsize) && return moves[1:end - 1]
        lm = get(lydiamoves, (x, y), nothing)

        move_s() = solve_maze(x, y + 1, vcat(moves, 'S'))
        move_e() = solve_maze(x + 1, y, vcat(moves, 'E'))

        lm === 'E' && return move_s()
        lm === 'S' && return move_e()

        ret = move_e()
        is_solved(ret) && return ret
        return move_s()
    end

    return join(solve_maze(1, 1, Char[]), "")
end

function main()
    numtasks = parse(Int, readline())

    for i in 1:numtasks
        bsize = parse(Int, readline())
        lydiastr = readline()
        sol = solve_task(bsize, lydiastr)
        @printf "Case #%d: %s\n" i sol
    end
end
end

isinteractive() || Ownway.main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# =========================================================================================
"""
    struct TessParam{T}
        name::String
        default::T
        desc::String
        debug::Bool
    end

Holds details about the default value of a parameter.

__Values:__

| Name    | Description
| :------ | :----------
| name    | The name/ID of the parameter.
| default | The default value of the parameter.
| desc    | The description of the variable.
| debug   | True if the value is a debug parameter.

__Details:__

This structure currently comes in 3 flavors, T may be Float64, Int32, or a String based on
the default value.

See also: [`tess_params_parsed`](@ref), [`tess_get_param`](@ref), [`tess_set_param`](@ref)
"""
struct TessParam{T}
    name::String
    default::T
    desc::String
    debug::Bool
end

# =========================================================================================
"""
    TessParam(
        name::AbstractString,
        value::AbstractString,
        help::AbstractString
    )::TessParam{T} where T

Construct a new instance of the [`TessParam`](@ref) structure from the provided values.

__Arguments:__

| T | Name  | Default | Description
|:--| :---- | :------ | :----------
| R | name  |         | The name/ID of the parameter.
| R | value |         | The default value of the parameter as a string.
| R | help  |         | Help text about the parameter.

__Details:__

This method looks at the contents of `value` and determines the correct [`TessParam`](@ref)
type to create.

See also: [`tess_params_parsed`](@ref)
"""
function TessParam(
            name::AbstractString,
            value::AbstractString,
            help::AbstractString
        )::TessParam
    local debug = findfirst("debug", name) !== nothing ||
                  findfirst("display", name) !== nothing

    if match(r"^-?[0-9]+\.[0-9]+$", value) != nothing
        return TessParam{Float64}(name, parse(Float64, value), help, debug)
    elseif match(r"^-?[0-9]+$", value) != nothing
        return TessParam{Int32}(name, parse(Int32, value), help, debug)
    else
        return TessParam{String}(name, value, help, debug)
    end
end

# =========================================================================================
"""
    tess_params_parsed(
            inst::TessInst
        )::Vector{TessParam}

Retrieved all the Tesseract parameters with their valuse as an array of [`TessParam`](@ref)
objects.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | inst |         | The Tesseract instance to get the parameters from.

__Details:__

Parses the result of `[`tess_params`](@ref) into something more easily digested by a
computer.  Each line of text is split into 3 values:

* The name of the parameter.
* The default value of the parameter (may be an empty string).
* Text describing the parameter.

Each value is separated by a tab and the description is terminated by a new line.

See also: [`TessParam`](@ref), [`tess_params`](@ref), [`tess_get_param`](@ref),
          [`tess_set_param`](@ref)
"""
function tess_params_parsed(
            inst::TessInst
        )::Vector{TessParam}
    local params = Vector{TessParam}()
    local data   = tess_params(inst)

    if data != nothing
        for m in eachmatch(r"([a-z0-9_]+)\t([^\t]*)\t([^\n]+)"sm, data)
            local name, value, help = m.captures
            push!(params, TessParam(name, value, help))
        end
    end

    return params
end


text:
using BinaryVectorProbability
using DataFrames
using Base.Test

X = rand(10,10000); p = collect(0.2:0.05:0.65); X = X .< p;
P = first_order_model(X); P2 = first_order_model(X);
P3 = BernoulliCodeDistribution(P.p)

@test isequal(P, P2)
@test !isequal(P, P3)

@test hash(P) == hash(P2)
@test hash(P) != hash(P3)

Base.Random.srand(17)

N_neurons = 10
Jtrue = rand(N_neurons, N_neurons); Jtrue = (Jtrue + Jtrue')/2
P2 = IsingDistribution(Jtrue; comment="true distribution")
X = random(P2, 10000)

P2hat = second_order_model(X, verbose=2)

max_stars_repo_name: sekunder/BinaryVectorProbability

max_issues_repo_name: sekunder/BinaryVectorProbability

max_forks_repo_name: sekunder/BinaryVectorProbability


text:
@testset "OUUnknown" begin
    θ₀ = [2.0,0.4]
    model_known_acv = OU(θ₀)
    model_unknown_acv = OUUnknown{5}(θ₀)
    ω = 0.6
    @testset "Function sdf" begin
        @test sdf(model_known_acv, ω)      == sdf(model_unknown_acv, ω)
    end
    @testset "Gradient sdf" begin
        @test grad_sdf(model_known_acv, ω) == grad_sdf(model_unknown_acv, ω)
    end
    @testset "Hessian sdf" begin
        @test hess_sdf(model_known_acv, ω) == hess_sdf(model_unknown_acv, ω)
    end
    @testset "Error handling" begin
        @test_throws ArgumentError OUUnknown{2}(-1.0,1.0)
        @test_throws ArgumentError OUUnknown{2}(1.0,-1.0)
        @test_throws ArgumentError OUUnknown{2}(ones(1))
        @test_throws ArgumentError OUUnknown{2}(ones(3))
    end
end


text:
#=

This file contains RT doubling-related functions

=#

"""
    $(FUNCTIONNAME)(pol_type, SFI, expk, ndoubl::Int, added_layer::AddedLayer, I_static::AbstractArray{FT},
                    architecture) where {FT}

Compute homogenous layer matrices from its elemental layer using Doubling
"""
function doubling_helper!(pol_type,
                          SFI,
                          expk,
                          ndoubl::Int,
                          added_layer::AddedLayer,
                          I_static::AbstractArray{FT},
                          architecture) where {FT}

    # Unpack the added layer
    @unpack r⁺⁻, r⁻⁺, t⁻⁻, t⁺⁺, J₀⁺, J₀⁻ = added_layer

    # Device architecture
    dev = devi(architecture)

    # Note: short-circuit evaluation => return nothing evaluated iff ndoubl == 0
    ndoubl == 0 && return nothing

    # Geometric progression of reflections (1-RR)⁻¹
    gp_refl      = similar(t⁺⁺)
    tt⁺⁺_gp_refl = similar(t⁺⁺)

    if SFI
        # Dummy for source
        J₁⁺ = similar(J₀⁺)
        # Dummy for J
        J₁⁻ = similar(J₀⁻)
    end

    # Loop over number of doublings
    for n = 1:ndoubl

        # T⁺⁺(λ)[I - R⁺⁻(λ)R⁻⁺(λ)]⁻¹, for doubling R⁺⁻,R⁻⁺ and T⁺⁺,T⁻⁻ is identical
        batch_inv!(gp_refl, I_static .- r⁻⁺ ⊠ r⁻⁺)
        tt⁺⁺_gp_refl[:] = t⁺⁺ ⊠ gp_refl

        if SFI

            # J⁺₂₁(λ) = J⁺₁₀(λ).exp(-τ(λ)/μ₀)
            J₁⁺[:,1,:] = J₀⁺[:,1,:] .* expk'

            # J⁻₁₂(λ)  = J⁻₀₁(λ).exp(-τ(λ)/μ₀)
            J₁⁻[:,1,:] = J₀⁻[:,1,:] .* expk'

            # J⁻₀₂(λ) = J⁻₀₁(λ) + T⁻⁻₀₁(λ)[I - R⁻⁺₂₁(λ)R⁺⁻₀₁(λ)]⁻¹[J⁻₁₂(λ) + R⁻⁺₂₁(λ)J⁺₁₀(λ)] (see Eqs.8 in Raman paper draft)
            J₀⁻[:] = J₀⁻ + (tt⁺⁺_gp_refl ⊠ (J₁⁻ + r⁻⁺ ⊠ J₀⁺))

            # J⁺₂₀(λ) = J⁺₂₁(λ) + T⁺⁺₂₁(λ)[I - R⁺⁻₀₁(λ)R⁻⁺₂₁(λ)]⁻¹[J⁺₁₀(λ) + R⁺⁻₀₁(λ)J⁻₁₂(λ)] (see Eqs.8 in Raman paper draft)
            J₀⁺[:] = J₁⁺ + (tt⁺⁺_gp_refl ⊠ (J₀⁺ + r⁻⁺ ⊠ J₁⁻))
            expk[:] = expk.^2
        end

        # R⁻⁺₂₀(λ) = R⁻⁺₁₀(λ) + T⁻⁻₀₁(λ)[I - R⁻⁺₂₁(λ)R⁺⁻₀₁(λ)]⁻¹R⁻⁺₂₁(λ)T⁺⁺₁₀(λ) (see Eqs.8 in Raman paper draft)
        r⁻⁺[:]  = r⁻⁺ + (tt⁺⁺_gp_refl ⊠ r⁻⁺ ⊠ t⁺⁺)

        # T⁺⁺₂₀(λ) = T⁺⁺₂₁(λ)[I - R⁺⁻₀₁(λ)R⁻⁺₂₁(λ)]⁻¹T⁺⁺₁₀(λ) (see Eqs.8 in Raman paper draft)
        t⁺⁺[:]  = tt⁺⁺_gp_refl ⊠ t⁺⁺
    end

    # After doubling, revert D(DR)->R, where D = Diagonal{1,1,-1,-1}
    # For SFI, after doubling, revert D(DJ₀⁻)->J₀⁻

    synchronize_if_gpu()

    apply_D_matrix!(pol_type.n, added_layer.r⁻⁺, added_layer.t⁺⁺, added_layer.r⁺⁻, added_layer.t⁻⁻)

    SFI && apply_D_matrix_SFI!(pol_type.n, added_layer.J₀⁻)

    return nothing

end

function doubling!(pol_type, SFI, expk,
                    ndoubl::Int,
                    added_layer::AddedLayer,#{FT},
                    I_static::AbstractArray{FT},
                    architecture) where {FT}

    doubling_helper!(pol_type, SFI, expk, ndoubl, added_layer, I_static, architecture)
    synchronize_if_gpu()
end

@kernel function apply_D!(n_stokes::Int,  r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻)
    iμ, jμ, n = @index(Global, NTuple)
    i = mod(iμ, n_stokes)
    j = mod(jμ, n_stokes)

    if (i > 2)
        r⁻⁺[iμ,jμ,n] = - r⁻⁺[iμ, jμ,n]
    end

    if ((i <= 2) & (j <= 2)) | ((i > 2) & (j > 2))
        r⁺⁻[iμ,jμ,n] = r⁻⁺[iμ,jμ,n]
        t⁻⁻[iμ,jμ,n] = t⁺⁺[iμ,jμ,n]
    else
        r⁺⁻[iμ,jμ,n] = - r⁻⁺[iμ,jμ,n]
        t⁻⁻[iμ,jμ,n] = - t⁺⁺[iμ,jμ,n]
    end

end

@kernel function apply_D_SFI!(n_stokes::Int, J₀⁻)
    iμ, _, n = @index(Global, NTuple)
    i = mod(iμ, n_stokes)

    if (i > 2)
        J₀⁻[iμ, 1, n] = - J₀⁻[iμ, 1, n]
    end
end

function apply_D_matrix!(n_stokes::Int, r⁻⁺::CuArray{FT,3}, t⁺⁺::CuArray{FT,3}, r⁺⁻::CuArray{FT,3}, t⁻⁻::CuArray{FT,3}) where {FT}
    if n_stokes == 1
        r⁺⁻[:] = r⁻⁺
        t⁻⁻[:] = t⁺⁺

        return nothing
    else
        device = devi(architecture(r⁻⁺))
        applyD_kernel! = apply_D!(device)
        event = applyD_kernel!(n_stokes, r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻, ndrange=size(r⁻⁺));
        wait(device, event);
        synchronize_if_gpu();
        return nothing
    end
end

function apply_D_matrix!(n_stokes::Int, r⁻⁺::Array{FT,3}, t⁺⁺::Array{FT,3}, r⁺⁻::Array{FT,3}, t⁻⁻::Array{FT,3}) where {FT}
    if n_stokes == 1
        r⁺⁻[:] = r⁻⁺
        t⁻⁻[:] = t⁺⁺

        return nothing
    else
        device = devi(Architectures.CPU())
        applyD_kernel! = apply_D!(device)
        event = applyD_kernel!(n_stokes, r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻, ndrange=size(r⁻⁺));
        wait(device, event);
        return nothing
    end
end

function apply_D_matrix_SFI!(n_stokes::Int, J₀⁻::CuArray{FT,3}) where {FT}

    n_stokes == 1 && return nothing

    device = devi(architecture(J₀⁻))
    applyD_kernel! = apply_D_SFI!(device)
    event = applyD_kernel!(n_stokes, J₀⁻, ndrange=size(J₀⁻));
    wait(device, event);
    synchronize();

    return nothing
end

function apply_D_matrix_SFI!(n_stokes::Int, J₀⁻::Array{FT,3}) where {FT}

    n_stokes == 1 && return nothing

    device = devi(architecture(J₀⁻))
    applyD_kernel! = apply_D_SFI!(device)
    event = applyD_kernel!(n_stokes, J₀⁻, ndrange=size(J₀⁻));
    wait(device, event);

    return nothing
end


text:
"""
    CenterPoint(dims) -> cp

Create a fixed point which can be used in `imresize` and `imrotate`
functions in order to keep the value in this point the same, i.e.,

```julia
img[cp] == imgr[cp]
```
"""
struct CenterPoint{N}
    p::CartesianIndex{N}
end

CenterPoint(dims::Dims{N}) where N = CenterPoint{N}(CartesianIndex(dims))
CenterPoint(dims::Int64...) = CenterPoint(Tuple(dims))


text:
import QuantumControlBase
using QuantumControlBase: getcontrols, getcontrolderivs, discretize_on_midpoints, evalcontrols
using QuantumControlBase: GradVector, TimeDependentGradGenerator
using QuantumPropagators: init_storage, initpropwrk
using ConcreteStructs

# GRAPE workspace (for internal use)
@concrete terse struct GrapeWrk

    # a copy of the objectives
    objectives

    # the adjoint objectives, containing the adjoint generators for the
    # backward propagation
    adjoint_objectives

    # The kwargs from the control problem
    kwargs

    # Tuple of the original controls (probably functions)
    controls

    # TODO: pulsevals0 and pulsevals1
    pulsevals :: Vector{Float64}

    gradient :: Vector{Float64}  # gradient for guess in iterations

    searchdirection :: Vector{Float64}  # search-direction for guess in iterations

    # Result object
    result

    #################################
    # scratch objects, per objective:

    # backward-propagated states
    # note: storage for fw-propagated states is in result.states
    bw_states

    # forward-propagated states (functional evaluation only)
    fw_states

    # foward-propagated grad-vectors
    fw_grad_states

    # gradients ∂τₖ/ϵₗ(tₙ)
    tau_grads :: Vector{Matrix{ComplexF64}}

    # dynamical generator (normal propagation) at a particular point in time
    G

    # dynamica generator for grad-propagation, time-dependent
    TDgradG

    # dynamical generator for grad-propagation at a particular point in time
    gradG

    control_derivs :: Vector{Vector{Union{Function, Nothing}}}

    vals_dict

    bw_storage # backward storage array (per objective)

    fw_prop_wrk # for normal forward propagation

    bw_prop_wrk # for normal propagation

    grad_prop_wrk # for gradient forward propagation

    use_threads :: Bool

end

function GrapeWrk(problem::QuantumControlBase.ControlProblem)
    use_threads = get(problem.kwargs, :use_threads, false)
    objectives = [obj for obj in problem.objectives]
    adjoint_objectives = [adjoint(obj) for obj in problem.objectives]
    controls = getcontrols(objectives)
    control_derivs = [
        getcontrolderivs(obj.generator, controls) for obj in objectives
    ]
    tlist = problem.tlist
    # interleave the pulse values as [ϵ₁(t̃₁), ϵ₂(t̃₁), ..., ϵ₁(t̃₂), ϵ₂(t̃₂), ...]
    # to allow access as reshape(pulsevals0, L :)[l, n] where l is the control
    # index and n is the time index
    pulsevals = convert(
        Vector{Float64},
        reshape(transpose(hcat(
            [discretize_on_midpoints(control, tlist)
            for control in controls]...
        )), :)
    )
    kwargs = Dict(problem.kwargs)
    pulse_options = problem.pulse_options
    if haskey(kwargs, :continue_from)
        @info "Continuing previous optimization"
        result = kwargs[:continue_from]
        if !(result isa GrapeResult)
            # account for continuing from a different optimization method
            result = convert(GrapeResult, result)
        end
        result.iter_stop = get(problem.kwargs, :iter_stop, 5000)
        result.converged = false
        result.start_local_time = now()
        result.message = "in progress"
        pulsevals = convert(
            Vector{Float64},
            reshape(transpose(hcat(
                [discretize_on_midpoints(control, wrk.result.tlist)
                for control in result.optimized_controls]...
            )), :)
        )
    else
        result = GrapeResult(problem)
    end
    gradient = zeros(length(pulsevals))
    searchdirection = zeros(length(pulsevals))
    bw_states = [similar(obj.initial_state) for obj in objectives]
    fw_states = [similar(obj.initial_state) for obj in objectives]
    fw_grad_states = [
        GradVector(obj.initial_state, length(controls))
        for obj in objectives
    ]
    dummy_vals = IdDict(control => 1.0 for (i, control) in enumerate(controls))
    G = [evalcontrols(obj.generator, dummy_vals) for obj in objectives]
    vals_dict = [copy(dummy_vals) for _ in objectives]
    bw_storage = [init_storage(obj.initial_state, tlist) for obj in objectives]
    fw_prop_method = [
        Val(
            QuantumControlBase.get_objective_prop_method(
                obj, :fw_prop_method, :prop_method; problem.kwargs...
            )
        )
        for obj in objectives
    ]
    bw_prop_method = [
        Val(
            QuantumControlBase.get_objective_prop_method(
                obj, :bw_prop_method, :prop_method; problem.kwargs...
            )
        )
        for obj in objectives
    ]
    grad_prop_method = [
        Val(
            QuantumControlBase.get_objective_prop_method(
                obj, :grad_prop_method, :fw_prop_method, :prop_method;
                problem.kwargs...
            )
        )
        for obj in objectives
    ]

    fw_prop_wrk = [
        QuantumControlBase.initobjpropwrk(obj, tlist, fw_prop_method[k];
                                          initial_state=obj.initial_state,
                                          kwargs...)
        for (k, obj) in enumerate(objectives)
    ]
    bw_prop_wrk = [
        QuantumControlBase.initobjpropwrk(obj, tlist, bw_prop_method[k];
                                          initial_state=obj.initial_state,
                                          kwargs...)
        for (k, obj) in enumerate(objectives)
    ]
    TDgradG = [TimeDependentGradGenerator(obj.generator) for obj in objectives]
    gradG = [evalcontrols(G̃_of_t, dummy_vals) for G̃_of_t ∈ TDgradG]
    grad_prop_wrk = [
        initpropwrk(Ψ̃, tlist, grad_prop_method[k], gradG[k]; kwargs...)
        for (k, Ψ̃) in enumerate(fw_grad_states)
    ]
    tau_grads :: Vector{Matrix{ComplexF64}} = [
        zeros(ComplexF64, length(controls), length(tlist)-1)
        for _ in objectives
    ]
    GrapeWrk(
        objectives,
        adjoint_objectives,
        kwargs,
        controls,
        pulsevals,
        gradient,
        searchdirection,
        result,
        bw_states,
        fw_states,
        fw_grad_states,
        tau_grads,
        G,
        TDgradG,
        gradG,
        control_derivs,
        vals_dict,
        bw_storage,
        fw_prop_wrk,
        bw_prop_wrk,
        grad_prop_wrk,
        use_threads
    )
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct ResourceHandleProto <: ProtoType
    device::AbstractString
    container::AbstractString
    name::AbstractString
    hash_code::UInt64
    maybe_type_name::AbstractString
    ResourceHandleProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ResourceHandleProto

export ResourceHandleProto


text:
function plot_sol(i_exp, sol, exp_data, Tlist, cap, sol0=nothing)
    T0, beta, ocen = l_exp_info[i_exp, :]
    ts = sol.t / 60.0
    ind = length(ts)
    plt = plot(
        exp_data[:, 1] / 60.0,
        exp_data[:, 3],
        seriestype=:scatter,
        label="Exp",
    )

    plot!(
        plt,
        ts,
        sum(clamp.(sol[1:end - 1, :], 0, Inf), dims=1)',
        lw=3,
        legend=:left,
        label="CRNN",
    )

    if sol0 !== nothing
        plot!(
            plt,
            sol0.t / 60,
            sum(sol0[1:end - 1, :], dims=1)',
            label="initial model",
        )
    end
    xlabel!(plt, "Time [min]")
    ylabel!(plt, "Mass")
    title!(plt, cap)
    exp_cond = string(
        @sprintf(
            "T0 = %.1f K \n beta = %.2f K/min \n [O2] = %.2f",
            T0,
            beta,
            exp(ocen) * 100.0
        )
    )
    annotate!(plt, exp_data[end, 1] / 60.0 * 0.85, 0.4, exp_cond)

    plt2 = twinx()
    plot!(
        plt2,
        exp_data[1:ind, 1] / 60,
        Tlist,
        lw=2,
        ls=:dash,
        legend=:topright,
        label="T",
    )
    ylabel!(plt2, "Temp")

    p2 = plot(ts, sol[1, :], lw=2, legend=:right, label="Cellulose")
    for i = 2:ns
        plot!(p2, ts, sol[i, :], lw=2, label="S$i")
    end
    xlabel!(p2, "Time [min]")
    ylabel!(p2, "Mass")

    plt = plot(plt, p2, framestyle=:box, layout=@layout [a; b])
    plot!(plt, size=(800, 800))
    return plt
end

cbi = function (p, i_exp)
    exp_data = l_exp_data[i_exp]
    sol = pred_n_ode(p, i_exp, exp_data)
    Tlist = similar(sol.t)
    T0, beta, ocen = l_exp_info[i_exp, :]
    for (i, t) in enumerate(sol.t)
        Tlist[i] = getsampletemp(t, T0, beta)
    end
    value = l_exp[i_exp]
    plt = plot_sol(i_exp, sol, exp_data, Tlist, "exp_$value")
    png(plt, string(fig_path, "/conditions/pred_exp_$value"))
    return false
end

function plot_loss(l_loss_train, l_loss_val; yscale=:log10)
    plt_loss = plot(l_loss_train, yscale=yscale, label="train")
    plot!(plt_loss, l_loss_val, yscale=yscale, label="val")
    plt_grad = plot(list_grad, yscale=yscale, label="grad_norm")
    xlabel!(plt_loss, "Epoch")
    ylabel!(plt_loss, "Loss")
    xlabel!(plt_grad, "Epoch")
    ylabel!(plt_grad, "Gradient Norm")
    # ylims!(plt_loss, (-Inf, 1e0))
    # ylims!(plt_grad, (-Inf, 1e3))
    plt_all = plot([plt_loss, plt_grad]..., legend=:top, framestyle=:box)
    plot!(
        plt_all,
        size=(1000, 450),
        xtickfontsize=11,
        ytickfontsize=11,
        xguidefontsize=12,
        yguidefontsize=12,
    )
    png(plt_all, string(fig_path, "/loss_grad"))
end

l_loss_train = []
l_loss_val = []
list_grad = []
iter = 1
cb = function (p, loss_train, loss_val, g_norm)
    global l_loss_train, l_loss_val, list_grad, iter
    push!(l_loss_train, loss_train)
    push!(l_loss_val, loss_val)
    push!(list_grad, g_norm)

    if iter % n_plot == 0
        display_p(p)
        list_exp = randperm(n_exp)[1]
        @printf(
            "Min Loss train: %.2e val: %.2e",
            minimum(l_loss_train),
            minimum(l_loss_val)
        )
        println("\n update plot ", l_exp[list_exp], "\n")
        for i_exp in list_exp
            cbi(p, i_exp)
        end

        plot_loss(l_loss_train, l_loss_val; yscale=:log10)

        @save string(ckpt_path, "/mymodel.bson") p opt l_loss_train l_loss_val list_grad iter
    end
    iter += 1
end

if is_restart
    @load string(ckpt_path, "/mymodel.bson") p opt l_loss_train l_loss_val list_grad iter
    iter += 1
end

max_stars_repo_name: DENG-MIT/CRNN-Pyrolysis

max_issues_repo_name: DENG-MIT/CRNN-Pyrolysis

max_forks_repo_name: DENG-MIT/CRNN-Pyrolysis