text:
module Frameworks

using Printf: @printf, @sprintf
using StaticArrays: SVector
using LaTeXStrings: latexstring
using Serialization: serialize
using TimerOutputs: TimerOutputs, TimerOutput, @timeit
using RecipesBase: RecipesBase, @recipe, @series
using ..QuantumOperators: Operator, Operators, Transformation, optype, idtype, identity
using ..Spatials: Bonds, AbstractLattice, acrossbonds, isintracell
using ..DegreesOfFreedom: Hilbert, Table, Term, ismodulatable
using ...Prerequisites: Float, atol, rtol, decimaltostr
using ...Prerequisites.Traits: efficientoperations, reparameter, commontype

import ...Interfaces: id, add!, expand, expand!
import ...Essentials: update, update!, reset!
import ...Prerequisites.Traits: contentnames, getcontent

export Parameters, Boundary, plain
export Engine, AbstractGenerator, Formulation, Entry, CompositeGenerator, Generator, Image, Action, Assignment, Algorithm
export prepare!, run!, save, rundependences!

"""
    Parameters{Names}(values::Number...) where Names

A NamedTuple that contain the key-value pairs.
"""
const Parameters{Names} = NamedTuple{Names, <:Tuple{Vararg{Number}}}
@inline Parameters{Names}(values::Number...) where {Names} = NamedTuple{Names}(values)
@inline @generated function update(params::NamedTuple; parameters...)
    names = fieldnames(params)
    values = Expr(:tuple, [:(get(parameters, $name, getfield(params, $name))) for name in QuoteNode.(names)]...)
    return :(NamedTuple{$names}($values))
end

"""
    match(params₁::Parameters, params₂::Parameters; atol=atol, rtol=rtol) -> Bool

Judge whether the second set of parameters matches the first.
"""
function Base.match(params₁::Parameters, params₂::Parameters; atol=atol, rtol=rtol)
    for name in keys(params₂)
        haskey(params₁, name) && !isapprox(getfield(params₁, name), getfield(params₂, name); atol=atol, rtol=rtol) && return false
    end
    return true
end

"""
    Boundary{Names}(values::AbstractVector{<:Number}, vectors::AbstractVector{<:AbstractVector{<:Number}}) where Names

Boundary twist of operators.
"""
struct Boundary{Names, D<:Number, V<:AbstractVector} <: Transformation
    values::Vector{D}
    vectors::Vector{V}
    function Boundary{Names}(values::AbstractVector{<:Number}, vectors::AbstractVector{<:AbstractVector{<:Number}}) where Names
        @assert length(Names)==length(values)==length(vectors) "Boundary error: mismatched names, values and vectors."
        datatype = promote_type(eltype(values), Float)
        new{Names, datatype, eltype(vectors)}(convert(Vector{datatype}, values), vectors)
    end
end
@inline Base.valtype(::Type{<:Boundary}, M::Type{<:Operator}) = reparameter(M, :value, promote_type(Complex{Int}, valtype(M)))
@inline Base.valtype(B::Type{<:Boundary}, MS::Type{<:Operators}) = (M = valtype(B, eltype(MS)); Operators{M, idtype(M)})

"""
    keys(bound::Boundary) -> Tuple{Vararg{Symbol}}
    keys(::Type{<:Boundary{Names}}) where Names -> Names

Get the names of the boundary parameters.
"""
@inline Base.keys(bound::Boundary) = keys(typeof(bound))
@inline Base.keys(::Type{<:Boundary{Names}}) where Names = Names

"""
    (bound::Boundary)(operator::Operator; origin::Union{AbstractVector, Nothing}=nothing) -> Operator

Get the boundary twisted operator.
"""
@inline function (bound::Boundary)(operator::Operator; origin::Union{AbstractVector, Nothing}=nothing)
    values = isnothing(origin) ? bound.values : bound.values-origin
    return replace(operator, operator.value*exp(1im*mapreduce(u->angle(u, bound.vectors, values), +, id(operator))))
end

"""
    update!(bound::Boundary; parameters...) -> Boundary

Update the values of the boundary twisted phase.
"""
@inline @generated function update!(bound::Boundary; parameters...)
    exprs = []
    for (i, name) in enumerate(QuoteNode.(keys(bound)))
        push!(exprs, :(bound.values[$i] = get(parameters, $name, bound.values[$i])))
    end
    return Expr(:block, exprs..., :(return bound))
end

"""
    merge!(bound::Boundary, another::Boundary) -> typeof(bound)

Merge the values and vectors of the twisted boundary condition from another one.
"""
@inline function Base.merge!(bound::Boundary, another::Boundary)
    @assert keys(bound)==keys(another) "merge! error: mismatched names of boundary parameters."
    bound.values .= another.values
    bound.vectors .= another.vectors
    return bound
end

"""
    Parameters(bound::Boundary)

Get the parameters of the twisted boundary condition.
"""
@inline Parameters(bound::Boundary) = NamedTuple{keys(bound)}(ntuple(i->bound.values[i], Val(fieldcount(typeof(keys(bound))))))

"""
    plain

Plain boundary condition without any twist.
"""
const plain = Boundary{()}(Float[], SVector{0, Float}[])
@inline Base.valtype(::Type{typeof(plain)}, M::Type{<:Operator}) = M
@inline Base.valtype(::Type{typeof(plain)}, M::Type{<:Operators}) = M
@inline (::typeof(plain))(operator::Operator; kwargs...) = operator

"""
    Engine

Abstract type for all engines.

An engine is the core to generate/update/manipulate (representations of) the quantum operators of a quantum lattice system.
"""
abstract type Engine end
@inline Base.:(==)(engine₁::Engine, engine₂::Engine) = ==(efficientoperations, engine₁, engine₂)
@inline Base.isequal(engine₁::Engine, engine₂::Engine) = isequal(efficientoperations, engine₁, engine₂)
@inline Base.repr(engine::Engine) = String(nameof(typeof(engine)))
@inline Base.show(io::IO, engine::Engine) = @printf io "%s" nameof(typeof(engine))
@inline Base.valtype(engine::Engine) = valtype(typeof(engine))
@inline update!(engine::Engine; kwargs...) = error("update! error: not implemented for $(nameof(typeof(engine))).")
@inline Parameters(engine::Engine) = error("Parameters error: not implemented for $(nameof(typeof(engine))).")

"""
    AbstractGenerator <: Engine

Abstract type of the generator of (representations of) the quantum operators of a quantum lattice system.
"""
abstract type AbstractGenerator <: Engine end
@inline Base.eltype(gen::AbstractGenerator) = eltype(typeof(gen))
@inline Base.IteratorSize(::Type{<:AbstractGenerator}) = Base.SizeUnknown()
@inline function Base.iterate(gen::AbstractGenerator)
    ops = expand(gen)
    index = iterate(ops)
    isnothing(index) && return nothing
    return index[1], (ops, index[2])
end
@inline function Base.iterate(gen::AbstractGenerator, state)
    index = iterate(state[1], state[2])
    isnothing(index) && return nothing
    return index[1], (state[1], index[2])
end
@inline Base.show(io::IO, ::MIME"text/latex", gen::AbstractGenerator) = show(io, MIME"text/latex"(), latexstring(latexstring(expand(gen))))

"""
    expand(gen::AbstractGenerator) -> valtype(gen)
    expand!(result, gen::AbstractGenerator) -> typeof(result)

Expand the generator, that is, get the (representations of the) quantum operators of a quantum lattice system (or some part of it).
"""
@inline expand(gen::AbstractGenerator) = expand!(zero(valtype(gen)), gen)
@inline expand!(result, gen::AbstractGenerator) = error("expand! error: not implemented for $(nameof(typeof(gen))).")

"""
    Formulation{F<:Function, P<:Parameters} <: AbstractGenerator

Generator of (representations of the) quantum operators of a quantum lattice system by analytical expressions.
"""
mutable struct Formulation{F<:Function, P<:Parameters} <: AbstractGenerator
    formula::F
    parameters::P
end
@inline function update!(formulation::Formulation; parameters...)
    formulation.parameters = update(formulation.parameters; parameters...)
    return formulation
end
@inline Parameters(formulation::Formulation) = formulation.parameters
@inline (formulation::Formulation)(; kwargs...) = formulation.formula(values(formulation.parameters)...; kwargs...)

"""
    Entry{C, A<:NamedTuple, B<:NamedTuple, P<:Parameters, D<:Boundary} <: AbstractGenerator

An entry of quantum operators (or representations of quantum operators) related to (part of) a quantum lattice system.
"""
mutable struct Entry{C, A<:NamedTuple, B<:NamedTuple, P<:Parameters, D<:Boundary} <: AbstractGenerator
    constops::C
    alterops::A
    boundops::B
    parameters::P
    boundary::D
end
@inline Entry(entry::Entry) = entry
@inline Base.eltype(E::Type{<:Entry}) = eltype(valtype(E))

"""
    Entry(terms::Tuple{Vararg{Term}}, bonds::Bonds, hilbert::Hilbert;
        half::Bool=false,
        table::Union{Nothing, Table}=nothing,
        boundary::Boundary=plain
        )

Construct an entry of quantum operators based on the input terms, bonds, Hilbert space and (twisted) boundary condition.
"""
function Entry(terms::Tuple{Vararg{Term}}, bonds::Bonds, hilbert::Hilbert;
        half::Bool=false,
        table::Union{Nothing, Table}=nothing,
        boundary::Boundary=plain
        )
    constterms, alterterms, choosedterms = termclassifier(terms)
    innerbonds = filter(acrossbonds, bonds, Val(:exclude))
    boundbonds = filter(acrossbonds, bonds, Val(:include))
    constops = Operators{mapreduce(term->optype(typeof(term), typeof(hilbert), eltype(bonds)), promote_type, choosedterms)}()
    map(term->expand!(constops, term, innerbonds, hilbert, half=half, table=table), constterms)
    alterops = NamedTuple{map(id, alterterms)}(map(term->expand(one(term), innerbonds, hilbert, half=half, table=table), alterterms))
    boundops = NamedTuple{map(id, terms)}(
        map(term->map!(
            boundary,
            expand!(Operators{valtype(typeof(boundary), optype(typeof(term), typeof(hilbert), eltype(bonds)))}(),one(term), boundbonds, hilbert, half=half, table=table)),
            terms
        )
    )
    parameters = NamedTuple{map(id, terms)}(map(term->term.value, terms))
    return Entry(constops, alterops, boundops, parameters, boundary)
end
@generated function termclassifier(terms::Tuple{Vararg{Term}})
    constterms, alterterms = [], []
    for (i, term) in enumerate(fieldtypes(terms))
        ismodulatable(term) ? push!(alterterms, :(terms[$i])) : push!(constterms, :(terms[$i]))
    end
    constterms, alterterms = Expr(:tuple, constterms...), Expr(:tuple, alterterms...)
    return Expr(:tuple, constterms, alterterms, (length(constterms.args)>0 ? constterms : alterterms))
end

"""
    (transformation::Transformation)(entry::Entry; kwargs...) -> Entry

Get the transformed entry of (representations of) quantum operators.
"""
function (transformation::Transformation)(entry::Entry; kwargs...)
    wrapper(m) = transformation(m; kwargs...)
    constops = wrapper(entry.constops)
    alterops = NamedTuple{keys(entry.alterops)}(map(wrapper, values(entry.alterops)))
    boundops = NamedTuple{keys(entry.boundops)}(map(wrapper, values(entry.boundops)))
    return Entry(constops, alterops, boundops, entry.parameters, deepcopy(entry.boundary))
end

"""
    valtype(::Type{<:Entry})

Get the valtype of an entry of (representations of) quantum operators.
"""
@inline @generated function Base.valtype(::Type{<:Entry{C, A, B}}) where {C, A<:NamedTuple, B<:NamedTuple}
    optp = C
    (fieldcount(A) > 0) && (optp = reduce(promote_type, fieldtypes(A), init=optp))
    (fieldcount(B) > 0) && (optp = reduce(promote_type, fieldtypes(B), init=optp))
    return optp
end

"""
    expand!(result, entry::Entry) -> typeof(result)

Expand the contents of an entry.
"""
@generated function expand!(result, entry::Entry)
    exprs = [:(add!(result, entry.constops))]
    for name in QuoteNode.(fieldnames(fieldtype(entry, :alterops)))
        push!(exprs, :(add!(result, get(entry.parameters, $name, 1)*getfield(entry.alterops, $name))))
    end
    for name in QuoteNode.(fieldnames(fieldtype(entry, :boundops)))
        push!(exprs, :(add!(result, get(entry.parameters, $name, 1)*getfield(entry.boundops, $name))))
    end
    push!(exprs, :(return result))
    return Expr(:block, exprs...)
end

