text:
#
# quasinewton.jl --
#
# Limited memory quasi-Newton methods for OptimPackNextGen.
#
# -----------------------------------------------------------------------------
#
# This file is part of OptimPackNextGen.jl which is licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2015-2021, Éric Thiébaut.
# <https://github.com/emmt/OptimPackNextGen.jl>.
#

# Improvements:
# - skip updates such that rho <= 0;
# - check for sufficient descent condition and implement restarts;
# - add other convergence criteria;
# - in case of early stop, revert to best solution so far;
# - better initial step than 1/vnorm2(g);

module QuasiNewton

export
    vmlmb,
    vmlmb!,
    EMULATE_BLMVM

using Printf, Requires

using LazyAlgebra
using ...OptimPackNextGen
using OptimPackNextGen.LineSearches
using OptimPackNextGen.SimpleBounds

# Use the same floating point type for scalars as in OptimPackNextGen.
import OptimPackNextGen.Float

const EMULATE_BLMVM = UInt(1)

# All scalar computations are done in double precision.  Thus manage to have
# `vdot` and `vnorm2` return a double precision result.

function vdot(x::AbstractArray{Float,N},
              y::AbstractArray{Float,N}) where {N}
    return LazyAlgebra.vdot(x, y)
end

function vdot(x::AbstractArray{<:Real,N},
              y::AbstractArray{<:Real,N}) where {N}
    return Float(LazyAlgebra.vdot(x, y))
end

function vdot(w::Union{AbstractArray{Float,N},AbstractVector{Int}},
              x::AbstractArray{Float,N},
              y::AbstractArray{Float,N}) where {N}
    return LazyAlgebra.vdot(w, x, y)
end

function vdot(w::Union{AbstractArray{<:Real,N},AbstractVector{Int}},
              x::AbstractArray{<:Real,N},
              y::AbstractArray{<:Real,N}) where {N}
    return Float(LazyAlgebra.vdot(w, x, y))
end

# FIXME: clean this!
vnorm2(x::AbstractArray{Float,N}) where {N} = LazyAlgebra.vnorm2(x)
vnorm2(x::AbstractArray{<:Real,N}) where {N} = Float(LazyAlgebra.vnorm2(x))

"""
## VMLMB: limited memory BFGS method with optional bounds

    x = vmlmb(fg!, x0; mem=..., lower=..., upper=..., ftol=..., fmin=...)

computes a local minimizer of a function of several variables by a limited
memory variable metric method.  The caller provides a function `fg!` to compute
the value and the gradient of the function as follows:

    fx = fg!(x, gx)

where `x` are the current variables, `fx` is the value of the function at `x`
and the contents of `gx` has to be overwritten with the gradient at `x` (when
`fg!` is called, `gx` is already allocated as `gx = vcreate(x0)`).  Argument `x0`
gives the initial approximation of the variables (its contents is left
unchanged).  The best solution found so far is returned in `x`.

If `Zygote` package has been loaded, the gradient of the function may be
automatically computed and the first argument to `vmlmb` can be a simpler
function, say `f`, that takes the variables `x` as a single argument and returns
the function value:

    fx = f(x)

It is also possible to extend the method
[`OptimPackNextGen.QuasiNewton.auto_differentiate!](@ref) to compute the gradient
for a specific function `f`.

The following keywords are available:

* `mem` specifies the amount of storage.

* `autodiff` is a boolean specifying whether to rely on
  automatic-differentiation by `Zygote` to compute the gradient of the
  objective function.  If not specified, the decision is based on whether a
  call like `fg!(x,g)` to compute the objective function and its gradient is
  applicable.

* `xtol` is a tuple of two nonnegative reals specifying respectively the
  absolute and relative tolerances for deciding convergence on the variables.
  Convergence occurs if the Euclidean norm of the the difference between
  successive iterates is less or equal `max(xtol[1], xtol[2]*vnorm2(x))`.  By
  default, `xtol = (0.0,1e-7)`.

* `ftol` is a tuple of two nonnegative reals specifying respectively the
  absolute and relative errors desired in the function.  Convergence occurs if
  the absolute error between `f(x)` and `f(xsol)` is less than `ftol[1]` or if
  the estimate of the relative error between `f(x)` and `f(xsol)`, where `xsol`
  is a local minimizer, is less than `ftol[2]`.  By default, `ftol =
  (0.0,1e-8)`.

* `gtol` is a tuple of two nonnegative reals specifying the absolute and a
  relative thresholds for the norm of the gradient, convergence is assumed as
  soon as:

      ||g(x)|| <= hypot(gtol[1], gtol[2]*||g(x0)||)

  where `||g(x)||` is the Euclidean norm of the gradient at the current
  solution `x`, `||g(x0)||` is the Euclidean norm of the gradient at the
  starting point `x0`.  By default, `gtol = (0.0,1e-6)`.

* `fmin` specifies a lower bound for the function.  If provided, `fmin` is used
  to estimate the steepest desxecnt step length this value.  The algorithm
  exits with a warning if `f(x) < fmin`.

* `maxiter` specifies the maximum number of iterations.

* `maxeval` specifies the maximum number of calls to `fg!`.

* `verb` specifies the verbosity level. It can be a boolean to specify
  whether to print information at every iteration or an integer to print
  information every `verb` iteration(s).  No information is printed if
  `verb` is less or equal zero. The default is `verb = false`.

* `printer` can be set with a user defined function to print iteration
  information, its signature is:

      printer(io::IO, iter::Integer, eval::Integer, rejects::Integer,
              f::Real, gnorm::Real, stp::Real)

  where `io` is the output stream, `iter` the iteration number (`iter = 0` for
  the starting point), `eval` is the number of calls to `fg!`, `rejects` is the
  number of times the computed direction was rejected, `f` and `gnorm` are the
  value of the function and norm of the gradient at the current point, `stp` is
  the length of the step to the current point.

* `output` specifies the output stream for printing information (`stdout` is
  used by default).

* `lnsrch` specifies the method to use for line searches (the default
   line search is `MoreThuenteLineSearch`).

* `lower` and `upper` specify the lower and upper bounds for the variables.
   The bound can be a scalar to indicate that all variables have the same bound
   value.  If the lower (resp. upper) bound is unspecified or set to `±∞`, the
   variables are assumed to be unbounded below (resp. above).  If no bounds are
   set, VMLMB amounts to an unconstrained limited memory BFGS method (L-BFGS).

* `blmvm` can be set true to emulate the BLMVM algorithm of Benson and Moré.
  This option has no effects for an uncostrained problem.

### History

The VMLMB algorithm in
[OptimPackNextGen](https://github.com/emmt/OptimPackNextGen.jl) provides a pure
Julia implementation of the original method (Thiébaut, 2002) with some
improvements and the capability to emulate L-BFGS and BLMVM methods.

The limited memory BFGS method (L-BFGS) was first described by Nocedal (1980)
who dubbed it SQN.  The method is implemented in MINPACK-2 (1995) by the
FORTRAN routine VMLM.  The numerical performances of L-BFGS have been studied
by Liu and Nocedal (1989) who proved that it is globally convergent for
unfiformly convex problems with a R-linear rate of convergence.  They provided
the FORTRAN code LBFGS.  The BLMVM and VMLMB algorithms were proposed by Benson
and Moré (2001) and Thiébaut (2002) to account for separable bound constraints
on the variables.  These two latter methods are rather different than L-BFGS-B
by Byrd at al. (1995) which has more overheads and is slower.

* J. Nocedal, "*Updating Quasi-Newton Matrices with Limited Storage*" in
  Mathematics of Computation, vol. 35, pp. 773-782 (1980).

* D.C. Liu & J. Nocedal, "*On the limited memory BFGS method for large scale
  optimization*" in Mathematical programming, vol. 45, pp. 503-528
  (1989).

* R.H. Byrd, P. Lu, J. Nocedal, & C. Zhu, "*A limited memory algorithm for
  bound constrained optimization*" in SIAM Journal on Scientific Computing,
  vol. 16, pp. 1190-1208 (1995).

* S.J. Benson & J.J. Moré, "*A limited memory variable metric method in
  subspaces and bound constrained optimization problems*" in Subspaces and
  Bound Constrained Optimization Problems (2001).

* É. Thiébaut, "*Optimization issues in blind deconvolution algorithms*" in
  Astronomical Data Analysis II, Proc. SPIE 4847, pp. 174-183 (2002).

"""
vmlmb(fg!::Function, x0; kwds...) = vmlmb!(fg!, vcopy(x0); kwds...)

"""
`vmlmb!` is the in-place version of `vmlmb` (which to see):

     vmlmb!(fg!, x; mem=..., lower=..., upper=..., ftol=..., fmin=...) -> x

finds a local minimizer of `f(x)` starting at `x` and stores the best solution
in `x`.
"""
function vmlmb!(fg!::Function, x::T;
                mem::Integer = min(5, length(x)),
                lower::Union{Real,T} = -Inf,
                upper::Union{Real,T} = +Inf,
                autodiff::Bool = !applicable(fg!, x, x),
                blmvm::Bool = false,
                fmin::Real = -Inf,
                maxiter::Integer = typemax(Int),
                maxeval::Integer = typemax(Int),
                xtol::NTuple{2,Real} = (0.0, 1e-7),
                ftol::NTuple{2,Real} = (0.0, 1e-8),
                gtol::NTuple{2,Real} = (0.0, 1e-6),
                epsilon::Real = 0.0,
                verb::Integer = false,
                printer::Function = print_iteration,
                output::IO = stdout,
                lnsrch::Union{LineSearch{Float},Nothing} = nothing) where {T}
    # Determine which options are used.
    flags::UInt = (blmvm ? EMULATE_BLMVM : 0)

    # Determine the type of bounds.
    bounds::UInt = 0
    if isa(lower, Real)
        lo = Float(lower)
        if lower > -Inf
            bounds |= 1
        end
    elseif isa(lower, T)
        lo = lower
        bounds |= 1
    else
        error("invalid lower bound type")
    end
    if isa(upper, Real)
        hi = Float(upper)
        if upper < +Inf
            bounds |= 2
        end
    elseif isa(upper, T)
        hi = upper
        bounds |= 2
    else
        error("invalid upper bound type")
    end

    # Determine the optimization method (0 for L-BFGS, 1 for BLMVM or 2 for
    # VMLMB).
    method = (bounds == 0 ? 0 : (flags & EMULATE_BLMVM) != 0 ? 1 : 2)

    # Provide a default line search method if needed.
    if isa(lnsrch, Nothing)
        if method == 0
            ls = MoreThuenteLineSearch(Float; ftol=1e-3, gtol=0.9, xtol=0.1)
        else
            ls = MoreToraldoLineSearch(Float; ftol=1e-3, gamma=(0.1,0.5))
        end
    else
        ls = lnsrch
    end

    # Call the real method.
    _vmlmb!(fg!, x, Int(mem), flags, lo, hi, bounds, autodiff, method,
            Float(fmin), Int(maxiter), Int(maxeval),
            Float(xtol[1]), Float(xtol[2]),
            Float(ftol[1]), Float(ftol[2]),
            Float(gtol[1]), Float(gtol[2]),
            Float(epsilon), Int(verb), printer, output, ls)
end

