text:
# The most important type of the package, GenericMERA, on which all specific MERA
# implementations (binary, ternary, ...) are built. Methods and functions that can be
# implemented on this level of abstraction, without having to know the details of the
# specific MERA type.
# To be included in MERAKit.jl.

"""
    GenericMERA{N, LT <: Layer, OT}

A `GenericMERA` is a collection of `Layer`s. The type of these layers then determines
whether the MERA is binary, ternary, etc.

On conventions and terminology:
* The physical indices of the MERA are at the "bottom", the scale invariant part at the
  "top".
* The counting of layers starts from the bottom, so the layer with physical indices is
  layer 1. The last layer is the scale invariant one, that then repeats upwards to infinity.
* Each layer is thought of as a linear map from its top, or input space to its bottom, or
  output space.

The type parameters are
* `N`:  The number of distinct layers (N-1 transition layers and one scale invariant one).
* `LT`: Layer type.
* `OT`: Operator type. The type of ascended and descended operators for this MERA. Determined
      from `LT`. Typically a `TensorMap` with input and output indices matching the causal
      cone width.

`GenericMERA` is immutable, and the layers can not be changed after construction. All
functions that modify a `GenericMERA` return new objects.
"""
struct GenericMERA{N, LT <: Layer, OT}
    # LT stands for Layer Type, OT for Operator Type.
    layers::NTuple{N, LT}
    cache::MERACache{N, LT, OT}

    function GenericMERA{N, LT, OT}(
        layers::NTuple{N}, cache::MERACache{N, LT, OT}
    ) where {N, LT, OT}
        # Note that this prevents the creation of types like GenericMERA{3, SimpleLayer}:
        # The second type parameter must be exactly the element type of layers, specified at
        # as a concrete type. This is intentional, to avoid accidentally creating
        # unnecessarily abstract types that would hamper inference.
        @assert LT === eltype(typeof(layers))
        @assert isconcretetype(LT)
        @assert OT === operatortype(LT)
        return new{N, LT, OT}(layers, cache)
    end
end

# # # Constructors

function GenericMERA(layers::NTuple{N}, cache::MERACache{N}) where {N}
    LT = eltype(typeof(layers))
    OT = operatortype(LT)
    cache::MERACache{N, LT, OT}
    return GenericMERA{N, LT, OT}(layers, cache)
end

function GenericMERA(layers::NTuple{N}) where {N}
    LT = eltype(typeof(layers))
    OT = operatortype(LT)
    cache = MERACache{N, LT, OT}()
    return GenericMERA(layers, cache)
end

GenericMERA(layers) = GenericMERA(tuple(layers...))

function GenericMERA{N, LT, OT}(layers::NTuple{N}) where {N, LT, OT}
    cache = MERACache{N, LT, OT}()
    return GenericMERA{N, LT, OT}(layers, cache)
end

function GenericMERA{N, LT, OT}(layers) where {N, LT, OT}
    return GenericMERA{N, LT, OT}(tuple(layers...))
end

function (::Type{GenericMERA{M, LT, OT} where M})(layers::NTuple{N}) where {N, LT, OT}
    return GenericMERA(layers)
end

function (::Type{GenericMERA{M, LT, OT} where M})(layers) where {LT, OT}
    return (GenericMERA where {N})(tuple(layers...))
end

# # # Basic utility functions

"""
    layertype(m::GenericMERA)
    layertype(::Type{<:GenericMERA})

Return the type of the layers of the MERA.
"""
layertype(::Type{GenericMERA{N, LT, OT} where N}) where {LT, OT} = LT
layertype(::Type{GenericMERA{N, LT, OT}}) where {N, LT, OT} = LT

"""
    baselayertype(m::GenericMERA)
    baselayertype(::Type{<:GenericMERA})

Return the generic type of the layers of the MERA, without specific type parameters
"""
baselayertype(::Type{GenericMERA{N, LT, OT} where N}) where {LT, OT} = baselayertype(LT)
baselayertype(::Type{GenericMERA{N, LT, OT}}) where {N, LT, OT} = baselayertype(LT)

"""
    operatortype(m::GenericMERA)
    operatortype(::Type{<: GenericMERA})

Return the type of operator associated with this MERA or MERA type. This means the type of
operator that fits in the causal cone, and naturally emerges as one ascends local operators.
"""
operatortype(::Type{GenericMERA{N, LT, OT} where N}) where {LT, OT} = OT
operatortype(::Type{GenericMERA{N, LT, OT}}) where {N, LT, OT} = OT
# we could also just do:
# operatortype(M::Type{<:GenericMERA}) = operatortype(layertype(M))

"""
    scalefactor(::Type{<: GenericMERA})

The ratio by which the number of sites changes when one descends by one layer, e.g. 2 for
binary MERA, 3 for ternary.
"""
scalefactor(M::Type{<:GenericMERA}) = scalefactor(layertype(M))

"""
    causal_cone_width(::Type{<: GenericMERA})

Return the width of the stable causal cone for this MERA type.
"""
causal_cone_width(M::Type{<:GenericMERA}) = causal_cone_width(layertype(M))

Base.eltype(M::Type{<:GenericMERA}) = eltype(layertype(M))

"""
    num_translayers(m::GenericMERA)

Return the number of transition layers, i.e. layers below the scale invariant one, in the
MERA.
"""
num_translayers(::Type{<:GenericMERA{N}}) where {N} = N-1
num_translayers(m::GenericMERA) = num_translayers(typeof(m))

"""
    getlayer(m::GenericMERA, depth)

Return the layer at the given depth. 1 is the lowest layer, i.e. the one with physical
indices.
"""
function getlayer(m::GenericMERA, depth)
    return (depth > num_translayers(m) ?  m.layers[end] : m.layers[depth])
end

"""
    replace_layer(m::GenericMERA, layer, depth; check_invar = true)

Replace the layer at `depth` in `m` with `layer`. If `check_invar = true`, check that the
indices match afterwards.
"""
function replace_layer(
    m::GenericMERA{N, LT}, newlayer::LT, depth; check_invar = true
) where {N, LT}
    index = min(num_translayers(m)+1, depth)
    new_layers = Base.setindex(m.layers, newlayer, index)
    new_cache = replace_layer(m.cache, depth)
    new_m = GenericMERA(new_layers, new_cache)
    check_invar && space_invar(new_m)
    return new_m
end

"""
    release_transitionlayer(m::GenericMERA)

Add one more transition layer at the top of the MERA, by taking the lowest of the scale
invariant ones and releasing it to vary independently.
"""
function release_transitionlayer(m::GenericMERA)
    new_layers = (m.layers..., m.layers[end])
    new_cache = release_transitionlayer(m.cache)
    new_m = GenericMERA(new_layers, new_cache)
    return new_m
end

"""
    projectisometric(m::GenericMERA)

Project all the tensors of the MERA to respect the isometricity condition.
"""
function TensorKitManifolds.projectisometric(m::GenericMERA)
    return typeof(m)(map(projectisometric, m.layers))
end

"""
    projectisometric!(m::GenericMERA)

Project all the tensors of the MERA to respect the isometricity condition, in place.
"""
function TensorKitManifolds.projectisometric!(m::GenericMERA)
    return typeof(m)(map(projectisometric!, m.layers))
end

"""
    outputspace(m::GenericMERA, depth)

Given a MERA and a depth, return the vector space of the downwards-pointing (towards the
physical level) indices of the layer at that depth.

See also: [`inputspace`](@ref)
"""
outputspace(m::GenericMERA, depth) = outputspace(getlayer(m, depth))

"""
    inputspace(m::GenericMERA, depth)

Given a MERA and a depth, return the vector space of the upwards-pointing (towards scale
invariance) indices of the layer at that depth.

See also: [`outputspace`](@ref)
"""
inputspace(m::GenericMERA, depth) = inputspace(getlayer(m, depth))

"""
    densitymatrix_entropies(m::GenericMERA)

Return a vector of entropies for the density matrices in the MERA. The first one is for the
density matrix at the physical level, the last one is the scale invariant density matrix.
"""
function densitymatrix_entropies(m::GenericMERA)
    return [densitymatrix_entropy(x) for x in densitymatrices(m)]
end

"""
    reset_storage(m::GenericMERA)

Reset cached operators, so that they will be recomputed when they are needed.
"""
reset_storage(m::GenericMERA) = typeof(m)(m.layers, typeof(m.cache)())

# # # Generating random MERAs

"""
    random_MERA(::Type{T <: GenericMERA}, ::Type{ET}, Vouts, Vints = Vouts; kwargs...)

Generate a random MERA of type `T`.

`ET` is the element type of the MERA, e.g. `Float64`. `Vouts` are the vector spaces between
the various layers. Number of layers will be the length of `Vouts`. `Vouts[1]` will be the
physical index space, and `Vs[end]` will be the one at the scale invariant layer. `Vints`
can be a vector/tuple of the same length as `Vouts`, and include vector spaces for the
internal indices of each layer. Alternatively, it can hold any other extra parameters
specific to each layer, as it is simply passed on to the function `randomlayer`. Also passed
to `randomlayer` will be any additional keyword arguments, but these will all be the same
for each layer.

See also: [`randomlayer`](@ref)
"""
function random_MERA(::Type{T}, ::Type{ET}, Vouts, Vints = Vouts; kwargs...) where {
    T <: GenericMERA, ET
}
    L = layertype(T)
    Vins = tuple(Vouts[2:end]..., Vouts[end])
    layers = ntuple(length(Vouts)) do i
        randomlayer(L, ET, Vins[i], Vouts[i], Vints[i]; kwargs...)
    end
    m = GenericMERA(layers)
    return m
end

# TODO Should the numbering be changed, so that the bond at `depth` would be the
# output bond of layer `depth`, instead of input? Would maybe be more consistent.
"""
    expand_bonddim(m::GenericMERA, depth, V; check_invar = true)

Expand the bond dimension of the MERA at the given depth.

The indices to expand are the input indices of the layer at `depth`, i.e. `depth = 1` means
the lowest virtual indices. The new bond dimension is given by `V`, which can be one of
three different things:
* `V` can be the new vector space of the indices in question.
* For a non-symmetric MERA `V` can be just a number, the bond dimension.
* For a symmetric MERA `V` can be dictionary of `irrep => block dimension`. Not all irreps
  for a bond need to be listed, the ones left out are left untouched.

The expansion is done by padding tensors with zeros. Note that this breaks isometricity of
the individual tensors. This is however of no consequence, since the MERA as a state remains
exactly the same. A round of optimization on the MERA will restore isometricity of each
tensor, or `projectisometric` can be called to do so explicitly.

If `check_invar = true` the function checks that the bond dimensions of various layers match
after the expansion.
"""
function expand_bonddim(m::GenericMERA, depth, V; check_invar = true)
    V_new = expand_vectorspace(inputspace(m, depth), V)
    return expand_bonddim(m, depth, V_new; check_invar=check_invar)
end

function expand_bonddim(m::GenericMERA, depth, V::IndexSpace; check_invar = true)
    if depth > num_translayers(m)
        msg = "expand_bonddim called with too large depth. To change the scale invariant bond dimension, use depth = num_translayers(m)."
        throw(ArgumentError(msg))
    end
    l = getlayer(m, depth)
    l = expand_inputspace(l, V)
    next_l = getlayer(m, depth+1)
    next_l = expand_outputspace(next_l, V)
    if depth == num_translayers(m)
        # next_layer is the scale invariant part, so we need to change its top index too
        # since we changed the bottom.
        next_l = expand_inputspace(next_l, V)
    end
    m = replace_layer(m, l, depth; check_invar = false)
    m = replace_layer(m, next_l, depth+1; check_invar = check_invar)
    expand_bonddim!(m.cache, depth, V)
    return m
end

"""
    expand_internal_bonddim(m::GenericMERA, depth, V; check_invar = true)

Expand the bond dimension of the layer-internal indices of the MERA at the given depth.

The new bond dimension is given by `V`, which can be one of
three different things:
* `V` can be the new vector space of the indices in question.
* For a non-symmetric MERA `V` can be just a number, the bond dimension.
* For a symmetric MERA `V` can be dictionary of `irrep => block dimension`. Not all irreps
  for a bond need to be listed, the ones left out are left untouched.

The expansion is done by padding tensors with zeros. Note that this breaks isometricity of
the individual tensors. This is however of no consequence, since the MERA as a state remains
exactly the same. A round of optimization on the MERA will restore isometricity of each
tensor, or `projectisometric` can be called to do so explicitly.

Note that not all MERAs have an internal bond dimension, and some may have several, so this
function will not make sense for all MERA types. Implementation relies on
the function `expand_internalspace`, defined for each `Layer` type.
"""
function expand_internal_bonddim(m::GenericMERA, depth, V; check_invar = true)
    V_new = expand_vectorspace(internalspace(m, depth), V)
    return expand_internal_bonddim(m, depth, V_new; check_invar=check_invar)
end

function expand_internal_bonddim(m::GenericMERA, depth, V::IndexSpace; check_invar = true)
    l = getlayer(m, depth)
    l = expand_internalspace(l, V)
    m = replace_layer(m, l, depth; check_invar = check_invar)
    return m
end

"""
    remove_symmetry(m::GenericMERA)

Given a MERA which may possibly be built of symmetry preserving `TensorMap`s, return
another, equivalent MERA that has the symmetry structure stripped from it, and all tensors
are dense.
"""
remove_symmetry(m::GenericMERA) = GenericMERA(map(remove_symmetry, m.layers))