"""
    update!(entry::Entry{<:Operators}; parameters...) -> Entry

Update the parameters (including the boundary parameters) of an entry of quantum operators.

!!! Note
    The coefficients of `boundops` are also updated due to the change of the boundary parameters.
"""
function update!(entry::Entry{<:Operators}; parameters...)
    entry.parameters = update(entry.parameters; parameters...)
    if !match(Parameters(entry.boundary), NamedTuple{keys(parameters)}(values(parameters)))
        old = copy(entry.boundary.values)
        update!(entry.boundary; parameters...)
        map(ops->map!(op->entry.boundary(op, origin=old), ops), values(entry.boundops))
    end
    return entry
end

"""
    update!(entry::Entry, transformation::Transformation, source::Entry{<:Operators}; kwargs...) -> Entry

Update the contents of an `Entry` by its source entry of quantum operators and the corresponding transformation.
"""
function update!(entry::Entry, transformation::Transformation, source::Entry{<:Operators}; kwargs...)
    entry.parameters = update(entry.parameters; source.parameters...)
    if !match(Parameters(entry.boundary), Parameters(source.boundary))
        update!(entry.boundary; Parameters(source.boundary)...)
        map((dest, ops)->map!(op->transformation(op; kwargs...), empty!(dest), ops), values(entry.boundops), values(source.boundops))
    end
    return entry
end

"""
    Parameters(entry::Entry)

Get the complete set of parameters of an entry of (representations of) quantum operators.
"""
@inline Parameters(entry::Entry) = merge(entry.parameters, Parameters(entry.boundary))

"""
    empty(entry::Entry) -> Entry
    empty!(entry::Entry) -> Entry

Get an empty copy entry or empty an entry of (representations of) quantum operators.
"""
@inline function Base.empty(entry::Entry)
    constops = empty(entry.constops)
    alterops = NamedTuple{keys(entry.alterops)}(map(empty, values(entry.alterops)))
    boundops = NamedTuple{keys(entry.boundops)}(map(empty, values(entry.boundops)))
    return Entry(constops, alterops, boundops, entry.parameters, deepcopy(entry.boundary))
end
@inline function Base.empty!(entry::Entry)
    empty!(entry.constops)
    map(empty!, values(entry.alterops))
    map(empty!, values(entry.boundops))
    return entry
end

"""
    merge(entry::Entry, another::Entry) -> Entry
    merge!(entry::Entry, another::Entry) -> Entry

1) Get a new entry of (representations of) quantum operators by merging two.
2) Merge an entry of (representations of) quantum operators by another one.
"""
@inline Base.merge(entry::Entry, another::Entry) = merge!(empty(entry), another)
function Base.merge!(entry::Entry, another::Entry)
    merge!(entry.constops, another.constops)
    for name in fieldnames(typeof(entry.alterops))
        merge!(getfield(entry.alterops, name), getfield(another.alterops, name))
    end
    for name in fieldnames(typeof(entry.boundops))
        merge!(getfield(entry.boundops, name), getfield(another.boundops, name))
    end
    entry.parameters = another.parameters
    merge!(entry.boundary, another.boundary)
    return entry
end

"""
    reset!(entry::Entry{<:Operators}, terms::Tuple{Vararg{Term}}, bonds::Bonds, hilbert::Hilbert;
        half::Bool=false,
        table::Union{Nothing, Table}=nothing,
        boundary::Boundary=entry.boundary
        ) -> Entry

Reset an entry of quantum operators by the new terms, bonds, Hilbert space and (twisted) boundary condition.
"""
function reset!(entry::Entry{<:Operators}, terms::Tuple{Vararg{Term}}, bonds::Bonds, hilbert::Hilbert;
        half::Bool=false,
        table::Union{Nothing, Table}=nothing,
        boundary::Boundary=entry.boundary
        )
    empty!(entry)
    constterms, alterterms, _ = termclassifier(terms)
    innerbonds = filter(acrossbonds, bonds, Val(:exclude))
    boundbonds = filter(acrossbonds, bonds, Val(:include))
    map(term->expand!(entry.constops, term, innerbonds, hilbert, half=half, table=table), constterms)
    map(term->expand!(getfield(entry.alterops, id(term)), one(term), innerbonds, hilbert, half=half, table=table), alterterms)
    map(term->map!(boundary, expand!(getfield(entry.boundops, id(term)), one(term), boundbonds, hilbert, half=half, table=table)), terms)
    entry.parameters = NamedTuple{map(id, terms)}(map(term->term.value, terms))
    merge!(entry.boundary, boundary)
    return entry
end

"""
    reset!(entry::Entry, transformation::Transformation, source::Entry{<:Operators}; kwargs...)

Reset the contents of an entry by its source entry of quantum operators and the corresponding transformation.
"""
function reset!(entry::Entry, transformation::Transformation, source::Entry{<:Operators}; kwargs...)
    map!(op->transformation(op; kwargs...), empty!(entry.constops), source.constops)
    map((dest, ops)->map!(op->transformation(op; kwargs...), empty!(dest), ops), values(entry.alterops), values(source.alterops))
    map((dest, ops)->map!(op->transformation(op; kwargs...), empty!(dest), ops), values(entry.boundops), values(source.boundops))
    entry.parameters = update(entry.parameters; source.parameters...)
    merge!(entry.boundary, source.boundary)
    return entry
end

"""
    CompositeGenerator{E<:Entry, T<:Union{Table, Nothing}}

Abstract type for composite generator of (representations of) quantum operators.

By protocol, it must have the following predefined contents:
* `operators::E`: the entry for the generated (representations of) quantum operators
* `table::T`: the index-sequence table if it is not nothing
"""
abstract type CompositeGenerator{E<:Entry, T<:Union{Table, Nothing}} <: AbstractGenerator end
@inline contentnames(::Type{<:CompositeGenerator}) = (:operators, :table)
@inline Base.valtype(::Type{<:CompositeGenerator{E}}) where {E<:Entry} = valtype(E)
@inline Base.eltype(::Type{<:CompositeGenerator{E}}) where {E<:Entry} = eltype(E)
@inline expand!(result, gen::CompositeGenerator) = expand!(result, getcontent(gen, :operators))
@inline Parameters(gen::CompositeGenerator) = Parameters(getcontent(gen, :operators))
@inline Entry(gen::CompositeGenerator) = getcontent(gen, :operators)

"""
    Generator{E<:Entry{<:Operators}, TS<:Tuple{Vararg{Term}}, BS<:Bonds, H<:Hilbert, T<:Union{Table, Nothing}} <: CompositeGenerator{E, T}

A generator of operators based on terms, bonds and Hilbert space.
"""
struct Generator{E<:Entry{<:Operators}, TS<:Tuple{Vararg{Term}}, BS<:Bonds, H<:Hilbert, T<:Union{Table, Nothing}} <: CompositeGenerator{E, T}
    operators::E
    terms::TS
    bonds::BS
    hilbert::H
    half::Bool
    table::T
end
@inline contentnames(::Type{<:Generator}) = (:operators, :terms, :bonds, :hilbert, :half, :table)

"""
    Generator(terms::Tuple{Vararg{Term}}, bonds::Bonds, hilbert::Hilbert;
        half::Bool=false,
        table::Union{Table,Nothing}=nothing,
        boundary::Boundary=plain
        )

Construct a generator of operators.
"""
@inline function Generator(terms::Tuple{Vararg{Term}}, bonds::Bonds, hilbert::Hilbert;
        half::Bool=false,
        table::Union{Table,Nothing}=nothing,
        boundary::Boundary=plain
        )
    return Generator(Entry(terms, bonds, hilbert; half=half, table=table, boundary=boundary), terms, bonds, hilbert, half, table)
end

"""
    update!(gen::Generator; parameters...) -> typeof(gen)

Update the coefficients of the terms in a generator.
"""
@inline function update!(gen::Generator; parameters...)
    update!(gen.operators; parameters...)
    map(term->(ismodulatable(term) ? update!(term; parameters...) : term), gen.terms)
    return gen
end

"""
    empty(gen::Generator) -> Generator
    empty!(gen::Generator) -> Generator

1) Get an empty copy of a generator.
2) Empty the `:bonds`, `:hilbert`, `:table` and `:operators` attributes of a generator.
"""
@inline function Base.empty(gen::Generator)
    Generator(
        empty(gen.operators),
        gen.terms,
        empty(gen.bonds),
        empty(gen.hilbert),
        gen.half,
        isnothing(gen.table) ? nothing : empty(gen.table),
        )
end
function Base.empty!(gen::Generator)
    empty!(gen.bonds)
    empty!(gen.hilbert)
    isnothing(gen.table) || empty!(gen.table)
    empty!(gen.operators)
    return gen
end

"""
    reset!(gen::Generator, lattice::AbstractLattice, boundary::Boundary=gen.operators.boundary) -> Generator

Reset a generator by a new lattice and a new twisted boundary condition.
"""
function reset!(gen::Generator, lattice::AbstractLattice, boundary::Boundary=gen.operators.boundary)
    reset!(gen.bonds, lattice)
    reset!(gen.hilbert, lattice.pids)
    isnothing(gen.table) || reset!(gen.table, gen.hilbert)
    reset!(gen.operators, gen.terms, gen.bonds, gen.hilbert, half=gen.half, table=gen.table, boundary=boundary)
    return gen
end

"""
    expand(gen::Generator, name::Symbol) -> Operators
    expand(gen::Generator, i::Int) -> Operators
    expand(gen::Generator, name::Symbol, i::Int) -> Operators

Expand the operators of a generator:
1) the operators of a specific term;
2) the operators on a specific bond;
3) the operators of a specific term on a specific bond.
"""
function expand(gen::Generator, name::Symbol)
    result = zero(valtype(gen))
    term = get(gen.terms, Val(name))
    if !ismodulatable(term)
        expand!(result, term, filter(acrossbonds, gen.bonds, Val(:exclude)), gen.hilbert, half=gen.half, table=gen.table)
    else
        for opt in getfield(gen.operators.alterops, name)
            add!(result, opt*term.value)
        end
    end
    for opt in getfield(gen.operators.boundops, name)
        add!(result, opt*term.value)
    end
    return result
end
function expand(gen::Generator, i::Int)
    bond = gen.bonds[i]
    result = zero(valtype(gen))
    map(term->expand!(result, term, bond, gen.hilbert, half=gen.half, table=gen.table), gen.terms)
    isintracell(bond) || for opt in result
        result[id(opt)] = gen.operators.boundary(opt)
    end
    return result
end
function expand(gen::Generator, name::Symbol, i::Int)
    bond = gen.bonds[i]
    term = get(gen.terms, Val(name))
    result = expand!(zero(valtype(gen)), term, bond, gen.hilbert, half=gen.half, table=gen.table)
    isintracell(bond) || for opt in result
        result[id(opt)] = gen.operators.boundary(opt)
    end
    return result
end
@inline @generated function Base.get(terms::Tuple{Vararg{Term}}, ::Val{Name}) where Name
    i = findfirst(isequal(Name), map(id, fieldtypes(terms)))::Int
    return :(terms[$i])
end

"""
    Image{E<:Entry, H<:Transformation, T<:Union{Table, Nothing}} <: CompositeGenerator{E, T}

The image of a transformation on an generator of (representation of) quantum operators.
"""
mutable struct Image{E<:Entry, H<:Transformation, T<:Union{Table, Nothing}} <: CompositeGenerator{E, T}
    operators::E
    transformation::H
    table::T
    sourceid::UInt
end
@inline contentnames(::Type{<:Image}) = (:operators, :transformation, :table, :sourceid)

"""
    (transformation::Transformation)(gen::AbstractGenerator; table::Union{Table, Nothing}=nothing, kwargs...) -> Image

Get the image of a transformation on an instance of `AbstractGenerator`.
"""
@inline function (transformation::Transformation)(gen::AbstractGenerator; table::Union{Table, Nothing}=nothing, kwargs...)
    return Image(transformation(Entry(gen); kwargs...), transformation, table, objectid(gen))
end

"""
    empty(gen::Image) -> Image
    empty!(gen::Image) -> Image

1) Get an empty copy of the image of a transformation on a generator.
2) Empty the `:operators` and `table` attributes of the image of a transformation on a generator.
"""
@inline Base.empty(gen::Image) = Image(empty(gen.operators), gen.transformation, isnothing(gen.table) ? nothing : empty(gen.table), gen.sourceid)
@inline function Base.empty!(gen::Image)
    empty!(gen.operators)
    isnothing(gen.table) || empty!(gen.table)
    return gen
end

"""
    update!(gen::Image; parameters...) -> typeof(gen)

Update the parameters of the image of a transformation on a generator.
"""
@inline function update!(gen::Image; parameters...)
    update!(gen.operators; parameters...)
    return gen