# The real worker.
function _vmlmb!(fg!::Function, x::T, mem::Int, flags::UInt,
                 lo::Union{Float, T},
                 hi::Union{Float, T},
                 bounds::UInt,
                 autodiff::Bool,
                 method::Int,
                 fmin::Float, maxiter::Int, maxeval::Int,
                 xatol::Float, xrtol::Float,
                 fatol::Float, frtol::Float,
                 gatol::Float, grtol::Float,
                 epsilon::Float,
                 verb::Int, printer::Function, output::IO,
                 lnsrch::LineSearch{Float}) where {T}

    @assert mem ≥ 1
    @assert maxiter ≥ 0
    @assert maxeval ≥ 1
    @assert xatol ≥ 0
    @assert xrtol ≥ 0
    @assert fatol ≥ 0
    @assert frtol ≥ 0
    @assert gatol ≥ 0
    @assert grtol ≥ 0
    @assert 0 ≤ epsilon < 1

    STPMIN = Float(1e-20)
    STPMAX = Float(1e+20)

    sel = Array{Int}(undef, 0)
    if bounds != 0
        sizehint!(sel, length(x))
    end
    mem = min(mem, length(x))
    reason::AbstractString = ""
    fminset::Bool = (! isnan(fmin) && fmin > -Inf)

    mark::Int = 1    # index of most recent step and gradient difference
    m::Int = 0       # number of memorized steps
    iter::Int = 0    # number of algorithm iterations
    eval::Int = 0    # number of objective function and gradient evaluations
    rejects::Int = 0 # number of rejected search directions

    f::Float = 0
    f0::Float = 0
    gd::Float = 0
    gd0::Float = 0
    stp::Float = 0
    stpmin::Float = 0
    stpmax::Float = 0
    smin::Float = 0
    smax::Float = 0
    gamma::Float = 1
    gnorm::Float = 0
    gtest::Float = 0

    # Variables for saving information about best point so far.
    beststp::Float = 0
    bestf::Float = 0
    bestgnorm::Float = 0

    # Allocate memory for the limited memory BFGS approximation of the inverse
    # Hessian.
    g = vcreate(x) # ------------> gradient
    if method > 0
        p = vcreate(x) # --------> projected gradient
    end
    d = vcreate(x) # ------------> search direction
    S = Array{T}(undef, mem) # --> memorized steps
    Y = Array{T}(undef, mem) # --> memorized gradient differences
    s = vcreate(x) # ------------> for effective step
    for k in 1:mem
        S[k] = vcreate(x)
        Y[k] = vcreate(x)
    end
    rho = Array{Float}(undef, mem)
    alpha = Array{Float}(undef, mem)

    # Variable used to control the stage of the algorithm:
    #   * stage = 0 at start or restart;
    #   * stage = 1 during a line search;
    #   * stage = 2 if line search has converged;
    #   * stage = 3 if algorithm has converged;
    #   * stage = 4 if algorithm is terminated with a warning.
    stage::Int = 0

    while true

        # Compute value of function and gradient, register best solution so far
        # and check for convergence based on the gradient norm.
        if method > 0
            project_variables!(x, x, lo, hi)
        end
        # if fg! takes a single argument (x), its derivative are automatically computed using Zygotes.jl
        if autodiff
            f = auto_differentiate!(fg!, x, g)
        else
            f = fg!(x, g)
        end

        eval += 1
        if method > 0
            project_direction!(p, x, lo, hi, -1, g)
        end
        if eval == 1 || f < bestf
            gnorm = vnorm2((method > 0 ? p : g))
            beststp = stp
            bestf = f
            bestgnorm = gnorm
            if eval == 1
                gtest = max(gatol, grtol*gnorm)
            end
            if gnorm ≤ gtest
                stage = 3
                if gnorm == 0
                    reason = "a stationary point has been found!"
                elseif method > 0
                    reason = "projected gradient sufficiently small"
                else
                    reason = "gradient sufficiently small"
                end
                if eval > 1
                    iter += 1 # FIXME: do this here???
                end
            end
        end
        if fminset && f < fmin
            stage = 4
            reason = "f < fmin"
        end

        if stage == 1
            # Compute effective step.
            vcombine!(s, 1, x, -1, S[mark])

            # Line search is in progress.
            if usederivatives(lnsrch)
                if method > 0
                    gd = vdot(g, s)/stp
                else
                    gd = -vdot(g, d)
                end
            end
            task = iterate!(lnsrch, stp, f, gd)
            if task == :SEARCH
                stp = getstep(lnsrch)
            elseif task == :CONVERGENCE
                # Line search has converged.  Increment iteration counter
                # and check for stopping conditions.
                iter += 1
                xtest = max(xatol, zero(Float))
                if xrtol > 0
                    xtest = max(xtest, xrtol*vnorm2(x))
                end
                if vnorm2(s) ≤ xtest
                    stage = 3
                    reason = "X test satisfied"
                else
                    delta = max(abs(f - f0), stp*abs(gd0))
                    if delta ≤ fatol
                        stage = 3
                        reason = "fatol test satisfied"
                    elseif delta ≤ frtol*abs(f0)
                        stage = 3
                        reason = "frtol test satisfied"
                    elseif iter ≥ maxiter
                        stage = 4
                        reason = "too many iterations"
                    elseif eval ≥ maxeval
                        stage = 4
                        reason = "too many evaluations"
                    else
                        stage = 2
                    end
                end
            else
                # Line seach terminated with a warningor an error.
                stage = 4
                reason = getreason(lnsrch)
            end
        end
        if stage < 2 && eval ≥ maxeval
            stage = 4
            reason = "too many evaluations"
        end

        if stage != 1
            # Initial step or line search has converged.

            # Make sure the gradient norm is correct or revert to best solution
            # so far if something wrong occured.
            if stp != beststp
                if stage ≥ 3
                    stp = beststp
                    f = bestf
                    gnorm = bestgnorm
                    vcombine!(x, 1, S[mark], -stp, d)
                    if method > 0
                        project_variables!(x, x, lo, hi)
                    end
                else
                    gnorm = vnorm2((method > 0 ? p : g))
                end
            end

            # Print some information if requested and terminate algorithm if
            # stage ≥ 3.
            if verbose(verb, iter)
                printer(output, iter, eval, rejects, f, gnorm, stp)
            end
            if stage ≥ 3
                break
            end

            # Compute next search direction.
            if stage == 2
                # Update limited memory BFGS approximation of the inverse
                # Hessian with the effective step and gradient change.
                vupdate!(S[mark], -1, x) # FIXME: vcopy!(S[mark], s)?
                vupdate!(Y[mark], -1, (method == 1 ? p : g))
                if method < 2
                    rho[mark] = vdot(Y[mark], S[mark])
                    if rho[mark] > 0
                        # The update is acceptable, compute the scale.
                        gamma = rho[mark]/vdot(Y[mark], Y[mark])
                    end
                end
                m = min(m + 1, mem)

                # Compute search direction.
                if method < 2
                    vcopy!(d, g)
                    change = apply_lbfgs!(S, Y, rho, gamma, m, mark, d, alpha)
                else
                    vcopy!(d, p)
                    change = apply_lbfgs!(S, Y, rho, m, mark, d, alpha,
                                          get_free_variables!(sel, d))
                end
                # FIXME: speedup project_direction with the free vars.?
                if change
                    if method > 0
                        project_direction!(d, x, lo, hi, -1, d)
                    end
                    gd = -vdot(g, d)
                    reject = ! sufficient_descent(gd, epsilon, gnorm, d)
                else
                    reject = true
                end
                if reject
                    # The steepest descent will be used.
                    stage = 0
                    rejects += 1
                end

                # Circularly move the mark to the next slot.
                mark = (mark < mem ? mark + 1 : 1)
            end
            if stage == 0
                # At the initial point or computed direction was rejected.  Use
                # the steepest descent which is the projected gradient.
                vcopy!(d, (method > 0 ? p : g))
                gd = -gnorm^2
            end

            # Save function value, variables and (projected) gradient at start
            # of line search.
            f0 = f
            gd0 = gd
            vcopy!(S[mark], x)
            vcopy!(Y[mark], (method == 1 ? p : g))

            # Choose initial step length.
            if stage == 2
                stp = 1
            else
                if fminset && fmin < f0
                    stp = 2*(fmin - f0)/gd0
                else
                    stp = 1/gnorm # FIXME: use a better scale
                end
            end
            if method > 0
                # Make sure the step is not longer than necessary.
                (smin, smax) = step_limits(x, lo, hi, -1, d)
                stp = min(stp, smax)
                stpmin = STPMIN*stp
                stpmax = min(STPMAX*stp, smax)
            else
                stpmin = STPMIN*stp
                stpmax = STPMAX*stp
            end

            # Initialize the line search.
            task = start!(lnsrch, f0, gd0, stp; stpmin=stpmin, stpmax=stpmax)
            if task != :SEARCH
                # Something wrong happens.
                stage = 4
                reason = getreason(lnsrch)
                break
            end
            stage = 1 # line search is in progress

        end

        # Compute the new iterate.
        vcombine!(x, 1, S[mark], -stp, d)

    end

    # Algorithm finished.
    if verbose(verb, 0) #always print last line if verb>0
        color = (stage > 3 ? :red : :green)
        prefix = (stage > 3 ? "WARNING: " : "CONVERGENCE: ")
        printstyled(output, "# ", prefix, reason, "\n"; color=color)
        #elseif stage > 3
        #@warn(reason)
    end
    return x
end

"""
    sufficient_descent(gd, ε, gnorm, d)

checks whether `d` is a sufficient descent direction (Zoutenjdik condition).
Argument `gd` is the scalar product `vdot(g,d)` between the gradient `g` and
the direction `d`, `ε` is a nonnegative small value in the range `[0,1)` and
`gnorm` is the Euclidean norm of the gradient `g`.

If the Euclidean norm of `d` has already been computed, then

    sufficient_descent(gd, ε, gnorm, dnorm)

should be used instead with `dnorm = vnorm2(d)`.
"""
sufficient_descent(gd::Real, ε::Real, gnorm::Real, d) =
    sufficient_descent(gd, ε, gnorm, vnorm2(d))

sufficient_descent(gd::Real, ε::Real, gnorm::Real, dnorm::Real) =
    ε > 0 ? (gd ≤ -ε*dnorm*gnorm) : gd < 0

"""
    verbose(verb, iter)

yields whether to print information at iteration `iter` with verbose level `verb`.

"""
verbose(verb::Integer,iter::Integer) =
    (verb > 0 && iter%verb == 0)

function print_iteration(iter::Integer, eval::Integer, rejects::Integer,
                         f::Real, gnorm::Real, step::Real)
    print_iteration(stdout, iter, eval, rejects, f, gnorm, step)
end

function print_iteration(io::IO, iter::Integer, eval::Integer, rejects::Integer,
                         f::Real, gnorm::Real, step::Real)
    if iter == 0
        @printf(io, "#%s%s\n#%s%s\n",
                " ITER   EVAL   REJECTS",
                "          F(X)           ||G(X)||    STEP",
                "----------------------",
                "-------------------------------------------")
    end
    @printf(io, " %5d  %5d  %5d  %24.16E %9.2E %9.2E\n",
            iter, eval, rejects, f, gnorm, step)
end

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

"""
One of the calls:

    modif = apply_lbfgs!(S, Y, rho, gamma, m, mark, v, alpha)
    modif = apply_lbfgs!(S, Y, rho, d,     m, mark, v, alpha)
    modif = apply_lbfgs!(S, Y, rho, H0!,   m, mark, v, alpha)

computes the matrix-vector product `H*v` where `H` is the limited memory
inverse BFGS approximation.  Operation is done *in-place*: on entry, argument
`v` contains the input vector `v`; on exit, argument `v` contains the
matrix-vector product `H*v`.  The returned value, `modif`, is a boolean
indicating whether the initial vector was modified (BFGS updates are skipped if
`rho[k] ≤ 0`).

The linear operator `H` depends on an initial approximation of the inverse
Hessian (see below), `m` steps `S[...]` and `m` gradient differences `Y[...]`.
The most recent step and gradient difference are stored in `S[mark]` and
`Y[mark]`, respectively.

Argument `rho` contains the inner products of the steps and the gradient
differences: `rho[k] = vdot(S[k],Y[k])`.  On exit rho is unchanged.

Argument `alpha` is a work vector of length at least `m`.

The initial approximation of the inverse Hessian can be specified in several
ways:

* by a scalar `gamma > 0` to assume that the initial approximation of the
  inverse Hessia is `gamma*I` where `I` is the identity;

* by a some instance of the same type as `v` (whose components are all strictly
  positive) to assume that the initial approximation of the inverse Hessian is
  a diagonal matrix whose diagonal terms are given by `d`;

* by a function `H0!` to assume that `H0!(u)` (for any instance `u` of the same
  type as `v`) applies the initial approximation of the inverse Hessian to `u`
  and stores the result in `u` (in-place operation).

"""
function apply_lbfgs!(S::Vector{T}, Y::Vector{T}, rho::Vector{Float},
                      gamma::Float, m::Int, mark::Int, v::T,
                      alpha::Vector{Float}) where {T}
    @assert gamma > 0
    apply_lbfgs!(S, Y, rho, u -> vscale!(u, gamma), m, mark, v, alpha)
end

function apply_lbfgs!(S::Vector{T}, Y::Vector{T}, rho::Vector{Float},
                      d::T, m::Int, mark::Int, v::T,
                      alpha::Vector{Float}) where {T}
    apply_lbfgs!(S, Y, rho, u -> vproduct!(u, d, u), m, mark, v, alpha)
end

function apply_lbfgs!(S::Vector{T}, Y::Vector{T}, rho::Vector{Float},
                      H0!::Function, m::Int, mark::Int, v::T,
                      alpha::Vector{Float}) where {T}
    mem = min(length(S), length(Y), length(rho), length(alpha))
    @assert 1 ≤ m ≤ mem
    @assert 1 ≤ mark ≤ mem
    modif::Bool = false
    @inbounds begin
        k::Int = mark + 1
        for i in 1:m
            k = (k > 1 ? k - 1 : mem)
            if rho[k] > 0
                alpha[k] = vdot(S[k], v)/rho[k]
                vupdate!(v, -alpha[k], Y[k])
                modif = true
            end
        end
        if modif
            H0!(v)
            for i in 1:m
                if rho[k] > 0
                    beta::Float = vdot(Y[k], v)/rho[k]
                    vupdate!(v, alpha[k] - beta, S[k])
                end
                k = (k < mem ? k + 1 : 1)
            end
        end
    end
    return modif
end

function apply_lbfgs!(S::Vector{T}, Y::Vector{T}, rho::Vector{Float},
                      m::Int, mark::Int, v::T,
                      alpha::Vector{Float}, sel) where {T}
    mem = min(length(S), length(Y), length(rho), length(alpha))
    @assert 1 ≤ m ≤ mem
    @assert 1 ≤ mark ≤ mem
    gamma::Float = 0
    @inbounds begin
        k::Int = mark + 1
        for i in 1:m
            k = (k > 1 ? k - 1 : mem)
            rho[k] = vdot(sel, Y[k], S[k])
            if rho[k] > 0
                alpha[k] = vdot(sel, S[k], v)/rho[k]
                vupdate!(v, sel, -alpha[k], Y[k])
                if gamma == 0
                    gamma = rho[k]/vdot(sel, Y[k], Y[k])
                end
            end
        end
        if gamma != 0
            vscale!(v, gamma)
            for i in 1:m
                if rho[k] > 0
                    beta::Float = vdot(sel, Y[k], v)/rho[k]
                    vupdate!(v, sel, alpha[k] - beta, S[k])
                end
                k = (k < mem ? k + 1 : 1)
            end
        end
    end
    return (gamma != 0)
end

"""
    OmptimPackNextGen.QuasiNewton.auto_differentiate!(f, x, g) -> fx

yields `fx = f(x)` for a given function `f` and variables `x` and overwrites the
contents of `g` with the gradient of the function at `x`.

This method may be extended to compute the gradient and function value for specific
`typeof(f)` or to automatically compute the gradient as can be done by the `Zygote`
package if it is loaded.

"""
auto_differentiate!(arg...; kwds...) = error("Zygote package must be loaded first")

function __init__()
    @require Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" auto_differentiate!(f, x, g) = begin
        vcopy!(g, Zygote.gradient(f, x)[1]);
        return f(x)
    end
end

end # module


text:
using Neuroimaging
using Test
using Logging
using MAT
using Plots
using Glob
using DataDeps
using LinearAlgebra
using Unitful

logger = SimpleLogger(stdout, Logging.Warn)

#
# Run all tests
#

tests = glob("**/*.jl")
append!(tests, glob("**/**/*.jl"))
append!(tests, glob("**/**/**/*.jl"))
tests = tests[.~(tests .== "runtests.jl")]

with_logger(logger) do
    for t in tests
        include(t)
    end
end

ext: jl


text:
@diffusion_process Lorenz begin
    :dimensions
    process --> 3
    wiener --> 3

    :parameters
    _ --> (3, Float64)
    σ --> Float64

    :additional
    constdiff --> true
end

function b(t, x, P::Lorenz)
    @SVector [
        P.p1*(x[2]-x[1]),
        P.p2*x[1] - x[2] - x[1]*x[3],
        x[1]*x[2] - P.p3*x[3]
    ]
end

σ(t, x, P::Lorenz) = SDiagonal(P.σ, P.σ, P.σ)


text:
include("../src/components/PermafrostSiBCASA.jl")
include("../src/components/PermafrostJULES.jl")
include("../src/components/PermafrostTotal.jl")