# # # Pseudo(de)serialization
# "Pseudo(de)serialization" refers to breaking the MERA down into types in Julia Base, and
# constructing it back. This can be used for storing MERAs on disk.
#
# Note that pseudoserialization discards the cache, forcing recomputation.

"""
    pseudoserialize(x)

Return a tuple of objects that can be used to reconstruct `x`, and that are all of Julia
base types.

The name refers to how this isn't quite serialization, since it doesn't break objects down
to bit strings, but kinda serves the same purpose: pseudoserialised objects can easily be
written and read from e.g. disk, without having to worry about quirks of the type system
(like with JLD).

See also: [`depseudoserialize`](@ref)
"""
function pseudoserialize(m::GenericMERA)
    return (repr(typeof(m)), map(pseudoserialize, m.layers))
end

"""
    depseudoserialize(::Type{T}, args) where T <: GenericMERA

Reconstruct an object given the output of `pseudoserialize`:
`x -> depseudoserialize(pseudoserialize(x)...)` should be an effective noop.

See also: [`pseudoserialize`](@ref)
"""
function depseudoserialize(::Type{T}, args) where T <: GenericMERA
    return GenericMERA(map(d->depseudoserialize(d...), args))
end

# # # Invariants

"""
    space_invar(m::GenericMERA)

Check that the indices of the various tensors in `m` are compatible with each other. If not,
throw an `ArgumentError`. If yes, return `true`.

This relies on two checks, `space_invar_intralayer` for checking indices within a layer and
`space_invar_interlayer` for checking indices between layers. These should be implemented
for each subtype of `Layer`.
"""
function space_invar(m::GenericMERA)
    l = getlayer(m, 1)
    # We go to num_translayers(m)+2, to check that the scale invariant layer is consistent
    # with itself.
    for i in 2:(num_translayers(m)+2)
        next_l = getlayer(m, i)
        if applicable(space_invar_intralayer, l)
            if !space_invar_intralayer(l)
                errmsg = "Mismatching bonds in MERA within layer $(i-1)."
                throw(ArgumentError(errmsg))
            end
        else
            msg = "space_invar_intralayer has no method for type $(baselayertype(m)). Please consider writing one, to enable checking for space mismatches when assigning tensors."
            @warn(msg)
        end

        if applicable(space_invar_interlayer, l, next_l)
            if !space_invar_interlayer(l, next_l)
                errmsg = "Mismatching bonds in MERA between layers $(i-1) and $i."
                throw(ArgumentError(errmsg))
            end
        else
            msg = "space_invar_interlayer has no method for type $(baselayertype(m)). Please consider writing one, to enable checking for space mismatches when assigning tensors."
            @warn(msg)
        end
        l = next_l
    end
    return true
end

# # # Scaling superoperators

"""
    ascend(op, m::GenericMERA, endscale = num_translayers(m)+1, startscale = 1)

Ascend the operator `op`, that lives on the layer `startscale`, through the MERA `m` to the
layer `endscale`. Living "on" a layer means living on the indices right below it, so
`startscale = 1` refers to the physical indices, and `endscale = num_translayers(m)+1` to
the indices just below the first scale invariant layer.
"""
function ascend(op, m::GenericMERA, endscale = num_translayers(m)+1, startscale = 1)
    if endscale < startscale
        throw(ArgumentError("endscale < startscale"))
    elseif endscale > startscale
        op_pre = ascend(op, m, endscale-1, startscale)
        l = getlayer(m, endscale-1)
        op_asc = ascend(op_pre, l)
    else
        # op_asc = convert(operatortype(m), op)
        # what if it is a charged operator, or one living on a smaller width?
        op_asc = op
    end
    return op_asc
end

"""
    descend(op, m::GenericMERA, endscale = 1, startscale = num_translayers(m)+1)

Descend the operator `op`, that lives on the layer `startscale`, through the MERA `m` to the
layer `endscale`. Living "on" a layer means living on the indices right below it, so
`endscale = 1` refers to the physical indices, and `startscale = num_translayers(m)+1` to
the indices just below the first scale invariant layer.
"""
function descend(op, m::GenericMERA, endscale = 1, startscale = num_translayers(m)+1)
    if endscale > startscale
        throw(ArgumentError("endscale > startscale"))
    elseif endscale < startscale
        op_pre = descend(op, m, endscale+1, startscale)
        l = getlayer(m, endscale)
        op_desc = descend(op_pre, l)
    else
        op_desc = convert(operatortype(m), op)
    end
    return op_desc
end