end

"""
    update!(gen::Image, source::AbstractGenerator; kwargs...) -> Image

Update the image of an transformation by its updated source generator.
"""
@inline function update!(gen::Image, source::AbstractGenerator; kwargs...)
    @assert gen.sourceid==objectid(source) "update! error: mismatched simplified generator, transformation and source generator."
    update!(gen.operators, gen.transformation, Entry(source); kwargs...)
    return gen
end

"""
    reset!(gen::Image, source::AbstractGenerator;
        transformation::Transformation=gen.transformation,
        table::Union{Table, Nothing}=nothing,
        kwargs...
        ) -> Image

Reset the image of a transformation on a generator.
"""
@inline function reset!(gen::Image, source::AbstractGenerator;
        transformation::Transformation=gen.transformation,
        table::Union{Table, Nothing}=nothing,
        kwargs...
        )
    @assert gen.sourceid==objectid(source) "reset! error: mismatched image and source generator."
    reset!(gen.operators, transformation, Entry(source); kwargs...)
    gen.transformation = transformation
    isnothing(table) && (table = isa(source, CompositeGenerator) ? getcontent(source, :table) : gen.table)
    isnothing(table) || gen.table===table || merge!(empty!(gen.table), table)
    return gen
end

"""
    Action

Abstract type for all actions.
"""
abstract type Action end
@inline Base.:(==)(action₁::Action, action₂::Action) = ==(efficientoperations, action₁, action₂)
@inline Base.isequal(action₁::Action, action₂::Action) = isequal(efficientoperations, action₁, action₂)
@inline prepare!(action::Action, engine::Engine) = nothing
@inline update!(action::Action; parameters...) = action

"""
    Assignment{A<:Action, P<:Parameters, M<:Function, D<:Tuple{Vararg{Symbol}}, R} <: Function

An assignment associated with an action.
"""
mutable struct Assignment{A<:Action, P<:Parameters, M<:Function, D<:Tuple{Vararg{Symbol}}, R} <: Function
    id::Symbol
    action::A
    parameters::P
    map::M
    dependences::D
    data::R
    ismatched::Bool
    save::Bool
end
@inline Base.:(==)(assign₁::Assignment, assign₂::Assignment) = ==(efficientoperations, assign₁, assign₂)
@inline Base.isequal(assign₁::Assignment, assign₂::Assignment) = isequal(efficientoperations, assign₁, assign₂)
@inline Parameters(assignment::Assignment) = assignment.parameters

"""
    valtype(assign::Assignment)
    valtype(::Type{<:Assignment})

The type of the data(result) of an assignment.
"""
@inline Base.valtype(assign::Assignment) = valtype(typeof(assign))
@inline Base.valtype(::Type{<:Assignment{<:Action, <:Parameters, <:Function, <:Tuple{Vararg{Symbol}}, R}}) where R = R

"""
    update!(assign::Assignment; parameters...) -> Assignment

Update the parameters of an assignment and the status of its associated action.
"""
function update!(assign::Assignment; parameters...)
    if length(parameters)>0
        assign.parameters = update(assign.parameters; parameters...)
        update!(assign.action; assign.map(assign.parameters)...)
    end
    return assign
end

"""
    Algorithm{E<:Engine, P<:Parameters, M<:Function} <: Function

An algorithm associated with an engine.
"""
mutable struct Algorithm{E<:Engine, P<:Parameters, M<:Function} <: Function
    name::Symbol
    engine::E
    din::String
    dout::String
    parameters::P
    map::M
    assignments::Dict{Symbol, Assignment}
    timer::TimerOutput
end
@inline run!(alg::Algorithm, assign::Assignment) = nothing
@inline function Base.:(==)(algorithm₁::Algorithm, algorithm₂::Algorithm)
    return ==(
        (algorithm₁.name, algorithm₁.engine, algorithm₁.din, algorithm₁.dout, algorithm₁.parameters, algorithm₁.map, algorithm₁.assignments),
        (algorithm₂.name, algorithm₂.engine, algorithm₂.din, algorithm₂.dout, algorithm₂.parameters, algorithm₂.map, algorithm₂.assignments),
    )
end
@inline function Base.isequal(algorithm₁::Algorithm, algorithm₂::Algorithm)
    return isequal(
        (algorithm₁.name, algorithm₁.engine, algorithm₁.din, algorithm₁.dout, algorithm₁.parameters, algorithm₁.map, algorithm₁.assignments),
        (algorithm₂.name, algorithm₂.engine, algorithm₂.din, algorithm₂.dout, algorithm₂.parameters, algorithm₂.map, algorithm₂.assignments),
    )
end
function Base.show(io::IO, alg::Algorithm)
    @printf io "%s(%s)" alg.name alg.engine
    for (name, value) in pairs(alg.parameters)
        @printf io "_%s" decimaltostr(value, 10)
    end
end
@inline Parameters(algorithm::Algorithm) = algorithm.parameters

"""
    Algorithm(name::Symbol, engine::Engine; din::String=".", dout::String=".", parameters::Parameters=Parameters(engine), map::Function=identity)

Construct an algorithm.
"""
@inline function Algorithm(name::Symbol, engine::Engine; din::String=".", dout::String=".", parameters::Parameters=Parameters(engine), map::Function=identity)
    return Algorithm(name, engine, din, dout, parameters, map, Dict{Symbol, Assignment}(), TimerOutput())
end

"""
    update!(alg::Algorithm; parameters...) -> Algorithm

Update the parameters of an algorithm and its associated engine.
"""
function update!(alg::Algorithm; parameters...)
    if length(parameters)>0
        alg.parameters = update(alg.parameters; parameters...)
        update!(alg.engine; alg.map(alg.parameters)...)
    end
    return alg
end

"""
    repr(alg::Algorithm, f::Function=param->true; ndecimal::Int=10) -> String

Get the repr representation of an algorithm.

Optionally, some parameters of the algorithm can be filtered by specifying the `f` function.
Besides, the maximum number of decimals of the parameters can also be specified.
"""
function Base.repr(alg::Algorithm, f::Function=param->true; ndecimal::Int=10)
    result = [@sprintf "%s(%s)" alg.name alg.engine]
    for (name, value) in pairs(alg.parameters)
        f(name) && push!(result, decimaltostr(value, ndecimal))
    end
    return join(result, "_")
end

"""
    summary(alg::Algorithm)

Provide a summary of an algorithm.
"""
function Base.summary(alg::Algorithm)
    @info "Summary of $(alg.name)($(nameof(typeof(alg.engine)))):"
    @info string(alg.timer)
end

"""
    nameof(alg::Algorithm, assign::Assignment) -> String

Get the name of the combination of an algorithm and an assignment.
"""
@inline Base.nameof(alg::Algorithm, assign::Assignment) = @sprintf "%s_%s" repr(alg) assign.id

"""
    add!(alg::Algorithm, id::Symbol, action::Action; parameters::Parameters=Parameters{()}(), kwargs...) -> Tuple{Algorithm, Assignment}

Add an assignment on a algorithm by providing the contents of the assignment without the execution of it.
"""
function add!(alg::Algorithm, id::Symbol, action::Action; parameters::Parameters=Parameters{()}(), kwargs...)
    parameters = merge(alg.parameters, parameters)
    map = get(kwargs, :map, identity)
    dependences = get(kwargs, :dependences, ())
    data = prepare!(action, alg.engine)
    save = get(kwargs, :save, false)
    assign = Assignment(id, action, parameters, map, dependences, data, false, save)
    alg.assignments[id] = assign
    return (alg, assign)
end

"""
    (alg::Algorithm)(assign::Assignment) -> Tuple{Algorithm, Assignment}
    (assign::Assignment)(alg::Algorithm) -> Tuple{Algorithm, Assignment}

Run an assignment based on an algorithm.

The difference between these two methods is that the first uses the parameters of `assign` as the current parameters while the second uses those of `alg`.
"""
function (alg::Algorithm)(assign::Assignment)
    ismatched = match(assign.parameters, alg.parameters)
    if !(assign.ismatched && ismatched)
        ismatched || update!(alg; assign.parameters...)
        run!(alg, assign)
        assign.ismatched = true
    end
    return (alg, assign)
end
function (assign::Assignment)(alg::Algorithm)
    ismatched = match(alg.parameters, assign.parameters)
    if !(assign.ismatched && ismatched)
        ismatched || update!(assign; alg.parameters...)
        run!(alg, assign)
        assign.ismatched = true
    end
    return (alg, assign)
end

"""
    (alg::Algorithm)(id::Symbol, action::Action; info::Bool=true, parameters::Parameters=Parameters{()}(), kwargs...) -> Tuple{Algorithm, Assignment}

Add an assignment on a algorithm by providing the contents of the assignment, and run this assignment.
"""
@inline function (alg::Algorithm)(id::Symbol, action::Action; info::Bool=true, parameters::Parameters=Parameters{()}(), kwargs...)
    add!(alg, id, action; parameters=parameters, kwargs...)
    alg(id, info=info)
end

"""
    (alg::Algorithm)(id::Symbol; info::Bool=true, parameters::Parameters=Parameters{()}()) -> Tuple{Algorithm, Assignment}

Run an assignment specified by its id.

Optionally, the time of the run process can be informed by setting the `info` argument to be `true`.
"""
function (alg::Algorithm)(id::Symbol; info::Bool=true, parameters::Parameters=Parameters{()}())
    assign = alg.assignments[id]
    update!(assign; parameters...)
    @timeit alg.timer string(id) alg(assign)
    info && @info "Action $id($(nameof(assign.action|>typeof))): time consumed $(TimerOutputs.time(alg.timer[string(id)]) / 10^9)s."
    assign.save && save(alg, assign)
    return (alg, assign)
end

"""
    save(alg::Algorithm, assign::Assignment) -> Tuple{Algorithm, Assignment}

Save the data of an assignment registered on an algorithm.
"""
function save(alg::Algorithm, assign::Assignment)
    serialize(@sprintf("%s/%s.dat", alg.dout, nameof(alg, assign)), data)
    return (alg, assign)
end

"""
    rundependences!(alg::Algorithm, assign::Assignment, f::Function=assign->true) -> Tuple{Algorithm, Assignment}

Run the dependences of an assignment.

Optionally, some dependences can be filtered by specifying the `f` function.
"""
function rundependences!(alg::Algorithm, assign::Assignment, f::Function=assign->true)
    for id in assign.dependences
        dependence = alg.assignments[id]
        f(dependence) && dependence(alg)
    end
    return (alg, assign)
end

"""
    @recipe plot(pack::Tuple{Algorithm, Assignment})

Define the recipe for the visualization of an assignment of an algorithm.
"""
@recipe function plot(pack::Tuple{Algorithm, Assignment})
    title --> nameof(pack[1], pack[2])
    titlefontsize --> 10
    legend --> false
    seriestype --> (isa(pack[2].data, Tuple{Any, Any, Any}) ? :heatmap : :path)
    pack[2].data
end

end  # module


text:
export StingerConfig, stingerconfig, generateconfig

"""
Used to set the configuration of the STINGER data structure.
"""
immutable StingerConfig
    nv::Int64
    nebs::Int64
    netypes::Int64
    nvtypes::Int64
    memory_size::Int64
    no_map_none_etype::Int8
    no_map_none_vtype::Int8
    no_resize::Int8
end

"""
    stingerconfig(
        nv, nebs=0, netypes=0, nvtypes=0, memory_size=0, no_map_none_etype=0,
        no_map_none_vtype=0, no_resize=0)

Creates a `StingerConfig`.
"""
function stingerconfig(
    nv::Int64; nebs::Int64=0, netypes::Int64=0, nvtypes::Int64=0,
    memory_size::Int64=0, no_map_none_etype::Int64=0, no_map_none_vtype::Int64=0,
    no_resize::Int64=0)
    StingerConfig(nv, nebs, netypes, nvtypes, memory_size, no_map_none_etype, no_map_none_vtype, no_resize)
end

function maxmemsize()
    return ccall(dlsym(stinger_core_lib, "stinger_max_memsize"), UInt64, ())
end

function verticessize(nv::Int64)
    return ccall(dlsym(stinger_core_lib, "stinger_vertices_size"), UInt64, (Int64,), nv)
end

function physmapsize(nv::Int64)
    return ccall(dlsym(stinger_core_lib, "stinger_physmap_size"), UInt64, (Int64,), nv)
end

function namessize(ntypes::Int64)
    return ccall(dlsym(stinger_core_lib, "stinger_names_size"), UInt64, (Int64,), ntypes)
end