using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = page_model()

    permafrost_sibcasa = add_comp!(m, PermafrostSiBCASA)
    permafrost_jules = add_comp!(m, PermafrostJULES)
    permafrost_total = add_comp!(m, PermafrostTotal)

    permafrost_total[:perm_sib_ce_co2] = permafrost_sibcasa[:perm_sib_ce_co2]
    permafrost_total[:perm_sib_e_co2] = permafrost_sibcasa[:perm_sib_e_co2]
    permafrost_total[:perm_sib_ce_ch4] = permafrost_sibcasa[:perm_sib_ce_ch4]
    permafrost_total[:perm_jul_ce_co2] = permafrost_jules[:perm_jul_ce_co2]
    permafrost_total[:perm_jul_e_co2] = permafrost_jules[:perm_jul_e_co2]
    permafrost_total[:perm_jul_ce_ch4] = permafrost_jules[:perm_jul_ce_ch4]

    set_param!(m, :PermafrostSiBCASA, :rt_g, readpagedata(m,"test/validationdata/$valdir/rt_g_globaltemperature.csv"))
    set_param!(m, :PermafrostJULES, :rt_g, readpagedata(m,"test/validationdata/$valdir/rt_g_globaltemperature.csv"))

    p = load_parameters(m)
    p["y_year_0"] = 2015.
    p["y_year"] = Mimi.dim_keys(m.md, :time)
    set_leftover_params!(m, p)

    run(m)

    permjcco2 = m[:PermafrostJULES,  :perm_jul_ce_c_co2]
    permjcco2_compare = readpagedata(m,"test/validationdata/$valdir/perm_jul_ce_c_co2.csv")
    @test permjcco2 ≈ permjcco2_compare rtol=1e-4

    permsco2 = m[:PermafrostSiBCASA,  :perm_sib_ce_co2]
    permsco2_compare = readpagedata(m,"test/validationdata/$valdir/perm_sib_ce_co2.csv")
    @test permsco2 ≈ permsco2_compare rtol=1e-4

    permco2 = m[:PermafrostTotal,  :perm_tot_e_co2]
    permco2_compare = readpagedata(m,"test/validationdata/$valdir/perm_tot_e_co2.csv")
    @test permco2 ≈ permco2_compare rtol=1e-4

    permch4 = m[:PermafrostTotal,  :perm_tot_ce_ch4]
    permch4_compare = readpagedata(m,"test/validationdata/$valdir/perm_tot_ce_ch4.csv")
    @test permch4 ≈ permch4_compare rtol=1e-4
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
abstract HopfieldNet

function energy(net::HopfieldNet)
    e = 0.0
    n = length(net.s)
    for i in 1:n
        for j in 1:n
            e += net.W[i, j] * net.s[i] * net.s[j]
        end
    end
    e *= -0.5
    for i in 1:n
        e += net.W[i] * net.s[i]
    end
    return e
end

function settle!(net::HopfieldNet,
                 iterations::Integer = 1_000,
                 trace::Bool = false)
    for i in 1:iterations
        update!(net)
        if trace
            @printf "%5.0d: %.4f\n" i energy(net)
        end
    end
    return
end

function associate!{T <: Real}(net::HopfieldNet,
                               pattern::Vector{T};
                               iterations::Integer = 1_000,
                               trace::Bool = false)
    copy!(net.s, pattern)
    settle!(net, iterations, trace)
    # TODO: Decide if this should really be a copy
    return copy(net.s)
end

# Hebbian learning steps w/ columns as patterns
function train!{T <: Real}(net::HopfieldNet, patterns::Matrix{T})
    n = length(net.s)
    p = size(patterns, 2)
    # Could use outer products here
    # (1 / p) * (patterns[:, mu] * patterns[:, mu]')
    for i in 1:n
        for j in (i + 1):n
            s = 0.0
            for mu in 1:p
                s += patterns[i, mu] * patterns[j, mu]
            end
            s = s / p # May need to be careful here
            net.W[i, j] += s
            net.W[j, i] += s
        end
    end
    return
end

function h{T <: Real}(i::Integer, j::Integer, mu::Integer, n::Integer,
                      W::Matrix{Float64}, patterns::Matrix{T})
    res = 0.0
    for k in 1:n
        if k != i && k != j
            res += W[i, k] * patterns[k, mu]
        end
    end
    return res
end

# Storkey learning steps w/ columns as patterns
function storkeytrain!{T <: Real}(net::HopfieldNet, patterns::Matrix{T})
    n = length(net.s)
    p = size(patterns, 2)
    for i in 1:n
        for j in (i + 1):n
            for mu in 1:p
                s = patterns[i, mu] * patterns[j, mu]
                s -= patterns[i, mu] * h(j, i, mu, n, net.W, patterns)
                s -= h(i, j, mu, n, net.W, patterns) * patterns[j, mu]
                s *= 1 / n
                net.W[i, j] += s
                net.W[j, i] += s
            end
        end
    end
    return
end


text:
function sma{T,N}(ta::TimeArray{T,N}, n::Int)
  tstamps = ta.timestamp[n:end]
  vals    = zeros(length(ta) - (n-1))
  cname   = ["sma$n"]
  for i in 1:length(ta) - (n-1)
    vals[i] =  mean(ta.values[i:i+(n-1)])
  end
  TimeArray(tstamps, vals, cname, ta.meta)
end

function ema{T,N}(ta::TimeArray{T,N}, n::Int; wilder=false)

  if  wilder
    k  = 1/n
  else
    k = 2/(n+1)
  end

  tstamps = ta.timestamp[n:end]
  vals    = ones(length(ta))
  vals[n] = sma(ta, n).values[1] # seed with first value an sma value
  for i = n+1:length(ta)
    vals[i] =  ta.values[i] * k + vals[i-1] * (1-k)
  end
  cname   = ["ema$n"]
  TimeArray(tstamps, vals[n:length(ta)], cname, ta.meta)
end

# Array dispatch for use by ta algorithms

function sma(a::Array{Float64,1}, n::Int)
  vals    = zeros(length(a) - (n-1))
  for i in 1:length(a) - (n-1)
    vals[i] =  mean(a[i:i+(n-1)])
  end
  vals
end

function ema(a::Array{Float64,1}, n::Int; wilder=false)

  if  wilder
    k  = 1/n
  else
    k = 2/(n+1)
  end

  vals    = ones(length(a))
  vals[n] = sma(a, n)[1] # seed with first value an sma value

  for i = n+1:length(a)
    vals[i] =  a[i] * k + vals[i-1] * (1-k)
  end
  vals[n:end]
end


text:
using PkgBenchmark
benchmarkpkg(
    dirname(@__DIR__),
    BenchmarkConfig(env = Dict("JULIA_NUM_THREADS" => "1", "OMP_NUM_THREADS" => "1")),
    resultfile = joinpath(@__DIR__, "result.json"),
)