"""
    fixedpoint_densitymatrix(m::GenericMERA, pars::NamedTuple = (;))

Find the fixed point density matrix of the scale invariant part of the MERA.

To find the fixed point, we use an iterative Krylov solver. The options for the solver
should be in `pars.scaleinvariant_krylovoptions`, they will be passed to
`KrylovKit.eigsolve`.
"""
function fixedpoint_densitymatrix(m::GenericMERA, pars = (;))
    f(x) = descend(x, getlayer(m, num_translayers(m)+1))
    # If we have stored the previous fixed point density matrix, and it has the right
    # dimensions, use that as the initial guess. Else, use a thermal density matrix.
    x0 = thermal_densitymatrix(m, Inf)
    old_rho = m.cache.previous_fixedpoint_densitymatrix
    if old_rho !== nothing && space(x0) == space(old_rho)
        x0 = old_rho
    end
    eigsolve_pars = get(pars, :scaleinvariant_krylovoptions, (;))
    solver = Arnoldi(; eager = true, eigsolve_pars...)
    _, vecs, vals, info = schursolve(f, x0, 1, :LM, solver)
    rho = vecs[1]
    # We know the result should always be Hermitian, and scaled to have trace 1.
    rho = (rho + rho')  # probably not even necessary
    rho /= tr(rho)
    m.cache.previous_fixedpoint_densitymatrix = rho
    if :verbosity in keys(pars) && pars[:verbosity] > 3
        msg = "Used $(info.numops) superoperator invocations to find the fixed point density matrix."
        @info(msg)
    end
    return rho
end

"""
    thermal_densitymatrix(m::GenericMERA, depth)

Return the thermal density matrix for the indices right below the layer at `depth`. Used as
an initial guess for the fixed-point density matrix.

See also: [`fixedpoint_densitymatrix`](@ref)
"""
function thermal_densitymatrix(m::GenericMERA, depth)
    width = causal_cone_width(typeof(m))
    V = ⊗(ntuple(n->inputspace(m, depth), Val(width))...)
    return convert(operatortype(m), id(storagetype(operatortype(m)), V))
end

"""
    densitymatrix(m::GenericMERA, depth, pars = (;))

Return the density matrix right below the layer at `depth`.

`pars` maybe a `NamedTuple` of options, passed on to the function
`fixedpoint_densitymatrix`.

This function utilises the cache, to avoid recomputation.

See also: [`fixedpoint_densitymatrix`](@ref), [`densitymatrices`](@ref)
"""
function densitymatrix(m::GenericMERA, depth, pars = (;))
    if !has_densitymatrix_stored(m.cache, depth)
        # If we don't find rho in storage, generate it.
        if depth > num_translayers(m)
            rho = fixedpoint_densitymatrix(m, pars)
        else
            rho_above = densitymatrix(m, depth+1, pars)
            rho = descend(rho_above, m, depth, depth+1)
        end
        # Store this density matrix for future use.
        set_stored_densitymatrix!(m.cache, rho, depth)
    end
    return get_stored_densitymatrix(m.cache, depth)
end

"""
    densitymatrices(m::GenericMERA, pars = (;))

Return all the distinct density matrices of the MERA, starting with the one at the physical
level, and ending with the scale invariant one.

`pars` maybe a `NamedTuple` of options, passed on to the function
`fixedpoint_densitymatrix`.

See also: [`fixedpoint_densitymatrix`](@ref), [`densitymatrix`](@ref)
"""
function densitymatrices(m::GenericMERA, pars = (;))
    rhos = [densitymatrix(m, depth, pars) for depth in 1:num_translayers(m)+1]
    return rhos
end

"""
    ascended_operator(op, m::GenericMERA, depth)

Return the operator `op` ascended from the physical level to `depth`.

This function utilises the cache, to avoid recomputation.

See also: [`scale_invariant_operator_sum`](@ref)
"""
function ascended_operator(op, m::GenericMERA, depth)
    # Note that if depth = 1, has_operator_stored always returns true, as it initializes
    # storage for this operator.
    if !has_operator_stored(m.cache, op, depth)
        op_below = ascended_operator(op, m, depth-1)
        opasc = ascend(op_below, m, depth, depth-1)
        # Store this density matrix for future use.
        set_stored_operator!(m.cache, opasc, op, depth)
    end
    return get_stored_operator(m.cache, op, depth)
end

"""
    scale_invariant_operator_sum(op, m::GenericMERA, pars)

Return the sum of the ascended versions of `op` in the scale invariant part of the MERA.

To be more precise, this sum is of course infinite, and what we return is the component of
it orthogonal to the dominant eigenoperator of the ascending superoperator (typically the
identity). This component converges to a finite result like a geometric series, since all
non-dominant eigenvalues of the ascending superoperator are smaller than 1.

To approximate the converging series, we use an iterative Krylov solver. The options for the
solver should be in `pars.scaleinvariant_krylovoptions`, they will be passed
to `KrylovKit.linsolve`.
"""
function scale_invariant_operator_sum(
    op, m::GenericMERA{N, LT, OT}, pars::NamedTuple = (;)
) where {N, LT, OT}
    nt = num_translayers(m)
    # fp is the dominant eigenvector of the ascending superoperator. We are not interested
    # in contributions to the sum along fp, since they will just be fp * infty, and fp is
    # merely the representation of the identity operator.
    fp = ascending_fixedpoint(getlayer(m, nt+1))
    rhop = densitymatrix(m, nt+1, pars)
    f = let fp = fp, rhop = rhop, lscaleinv = getlayer(m, nt+1)
        function (x)
            xasc = ascend(x, lscaleinv)
            return axpy!(-dot(rhop, xasc), fp, xasc)
        end
    end
    op_top = ascended_operator(op, m, nt+1)
    x0 = op_top
    old_opsum = m.cache.previous_operatorsum
    if old_opsum !== nothing && space(x0) == space(old_opsum)
        x0 = old_opsum
    end
    linsolve_pars = get(pars, :scaleinvariant_krylovoptions, (;))
    one_ = one(eltype(m))
    opsum, info = linsolve(f, op_top, x0, one_, -one_; linsolve_pars...)
    # We know the result should always be Hermitian.
    opsum = (opsum + opsum') / 2
    m.cache.previous_operatorsum = opsum
    # We are not interested in the component along fp.
    opsum = opsum - fp * dot(rhop, opsum)
    if :verbosity in keys(pars) && pars[:verbosity] > 3
        msg = "Used $(info.numops) superoperator invocations to find the scale invariant operator sum."
        @info(msg)
    end
    # TODO The type assertion should be removed once 1.6 is out. This is inferred fine
    # without it in 1.6, but not in 1.5, and I just use this to silence the failure in
    # tests.
    return opsum::operatortype(m)
end

"""
    environment(op, m::GenericMERA, depth, pars)

Return a `Layer` consisting of the environments of the various tensors of `m` at `depth`,
with respect to the expectation value of `op`. Note the return value isn't really a proper
MERA layer, i.e. the tensors are not isometric, it just has the same structure, and hence
the same data structure is used.

`pars` are parameters that are passed to `scale_invariant_operator_sum` and
`fixedpoint_densitymatrix`. `pars[:vary_disentanglers] = false` gives the option of
computing the environments only for the isometries, setting the environments of the
disentanglers to zero.

This function utilises the cache to avoid recomputation.

See also: [`fixedpoint_densitymatrix`](@ref), [`scale_invariant_operator_sum`](@ref)
"""
function environment(op, m::GenericMERA{N, LT, OT}, depth, pars) where {N, LT, OT}
    if !has_environment_stored(m.cache, op, depth)
        if depth <= num_translayers(m)
            op_below = ascended_operator(op, m, depth)
        else
            op_below = scale_invariant_operator_sum(op, m, pars)
        end
        op_below = normalise_hamiltonian(op_below)
        rho_above = densitymatrix(m, depth+1, pars)
        l = getlayer(m, depth)
        env = environment(op_below, l, rho_above;
                          vary_disentanglers = pars[:vary_disentanglers])
        set_stored_environment!(m.cache, env, op, depth)
    end
    return get_stored_environment(m.cache, op, depth)
end

# # # Extracting CFT data

"""
    scalingdimensions(m::GenericMERA, howmany = 20)

Diagonalize the scale invariant ascending superoperator to compute the scaling dimensions of
the underlying CFT.

The return value is a dictionary, the keys of which are symmetry sectors for a possible
internal symmetry of the MERA (`Trivial()` if there is no internal symmetry), and values are
scaling dimensions in this symmetry sector.

`howmany` controls how many of lowest scaling dimensions are computed.
"""
function scalingdimensions(m::GenericMERA, howmany = 20)
    width = causal_cone_width(m)
    V = inputspace(m, Inf)
    interlayer_space = ⊗(ntuple(n->V, Val(width))...)
    scaldim_dict = Dict()
    l = getlayer(m, num_translayers(m)+1)
    # Scaling operators come with charges, depending on the symmetry that the MERA tensors
    # have. If the symmetry is anyonic, then each scaling operator is labeled by two irreps,
    # otherwise by only a single irrep. We handle these cases separately, and diagonalise in
    # each charge sector separately, collecting the results in scaldim_dict.
    if BraidingStyle(sectortype(spacetype(layertype(m)))) == Anyonic()
        sectors = values(sectortype(interlayer_space))
        for irrep1 in sectors, irrep2 in sectors
            # Don't even try to get more than half of the eigenvalues. Its too expensive,
            # and they are garbage anyway.
            d = sum(blockdim(interlayer_space, irrep) for irrep in irrep1 ⊗ irrep2)
            d == 0 && continue
            howmanysector = min(howmany, Int(div(d, 2, RoundUp)))
            x0 = scalingoperator_initialguess(l, irrep1, irrep2)
            scaldims = scalingdimensions_from_initialguess(m, x0, howmanysector)
            scaldim_dict[(irrep1, irrep2)] = scaldims
        end
    else
        sectors = blocksectors(interlayer_space)
        for irrep in sectors
            # Don't even try to get more than half of the eigenvalues. Its too expensive,
            # and they are garbage anyway.
            d = blockdim(interlayer_space, irrep)
            howmanysector = min(howmany, Int(div(d, 2, RoundUp)))
            x0 = scalingoperator_initialguess(l, irrep)
            scaldims = scalingdimensions_from_initialguess(m, x0, howmanysector)
            scaldim_dict[irrep] = scaldims
        end
    end
    return scaldim_dict
end

function scalingdimensions_from_initialguess(m, x0, howmany)
    # Define a function that takes an operator and ascends it once through the scale
    # invariant layer.
    # Closures are more type stable if they only depend on arguments of the function.
    f(x) = ascend(x, getlayer(m, num_translayers(m)+1))
    S, U, info = eigsolve(f, x0, howmany, :LM)
    # sfact is the ratio by which the number of sites changes at each coarse-graining.
    sfact = scalefactor(typeof(m))
    scaldims = sort(-log.(sfact, abs.(real(S))))
    return scaldims
end

# # # Evaluation

"""
    expect(op, m::GenericMERA, pars = (;), opscale = 1, evalscale = 1)

Return the expecation value of operator `op` for the MERA `m`.

The layer on which `op` lives is set by `opscale`, which by default is the physical one.
`evalscale` can be used to set whether the operator is ascended through the network or the
density matrix is descended (the default). `pars` can hold additional options that are
further down the line passed on to `fixedpoint_densitymatrix`.

See also: [`fixedpoint_densitymatrix`](@ref)
"""
function expect(op, m::GenericMERA, pars = (;), opscale = 1, evalscale = 1)
    rho = densitymatrix(m, evalscale, pars)
    op = ascended_operator(op, m, evalscale)
    value = dot(rho, op)
    if abs(imag(value)/(norm(op)*norm(rho))) > 1e-13
        @warn("Non-real expectation value: $value")
    end
    return real(value)
end

# # # Optimization

const default_pars = (
    method = :lbfgs,
    retraction = :exp,
    transport = :exp,
    metric = :euclidean,
    precondition = true,
    vary_disentanglers = true,
    gradient_delta = 1e-12,
    isometries_only_iters = 0,
    maxiter = 2000,
    ev_layer_iters = 1,
    ls_epsilon = 1e-6,
    lbfgs_m = 8,
    cg_flavor = :HagerZhang,
    verbosity = 2,
    scaleinvariant_krylovoptions = (
        tol = 1e-13,
        verbosity = 0,
        maxiter = 20,
    ),
)

"""
    minimize_expectation(m::GenericMERA, h, pars = (;); finalize! = OptimKit._finalize!)

Return a MERA optimized to minimize the expectation value of operator `h`, starting with `m`
as the initial guess.

`pars` is a `NamedTuple` of parameters for the optimisation. They are,
* `method`: A `Symbol` that chooses which optimisation method to use. Options are `:lbfgs`
  for L-BFGS (default), `:ev` for Evenbly-Vidal, `:cg` for conjugate gradient, and `:gd` for
  gradient descent. `:lbfgs`, `:cg`, and `:gd` are together known as the gradient methods.
* `maxiter`: Maximum number of iterations to use. 2000 by default.
* `gradient_delta`: Convergence threshold, as measured by the norm of the gradient. `1e-14`
  by default.
* `precondition`: Whether to apply preconditioning with the physical Hilbert space inner
  product. `true` by default. See https://arxiv.org/abs/2007.03638 for more details.
* `verbosity`: How much output to log. 2 by default.
* `isometries_only_iters`: An integer for how many iterations should at first be done
  optimising only the isometries, and leaving the disentanglers be. 0 by default.
* `scaleinvariant_krylovoptions`: A `NamedTuple` of keyword arguments passed to
  `KrylovKit.linsolve` and `KrylovKit.eigsolve`, when solving for the fixed-point density
  matrix and the scale invariant operator sum. The default is
  `(tol = 1e-13, verbosity = 0, maxiter = 20)`.
* `retraction`: Which retraction method to use. Options are `:exp` for geodesics
  (default), and `:cayley` for Cayley transforms. Only affects gradient methods.
* `transport`: Which vector transport method to use. Currently each `retraction` method only
  comes with a single compatible transport, so one should always set `transport` to be the
  same as `retraction`. This may change. Only affects gradient methods.
* `metric`: Which metric to use for Stiefel manifold. Options are `:euclidean`
  (default) and `:canonical`. Only affects gradient methods.
* `ls_epsilon`: The ϵ parameter for the Hager-Zhang line search. `1e-6` be default. Only
  affects gradient methods.
* `lbfgs_m`: The rank of the approximation of the inverse Hessian in L-BFGS. 8 by default.
  Only affects the `:lbfgs` method.
* `cg_flavor`: The "flavor" of conjguate gradient to use. The options are `:HagerZhang` (the
  default), `:HestenesStiefel`, `:PolakRibierePolyak`, `:DaiYuan`, `:FletcherReeves`. Only
  affects the `:cg` method.
* `ev_layer_iters`: How many times a single layer is optimised before moving to the next
  layer in the Evenbly-Vidal algorithm. `1` by default. Only affects the `:ev` method.
* `vary_disentanglers`: If `false`, the optimisation is run only for the
  isometries, leaving the disentanglers as they are. `true` by default.
If any of these are specified in `pars`, the specified values override the defaults.

`finalize!` is a function that will be called at every iteration. It can be used to for
instance log the development of some quantity during the optimisation, or modify the MERA in
some custom way (although undefined behavior may follow depending on how the state is
changed).  Its signature is `finalize!(m, f, g, counter)` where `m` is the current MERA, `f`
is its expectation value for `h`, `g` is the gradient MERA at `m`, and `counter` is the
current iteration number. It should return the possibly modified `m`, `f`, and `g`. If
`method = :ev`, then it should also be able to handle the case `g = nothing`, since the
Evenbly-Vidal algorithm does not use gradients.

"""
function minimize_expectation(
    m::GenericMERA, h, pars = (;); finalize! = OptimKit._finalize!
)
    pars = merge(default_pars, pars)
    # If pars[:isometries_only_iters] is set, but the optimization on the whole is supposed
    # to vary_disentanglers too, then first run a pre-optimization without touching the
    # disentanglers, with pars[:isometries_only_iters] as the maximum iteration count,
    # before moving on to the main optimization with all tensors varying.
    if pars[:vary_disentanglers] && pars[:isometries_only_iters] > 0
        temp_pars = merge(pars, (maxiter = pars[:isometries_only_iters],
                                 vary_disentanglers = false))
        m = minimize_expectation(m, h, temp_pars; finalize! = finalize!)
    end

    method = pars[:method]
    if method in (:cg, :conjugategradient, :gd, :gradientdescent, :lbfgs)
        return minimize_expectation_grad(m, h, pars; finalize! = finalize!)
    elseif method == :ev || method == :evenblyvidal
        return minimize_expectation_ev(m, h, pars; finalize! = finalize!)
    else
        msg = "Unknown optimization method $(method)."
        throw(ArgumentError(msg))
    end
end

"""
    minimize_expectation_ev(m::GenericMERA, h, pars; finalize! = OptimKit._finalize!)

Return a MERA optimized with the Evenbly-Vidal method to minimize the expectation value of
operator `h`, starting with `m` as the initial guess. See [`minimize_expectation`](@ref) for
details.
"""
function minimize_expectation_ev(m::GenericMERA, h, pars; finalize! = OptimKit._finalize!)
    nt = num_translayers(m)
    rhos = densitymatrices(m, pars)
    expectation = expect(h, m, pars)
    rhos_maxchange = Inf
    gradnorm = Inf
    counter = 0
    if pars[:verbosity] >= 2
        @info(@sprintf("E-V: initializing with f = %.12f,", expectation))
    end
    while gradnorm > pars[:gradient_delta] && counter < pars[:maxiter]
        counter += 1
        old_rhos = rhos
        old_expectation = expectation
        gradnorm_sq = 0.0

        for i in 1:nt+1
            local env, l
            for j in 1:pars[:ev_layer_iters]
                env = environment(h, m, i, pars)
                l = getlayer(m, i)
                new_l = minimize_expectation_ev(
                    l, env; vary_disentanglers = pars[:vary_disentanglers]
                )
                m = replace_layer(m, new_l, i)
            end
            # We use the latest env and the corresponding layer to compute the norm of the
            # gradient. This isn't quite the gradient at the end point, which is what we
            # would want, but close enough.
            gradnorm_sq += gradient_normsq(l, env; metric = pars[:metric])
        end

        gradnorm = sqrt(gradnorm_sq)
        rhos = densitymatrices(m, pars)
        expectation = expect(h, m, pars)
        rho_diffs = [norm(r - or) for (r, or) in zip(rhos, old_rhos)]
        rhos_maxchange = maximum(rho_diffs)
        # The nothing is for the gradient, which we don't use, but OptimKit's format of
        # finalize! expects.
        m = finalize!(m, expectation, nothing, counter)[1]
        if pars[:verbosity] >= 2
            @info(@sprintf(
                "E-V: iter %4d: f = %.12f, ‖∇f‖ = %.4e, max‖Δρ‖ = %.4e",
                counter, expectation, gradnorm, rhos_maxchange
            ))
        end
    end
    if pars[:verbosity] > 0
        if gradnorm <= pars[:gradient_delta]
            @info(@sprintf(
                "E-V: converged after %d iterations: f = %.12f, ‖∇f‖ = %.4e, max‖Δρ‖ = %.4e",
                counter, expectation, gradnorm, rhos_maxchange
            ))
        else
            @warn(@sprintf(
                "E-V: not converged to requested tol: f = %.12f, ‖∇f‖ = %.4e, max‖Δρ‖ = %.4e",
                expectation, gradnorm, rhos_maxchange
            ))
        end
    end
    return m
end

"""
    normalise_hamiltonian(h)

Change the additive normalisation of a local Hamiltonian term `h` to make it suitable for
computing an Evenbly-Vidal update environment.

More specifically, return `h - ub*one(h)`, where `one(h)` is the identity operator, and
`ub` is an upper bound for the largest eigenvalue of `h`
"""
function normalise_hamiltonian(h)
    lb, ub = gershgorin_bounds(h)
    return h - ub*one(h)
end

# # # Gradient optimization

"""
    tensorwise_scale(m::GenericMERA, alpha::Number)

Scale all the tensors of `m` by `alpha`.
"""
function tensorwise_scale(m::GenericMERA, alpha::Number)
    return GenericMERA(tensorwise_scale.(m.layers, alpha))
end

"""
    tensorwise_sum(m1::T, m2::T) where T <: GenericMERA

Return a MERA for which each tensor is the sum of the corresponding tensors of `m1` and
`m2`.
"""
function tensorwise_sum(m1::T, m2::T) where T <: GenericMERA
    n = max(num_translayers(m1), num_translayers(m2)) + 1
    layers = ntuple(Val(n)) do i
        tensorwise_sum(getlayer(m1, i), getlayer(m2, i))
    end
    return GenericMERA(layers)
end

"""
    inner(m::GenericMERA, m1::GenericMERA, m2::GenericMERA; metric = :euclidean)

Given two tangent MERAs `m1` and `m2`, both at base point `m`, compute their inner product.
This means the sum of the inner products of the individual tensors.

See the documentation for `TensorKitManifolds` for more details.
"""
function TensorKitManifolds.inner(
    m::GenericMERA, m1::GenericMERA, m2::GenericMERA; metric = :euclidean
)
    n = max(num_translayers(m1), num_translayers(m2)) + 1
    res = sum([
        inner(getlayer(m, i), getlayer(m1, i), getlayer(m2, i); metric = metric)
        for i in 1:n
    ])
    return res
end

"""
    gradient(h, m::GenericMERA, pars::NamedTuple)

Compute the gradient of the expectation value of `h` at the point `m`.

`pars` should have `pars.metric` that specifies whether to use the `:euclidean` or
`:canonical` metric for the Stiefel manifold, and `pars.scaleinvariant_krylovoptions` that
is passed on to [`environment`](@ref). It should also have `pars.vary_disentanglers`, which
being `false` means that we should only compute gradients for the isometries, and set the
gradients for the disentanglers to zero.

The return value is a "tangent MERA": An object of a similar type as `m`, but instead of
regular layers with tensors that have isometricity constraints, each layer holds the
corresponding gradients for each tensor.
"""
function gradient(h, m::GenericMERA{N}, pars::NamedTuple) where {N}
    nt = num_translayers(m)
    layers = ntuple(Val(nt+1)) do i
        l = getlayer(m, i)
        env = environment(h, m, i, pars)
        gradient(l, env; metric = pars[:metric])
    end
    g = GenericMERA(layers)
    return g
end

"""
    precondition_tangent(m::GenericMERA, tan::GenericMERA, pars::NamedTuple)

Precondition the gradient `tan`, living at the base point `m`, using the physical Hilbert
space inner product.

For details on how the preconditioning is done, see https://arxiv.org/abs/2007.03638.

`pars` is passed on to [`densitymatrix`](@ref).
"""
function precondition_tangent(m::GenericMERA, tan::GenericMERA, pars::NamedTuple)
    nt = num_translayers(m)
    tanlayers_prec = ntuple(Val(nt+1)) do i
        l = getlayer(m, i)
        tanl = getlayer(tan, i)
        rho = densitymatrix(m, i+1, pars)
        precondition_tangent(l, tanl, rho)
    end
    tan_prec = GenericMERA(tanlayers_prec)
    return tan_prec
end

"""
    retract(m::GenericMERA, mtan::GenericMERA, alpha::Real; kwargs...)

Given a "tangent MERA" `mtan`, at base point `m`, retract in the direction of `mtan` by
distance `alpha`. This is done tensor-by-tensor, i.e. each tensor is retracted along its
respective Stiefel/Grassmann tangent.

The additional keyword argument are passed on to the respective `TensorKitManifolds`
function.

See the documentation for `TensorKitManifolds` for more details.

See also: [`transport!`](@ref)
"""
function TensorKitManifolds.retract(
    m::GenericMERA{N}, mtan::GenericMERA{N}, alpha::Real; kwargs...
) where N
    layers_and_layers_tan = retract.(m.layers, mtan.layers, alpha; kwargs...)
    layers = first.(layers_and_layers_tan)
    layers_tan = last.(layers_and_layers_tan)
    return GenericMERA(layers), GenericMERA(layers_tan)
end

"""
    transport!(mvec::GenericMERA, m::GenericMERA, mtan::GenericMERA, alpha::Real,
               mend::GenericMERA; kwargs...)

Given "tangent MERAs" `mtan` and `mvec`, at base point `m`, transport `mvec` in the
direction of `mtan` by distance `alpha`. This is done tensor-by-tensor, i.e. each tensor is
transported along its respective Stiefel/Grassmann tangent.

`mend` is the endpoint on the manifold, i.e. the result of retracting by `alpha` in the
direction of `mtan`. The additional keyword argument are passed on to the respective
`TensorKitManifolds` function.

See the documentation for `TensorKitManifolds` for more details.

See also: [`retract`](@ref)
"""
function TensorKitManifolds.transport!(
    mvec::GenericMERA, m::GenericMERA, mtan::GenericMERA, alpha::Real, mend::GenericMERA;
    kwargs...
)
    layers = transport!.(mvec.layers, m.layers, mtan.layers, alpha, mend.layers; kwargs...)
    return GenericMERA(layers)
end

"""
minimize_expectation_grad(m, h, pars; finalize! = OptimKit._finalize!)

Return a MERA optimized with one of the gradient methods to minimize the expectation value
of operator `h`, starting with `m` as the initial guess. See [`minimize_expectation`](@ref)
for details.
"""
function minimize_expectation_grad(m, h, pars; finalize! = OptimKit._finalize!)
    function fg(x)
        f = expect(h, x, pars)
        g = gradient(h, x, pars)
        return f, g
    end

    rtrct(args...) = retract(args...; alg = pars[:retraction])
    trnsprt!(args...) = transport!(args...; alg = pars[:transport])
    innr(args...) = inner(args...; metric = pars[:metric])
    scale(vec, beta) = tensorwise_scale(vec, beta)
    add(vec1, vec2, beta) = tensorwise_sum(vec1, scale(vec2, beta))
    linesearch = HagerZhangLineSearch(; ϵ = pars[:ls_epsilon])
    if pars[:precondition]
        precondition(x, g) = precondition_tangent(x, g, pars)
    else
        # The default that does nothing.
        precondition = OptimKit._precondition
    end

    algkwargs = (
        maxiter = pars[:maxiter],
        linesearch = linesearch,
        verbosity = pars[:verbosity],
        gradtol = pars[:gradient_delta]
    )
    local alg
    if pars[:method] == :cg || pars[:method] == :conjugategradient
        if pars[:cg_flavor] == :HagerZhang
            flavor = HagerZhang()
        elseif pars[:cg_flavor] == :HestenesStiefel
            flavor = HestenesStiefel()
        elseif pars[:cg_flavor] == :PolakRibierePolyak
            flavor = PolakRibierePolyak()
        elseif pars[:cg_flavor] == :DaiYuan
            flavor = DaiYuan()
        elseif pars[:cg_flavor] == :FletcherReeves
            flavor = FletcherReeves()
        else
            msg = "Unknown conjugate gradient flavor $(pars[:cg_flavor])"
            throw(ArgumentError(msg))
        end
        alg = ConjugateGradient(; flavor = flavor, algkwargs...)
    elseif pars[:method] == :gd || pars[:method] == :gradientdescent
        alg = GradientDescent(; algkwargs...)
    elseif pars[:method] == :lbfgs
        alg = LBFGS(pars[:lbfgs_m]; algkwargs...)
    else
        msg = "Unknown optimization method $(pars[:method])."
        throw(ArgumentError(msg))
    end
    res = optimize(
        fg, m, alg;
        scale! = scale,
        add! = add,
        retract = rtrct,
        inner = innr,
        transport! = trnsprt!,
        isometrictransport = true,
        precondition = precondition,
        finalize! = finalize!
    )
    m, expectation, normgrad, normgradhistory = res
    if pars[:verbosity] > 0
        @info("Gradient optimization done. Expectation = $(expectation).")
    end
    return m
end


text:
using ProblemParser
using Utils

using Chain
using DataStructures
using Dictionaries
using StaticArrays

const file = find_input(@__FILE__)
const input = parse(Rectangular(Convert(Int8)), slurp(file))

const NSEW = SVector{4,CartesianIndex}(CartesianIndex(lr...) for lr in ((0,-1),(0,1),(-1,0),(1,0)))

function part1(input)::Int
    pq = @chain begin
        input
        size
        CartesianIndices
        zip(_, Iterators.repeated(Inf))
        PriorityQueue
    end
    pq[CartesianIndex(1,1)] = input[1,1]
    prev = Dictionary{CartesianIndex,CartesianIndex}()

    while !isempty(pq)
        u,dist = dequeue_pair!(pq)

        for v in intersect(NSEW .+ [u], keys(pq))
            alt = dist + input[v]
            if alt < pq[v]
                pq[v] = alt
                set!(prev, v, u)
            end
        end
    end
    risk::Int = 0
    backtrace = CartesianIndex(size(input))
    while backtrace != CartesianIndex(1,1)
        risk += input[backtrace]
        backtrace = prev[backtrace]
    end
    risk
end

assertequal(part1(input), 498)

function part2(input)
    repeat = 5
    extended = hvcat(repeat, [mod1.(input .+ (i + j), 9) for i=0:repeat - 1 for j=0:repeat - 1]...)
    part1(extended)
end

assertequal(part2(input), 2901)

ext: jl

max_stars_repo_name: CmdQ/AoC


text:
mutable struct Texture
    sampler::Union{Nothing,Int}
    source::Union{Nothing,Int}
    name::Union{Nothing,String}
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function Texture()
        obj = new()
        obj.sampler = nothing
        obj.source = nothing
        obj.name = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

function Base.setproperty!(obj::Texture, sym::Symbol, x)
    if sym === :sampler && x !== nothing
        x ≥ 0 || throw(ArgumentError("the index of the sampler used by this texture should be ≥ 0"))
    elseif sym === :source && x !== nothing
        source ≥ 0 || throw(ArgumentError("the index of the image used by this texture should be ≥ 0"))
    end
    setfield!(obj, sym, x)
end

JSON3.StructType(::Type{Texture}) = JSON3.Mutable()
JSON3.omitempties(::Type{Texture}) = (:sampler, :source, :name, :extensions, :extras)


text:
#=
    A note on how we represent polynomials

    First, coefficients, exponents and polynomial ring metainfo
    are extracted from input polynomials with `convert_to_internal`.

    Inside the algorithm all exponents are hashed without collisions,
    so that an integer represents a single monomial.
    That drastically decreases memory consumption.

    For each monomial we also assign a divmask -- a compressed representation
    of exponent vector used to speed up divisibility checks.

    Hence, a polynomial is represented with coefficients vector
    together with vector of hashtable indices.

    Finally, the resulting basis in internal representation is dehashed,
    and passed to `convert_to_output`
=#

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

"""
    Contains info about polynomial ring
"""
mutable struct PolyRing
    # number of variables
    nvars::Int
    # raw length of exponent vector
    # (should be always equal to nvars+1)
    explen::Int
    # ring monomial ordering,
    # possible are :lex and :degrevlex
    ord::Symbol
    # characteristic of coefficient field
    ch::UInt64
    # information about the original ring of input. Options are:
    #    :abstract for AbstractAlgebra,
    #    :multivariate for MultivariatePolynomials, e.g, DynamicPolynomials,
    #    :hasparent for polynomials constructed with parent ring, e.g., Nemo
    origring::Symbol
end

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

"""
    Converts input polynomials to internal representation used by the algorithm.
    Extracts base ring information, exponents, and coefficients.

    This is the most general implementation.
    It happened to work for polynomials from
        . `Nemo`

    Currently, there are more efficient specializations for:
        . `AbstractAlgebra.MPoly`
        . `MultivariatePolynomials.AbstractPolynomial`
"""
function convert_to_internal(
            orig_polys::Vector{T},
            ordering::Symbol) where {T}
    isempty(orig_polys) && error("Empty input")
    ordering in (:input, :lex, :degrevlex, :deglex) || error("Not supported ordering $ordering")

    if hasmethod(parent, Tuple{typeof(first(orig_polys))})
        convert_to_internal(orig_polys, ordering, Val(:hasparent))
    else
        error("Sorry, we don't work with this type of polynomials yet. Feel free to open a github issue")
    end
end

function extract_ring(R::T) where {T}
    @assert hasmethod(nvars, Tuple{T})
    # @assert hasmethod(ordering, Tuple{typeof(R)})
    @assert hasmethod(characteristic, Tuple{T})

    nv     = nvars(R)
    explen = nv + 1
    ord = hasmethod(ordering, Tuple{T}) ? ordering(R) : :lex
    ch     = characteristic(R)

    @assert nv + 1 == explen
    @assert ord in (:lex, :degrevlex, :deglex)
    @assert 0 <= ch < 2^32

    PolyRing(nv, explen, ord, UInt64(ch), :hasparent)
end

function extract_coeffs(ring::PolyRing, orig_polys::Vector{T}) where {T}
    if ring.ch > 0
        extract_coeffs_ff(ring, orig_polys)
    else
        extract_coeffs_qq(ring, orig_polys)
    end
end

function extract_coeffs_ff(ring::PolyRing, poly::Poly)
    reverse(map(CoeffFF ∘ data, filter(!iszero, collect(coefficients(poly)))))
end

function extract_coeffs_qq(ring::PolyRing, poly::Poly)
    reverse(map(Rational, filter(!iszero, collect(coefficients(poly)))))
end

function extract_coeffs_ff(ring::PolyRing, poly)
    map(CoeffFF ∘ data, coefficients(poly))
end

function extract_coeffs_qq(ring::PolyRing, poly)
    map(Rational, coefficients(poly))
end

function extract_coeffs_ff(ring::PolyRing, orig_polys::Vector{T}) where {T}
    npolys = length(orig_polys)
    coeffs = Vector{Vector{CoeffFF}}(undef, npolys)
    for i in 1:npolys
        coeffs[i] = extract_coeffs_ff(ring, orig_polys[i])
    end
    coeffs
end

function extract_coeffs_qq(ring::PolyRing, orig_polys::Vector{T}) where {T}
    npolys = length(orig_polys)
    coeffs = Vector{Vector{CoeffQQ}}(undef, npolys)
    for i in 1:npolys
        coeffs[i] = extract_coeffs_qq(ring, orig_polys[i])
    end
    coeffs
end

function extract_exponents(ring, poly)
    exps = Vector{ExponentVector}(undef, length(poly))
    @inbounds for j in 1:length(poly)
        exps[j] = ExponentVector(undef, ring.explen)
        exps[j][1:ring.nvars] .= exponent_vector(poly, j)
        exps[j][end] = sum(exps[j][k] for k in 1:ring.nvars)
    end
    exps
end

function extract_exponents(ring, poly::Poly)
    # Why define length of univeriate polynomial to be the dense length??
    exps = Vector{ExponentVector}(undef, 0)
    @inbounds while !iszero(poly)
        push!(exps, ExponentVector(undef, ring.explen))
        exps[end][1] = degree(poly)
        exps[end][end] = degree(poly)
        poly = AbstractAlgebra.tail(poly)
    end
    exps
end

function extract_exponents(ring::PolyRing, orig_polys::Vector{T}) where {T}
    npolys = length(orig_polys)
    exps = Vector{Vector{ExponentVector}}(undef, npolys)
    for i in 1:npolys
        poly = orig_polys[i]
        exps[i] = extract_exponents(ring, poly)
    end
    exps
end

function extract_polys(ring::PolyRing, orig_polys::Vector{T}) where {T}
    cfs = extract_coeffs(ring, orig_polys)
    exps = extract_exponents(ring, orig_polys)
    exps, cfs
end

function convert_to_internal(
            orig_polys::Vector{T},
            ordering::Symbol,
            ::Val{:hasparent}) where {T}

    R = parent(first(orig_polys))
    ring = extract_ring(R)
    exps, cfs = extract_polys(ring, orig_polys)
    ring, exps, cfs
end

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

function extract_ring(
        orig_polys::Vector{<:AbstractPolynomialLike{T}}) where {T}

    f = first(orig_polys)

    nv = Groebner.MultivariatePolynomials.nvariables(orig_polys)
    explen = nv + 1
    ord    = :deglex
    ch     = T <: Union{Integer, Rational} ? 0 : characteristic(parent(first(coefficients(f))))

    @assert nv + 1 == explen
    @assert ord in (:lex, :degrevlex, :deglex)
    @assert 0 <= ch < 2^31

    PolyRing(nv, explen, ord, UInt64(ch), :multivariate)
end

function extract_coeffs_qq(
            ring::PolyRing,
            orig_polys::Vector{T}) where {T<:AbstractPolynomialLike{U}} where {U}
    npolys = length(orig_polys)
    coeffs = Vector{Vector{CoeffQQ}}(undef, npolys)
    for i in 1:npolys
        poly = orig_polys[i]
        coeffs[i] = map(Rational, MultivariatePolynomials.coefficients(poly))
    end
    coeffs
end

function exponents_wrt_vars(t, var2idx)
    exp = zeros(UInt16, length(var2idx))
    @inbounds for (v, p) in Groebner.MultivariatePolynomials.powers(t)
        exp[var2idx[v]] = p
    end
    exp
end

multivariate_length(p::MultivariatePolynomials.AbstractMonomialLike) = 1
multivariate_length(p::MultivariatePolynomials.AbstractTermLike) = 1
multivariate_length(p::AbstractPolynomialLike) = length(p)

function extract_exponents(
            ring::PolyRing,
            orig_polys::Vector{T}) where {T<:AbstractPolynomialLike{U}} where {U}

    npolys = length(orig_polys)
    exps = Vector{Vector{ExponentVector}}(undef, npolys)
    vars = MultivariatePolynomials.variables(orig_polys)
    @assert issorted(vars, rev=true)

    var2idx = Dict(vars[i] => i for i in 1:length(vars))
    for i in 1:npolys
        poly = orig_polys[i]
        exps[i] = Vector{ExponentVector}(undef, multivariate_length(poly))
        @inbounds for (j, t) in enumerate(MultivariatePolynomials.monomials(poly))
            exps[i][j] = ExponentVector(undef, ring.explen)
            exps[i][j][end] = zero(exps[i][j][end])
            et = exponents_wrt_vars(t, var2idx)
            @inbounds for ei in 1:ring.nvars
                exps[i][j][ei] = et[ei]
                exps[i][j][end] += exps[i][j][ei]
            end
        end
    end
    exps
end

"""
    `MultivariatePolynomials.AbstractPolynomial` conversion specialization
"""
function convert_to_internal(
        orig_polys::Vector{T},
        ordering::Symbol) where {T<:AbstractPolynomialLike{U}} where {U}

    isempty(orig_polys) && error("Empty input")
    ordering in (:input, :lex, :degrevlex, :deglex) || error("Not supported ordering $ordering")

    ring = extract_ring(orig_polys)
    exps, cfs = extract_polys(ring, orig_polys)
    ring, exps, cfs
end

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

function extract_exponents(
        ring::PolyRing,
        orig_polys::Vector{T},
        ::Ord) where {Ord<:Union{Val{:lex}, Val{:deglex}}} where {T}

    npolys = length(orig_polys)
    exps   = Vector{Vector{ExponentVector}}(undef, npolys)
    for i in 1:npolys
        poly = orig_polys[i]
        exps[i] = Vector{ExponentVector}(undef, length(poly))
        @inbounds for j in 1:length(poly)
            exps[i][j] = ExponentVector(undef, ring.explen)
            exps[i][j][end] = zero(exps[i][j][end])
            @inbounds for ei in 1:ring.nvars
                exps[i][j][ei] = poly.exps[ring.nvars - ei + 1, j]
                exps[i][j][end] += exps[i][j][ei]
            end
        end
    end
    exps
end

function extract_exponents(
        ring::PolyRing,
        orig_polys::Vector{T},
        ::Val{:degrevlex}) where {T}

    npolys = length(orig_polys)
    exps   = Vector{Vector{ExponentVector}}(undef, npolys)
    for i in 1:npolys
        poly = orig_polys[i]
        exps[i] = Vector{ExponentVector}(undef, length(poly))
        @inbounds for j in 1:length(poly)
            exps[i][j] = poly.exps[:, j]
        end
    end
    exps
end

function extract_polys(
        ring::PolyRing,
        orig_polys::Vector{T},
        ord::Symbol) where {T}
    cfs = extract_coeffs(ring, orig_polys)
    exps = extract_exponents(ring, orig_polys, Val(ord))
    exps, cfs
end

"""
    `AbstractAlgebra.MPoly` conversion specialization
"""
function convert_to_internal(
        orig_polys::Vector{MPoly{T}},
        ordering::Symbol) where {T}
    isempty(orig_polys) && error("Empty input")

    R = parent(first(orig_polys))
    ring = extract_ring(R)
    exps, cfs = extract_polys(ring, orig_polys, ring.ord)
    ring, exps, cfs
end

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

"""
    Converts internal polynomials for export as elements of `origring`.

    This is the most general implementation.
    It happened to work for polynomials from
        . `Nemo`

    There are also more efficient specializations for:
        . `AbstractAlgebra.MPoly`
        . `MultivariatePolynomials.AbstractPolynomial`
"""
function convert_to_output(
            ring::PolyRing,
            origpolys::Vector{P},
            gbexps::Vector{Vector{ExponentVector}},
            gbcoeffs::Vector{Vector{I}},
            metainfo::GroebnerMetainfo) where {P, I<:Coeff}

    if ring.origring == :abstract
        convert_to_output(ring, parent(first(origpolys)), gbexps, gbcoeffs, metainfo)
    elseif ring.origring == :multivariate
        convert_to_output(ring, origpolys, gbexps, gbcoeffs, metainfo)
    elseif ring.origring == :hasparent
        convert_to_output(ring, parent(first(origpolys)), gbexps, gbcoeffs, metainfo)
    else
        # this actually never happens
        error("Sorry, unknown polynomial ring.")
    end
end

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

checkexact(c, T::Type{BigInt}) = true
checkexact(c, T::Type{Rational{U}}) where {U} = checkexact(numerator(c), U) && checkexact(denominator(c), U)
function checkexact(c, T)
    if typemin(T) <= c <= typemax(T)
        return true
    else
        throw(DomainError(c, "Coefficient $c in the final basis cannot be converted exactly to $T. Using big arithmetic in input should fix this."))
        return false
    end
end

function check_and_convert_coeffs(coeffs_zz, T)
    cfs = Vector{T}(undef, length(coeffs_zz))
    for i in 1:length(coeffs_zz)
        checkexact(coeffs_zz[i], T)
        cfs[i] = coeffs_zz[i]
    end
    cfs
end

function convert_coeffs_to_output(
        polycoeffs::Vector{CoeffQQ},
        ::Type{T}) where {T<:Rational}
    check_and_convert_coeffs(polycoeffs, T)
end

function convert_coeffs_to_output(
        polycoeffs::Vector{CoeffQQ},
        ::Type{T}) where {T<:Integer}
    coeffs_zz = scale_denominators(polycoeffs)
    check_and_convert_coeffs(coeffs_zz, T)
end

"""
    `multivariate` conversion specialization
"""
function convert_to_output(
            ring::PolyRing,
            origpolys::Vector{P},
            gbexps::Vector{Vector{ExponentVector}},
            gbcoeffs::Vector{Vector{I}},
            metainfo::GroebnerMetainfo) where {P<:AbstractPolynomialLike{J}, I<:Coeff} where {J}

    # TODO: hardcoded
    (metainfo.targetord != :deglex) && @warn "Input polynomial type does not support ordering $(metainfo.targetord). Computed basis is in $(metainfo.targetord), but terms are ordered in $(:deglex) in output"

    origvars = MultivariatePolynomials.variables(origpolys)
    # xd
    T = typeof(origpolys[1] + origpolys[1])
    exported = Vector{T}(undef, length(gbexps))
    for i in 1:length(gbexps)
        cfs::Vector{J} = convert_coeffs_to_output(gbcoeffs[i], J)
        expvectors = [map(Int, gbexps[i][j][1:end-1]) for j in 1:length(gbexps[i])]
        expvars = map(t -> t[1]*prod(map(^, origvars, t[2])), zip(cfs, expvectors))
        exported[i] = sum(expvars)
    end
    exported
end

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

"""
    `hasparent` conversion specialization
"""
function convert_to_output(
            ring::PolyRing,
            origring::M,
            gbexps::Vector{Vector{ExponentVector}},
            gbcoeffs::Vector{Vector{I}},
            metainfo::GroebnerMetainfo) where {M, T, I}

    @assert hasmethod(base_ring, Tuple{typeof(origring)})

    # TODO: hardcoded
    if hasmethod(ordering, Tuple{M})
        (metainfo.targetord != ordering(origring)) && @warn "Unknown polynomial type. Computed basis is in $(metainfo.targetord), but terms are ordered in $(ordering(origring)) in output"
    end

    etype = elem_type(base_ring(origring))
    # rather weak but okay for now
    if etype <: Integer
        coeffs_zz = scale_denominators(gbcoeffs)
        convert_to_output(origring, gbexps, coeffs_zz, metainfo)
    else
        convert_to_output(origring, gbexps, gbcoeffs, metainfo)
    end
end

function convert_to_output(
            origring::M,
            gbexps::Vector{Vector{ExponentVector}},
            gbcoeffs::Vector{Vector{I}},
            metainfo::GroebnerMetainfo) where {M<:AbstractAlgebra.Generic.PolyRing, I}

    ground   = base_ring(origring)
    exported = Vector{elem_type(origring)}(undef, length(gbexps))
    @inbounds for i in 1:length(gbexps)
        cfs = Vector{}
        cfs    = zeros(ground, gbexps[i][1][1] + 1)
        for (idx, j) in enumerate(gbexps[i])
            cfs[j[1] + 1] = ground(gbcoeffs[i][idx])
        end
        exported[i] = origring(cfs)
    end
    exported
end

function convert_to_output(
            origring::M,
            gbexps::Vector{Vector{ExponentVector}},
            gbcoeffs::Vector{Vector{I}},
            metainfo::GroebnerMetainfo) where {M, I}

    ground   = base_ring(origring)
    exported = Vector{elem_type(origring)}(undef, length(gbexps))
    @inbounds for i in 1:length(gbexps)
        cfs    = map(ground, gbcoeffs[i])
        exps   = [Int.(gbexps[i][j][1:end-1]) for j in 1:length(gbexps[i])]
        exported[i] = origring(cfs, exps)
    end
    exported
end

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

function create_polynomial(
            origring::MPolyRing{T}, coeffs::Vector{T}, exps) where {T}
    ground   = base_ring(origring)
    if !isempty(coeffs)
        MPoly{elem_type(ground)}(origring, coeffs, exps)
    else
        MPoly{elem_type(ground)}(origring)
    end
end

"""
    `AbstractAlgebra.MPolyRing` conversion specialization
"""
function convert_to_output(
            ring::PolyRing,
            origring::MPolyRing{T},
            gbexps::Vector{Vector{ExponentVector}},
            gbcoeffs::Vector{Vector{I}},
            metainfo::GroebnerMetainfo) where {T, I}

    ord = ordering(origring)
    if metainfo.targetord != ord
        origring, _ = PolynomialRing(base_ring(origring), symbols(origring), ordering=metainfo.targetord)
    end

    if elem_type(base_ring(origring)) <: Integer
        coeffs_zz = scale_denominators(gbcoeffs)
        convert_to_output(origring, gbexps, coeffs_zz, Val(ordering(origring)))
    else
        convert_to_output(origring, gbexps, gbcoeffs, Val(ordering(origring)))
    end
end

"""
    Rational, Integer, and Finite field :degrevlex
    `AbstractAlgebra.MPolyRing` conversion specialization
"""
function convert_to_output(
            origring::MPolyRing{U},
            gbexps::Vector{Vector{ExponentVector}},
            gbcoeffs::Vector{Vector{T}},
            ::Val{:degrevlex}) where  {T<:Coeff, U}

    nv = nvars(origring)
    ground   = base_ring(origring)
    exported = Vector{elem_type(origring)}(undef, length(gbexps))
    for i in 1:length(gbexps)
        cfs    = map(ground, gbcoeffs[i])
        exps   = Matrix{UInt64}(undef, nv + 1, length(gbcoeffs[i]))
        @inbounds for jt in 1:length(gbcoeffs[i])
            for je in 1:nv
                exps[je, jt] = gbexps[i][jt][je]
            end
            exps[nv + 1, jt] = gbexps[i][jt][end]
        end
        exported[i] = create_polynomial(origring, cfs, exps)
    end
    exported
end

"""
    Rational, Integer, and Finite field :lex
    `AbstractAlgebra.MPolyRing` conversion specialization
"""
function convert_to_output(
            origring::MPolyRing{U},
            gbexps::Vector{Vector{ExponentVector}},
            gbcoeffs::Vector{Vector{T}},
            ::Val{:lex}) where {T<:Coeff, U}

    nv = nvars(origring)
    ground   = base_ring(origring)
    exported = Vector{elem_type(origring)}(undef, length(gbexps))
    for i in 1:length(gbexps)
        cfs    = map(ground, gbcoeffs[i])
        exps   = Matrix{UInt64}(undef, nv, length(gbcoeffs[i]))
        @inbounds for jt in 1:length(gbcoeffs[i])
            for je in 1:nv
                exps[je, jt] = gbexps[i][jt][nv - je + 1]
            end
        end
        # exps   = UInt64.(hcat(map(x -> x[end-1:-1:1], gbexps[i])...))
        exported[i] = create_polynomial(origring, cfs, exps)
    end
    exported
end

"""
    Rational, Integer, and Finite field :deglex
    `AbstractAlgebra.MPolyRing` conversion specialization
"""
function convert_to_output(
            origring::MPolyRing{U},
            gbexps::Vector{Vector{ExponentVector}},
            gbcoeffs::Vector{Vector{T}},
            ::Val{:deglex}) where {T<:Coeff, U}

    nv = nvars(origring)
    ground   = base_ring(origring)
    exported = Vector{elem_type(origring)}(undef, length(gbexps))
    for i in 1:length(gbexps)
        cfs    = map(ground, gbcoeffs[i])
        exps   = Matrix{UInt64}(undef, nv + 1, length(gbcoeffs[i]))
        @inbounds for jt in 1:length(gbcoeffs[i])
            for je in 1:nv
                exps[je, jt] = gbexps[i][jt][nv - je + 1]
            end
            exps[nv + 1, jt] = gbexps[i][jt][end]
        end
        exported[i] = create_polynomial(origring, cfs, exps)
    end
    exported
end

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


text:
using TSeriesForecast
using Test

tests = [
    "simpleexponentialsmoothing",
    "holts_trend_method",
    "seasonality_exponential_smoothing"
]

@testset "TSeriesForecast.jl" begin
    @info("Running tests:")

    for test ∈ tests
        @info("\t* $test ...")
        include("./$test.jl")
    end
end


text:
# Helper functions for spherical geometry

using Rotations

"""
    c2s(x, y, z) -> [θ, ϕ]
    c2s(xyz::AbstractVector) -> [θ, ϕ]
    c2s(xyz::AbstractMatrix) -> [θs; ϕs]

Convert Cartesian vector `[x, y, z]` (or `xyz`) to spherical unit vector `[θ,
ϕ]`.  For AbstractMatrix `xyz`, each column is converted and the results are
reassmebled into a Matrix using `hcat`.
"""
function c2s(x,y,z)
    θ=atan(y,x)
    ϕ=atan(z, hypot(x,y))
    [θ, ϕ]
end

function c2s(xyz::AbstractVector{<:Real})
    c2s(xyz[1], xyz[2], xyz[3])
end

function c2s(xyz::AbstractMatrix{<:Real})
    @views reduce(hcat, c2s.(xyz[1,:], xyz[2,:], xyz[3,:]))
end

"""
    s2c(θ, ϕ) -> [x, y, z)
    s2c(θϕ::AbstractVector) -> [x, y, z]
    s2c(θϕ::AbstractMaxtrix) -> [xs; ys; zs]

Convert spherical unit vector `[θ, ϕ]` (or `θϕ`) to Cartesian unit vector `[x,
y, z]`.  For AbstractMatrix `θϕ`, each column is converted and the results are
reassmebled into a Matrix using `hcat`.
"""
function s2c(θ,ϕ)
    z, r = sincos(ϕ)
    y, x = sincos(θ) .* r
    [x, y, z]
end

function s2c(θϕ::AbstractVector{<:Real})
    s2c(θϕ[1], θϕ[2])
end

function s2c(θϕ::AbstractMatrix{<:Real})
    @views reduce(hcat, s2c.(θϕ[1,:], θϕ[2,:]))
end

"""
    beamrings(θ, ϕ; nrings=4, dϕ=deg2rad(10/3600)) -> Matrix

Compute the centers of beams that are arranged in concentric rings around the
spherical coordinates `(θ, ϕ)`, both in radians with `θ` being the longitudinal
coordinate and `ϕ` being the latitudinal coordinate.  The number of rings for
which beam centers are returned is given by `nrings`, which defaults to 4.  The
rings are radially separated by `dϕ` radians, which defaults to the equivalent
of 10 arcseconds.  The beams returned do not form a packed hexagon because that
is not possible on a sperical sky, but the number of beams per ring does follow
the hexagonal packing convention of `6n` beams for ring `n` (when `n>0`).  The
so-called boresight beam center at `(θ, ϕ)` is included as well, so 61 beam
centers are returned for `nrings==4`.  More generally,
`nbeams=6*sum(1:nrings)+1` beam centers are returned in a `2xnbeams` Matrix,
where the first row contains the `θ_beam` coordinates and the second row
contains the `ϕ_beam` coordinates.
"""
function beamrings(θ, ϕ; nrings=4, dϕ=deg2rad(10/3600))
    nbeams = 6 * sum(1:nrings) + 1
    θs = Vector{Float64}(undef, nbeams)
    ϕs = Vector{Float64}(undef, nbeams)
    θs[1] = 0.0
    ϕs[1] = π/2
    b=2
    for r in 1:nrings
        nθ = 6r
        dθ = 2π/nθ
        ϕr = π/2 - r * dϕ
        for i in 0:nθ-1
            θs[b] = i*dθ
            ϕs[b] = ϕr
            b += 1
        end
    end
    reduce(hcat, c2s.(Ref(RotZY(θ, π/2-ϕ)) .* s2c.(θs, ϕs)))
end


text:
using Documenter, CreditMetrics

makedocs(sitename="CreditMetrics.jl",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    modules = [CreditMetrics],
    pages = ["Home" => "index.md", "Reference" => "reference.md"]
)

deploydocs(
    repo = "github.com/kasperrisager/CreditMetrics.jl.git",
)


text:
using IterativeSolvers
abstract type ForcingSequence end

struct FixedForceTerm{T} <: ForcingSequence
    η::T
end
η(fft::FixedForceTerm, info) = fft.η

struct DemboSteihaug <: ForcingSequence
end
η(fft::DemboSteihaug, info) = min(1/(info.k+2), info.normFx)

struct EisenstatWalkerA{T} <: ForcingSequence
    α::T
end
function η(fft::EisenstatWalkerA, info)
    α = fft.α

    ξ = (info.normFx - info.residual_old)/info.normFx_old
    β = info.η_old^(1+sqrt(5)/2)
    η = β ≦ α ? η : max(ξ, β)
end

struct EisenstatWalkerB{T} <: ForcingSequence
    α::T
    γ::T
    ω::T
    η₀::T
end
function η(fft::EisenstatWalkerB, info)
    γ, ω, α = fft.γ, fft.ω, fft.α
    T = typeof(info.normFx)

    ρ = info.normFx/info.normFx_old
    ξ = γ*ρ^ω
    β = info.η_old^(1+sqrt(5)/2)
    η = β ≦ α ? η : max(ξ, β)
end

struct ResidualKrylov{ForcingType<:ForcingSequence, Tη}
    force_seq::ForcingType
    η₀::Tη
    itermax::Int
end
ResidualKrylov(; force_seq=FixedForceTerm(1e-4), eta0 = 1e-4, itermax=300)=ResidualKrylov(force_seq, eta0, itermax)
# map from method to forcing sequence
η(fft::ResidualKrylov, info) = η(fft.force_seq, info)

# Default to EisenstatA-values. The paper [[[]]] suggests that γ ≥ 0.7 and
# ω ≥ (1+sqrt(5))/2
EisenstatWalkerB() = EisenstatWalkerB(0.1, 1.0, 1+sqrt(5)/2, 0.5)

struct ResidualKrylovOp{T1, T2}
    F::T1
    JacOp::T2
end
function ResidualKrylovOp(F; seed, autodiff=false)
    JacOp = JacVec(F, seed; autodiff=false)
    ResidualKrylovOp(F, JacOp)
end

function nlsolve!(kp::ResidualKrylovOp, x, method::ResidualKrylov; f_abstol=1e-8, f_reltol=1e-12)
    xp, Fx = copy(x), copy(x)

    Fx = value!(kp, Fx, x)
    normFx = norm(Fx, 2)

    stoptol = T(f_reltol)*normFx + T(f_abstol)
    force_info = (k = 1, normFx=normFx, normFx_old=nothing, η_old=nothing)

    for i = 1:method.itermax
        # Update the point of evaluation
        jacvec_fn(kp).u .= x

        # Refactor this
        ηₖ = i == 1 && !isa(method.force_seq, FixedForceTerm) ? method.η₀ : η(method, force_info)

        krylov_iter = IterativeSolvers.gmres_iterable!(xp, jacvec_fn(kp), Fx; maxiter=20)

        rhs = ηₖ*norm(Fx)
        for item in krylov_iter
            res = krylov_iter.residual.current
            if res <= rhs
                break
            end
        end
        x .= x .- xp
        Fx = value!(kp, Fx, x)

        if norm(Fx, 2) < stoptol
            break
        end
        normFx_old = force_info.normFx
        η_old = ηₖ
        normFx = norm(Fx, 2)
        force_info = (k = i, normFx=normFx, normFx_old=normFx_old, η_old=η_old, residual_old=res)
    end
    return x, Fx
end
value!(kp::ResidualKrylovOp, Fx, x) = kp.F(Fx, x)
value_fn(kp) = kp.F
jacvec_fn(kp) = kp.JacOp


text:
export UnspecifiedFormConstraint

"""
    UnspecifiedFormConstraint

One of the form constraint objects. Does not imply any form constraints and simply returns the same object as receives.
However it does not allow `DistProduct` to be a valid functional form in the inference backend.

See also: [`constrain_form`](@ref), [`DistProduct`](@ref)
"""
struct UnspecifiedFormConstraint end

default_form_check_strategy(::UnspecifiedFormConstraint) = FormConstraintCheckLast()

is_point_mass_form_constraint(::UnspecifiedFormConstraint) = false

constrain_form(::UnspecifiedFormConstraint, something)                            = something
constrain_form(::UnspecifiedFormConstraint, something::Message{ <: DistProduct }) = error("`DistProduct` object cannot be used as a functional form in inference backend. Use form constraints to restrict the functional form of marginal posteriors.")


text:
function sign_change(x, sgn::Pair{Symbol,Symbol})
    if sgn == (:- => :+)
        idx = sign_change(x, true)
    else
        idx = sign_change(x, false)
    end

    # Adjust to closest value
    inext = idx == lastindex(x) ? firstindex(x) : idx+1
    iprev = idx == firstindex(x) ? lastindex(x) : idx-1
    if abs(x[inext]) < abs(x[idx])
        idx = inext
    elseif abs(x[iprev]) < abs(x[idx])
        idx = iprev
    end

    return idx
end

function sign_change(x, sgn)
    if signbit(x[end]) == sgn
        # we consider peroidic boundary conditions
        start_idx = firstindex(x)
        xprev = x[end]
    elseif signbit(x[begin]) == sgn
        # we start on the correct sign
        start_idx = firstindex(x) + 1
        xprev = x[begin]
    else
        start_idx = firstindex(x)
        while signbit(x[start_idx]) ≠ sgn
            start_idx +=1
        end
        xprev = x[start_idx]
        start_idx += 1
    end
    sprev = signbit(xprev)
    @assert sprev == sgn
    idx = start_idx
    offset_x = OffsetArray(x[start_idx:end], start_idx:lastindex(x))
    for i in eachindex(offset_x)
        current_sign = signbit(x[i])
        if (current_sign ≠ sprev)
            idx = i
            break
        end
        xprev = x[i]
        sprev = current_sign
    end
    return idx
end

function integration_domain(x, y, quadrant, direction)
    if direction == :trigonometric
        plus = :+
        minus = :-
    else
        plus = :-
        minus = :+
    end

    if quadrant == 1
        a_idx = sign_change(y, minus => plus)
        b_idx = sign_change(x, plus => minus)
    elseif quadrant == 2
        a_idx = sign_change(x, plus => minus)
        b_idx = sign_change(y, plus => minus)
    elseif quadrant == 3
        a_idx = sign_change(y, plus => minus)
        b_idx = sign_change(x, minus => plus)
    elseif quadrant == 4
        a_idx = sign_change(x, minus => plus)
        b_idx = sign_change(y, minus => plus)
    else
        @error "Unkown quadrant $quadrant"
        a_idx = 0
        b_idx = 0
    end

    return direction == :trigonometric ? (a_idx, b_idx) : (b_idx, a_idx)
end

function discharge_area(df, quadrant, fixed_ΔV)
    I, V = df[!, Symbol("Current (mA)")], df[!, Symbol("Potential (V)")]
    t = df[!, Symbol("Time (s)")]

    a_idx, b_idx = integration_domain(V, I, quadrant, :orar)
    @assert a_idx < b_idx

    I, V = I[a_idx:b_idx], V[a_idx:b_idx]
    # account for integrating "in reverse"
    sgn = V[begin] > V[end] ? -1 : 1
    ∫IdV = sgn * integrate(V, I) * u"V*A"
    Δt = -sgn*(t[b_idx] - t[a_idx]) * u"s"
    ΔV = !isnothing(fixed_ΔV) ? sgn*fixed_ΔV : (V[end] - V[begin]) * u"V"

    return Δt, ΔV, ∫IdV
end


text:
using ZygoteFFTs
using Test
using FillArrays
using FFTW
using Zygote

function ngradient(f, xs::AbstractArray...)
  grads = zero.(xs)
  for (x, Δ) in zip(xs, grads), i in 1:length(x)
    δ = sqrt(eps())
    tmp = x[i]
    x[i] = tmp - δ/2
    y1 = f(xs...)
    x[i] = tmp + δ/2
    y2 = f(xs...)
    x[i] = tmp
    Δ[i] = (y2-y1)/δ
  end
  return grads
end

gradcheck(f, xs...) =
  all(isapprox.(ngradient(f, xs...),
                gradient(f, xs...), rtol = 1e-5, atol = 1e-5))

gradtest(f, xs::AbstractArray...) = gradcheck((xs...) -> sum(sin.(f(xs...))), xs...)
gradtest(f, dims...) = gradtest(f, rand.(Float64, dims)...)

# utilities for using gradcheck with complex matrices
_splitreim(A) = (real(A),)
_splitreim(A::AbstractArray{<:Complex}) = reim(A)

_joinreim(A, B) = complex.(A, B)
_joinreim(A) = A

function _dropimaggrad(A)
  back(Δ) = real(Δ)
  back(Δ::Nothing) = nothing
  return Zygote.hook(back, A)
end

@testset "ZygoteFFTs.jl" begin
    findicateMat(i,j,n1,n2) = [(k==i) && (l==j) ? 1.0 : 0.0 for k=1:n1,
                               l=1:n2]
    mirrorIndex(i,N) = i - 2*max(0,i - (N>>1+1))
    for sizeX in [(2,3), (10,10), (13,15)]
        X = randn(sizeX)
        X̂r = rfft(X)
        X̂ = fft(X)
        N = prod(sizeX)
        for i=1:size(X,1), j=1:size(X,2)
            indicateMat = [(k==i) && (l==j) ? 1.0 : 0.0 for k=1:size(X, 1),
                           l=1:size(X,2)]
            # gradient of ifft(fft) must be (approximately) 1 (for various cases)
            @test gradient((X)->real.(ifft(fft(X))[i, j]), X)[1] ≈ indicateMat
            # same for the inverse
            @test gradient((X̂)->real.(fft(ifft(X̂))[i, j]), X̂)[1] ≈ indicateMat
            # same for rfft(irfft)
            @test gradient((X)->real.(irfft(rfft(X),
                                            size(X,1)))[i, j], X)[1] ≈
                                                real.(indicateMat)
            # rfft isn't actually surjective, so rffft(irfft) can't really
            # be tested this way.

            # the gradients are actually just evaluating the inverse
            # transform on the indicator matrix
            mirrorI = mirrorIndex(i,sizeX[1])
            FreqIndMat = findicateMat(mirrorI, j, size(X̂r,1), sizeX[2])
            listOfSols = [(fft, bfft(indicateMat), bfft(indicateMat*im),
                           plan_fft(X), i, X),
                          (ifft, 1/N*fft(indicateMat), 1/N*fft(indicateMat*im),
                           plan_fft(X), i, X),
                          (bfft, fft(indicateMat), fft(indicateMat*im),
                           nothing, i, X),
                          (rfft, real.(brfft(FreqIndMat, sizeX[1])),
                           real.(brfft(FreqIndMat*im, sizeX[1])), plan_rfft(X),
                           mirrorI, X),
                          ((K)->(irfft(K,sizeX[1])), 1/N * rfft(indicateMat),
                           zeros(size(X̂r)), plan_rfft(X), i, X̂r)]
            for (trans, solRe, solIm, P, mI, evalX) in listOfSols
                @test gradient((X)->real.(trans(X))[mI, j], evalX)[1] ≈
                    solRe
                @test gradient((X)->imag.(trans(X))[mI, j], evalX)[1] ≈
                    solIm
                if typeof(P) <:AbstractFFTs.Plan && maximum(trans .== [fft,rfft])
                    @test gradient((X)->real.(P * X)[mI, j], evalX)[1] ≈
                        solRe
                    @test gradient((X)->imag.(P * X)[mI, j], evalX)[1] ≈
                        solIm
                elseif typeof(P) <: AbstractFFTs.Plan
                    @test gradient((X)->real.(P \ X)[mI, j], evalX)[1] ≈
                        solRe
                    # for whatever reason the rfft_plan doesn't handle this
                    # case well, even though irfft does
                    if eltype(evalX) <: Real
                        @test gradient((X)->imag.(P \ X)[mI, j], evalX)[1] ≈
                            solIm
                    end
                end
            end
        end
    end

    x = [-0.353213 -0.789656 -0.270151; -0.95719 -1.27933 0.223982]
    for trans in (fft, ifft, bfft)
        @test gradient((x)->sum(abs.(trans(x))), x)[1] ≈
            gradient( (x) -> sum(abs.(trans(trans(x,1),2))),  x)[1]
        # switch sum abs order
        @test gradient((x)->abs(sum((trans(x)))),x)[1] ≈
            gradient( (x) -> abs(sum(trans(trans(x,1),2))),  x)[1]
        # dims parameter for the function
        @test gradient((x, dims)->sum(abs.(trans(x,dims))), x, (1,2))[1] ≈
            gradient( (x) -> sum(abs.(trans(x))), x)[1]
        # (1,2) should be the same as no index
        @test gradient( (x) -> sum(abs.(trans(x,(1,2)))), x)[1] ≈
            gradient( (x) -> sum(abs.(trans(trans(x,1),2))), x)[1]
        @test gradcheck(x->sum(abs.(trans(x))), x)
        @test gradcheck(x->sum(abs.(trans(x, 2))), x)
    end

    @test gradient((x)->sum(abs.(rfft(x))), x)[1] ≈
        gradient( (x) -> sum(abs.(fft(rfft(x,1),2))),  x)[1]
    @test gradient((x, dims)->sum(abs.(rfft(x,dims))), x, (1,2))[1] ≈
        gradient( (x) -> sum(abs.(rfft(x))), x)[1]
end


text:
using Pkg

Pkg.develop(path = "..")

using Publish
using Artifacts, LazyArtifacts
using Metalhead

# override default theme
cp(artifact"flux-theme", "../_flux-theme"; force = true)

p = Publish.Project(Metalhead)

function build_and_deploy(label)
    rm(label; recursive = true, force = true)
    deploy(Metalhead; root = "/Metalhead.jl", label = label)
end


text:
using LoopVectorization, OffsetArrays, Test
using LoopVectorization.VectorizationBase: StaticUnitRange
T = Float64
# T = Float32

@testset "OffsetArrays" begin

     function old2d!(out::AbstractMatrix, A::AbstractMatrix, kern)
        rng1k, rng2k = axes(kern)
        rng1,  rng2  = axes(out)
        for j in rng2, i in rng1
            tmp = zero(eltype(out))
            @inbounds for jk in rng2k, ik in rng1k
                tmp += oftype(tmp, A[i+ik,j+jk])*kern[ik,jk]
            end
            @inbounds out[i,j] = tmp
        end
        out
    end

    # out = out1;
    # R=CartesianIndices(out);
    # z=zero(eltype(out));
    # rng1k, rng2k = axes(skern);
    # rng1,  rng2  = R.indices;
    # tmp = z; i = 2; j = 2;
    # ls1 = LoopVectorization.@avx_debug for jk in rng2k, ik in rng1k
    #     tmp += A[i+ik,j+jk]*skern[ik,jk]
    # end;
    # ls1

    # out = out1;
    # rng1,  rng2  = CartesianIndices(out1).indices;
    # rng1k, rng2k = axes(kern);
    # # C = At';
    # ls2d = LoopVectorization.@avx_debug for j in rng2, i in rng1
    #         tmp = zero(eltype(out))
    #         for jk in rng2k, ik in rng1k
    #             tmp += A[i+ik,j+jk]*kern[ik,jk]
    #         end
    #         out1[i,j] = tmp
    # end;
    # LoopVectorization.choose_order(ls2d)

    # q2d = :(for j in rng2, i in rng1
    #         tmp = zero(eltype(out))
    #         for jk in rng2k, ik in rng1k
    #             tmp += A[i+ik,j+jk]*kern[ik,jk]
    #         end
    #         out[i,j] = tmp
    #        end);
    # lsq2d = LoopVectorization.LoopSet(q2d); LoopVectorization.choose_order(lsq2d)

    # oq2 = :(for j in rng2, i in rng1
    #         tmp = zero(eltype(out))
    #         for jk in -1:1, ik in -1:1
    #             tmp += A[i+ik,j+jk]*kern[ik,jk]
    #         end
    #         out[i,j] = tmp
    #        end);
    # lsoq = LoopVectorization.LoopSet(oq2);
    # LoopVectorization.choose_order(lsoq)

    function avx2d!(out::AbstractMatrix, A::AbstractMatrix, kern)
        rng1k, rng2k = axes(kern)
        rng1,  rng2  = axes(out)
        for j in rng2, i in rng1
            tmp = zero(eltype(out))
            @avx for jk in rng2k, ik in rng1k
                tmp += A[i+ik,j+jk]*kern[ik,jk]
            end
            out[i,j] = tmp
        end
        out
    end
    function avx2douter!(out::AbstractMatrix, A::AbstractMatrix, kern)
        rng1k, rng2k = axes(kern)
        rng1,  rng2  = axes(out)
        @avx for j in rng2, i in rng1
            tmp = zero(eltype(out))
            for jk in rng2k, ik in rng1k
                tmp += A[i+ik,j+jk]*kern[ik,jk]
            end
            out[i,j] = tmp
        end
        out
    end

    struct SizedOffsetMatrix{T,LR,UR,LC,RC} <: DenseMatrix{T}
        data::Matrix{T}
    end
    Base.axes(::SizedOffsetMatrix{T,LR,UR,LC,UC}) where {T,LR,UR,LC,UC} = (StaticUnitRange{LR,UR}(),StaticUnitRange{LC,UC}())
    Base.parent(A::SizedOffsetMatrix) = A.data
    @generated function LoopVectorization.stridedpointer(A::SizedOffsetMatrix{T,LR,UR,LC,RC}) where {T,LR,UR,LC,RC}
        quote
            $(Expr(:meta,:inline))
            LoopVectorization.OffsetStridedPointer(
                LoopVectorization.StaticStridedPointer{$T,Tuple{1,$(UR-LR+1)}}(pointer(parent(A))),
                ($(LR-1), $(LC-1))
            )
        end
    end
    # Base.size(A::SizedOffsetMatrix{T,LR,UR,LC,UC}) where {T,LR,UR,LC,UC} = (1 + UR-LR, 1 + UC-LC)
    # Base.CartesianIndices(::SizedOffsetMatrix{T,LR,UR,LC,UC}) where {T,LR,UR,LC,UC} = CartesianIndices((LR:UR,LC:UC))
    Base.getindex(A::SizedOffsetMatrix, i, j) = LoopVectorization.vload(LoopVectorization.stridedpointer(A), (i-1,j-1))
    function avx2dunrolled!(out::AbstractMatrix, A::AbstractMatrix, kern::SizedOffsetMatrix{T,-1,1,-1,1}) where {T}
        rng1,  rng2  = axes(out)
        Base.Cartesian.@nexprs 3 jk -> Base.Cartesian.@nexprs 3 ik -> kern_ik_jk = kern[ik-2,jk-2]
        # Manually unpack the OffsetArray
        @avx for j in rng2, i in rng1
            tmp_0 = zero(eltype(out))
            Base.Cartesian.@nexprs 3 jk -> Base.Cartesian.@nexprs 3 ik -> tmp_{ik+(jk-1)*3} = A[i+(ik-2),j+(jk-2)] * kern_ik_jk + tmp_{ik+(jk-1)*3-1}
            out[i,j] = tmp_9
        end
        out
    end
    function avx2dunrolled2x2!(out::AbstractMatrix, A::AbstractMatrix, kern::SizedOffsetMatrix{T,-1,1,-1,1}) where {T}
        rng1,  rng2  = axes(out)
        Base.Cartesian.@nexprs 3 jk -> Base.Cartesian.@nexprs 3 ik -> kern_ik_jk = kern[ik-2,jk-2]
        # Manually unpack the OffsetArray
        @avx unroll=(2,2) for j in rng2, i in rng1
            tmp_0 = zero(eltype(out))
            Base.Cartesian.@nexprs 3 jk -> Base.Cartesian.@nexprs 3 ik -> tmp_{ik+(jk-1)*3} = A[i+(ik-2),j+(jk-2)] * kern_ik_jk + tmp_{ik+(jk-1)*3-1}
            out[i,j] = tmp_9
        end
        out
    end
    function avx2dunrolled3x3!(out::AbstractMatrix, A::AbstractMatrix, kern::SizedOffsetMatrix{T,-1,1,-1,1}) where {T}
        rng1,  rng2  = axes(out)
        Base.Cartesian.@nexprs 3 jk -> Base.Cartesian.@nexprs 3 ik -> kern_ik_jk = kern[ik-2,jk-2]
        # Manually unpack the OffsetArray
        @avx unroll=(3,3) for j in rng2, i in rng1
            tmp_0 = zero(eltype(out))
            Base.Cartesian.@nexprs 3 jk -> Base.Cartesian.@nexprs 3 ik -> tmp_{ik+(jk-1)*3} = A[i+(ik-2),j+(jk-2)] * kern_ik_jk + tmp_{ik+(jk-1)*3-1}
            out[i,j] = tmp_9
        end
        out
    end
    # uq = :(for j in rng2, i in rng1
    #         tmp_0 = zero(eltype(out))
    #         Base.Cartesian.@nexprs 3 jk -> Base.Cartesian.@nexprs 3 ik -> tmp_{ik+(jk-1)*3} = A[i+(ik-2),j+(jk-2)] * kern_ik_jk + tmp_{ik+(jk-1)*3-1}
    #         out[i,j] = tmp_9
    #        end);
    # lsuq = LoopVectorization.LoopSet(macroexpand(Base, uq));
    # LoopVectorization.choose_order(lsuq)
    # lsuq

    # using LoopVectorization, OffsetArrays
    # T = Float64
    # A = rand(T, 100, 100);
    # kern = OffsetArray(rand(T, 3, 3), -1:1, -1:1);
    # out = OffsetArray(similar(A, size(A).-2), 1, 1);   # stay away from the edges of A
    # lsgeneric = LoopVectorization.@avx_debug for I in CartesianIndices(out)
    #        tmp = zero(eltype(out))
    #        for J in CartesianIndices(kern)
    #            tmp += A[I+J]*kern[J]
    #        end
    #        out[I] = tmp
    #    end;
    # LoopVectorization.choose_order(lsgeneric)
    # # out = out1;
#     lsgenerics = LoopVectorization.@avx_debug for I in CartesianIndices(out)
#            tmp = zero(eltype(out))
#            for J in CartesianIndices(skern)
#                tmp += A[I+J]*kern[J]
#            end
#            out[I] = tmp
#        end;
#     LoopVectorization.choose_order(lsgenerics)
# @macroexpand @avx for I in R
#            tmp = z
#            for J in Rk
#                tmp += A[I+J]*kern[J]
#            end
#            out[I] = tmp
#        end

    function avxgeneric!(out, A, kern, R=CartesianIndices(out), z=zero(eltype(out)))
       Rk = CartesianIndices(kern)
       @avx for I in R
           tmp = z
           for J in Rk
               tmp += A[I+J]*kern[J]
           end
           out[I] = tmp
       end
       out
    end
    function avxgeneric2!(out, A, kern)
       @avx for I in CartesianIndices(out)
           tmp = zero(eltype(out))
           for J in CartesianIndices(kern)
               tmp += A[I+J]*kern[J]
           end
           out[I] = tmp
       end
       out
   end

    for T ∈ (Float32, Float64)
        @show T, @__LINE__
        A = rand(T, 100, 100); At = copy(A');
        kern = OffsetArray(rand(T, 3, 3), -1:1, -1:1);
        out1 = OffsetArray(view(similar(A, size(A) .+ 32), (1:98) .+ 32, (1:98) .+ 32), 1, 1);   # stay away from the edges of A
        # out1 = OffsetArray(similar(A, size(A).-2), 1, 1);   # stay away from the edges of A
        out2 = similar(out1); out3 = similar(out1); out4 = similar(out1);
        skern = SizedOffsetMatrix{T,-1,1,-1,1}(parent(kern));

        old2d!(out1, A, kern);
        avx2d!(out2, A, kern);
        @test out1 ≈ out2

        avx2douter!(out3, A, kern);
        @test out1 ≈ out3

        fill!(out2, NaN); avx2d!(out2, A, skern);
        @test out1 ≈ out2

        fill!(out2, NaN); avx2douter!(out2, At', kern);
        @test out1 ≈ out2

        fill!(out2, NaN); avx2douter!(out2', A, kern);
        @test out1 ≈ out2'

        fill!(out2, NaN); avx2douter!(out2', At', kern);
        @test out1 ≈ out2'

        fill!(out3, NaN); avx2douter!(out3, A, skern);
        @test out1 ≈ out3

        fill!(out3, NaN); avx2dunrolled!(out3, A, skern);
        @test out1 ≈ out3

        fill!(out3, NaN); avx2dunrolled2x2!(out3, A, skern);
        @test out1 ≈ out3

        fill!(out3, NaN); avx2dunrolled3x3!(out3, A, skern);
        @test out1 ≈ out3

        @test avxgeneric!(out4, A, kern) ≈ out1
        fill!(out4, NaN);
        @test avxgeneric!(out4, A, skern) ≈ out1

        fill!(out4, NaN); @test avxgeneric2!(out4, A, kern) ≈ out1
        fill!(out4, NaN); @test avxgeneric2!(out4, A, skern) ≈ out1
        fill!(out4, NaN); @test avxgeneric2!(out4, At', kern) ≈ out1
        fill!(out4, NaN); @test avxgeneric2!(out4, At', skern) ≈ out1
        fill!(out4, NaN); @test avxgeneric2!(out4', A, kern) ≈ out1
        fill!(out4, NaN); @test avxgeneric2!(out4', A, skern) ≈ out1
        fill!(out4, NaN); @test avxgeneric2!(out4', At', kern) ≈ out1
        fill!(out4, NaN); @test avxgeneric2!(out4', At', skern) ≈ out1
    end

end


text:
using Documenter, ReverseDiffExpressions

makedocs(;
    modules=[ReverseDiffExpressions],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/chriselrod/ReverseDiffExpressions.jl/blob/{commit}{path}#L{line}",
    sitename="ReverseDiffExpressions.jl",
    authors="Chris Elrod",
    assets=String[],
)

deploydocs(;
    repo="github.com/chriselrod/ReverseDiffExpressions.jl",
)


text:
module PoroTomo

using AWS, Dates, SeisIO, SeisIO.Nodal
using AWS: @service
@service S3
export get_porotomo_files, porotomo_stream, porotomo_download
aws_config = AWSConfig(region="us-west-2")

"""
  get_porotomo_files(startdate; enddate=endate, filetype=filetype,orientation=orientation)

Query PoroTomo file availability on AWS.

# Arguments
- `startdate::TimeType`: Starting time to query PoroTomo files.

# Keywords
- `enddate`: Ending time to query PoroTomo files.
- `filetype`: Return file type. "SEG-Y" is only currently supported.
- `orientation`: DAS orientation, e.g. "H" for horizontal and "V" for vertical.
"""
function get_porotomo_files(
    startdate::TimeType;
    enddate=nothing,
    filetype="SEG-Y",
    orientation="H",
)
    if isnothing(enddate)
        enddate = Date(startdate + Day(1))
    end
    @assert startdate <= enddate
    @assert startdate > Date(2016,3,8)
    @assert enddate < Date(2016,3,27)
    orientation = uppercase(orientation)

    # get filelist for each date
    filelist = Array{String}(undef,0)
    daterange = Date(startdate):Day(1):Date(enddate)
    for d in daterange
        datestr = Dates.format(d,"yyyymmdd")
        path = "DAS/$filetype/DAS$orientation/$datestr/"
        dlist = s3_list_objects("nrel-pds-porotomo",path)
        append!(filelist,dlist)
    end

    # cull by datetime
    times = Array{DateTime}(undef,length(filelist))
    for ii = 1:length(filelist)
        datestr = "20"*String(split(basename(replace(filelist[ii],".sgy"=>"")),"_")[3])
        times[ii] = DateTime(datestr,"yyyymmddHHMMSS")
    end
    ind = findall(startdate .<= times .<= enddate)
    return filelist[ind]
end

function s3_list_objects(bucket,path)
    s3dict = S3.list_objects_v2(bucket,Dict("prefix"=>path),aws_config=aws_config)
    if haskey(s3dict,"Contents")
        keycount = parse(Int,s3dict["KeyCount"])
        filelist = Array{String}(undef,keycount)
        for ii = 1:keycount
            filelist[ii] = s3dict["Contents"][ii]["Key"]
        end

        if s3dict["IsTruncated"] == "true"
            while s3dict["IsTruncated"] == "true"
                token = s3dict["NextContinuationToken"]
                s3dict = S3.list_objects_v2(
                    bucket,
                    Dict("prefix"=>path,"continuation-token"=>token),
                    aws_config=aws_config,
                )
                keycount = parse(Int,s3dict["KeyCount"])
                newfiles = Array{String}(undef,keycount)
                for ii = 1:keycount
                    newfiles[ii] = s3dict["Contents"][ii]["Key"]
                end
                append!(filelist,newfiles)
            end
        end
    else
        ErrorException("No files for bucket $bucket path $path")
    end
    return filelist
end

"""
  porotomo_stream(path)

Stream file from AWS to SeisIO.NodalData structure.

# Arguments
- `path`: Filepath in `nrel-pds-porotomo` bucket to pull from S3.

# Examples
```julia-repl
julia> N = porotomo_stream("SEG-Y/DASH/20160325/PoroTomo_iDAS16043_160325235445.sgy")
NodalData with 8721 channels (4 shown)
```
"""
function porotomo_stream(path)
    stream = S3.get_object("nrel-pds-porotomo",path,aws_config=aws_config)
    return stream_nodal_segy(path,IOBuffer(stream))
end

"""
  porotomo_download(path,outdir)

# Arguments
- `path`: Filepath in `nrel-pds-porotomo` bucket to download from S3.
- `OUTDIR::String`: The output directory on EC2 instance.
"""
function porotomo_download(path,OUTDIR)

    OUTDIR = expanduser(OUTDIR)
	outfile = joinpath(OUTDIR,path)
	filedir = dirname(outfile)
	if !isdir(filedir)
	    mkpath(filedir)
    end
    stream = S3.get_object("nrel-pds-porotomo",path,aws_config=aws_config)

    open(outfile, "w") do file
        write(file, stream)
    end
    return nothing
end

function stream_nodal_segy(fname,f)

    chans = Int64[]
    nn = "N0"
    s = "0001-01-01T00:00:00"
    t = "9999-12-31T12:59:59"

    # Preprocessing
    (d0, d1) = SeisIO.parsetimewin(s, t)
    t0 = DateTime(d0).instant.periods.value*1000 - SeisIO.dtconst
    trace_fh = Array{Int16, 1}(undef, 3)
    shorts  = getfield(SeisIO.BUF, :int16_buf)
    fhd = Dict{String,Any}()
    # ww = Array{String, 1}(undef, 0)

    # ==========================================================================
    # Read file header
    fhd["filehdr"]  = SeisIO.fastread(f, 3200)
    fhd["jobid"]    = bswap(SeisIO.fastread(f, Int32))
    fhd["lineid"]   = bswap(SeisIO.fastread(f, Int32))
    fhd["reelid"]   = bswap(SeisIO.fastread(f, Int32))
    SeisIO.fast_readbytes!(f, SeisIO.BUF.buf, 48)
    SeisIO.fillx_i16_be!(shorts, SeisIO.BUF.buf, 24, 0)
    SeisIO.fastskip(f, 240)
    for i in 25:27
      shorts[i] = read(f, Int16)
    end
    SeisIO.fastskip(f, 94)

    # check endianness; can be inconsistent; 0x0400 is a kludge
    # (stands for SEG Y rev 10.24, current is rev 2.0)
    # if (unsigned(shorts[25]) > 0x0400) || (shorts[26] < zero(Int16)) || (shorts[26] > one(Int16)) || (shorts[27] < zero(Int16))
    #   shorts[25:27] .= bswap.(shorts[25:27])
    #   push!(ww, "Inconsistent file header endianness")
    # end
    # fastskip(f, 94)

    # Process
    nh = max(zero(Int16), getindex(shorts, 27))
    fhd["exthdr"] = Array{String,1}(undef, nh)
    [fhd["exthdr"][i] = SeisIO.fastread(f, 3200) for i in 1:nh]
    for (j,i) in enumerate(String[ "ntr", "naux", "filedt", "origdt", "filenx",
                                   "orignx", "fmt", "cdpfold", "trasort", "vsum",
                                   "swst", "swen0", "swlen", "swtyp", "tapnum",
                                   "swtapst", "swtapen", "taptyp", "corrtra", "bgainrec",
                                   "amprec", "msys", "zupdn", "vibpol", "segyver",
                                   "isfixed", "n_exthdr" ])
      fhd[i] = shorts[j]
    end

    nt = getindex(shorts, 1)
    trace_fh[1] = getindex(shorts,3)
    trace_fh[2] = getindex(shorts,5)
    trace_fh[3] = getindex(shorts,7)

    # ==========================================================================
    # initialize NodalData container; set variables
    if isempty(chans)
      chans = 1:nt
    end
    fs = SeisIO.sμ / Float64(trace_fh[1])
    data = Array{Float32, 2}(undef, trace_fh[2], nt)
    S = NodalData(data, fhd, chans, t0)
    net = nn * "."
    cha = string("..O", getbandcode(fs), "0")

    # ==========================================================================
    # Read traces
    j = 0

    for i = 1:nt
      C = SeisIO.do_trace(f, false, true, 0x00, true, trace_fh)
      if i in chans
        j += 1
        S.id[j]       = string(net, lpad(i, 5, '0'), cha)
        S.name[j]     = string(C.misc["rec_no"], "_", i)
        S.fs[j]       = C.fs
        S.gain[j]     = C.gain
        S.misc[j]     = C.misc
        S.t[j]        = C.t
        S.data[:, j] .= C.x
      end
    end

    # TO DO: actually use s, t here

    # ==========================================================================
    # Cleanup
    close(f)
    resize!(SeisIO.BUF.buf, 65535)
    fill!(S.fs, fs)
    fill!(S.src, fname)
    fill!(S.units, "m/m")

    # Output warnings to STDOUT
    # if !isempty(ww)
    #   for i = 1:length(ww)
    #     @warn(ww[i])
    #   end
    # end
    # S.info["Warnings"] = ww
    return S
end

end