"""
    generateconfig(nv::Int64; netypes::Int64=1, nvtypes::Int64=1)

Generates a config for the specified number and types of vertices and the number
of edge types. The generated config attempts to maximize the number of edge
blocks that can be allocated.
"""
function generateconfig(nv::Int64; netypes::Int64=1, nvtypes::Int64=1)
    sz = maxmemsize() * 0.5 - (verticessize(nv)  + physmapsize(nv) +
        namessize(netypes) + namessize(nvtypes) + 16 + 16)
    nebs = floor(Int, sz / (512 + 8))
    return stingerconfig(nv, nebs=nebs, netypes=netypes, nvtypes=nvtypes, no_resize=1)
end

ext: jl


text:
using MacroTools
using Logging
using Base.Threads

export @trace

abstract type RFTContext end

struct RFTFunc <: RFTContext end
struct RFTCall <: RFTContext end
struct RFTParam <: RFTContext end
struct RFTKw <: RFTContext end

rewrite_for_trace(def::Expr) = rewrite_for_trace(RFTFunc(), def)

rewrite_for_trace(ctx::RFTContext, def::Expr) =
    # Dispatch on def.head
    rewrite_for_trace(ctx, Val(def.head), def)

rewrite_for_trace(ctx::RFTFunc, head::Val{:function}, def::Expr) =
    Expr(:quote, rewrite_for_trace(RFTCall(), def.args[1]))

rewrite_for_trace(ctx::RFTCall, head::Val{:call}, expr::Expr) =
    Expr(:call, map(expr.args) do param
             rewrite_for_trace(RFTParam(), param)
         end...)

rewrite_for_trace(ctx::RFTParam, pexp::Symbol) = Expr(:$, pexp)

rewrite_for_trace(ctx::RFTContext, head::Val{Symbol("::")}, pexp::Expr) =
    rewrite_for_trace(ctx, pexp.args[1])  # Discard the specializer

rewrite_for_trace(ctx::RFTParam, head::Val{:parameters}, pexp::Expr) =
    Expr(:parameters, map(pexp.args) do p
             rewrite_for_trace(ctx, p)
             end...)

rewrite_for_trace(ctx::RFTParam, head::Val{:kw}, pexp::Expr) =
    Expr(:kw,
         rewrite_for_trace(RFTKw(), pexp.args[1]),
         Expr(:$, pexp.args[1]))

rewrite_for_trace(ctx::RFTKw, kw::Symbol) = kw

const TRACE_ENTER = "Trace enter"
const TRACE_EXIT = "Trace exit"

trace_enter(id, fcall) =
    @info(TRACE_ENTER,
          id=id,
          threadid=Threads.threadid(),
          fcall=string(fcall))

trace_exit(id, result) =
    @info(TRACE_EXIT,
          id=id,
          threadid=Threads.threadid(),
          result=result)

_trace_counter = Threads.Atomic{Int}(1)

"""
    @trace(global_flag, definition)
Cause the call arguments and return values of the function defined by
`definition` to be logged if `global_flag` is true at run time.
`definition` should define a method.
"""
macro trace(global_flag, definition)
    definition = longdef(definition)
    pieces = splitdef(definition)
    result = gensym("result")
    this_trace = gensym("this_trace")
    bodyfunction = gensym("bodyfunction")
    Expr(:escape,
         Expr(definition.head,
              definition.args[1],   # function signature
              Expr(:block,
                   Expr(:function,
                        Expr(:call, bodyfunction),
                        pieces[:body]),
                   Expr(:if, global_flag,
                        Expr(:block,
                             Expr(Symbol("="), this_trace,
                                  Expr(:call, Base.Threads.atomic_add!,
                                       :(NahaJuliaLib._trace_counter), 1)),
                             Expr(:call, :(NahaJuliaLib.trace_enter), this_trace,
                                  rewrite_for_trace(definition)))),
                   Expr(:(=), result, Expr(:call, bodyfunction)),
                   Expr(:if, global_flag,
                        Expr(:call, :(NahaJuliaLib.trace_exit), this_trace, result)),
                   result)))
end


text:
function ispossemidef(A)
	eigs = eigvals(A)
	if any(real(eigs) .< 0)
		return false
	else
		return true
	end
end

struct NotImplemented <: Exception
	fun::Symbol
	type::Symbol
end

@inline NotImplemented(fun::Symbol, type::DataType) = NotImplemented(fun, Symbol(DataType))
@inline NotImplemented(fun::Symbol, type) = NotImplemented(fun, Symbol(typeof(type)))

Base.showerror(io::IO, ex::NotImplemented) =
	print(io, "Not Implemented Error: ", ex.fun, " not implemented for type ", ex.type)

function gen_zinds(n::Int, m::Int, N::Int, isequal::Bool=false)
	Nu = isequal ? N : N-1
	zinds = [(k-1)*(n+m) .+ (1:n+m) for k = 1:Nu]
	if !isequal
		push!(zinds, (N-1)*(n+m) .+ (1:n))
	end
	return zinds
end


text:
function reconstruct_obs!(g::GGLRM, XY::Array{VecOrMat{Float64}}; X = g.X, Y = g.Y)
  yidxs = get_yidxs(g.losses)
  obsex = g.observed_examples
  for j in 1:length(g.losses)
    @inbounds Yj = Y[:,yidxs[j]]
    @inbounds Xisj = X[:, obsex[j]]
    #Is a column of the XY matrix proper
    At_mul_B!(XY[j], Xisj, Yj)
  end
end

function reconstruct_obs(g::GGLRM)
  yidxs = get_yidxs(g.losses)
  obsex = g.observed_examples
  XY = VecOrMat{Float64}[zeros(length(obsex[j]), length(yidxs[j]))
                                        for j in 1:length(obsex)]
  for j in 1:length(XY)
    if size(XY[j], 2) == 1
      XY[j] = vec(XY[j])
    end
  end
  reconstruct_obs!(g, XY)
  XY
end

#Evaluates the loss function over all the observed values, optimized for few observations
function sparse_loss_objective(g::GGLRM, XY::Array{VecOrMat{Float64}})
  yidxs = get_yidxs(g.losses)
  obj = zeros(Threads.nthreads())
  Threads.@threads for j in 1:length(g.losses)
    obsex = g.observed_examples[j]
    #A may be a DataFrame or DataArray, but native Arrays are fastest
    @inbounds Aj = convert(Array, g.A[obsex, j])
    #This assumes a reconstruction over only the observations as done in reconstruct_obs!
    @inbounds XYj = XY[j]
    obj[Threads.threadid()] += evaluate(g.losses[j], XYj, Aj)
  end
  sum(obj)
end

#Calculates the whole objective of the GLRM
function sparse_whole_objective(g::GGLRM, XY::Array{VecOrMat{Float64}};
                        X::Matrix{Float64}=g.X, Y::Matrix{Float64}=g.Y)
  sparse_loss_objective(g, XY) + evaluate(g.rx, X) + evaluate(g.ry, Y)
end

#Makes some performance optimizations
#Finds all of the gradients in a column so that the size of the column gradient is consistent
#as opposed to the row gradient which has column-chunks and whatnot
@inline function _threadedupdateGradX!(g::AbstractGLRM, XY::Array{VecOrMat{Float64}}, gx::Matrix{Float64})
  yidxs = get_yidxs(g.losses)
  scale!(gx,0)

  #Update the gradient, go by column then by row
  Threads.@threads for j in 1:length(g.losses)
    #Yj for computing gradient
    @inbounds Yj = view(g.Y, :, yidxs[j])
    obsex = g.observed_examples[j]

    #Take whole columns of XY and A and take the gradient of those
    @inbounds Aj = convert(Array, g.A[obsex, j])
    @inbounds XYj = XY[j]
    grads = grad(g.losses[j], XYj, Aj)

    #Single dimensional losses
    if isa(grads, Vector)
      for e in 1:length(obsex)
        #i = obsex[e], so update that portion of gx
        @inbounds gxi = view(gx, :, obsex[e])
        axpy!(grads[e], Yj, gxi)
      end
    else
      for e in 1:length(obsex)
        @inbounds gxi = view(gx, :, obsex[e])
        gemm!('N','N',1.0,Yj, grads[e,:], 1.0, gxi)
      end
    end
  end
end

@inline function _threadedupdateGradY!(g::AbstractGLRM, XY::Array{VecOrMat{Float64}}, gy::Matrix{Float64})
  yidxs = get_yidxs(g.losses)
  #scale y gradient to zero
  scale!(gy, 0)

  #Update the gradient
  Threads.@threads for j in 1:length(g.losses)
    @inbounds gyj = view(gy, :, yidxs[j])
    obsex = g.observed_examples[j]
    #Take whole columns of XY and A and take the gradient of those
    @inbounds Aj = convert(Array, g.A[obsex, j])
    @inbounds XYj = XY[j]
    grads = grad(g.losses[j], XYj, Aj)
    #Single dimensional losses
    if isa(grads, Vector)
      for e in 1:length(obsex)
        #i = obsex[e], so use that for Xi
        @inbounds Xi = view(g.X, :, obsex[e])
        axpy!(grads[e], Xi, gyj)
      end
    else
      for e in 1:length(obsex)
        @inbounds Xi = view(g.X, :, obsex[e])
        gemm!('N','T',1.0, Xi, grads[e,:], 1.0, gyj)
      end
    end
  end
end

#Does a line search for the step size for X, returns the new step size
@inline function _threadedproxStepX!(g::AbstractGLRM, params::ProxGradParams,
                            newX::Matrix{Float64}, gx::Matrix{Float64},
                            XY::Array{VecOrMat{Float64}}, newXY::Array{VecOrMat{Float64}},
                            αx::Number)
  #l = 1.5
  l = maximum(map(length, g.observed_features))+1#(mapreduce(length,+,g.observed_features) + 1)

  obj = sparse_loss_objective(g,XY) + evaluate(g.rx, g.X)
  newobj = NaN
  while αx > params.min_stepsize #Linesearch to find the new step size
    stepsize = αx/l

    axpy!(-stepsize, gx, newX)
    prox!(g.rx, newX, stepsize)
    reconstruct_obs!(g, newXY, X=newX)#At_mul_B!(newXY, newX, g.Y)
    newobj = sparse_loss_objective(g, newXY) + evaluate(g.rx, newX)
    #newobj = threaded_objective(g, newXY)
    if newobj < obj
      copy!(g.X, newX)
      αx *= 1.03#(1.05)
      break
    else #Try again with smaller step-size
      copy!(newX, g.X)
      αx *= 0.9 #0.7
      if αx < params.min_stepsize
        αx = params.min_stepsize * 1.1
        break
      end #if
    end #if else
  end #while
  αx, newobj
end

@inline function _threadedproxStepY!(g::AbstractGLRM, params::ProxGradParams,
                              newY::Matrix{Float64}, gy::Matrix{Float64},
                              XY::Array{VecOrMat{Float64}}, newXY::Array{VecOrMat{Float64}},
                              αy::Number)
  #l = 1.5
  l = maximum(map(length, g.observed_examples)) + 1#(mapreduce(length,+,g.observed_features) + 1)
  #obj = threaded_objective(g,XY)
  obj = sparse_loss_objective(g, XY) + evaluate(g.ry, g.Y)
  newobj = NaN
  while αy > params.min_stepsize #Linesearch to find the new step size
    stepsize = αy/l
    axpy!(-stepsize, gy, newY)
    prox!(g.ry, newY, stepsize)
    reconstruct_obs!(g, newXY, Y=newY)#At_mul_B!(newXY, g.X, newY)
    newobj = sparse_loss_objective(g, newXY) + evaluate(g.ry, newY)
    #newobj = threaded_objective(g, newXY)
    if newobj < obj
      copy!(g.Y, newY)
      αy *= 1.03#1.05
      break
    else #Try again with smaller step-size
      copy!(newY, g.Y)
      αy *= 0.9 #0.7
      if αy < params.min_stepsize
        αy = params.min_stepsize * 1.1
        break
      end #if
    end #if else
  end #while
  αy, newobj
end