text:
 """
```
pseudo_measurement(m::Model1002{T},
    TTT::Matrix{T}, RRR::Matrix{T}, CCC::Vector{T}) where {T<:AbstractFloat}
```

Assign pseudo-measurement equation (a linear combination of states):

```
x_t = ZZ_pseudo*s_t + DD_pseudo
```
"""
function pseudo_measurement(m::Model1002{T},
                            TTT::Matrix{T},
                            RRR::Matrix{T},
                            CCC::Vector{T}) where {T<:AbstractFloat}

    endo      = m.endogenous_states
    endo_addl = m.endogenous_states_augmented
    pseudo    = m.pseudo_observables

    _n_states = n_states_augmented(m)
    _n_pseudo = n_pseudo_observables(m)

    # Compute TTT^10, used for Expected10YearRateGap, Expected10YearRate, and Expected10YearNaturalRate
    TTT10 = (1/40)*((UniformScaling(1.) - TTT)\(TTT - TTT^41))

    # Initialize pseudo ZZ and DD matrices
    ZZ_pseudo = zeros(_n_pseudo, _n_states)
    DD_pseudo = zeros(_n_pseudo)

    ##########################################################
    ## PSEUDO-OBSERVABLE EQUATIONS
    ##########################################################

    ## Output
    ZZ_pseudo[pseudo[:y_t],endo[:y_t]] = 1.

    ## Flexible Output
    ZZ_pseudo[pseudo[:y_f_t],endo[:y_f_t]] = 1.

    ## Natural Rate
    ZZ_pseudo[pseudo[:NaturalRate],endo[:r_f_t]] = 1.
    DD_pseudo[pseudo[:NaturalRate]]              = 100.0*(m[:rstar]-1.0)

    ## π_t
    ZZ_pseudo[pseudo[:π_t],endo[:π_t]] = 1.
    DD_pseudo[pseudo[:π_t]]            = 100*(m[:π_star]-1);

    ## Output Gap
    ZZ_pseudo[pseudo[:OutputGap],endo[:y_t]] = 1;
    ZZ_pseudo[pseudo[:OutputGap],endo[:y_f_t]] = -1;

    ## Ex Ante Real Rate
    ZZ_pseudo[pseudo[:ExAnteRealRate],endo[:R_t]]  = 1;
    ZZ_pseudo[pseudo[:ExAnteRealRate],endo[:Eπ_t]] = -1;
    DD_pseudo[pseudo[:ExAnteRealRate]]             = m[:Rstarn] - 100*(m[:π_star]-1);

    ## Long Run Inflation
    ZZ_pseudo[pseudo[:LongRunInflation],endo[:π_star_t]] = 1.
    DD_pseudo[pseudo[:LongRunInflation]]                 = 100. *(m[:π_star]-1.)

    ## Marginal Cost
    ZZ_pseudo[pseudo[:MarginalCost],endo[:mc_t]] = 1.

    ## Wages
    ZZ_pseudo[pseudo[:Wages],endo[:w_t]] = 1.

    ## Flexible Wages
    ZZ_pseudo[pseudo[:FlexibleWages],endo[:w_f_t]] = 1.

    ## Hours
    ZZ_pseudo[pseudo[:Hours],endo[:L_t]] = 1.

    ## Flexible Hours
    ZZ_pseudo[pseudo[:FlexibleHours],endo[:L_f_t]] = 1.

    ## z_t
    ZZ_pseudo[pseudo[:z_t], endo[:z_t]] = 1.

    ## Expected 10-Year Rate Gap
    ZZ_pseudo[pseudo[:Expected10YearRateGap], :] = TTT10[endo[:R_t], :] - TTT10[endo[:r_f_t], :] - TTT10[endo[:Eπ_t], :]

    ## Nominal FFR
    ZZ_pseudo[pseudo[:NominalFFR], endo[:R_t]] = 1.
    DD_pseudo[pseudo[:NominalFFR]] = m[:Rstarn]

    ## Expected 10-Year Interest Rate
    ZZ_pseudo[pseudo[:Expected10YearRate], :] = TTT10[endo[:R_t], :]
    DD_pseudo[pseudo[:Expected10YearRate]]    = m[:Rstarn]

    ## Expected 10-Year Natural Rate
    ZZ_pseudo[pseudo[:Expected10YearNaturalRate], :] = TTT10[endo[:r_f_t], :] + TTT10[endo[:Eπ_t], :]
    DD_pseudo[pseudo[:Expected10YearNaturalRate]]    = m[:Rstarn]

    ## Expected Nominal Natural Rate
    ZZ_pseudo[pseudo[:ExpectedNominalNaturalRate], endo[:r_f_t]] = 1.
    ZZ_pseudo[pseudo[:ExpectedNominalNaturalRate], endo[:Eπ_t]]  = 1.
    DD_pseudo[pseudo[:ExpectedNominalNaturalRate]]               = m[:Rstarn]

    ## Nominal Rate Gap
    ZZ_pseudo[pseudo[:NominalRateGap], endo[:R_t]]   = 1.
    ZZ_pseudo[pseudo[:NominalRateGap], endo[:r_f_t]] = -1.
    ZZ_pseudo[pseudo[:NominalRateGap], endo[:Eπ_t]]  = -1.

    ## Labor Productivity Growth
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo[:y_t]]           = 1.
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo_addl[:y_t1]]     = -1.
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo[:z_t]]           = 1.
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo_addl[:e_gdp_t]]  = 1.
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo_addl[:e_gdp_t1]] = -m[:me_level]
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo[:L_t]]           = -1
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo_addl[:L_t1]]     = 1.
    DD_pseudo[pseudo[:LaborProductivityGrowth]]                       = 100*(exp(m[:z_star]) - 1)

    ## u_t
    ZZ_pseudo[pseudo[:u_t], endo[:u_t]] = 1.

    ## Fundamental inflation related pseudo-obs
    if subspec(m) in ["ss13", "ss14", "ss15", "ss16", "ss17", "ss18", "ss19"]
        # Compute coefficient on Sinf
        betabar = exp((1-m[:σ_c] ) * m[:z_star]) * m[:β]
        κ = ((1 - m[:ζ_p]*m[:β]*exp((1 - m[:σ_c])*m[:z_star]))*
             (1 - m[:ζ_p]))/(m[:ζ_p]*((m[:Φ]- 1)*m[:ϵ_p] + 1))/
        (1 + m[:ι_p]*m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
        κcoef = κ * (1 + m[:ι_p] * betabar)

        ZZ_pseudo[pseudo[:Sinf_t], endo_addl[:Sinf_t]] = 1.
        ZZ_pseudo[pseudo[:Sinf_w_coef_t], endo_addl[:Sinf_t]] = κcoef
        DD_pseudo[pseudo[:ι_p]] = m[:ι_p]
        ZZ_pseudo[pseudo[:πtil_t], endo_addl[:πtil_t]] = 1.
        DD_pseudo[pseudo[:πtil_t]] = 100 * (m[:π_star] - 1)
        ZZ_pseudo[pseudo[:e_tfp_t], endo_addl[:e_tfp_t]] = 1.
        if subspec(m) in ["ss14", "ss15", "ss16", "ss18", "ss19"]
            ZZ_pseudo[pseudo[:e_tfp_t1], endo_addl[:e_tfp_t1]] = 1.
        end
    end

    ## Exogenous processes
    if subspec(m) == "ss12"
        to_add = [:g_t, :b_t, :μ_t, :z_t, :λ_f_t, :λ_w_t, :rm_t, :σ_ω_t, :μ_e_t,
                  :γ_t, :π_star_t]
        to_add_addl = [:e_lr_t, :e_tfp_t, :e_gdpdef_t, :e_corepce_t, :e_gdp_t, :e_gdi_t]
        for i in to_add
            ZZ_pseudo[pseudo[i], endo[i]] = 1.
        end
        for i in to_add_addl
            ZZ_pseudo[pseudo[i], endo_addl[i]] = 1.
        end
    end

    return PseudoMeasurement(ZZ_pseudo, DD_pseudo)
end


text:
module FDM

    using Printf, LinearAlgebra

    const AV = AbstractVector

    include("methods.jl")
    include("numerics.jl")
    include("grad.jl")
end


text:
# Add MDP to Julia search path
push!(LOAD_PATH, MDPpath * "src/Julia");
push!(LOAD_PATH, MDPpath * "src/Julia/Gencode");
push!(LOAD_PATH, MDPpath * "src/Julia/Preprocessing");
push!(LOAD_PATH, MDPpath * "src/Julia/Postprocessing");
push!(LOAD_PATH, MDPpath * "src/Julia/Potential");
push!(LOAD_PATH, MDPpath * "src/Julia/Optimization");
push!(LOAD_PATH, MDPpath * "src/Julia/Integration");

# Set Julia's PATH enviroment variable so that MDP can call external programs
ENV["PATH"] = "/usr/local/bin:/usr/bin:/opt/local/bin:/bin:/usr/sbin:/sbin";

using Preprocessing, Potential, Optimization, Integration

# load Snap library
Potential.loadsnap(MDPpath)

# load integration library
Integration.loadintegration(MDPpath)

# initialize arrays with nothing
traindata = Array{Any}(nothing, 100)
testdata = Array{Any}(nothing, 100)
validdata = Array{Any}(nothing, 100)
descriptors = Array{Any}(nothing, 100)
potentials  = Array{Any}(nothing, 100)

# initialize app
app = Preprocessing.initializeapp(cdir,MDPpath);

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
## ---  Define a StatiCompiler- and LLVM-compatible statically-sized string type

    # Define the StaticString type, backed by an NTuple

    # Definition and constructors:
    """
    ```julia
    StaticString{N}
    ```
    A stringy type which should generally behave like a base Julia `String`, but
    is explicitly null-terminated, mutable, and standalone-StaticCompiler safe
    (does not require libjulia).

    Can be constructed with the `c"..."` string macro.

    Unlike base Julia `String`s, slicing does not create a copy, but rather a view.
    You are responsible for ensuring that any such views are null-terminated if you
    wish to pass them to any functions (including most system IO) that expect
    null-termination.

    ## Examples
    ```julia
    julia> s = c"Hello world!"
    c"Hello world!"

    julia> s[8:12] = c"there"; s
    c"Hello there!"

    julia> s[1:5]
    StringView: "Hello"

    julia> s[1:5] == "Hello"
    true

    julia> StaticString(s[1:5])
    c"Hello"
    ```
    """
    mutable struct StaticString{N} <: AbstractStaticString
        data::NTuple{N,UInt8}
        @inline function StaticString{N}(::UndefInitializer) where N
            new{N}()
        end
        @inline function StaticString(data::NTuple{N,UInt8}) where N
            new{N}(data)
        end
    end

    """
    ```julia
    StaticString{N}(undef)
    ```
    Construct an uninitialized `N`-byte `StaticString`

    ```julia
    StaticString(data::NTuple{N,UInt8})
    ```
    Construct a `StaticString` containing the `N` bytes specified by `data`.
    To yield a valid string, `data` must be null-terminated, i.e., end in `0x00`.

    ```julia
    StaticString(s::AbstractStaticString)
    ```
    Construct a `StaticString` containing the same data as the input string `s`.

    ## Examples
    ```julia
    julia> data = (0x48, 0x65, 0x6c, 0x6c, 0x6f, 0x20, 0x77, 0x6f, 0x72, 0x6c, 0x64, 0x21, 0x00);

    julia> s = StaticString(data)
    c"Hello world!"

    julia> StaticString(s[1:5])
    c"Hello"
    ```
    """
    @inline function StaticString(s::AbstractStaticString)
        c = StaticString{length(s)+1}(undef)
        c[1:length(s)] = s
        c[end] = 0x00
        return c
    end

    # String macro to create null-terminated `StaticString`s
    """
    ```julia
    @c_str -> StaticString
    ```
    Construct a `StaticString`, such as `c"Foo"`.

    A `StaticString` should generally behave like a base Julia `String`, but is
    explicitly null-terminated, mutable, and standalone-StaticCompiler safe (does
    not require libjulia).

    ## Examples
    ```julia
    julia> c"Hello there!"
    c"Hello there!"

    julia> c"foo" == "foo"
    true
    ```
    """
    macro c_str(s)
        n = _unsafe_unescape!(s)
        t = Expr(:tuple, codeunits(s)[1:n]..., 0x00)
        :(StaticString($t))
    end

    # Fundamentals -- where overriding AbstractStaticString defaults
    @inline Base.unsafe_convert(::Type{Ptr{T}}, m::StaticString) where {T} = Ptr{T}(pointer_from_objref(m))
    @inline Base.pointer(m::StaticString{N}) where {N} = Ptr{UInt8}(pointer_from_objref(m))
    @inline Base.length(s::StaticString{N}) where N = N-1
    @inline Base.:(==)(::StaticString, ::StaticString) = false
    @inline function Base.:(==)(a::StaticString{N}, b::StaticString{N}) where N
        GC.@preserve a b begin
            pa, pb = pointer(a), pointer(b)
            for n ∈ 0:N-1
                unsafe_load(pa + n) == unsafe_load(pb + n) || return false
            end
            return true
        end
    end

    # Custom replshow for interactive use (n.b. _NOT_ static-compilerable)
    function Base.show(io::IO, s::StaticString)
        Base.print(io, "c\"")
        Base.escape_string(io, Base.unsafe_string(pointer(s)))
        Base.print(io, "\"")
    end

    # Implement some of the AbstractArray interface -- where overriding AbstractStaticString defaults
    @inline Base.firstindex(s::StaticString) = 1
    @inline Base.lastindex(s::StaticString{N}) where {N} = N
    @inline Base.eachindex(s::StaticString{N}) where {N} = 1:N
    @inline Base.getindex(s::StaticString, i::Int) = unsafe_load(pointer(s)+(i-1))
    @inline Base.copy(s::StaticString) = StaticString(codeunits(s))

    # Implement some of the AbstractString interface -- where overriding AbstractStaticString defaults
    @inline Base.ncodeunits(s::StaticString{N}) where N = N
    @inline Base.codeunits(s::StaticString) = s.data
    @inline function Base.:*(a::AbstractStaticString, b::AbstractStaticString)  # Concatenation
        N = length(a) + length(b) + 1 # n.b. `length` excludes null-termination
        c = StaticString{N}(undef)
        c[1:length(a)] = a
        c[length(a)+1:length(a)+length(b)] = b
        c[end] = 0x00 # Null-terminate
        return c
    end
    @inline function Base.:^(s::AbstractStaticString, n::Integer)       # Repetition
        l = length(s) # Excluding null-termination
        N = n*l + 1
        c = StaticString{N}(undef)
        for i=1:n
            c[(l*(i-1) + 1):(l*i)] = s
        end
        c[end] = 0x00 # Null-terminate
        return c
    end

    # String macro to directly create null-terminated `ManualMemory.MemoryBuffer`s
    macro mm_str(s)
        n = _unsafe_unescape!(s)
        t = Expr(:tuple, codeunits(s)[1:n]...)
        :(MemoryBuffer($t))
    end


text:
using DataFrames: DataFrame
using Dates
using PlutoStaticHTML
using Pluto:
    Cell,
    Notebook,
    Pluto,
    ServerSession,
    SessionActions
using Test
using TimerOutputs: TimerOutput, @timeit

const PKGDIR = string(pkgdir(PlutoStaticHTML))::String
const NOTEBOOK_DIR = joinpath(PKGDIR, "docs", "src", "notebooks")

function pluto_notebook_content(code)
    return """
        ### A Pluto.jl notebook ###
        # v0.18.1

        using Markdown
        using InteractiveUtils

        # ╔═╡ a6dda572-3f2c-11ec-0eeb-69e2323a92de
        $(code)

        # ╔═╡ Cell order:
        # ╠═a6dda572-3f2c-11ec-0eeb-69e2323a92de
        """
end

function drop_cache_info(html::AbstractString)
    n = PlutoStaticHTML.n_cache_lines()
    sep = '\n'
    lines = split(html, sep)
    return join(lines[n:end], sep)
end

function drop_begin_end(html::AbstractString)
    sep = '\n'
    lines = split(html, sep)
    return join(lines[2:end-1], sep)
end

"Helper function to simply pass a `nb::Notebook` and run it."
function notebook2html_helper(
        nb::Notebook,
        hopts=HTMLOptions();
        use_distributed::Bool=true
    )
    tmpdir = mktempdir()
    tmppath = joinpath(tmpdir, "notebook.jl")
    Pluto.save_notebook(nb, tmppath)
    session = ServerSession()
    session.options.evaluation.workspace_use_distributed = use_distributed
    nb = PlutoStaticHTML.run_notebook!(tmppath, session)
    html = PlutoStaticHTML.notebook2html(nb, tmppath, hopts)

    has_begin_end = contains(html, PlutoStaticHTML.BEGIN_IDENTIFIER)
    without_begin_end = has_begin_end ? drop_begin_end(html) : html

    # Remove the caching information because it's not important for most tests.
    has_cache = contains(html, PlutoStaticHTML.STATE_IDENTIFIER)
    without_cache = has_cache ? drop_cache_info(without_begin_end) : html

    return (without_cache, nb)
end

# Credits to Tensors.jl/test/runtests.jl
macro timed_testset(str, block)
    return quote
        @timeit TIMEROUTPUT "$($(esc(str)))" begin
            @testset "$($(esc(str)))" begin
                $(esc(block))
            end
        end
    end
end

# Hide output when using `TestEnv.activate(); include("test/preliminaries.jl")`.
nothing


text:
using MultivariateFunctions
using Dates
tol = 1e-13

tol = 10*eps()
constant    = PE_Unit()
linear      = PE_Unit(0.0,2.0,1)
quadratic   = PE_Unit(0.0,1.0,2)
exponential = PE_Unit(2.0,2.0,0)
higher      = PE_Unit(2.0,2.0,2)

f = PE_Function(1.8, Dict([:w, :x, :y, :z] .=> [linear, quadratic, exponential, higher]))

degree = 2
dimensions = underlying_dimensions(f)
g = all_derivatives(f, degree, dimensions)
coordinates = Dict{Symbol,Float64}([:w, :x, :y, :z] .=> [1.1, 2.1, 3.3, 4.0])
abs(evaluate(g[Dict{Symbol,Int}([:w,:z] .=> [1,1])], coordinates) - evaluate(derivative(f, Dict{Symbol,Int}([:w,:z] .=> [1,1])), coordinates) ) < 1e-10

# Testing uniroot
step_size = 1.0
max_iters = 40
convergence_tol = 1e-10
print_reports = false
initial_guess = coordinates
root_coordinates, root_value, root_convergence  = uniroot(f, initial_guess; step_size = step_size, max_iters = max_iters, convergence_tol = convergence_tol, print_reports = print_reports)
evaluate(f, root_coordinates) < tol
initial_guess = Dict{Symbol,Float64}([:w, :x, :y, :z] .=> [10.0, 5.0, 0.2, 1.2])
root_coordinates2, root_value2, root_convergence2 = uniroot(f, initial_guess; step_size = step_size, max_iters = max_iters, convergence_tol = convergence_tol, print_reports = print_reports)
evaluate(f, root_coordinates2) < tol
# Test we have found a different root:
sum(values(merge(-,root_coordinates, root_coordinates2))) > 1.0

print_reports = false
maximise = true
initial_guess = coordinates
opt_coordinates , opt_value , opt_convergence, hessian_det_sign  = find_local_optima(f, initial_guess; step_size = step_size, max_iters = max_iters, convergence_tol = convergence_tol, print_reports = print_reports)
opt_convergence
coordinates2 = Dict{Symbol,Float64}([:w, :x, :y, :z] .=> [1.9, 4.1, -4.0, 4.0])
opt_coordinates2, opt_value2, opt_convergence2, hessian_det_sign2 = find_local_optima(f, coordinates2; step_size = step_size, max_iters = max_iters, convergence_tol = convergence_tol, print_reports = print_reports)
opt_convergence2


text:
function divedivedive(d::Int)
    try
        divedivedive(d+1)
    catch
        return d
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# ---
# title: 194. Transpose File
# id: problem194
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories:
# link: <https://leetcode.com/problems/transpose-file/description/>
# hidden: true
# ---
#
# Given a text file `file.txt`, transpose its content.
#
# You may assume that each row has the same number of columns and each field is
# separated by the `' '` character.
#
# **Example:**
#
# If `file.txt` has the following content:
#
#
#
#     name age
#     alice 21
#     ryan 30
#
#
# Output the following:
#
#
#
#     name alice ryan
#     age 21 30
#
#
#
#
#
## @lc code=start
using LeetCode

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


text:
# this file tests various aspects of the tracers routines
using MagmaThermoKinematics
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
using Plots
using LinearAlgebra
using SpecialFunctions
using Test

const CreatePlots = true      # easy way to deactivate plotting throughout

function test_TracerUpdate(Dimension="2D", InterpolationMethod="Linear")
  # test interpolation methods from grid to tracers in 2D and 3D

  if Dimension=="2D"
    # Model parameters
    W,H                     =   1.,  1.;                                    # Width, Length, Height

    # Define grid
    Nx, Nz                  =   33, 33;                                     # resolution of coarse grid
    dx,dz                   =   W/(Nx-1), H/(Nz-1);                         # grid size [m]
    x,z                     =   0:dx:((Nx-1)*dx), 0:dz:((Nz-1)*dz);         # 1D coordinate arrays
    coords                  =   collect(Iterators.product(x,z))             # generate coordinates from 1D coordinate vectors
    X,Z                     =   (x->x[1]).(coords), (x->x[2]).(coords);     # transfer coords to 3D arrays
    Grid,FullGrid,Spacing   =   (x,z), (X,Z), (dx,dz);

    # Define function on grid
    T                       =   cos.(pi.*X).*sin.(2*pi.*Z)

  elseif Dimension=="3D"
      # Model parameters
      W,L,H                 =   1., 1., 1.;                                    # Width, Length, Height

      # Define grid
      Nx, Ny, Nz              =   33,33, 33;                                                    # resolution of coarse grid
      dx,dy,dz                =   W/(Nx-1), L/(Ny-1), H/(Nz-1);                                 # grid size [m]
      x,y,z                   =   0:dx:((Nx-1)*dx),  0:dy:((Ny-1)*dy), 0:dz:((Nz-1)*dz);        # 1D coordinate arrays
      coords                  =   collect(Iterators.product(x,y,z))                             # generate coordinates from 1D coordinate vectors
      X,Y,Z                   =   (x->x[1]).(coords), (x->x[2]).(coords), (x->x[3]).(coords);   # transfer coords to 3D arrays
      Grid, FullGrid, Spacing =   (x,y,z), (X,Y,Z), (dx,dy,dz);

      # Define function on coarse grid
      T                       =   cos.(pi.*X).*sin.(2*pi.*Z).*sin.(2*pi.*Y)

  end

  # Create tracer structure that fill te full grid
  Tracers                     =   InitializeTracers(FullGrid,3, false);

  Phi     =   Z./H;

  # Perform interpolation from grid -> tracers
  Tracers =  UpdateTracers(Tracers, Grid, T, Phi, InterpolationMethod);

  # Compute error
  Tr_coord    =   Tracers.coord; Tr_coord = hcat(Tr_coord...)';       # extract array with coordinates of tracers

  if Dimension=="2D"
    Tanal       =   cos.(pi.*Tr_coord[:,1]).*sin.(2*pi.*Tr_coord[:,2]);
    Phi_anal    =    1 .- Tr_coord[:,end]/H;
    Error       =   (Tanal - Tracers.T)  .+ (Phi_anal-Tracers.Phi_melt);

    if CreatePlots
      p1          =   contourf(x, z,      T',       aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Grid",  dpi=300, levels=10)
      p2          =   scatter(Tr_coord[:,1], Tr_coord[:,2], zcolor = Tracers.T, m = (:inferno , 0.8, Plots.stroke(0.01, :black)), title="Tracers", markersize=2.0,dpi=300)

    #   p1          =   contourf(x, z,      (1.0-Phi)',       aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Grid",  dpi=150, levels=10)
    #   p2          =   scatter(Tr_coord[:,1], Tr_coord[:,2], zcolor = Tracers.Phi_melt, m = (:inferno , 0.8, Plots.stroke(0.01, :black)), title="Tracers", markersize=5.0)

      plot(p1,p2);

     png("TracerUpdate_2D_$InterpolationMethod")
    end

  elseif Dimension=="3D"
    Tanal       =   cos.(pi.*Tr_coord[:,1]).*sin.(2*pi.*Tr_coord[:,3]).*sin.(2*pi.*Tr_coord[:,2]);
    Phi_anal    =   1 .- Tr_coord[:,end]/H;
    Error       =   (Tanal - Tracers.T)  .+ (Phi_anal-Tracers.Phi_melt);

    if CreatePlots
    #  p1          =   contourf(x, z,      Data_coarse1[1][:,10,:]', aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Grid",  dpi=150, levels=10)
    #  p2          =   scatter(Tr_coord[:,1], Tr_coord[:,2], zcolor = Tracers.T, m = (:inferno , 0.8, Plots.stroke(0.01, :black)), title="Tracers", markersize=5.0)

   #   plot(p1,p2);
   #   png("TracerUpdate_2D_$InterpolationMethod")
    end

  end

  error = norm(Error[:],2)/length(Error[:]);
  return error;        # return error
end

function test_PhaseRatioFromTracers(Dimension="2D", InterpolationMethod="Linear", Method="TracersEverywhere")
  # test routines to go from tracers -> Grid

  if Dimension=="2D"
    dim                     = 2;
    # Model parameters
    W,H                     =   1.,  1.;                                    # Width, Length, Height

    # Define grid
    Nx, Nz                  =   65, 65;                                   # resolution of coarse grid
    dx,dz                   =   W/(Nx-1), H/(Nz-1);                         # grid size [m]
    x,z                     =   0:dx:((Nx-1)*dx), 0:dz:((Nz-1)*dz);         # 1D coordinate arrays
    coords                  =   collect(Iterators.product(x,z))             # generate coordinates from 1D coordinate vectors
    X,Z                     =   (x->x[1]).(coords), (x->x[2]).(coords);     # transfer coords to 3D arrays
    Grid,FullGrid,Spacing   =   (x,z), (X,Z), (dx,dz);

    # Define function on grid
    T                       =   cos.(pi.*X).*sin.(2*pi.*Z)

  elseif Dimension=="3D"
      dim                   = 3;
      # Model parameters
      W,L,H                 =   1., 1., 1.;                                    # Width, Length, Height

      # Define grid
      Nx, Ny, Nz              =   33,33, 33;                                                    # resolution of coarse grid
      dx,dy,dz                =   W/(Nx-1), L/(Ny-1), H/(Nz-1);                                 # grid size [m]
      x,y,z                   =   0:dx:((Nx-1)*dx),  0:dy:((Ny-1)*dy), 0:dz:((Nz-1)*dz);        # 1D coordinate arrays
      coords                  =   collect(Iterators.product(x,y,z))                             # generate coordinates from 1D coordinate vectors
      X,Y,Z                   =   (x->x[1]).(coords), (x->x[2]).(coords), (x->x[3]).(coords);   # transfer coords to 3D arrays
      Grid, FullGrid, Spacing =   (x,y,z), (X,Y,Z), (dx,dy,dz);

      # Define function on coarse grid
      T                       =   cos.(pi.*X).*sin.(2*pi.*Z).*sin.(2*pi.*Y)

  end

  # Create tracer structure that fill the full grid
  #println("InitializeTracers:")
  if      Method=="TracersEverywhere"
    Tracers                     =   InitializeTracers(FullGrid,3, false);   # tracers defined everywhere

  elseif  Method=="TracersLimitedRegion"
    # Tracers defined in a limited region only, with a default phase outside this area

    # we use the same routine to initialize tracers but in a smaller square region
    W_l, H_l      =   W/3.0, H/3.0;
    dx_l,dz_l     =   W_l/(Nx-1), H_l/(Nz-1);                         # grid size [m]
    xl,zl         =   0.3:dx_l:((Nx-1)*dx_l)+0.3, 0.3:dz_l:((Nz-1)*dz_l)+0.3;         # 1D coordinate arrays of limited region
    if dim==3
      L_l         =   L/3.0;
      dy_l        =   L_l/(Ny-1);
      yl          =   0.3:dy_l:((Ny-1)*dy_l)+0.3;
      coords      =   collect(Iterators.product(xl,yl,zl))
      X,Y,Z       =   (x->x[1]).(coords), (x->x[2]).(coords), (x->x[3]).(coords);   # transfer coords to 3D arrays
      FullGrid_l  =   (X,Y,Z)

    else
      coords      =   collect(Iterators.product(xl,zl))
      X,Z         =   (x->x[1]).(coords), (x->x[2]).(coords);   # transfer coords to 2D arrays
      FullGrid_l  =   (X,Z)
    end
    Tracers       =   InitializeTracers(FullGrid_l,2, false);   # tracers defined everywhere

  else
    error("unknown Method = $Method")
  end

  # Set different phases/temperatures on tracers (default Phase = 0)
  # Sphere 1 (Phase=2)
  if      dim==2; cen = [0.5; 0.5];
  elseif  dim==3; cen = [0.5; 0.5; 0.5]; end
  R = 0.1;
  [Tracers.Phase[i]=2  for i=1:length(Tracers) if sum( (Tracers.coord[i]-cen).^2)<R^dim ]   # julia iterator to set props on Tracers

  if      dim==2; cen = [0.15; 0.7];
  elseif  dim==3; cen = [0.15; 0.1; 0.7]; end
  R = 0.1;
  [Tracers.Phase[i]=3  for i=1:length(Tracers) if sum( (Tracers.coord[i]-cen).^2)<R^dim ]

  # Call main routine to compute phase fractions from particles
  #println("PhaseRatioFromTracers:")
  if      Method=="TracersEverywhere"
    PhaseRatio, NumTracers  =   PhaseRatioFromTracers(FullGrid, Grid, Tracers, InterpolationMethod, true);
  else
    PhaseRatio,  NumTracers  =   PhaseRatioFromTracers(FullGrid, Grid, Tracers, InterpolationMethod, true, BackgroundPhase=4);
  end

  # Also test the Rocktype routine
  RockType    =   RockAssemblage(PhaseRatio);

  Tr_coord    =   Tracers.coord; Tr_coord = hcat(Tr_coord...)';       # extract array with coordinates of tracers
  if Dimension=="2D"
    Data =   PhaseRatio[:,:,1];
    if CreatePlots
     # @show size(PhaseRatio)
      #p1          =   contourf(x, z,      PhaseRatio[:,:,3]',       aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Grid",  dpi=300, levels=10)
      #p1          =   heatmap(x, z,      PhaseRatio[:,:,3]',       aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Grid",  dpi=300)

       p1          =   heatmap(x, z,      RockType',       aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="RockType",  dpi=300)
      #p1          =   heatmap(x, z,      NumTracers',       aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="NumTracers",  dpi=300)

      #p2 = plot(X[:],Z[:],markershape = :plus, markersize=0.2)
      p2          =   scatter(Tr_coord[:,1], Tr_coord[:,2], zcolor = Tracers.T, m = (:inferno , 0.8, Plots.stroke(0.01, :black)), title="Tracers", markersize=1.0,dpi=300)

    #   p1          =   contourf(x, z,      (1.0-Phi)',       aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Grid",  dpi=150, levels=10)
    #   p2          =   scatter(Tr_coord[:,1], Tr_coord[:,2], zcolor = Tracers.Phi_melt, m = (:inferno , 0.8, Plots.stroke(0.01, :black)), title="Tracers", markersize=5.0)

      #plot(p1,p2);
      plot(p1);

      png("PhaseRatioFromTracers_2D_$(InterpolationMethod)_$(Method)")
    end

  elseif Dimension=="3D"
    Data =   PhaseRatio[:,:,:,1];

    if CreatePlots
    #  p1          =   contourf(x, z,      Data_coarse1[1][:,10,:]', aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Grid",  dpi=150, levels=10)
    #  p2          =   scatter(Tr_coord[:,1], Tr_coord[:,2], zcolor = Tracers.T, m = (:inferno , 0.8, Plots.stroke(0.01, :black)), title="Tracers", markersize=5.0)

   #   plot(p1,p2);
   #   png("Tracer2Grid_2D_$InterpolationMethod")
    end

  end

  return norm(RockType);        # return norm of data
end

function test_TracerToGrid(Dimension="2D")
  # tests interpolating TracersToGrid! routine, which interpolates, e.g. temperature from tracers to the grid

  if Dimension=="2D"
    # Model parameters
    W,H                     =   1.,  1.;                                    # Width, Length, Height

    # Define grid
    Nx, Nz                  =   33, 33;                                     # resolution of coarse grid
    dx,dz                   =   W/(Nx-1), H/(Nz-1);                         # grid size [m]
    x,z                     =   0:dx:((Nx-1)*dx), 0:dz:((Nz-1)*dz);         # 1D coordinate arrays
    coords                  =   collect(Iterators.product(x,z))             # generate coordinates from 1D coordinate vectors
    X,Z                     =   (x->x[1]).(coords), (x->x[2]).(coords);     # transfer coords to 3D arrays
    Grid,FullGrid,Spacing   =   (x,z), (X,Z), (dx,dz);

    # Define function on grid
    T                       =   cos.(pi.*X).*sin.(2*pi.*Z)

  elseif Dimension=="3D"
      # Model parameters
      W,L,H                 =   1., 1., 1.;                                    # Width, Length, Height

      # Define grid
      Nx, Ny, Nz              =   33,33, 33;                                                    # resolution of coarse grid
      dx,dy,dz                =   W/(Nx-1), L/(Ny-1), H/(Nz-1);                                 # grid size [m]
      x,y,z                   =   0:dx:((Nx-1)*dx),  0:dy:((Ny-1)*dy), 0:dz:((Nz-1)*dz);        # 1D coordinate arrays
      coords                  =   collect(Iterators.product(x,y,z))                             # generate coordinates from 1D coordinate vectors
      X,Y,Z                   =   (x->x[1]).(coords), (x->x[2]).(coords), (x->x[3]).(coords);   # transfer coords to 3D arrays
      Grid, FullGrid, Spacing =   (x,y,z), (X,Y,Z), (dx,dy,dz);

      # Define function on coarse grid
      T                       =   cos.(pi.*X).*sin.(2*pi.*Z).*sin.(2*pi.*Y)

  end

  # Create tracer structure that fill te full grid
  Tracers                     =   InitializeTracers(FullGrid,3, false);
  @time Tracers               =   InitializeTracers(FullGrid,3, false);

  Phi     =   Z./H;

  # Perform interpolation from grid -> tracers
  Tracers =  UpdateTracers(Tracers, Grid, T, Phi, "Linear");

  # Go back from tracers to a new grid
  Tnew        =   copy(T)*.0;
  NumTracers  =   TracersToGrid!(Tnew, FullGrid, Grid, Tracers, "T", "Constant", true);
  NumTracers  =   TracersToGrid!(Tnew, FullGrid, Grid, Tracers, "T", "Constant", true);   # do it a second time, to make sire it doesn't double

  # Compute error
  Tr_coord    =   Tracers.coord; Tr_coord = hcat(Tr_coord...)';       # extract array with coordinates of tracers

  if Dimension=="2D"
    Tanal       =   cos.(pi.*X).*sin.(2*pi.*Z);
    #Phi_anal    =    1 .- Tr_coord[:,end]/H;
    Error       =   (Tanal - Tnew);

    if CreatePlots
      p2          =   contourf(x, z,      Tnew',       aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Tnew",  dpi=300, levels=10)
      p1          =   scatter(Tr_coord[:,1], Tr_coord[:,2], zcolor = Tracers.T, m = (:inferno , 0.8, Plots.stroke(0.01, :black)), title="T on Tracers", markersize=2.0,dpi=300)

    #   p1          =   contourf(x, z,      (1.0-Phi)',       aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Grid",  dpi=150, levels=10)
    #   p2          =   scatter(Tr_coord[:,1], Tr_coord[:,2], zcolor = Tracers.Phi_melt, m = (:inferno , 0.8, Plots.stroke(0.01, :black)), title="Tracers", markersize=5.0)

      plot(p1,p2);

     png("TracerToGrid_2D")
    end

  elseif Dimension=="3D"
    Tanal       =   cos.(pi.*X).*sin.(2*pi.*Z).*sin.(2*pi.*Y);
    #Phi_anal    =   1 .- Tr_coord[:,end]/H;
    Error       =   (Tanal - Tnew);

    if CreatePlots
    #  p1          =   contourf(x, z,      Data_coarse1[1][:,10,:]', aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Grid",  dpi=150, levels=10)
    #  p2          =   scatter(Tr_coord[:,1], Tr_coord[:,2], zcolor = Tracers.T, m = (:inferno , 0.8, Plots.stroke(0.01, :black)), title="Tracers", markersize=5.0)

   #   plot(p1,p2);
   #   png("TracerUpdate_2D_$InterpolationMethod")
    end

  end

  error = norm(Error[:],2)/length(Error[:]);
  return error;        # return error
end

# ===================================================================================================
if 1==0
  @testset "Update Tracer" begin
    @test test_TracerUpdate("2D", "Linear") ≈  2.2699573946787955e-5  atol=1e-8;
    @test test_TracerUpdate("2D", "Cubic")  ≈  4.217098535399923e-7   atol=1e-8;
    @test test_TracerUpdate("3D", "Linear") ≈  2.9095783218223196e-6  atol=1e-8;
    @test test_TracerUpdate("3D", "Cubic")  ≈  3.252909873990004e-8   atol=1e-10;
  end;

  @testset "PhaseRatioFromTracers" begin
    @test test_PhaseRatioFromTracers("2D", "DistanceWeighted")  ≈  74.96665925596525    atol=1e-3;
    @test test_PhaseRatioFromTracers("2D", "Constant")          ≈  75.17978451685       atol=1e-3;
    @test test_PhaseRatioFromTracers("3D", "DistanceWeighted")  ≈  189.66285877841239   atol=1e-3;
    @test test_PhaseRatioFromTracers("3D", "Constant")          ≈  189.68658360569415   atol=1e-3;
    @test test_PhaseRatioFromTracers("2D", "DistanceWeighted","TracersLimitedRegion") ≈  246.4710936398019 atol=1e-3;
  end;

  @testset "TracerToGrid" begin
    @test test_TracerToGrid("2D")  ≈  0.0006078889344375056    atol=1e-8;
    @test test_TracerToGrid("3D")  ≈  0.00010053219987571875   atol=1e-8;
  end;

end

 test_TracerToGrid("2D")
 test_TracerToGrid("3D")


text:
# Test basic STRIPS functionality in a gripper domain
@testset "strips" begin

path = joinpath(dirname(pathof(PDDL)), "..", "test", "strips")

domain = load_domain(joinpath(path, "domain.pddl"))
@test domain.name == :gripper
@test convert(Term, domain.predicates[:room]) == pddl"(room ?r)"

problem = load_problem(joinpath(path, "problem.pddl"))
@test problem.name == Symbol("gripper-problem")
@test problem.objects == @pddl("rooma", "roomb", "ball1", "ball2", "left", "right")

state = initstate(domain, problem)
implementations = [
    "concrete interpreter" => (domain, state),
    "abstract interpreter" => abstracted(domain, state),
    "ground interpreter" => (ground(domain, state), state),
    "concrete compiler" => compiled(domain, state),
    "abstract compiler" => compiled(abstracted(domain), state)
]

@testset "strips ($name)" for (name, (domain, _)) in implementations

    # Test forward execution of plans
    state = initstate(domain, problem)
    state = execute(domain, state, pddl"(pick ball1 rooma left)")
    @test domain[state => pddl"(carry ball1 left)"] ≃ true
    state = execute(domain, state, pddl"(move rooma roomb)")
    @test domain[state => pddl"(robbyat roomb)"] ≃ true
    state = execute(domain, state, pddl"(drop ball1 roomb left)")
    @test domain[state => pddl"(at ball1 roomb)"] ≃ true

    @test satisfy(domain, state, problem.goal) ≃ true

    # Test action availability
    state = initstate(domain, problem)
    @test Set{Term}(available(domain, state)) == Set{Term}(@pddl(
        "(pick ball1 rooma right)", "(pick ball1 rooma left)",
        "(pick ball2 rooma right)", "(pick ball2 rooma left)",
        "(move rooma roomb)", "(move rooma rooma)"
    ))

end

# Test backward regression of plans
state = goalstate(domain, problem)
state = regress(domain, state, pddl"(drop ball1 roomb left)")
@test domain[state => pddl"(carry ball1 left)"] == true
state = regress(domain, state, pddl"(move rooma roomb)")
@test domain[state => pddl"(robbyat rooma)"] == true
state = regress(domain, state, pddl"(pick ball1 rooma left)")
@test domain[state => pddl"(at ball1 rooma)"] == true
@test issubset(state, initstate(domain, problem))

# Test action relevance
state = goalstate(domain, problem)
@test Set{Term}(relevant(domain, state)) == Set{Term}(@pddl(
    "(drop ball1 roomb left)", "(drop ball1 roomb right)",
    "(drop ball1 roomb ball1)", "(drop ball1 roomb ball2)",
    "(drop ball1 roomb rooma)", "(drop ball1 roomb roomb)"
))

end # strips


text:
using ITensors, Test, Random

include("util.jl")

function isingMPO(sites)::MPO
  H = MPO(sites)
  N = length(H)
  link = Vector{Index}(undef, N + 1)
  for n in 1:(N + 1)
    link[n] = Index(3, "Link,Ising,l=$(n-1)")
  end
  for n in 1:N
    s = sites[n]
    ll = link[n]
    rl = link[n + 1]
    H[n] = ITensor(dag(ll), dag(s), s', rl)
    H[n] += setelt(ll[1]) * setelt(rl[1]) * op(sites, "Id", n)
    H[n] += setelt(ll[3]) * setelt(rl[3]) * op(sites, "Id", n)
    H[n] += setelt(ll[2]) * setelt(rl[1]) * op(sites, "Sz", n)
    H[n] += setelt(ll[3]) * setelt(rl[2]) * op(sites, "Sz", n)
  end
  LE = ITensor(link[1])
  LE[3] = 1.0
  RE = ITensor(dag(link[N + 1]))
  RE[1] = 1.0
  H[1] *= LE
  H[N] *= RE
  return H
end

#! format: off
# JuliaFormatter tries to change this line to:
# heisenbergMPO(sites, h::Vector{Float64}; onsite::String="Sz")
# so turn it off for this line.
function heisenbergMPO(sites, h::Vector{Float64}, onsite::String="Sz")::MPO
#! format: on
  H = MPO(sites)
  N = length(H)
  link = Vector{Index}(undef, N + 1)
  for n in 1:(N + 1)
    link[n] = Index(5, "Link,Heis,l=$(n-1)")
  end
  for n in 1:N
    s = sites[n]
    ll = link[n]
    rl = link[n + 1]
    H[n] = ITensor(ll, s, s', rl)
    H[n] += setelt(ll[1]) * setelt(rl[1]) * op(sites, "Id", n)
    H[n] += setelt(ll[5]) * setelt(rl[5]) * op(sites, "Id", n)
    H[n] += setelt(ll[2]) * setelt(rl[1]) * op(sites, "S+", n)
    H[n] += setelt(ll[3]) * setelt(rl[1]) * op(sites, "S-", n)
    H[n] += setelt(ll[4]) * setelt(rl[1]) * op(sites, "Sz", n)
    H[n] += setelt(ll[5]) * setelt(rl[2]) * op(sites, "S-", n) * 0.5
    H[n] += setelt(ll[5]) * setelt(rl[3]) * op(sites, "S+", n) * 0.5
    H[n] += setelt(ll[5]) * setelt(rl[4]) * op(sites, "Sz", n)
    H[n] += setelt(ll[5]) * setelt(rl[1]) * op(sites, onsite, n) * h[n]
  end
  H[1] *= setelt(link[1][5])
  H[N] *= setelt(link[N + 1][1])
  return H
end

function NNheisenbergMPO(sites, J1::Float64, J2::Float64)::MPO
  H = MPO(sites)
  N = length(H)
  link = Vector{Index}(undef, N + 1)
  for n in 1:(N + 1)
    link[n] = Index(8, "Link,H,l=$(n-1)")
  end
  for n in 1:N
    s = sites[n]
    ll = link[n]
    rl = link[n + 1]
    H[n] = ITensor(ll, s, s', rl)
    H[n] += setelt(ll[1]) * setelt(rl[1]) * op(sites, "Id", n)
    H[n] += setelt(ll[8]) * setelt(rl[8]) * op(sites, "Id", n)

    H[n] += setelt(ll[2]) * setelt(rl[1]) * op(sites, "S-", n)
    H[n] += setelt(ll[5]) * setelt(rl[2]) * op(sites, "Id", n)
    H[n] += setelt(ll[8]) * setelt(rl[2]) * op(sites, "S+", n) * J1 / 2
    H[n] += setelt(ll[8]) * setelt(rl[5]) * op(sites, "S+", n) * J2 / 2

    H[n] += setelt(ll[3]) * setelt(rl[1]) * op(sites, "S+", n)
    H[n] += setelt(ll[6]) * setelt(rl[3]) * op(sites, "Id", n)
    H[n] += setelt(ll[8]) * setelt(rl[3]) * op(sites, "S-", n) * J1 / 2
    H[n] += setelt(ll[8]) * setelt(rl[6]) * op(sites, "S-", n) * J2 / 2

    H[n] += setelt(ll[4]) * setelt(rl[1]) * op(sites, "Sz", n)
    H[n] += setelt(ll[7]) * setelt(rl[4]) * op(sites, "Id", n)
    H[n] += setelt(ll[8]) * setelt(rl[4]) * op(sites, "Sz", n) * J1
    H[n] += setelt(ll[8]) * setelt(rl[7]) * op(sites, "Sz", n) * J2
  end
  H[1] *= setelt(link[1][8])
  H[N] *= setelt(link[N + 1][1])
  return H
end

function threeSiteIsingMPO(sites, h::Vector{Float64})::MPO
  H = MPO(sites)
  N = length(H)
  link = Vector{Index}(undef, N + 1)
  for n in 1:(N + 1)
    link[n] = Index(4, "Link,l=$(n-1)")
  end
  for n in 1:N
    s = sites[n]
    ll = link[n]
    rl = link[n + 1]
    H[n] = ITensor(ll, s, s', rl)
    H[n] += setelt(ll[1]) * setelt(rl[1]) * op(sites, "Id", n)
    H[n] += setelt(ll[4]) * setelt(rl[4]) * op(sites, "Id", n)
    H[n] += setelt(ll[2]) * setelt(rl[1]) * op(sites, "Sz", n)
    H[n] += setelt(ll[3]) * setelt(rl[2]) * op(sites, "Sz", n)
    H[n] += setelt(ll[4]) * setelt(rl[3]) * op(sites, "Sz", n)
    H[n] += setelt(ll[4]) * setelt(rl[1]) * op(sites, "Sx", n) * h[n]
  end
  H[1] *= setelt(link[1][4])
  H[N] *= setelt(link[N + 1][1])
  return H
end

function fourSiteIsingMPO(sites)::MPO
  H = MPO(sites)
  N = length(H)
  link = Vector{Index}(undef, N + 1)
  for n in 1:(N + 1)
    link[n] = Index(5, "Link,l=$(n-1)")
  end
  for n in 1:N
    s = sites[n]
    ll = link[n]
    rl = link[n + 1]
    H[n] = ITensor(ll, s, s', rl)
    H[n] += setelt(ll[1]) * setelt(rl[1]) * op(sites, "Id", n)
    H[n] += setelt(ll[5]) * setelt(rl[5]) * op(sites, "Id", n)
    H[n] += setelt(ll[2]) * setelt(rl[1]) * op(sites, "Sz", n)
    H[n] += setelt(ll[3]) * setelt(rl[2]) * op(sites, "Sz", n)
    H[n] += setelt(ll[4]) * setelt(rl[3]) * op(sites, "Sz", n)
    H[n] += setelt(ll[5]) * setelt(rl[4]) * op(sites, "Sz", n)
  end
  H[1] *= setelt(link[1][5])
  H[N] *= setelt(link[N + 1][1])
  return H
end

@testset "AutoMPO" begin
  N = 10

  @testset "Show MPOTerm" begin
    ampo = AutoMPO()
    add!(ampo, "Sz", 1, "Sz", 2)
    @test length(sprint(show, ITensors.data(ampo)[1])) > 1
  end

  @testset "Show AutoMPO" begin
    ampo = AutoMPO()
    add!(ampo, "Sz", 1, "Sz", 2)
    add!(ampo, "Sz", 2, "Sz", 3)
    @test length(sprint(show, ampo)) > 1
  end

  @testset "Single creation op" begin
    ampo = AutoMPO()
    add!(ampo, "Adagup", 3)
    sites = siteinds("Electron", N)
    W = MPO(ampo, sites)
    psi = makeRandomMPS(sites)
    cdu_psi = copy(psi)
    cdu_psi[3] = noprime(cdu_psi[3] * op(sites, "Adagup", 3))
    @test inner(psi, W, psi) ≈ inner(cdu_psi, psi)
  end

  @testset "Ising" begin
    ampo = AutoMPO()
    for j in 1:(N - 1)
      ampo += "Sz", j, "Sz", j + 1
    end
    sites = siteinds("S=1/2", N)
    Ha = MPO(ampo, sites)
    He = isingMPO(sites)
    psi = makeRandomMPS(sites)
    Oa = inner(psi, Ha, psi)
    Oe = inner(psi, He, psi)
    @test Oa ≈ Oe
  end

  @testset "Ising-Different Order" begin
    ampo = AutoMPO()
    for j in 1:(N - 1)
      ampo += "Sz", j, "Sz", j + 1
    end
    sites = siteinds("S=1/2", N)
    Ha = MPO(ampo, sites)
    He = isingMPO(sites)
    psi = makeRandomMPS(sites)
    Oa = inner(psi, Ha, psi)
    Oe = inner(psi, He, psi)
    @test Oa ≈ Oe
  end

  @testset "Heisenberg" begin
    ampo = AutoMPO()
    h = rand(N) #random magnetic fields
    for j in 1:(N - 1)
      ampo += "Sz", j, "Sz", j + 1
      ampo += 0.5, "S+", j, "S-", j + 1
      ampo += 0.5, "S-", j, "S+", j + 1
    end
    for j in 1:N
      ampo += h[j], "Sz", j
    end

    sites = siteinds("S=1/2", N)
    Ha = MPO(ampo, sites)
    He = heisenbergMPO(sites, h)
    psi = makeRandomMPS(sites)
    Oa = inner(psi, Ha, psi)
    Oe = inner(psi, He, psi)
    @test Oa ≈ Oe
  end

  @testset "Multiple Onsite Ops" begin
    sites = siteinds("S=1", N)
    ampo1 = AutoMPO()
    for j in 1:(N - 1)
      ampo1 += "Sz", j, "Sz", j + 1
      ampo1 += 0.5, "S+", j, "S-", j + 1
      ampo1 += 0.5, "S-", j, "S+", j + 1
    end
    for j in 1:N
      ampo1 += "Sz*Sz", j
    end
    Ha1 = MPO(ampo1, sites)

    ampo2 = AutoMPO()
    for j in 1:(N - 1)
      ampo2 += "Sz", j, "Sz", j + 1
      ampo2 += 0.5, "S+", j, "S-", j + 1
      ampo2 += 0.5, "S-", j, "S+", j + 1
    end
    for j in 1:N
      ampo2 += "Sz", j, "Sz", j
    end
    Ha2 = MPO(ampo2, sites)

    He = heisenbergMPO(sites, ones(N), "Sz*Sz")
    psi = makeRandomMPS(sites)
    Oe = inner(psi, He, psi)
    Oa1 = inner(psi, Ha1, psi)
    @test Oa1 ≈ Oe
    Oa2 = inner(psi, Ha2, psi)
    @test Oa2 ≈ Oe
  end

  @testset "Three-site ops" begin
    ampo = AutoMPO()
    # To test version of add! taking a coefficient
    add!(ampo, 1.0, "Sz", 1, "Sz", 2, "Sz", 3)
    @test length(ITensors.data(ampo)) == 1
    for j in 2:(N - 2)
      add!(ampo, "Sz", j, "Sz", j + 1, "Sz", j + 2)
    end
    h = ones(N)
    for j in 1:N
      add!(ampo, h[j], "Sx", j)
    end
    sites = siteinds("S=1/2", N)
    Ha = MPO(ampo, sites)
    He = threeSiteIsingMPO(sites, h)
    psi = makeRandomMPS(sites)
    Oa = inner(psi, Ha, psi)
    Oe = inner(psi, He, psi)
    @test Oa ≈ Oe
  end

  @testset "Four-site ops" begin
    ampo = AutoMPO()
    for j in 1:(N - 3)
      add!(ampo, "Sz", j, "Sz", j + 1, "Sz", j + 2, "Sz", j + 3)
    end
    sites = siteinds("S=1/2", N)
    Ha = MPO(ampo, sites)
    He = fourSiteIsingMPO(sites)
    psi = makeRandomMPS(sites)
    Oa = inner(psi, Ha, psi)
    Oe = inner(psi, He, psi)
    @test Oa ≈ Oe
  end

  @testset "Next-neighbor Heisenberg" begin
    ampo = AutoMPO()
    J1 = 1.0
    J2 = 0.5
    for j in 1:(N - 1)
      add!(ampo, J1, "Sz", j, "Sz", j + 1)
      add!(ampo, J1 * 0.5, "S+", j, "S-", j + 1)
      add!(ampo, J1 * 0.5, "S-", j, "S+", j + 1)
    end
    for j in 1:(N - 2)
      add!(ampo, J2, "Sz", j, "Sz", j + 2)
      add!(ampo, J2 * 0.5, "S+", j, "S-", j + 2)
      add!(ampo, J2 * 0.5, "S-", j, "S+", j + 2)
    end
    sites = siteinds("S=1/2", N)
    Ha = MPO(ampo, sites)

    He = NNheisenbergMPO(sites, J1, J2)
    psi = makeRandomMPS(sites)
    Oa = inner(psi, Ha, psi)
    Oe = inner(psi, He, psi)
    @test Oa ≈ Oe
    #@test maxlinkdim(Ha) == 8
  end

  @testset "Onsite Regression Test" begin
    sites = siteinds("S=1", 4)
    ampo = AutoMPO()
    add!(ampo, 0.5, "Sx", 1)
    add!(ampo, 0.5, "Sy", 1)
    H = MPO(ampo, sites)
    l = commonind(H[1], H[2])
    T = setelt(l[1]) * H[1]
    O = op(sites[1], "Sx") + op(sites[1], "Sy")
    @test norm(T - 0.5 * O) < 1E-8

    sites = siteinds("S=1", 2)
    ampo = AutoMPO()
    add!(ampo, 0.5im, "Sx", 1)
    add!(ampo, 0.5, "Sy", 1)
    H = MPO(ampo, sites)
    T = H[1] * H[2]
    O =
      im * op(sites[1], "Sx") * op(sites[2], "Id") + op(sites[1], "Sy") * op(sites[2], "Id")
    @test norm(T - 0.5 * O) < 1E-8
  end

  @testset "+ syntax" begin
    @testset "Single creation op" begin
      ampo = AutoMPO()
      ampo += "Adagup", 3
      sites = siteinds("Electron", N)
      W = MPO(ampo, sites)
      psi = makeRandomMPS(sites)
      cdu_psi = copy(psi)
      cdu_psi[3] = noprime(cdu_psi[3] * op(sites, "Adagup", 3))
      @test inner(psi, W, psi) ≈ inner(cdu_psi, psi)
    end

    @testset "Ising" begin
      ampo = AutoMPO()
      for j in 1:(N - 1)
        ampo += "Sz", j, "Sz", j + 1
      end
      sites = siteinds("S=1/2", N)
      Ha = MPO(ampo, sites)
      He = isingMPO(sites)
      psi = makeRandomMPS(sites)
      Oa = inner(psi, Ha, psi)
      Oe = inner(psi, He, psi)
      @test Oa ≈ Oe
    end

    @testset "Ising-Different Order" begin
      ampo = AutoMPO()
      for j in 1:(N - 1)
        ampo += "Sz", j + 1, "Sz", j
      end
      sites = siteinds("S=1/2", N)
      Ha = MPO(ampo, sites)
      He = isingMPO(sites)
      psi = makeRandomMPS(sites)
      Oa = inner(psi, Ha, psi)
      Oe = inner(psi, He, psi)
      @test Oa ≈ Oe
    end

    @testset "Heisenberg" begin
      ampo = AutoMPO()
      h = rand(N) #random magnetic fields
      for j in 1:(N - 1)
        ampo += "Sz", j, "Sz", j + 1
        ampo += 0.5, "S+", j, "S-", j + 1
        ampo += 0.5, "S-", j, "S+", j + 1
      end
      for j in 1:N
        ampo += h[j], "Sz", j
      end

      sites = siteinds("S=1/2", N)
      Ha = MPO(ampo, sites)
      He = heisenbergMPO(sites, h)
      psi = makeRandomMPS(sites)
      Oa = inner(psi, Ha, psi)
      Oe = inner(psi, He, psi)
      @test Oa ≈ Oe
    end

    @testset "Multiple Onsite Ops" begin
      sites = siteinds("S=1", N)
      ampo1 = AutoMPO()
      for j in 1:(N - 1)
        ampo1 += "Sz", j, "Sz", j + 1
        ampo1 += 0.5, "S+", j, "S-", j + 1
        ampo1 += 0.5, "S-", j, "S+", j + 1
      end
      for j in 1:N
        ampo1 += "Sz*Sz", j
      end
      Ha1 = MPO(ampo1, sites)

      ampo2 = AutoMPO()
      for j in 1:(N - 1)
        ampo2 += "Sz", j, "Sz", j + 1
        ampo2 += 0.5, "S+", j, "S-", j + 1
        ampo2 += 0.5, "S-", j, "S+", j + 1
      end
      for j in 1:N
        ampo2 += "Sz", j, "Sz", j
      end
      Ha2 = MPO(ampo2, sites)

      He = heisenbergMPO(sites, ones(N), "Sz*Sz")
      psi = makeRandomMPS(sites)
      Oe = inner(psi, He, psi)
      Oa1 = inner(psi, Ha1, psi)
      @test Oa1 ≈ Oe
      Oa2 = inner(psi, Ha2, psi)
      @test Oa2 ≈ Oe
    end

    @testset "Three-site ops" begin
      ampo = AutoMPO()
      # To test version of add! taking a coefficient
      ampo += 1.0, "Sz", 1, "Sz", 2, "Sz", 3
      @test length(ITensors.data(ampo)) == 1
      for j in 2:(N - 2)
        ampo += "Sz", j, "Sz", j + 1, "Sz", j + 2
      end
      h = ones(N)
      for j in 1:N
        ampo += h[j], "Sx", j
      end
      sites = siteinds("S=1/2", N)
      Ha = MPO(ampo, sites)
      He = threeSiteIsingMPO(sites, h)
      psi = makeRandomMPS(sites)
      Oa = inner(psi, Ha, psi)
      Oe = inner(psi, He, psi)
      @test Oa ≈ Oe
    end

    @testset "Four-site ops" begin
      ampo = AutoMPO()
      for j in 1:(N - 3)
        ampo += "Sz", j, "Sz", j + 1, "Sz", j + 2, "Sz", j + 3
      end
      sites = siteinds("S=1/2", N)
      Ha = MPO(ampo, sites)
      He = fourSiteIsingMPO(sites)
      psi = makeRandomMPS(sites)
      Oa = inner(psi, Ha, psi)
      Oe = inner(psi, He, psi)
      @test Oa ≈ Oe
    end

    @testset "Next-neighbor Heisenberg" begin
      ampo = AutoMPO()
      J1 = 1.0
      J2 = 0.5
      for j in 1:(N - 1)
        ampo += J1, "Sz", j, "Sz", j + 1
        ampo += J1 * 0.5, "S+", j, "S-", j + 1
        ampo += J1 * 0.5, "S-", j, "S+", j + 1
      end
      for j in 1:(N - 2)
        ampo += J2, "Sz", j, "Sz", j + 2
        ampo += J2 * 0.5, "S+", j, "S-", j + 2
        ampo += J2 * 0.5, "S-", j, "S+", j + 2
      end
      sites = siteinds("S=1/2", N)
      Ha = MPO(ampo, sites)

      He = NNheisenbergMPO(sites, J1, J2)
      psi = makeRandomMPS(sites)
      Oa = inner(psi, Ha, psi)
      Oe = inner(psi, He, psi)
      @test Oa ≈ Oe
      #@test maxlinkdim(Ha) == 8
    end

    #@testset "-= syntax" begin
    #  ampo = AutoMPO()
    #  ampo += (-1,"Sz",1,"Sz",2)
    #  ampo2 = AutoMPO()
    #  ampo2 -= ("Sz",1,"Sz",2)
    #  @test ampo == ampo2
    #end

    @testset "Onsite Regression Test" begin
      sites = siteinds("S=1", 4)
      ampo = AutoMPO()
      ampo += 0.5, "Sx", 1
      ampo += 0.5, "Sy", 1
      H = MPO(ampo, sites)
      l = commonind(H[1], H[2])
      T = setelt(l[1]) * H[1]
      O = op(sites[1], "Sx") + op(sites[1], "Sy")
      @test norm(T - 0.5 * O) < 1E-8

      sites = siteinds("S=1", 2)
      ampo = AutoMPO()
      ampo += 0.5im, "Sx", 1
      ampo += 0.5, "Sy", 1
      H = MPO(ampo, sites)
      T = H[1] * H[2]
      O =
        im * op(sites[1], "Sx") * op(sites[2], "Id") +
        op(sites[1], "Sy") * op(sites[2], "Id")
      @test norm(T - 0.5 * O) < 1E-8
    end
  end

  @testset ".+= and .-= syntax" begin

    #@testset ".-= syntax" begin
    #  ampo = AutoMPO()
    #  ampo .+= (-1,"Sz",1,"Sz",2)
    #  ampo2 = AutoMPO()
    #  ampo2 .-= ("Sz",1,"Sz",2)
    #  @test ampo == ampo2
    #end

    @testset "Single creation op" begin
      ampo = AutoMPO()
      ampo .+= "Adagup", 3
      sites = siteinds("Electron", N)
      W = MPO(ampo, sites)
      psi = makeRandomMPS(sites)
      cdu_psi = copy(psi)
      cdu_psi[3] = noprime(cdu_psi[3] * op(sites, "Adagup", 3))
      @test inner(psi, W, psi) ≈ inner(cdu_psi, psi)
    end

    @testset "Ising" begin
      ampo = AutoMPO()
      for j in 1:(N - 1)
        ampo .+= "Sz", j, "Sz", j + 1
      end
      sites = siteinds("S=1/2", N)
      Ha = MPO(ampo, sites)
      He = isingMPO(sites)
      psi = makeRandomMPS(sites)
      Oa = inner(psi, Ha, psi)
      Oe = inner(psi, He, psi)
      @test Oa ≈ Oe
    end

    @testset "Ising-Different Order" begin
      ampo = AutoMPO()
      for j in 1:(N - 1)
        ampo .+= "Sz", j + 1, "Sz", j
      end
      sites = siteinds("S=1/2", N)
      Ha = MPO(ampo, sites)
      He = isingMPO(sites)
      psi = makeRandomMPS(sites)
      Oa = inner(psi, Ha, psi)
      Oe = inner(psi, He, psi)
      @test Oa ≈ Oe
    end

    @testset "Heisenberg" begin
      ampo = AutoMPO()
      h = rand(N) #random magnetic fields
      for j in 1:(N - 1)
        ampo .+= "Sz", j, "Sz", j + 1
        ampo .+= 0.5, "S+", j, "S-", j + 1
        ampo .+= 0.5, "S-", j, "S+", j + 1
      end
      for j in 1:N
        ampo .+= h[j], "Sz", j
      end

      sites = siteinds("S=1/2", N)
      Ha = MPO(ampo, sites)
      He = heisenbergMPO(sites, h)
      psi = makeRandomMPS(sites)
      Oa = inner(psi, Ha, psi)
      Oe = inner(psi, He, psi)
      @test Oa ≈ Oe
    end

    @testset "Multiple Onsite Ops" begin
      sites = siteinds("S=1", N)
      ampo1 = AutoMPO()
      for j in 1:(N - 1)
        ampo1 .+= "Sz", j, "Sz", j + 1
        ampo1 .+= 0.5, "S+", j, "S-", j + 1
        ampo1 .+= 0.5, "S-", j, "S+", j + 1
      end
      for j in 1:N
        ampo1 .+= "Sz*Sz", j
      end
      Ha1 = MPO(ampo1, sites)

      ampo2 = AutoMPO()
      for j in 1:(N - 1)
        ampo2 .+= "Sz", j, "Sz", j + 1
        ampo2 .+= 0.5, "S+", j, "S-", j + 1
        ampo2 .+= 0.5, "S-", j, "S+", j + 1
      end
      for j in 1:N
        ampo2 .+= "Sz", j, "Sz", j
      end
      Ha2 = MPO(ampo2, sites)

      He = heisenbergMPO(sites, ones(N), "Sz*Sz")
      psi = makeRandomMPS(sites)
      Oe = inner(psi, He, psi)
      Oa1 = inner(psi, Ha1, psi)
      @test Oa1 ≈ Oe
      Oa2 = inner(psi, Ha2, psi)
      @test Oa2 ≈ Oe
    end

    @testset "Three-site ops" begin
      ampo = AutoMPO()
      # To test version of add! taking a coefficient
      ampo .+= 1.0, "Sz", 1, "Sz", 2, "Sz", 3
      @test length(ITensors.data(ampo)) == 1
      for j in 2:(N - 2)
        ampo .+= "Sz", j, "Sz", j + 1, "Sz", j + 2
      end
      h = ones(N)
      for j in 1:N
        ampo .+= h[j], "Sx", j
      end
      sites = siteinds("S=1/2", N)
      Ha = MPO(ampo, sites)
      He = threeSiteIsingMPO(sites, h)
      psi = makeRandomMPS(sites)
      Oa = inner(psi, Ha, psi)
      Oe = inner(psi, He, psi)
      @test Oa ≈ Oe
    end

    @testset "Four-site ops" begin
      ampo = AutoMPO()
      for j in 1:(N - 3)
        ampo .+= "Sz", j, "Sz", j + 1, "Sz", j + 2, "Sz", j + 3
      end
      sites = siteinds("S=1/2", N)
      Ha = MPO(ampo, sites)
      He = fourSiteIsingMPO(sites)
      psi = makeRandomMPS(sites)
      Oa = inner(psi, Ha, psi)
      Oe = inner(psi, He, psi)
      @test Oa ≈ Oe
    end

    @testset "Next-neighbor Heisenberg" begin
      ampo = AutoMPO()
      J1 = 1.0
      J2 = 0.5
      for j in 1:(N - 1)
        ampo .+= J1, "Sz", j, "Sz", j + 1
        ampo .+= J1 * 0.5, "S+", j, "S-", j + 1
        ampo .+= J1 * 0.5, "S-", j, "S+", j + 1
      end
      for j in 1:(N - 2)
        ampo .+= J2, "Sz", j, "Sz", j + 2
        ampo .+= J2 * 0.5, "S+", j, "S-", j + 2
        ampo .+= J2 * 0.5, "S-", j, "S+", j + 2
      end
      sites = siteinds("S=1/2", N)
      Ha = MPO(ampo, sites)

      He = NNheisenbergMPO(sites, J1, J2)
      psi = makeRandomMPS(sites)
      Oa = inner(psi, Ha, psi)
      Oe = inner(psi, He, psi)
      @test Oa ≈ Oe
      #@test maxlinkdim(Ha) == 8
    end

    @testset "Onsite Regression Test" begin
      sites = siteinds("S=1", 4)
      ampo = AutoMPO()
      ampo .+= 0.5, "Sx", 1
      ampo .+= 0.5, "Sy", 1
      H = MPO(ampo, sites)
      l = commonind(H[1], H[2])
      T = setelt(l[1]) * H[1]
      O = op(sites[1], "Sx") + op(sites[1], "Sy")
      @test norm(T - 0.5 * O) < 1E-8

      sites = siteinds("S=1", 2)
      ampo = AutoMPO()
      ampo .+= 0.5im, "Sx", 1
      ampo .+= 0.5, "Sy", 1
      H = MPO(ampo, sites)
      T = H[1] * H[2]
      O =
        im * op(sites[1], "Sx") * op(sites[2], "Id") +
        op(sites[1], "Sy") * op(sites[2], "Id")
      @test norm(T - 0.5 * O) < 1E-8
    end
  end

  @testset "Fermionic Operators" begin
    N = 5
    s = siteinds("Fermion", N)

    a1 = AutoMPO()
    a1 += "Cdag", 1, "C", 3
    M1 = MPO(a1, s)

    a2 = AutoMPO()
    a2 += -1, "C", 3, "Cdag", 1
    M2 = MPO(a2, s)

    a3 = AutoMPO()
    a3 += "Cdag", 1, "N", 2, "C", 3
    M3 = MPO(a3, s)

    p011 = productMPS(s, [1, 2, 2, 1, 1])
    p110 = productMPS(s, [2, 2, 1, 1, 1])

    @test inner(p110, M1, p011) ≈ -1.0
    @test inner(p110, M2, p011) ≈ -1.0
    @test inner(p110, M3, p011) ≈ -1.0

    p001 = productMPS(s, [1, 1, 2, 1, 1])
    p100 = productMPS(s, [2, 1, 1, 1, 1])

    @test inner(p100, M1, p001) ≈ +1.0
    @test inner(p100, M2, p001) ≈ +1.0
    @test inner(p100, M3, p001) ≈ 0.0

    #
    # Repeat similar test but
    # with Electron sites
    #

    s = siteinds("Electron", N; conserve_qns=true)

    a1 = AutoMPO()
    a1 += "Cdagup", 1, "Cup", 3
    M1 = MPO(a1, s)

    a2 = AutoMPO()
    a2 += -1, "Cdn", 3, "Cdagdn", 1
    M2 = MPO(a2, s)

    p0uu = productMPS(s, [1, 2, 2, 1, 1])
    puu0 = productMPS(s, [2, 2, 1, 1, 1])
    p0ud = productMPS(s, [1, 2, 3, 1, 1])
    pdu0 = productMPS(s, [3, 2, 1, 1, 1])
    p00u = productMPS(s, [1, 1, 2, 1, 1])
    pu00 = productMPS(s, [2, 1, 1, 1, 1])
    p00d = productMPS(s, [1, 1, 3, 1, 1])
    pd00 = productMPS(s, [3, 1, 1, 1, 1])

    @test inner(puu0, M1, p0uu) ≈ -1.0
    @test inner(pdu0, M2, p0ud) ≈ -1.0
    @test inner(pu00, M1, p00u) ≈ +1.0
    @test inner(pd00, M2, p00d) ≈ +1.0
  end

  @testset "Complex AutoMPO Coefs" begin
    N = 4

    for use_qn in [false, true]
      sites = siteinds("S=1/2", N; conserve_qns=use_qn)
      ampo = AutoMPO()
      for i in 1:(N - 1)
        ampo += +1im, "S+", i, "S-", i + 1
        ampo += -1im, "S-", i, "S+", i + 1
      end
      H = MPO(ampo, sites)
      psiud = productMPS(sites, [1, 2, 1, 2])
      psidu = productMPS(sites, [2, 1, 1, 2])
      @test inner(psiud, H, psidu) ≈ +1im
      @test inner(psidu, H, psiud) ≈ -1im
    end
  end

  @testset "Fermion AutoMPO Issue 514 Regression Test" begin
    N = 4
    s = siteinds("Electron", N; conserve_qns=true)
    ampo1 = AutoMPO()
    ampo2 = AutoMPO()

    ampo1 += "Nup", 1
    ampo2 += "Cdagup", 1, "Cup", 1

    M1 = MPO(ampo1, s)
    M2 = MPO(ampo2, s)

    H1 = M1[1] * M1[2] * M1[3] * M1[4]
    H2 = M2[1] * M2[2] * M2[3] * M2[4]

    @test norm(H1 - H2) ≈ 0.0
  end
end

nothing


text:
using HyperGraphs

## hyperedges and hypergraphs ##

# hyperedges: constructors
HyperEdge()
HyperEdge{Int8}()
HyperEdge([:a, :b])

# hyperedges: functions
vertices_e = [1, 2, 3]
e = HyperEdge(vertices_e)
@test vertices(e) == vertices_e
@test vertices(HyperEdge([1, :∅])) == [1, :∅]
vertices_es = [[1, 2, 3], [2, 3, 4]]
es = [HyperEdge(vertices_es[1]), HyperEdge(vertices_es[2])]
@test vertices(es) == vertices_es
@test eltype(e) == Int64
@test HyperEdge([1, 2, 3]) == e

# hypergraphs: constructors
HyperGraph()
HyperGraph{Int8}()
HyperGraph(HyperEdge([1, :∅]))
HyperGraph([HyperEdge([1, 1]), HyperEdge([2, 2])])
HyperGraph(HyperEdge([1, 1]))
@test HyperGraph([1], HyperEdge{Int}[]) == HyperGraph{Int}([1], HyperEdge[])
@test HyperGraph(1, HyperEdge([1, 1])) == HyperGraph([1], HyperEdge([1, 1])) == HyperGraph(1, [HyperEdge([1, 1])]) == HyperGraph([1], [HyperEdge([1, 1])])
@test HyperGraph(es) == HyperGraph(unique(vcat(vertices_es...)), es)

# hypergraphs: functions
@test vertices(HyperGraph(HyperEdge([1, :∅]))) == [1, :∅]
x = HyperGraph(es)
@test vertices(x) == unique(vcat(vertices_es...)) == [1, 2, 3, 4]
@test vertices(hyperedges(x)) == vertices_es
@test hyperedges(x) == es
@test eltype(x) == eltype(hyperedges(x)[1])

## chemical hyperedges and chemical hypergraphs ##

# incidence: constructors
SpeciesSet{Symbol}()
SpeciesSet()
@test SpeciesSet([1]) == SpeciesSet(1)
@test SpeciesSet(["X"]) == SpeciesSet("X") == SpeciesSet("X", 1) == SpeciesSet(["X"], [1])

# incidence: functions
s = SpeciesSet(:∅)
@test objects(s) == [:∅]
@test multiplicities(s) == [1]
@test multiplicities(SpeciesSet([1, 2, 3], [1, 1, 2])) == [1, 1, 2]
@test_throws MethodError SpeciesSet(["X"], [1.2]) # stoichiometries are defined on integers
@test species(s) == [:∅]
@test stoich(s) == [1]
@test length(s) == 1

# chemical hyperedges: constructors
ChemicalHyperEdge{String}()
ChemicalHyperEdge()
ChemicalHyperEdge(["X", "Y"], ["Z"], 2)
ChemicalHyperEdge(["X", "Y"], ["Z"], [[[2]]]) # is this undesirable?
ChemicalHyperEdge(["X"], ["Y"])
@test ChemicalHyperEdge(["X"], ["Y"], 1) == ChemicalHyperEdge(["X"], ["Y"]) # rate defaults to 1

# chemical hyperedges: functions
reactant = SpeciesSet(["X"])
product = SpeciesSet(["Y"], [2])
che = ChemicalHyperEdge(reactant, product)
@test src(che) == objects(reactant)
@test tgt(che) == objects(product)
@test vertices(che) == ["X", "Y"]
@test src(ChemicalHyperEdge(["X", "Y"], ["Z"])) == ["X", "Y"]
@test src_multiplicities(che) == [1]
@test tgt_multiplicities(che) == [2]
@test weight(ChemicalHyperEdge()) == 1
e = ChemicalHyperEdge([], [], 6)
@test weight(e) == 6
x = ChemicalHyperGraph([e, e])
@test weights(hyperedges(x)) == weights(x) == [6, 6]
srcs = [[1], [4]]; tgts = [[1, 2, 3], [2, 3, 4]]
ches = [ChemicalHyperEdge(src, tgt) for (src, tgt) in zip(srcs, tgts)]
@test src(ches) == srcs
@test tgt(ches) == tgts
@test vertices(ches) == [[1, 1, 2, 3], [4, 2, 3, 4]]
@test eltype(reactant) == eltype(che) == String
@test src_stoich(che) == inputs_stoich(che) == [1]
@test tgt_stoich(che) == outputs_stoich(che) == [2]

# chemical hypergraphs: constructors
ChemicalHyperGraph{String}()
ChemicalHyperGraph()

# chemical hypergraphs: functions
vertices_chx = [1, 2, 3, 4]
chx = ChemicalHyperGraph(vertices_chx, ches)
@test vertices(chx) == vertices_chx
hyperedges(chx)
@test ChemicalHyperGraph(ches) == chx
@test eltype(chx) == eltype(hyperedges(chx)[1])


text:
# This file was generated, do not modify it. # hide
figure(figsize=(8,6))

res = r.plotting

vals_sf = res.parameter_values[:, 1]
vals_bf = res.parameter_values[:, 2]

tricontourf(vals_sf, vals_bf, res.measurements)
xlabel("Number of sub-features", fontsize=14)
ylabel("Bagging fraction", fontsize=14)
xticks([1, 2, 3], fontsize=12)
yticks(fontsize=12)

savefig(joinpath(@OUTPUT, "A-model-tuning-hm.svg")) # hide


text:
function calculate_set(U::TaylorN, V::TaylorN, bounds)

    xs = []
    ys = []

    # iterate over 4 sides of initial box

    num_points = 30

    b = bounds[2].lo
    for a in linspace(bounds[1].lo, bounds[2].hi, num_points)
        push!(xs, U([a, b]))
        push!(ys, V([a, b]))
    end

    a = bounds[1].hi
    for b in linspace(bounds[2].lo, bounds[2].hi, num_points)
        push!(xs, U([a, b]))
        push!(ys, V([a, b]))
    end

    b = bounds[2].hi
    for a in linspace(bounds[2].hi, bounds[2].lo, num_points)
        push!(xs, U([a, b]))
        push!(ys, V([a, b]))
    end

    a = bounds[1].lo
    for b in linspace(bounds[2].hi, bounds[2].lo, num_points)
        push!(xs, U([a, b]))
        push!(ys, V([a, b]))
    end

    return mid.(xs), mid.(ys)

end