function fit_sparse!(g::GGLRM,
                      params::ProxGradParams=ProxGradParams(),
                      ch::ConvergenceHistory=ConvergenceHistory("ProxGradGLRM"),
                      verbose=true)
  X,Y = g.X, g.Y
  A = g.A
  losses, rx, ry = g.losses, g.rx, g.ry
  yidxs = get_yidxs(g.losses)

  #Initialize X*Y
  XY = reconstruct_obs(g)

  tm = 0
  update_ch!(ch, tm, sparse_whole_objective(g,XY))
  (objx, objy) = NaN, NaN
  #Step sizes
  αx = params.stepsize
  αy = params.stepsize

  # stopping criterion: stop when decrease in objective < tol, scaled by the number of observations
  scaled_abs_tol = params.abs_tol * mapreduce(length,+,g.observed_features)

  newX = copy(X)
  newY = copy(Y)
  newXY = deepcopy(XY) #It's a ragged array of matrices, so deepcopy is needed

  #Gradient matrices
  gx = zeros(X)
  gy = zeros(Y)

  if verbose println("Fitting GGLRM") end
  for t in 1:params.max_iter
    #X update-----------------------------------------------------------------
    for xit in 1:params.inner_iter_X
      #_threadedupdateGradX!(g,XY,gx)
      _threadedupdateGradX!(g,XY,gx)
      #Take a prox step with line search
      αx, objx = _threadedproxStepX!(g, params, newX, gx, XY, newXY, αx)
      reconstruct_obs!(g, XY) #Get the new XY matrix for objective
    end
    #Y Update---------------------------------------------------------------
    for yit in 1:params.inner_iter_Y
      #_threadedupdateGradY!(g,XY,gy)
      _threadedupdateGradY!(g,XY,gy)
      αy, objy = _threadedproxStepY!(g, params, newY, gy, XY, newXY, αy)
      reconstruct_obs!(g, XY) #Get the new XY matrix for objective
    end
    if t % 10 == 0
      if verbose
        println("Iteration $t, objective value: $(objy + evaluate(rx, g.X))")
      end
    end
    #Update convergence history
    obj = objy + evaluate(rx, g.X)
    tm = time() - tm
    update_ch!(ch, tm, obj)
    tm = time()
    #Check stopping criterion
    obj_decrease = ch.objective[end-1] - obj
    if t>10 && (obj_decrease < scaled_abs_tol || obj_decrease/obj < params.rel_tol)
      if verbose
        println("Iteration $t, objective value: $obj")
      end
      break
    end
  end #For
  ch
end


text:
using TypeClasses

"""
Iterables can be seen two ways. On the one hand, an iterable is mainly defined by its `iterate` method, which can be
thought of as a kind of TypeClass (similar to how `map`, `combine`, `Monad`, or `Monoid` refer to TypeClasses). In
this sense being iterable is regarded as a decisive characteristic which is usually checked via `Base.isiterable`.
As an example, the TypeClass FlipTypes has a default implementation which uses exactly this.

The alternativ semantics, which is implemented in this file, is that an iterable type can be seen as an (abstract)
DataType on top of which we can define TypeClasses themselves.

Because of this duality, we don't define TypeClasses for iterables directly by dispatching on `isiterable`, but we
provide a custom wrapper `TypeClasses.Iterable` on which all TypeClasses are defined.
I.e. if you want to use standard TypeClasses on top of your iterable Type, just wrap it within an `Iterable`:
```
myiterable = ...
Iterable(myiterable)
```
"""

# MonoidAlternative
# =================

TypeClasses.neutral(iter::Type{<:Iterable}) = Iterable()
TypeClasses.combine(it1::Iterable, it2::Iterable) = Iterable(chain(it1.iter, it2.iter))

# FunctorApplicativeMonad
# =======================

TypeClasses.pure(::Type{<:Iterable}, a) = Iterable(TypeClasses.DataTypes.Iterables.IterateSingleton(a))
TypeClasses.ap(fs::Iterable, it::Iterable) = Iterable(f(a) for f ∈ fs.iter for a ∈ it.iter)

# Base.map(f, iter1, iter2, iter3...) is already defined and uses zip semantics, hence we don't overload it

# Flattening Iterables works with everything being iterable itself (it is treated as iterable)
TypeClasses.flatten(it::Iterable) = Iterable(Iterators.flatten(it.iter))
TypeClasses.flatmap(f, it::Iterable) = flatten(map(f, it))

# FlipTypes
# =========

# we define flip_types for all iterable despite it only works if the underlying element defines `ap`
# as there is no other sensible definition for Iterable, an error that the element does not implement `ap`
# is actually the correct error
flip_types(iter::Iterable) = default_flip_types_having_pure_combine_apEltype(iter)


text:
macro linklibrary_div(n, es)
  path = normpath(Pkg.dir("FastSigmoid"),"c-src","libfastposit.so")
  fnname = QuoteNode(Symbol(:p, n, :e, es, :_div_e))

  posittype = Symbol(:UInt, n)

  esc(:(pdiv = (a, b) -> ccall( ($fnname, $path), $posittype, ($posittype, $posittype), a, b )))
end

#comprehensive integration testing
@linklibrary_div(8,0)
for x = 0x00:0xFF, y = 0x00:0xFF
  xinf = (x == 0x80)
  yinf = (y == 0x80)
  xzer = (x == 0x00)
  yzer = (y == 0x00)
  if !((xinf & yinf) | (xzer & yzer))
    @test (x,y,pdiv(x, y)) == (x,y,reinterpret(UInt8, (reinterpret(Posit{8,0}, x) / reinterpret(Posit{8,0}, y))))
  end
end

@linklibrary_div(8,1)
for x = 0x00:0xFF, y = 0x00:0xFF
  xinf = (x == 0x80)
  yinf = (y == 0x80)
  xzer = (x == 0x00)
  yzer = (y == 0x00)
  if !((xinf & yinf) | (xzer & yzer))
    @test (x,y,pdiv(x, y)) == (x,y,reinterpret(UInt8, (reinterpret(Posit{8,1}, x) / reinterpret(Posit{8,1}, y))))
  end
end

@linklibrary_div(8,2)
for x = 0x00:0xFF, y = 0x00:0xFF
  xinf = (x == 0x80)
  yinf = (y == 0x80)
  xzer = (x == 0x00)
  yzer = (y == 0x00)
  if !((xinf & yinf) | (xzer & yzer))
    @test (x,y,pdiv(x, y)) == (x,y,reinterpret(UInt8, (reinterpret(Posit{8,2}, x) / reinterpret(Posit{8,2}, y))))
  end
end

macro test_big_bittype_div(N, ES)
  quote
    @linklibrary_div($N, $ES)

    println("testing division on p$($N)e$($ES)")

    UType = UIntLookup[$N]

    for idx = 1:10000
      x = rand(UType)
      y = rand(UType)

      xinf = (x == topbits($N))
      yinf = (y == topbits($N))
      xzer = (x == zero($N))
      yzer = (y == zero($N))

      if !((xinf & yinf) | (xzer & yzer)) #ignore the cases that doesn't work.
        @test (x,y,pdiv(x, y)) == (x,y,reinterpret(UType, (reinterpret(Posit{$N,$ES}, x) / reinterpret(Posit{$N,$ES}, y))))
      end
    end
  end
end

@test_big_bittype_div(16,0)
@test_big_bittype_div(16,1)
@test_big_bittype_div(16,2)

@test_big_bittype_div(32,0)
@test_big_bittype_div(32,1)
@test_big_bittype_div(32,2)
@test_big_bittype_div(32,3)


text:
using Sunaasobi
using Test

@testset "Sunaasobi.jl" begin
    # Write your tests here.
end


text:
###
### Bipartite graph utilities
###

"""
    maximal_matching(s::SystemStructure, eqfilter=eq->true, varfilter=v->true) -> Matching

Find equation-variable maximal bipartite matching. `s.graph` is a bipartite graph.
"""
BipartiteGraphs.maximal_matching(s::SystemStructure, eqfilter=eq->true, varfilter=v->true) =
    maximal_matching(s.graph, eqfilter, varfilter)

function error_reporting(sys, bad_idxs, n_highest_vars, iseqs)
    io = IOBuffer()
    if iseqs
        error_title = "More equations than variables, here are the potential extra equation(s):\n"
        out_arr = equations(sys)[bad_idxs]
    else
        error_title = "More variables than equations, here are the potential extra variable(s):\n"
        out_arr = structure(sys).fullvars[bad_idxs]
    end

    Base.print_array(io, out_arr)
    msg = String(take!(io))
    neqs = length(equations(sys))
    if iseqs
        throw(ExtraEquationsSystemException(
            "The system is unbalanced. "
            * "There are $n_highest_vars highest order derivative variables "
            * "and $neqs equations.\n"
            * error_title
            * msg
        ))
    else
        throw(ExtraVariablesSystemException(
            "The system is unbalanced. "
            * "There are $n_highest_vars highest order derivative variables "
            * "and $neqs equations.\n"
            * error_title
            * msg
        ))
    end
end

###
### Structural check
###
function check_consistency(sys::AbstractSystem)
    s = structure(sys)
    @unpack graph, var_to_diff, fullvars = s
    n_highest_vars = count(v->length(outneighbors(s.var_to_diff, v)) == 0, vertices(s.var_to_diff))
    neqs = nsrcs(graph)
    is_balanced = n_highest_vars == neqs

    if neqs > 0 && !is_balanced
        varwhitelist = var_to_diff .== nothing
        var_eq_matching = maximal_matching(graph, eq->true, v->varwhitelist[v]) # not assigned
        # Just use `error_reporting` to do conditional
        iseqs = n_highest_vars < neqs
        if iseqs
            eq_var_matching = invview(complete(var_eq_matching)) # extra equations
            bad_idxs = findall(isnothing, @view eq_var_matching[1:nsrcs(graph)])
        else
            bad_idxs = findall(isequal(unassigned), var_eq_matching)
        end
        error_reporting(sys, bad_idxs, n_highest_vars, iseqs)
    end

    # This is defined to check if Pantelides algorithm terminates. For more
    # details, check the equation (15) of the original paper.
    extended_graph = (@set graph.fadjlist = Vector{Int}[graph.fadjlist; map(collect, edges(var_to_diff))])
    extended_var_eq_matching = maximal_matching(extended_graph)

    unassigned_var = []
    for (vj, eq) in enumerate(extended_var_eq_matching)
        if eq === unassigned
            push!(unassigned_var, fullvars[vj])
        end
    end

    if !isempty(unassigned_var) || !is_balanced
        io = IOBuffer()
        Base.print_array(io, unassigned_var)
        unassigned_var_str = String(take!(io))
        errmsg = "The system is structurally singular! " *
                 "Here are the problematic variables: \n" *
                 unassigned_var_str
        throw(InvalidSystemException(errmsg))
    end

    return nothing
end

###
### BLT ordering
###

"""
    find_var_sccs(g::BipartiteGraph, assign=nothing)

Find strongly connected components of the variables defined by `g`. `assign`
gives the undirected bipartite graph a direction. When `assign === nothing`, we
assume that the ``i``-th variable is assigned to the ``i``-th equation.
"""
function find_var_sccs(g::BipartiteGraph, assign=nothing)
    cmog = DiCMOBiGraph{true}(g, Matching(assign === nothing ? Base.OneTo(nsrcs(g)) : assign))
    sccs = Graphs.strongly_connected_components(cmog)
    foreach(sort!, sccs)
    return sccs
end

function sorted_incidence_matrix(sys, val=true; only_algeqs=false, only_algvars=false)
    var_eq_matching, var_scc = algebraic_variables_scc(sys)
    s = structure(sys)
    @unpack fullvars, graph = s
    g = graph
    varsmap = zeros(Int, ndsts(graph))
    eqsmap = zeros(Int, nsrcs(graph))
    varidx = 0
    eqidx = 0
    for vs in scc, v in vs
        eq = var_eq_matching[v]
        if eq !== unassigned
            eqsmap[eq] = (eqidx += 1)
            varsmap[var] = (varidx += 1)
        end
    end
    for i in diffvars_range(s)
        varsmap[i] = (varidx += 1)
    end
    for i in dervars_range(s)
        varsmap[i] = (varidx += 1)
    end
    for i in 1:nsrcs(graph)
        if eqsmap[i] == 0
            eqsmap[i] = (eqidx += 1)
        end
    end

    I = Int[]
    J = Int[]
    for eq in 𝑠vertices(g)
        only_algeqs && (isalgeq(s, eq) || continue)
        for var in 𝑠neighbors(g, eq)
            only_algvars && (isalgvar(s, var) || continue)
            i = eqsmap[eq]
            j = varsmap[var]
            (iszero(i) || iszero(j)) && continue
            push!(I, i)
            push!(J, j)
        end
    end
    #sparse(I, J, val, nsrcs(g), ndsts(g))
    sparse(I, J, val)
end

###
### Structural and symbolic utilities
###

function find_solvables!(sys)
    s = structure(sys)
    @unpack fullvars, graph, solvable_graph = s
    eqs = equations(sys)
    empty!(solvable_graph)
    for (i, eq) in enumerate(eqs)
        isdiffeq(eq) && continue
        term = value(eq.rhs - eq.lhs)
        for j in 𝑠neighbors(graph, i)
            isalgvar(s, j) || continue
            var = fullvars[j]
            isinput(var) && continue
            a, b, islinear = linear_expansion(term, var)
            a = unwrap(a)
            if islinear && (!(a isa Symbolic) && a isa Number && a != 0)
                add_edge!(solvable_graph, i, j)
            end
        end
    end
    s
end

# debugging use
function reordered_matrix(sys, torn_matching)
    s = structure(sys)
    @unpack graph = s
    eqs = equations(sys)
    nvars = ndsts(graph)
    max_matching = complete(maximal_matching(s))
    torn_matching = complete(torn_matching)
    sccs = find_var_sccs(s.graph, max_matching)
    I, J = Int[], Int[]
    ii = 0
    M = Int[]
    solved = BitSet(findall(torn_matching .!== unassigned))
    for vars in sccs
        append!(M, filter(in(solved), vars))
        append!(M, filter(!in(solved), vars))
    end
    M = invperm(vcat(M, setdiff(1:nvars, M)))
    for vars in sccs
        e_solved = [torn_matching[v] for v in vars if torn_matching[v] !== unassigned]
        for es in e_solved
            isdiffeq(eqs[es]) && continue
            ii += 1
            js = [M[x] for x in 𝑠neighbors(graph, es) if isalgvar(s, x)]
            append!(I, fill(ii, length(js)))
            append!(J, js)
        end

        e_residual = setdiff([max_matching[v] for v in vars if max_matching[v] !== unassigned], e_solved)
        for er in e_residual
            isdiffeq(eqs[er]) && continue
            ii += 1
            js = [M[x] for x in 𝑠neighbors(graph, er) if isalgvar(s, x)]
            append!(I, fill(ii, length(js)))
            append!(J, js)
        end
    end
    # only plot algebraic variables and equations
    sparse(I, J, true)
end

###
### Nonlinear equation(s) solving
###

@noinline nlsolve_failure(rc) = error("The nonlinear solver failed with the return code $rc.")

function numerical_nlsolve(f, u0, p)
    prob = NonlinearProblem{false}(f, u0, p)
    sol = solve(prob, NewtonRaphson())
    rc = sol.retcode
    rc === :DEFAULT || nlsolve_failure(rc)
    # TODO: robust initial guess, better debugging info, and residual check
    sol.u
end


text:
function within{T <: Real}(scenario::Vector{T}, gamma::Gamma, p::DeploymentProblem)
    # returns true if scenario ∈ gamma; false otherwise
    if any(scenario .> gamma._single)
        return false
    elseif sum(scenario) > gamma._global
        return false
    else
        for j in 1:p.nregions
            if sum(scenario[vec(p.adjacency[j,:])]) > gamma._local[j]
                return false
            end
        end
        for i in 1:p.nlocations
            if sum(scenario[p.coverage[:,i]]) > gamma._regional[i]
                return false
            end
        end
    end
    true
end

type Evaluation
    objvalue::Int
    shortfall::Vector{Int}
    dispatch::Matrix{Int}
end

function evaluate{T1, T2 <: Real}(x::Vector{T1}, scenario::Vector{T2}, p::DeploymentProblem;
    solver=GurobiSolver(OutputFlag=0,PrePasses=3))

    I = 1:p.nlocations
    J = 1:p.nregions

    m = JuMP.Model(solver=solver)
    JuMP.@variable(m, y[1:p.nlocations,1:p.nregions] >= 0, Int)
    JuMP.@variable(m, z[1:p.nregions] >= 0, Int)

    JuMP.@objective(m, Min, sum(z[j] for j in J))

    # flow constraints at each station
    for i in I
        JuMP.@expression(m, outflow, sum(y[i,j] for j in filter(j->p.coverage[j,i], J)))
        JuMP.@constraint(m, x[i] >= outflow)
    end
    # shortfall from satisfying demand/calls
    for j in J
        JuMP.@expression(m, inflow, sum(y[i,j] for i in filter(i->p.coverage[j,i], I)))
        JuMP.@constraint(m, z[j] >= scenario[j] - inflow)
    end
    status = JuMP.solve(m)
    @assert status == :Optimal

    Evaluation(Int(JuMP.getObjectiveValue(m)),
               [round(Int,z) for z in JuMP.getValue(z)],
               map(y->round(Int,y), JuMP.getValue(y)))
end

function evaluate{T <: Real}(x::Vector{T}, p::DeploymentProblem)
    Evaluation[evaluate(x,vec(p.demand[i,:]),p) for i in p.test]
end

type Result
    stoch_x
    robust_x
    bad_scenarios
    upperbounds
    lowerbounds
    stoch_timing
    robust_timing
end

function performance(p::DeploymentProblem, α::Float64; verbose=false)
    stoch_model = StochasticDeployment(p)
    tic()
    solve(stoch_model)
    stoch_timing = toq()
    robust_model = RobustDeployment(p, α=α, verbose=verbose, master_verbose=verbose)
    tic()
    solve(robust_model, p, verbose=verbose)
    robust_timing = toq()
    Result( deployment(stoch_model),
            deployment(robust_model),
            robust_model.scenarios,
            robust_model.upperbounds,
            robust_model.lowerbounds,
            stoch_timing,
            robust_timing)
end

function test_performance(p::DeploymentProblem; namb::Vector{Int}=[25:5:45], alpha::Vector{Float64}=[0.1,0.05,0.01,0.001,0.0001], verbose=false)
    results = Array(Result, (5,5))
    for (i,n) in enumerate(namb)
        print(n)
        p.nambulances = n
        stoch_model = StochasticDeployment(p)
        tic()
        solve(stoch_model)
        stoch_timing = toq()
        for (j,α) in enumerate(alpha)
            print(" $α")
            robust_model = RobustDeployment(p, α=α, verbose=verbose, master_verbose=verbose)
            tic()
            solve(robust_model, p, verbose=verbose)
            robust_timing = toq()
            results[i,j] = Result(  deployment(stoch_model),
                                    deployment(robust_model),
                                    robust_model.scenarios,
                                    robust_model.upperbounds,
                                    robust_model.lowerbounds,
                                    stoch_timing,
                                    robust_timing)
        end
        println("")
    end
    results
end


text:
using Documenter
using NoiseRobustDifferentiation

makedocs(;
    modules=[NoiseRobustDifferentiation],
    sitename="NoiseRobustDifferentiation.jl",
    authors="Adrian Hill",
    format=Documenter.HTML(; prettyurls=get(ENV, "CI", "false") == "true", assets=String[]),
    pages=["Home" => "index.md", "Examples" => "examples.md"],
)

deploydocs(;
    repo="github.com/adrhill/NoiseRobustDifferentiation.jl.git",
    devbranch="main",
)


text:
# a lightweight CUDA array type for testing purposes

## ctor & finalizer

mutable struct CuTestArray{T,N}
    ptr::CuPtr{T}
    shape::NTuple{N,Int}
    function CuTestArray{T,N}(shape::NTuple{N,Int}) where {T,N}
        len = prod(shape)
        buf = Mem.alloc(Mem.Device, len*sizeof(T))
        ptr = convert(CuPtr{T}, buf)

        obj = new{T,N}(ptr, shape)
        finalizer(obj) do a
            CUDAdrv.isvalid(buf.ctx) && Mem.free(buf)
        end
        return obj
    end
    function CuTestArray{T,N}(ptr::CuPtr{T}, shape::NTuple{N,Int}) where {T,N}
        new{T,N}(ptr, shape)
    end
end

Base.pointer(xs::CuTestArray) = xs.ptr

Base.length(xs::CuTestArray) = prod(xs.shape)

Base.unsafe_convert(::Type{<:CuPtr}, x::CuTestArray) = pointer(x)

## memory copy operations

function CuTestArray(src::Array{T,N}) where {T,N}
    dst = CuTestArray{T,N}(size(src))
    unsafe_copyto!(pointer(dst), pointer(src), length(src))
    return dst
end

function Base.Array(src::CuTestArray{T,N}) where {T,N}
    dst = Array{T,N}(undef, src.shape)
    unsafe_copyto!(pointer(dst), pointer(src), length(src))
    return dst
end

## conversions

using Adapt
function Adapt.adapt_storage(::CUDAnative.Adaptor, a::CuTestArray{T,N}) where {T,N}
    ptr = pointer(a)
    devptr = CUDAnative.DevicePtr{T,AS.Global}(ptr)
    CuDeviceArray{T,N,AS.Global}(a.shape, devptr)
end


text:
using QuasiCopula, LinearAlgebra, GLM
using Random, Distributions, DataFrames, ToeplitzMatrices
using Test, BenchmarkTools

BLAS.set_num_threads(1)
Threads.nthreads()

p = 3    # number of fixed effects, including intercept

# true parameter values
Random.seed!(12345)
βtrue = rand(Uniform(-2, 2), p)
# βtrue = 0.1 * ones(p)
σ2true = [0.5]
ρtrue = [0.5]

function get_V(ρ, n)
    vec = zeros(n)
    vec[1] = 1.0
    for i in 2:n
        vec[i] = ρ
    end
    V = ToeplitzMatrices.SymmetricToeplitz(vec)
    V
end

#simulation parameters
samplesize = 10000

st = time()
currentind = 1
d = Poisson()
link = LogLink()
D = typeof(d)
Link = typeof(link)
T = Float64

gcs = Vector{GLMCopulaCSObs{T, D, Link}}(undef, samplesize)

ni = 25 #  number of observations per individual
V = get_V(ρtrue[1], ni)

# true Gamma
Γ = σ2true[1] * V

# for reproducibility I will simulate all the design matrices here
Random.seed!(12345)
X_samplesize = [randn(ni, p - 1) for i in 1:samplesize]

for i in 1:samplesize
    X = [ones(ni) X_samplesize[i]]
    η = X * βtrue
    μ = exp.(η)
    vecd = Vector{DiscreteUnivariateDistribution}(undef, ni)
    for i in 1:ni
        vecd[i] = Poisson(μ[i])
    end
    nonmixed_multivariate_dist = NonMixedMultivariateDistribution(vecd, Γ)
    # simuate single vector y
    y = Vector{Float64}(undef, ni)
    res = Vector{Float64}(undef, ni)
    rand(nonmixed_multivariate_dist, y, res)
    # push!(Ystack, y)
    V = [Float64.(Matrix(I, ni, ni))]
    # V = [ones(ni, ni)]
    gcs[i] = GLMCopulaCSObs(y, X, d, link)
end

# form model
gcm = GLMCopulaCSModel(gcs)
# precompile
println("precompiling Poisson CS fit")
gcm2 = deepcopy(gcm);
QuasiCopula.fit!(gcm2, IpoptSolver(print_level = 0, max_iter = 20));

fittime = @elapsed QuasiCopula.fit!(gcm, IpoptSolver(print_level = 5, max_iter = 100, tol = 10^-8, accept_after_max_steps = 2, limited_memory_max_history = 50, hessian_approximation = "limited-memory"))

loglikelihood!(gcm, true, true)
@show fittime
@show gcm.β
@show gcm.σ2
@show gcm.ρ
@show gcm.∇β
@show gcm.∇σ2
@show gcm.∇ρ

vcov!(gcm)
@show QuasiCopula.confint(gcm)
mseβ, mseρ, mseσ2 = MSE(gcm, βtrue, ρtrue, σ2true)
@show mseβ
@show mseσ2
@show mseρ

@test mseβ < 0.01
@test mseσ2 < 1
@test mseρ < 1

println("checking memory allocation for Poisson CS")
# logl_gradient_memory = @benchmark loglikelihood!($gcm, true, false) # this will allocate for threads
logl_gradient_memory = @benchmark loglikelihood!($gcm.data[1], $gcm.β, $gcm.ρ[1], $gcm.σ2[1], true, false)
@test logl_gradient_memory.memory == 0.0


text:
# fallback methods
# ----------------

set_format_options!(docformat::WeaveFormat; _kwargs...) = return

function restore_header!(doc)
    (hasproperty(doc.format, :preserve_header) && doc.format.preserve_header) || return

    # only strips Weave headers
    delete!(doc.header, WEAVE_OPTION_NAME)
    if haskey(doc.header, WEAVE_OPTION_NAME_DEPRECATED)
        @warn "Weave: `options` key is deprecated. Use `weave_options` key instead." _id = WEAVE_OPTION_DEPRECATE_ID maxlog = 1
        delete!(doc.header, WEAVE_OPTION_NAME_DEPRECATED)
    end
    isempty(doc.header) && return

    # restore remained headers as `DocChunk`
    header_text = "---\n$(YAML.write(doc.header))---"
    pushfirst!(doc.chunks, DocChunk(header_text, 0, 0))
end

render_chunk(docformat::WeaveFormat, chunk::DocChunk) = join((render_inline(c) for c in chunk.content))

render_inline(inline::InlineText) = inline.content

function render_inline(inline::InlineCode)
    isempty(inline.rich_output) || return inline.rich_output
    isempty(inline.figures) || return inline.figures[end]
    return inline.output
end

function render_chunk(docformat::WeaveFormat, chunk::CodeChunk)

    # Fill undefined options with format specific defaults
    isnothing(chunk.options[:out_width]) && (chunk.options[:out_width] = docformat.out_width)
    isnothing(chunk.options[:fig_pos]) && (chunk.options[:fig_pos] = docformat.fig_pos)

    # Only use floats if chunk has caption or sets fig_env
    if !isnothing(chunk.options[:fig_cap]) && isnothing(chunk.options[:fig_env])
        (chunk.options[:fig_env] = docformat.fig_env)
    end

    hasproperty(docformat, :indent) && (chunk.content = indent(chunk.content, docformat.indent))

    chunk.content = render_code(docformat, chunk.content)

    if !chunk.options[:eval]
        return if chunk.options[:echo]
            string(docformat.codestart, '\n', chunk.content, docformat.codeend)
        else
            ""
        end
    end

    if chunk.options[:term]
        result = render_termchunk(docformat, chunk)
    else
        result = if chunk.options[:echo]
            # Convert to output format and highlight (html, tex...) if needed
            string(docformat.codestart, chunk.content, docformat.codeend, '\n')
        else
            ""
        end

        if (strip(chunk.output) ≠ "" || strip(chunk.rich_output) ≠ "") &&
           (chunk.options[:results] ≠ "hidden")
            if chunk.options[:results] ≠ "markup" && chunk.options[:results] ≠ "hold"
                strip(chunk.output) ≠ "" && (result *= "$(chunk.output)\n")
                strip(chunk.rich_output) ≠ "" && (result *= "$(chunk.rich_output)\n")
            else
                if chunk.options[:wrap]
                    chunk.output =
                        '\n' * wraplines(chunk.output, chunk.options[:line_width])
                    chunk.output = render_output(docformat, chunk.output)
                else
                    chunk.output = '\n' * rstrip(chunk.output)
                    chunk.output = render_output(docformat, chunk.output)
                end

                hasproperty(docformat, :indent) && (chunk.output = indent(chunk.output, docformat.indent))

                strip(chunk.output) ≠ "" && (
                    result *= "$(docformat.outputstart)$(chunk.output)\n$(docformat.outputend)\n"
                )
                strip(chunk.rich_output) ≠ "" && (result *= chunk.rich_output * '\n')
            end
        end
    end

    # Handle figures
    if chunk.options[:fig] && length(chunk.figures) > 0
        result *= render_figures(docformat, chunk)
    end

    return result
end

render_code(docformat::WeaveFormat, code) = code

indent(text, nindent) = join(map(x -> string(repeat(' ', nindent), x), split(text, '\n')), '\n')

function wraplines(text, line_width = 75)
    result = AbstractString[]
    lines = split(text, '\n')
    for line in lines
        if length(line) > line_width
            push!(result, wrapline(line, line_width))
        else
            push!(result, line)
        end
    end

    return strip(join(result, '\n'))
end

function wrapline(text, line_width = 75)
    result = ""
    while length(text) > line_width
        result *= first(text, line_width) * '\n'
        text = chop(text, head = line_width, tail = 0)
    end
    result *= text
end

render_output(docformat::WeaveFormat, output) = output

function render_termchunk(docformat::WeaveFormat, chunk)
    return if should_render(chunk)
        string(docformat.termstart, chunk.output, '\n', docformat.termend, '\n')
    else
        ""
    end
end

should_render(chunk) = chunk.options[:echo] && chunk.options[:results] ≠ "hidden"

render_doc(docformat, body, doc) = body

# utilities
# ---------

function clear_buffer_and_format!(io::IOBuffer, out::IOBuffer, render_function)
    text = take2string!(io)
    m = Markdown.parse(text, flavor = WeaveMarkdown.weavemd)
    write(out, string(render_function(m)))
end

addlines(op, inline) = inline.ctype === :line ? string('\n', op, '\n') : op

get_mustache_template(path::AbstractString) = Mustache.template_from_file(path)
get_mustache_template(tpl::Mustache.MustacheTokens) = tpl

get_highlight_stylesheet(mime, highlight_theme) =
    get_highlight_stylesheet(mime, get_highlight_theme(highlight_theme))
get_highlight_stylesheet(mime, highlight_theme::Type{<:Highlights.AbstractTheme}) =
    sprint((io, x) -> Highlights.stylesheet(io, mime, x), highlight_theme)

get_highlight_theme(::Nothing) = Highlights.Themes.DefaultTheme
get_highlight_theme(highlight_theme::Type{<:Highlights.AbstractTheme}) = highlight_theme

highlight_code(mime, code, highlight_theme) =
    highlight(mime, strip(code), Highlights.Lexers.JuliaLexer, highlight_theme)
highlight_term(mime, output, highlight_theme) =
    highlight(mime, strip(output), Highlights.Lexers.JuliaConsoleLexer, highlight_theme)
highlight(mime, output, lexer, theme = Highlights.Themes.DefaultTheme) =
    sprint((io, x) -> Highlights.highlight(io, mime, x, lexer, theme), output)


text:
using Revise
using NNFEM
using PyPlot
using LinearAlgebra
using ADCME
using ADCMEKit
using DelimitedFiles

NT = 100
Δt = 1/NT

n = 10
m = 2n
h = 1/n

# Create a very simple mesh
elements = FiniteStrainContinuum[]
prop = Dict("name"=> "PlaneStressIncompressibleRivlinSaunders", "rho"=> 1.522,  "C1"=>0.162, "C2"=>5.9e-3)
coords = zeros((m+1)*(n+1), 2)
for j = 1:n
    for i = 1:m
        idx = (m+1)*(j-1)+i
        elnodes = [idx; idx+1; idx+1+m+1; idx+m+1]
        ngp = 3
        nodes = [
            (i-1)*h (j-1)*h
            i*h (j-1)*h
            i*h j*h
            (i-1)*h j*h
        ]
        coords[elnodes, :] = nodes
        push!(elements, FiniteStrainContinuum(nodes, elnodes, prop, ngp))
    end
end

# impose force on the top
Edge_Traction_Data = Array{Int64}[]
for i = 1:m
  elem = elements[i]
  for k = 1:4
    if elem.coords[k,2]<0.001 && elem.coords[k+1>4 ? 1 : k+1,2]<0.001
      push!(Edge_Traction_Data, [i, k, 1])
    end
  end
end
Edge_Traction_Data = hcat(Edge_Traction_Data...)'|>Array

# fixed on the bottom, push on the right
EBC = zeros(Int64, (m+1)*(n+1), 2)
FBC = zeros(Int64, (m+1)*(n+1), 2)
g = zeros((m+1)*(n+1), 2)
f = zeros((m+1)*(n+1), 2)
for j = 1:n+1
    idx = (j-1)*(m+1) + m+1
    EBC[idx,:] .= -1 # time-independent boundary, right
end
for i = 1:m+1
    idx = n*(m+1) + i
    EBC[idx,:] .= -1 # fixed boundary, bottom
end

dimension = 2
domain = Domain(coords, elements, dimension, EBC, g, FBC, f, Edge_Traction_Data)

# Set initial condition: at reset
Dstate = zeros(domain.neqs)
state = zeros(domain.neqs)
velo = zeros(domain.neqs)
acce = zeros(domain.neqs)
EBC_func = nothing
FBC_func = nothing
Body_func = nothing
# Construct Edge_func
function Edge_func_hyperelasticity(x, y, t, idx)
  return [zeros(length(x)) 0.1*ones(length(x)) * sin(2π*t)]
end

globaldata = GlobalData(state, Dstate, velo, acce, domain.neqs, EBC_func, FBC_func,Body_func, Edge_func_hyperelasticity)

assembleMassMatrix!(globaldata, domain)
SolverInitial!(Δt, globaldata, domain)

updateStates!(domain, globaldata)
ω = EigenMode(Δt, globaldata, domain)
@show "stable time step is ", 0.8 * 2/ω, " current time step is ", Δt
for i = 1:NT
    @info i
    # global globaldata, domain = GeneralizedAlphaSolverStep(globaldata, domain, Δt)
    global globaldata, domain = ExplicitSolverStep(globaldata, domain, Δt)
end
d_ = hcat(domain.history["state"]...)'|>Array

# p = visualize_displacement(d_, domain)
# saveanim(p, "hyperelasticity_exp.gif")

close("all")
  plot(d_[:,1], "-", color="C1")
  plot(d_[:,1+domain.nnodes], color="C2")


text:
function dotcheck(T,n,tol)
    x = rand(T,n)
    y = rand(T,n)
    xb = big.(x)
    yb = big.(y)
    db = dot(xb,yb)
    e = eps(real(T))
    sxy = norm(xb)*norm(yb)
    @test abs(big(dot(x,y)) - db) / sxy < tol * n * e

    # check for correct handling of strides
    xs = zeros(T,2*n) .+ T(NaN)
    ys = zeros(T,3*n) .+ T(NaN)
    for j=1:n
        xs[2*j]=x[j]
        ys[3*j]=y[j]
    end
    vx = view(xs,2:2:2*n)
    vy = view(ys,3:3:3*n)
    @test abs(big(dot(x,vy)) - db) / sxy < tol * n * e
    @test abs(big(dot(vx,y)) - db) / sxy < tol * n * e
    @test abs(big(dot(vx,vy)) - db) / sxy < tol * n * e
end

@testset "dot $T" for T in (Double64, Double32, Complex{Double64})
    n = 67
    dotcheck(T,n,2)
end

function axpycheck(T,n,tol)
    x = rand(T,n)
    y = rand(T,n)
    xb = big.(x)
    yb = big.(y)
    if T <: Complex
        a = T(1.5 + 2.5im)
        ab = big(1.5) + 1im*big(2.5)
    else
        a = T(1.5)
        ab = big(1.5)
    end
    zb = LinearAlgebra.axpy!(ab,xb,yb)
    e = eps(real(T))
    @test norm(big.(axpy!(a,x,y)) - zb) / norm(zb) < tol * n * e
end

@testset "axpy $T" for T in (Double64, Double32, Complex{Double64})
    n = 67
    axpycheck(T,n,2)
end

function scalecheck(T,n,tol)
    x = rand(T,n)
    xb = big.(x)
    if T <: Complex
        a = T(1.5 + 2.5im)
        ab = big(1.5) + 1im*big(2.5)
    else
        a = T(1.5)
        ab = big(1.5)
    end
    e = eps(real(T))
    zb = copy(xb)
    LinearAlgebra.lmul!(ab,zb)
    z = copy(x)
    @test norm(big.(lmul!(a,z)) - zb) / norm(zb) < tol * n * e
    zb = copy(xb)
    LinearAlgebra.rmul!(zb,ab)
    z = copy(x)
    @test norm(big.(rmul!(z,a)) - zb) / norm(zb) < tol * n * e
end

@testset "scale $T" for T in (Double64, Double32, Complex{Double64})
    n = 67
    scalecheck(T,n,2)
end

max_stars_repo_name: RalphAS/DoubleBLAS


text:
"""
Collection of utilities for creating, validating and sending HTTP requests
"""
module Requests

using Base: @kwdef
using Chain: @chain
using Dates: Date, Day, today, days, unix2datetime
using HTTP: HTTP
using JSON3: JSON3
using Logging: @debug, @warn

using Stonks: UpdatableSymbol, RequestBuilderError, APIResponseError
using Stonks: split_tickers_in_batches, get_minimum_dates
using Stonks.APIClients: APIResource, APIClient
using Stonks.Models: AbstractStonksRecord, AssetInfo, AssetPrice, ExchangeRate
using Stonks.Parsers: AbstractContentParser, parse_content

@kwdef struct RequestParams
  url::String
  tickers::Vector{UpdatableSymbol}
  retries::Int = 0
  headers::Dict{String,String} = Dict()
  params::Dict{String,String} = Dict()
  query::String = ""
  from::Union{Date,Missing} = missing
  to::Union{Date,Missing} = missing
end

"""
    prepare_requests(tickers, resource; kwargs...)

  Decides how many requests should be created based upon `resource.max_batch_size`.
  After this step, for each request, the parameters with {param} are replaced by the value of the keyword argumets with the same name.
"""
function prepare_requests(
  tickers::Vector{UpdatableSymbol}, resource::APIResource; kwargs...
)::Union{Vector{RequestParams},Exception}
  requests = @chain begin
    tickers
    split_tickers_in_batches(_, resource.max_batch_size)
    map(t -> resolve_request_parameters(t, resource; kwargs...), _)
  end
  valid_requests = @chain requests filter(r -> isa(r, RequestParams), _)
  isempty(valid_requests) && return first(requests)
  length(valid_requests) != length(requests) &&
    return first(filter(r -> isa(r, RequestBuilderError), requests))
  #@debug "Prepared $(length(valid_requests)) requests for $(resource.url)), params -> $(map(r -> r.params, valid_requests))"
  return valid_requests
end

#"""
#    prepare_requests(tickers, provider; kwargs...)
#
#Prepares request for all all the API resources. It will run `prepare_requests(tickers, resource)` for each key of `provider.resources`.
#"""
# function prepare_requests(tickers::Vector{UpdatableSymbol}, provider::APIClient; kwargs...)::Union{Dict{String,Vector{RequestParams}}, Exception}
#   requests = Dict(k => prepare_requests(tickers, resource; kwargs...) for (k, resource) in provider.resources)
#   # Requests: Dict{String, Vector{RequestParams}}
#   valid_requests = Dict()
#   for (k, reqs) in requests
#     if all([r -> isa(r, RequestParams) for r in reqs])
#        valid_requests[k] = reqs
#     end
#   end
#   #valid_requests = Dict(k => reqs for (k, reqs) in requests if all(map(r -> isa(r, RequestParams), reqs)))
#   println(valid_requests)
#   if isempty(valid_requests)
#     @warn "No valid requests"
#   end
#   first_failure = first([ first(reqs) for (k, reqs) in requests if any(map(r -> isa(r, RequestBuilderError), reqs))])
#   isempty(valid_requests) && return first_failure
#   diff = setdiff(keys(requests), keys(valid_requests))
#   if !isempty(diff)
#     @warn("The following resources have errors and were removed: $(join(diff, ','))")
#   end
#   _req_per_resource = [k => length(v) for (k, v) in valid_requests]
#   _total_req = sum([v for (k, v) in _req_per_resource])
#   @debug "Prepared a total of $(_total_req) requests: $(_req_per_resource)"
#   return requests
# end

function materialize_request(rp::RequestParams, parser::AbstractContentParser; kwargs...)
  @chain begin
    rp
    send_request
    extract_content
    deserialize_content(_, parser; kwargs...)
  end
end

"""
  resolve_request_parameters(tickers, resource; kwargs...)

Replaces all template values from `resource.headers` and `resouurce.params` with the concrete values passed as keyword arguments.
Puts all the resolved fields  into a RequestParams struct.
"""
function resolve_request_parameters(
  tickers::Vector{UpdatableSymbol}, resource::APIResource; kwargs...
)::Union{RequestParams,Exception}
  symbol_value = join(map(x -> x.ticker, tickers), ",")
  args = Dict(kwargs)
  from, to = get_minimum_dates(
    tickers; from=get(args, :from, missing), to=get(args, :to, missing)
  )
  kwargs_to_add = Dict()
  kwargs_to_add[:from] = from
  kwargs_to_add[:to] = to
  if haskey(resource.query_params, resource.symbol_key) || contains(resource.url, resource.symbol_key)
    kwargs_to_add[Symbol(resource.symbol_key)] = symbol_value
  end
  if length(tickers) == 1
    ticker = first(tickers)
    if !ismissing(ticker.fx_pair) && !haskey(resource.query_params, resource.symbol_key)
      base, target = ticker.fx_pair
      kwargs_to_add[:base] = base
      kwargs_to_add[:target] = target
    end
  end
  kwargs_new = @chain merge(args, kwargs_to_add) filter(kw -> kw[1] !== :symbol_key, _)
  url_part = resolve_url_params(resource.url; symbol_key=symbol_value, kwargs_new...)
  query_params = resolve_query_params(resource.query_params, resource.url; kwargs_new...)
  validation = validate_parameter_substitution(url_part, query_params)
  isa(validation, RequestBuilderError) && return validation
  return RequestParams(;
    url=url_part,
    params=query_params,
    retries=resource.max_retries,
    tickers=tickers,
    headers=resource.headers,
    query=join(["$k=$v" for (k, v) in query_params], "&"),
    from=from,
    to=to,
  )
end

function validate_parameter_substitution(url::String, query_params::Dict{String,String})
  pattern = r"(?<=\{).+?(?=\})"
  url_params_unsolved = match(pattern, url)
  if url_params_unsolved !== nothing
    return RequestBuilderError(
      "Couldn't substitute url parameter at $(url_params_unsolved.match)"
    )
  end
  query_params_unsolved = @chain begin
    [match(pattern, v) for (k, v) in query_params]
    filter(!isnothing, _)
    !isempty(_) ? map(m -> m.match, _) : []
  end
  if !isempty(query_params_unsolved)
    qp = join(query_params_unsolved, ", ")
    return RequestBuilderError("Couldn't substitute query parameters: '$qp'")
  end
end

function send_request(req::RequestParams)::Union{HTTP.Response,Exception}
  try
    resp = HTTP.request(
      "GET", req.url; headers=req.headers, query=req.query, retries=req.retries
    )
    if resp.status != 200
      return APIResponseError("Code: $(resp.status)")
    end
    return resp
  catch err
    if isa(err, HTTP.ExceptionRequest.StatusError)
      msg = String(err.response.body)
      return APIResponseError("Code: $(err.status); msg: $msg")
    end
    return APIResponseError(String(err))
  end
end

function extract_content(
  resp::Union{HTTP.Response,Exception}
)::Union{String,Exception}
  !isa(resp, HTTP.Response) && return resp
  try
    content = String(resp.body)
    # @TODO: Add checks on content length
    return content
  catch err
    return err
  end
end

function deserialize_content(
  c::Union{String,Exception}, parser::AbstractContentParser; kwargs...
)
  !isa(c, String) && return c
  content = c
  try
    return parse_content(parser, content; kwargs...)
  catch err
    return APIResponseError(err.msg)
  end
end

function resolve_query_params(
  params::Dict{String,String}, url::AbstractString=""; kwargs...
)::Dict{String,String}
  pattern = r"(?<=\{).+?(?=\})"
  keys_to_resolve = [k for (k, v) in params if match(pattern, v) !== nothing]
  isempty(keys_to_resolve) && return params
  params_new = Dict()
  # insert values that don't need to be resolved
  [params_new[k] = v for (k, v) in params if match(pattern, v) === nothing]
  args = Dict(kwargs)
  for key in keys_to_resolve
    mutated_value = params[key]
    for m in eachmatch(pattern, params[key])
      kwarg_key = replace(m.match, r"[{}]" => "") # {param} => param
      kwarg_value = get(args, Symbol(kwarg_key), missing)
      if !ismissing(kwarg_value)
        new_value = replace_params_with_kwargs(
          kwarg_value, pattern; Dict(Symbol(kwarg_key) => kwarg_value)...
        )
        mutated_value = replace(mutated_value, "{$kwarg_key}" => new_value)
      end
    end
    params_new[key] = mutated_value
  end
  args = Dict(kwargs)
  specific_params = resolve_specific_query_params(
    url; from=get(args, :from, missing), to=get(args, :to, missing)
  )
  return merge(params_new, specific_params)
end

function resolve_url_params(url::AbstractString; kwargs...)
  return replace_params_with_kwargs(url, r"(?<=\{).+?(?=\})"; kwargs...)
end

function resolve_specific_query_params(url::AbstractString; from=missing, to=missing)
  specific_params = Dict()
  default_date = Date(unix2datetime(0)) # or Dates.today() - Dates.Day(28) ?
  if contains(url, "alphavantage")
    dt_from = isa(from, Date) ? from : default_date
    specific_params["outputsize"] = dt_from > today() - Day(99) ? "compact" : "full"
  end
  if contains(url, "yfapi")
    max_date = isa(from, Date) ? from : default_date
    specific_params["range"] = convert_date_to_range_expr(max_date)
  end
  return specific_params
end

function replace_params_with_kwargs(s::AbstractString, pattern::Regex; kwargs...)
  new_s = s
  for (k, v) in kwargs
    for m in eachmatch(pattern, s)
      if m.match == String(k)
        new_s = replace(new_s, "{$(m.match)}" => v)
      end
    end
  end
  return new_s
end

function convert_date_to_range_expr(end_date::Date)
  days_delta = days(today() - end_date)
  in(days_delta, 0:1) && return "1d"
  in(days_delta, 2:5) && return "5d"
  in(days_delta, 6:28) && return "1mo"
  in(days_delta, 29:88) && return "3mo"
  in(days_delta, 89:179) && return "6mo"
  in(days_delta, 180:364) && return "1y"
  in(days_delta, 365:1820) && return "5y"
  return "5y"
end

"""
What to do when there are multiple requests, some succesfull, but some failed.
  - optimistic approach => if at least one request succedes, return a Vector{AbstractStonksRecord},
    but print warnings for all requests which failed
"""
function optimistic_request_resolution(
  ::Type{T}, c::Channel
)::Union{Vector{T},Exception} where {T<:AbstractStonksRecord}
  result = T[]
  failures = Exception[]
  for response in c
    if response.err !== nothing
      @warn "Request to: $(response.params.url) failed. $(response.err)"
      push!(failures, response.err)
    else
      append!(result, response.value)
    end
  end
  if !isempty(result)
    return result
  end
  if !isempty(failures)
    return first(failures)
  end
  return result
end

end


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie
CairoMakie.activate!() # HIDE
let
    x = LinRange(0,2π,50)
    fig = Figure(resolution = (700, 450), font =:sans)
    ax = Axis(fig, xlabel = "x", ylabel = "")
    line1 = lines!(x, sin.(x), color = :black)
    pnts1 = scatter!(x, sin.(x), color = :black, marker = '◆', markersize = 10)
    line2 = lines!(x, cos.(x), color = :black)
    pnts2 = scatter!(x, cos.(x), color = :red, strokecolor = :red, marker = '■',
                    markersize = 10)
    leg = Legend(fig, [[line1, pnts1], [line2, pnts2]], ["sin(x)", "cos(x)"],
        halign = :center, valign = :center, markersize = 7,tellheight = false,
        tellwidth = false, bgcolor = :white, framecolor=:white)
    fig[1, 1] = ax
    fig[1, 1] = leg
    fig
    #save("FigLineScatter.png", fig, px_per_unit = 2)
    save(joinpath(@__DIR__, "output", "FigLinesScatters.png"), fig, px_per_unit = 2) # HIDE
end
using Pkg # HIDE
Pkg.status("CairoMakie") # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using Test

@testset "FrameFunDomains" begin
    include("test_framefundomains.jl")
end
@testset "Platforms" begin
    include("test_platforms.jl")
end
@testset "ExtensionFrames" begin
    include("test_extensionframes.jl")
end
@testset "FrameFunInterface" begin
    include("test_framefuninterface.jl")
end
@testset "ExtensionFramePlatforms" begin
    include("test_extensionframeplatforms.jl")
end
@testset "Other Frame Platforms" begin
    include("test_otherframeplatforms.jl")
end
@testset "FrameFunPlatforms" begin
    include("test_framefunplatforms.jl")
end
@testset "FrameFunInterfaceExtension" begin
    include("test_framefuninterfaceextension.jl")
end
@testset "Adaptivity" begin
    include("test_adaptivity.jl")
end
include("test_scenariolist.jl")
include("test_suite.jl")
include("test_suite_applications.jl")
include("test_specificplatforms.jl")
include("test_plots.jl")

# include("test_notebooks.jl")
# include("create_readme.jl")

println("All tests succeeded!!!")


text:
using MyExample
using Documenter

makedocs(;
    modules=[MyExample],
    authors="Spencer Gibbs",
    repo="https://github.com/sgibbs409/MyExample.jl/blob/{commit}{path}#L{line}",
    sitename="MyExample.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://sgibbs409.github.io/MyExample.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/sgibbs409/MyExample.jl",
)


text:
using Documenter, DoubleFloats

makedocs(
    modules = [DoubleFloats],
    format = :html,
    sitename = "DoubleFloats.jl",
    authors = "Jeffrey Sarnoff, Simon Byrne, Sascha Timme, and other contributors",
    pages = [
        "Home" => "index.md",
        "Overview" => "appropriate.md",
        "Names" => "naming.md",
        "Capabilities" => "capabilities.md",
        "Refs" => "references.md"
    ],
)

deploydocs(
    julia = "1.0",
    repo = "github.com/JuliaMath/DoubleFloats.jl.git",
    target = "build",
    osname = "linux",
    deps   = nothing,
    make   = nothing
)


text:
using Test
using BSDiff
using PkgFetcher
using Pkg.Artifacts

@testset "BSPatch" begin
    registry_data = joinpath(artifact"test_data", "registry")
    old = joinpath(registry_data, "before.tar")
    new = joinpath(registry_data, "after.tar")
    patch = bsdiff(old, new)
    new′ = tempname()
    PkgFetcher.bspatch(old, new′, patch)
    @test read(new) == read(new′)
end

@testset "PkgFetcher.jl" begin
    # Write your tests here.
end