text:
import Base: eltype, getindex, isvalid, ndims, show, size

#### Interpolation of evenly-spaced data ####

# This implements "generalized interpolation." See, for example,
# P. Thevenaz, T. Blu, and M. Unser (2000). "Interpolation Revisited."
# IEEE Transactions on Medical Imaging, 19: 739-758.

# This type manages temporary storage needed for efficient
# interpolation on a grid, so that once created it is possible to
# calculate many interpolated values without additional
# temporaries. It is designed to support both single- and multi-valued
# interpolation, and computation of the gradient with respect to the
# position of evaluation.
# An example of multi-valued interpolation is sub-pixel interpolation
# within an RGB image, and getting the 3 color values out.
type InterpGridCoefs{T<:FloatingPoint, IT<:InterpType}
    coord1d::Vector{Vector{Int}}  # for 1-d positions
    coef1d::Vector{Vector{T}}     # for 1-d coefficients
    gcoef1d::Vector{Vector{T}}    # for 1-d coefficients of the gradient
    c1d::Vector{Vector{T}}        # temp space for switch btw value & gradient
    dims::Vector{Int}
    strides::Vector{Int}
    offset::Vector{Int}  # the default N-d offset when there is no wrapping
    offset_base::Int
    index::Vector{Int}   # the alternative N-d index, needed when wrapping
    wrap::Bool           # which of these two indices to use
    valid::Bool          # some BCs return NaN when wrapping
    coef::Vector{T}      # the N-d coefficient of each neighbor (val & gradient)
end

#### InterpGrid ####
# This is the high-level interface. Use this for scalar-valued
# interpolation, and for calculation of the gradient with respect to
# the location of the evaluation point.
# For efficient multi-valued interpolation, use the low-level
# interface below. (getindex and valgrad provide examples of how to use the
# low-level interface.)
type InterpGrid{T<:FloatingPoint, N, BC<:BoundaryCondition, IT<:InterpType} <: AbstractArray{T,N}
    coefs::Array{T,N}
    ic::InterpGridCoefs{T, IT}
    x::Vector{T}
    fillval::T  # used only for BCfill (if ever)
end
function InterpGrid{T<:FloatingPoint, BC<:BoundaryCondition, IT<:InterpType}(A::Array{T}, ::Type{BC}, ::Type{IT})
    if BC == BCfill
        error("Construct BCfill InterpGrids by supplying the fill value")
    end
    coefs = copy(A)
    interp_invert!(coefs, BC, IT, 1:ndims(A))
    ic = InterpGridCoefs(coefs, IT)
    x = zeros(T, ndims(A))
    InterpGrid{T, ndims(A), BC, IT}(coefs, ic, x, nan(T))
end
function InterpGrid{T<:FloatingPoint, IT<:InterpType}(A::Array{T}, f::Number, ::Type{IT})
    coefs = pad1(A, f, 1:ndims(A))
    interp_invert!(coefs, BCnearest, IT, 1:ndims(A))
    ic = InterpGridCoefs(coefs, IT)
    x = zeros(T, ndims(A))
    InterpGrid{T, ndims(A), BCfill, IT}(coefs, ic, x, convert(T, f))
end

setx{T}(G::InterpGrid{T,1}, x::Real) = G.x[1] = x
function setx{T}(G::InterpGrid{T,2}, x::Real, y::Real)
    xG = G.x
    xG[1] = x
    xG[2] = y
end
function setx{T}(G::InterpGrid{T,3}, x::Real, y::Real, z::Real)
    xG = G.x
    xG[1] = x
    xG[2] = y
    xG[3] = z
end
function setx{T,N}(G::InterpGrid{T,N}, x::Real...)
    if length(x) != N
        error("Incorrect number of dimensions supplied")
    end
    xG = G.x
    for idim = 1:N
        xG[idim] = x[idim]
    end
end
# a version that corrects for the padding of BCfill types
setx{T}(G::InterpGrid{T,1,BCfill}, x::Real) = G.x[1] = x+1
function setx{T}(G::InterpGrid{T,2,BCfill}, x::Real, y::Real)
    xG = G.x
    xG[1] = x+1
    xG[2] = y+1
end
function setx{T}(G::InterpGrid{T,3,BCfill}, x::Real, y::Real, z::Real)
    xG = G.x
    xG[1] = x+1
    xG[2] = y+1
    xG[3] = z+1
end
function setx{T,N}(G::InterpGrid{T,N,BCfill}, x::Real...)
    if length(x) != N
        error("Incorrect number of dimensions supplied")
    end
    xG = G.x
    for idim = 1:N
        xG[idim] = x[idim]+1
    end
end

## Evaluation at single points
function _getindex{T}(G::InterpGrid{T})
    set_position(G.ic, boundarycondition(G), false, G.x)
    interp(G.ic, G.coefs)
end
function getindex(G::InterpGrid, x::Real...)
    setx(G, x...)
    _getindex(G)
end
function _valgrad{T,N}(g::Vector{T}, G::InterpGrid{T,N})
    if length(g) != N
        error("Wrong number of components for the gradient")
    end
    ic = G.ic
    coefs = G.coefs
    set_position(ic, boundarycondition(G), true, G.x)
    val = interp(ic, coefs)
    for idim = 1:N
        set_gradient_coordinate(ic, idim)
        g[idim] = interp(ic, coefs)
    end
    return val
end
function _valgrad{T}(G::InterpGrid{T,1})
    ic = G.ic
    coefs = G.coefs
    set_position(ic, boundarycondition(G), true, G.x)
    val = interp(ic, coefs)
    set_gradient_coordinate(ic, 1)
    g = interp(ic, coefs)
    return val, g
end
function valgrad{T}(G::InterpGrid{T,1}, x::Real)
    setx(G, x)
    _valgrad(G)
end
function valgrad{T}(G::InterpGrid{T}, x::Real...)
    setx(G, x...)
    g = Array(T, length(x))
    val = _valgrad(g, G)
    return val, g
end
function valgrad{T}(g::Vector{T}, G::InterpGrid{T}, x::Real...)
    setx(G, x...)
    _valgrad(g, G)
end

## Vectorized evaluation at multiple points
function getindex{T,R<:Real}(G::InterpGrid{T,1}, x::AbstractVector{R})
    n = length(x)
    v = Array(T, n)
    for i = 1:n
        setx(G, x[i])
        v[i] = _getindex(G)
    end
    v
end
getindex{T,N,R<:Real}(G::InterpGrid{T,N}, x::AbstractVector{R}) = error("Linear indexing not supported")
function getindex{T,R<:Real}(G::InterpGrid{T,2}, x::AbstractVector{R}, y::AbstractVector{R})
    nx, ny = length(x), length(y)
    v = Array(T, nx, ny)
    for i = 1:nx
        for j = 1:ny
            setx(G, x[i], y[j])
            v[i,j] = _getindex(G)
        end
    end
end
function getindex{T,N,R<:Real}(G::InterpGrid{T,N}, x::AbstractVector{R}, xrest::AbstractVector{R}...)
    if length(xrest) != N-1
        error("Dimensionality mismatch")
    end
    nx = length(x)
    nrest = [length(y) for y in xrest]
    v = Array(T, nx, nrest...)
    for c in Counter(nrest)
        for i = 1:nx
            setx(G, x[i], ntuple(N-1, i->xrest[i][c[i]])...)  # FIXME performance?? May not matter...
            v[i,c...] = _getindex(G)
        end
    end
end

# This should work for the convert in the InterpIrregular constructor below,
# but for some reason it does not
#convert{T,S}(::Type{Array{T,1}}, r::Ranges{S}) = [convert(T, x) for x in r]

#### Non-uniform grid interpolation ####

# Currently supports only 1d, nearest-neighbor or linear
# Consequently, the internal representation may change in the future
# BCperiodic and BCreflect not supported
type InterpIrregular{T<:FloatingPoint, N, BC<:BoundaryCondition, IT<:InterpType} <: AbstractArray{T,N}
    grid::Vector{Vector{T}}
    coefs::Array{T,N}
    x::Vector{T}
    fillval::T  # used only for BCfill (if ever)
end
InterpIrregular{T<:FloatingPoint, BC<:BoundaryCondition, IT<:Union(InterpNearest,InterpLinear)}(grid::Vector{T}, A::AbstractVector, ::Type{BC}, ::Type{IT}) =
    InterpIrregular(Vector{T}[grid], A, BC, IT) # special 1d syntax
InterpIrregular{T<:FloatingPoint, BC<:BoundaryCondition, IT<:Union(InterpNearest,InterpLinear)}(grid::(Vector{T}...), A::AbstractVector, ::Type{BC}, ::Type{IT}) =
    InterpIrregular(Vector{T}[grid...], A, BC, IT)
function InterpIrregular{T<:FloatingPoint, BC<:BoundaryCondition, IT<:Union(InterpNearest,InterpLinear)}(grid::Vector{Vector{T}}, A::AbstractArray, ::Type{BC}, ::Type{IT})
    if length(grid) != 1
        error("Sorry, for now only 1d is supported")
    end
    if BC == BCreflect || BC == BCperiodic
        error("Sorry, reflecting or periodic boundary conditions not yet supported")
    end
    for i = 1:length(grid)
        if !issorted(grid[i])
            error("Coordinates must be supplied in increasing order")
        end
    end
    grid = copy(grid)
    coefs = Array(T, size(A))
    copy!(coefs, A)
    x = zeros(T, ndims(A))
    InterpIrregular{T, ndims(A), BC, IT}(grid, coefs, x, nan(T))
end
function InterpIrregular{T<:FloatingPoint, IT<:InterpType}(grid, A::Array{T}, f::Number, ::Type{IT})
    iu = InterpIrregular(grid, A, BCfill, IT)
    iu.fillval = f
    iu
end

function _getindexii{T,BC<:Union(BCfill,BCna,BCnan)}(G::InterpIrregular{T,1,BC}, x::Real)
    g = G.grid[1]
    i = (x == g[1]) ? 2 : searchsortedfirst(g, x)
    (i == 1 || i == length(g)+1) ? G.fillval : _interpu(x, g, i, G.coefs, interptype(G))
end
function _getindexii{T}(G::InterpIrregular{T,1,BCnil}, x::Real)
    g = G.grid[1]
    i = (x == g[1]) ? 2 : searchsortedfirst(g, x)
    (i == 1 || i == length(g)+1) ? error(BoundsError) : _interpu(x, g, i, interptype(G))
end
# This next is necessary for precedence
getindex(G::InterpIrregular, x::Real) = _getindexii(G, x)

_interpu(x, g, i, coefs, ::Type{InterpNearest}) = (x-g[i-1] < g[i]-x) ? coefs[i-1] : coefs[i]
function _interpu(x, g, i, coefs, ::Type{InterpLinear})
    f = (x-g[i-1])/(g[i]-g[i-1])
    (1-f)*coefs[i-1] + f*coefs[i]
end

# Low-level interpolation interface

# Low-level interpolation. The process is split into two parts:
#   1. Identify the neighbors that will contribute and calculate the
#      interpolation coefficients
#   2. Use the neighbor indices and coefficients to produce
#      interpolated values
# Splitting it this way helps support efficient multi-valued
# interpolation and gradient computation.
#
# set_position handles step 1 of this process.
function set_position{T,BC<:BoundaryCondition,IT<:InterpType}(ic::InterpGridCoefs{T,IT}, ::Type{BC}, calc_grad::Bool, x::Vector{T})
    N = ndims(ic)
    if length(x) != N
        error("Dimensionality mismatch")
    end
    valid = true
    for idim = 1:N
        if !isvalid(BC, IT, x[idim], ic.dims[idim])
            valid = false
            break
        end
    end
    ic.valid = valid
    wrap = false
    if valid
        ib::Int = 0
        for idim = 1:N
            ix::Int, dx::T, twrap::Bool = interp_coords_1d(ic.coord1d[idim], BC, IT, x[idim], ic.dims[idim])
            wrap = wrap | twrap
            ib += (ix-1)*ic.strides[idim]
            interp_coefs_1d(ic.coef1d[idim], BC, IT, dx)
            if calc_grad
                interp_gcoefs_1d(ic.gcoef1d[idim], BC, IT, dx)
            end
        end
        if wrap
            interp_index_coef(ic.index, ic.coef, ic.coord1d, ic.coef1d, ic.strides)
        else
            ic.offset_base = ib
            interp_coef(ic.coef, ic.coef1d)
        end
    end
    ic.wrap = wrap
end

# Once you're done calculating interpolated _values_, call this function
# for each component of the _gradient_.
# It effectively re-does step #1 above for the gradient component,
# although it skips over the calculation of the neighbor indices. You
# then call "interp" below to evaluate the gradient, just as you to
# evaluate the value. Put each component in a loop to get all the
# gradient components.
# Naturally, if you're interpolating a multi-valued function, it's far
# better to evaluate all values, then all first components of the
# gradient, then all second components of the gradient, etc, than it
# is to recalculate the coefficients for each of the multiple values.
function set_gradient_coordinate(ic::InterpGridCoefs, i::Int)
    if !ic.valid
        return
    end
    N = ndims(ic)
    if i < 1 || i > N
        error("Wrong dimension index")
    end
    for idim = 1:N
        if idim == i
            ic.c1d[idim] = ic.gcoef1d[idim]
        else
            ic.c1d[idim] = ic.coef1d[idim]
        end
    end
    interp_coef(ic.coef, ic.c1d)# , npoints(IT), N)
end

# "interp" evaluates the interpolation. Call "set_position" first. If
# you want to also evaluate the gradient, call
# "set_gradient_coordinate" followed by "interp" for each component.
function interp{T}(ic::InterpGridCoefs{T}, A::AbstractArray, index::Int)
    if !ic.valid
        return nan(T)
    end
    coef = ic.coef
    if ic.wrap
        offset = ic.index
        index -= 1
    else
        offset = ic.offset
        index += ic.offset_base
    end
    val = zero(T)
    for i = 1:length(coef)
        c = coef[i]
        val += c == zero(T) ? zero(T) : c*convert(T,A[offset[i]+index])
    end
    return convert(T, val)
end
interp(ic::InterpGridCoefs, A::AbstractArray) = interp(ic, A, 1)

# Change the array dimensions
function set_size(ic::InterpGridCoefs, dims, strides)
    N = ndims(ic)
    if length(strides) != N
        error("Strides do not have the correct dimensionality")
    end
    if length(dims) != N
        error("Dimensions do not have the correct dimensionality")
    end
    d = [dims...]
    s = [strides...]
    for idim = 1:N
        ic.dims[idim] = d[idim]
        ic.strides[idim] = s[idim]
        interp_coords_1d(ic.coord1d[idim], interptype(ic))
    end
    interp_index(ic.offset, ic.coord1d, s)
end
function set_size(ic::InterpGridCoefs, dims)
    N = length(dims)
    s = Array(Int, N)
    s[1] = 1
    for i = 1:N-1
        s[i+1] = s[i]*dims[i]
    end
    set_size(ic, dims, s)
end

# strides is supplied separately to allow working with subarrays,
# and/or specific dimensions. For example, an RGB image might be
# interpolated with respect to the spatial dimensions but not the
# dimension that holds color.
function InterpGridCoefs{T<:FloatingPoint,IT<:InterpType}(::Type{T}, ::Type{IT}, dims::Union(Dims,Vector{Int}), strides::Union(Dims,Vector{Int}))
    N = length(strides)
    if length(dims) != N
        error("Length of dims and strides must match")
    end
    coord1d = Array(Vector{Int}, N)
    coef1d = Array(Vector{T}, N)
    gcoef1d = Array(Vector{T}, N)
    c1d = Array(Vector{T}, N)
    l = npoints(IT)
    for idim = 1:N
        coord1d[idim] = Array(Int, l)
        coef1d[idim] = Array(T, l)
        gcoef1d[idim] = Array(T, l)
        # do not allocate entries for c1d
    end
    n_coef = l^N
    offset = Array(Int, n_coef)
    index = Array(Int, n_coef)
    coef = Array(T, n_coef)
    # Pre-calculate the default offset from the strides
    for idim = 1:N
        interp_coords_1d(coord1d[idim], IT)
    end
    interp_index(offset, coord1d, strides)
    InterpGridCoefs{T,IT}(coord1d,coef1d,gcoef1d,c1d,[dims...],[strides...],offset,0,index,false,false,coef)
end
InterpGridCoefs{IT<:InterpType}(A::Array, ::Type{IT}) = InterpGridCoefs(eltype(A), IT, [size(A)...], [strides(A)...])

# Interpolation support routines

npoints(::Type{InterpNearest}) = 1
npoints(::Type{InterpLinear}) = 2
npoints(::Type{InterpQuadratic}) = 3

# Test whether a given coordinate will yield an interpolated result
isvalid{BC<:BoundaryCondition,IT<:InterpType}(::Type{BC},::Type{IT}, x, len::Int) = true
isvalid{BC<:Union(BCnil,BCnan,BCna)}(::Type{BC}, ::Type{InterpNearest}, x, len::Int) = x >= 0.5 && x-0.5 <= len
isvalid{BC<:Union(BCnil,BCnan,BCna)}(::Type{BC}, ::Type{InterpLinear}, x, len::Int) = x >= 1 && x <= len
isvalid{BC<:Union(BCnil,BCnan,BCna)}(::Type{BC}, ::Type{InterpQuadratic}, x, len::Int) = x >= 1 && x <= len

# version for offsets
function interp_coords_1d(coord1d::Vector{Int}, ::Type{InterpNearest})
    coord1d[1] = 0
end
# version for indices
function interp_coords_1d{BC<:BoundaryCondition}(coord1d::Vector{Int}, ::Type{BC}, ::Type{InterpNearest}, x, len::Int)
    ix = wrap(BC, iround(x), len)
    coord1d[1] = ix
    return ix, zero(typeof(x)), false
end

# version for offsets
function interp_coords_1d(coord1d::Vector{Int}, ::Type{InterpLinear})
    coord1d[1] = 0
    coord1d[2] = 1
end
# version for indices
function interp_coords_1d{T,BC<:BoundaryCondition}(coord1d::Vector{Int}, ::Type{BC}, ::Type{InterpLinear}, x::T, len::Int)
    ix::Int = wrap(BC, ifloor(x), len)
    dx::T = x-trunc(x)
    coord1d[1] = ix
    iswrap::Bool = (ix == len && dx > 0)
    coord1d[2] = wrap(BC, ix+1, len)
    return ix, dx, iswrap
end
function interp_coords_1d{T}(coord1d::Vector{Int}, ::Type{BCreflect}, ::Type{InterpLinear}, x::T, len::Int)
    ix = mod(ifloor(x)-1, 2*len)
    dx = x-trunc(x)
    if ix < len
        ix += 1
        ixp = ix+1
        if ixp > len
            ixp = len
        end
    else
        ix = 2*len-ix-1
        dx = one(T)-dx
        ixp = ix+1
        if ix == 0
            ix = 1
        end
    end
    coord1d[1] = ix
    coord1d[2] = wrap(BCreflect, ixp, len)
    iswrap = (ix == ixp && dx > 0)
    return ix, dx, iswrap
end

# version for offsets
function interp_coords_1d(coord1d::Vector{Int}, ::Type{InterpQuadratic})
    coord1d[1] = -1
    coord1d[2] = 0
    coord1d[3] = 1
end
# versions for indices
function interp_coords_1d{BC<:BoundaryCondition}(coord1d::Vector{Int}, ::Type{BC}, ::Type{InterpQuadratic}, x, len::Int)
    ix = iround(x)
    dx = x-ix
    ix = wrap(BC, ix, len)
    coord1d[2] = ix
    iswrap = (ix == 1 || ix == len)
    if iswrap
        coord1d[1] = wrap(BC, ix-1, len)
        coord1d[3] = wrap(BC, ix+1, len)
    else
        coord1d[1] = ix-1
        coord1d[3] = ix+1
    end
    return ix, dx, iswrap
end
# for InterpQuadratic, several require special handling
# BCnil, Bnan, BCna: for 1 <= x <= 1.5, continue the quadratic centered at x=2
function interp_coords_1d{BC<:Union(BCnil,BCnan,BCna)}(coord1d::Vector{Int}, ::Type{BC}, ::Type{InterpQuadratic}, x, len::Int)
    if x > 1.5 && x+0.5 < len
        ix = iround(x)
    elseif x <= 1.5
        ix = 2
    else
        ix = len-1
    end
    dx = x-ix
    coord1d[1] = ix-1
    coord1d[2] = ix
    coord1d[3] = ix+1
    return ix, dx, false
end
# BCnearest & BCfill: for 1 <= x <= 1.5, ensure the slope tends to 0 at x=1
function interp_coords_1d{T,BC<:Union(BCnearest,BCfill)}(coord1d::Vector{Int}, ::Type{BC}, ::Type{InterpQuadratic}, x::T, len::Int)
    if x > 1.5 && x+0.5 < len
        ix = iround(x)
        coord1d[1] = ix-1
        coord1d[2] = ix
        coord1d[3] = ix+1
        iswrap = false
        dx = x-ix
    elseif x <= 1.5
        ix = one(T)
        coord1d[1] = 2
        coord1d[2] = 1
        coord1d[3] = 2
        iswrap = true
        if x < 1
            dx = zero(T)
        else
            dx = x-ix
        end
    else
        ix = convert(T, len)
        coord1d[1] = len-1
        coord1d[2] = len
        coord1d[3] = len-1
        iswrap = true
        if x > len
            dx = zero(T)
        else
            dx = x-ix
        end
    end
    return ix, dx, iswrap
end
function interp_coords_1d(coord1d::Vector{Int}, ::Type{BCreflect}, ::Type{InterpQuadratic}, x, len::Int)
    ix = iround(x)
    dx = x-ix
    ix = mod(ix-1, 2*len)
    if ix < len
        ix += 1
    else
        dx = -dx
        ix = 2*len-ix
    end
    coord1d[2] = ix
    iswrap = (ix == 1 || ix == len)
    if iswrap
        coord1d[1] = wrap(BCreflect, ix-1, len)
        coord1d[3] = wrap(BCreflect, ix+1, len)
    else
        coord1d[1] = ix-1
        coord1d[3] = ix+1
    end
    return ix, dx, iswrap
end

function interp_coefs_1d{T,BC<:BoundaryCondition}(coef1d::Vector{T}, ::Type{BC}, ::Type{InterpNearest}, dx::T)
    coef1d[1] = one(T)
end
function interp_gcoefs_1d{T,BC<:BoundaryCondition}(coef1d::Vector{T}, ::Type{BC}, ::Type{InterpQuadratic}, dx::T)
    coef1d[1] = zero(T)
end

function interp_coefs_1d{T,BC<:BoundaryCondition}(coef1d::Vector{T}, ::Type{BC}, ::Type{InterpLinear}, dx::T)
    coef1d[1] = 1.0-dx
    coef1d[2] = dx
end
function interp_gcoefs_1d{T,BC<:BoundaryCondition}(coef1d::Vector{T}, ::Type{BC}, ::Type{InterpQuadratic}, dx::T)
    coef1d[1] = -one(T)
    coef1d[2] = one(T)
end

function interp_coefs_1d{T,BC<:BoundaryCondition}(coef1d::Vector{T}, ::Type{BC}, ::Type{InterpQuadratic}, dx::T)
    coef1d[1] = (dx-0.5)^2/2
    coef1d[2] = 0.75-dx^2
    coef1d[3] = (dx+0.5)^2/2
end
function interp_gcoefs_1d{T,BC<:BoundaryCondition}(coef1d::Vector{T}, ::Type{BC}, ::Type{InterpQuadratic}, dx::T)
    coef1d[1] = dx-0.5
    coef1d[2] = -2dx
    coef1d[3] = dx+0.5
end

eltype{T, N, BC, IT}(G::InterpGrid{T, N, BC, IT}) = T
ndims{T, N, BC, IT}(G::InterpGrid{T, N, BC, IT}) = N
boundarycondition{T, N, BC, IT}(G::InterpGrid{T, N, BC, IT}) = BC
interptype{T, N, BC, IT}(G::InterpGrid{T, N, BC, IT}) = IT
size(G::InterpGrid) = size(G.coefs)
size(G::InterpGrid, i::Integer) = size(G.coefs, i)

eltype{T, N, BC, IT}(G::InterpIrregular{T, N, BC, IT}) = T
ndims{T, N, BC, IT}(G::InterpIrregular{T, N, BC, IT}) = N
boundarycondition{T, N, BC, IT}(G::InterpIrregular{T, N, BC, IT}) = BC
interptype{T, N, BC, IT}(G::InterpIrregular{T, N, BC, IT}) = IT
size(G::InterpIrregular) = size(G.coefs)
size(G::InterpIrregular, i::Integer) = size(G.coefs, i)

eltype{T,IT<:InterpType}(ic::InterpGridCoefs{T,IT}) = T
interptype{IT<:InterpType,T}(ic::InterpGridCoefs{T,IT}) = IT
ndims(ic::InterpGridCoefs) = length(ic.coord1d)
show(io::IO, ic::InterpGridCoefs) = print(io, "InterpGridCoefs{", eltype(ic), ",", interptype(ic), "}")

# Generalized interpolation of higher order than InterpLinear requires
# inversion of the interpolation operator. See, for example, the
# Thevanez citation above.
const Q3inv = [7/8 1/4 -1/8; -1/8 5/4 -1/8; -1/8 1/4 7/8] # for handling "snippets" of size 3 along each dimension (InterpQuadratic)
# This works in place. If instead it allocated the output for you, then
# calling multiple times (e.g., to apply inversions for different
# interptypes along different dimensions) would result in unnecessary
# allocations.
interp_invert!{BC<:BoundaryCondition}(A::Array, ::Type{BC}, ::Type{InterpNearest}, dimlist) = A
interp_invert!{BC<:BoundaryCondition}(A::Array, ::Type{BC}, ::Type{InterpLinear}, dimlist) = A
function interp_invert!{BC<:BoundaryCondition}(A::Array, ::Type{BC}, ::Type{InterpQuadratic}, dimlist)
    sizeA = [size(A)...]
    stridesA = [strides(A)...]
    for idim = dimlist
        n = size(A, idim)
        # Set up the tridiagonal system
        du = fill(convert(eltype(A), 1/8), n-1)
        d = fill(convert(eltype(A), 3/4), n)
        dl = copy(du)
        M = _interp_invert_matrix(BC, InterpQuadratic, dl, d, du)
        sizeA[idim] = 1  # don't iterate over the dimension we're solving on
        for cc in Counter(sizeA)
            rng = Range(coords2lin(cc, stridesA), stridesA[idim], n)
            Base.LinAlg.solve(A, rng, M, A, rng) # in-place
        end
        sizeA[idim] = size(A, idim)
    end
end
interp_invert!{BC<:Union(BoundaryCondition,Number)}(A::Array, ::Type{BC}, IT) = interp_invert!(A, BC, IT, 1:ndims(A))
interp_invert!{BC<:Union(BoundaryCondition,Number)}(A::Array, ::Type{BC}, IT, dimlist...) = interp_invert!(A, BC, IT, dimlist)

function _interp_invert_matrix{BC<:Union(BCnil,BCnan,BCna),T}(::Type{BC}, ::Type{InterpQuadratic}, dl::Vector{T}, d::Vector{T}, du::Vector{T})
    # For these, the quadratic centered on x=2 is continued down to
    # 1 rather than terminating at 1.5
    n = length(d)
    d[1] = d[n] = 9/8
    dl[n-1] = du[1] = -1/4
    MT = Tridiagonal(dl, d, du)
    # Woodbury correction to add 1/8 for row 1, col 3 and row n, col n-2
    U = zeros(T, n, 2)
    V = zeros(T, 2, n)
    C = zeros(T, 2, 2)
    C[1,1] = C[2,2] = 1/8
    U[1,1] = U[n,2] = 1
    V[1,3] = V[2,n-2] = 1
    M = Woodbury(MT, U, C, V)
end
function _interp_invert_matrix{T}(::Type{BCreflect}, ::Type{InterpQuadratic}, dl::Vector{T}, d::Vector{T}, du::Vector{T})
    n = length(d)
    d[1] += 1/8
    d[n] += 1/8
    M = Tridiagonal(dl, d, du)
end
function _interp_invert_matrix{T}(::Type{BCperiodic}, ::Type{InterpQuadratic}, dl::Vector{T}, d::Vector{T}, du::Vector{T})
    n = length(d)
    MT = Tridiagonal(dl, d, du)
    # Woodbury correction to wrap around
    U = zeros(T, n, 2)
    V = zeros(T, 2, n)
    C = zeros(T, 2, 2)
    C[1,1] = C[2,2] = 1/8
    U[1,1] = U[n,2] = 1
    V[1,n] = V[2,1] = 1
    M = Woodbury(MT, U, C, V)
end
function _interp_invert_matrix{T,BC<:Union(BCnearest,BCfill)}(::Type{BC}, ::Type{InterpQuadratic}, dl::Vector{T}, d::Vector{T}, du::Vector{T})
    n = length(d)
    du[1] += 1/8
    dl[n-1] += 1/8
    M = Tridiagonal(dl, d, du)
end

# We have 3 functions for computing indices and coefficients:
#   interp_coef: fast, computes just the coefficients
#   interp_index_coef: fast, computes both indices and coefficients
#   interp_index: slow, computes just the indices
# The latter is rarely used (presumably, only upon InterpGridCoefs
# allocation), and hence should not be performance-critical, so the
# code is kept simple.
# The motivation for this organization is that in the "generic" case
# (see below), the overhead of looping is such that it's better to
# compute the indices and coefficients at the same time. However, when
# you don't need the indices, for fixed dimensions it's a significant
# savings (roughly 2x, at least for 3d) to just compute the
# coefficients.
function interp_index(index::Vector{Int}, coord1d::Vector{Vector{Int}}, s)
    N = length(coord1d)
    l = map(length, coord1d)
    i = 1
    for c in Counter(l)
        indx = s[1]*coord1d[1][c[1]]
        for idim = 2:N
            indx += s[idim]*coord1d[idim][c[idim]]
        end
        index[i] = indx
        i += 1
    end
end
# This is the main bottleneck, so it pays handsomely to specialize on
# dimension.
function interp_index_coef{T}(indices::Array{Int}, coefs::Array{T}, coord1d::Vector{Vector{Int}}, coef1d::Vector{Vector{T}}, strides)
    n_dims = length(coef1d)
    l = length(coef1d[1])
    if n_dims == 1
        copy!(coefs, coef1d[1])
        s = strides[1]
        c = coord1d[1]
        for i = 1:l
            indices[i] = (c[i]-1)*s + 1
        end
    elseif n_dims == 2
        s1 = strides[1]
        s2 = strides[2]
        cd1 = coord1d[1]
        cd2 = coord1d[2]
        cf1 = coef1d[1]
        cf2 = coef1d[2]
        ind = 1
        for i2 = 1:l
            offset = (cd2[i2]-1)*s2
            p = cf2[i2]
            for i1 = 1:l
                indices[ind] = offset + (cd1[i1]-1)*s1 + 1
                coefs[ind] = p*cf1[i1]
                ind += 1
            end
        end
    elseif n_dims == 3
        s1 = strides[1]
        s2 = strides[2]
        s3 = strides[3]
        cd1 = coord1d[1]
        cd2 = coord1d[2]
        cd3 = coord1d[3]
        cf1 = coef1d[1]
        cf2 = coef1d[2]
        cf3 = coef1d[3]
        ind = 1
        for i3 = 1:l
            offset3 = (cd3[i3]-1)*s3
            p3 = cf3[i3]
            for i2 = 1:l
                offset2 = offset3 + (cd2[i2]-1)*s2
                p2 = p3*cf2[i2]
                for i1 = 1:l
                    indices[ind] = offset2 + (cd1[i1]-1)*s1 + 1
                    coefs[ind] = p2*cf1[i1]
                    ind += 1
                end
            end
        end
    elseif n_dims == 4
        s1 = strides[1]
        s2 = strides[2]
        s3 = strides[3]
        s4 = strides[4]
        cd1 = coord1d[1]
        cd2 = coord1d[2]
        cd3 = coord1d[3]
        cd4 = coord1d[4]
        cf1 = coef1d[1]
        cf2 = coef1d[2]
        cf3 = coef1d[3]
        cf4 = coef1d[4]
        ind = 1
        for i4 = 1:l
            offset4 = (cd4[i4]-1)*s4
            p4 = cf4[i4]
            for i3 = 1:l
                offset3 = offset4 + (cd3[i3]-1)*s3
                p3 = p4*cf3[i3]
                for i2 = 1:l
                    offset2 = offset3 + (cd2[i2]-1)*s2
                    p2 = p3*cf2[i2]
                    for i1 = 1:l
                        indices[ind] = offset2 + (cd1[i1]-1)*s1 + 1
                        coefs[ind] = p2*cf1[i1]
                        ind += 1
                    end
                end
            end
        end
    else
        interp_index_coef_generic(indices, coefs, coord1d, coef1d, strides)
    end
end
function interp_coef{T}(coefs::Array{T}, coef1d::Vector{Vector{T}})
    n_dims = length(coef1d)
    l = length(coef1d[1])
    if n_dims == 1
        copy!(coefs, coef1d[1])
    elseif n_dims == 2
        cf1 = coef1d[1]
        cf2 = coef1d[2]
        ind = 1
        for i2 = 1:l
            p = cf2[i2]
            for i1 = 1:l
                coefs[ind] = p*cf1[i1]
                ind += 1
            end
        end
    elseif n_dims == 3
        cf1 = coef1d[1]
        cf2 = coef1d[2]
        cf3 = coef1d[3]
        ind = 1
        for i3 = 1:l
            p3 = cf3[i3]
            for i2 = 1:l
                p2 = p3*cf2[i2]
                for i1 = 1:l
                    coefs[ind] = p2*cf1[i1]
                    ind += 1
                end
            end
        end
    elseif n_dims == 4
        cf1 = coef1d[1]
        cf2 = coef1d[2]
        cf3 = coef1d[3]
        cf4 = coef1d[4]
        ind = 1
        for i4 = 1:l
            p4 = cf4[i4]
            for i3 = 1:l
                p3 = p4*cf3[i3]
                for i2 = 1:l
                    p2 = p3*cf2[i2]
                    for i1 = 1:l
                        coefs[ind] = p2*cf1[i1]
                        ind += 1
                    end
                end
            end
        end
    else
        interp_coef_generic(coefs, coef1d)
    end
end
# Compared to the dimension-specialized versions, these next two will
# be slow due to having many more arrayref operations. But at least it
# avoids an extra loop compared to the naive implementation using
# Counter, so it could be a lot worse.
# Consider making something like make_arrayind_loop_nest?
function interp_index_coef_generic{T}(indices::Array{Int}, coefs::Array{T}, coord1d::Vector{Vector{Int}}, coef1d::Vector{Vector{T}}, strides)
    n_dims = length(coord1d)
    l = length(coord1d[1])
    c = fill(1, n_dims)
    offset = Array(Int, n_dims)
    p = Array(T, n_dims)
    offset[end] = (coord1d[end][1]-1)*strides[end]
    p[end] = coef1d[end][1]
    for idim = n_dims-1:-1:1
        offset[idim] = offset[idim+1] + (coord1d[idim][1]-1)*strides[idim]
        p[idim] = p[idim+1] * coef1d[idim][1]
    end
    ind = 1
    while c[n_dims] <= l
        indices[ind] = offset[1]+1
        coefs[ind] = p[1]
        ind += 1
        c[1] += 1
        thisc = c[1]
        if thisc <= l
            if n_dims > 1
                p[1] = p[2] * coef1d[1][thisc]
                offset[1] = offset[2] + (coord1d[1][thisc]-1)*strides[1]
            else
                p[1] = coef1d[1][thisc]
                offset[1] = (coord1d[1][thisc]-1)*strides[1]
            end
        else
            idim = 1
            while c[idim] > l && idim < n_dims
                c[idim] = 1
                c[idim+1] += 1
                idim += 1
            end
            if idim == n_dims
                thisc = c[idim]
                if thisc <= l
                    p[idim] = coef1d[idim][thisc]
                    offset[idim] = (coord1d[idim][thisc]-1)*strides[idim]
                    idim -= 1
                end
            end
            while idim > 0 && idim < n_dims
                thisc = c[idim]
                offset[idim] = offset[idim+1] + (coord1d[idim][thisc]-1)*strides[idim]
                p[idim] = p[idim+1] * coef1d[idim][thisc]
                idim -= 1
            end
        end
    end
end
function interp_coef_generic{T}(indices::Array{Int}, coefs::Array{T}, coord1d::Vector{Vector{Int}}, coef1d::Vector{Vector{T}}, strides)
    n_dims = length(coord1d)
    l = length(coord1d[1])
    c = ones(Int, n_dims)
    offset = Array(Int, n_dims)
    p = Array(T, n_dims)
    p[end] = coef1d[end][1]
    for idim = n_dims-1:-1:1
        p[idim] = p[idim+1] * coef1d[idim][1]
    end
    ind = 1
    while c[n_dims] <= l
        coefs[ind] = p[1]
        ind += 1
        c[1] += 1
        thisc = c[1]
        if thisc <= l
            if n_dims > 1
                p[1] = p[2] * coef1d[1][thisc]
            else
                p[1] = coef1d[1][thisc]
            end
        else
            idim = 1
            while c[idim] > l && idim < n_dims
                c[idim] = 1
                c[idim+1] += 1
                idim += 1
            end
            if idim == n_dims
                thisc = c[idim]
                if thisc <= l
                    p[idim] = coef1d[idim][thisc]
                    idim -= 1
                end
            end
            while idim > 0 && idim < n_dims
                thisc = c[idim]
                p[idim] = p[idim+1] * coef1d[idim][thisc]
                idim -= 1
            end
        end
    end
end

# Padding an array by 1 element on each boundary of chosen dimensions
# This is needed to support BCfill
function pad1_index(szt::Tuple, dimpad)
    szv = [szt...]
    for i in dimpad
        szv[i] += 2
    end
    N = length(szt)
    ind = cell(N)
    for i in 1:N
        if szv[i] > szt[i]
            ind[i] = 2:szv[i]-1
        else
            ind[i] = 1:szv[i]
        end
    end
    return szv, ind
end
pad1{BC<:BoundaryCondition}(Ain::AbstractArray, ::Type{BC}, dimpad...) = copy(Ain)  # no padding needed for most types
pad1(Ain::AbstractArray, ::Type{BCfill}, dimpad...) = error("For BCfill, specify the fill value instead of the type BCfill")
function pad1(Ain::AbstractArray, f::Number, dimpad...)
    sz, ind = pad1_index(size(Ain), dimpad)
    A = fill(convert(eltype(Ain), f), sz...)
    A[ind...] = Ain
    return A
end


text:
const SRange{T} = Union{<:UnitSRange{T},<:StepSRange{T},<:StepSRangeLen{T},<:LinSRange{T},<:OffsetSRange{T}}
const MRange{T} = Union{<:UnitMRange{T},<:StepMRange{T},<:StepMRangeLen{T},<:LinMRange{T},<:OffsetMRange{T}}

abstract type StaticRangeWrapper{T,R<:StaticRange{T}} <: AbstractRange{T} end

@pure parent(r::StaticRangeWrapper{T,R}) where {T,R<:SRange} = R()::R

@inline first(r::StaticRangeWrapper) = first(parent(r))
@inline step(r::StaticRangeWrapper) = step(parent(r))
@inline last(r::StaticRangeWrapper) = last(parent(r))

@inline firstindex(r::StaticRangeWrapper) = firstindex(parent(r))
@inline lastindex(r::StaticRangeWrapper) = lastindex(parent(r))

@inline size(r::StaticRangeWrapper) = firstindex(parent(r))
@inline length(r::StaticRangeWrapper) = firstindex(parent(r))

ext: jl


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

struct ConstDensity{B<:ParamVolumeBounds,T<:Real} <: AbstractDensity
    bounds::B
    log_value::T
end

export ConstDensity

ConstDensity(bounds::ParamVolumeBounds{T,V}, ::typeof(one))  where {T,V}=
    ConstDensity(bounds, zero(T))

ConstDensity(bounds::ParamVolumeBounds, ::typeof(normalize)) =
    ConstDensity(bounds, -log_volume(spatialvolume(bounds)))

param_bounds(density::ConstDensity) = density.bounds

nparams(density::ConstDensity) = nparams(density.bounds)

function unsafe_density_logval(
    density::ConstDensity,
    params::AbstractVector{<:Real},
    exec_context::ExecContext = ExecContext()
)
    density.log_value
end

@inline exec_capabilities(::typeof(unsafe_density_logval), density::ConstDensity, args...) =
    ExecCapabilities(0, true, 0, true)

function unsafe_density_logval!(
    r::AbstractArray{<:Real},
    density::ConstDensity,
    params::AbstractMatrix{<:Real},
    exec_context::ExecContext
)
    fill!(r, density.log_value)
end

@inline exec_capabilities(::typeof(unsafe_density_logval!), r::AbstractArray{<:Real}, density::ConstDensity, args...) =
    ExecCapabilities(0, true, 0, true)

Distributions.sampler(density::ConstDensity) = spatialvolume(param_bounds(density))


text:
export  Conn,
        execute,
        dbconnection,
        dbname,
        gettables,
        getcolumnnames

abstract type SQLStatement end

"""
`Conn` is a shorthand for `LibPQ.Connection`, the Julia type that represents a connection to a PostgreSQL database.
"""
const Conn = LibPQ.Connection

"""
```
execute(
    conn::Conn,
    sqls::SQLStatement;
    kwargs... )
```
This function executes the SQL statement `sqls` on the database that `conn` connects to. Extra keyword arguments can be passed if necessary.

```
execute( conn::Conn; kwargs... )
```
This variant generates a function that takes a SQL statement and executes it.
"""
function execute( conn::Conn, sqls::SQLStatement; kwargs... )
    tname = split( typeof(sqls) |> string, "." )[end] |> lowercase
    fn = Core.eval( @__MODULE__, tname |> Meta.parse )
    sqls |> fn( conn, kwargs... )
end

execute( conn::Conn; kwargs... ) =
    sqls::SQLStatement -> execute( conn, sqls, kwargs... )

"""
```
dbconnection(
    dbname::AbstractString;
    host::AbstractString="localhost",
    user::AbstractString="postgres",
    password::AbstractString="postgres" )
```
This function creates a connection to a PostgreSQL database with name `dbname`. The database is hosted on the server with address `host`, and connects with the given `user` and `password`.

Because the user credentials are transmitted in plaintext, this function should **ONLY** be used to connect to a local database!
"""
dbconnection( dbname::AbstractString; host::AbstractString="localhost", user::AbstractString="postgres", password::AbstractString="postgres" ) =
    string( "host=", host == "localhost" ? "127.0.0.1" : host, " dbname=$dbname user=$user password=$password" ) |> Conn

function makeunique!( v::Vector{T} ) where T <: AbstractString
    vunique = unique(v)
    length(vunique) == length(v) && return v

    for el in vunique
        inds = findall( v .== el )

        for ii in 1:(length(inds) - 1)
            v[inds[ii+1]] = string( v[inds[ii+1]], '_', ii )
        end
    end

    v
end

function makeunique!( res::LibPQ.Result )
    makeunique!(res.column_names)
    res
end

makeunique( v::Vector{T} ) where T <: AbstractString =
    deepcopy(v) |> makeunique!

"""
```
tableresult(
    res::LibPQ.Result{false};
    df::Bool=true )
```
This function converts the result `res` of a PostgreSQL query to either a `DataFrame` if `df` is `true`, and to a `NamedTuple` (using the function `columntables` from `Tables.jl`) otherwise.

```
tableresult( df::Bool=true )
```
This variant generates a function that takes a SQL statement and executes it.
"""
function tableresult( res::LibPQ.Result{false}; df::Bool=true )
    tmpres = deepcopy(res)
    makeunique!(tmpres.column_names)
    tmpres |> (df ? DataFrame : columntable)
end

tableresult( df::Bool=true ) =
    res::LibPQ.Result{false} -> tableresult( res, df=df )

"""
```
dbname( conn::Conn )
```
This function returns the name of the database that `conn` is connecting to.
"""
function dbname( conn::Conn )
    copts = filter( copt -> copt.keyword == "dbname", LibPQ.conninfo(conn) )
    copts[1].val
end

"""
```
gettables( conn::Conn )
```
This function returns a vector with the names of all the user-defined tables in the database that `conn` is connecting to. System tables are left out of this list.
"""
function gettables( conn::Conn )
    res = Select( "information_schema.tables" ) |>
        addqueryfield!( "table_name" ) |>
        setwhereclause!( SQLAnd(
            SQLeq( "table_schema", "public" |> SQLString ),
            SQLeq( "table_catalog", dbname(conn) |> SQLString ) ) ) |>
        select(conn)
    isempty(res) && return String[]
    res[:, :table_name] |> Vector{String}
end

"""
```
getcolumnnames(
    conn::Conn,
    tablename::AbstractString )
```
This function returns the names of the columns in the table with name `tablename` of the database that `conn` is connecting to. If the table is not a user-defined table or doesn't exist, this function gives a warning returns an empty vector.
"""
function getcolumnnames( conn::Conn, tablename::AbstractString )
    if tablename ∉ gettables(conn)
        @warn "Database does not contain a table '$tablename'."
        return String[]
    end

    res = Select( "information_schema.columns" ) |>
        addqueryfield!( "column_name" ) |>
        setwhereclause!( SQLAnd(
            SQLeq( "table_catalog", dbname(conn) |> SQLString ),
            SQLeq( "table_name", tablename |> SQLString ) ) ) |>
        addorderby!( "ordinal_position" ) |>
        select(conn)
    isempty(res) && return String[]
    res[:, :column_name] |> Vector{String}
end


text:
"""
    lif(t::Real, I, V; R, tau)
    lif!(t::AbstractArray{<:Real}, I::AbstractArray{<:Real}, V::AbstractArray{<:Real}; R::AbstractArray{<:Real}, tau::AbstractArray{<:Real})
    lif!(t::CuVector{<:Real}, I::CuVector{<:Real}, V::CuVector{<:Real}; R::CuVector{<:Real}, tau::CuVector{<:Real})

Evaluate a leaky integrate-and-fire neuron.
Use `CuVector` instead of `Vector` to evaluate on GPU.

# Fields
- `t`: time since last evaluation in seconds
- `I`: external current
- `V`: current membrane potential
- `R`: resistance constant
- `vreset`: reset potential
- `tau`: time constant
"""
function lif(t::Real, I, V; R, tau)
    V *= exp(-t / tau)
    V += I * (R / tau)

    return V
end
function lif!(t::AbstractArray{<:Real}, I::AbstractArray{<:Real}, V::AbstractArray{<:Real}; R::AbstractArray{<:Real}, tau::AbstractArray{<:Real})
    # account for leakage
    @. V *= exp(-t / tau)

    # apply update step
    @. V += I * (R / tau)

    return V
end
function lif!(t::CuVector{<:Real}, I::CuVector{<:Real}, V::CuVector{<:Real}; R::CuVector{<:Real}, tau::CuVector{<:Real})
    # account for leakage
    @. V *= exp(-t / tau)

    # apply update step
    @. V += I * (R / tau)

    return V
end


text:
include("test_gage_RotationAngle.jl")
include("test_gage_ProjectedRotationAngle.jl")


text:
win = GtkWindow("hello", 100, 100)
fullscreen(win)
sleep(10)
println(width(win), " ", height(win))
destroy(win)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module SuperResModels
using SparseInverseProblemsMod
using SparseInverseProblemsMod.Util
using NLopt
import SparseInverseProblemsMod: psi, dpsi, getStartingPoint, parameterBounds, phi, computeGradient
export SuperRes, DynamicSuperRes, Conv1d, DynamicConv1d, Fourier1d, DynamicFourier1d, Fourier2d, DynamicFourier2d, is_in_bounds, dim,setBounds,psi, Conv2d, Conv2dParameters, DynamicConv2d, DynamicConv2dParameters

include("SuperRes.jl")
include("DynamicSuperRes.jl")
include("Conv1d.jl")
include("Conv2d.jl")
include("Fourier1d.jl")
include("Fourier2d.jl")
end


text:
include("Constants.jl")
using LightXML, Printf

cf = Dict("M"=>1,"CM"=>0.01,"KJ/MOL"=>1000)
convert2si(u::String) = cf[strip(uppercase(u))]

"""
Get the index of species in the gasphase limits
#   Usage
    get_index(sp,ig)
-   'sp::String' : species name
-   'ig::Gasphase' : gasphase object
"""
function get_index(sp::String,splist::Array{String,1})
    if uppercase(strip(sp)) in splist
        return findfirst(x-> x == uppercase(strip(sp)),splist)
    else
        throw(ErrorException("species $sp not present in $splist"))
    end
end

"""
Function to parse the composition string
#   Usage:
    parse_composition(composition_data,splist)
-   composition_data::String : composition specification String
-   sp_list:Array{String,1}: Array of species names
"""
function parse_composition(composition_data::String, sp_list::Array{String,1})
    species_composition = split(composition_data,",")
    composition = zeros(length(sp_list))
    for sc in species_composition
        species,val = split(sc,"=")
        composition[get_index(String(species), sp_list)] = parse(Float64,val)
    end
    return composition
end

"""
Function to convert massfractions to mole fractions
#   Usage:
    massfrac_to_molefrac!(massfracs,mol_wt,molefracs)
-   massfracs::Array{Float64,1} : Vector of massfractions
-   mol_wt::Array{Float64,1} : Vector of molecular weights
-   molefracs::Array{Float64,1} : Vector of molefractions
"""
function massfrac_to_molefrac!(molefracs::T, massfracs::T, mol_wt::T) where T <: Array{Float64,1}
    moles = massfracs ./ mol_wt
    molefracs .= moles/sum(moles)
end

"""
Function to convert massfractions to mole fractions
#   Usage:
    massfrac_to_molefrac(massfracs,mol_wt,molefracs)
-   massfracs::Array{Float64,1} : Vector of massfractions
-   mol_wt::Array{Float64,1} : Vector of molecular weights
"""
function massfrac_to_molefrac(massfracs::T,mol_wt::T) where T <: Array{Float64,1}
    moles = massfracs ./ mol_wt
    return moles/sum(moles)
end

"""
Function to convert molefractions to massfractions
#   Usage:
    molefrac_to_massfrac!(massfracs,molefracs,mol_wt)
-   molefracs::Array{Float64,1} : Vector of molefractions
-   mol_wt::Array{Float64,1} : Vector of molecular weights
"""
function molefrac_to_massfrac!(massfracs::T, molefracs::T,mol_wt::T) where T <: Array{Float64,1}
    wt = molefracs .* mol_wt
    massfracs .= wt/sum(wt)
end

"""
Average molecular weight calculation
"""
average_molwt(molfracs::T,molwt::T) where T = sum(molfracs .* molwt)
"""
density calculation
"""
density(avg_molwt::T1, T::T1, p::T1) where T1 = p*avg_molwt/(R*T)
density(molefrac::Array{T1,1},molwt::Array{T1,1}, T::T1, p::T1) where T1 = p*average_molwt(molefrac,molwt)/(R*T)

"""
Function to get floating point value from an  xml tag
#   Usage:
    get_value_from_xml(xmlroot,tag)
-   xmlroot::XMLElement : The xml root node
-   tag::String : name of the tag
The function returns the floating point value contained in the tag
"""
function get_value_from_xml(xmlroot::XMLElement,tag::String)
    tag_node = get_elements_by_tagname(xmlroot,tag)
    if is_elementnode(tag_node[1])
        tag_content = String(content(tag_node[1]))
        return parse(Float64,tag_content)
    else
        throw(error("Unable to find the node $tag in the xmlroot\n"))
    end
end

"""
Function to get the file path from an xml tag
#   Usage:
    get_text_from_xml(xmlroot,tag)
-   xmlroot::XMLElement : The xml root node
-   tag::String : name of the tag
The function returns string of file path
"""
function get_text_from_xml(xmlroot::XMLElement,tag::String)
    tag_node = get_elements_by_tagname(xmlroot,tag)
    if is_elementnode(tag_node[1])
        return String(strip(content(tag_node[1])))
    else
        throw(error("Unable to find the node $tag in the xmlroot\n"))
    end
end

"""
Function to get vector of strings contained in a tag.
Ideal for getting the species names
#   Usage:
    get_collection_from_xml(xmlroot,tag)
-   xmlroot::XMLElement : The xml root node
-   tag::String : name of the tag
The function returns vector of strings.

"""
function get_collection_from_xml(xmlroot::XMLElement,tag::String)
    tag_node = get_elements_by_tagname(xmlroot,tag)
    if is_elementnode(tag_node[1])
        return map(x->String(uppercase(x)),split(content(tag_node[1])))
    else
        throw(error("Unable to find the node $tag in the xmlroot\n"))
    end
end

"""
Function to get the mole fractions of species
#   Usage:
    get_molefraction_from_xml(xmlroot,molwt,gasphase)
-   xmlroot::XMLElement : The xml root node
-   molwt::Array{Float64} :  Vector of molecular weights
"""
function get_molefraction_from_xml(xmlroot::XMLElement, molwt::Array{Float64},gasphase::Array{String,1})
    tag_molefraction = get_elements_by_tagname(xmlroot,"molefractions")
    tag_massfraction = get_elements_by_tagname(xmlroot,"massfractions")
    local mass_fracs = Array{Float64,1}()
    local mole_fracs = Array{Float64,1}()
    if length(tag_massfraction) > 0
        massfraction = content(tag_massfraction[1])
        mass_fracs = parse_composition(massfraction,gasphase)
        mole_fracs = massfrac_to_molefrac(mass_fracs,molwt)
        return mole_fracs
    end

    if length(tag_molefraction) > 0
        molefraction = content(tag_molefraction[1])
        mole_fracs = parse_composition(molefraction,gasphase)
        return mole_fracs
    end
end

function create_header(file_stream, args...)
    for i in eachindex(args)
        for k in args[i]
            @printf(file_stream,"%10s\t",k)
        end
    end
    @printf(file_stream,"\n")
end

function write_to_file(file_stream, args...)
    for i in eachindex(args)
        for k in args[i]
            @printf(file_stream,"%.4e\t",k)
        end
    end
    @printf(file_stream,"\n")
end


text:
using Reduce

# @testset "simplify 1D L=2 Spinless Nv = 1" begin
#     L = 4
#     fos = FOS()

#     for i = 1:L
#         push!(fos.fos,PFO(FO(:α,(i%L)+1,Sl,false),symbols(Meta.parse("b$(i)"))))
#         push!(fos.fos,PFO(FO(:β,i,Sl,false),symbols(Meta.parse("b$(i)"))))
#     end

#     for i = 1:L
#         push!(fos.fos,PFO(FO(:α,i,Sl,true),symbols(Meta.parse("i$(i)"))))
#         push!(fos.fos,PFO(FO(:σ,i,Sl,true),symbols(Meta.parse("s$(i)"))))
#         push!(fos.fos,PFO(FO(:β,i,Sl,true),symbols(Meta.parse("j$(i)"))))
#     end

#     bra = FOBra([[FO(:α,i,Sl,true) for i =1:L]...,[FO(:β,i,Sl,true) for i =1:L]...])
#     ket = FOKet([[FO(:α,i,Sl,false) for i =1:L]...,[FO(:β,i,Sl,false) for i =1:L]...])

#     (expr,redundant,delta) = simplify(fos,bra,ket)
#     display(expr)
#     print("\n")
#     display(redundant)
#     print("\n")
#     display(delta)
#     print("\n")
# end

# @testset "simplify 1D L=2 Spinless Nv = 2" begin
#     L = 4
#     fos = FOS()

#     for i = 1:L
#         push!(fos.fos,PFO(FO(:α1,(i%L)+1,Sl,false),symbols(Meta.parse("b$(i)n1"))))
#         push!(fos.fos,PFO(FO(:β1,i,Sl,false),symbols(Meta.parse("b$(i)n1"))))

#         push!(fos.fos,PFO(FO(:α2,(i%L)+1,Sl,false),symbols(Meta.parse("b$(i)n2"))))
#         push!(fos.fos,PFO(FO(:β2,i,Sl,false),symbols(Meta.parse("b$(i)n2"))))
#     end

#     for i = 1:L
#         push!(fos.fos,PFO(FO(:α1,i,Sl,true),symbols(Meta.parse("i$(i)1"))))
#         push!(fos.fos,PFO(FO(:α2,i,Sl,true),symbols(Meta.parse("i$(i)2"))))
#         push!(fos.fos,PFO(FO(:σ,i,Sl,true),symbols(Meta.parse("s$(i)"))))
#         push!(fos.fos,PFO(FO(:β1,i,Sl,true),symbols(Meta.parse("j$(i)1"))))
#         push!(fos.fos,PFO(FO(:β2,i,Sl,true),symbols(Meta.parse("j$(i)2"))))
#     end

#     bra = FOBra([[FO(Meta.parse("$(s)$(nv)"),i,Sl,true) for i =1:L for s in ["α";"β"] for nv = 1:2]...])
#     ket = FOKet([[FO(Meta.parse("$(s)$(nv)"),i,Sl,false) for i =1:L for s in ["α";"β"] for nv = 1:2]...])

#     (expr,redundant,delta) = simplify(fos,bra,ket)
#     display(expr)
#     print("\n")
#     display(redundant)
#     print("\n")
#     display(delta)
#     print("\n")
# end

# @testset  "simplify 2D Lx,Ly = arbitary; Spinless; Nv = arbitary" begin
#     Lx,Ly = 2,2
#     Nv = 2
#     fos = FOS()

#     for nv = 1:Nv
#         for dir = ["r";"d"]
#             for iy = 1:Ly
#                 for ix = 1:Lx
#                     if dir == "r"
#                         ox,oy = (ix%Lx)+1,iy
#                         si,so = iy+Ly*(ix-1),oy+Ly*(ox-1)
#                         push!(fos.fos,PFO(FO(Meta.parse("l$(nv)"),si,Sl,false),symbols(Meta.parse("b$(si)_$(so)nv$(nv)"))))
#                         push!(fos.fos,PFO(FO(Meta.parse("r$(nv)"),so,Sl,false),symbols(Meta.parse("b$(si)_$(so)nv$(nv)"))))
#                     else
#                         ox,oy = ix,(iy%Ly)+1
#                         si,so = iy+Ly*(ix-1),oy+Ly*(ox-1)
#                         push!(fos.fos,PFO(FO(Meta.parse("u$(nv)"),si,Sl,false),symbols(Meta.parse("b$(si)_$(so)nv$(nv)"))))
#                         push!(fos.fos,PFO(FO(Meta.parse("d$(nv)"),so,Sl,false),symbols(Meta.parse("b$(si)_$(so)nv$(nv)"))))
#                     end
#                 end
#             end
#         end
#     end

#     for lx = 1:Lx
#         for ly = 1:Ly
#             site = ly+Ly*(lx-1)
#             for dir in ["l";"u";"r";"d"]
#                 for nv = 1:Nv
#                     push!(fos.fos,PFO(FO(Meta.parse("$(dir)$(nv)"),site,Sl,true),symbols(Meta.parse("$(dir)$(site)nv$(nv)"))))
#                 end
#             end
#             push!(fos.fos,PFO(FO(:σ,site,Sl,true),symbols(Meta.parse("s$(site)"))))
#         end
#     end

#     bra = FOBra([[FO(Meta.parse("$(s)$(nv)"),i,Sl,true) for i =1:Lx*Ly for s in ["l";"r";"d";"u"] for nv = 1:Nv]...])
#     ket = FOKet([[FO(Meta.parse("$(s)$(nv)"),i,Sl,false) for i =1:Lx*Ly for s in ["l";"r";"d";"u"] for nv = 1:Nv]...])

#     (expr,redundant,delta) = simplify(fos,bra,ket)
#     display(Reduce.Algebra.factor(expr))
#     print("\n")
#     display(redundant)
#     print("\n")
#     display(delta)
#     print("\n")
# end

@testset  "simplify 2D Lx,Ly = arbitary; Spinless; Nv = 1" begin
    Lx,Ly = 2,2
    fos = FOS()

    for dir = ["r";"d"]
        for ix = 1:Lx
            for iy = 1:Ly
                if dir == "r"
                    ox,oy = (ix%Lx)+1,iy
                    si,so = ix+Ly*(iy-1),ox+Lx*(oy-1)
                    push!(fos.fos,PFO(FO(Meta.parse("l"),so,Sl,false),symbols(Meta.parse("b$(si)$(so)"))))
                    push!(fos.fos,PFO(FO(Meta.parse("r"),si,Sl,false),symbols(Meta.parse("b$(si)$(so)"))))
                else
                    ox,oy = ix,(iy%Ly)+1
                    si,so = ix+Ly*(iy-1),ox+Lx*(oy-1)
                    push!(fos.fos,PFO(FO(Meta.parse("d"),si,Sl,false),symbols(Meta.parse("b$(si)$(so)"))))
                    push!(fos.fos,PFO(FO(Meta.parse("u"),so,Sl,false),symbols(Meta.parse("b$(si)$(so)"))))
                end
            end
        end
    end

    for iy = 1:1:Ly
        for ix = 1:1:Lx
            site = ix+Lx*(iy-1)
            for dir in ["l";"u";"r";"d"]
                push!(fos.fos,PFO(FO(Meta.parse("$(dir)"),site,Sl,true),symbols(Meta.parse("$(dir)$(site)"))))
            end
            push!(fos.fos,PFO(FO(:σ,site,Sl,true),symbols(Meta.parse("S$(site)"))))

            # push!(fos.fos,PFO(FO(:σ,site,Up,true),symbols(Meta.parse("Sup$(site)"))))
            # push!(fos.fos,PFO(FO(:σ,site,Dn,true),symbols(Meta.parse("Sdn$(site)"))))
        end
    end

    bra = FOBra([[FO(Meta.parse("$(s)"),i,Sl,true) for i =1:Lx*Ly for s in ["l";"r";"d";"u"]]...])
    ket = FOKet([[FO(Meta.parse("$(s)"),i,Sl,false) for i =1:Lx*Ly for s in ["l";"r";"d";"u"]]...])

    display(fos)
    (expr,redundant,delta) = simplify(fos,bra,ket)
    display(Reduce.Algebra.factor(expr))
    print("\n")
    display(redundant)
    print("\n")
    display(delta)
    print("\n")
end

@testset "Steps 1 by steps contraction for Lx,Ly =2,2; Nv=1 PBC" begin
    Lx,Ly = 2,2
    fos = FOS()

    push!(fos.fos,PFO(FO(Meta.parse("l"),2,Sl,false),:b12))
    push!(fos.fos,PFO(FO(Meta.parse("r"),1,Sl,false),:b12))

    push!(fos.fos,PFO(FO(Meta.parse("l"),1,Sl,false),:b21))
    push!(fos.fos,PFO(FO(Meta.parse("r"),2,Sl,false),:b21))

    for site in [1;2]
        for dir in ["l";"u";"r";"d"]
            push!(fos.fos,PFO(FO(Meta.parse("$(dir)"),site,Sl,true),symbols(Meta.parse("$(dir)$(site)"))))
        end
        push!(fos.fos,PFO(FO(:σ,site,Sl,true),symbols(Meta.parse("S$(site)"))))
    end
    bra = FOBra([[FO(Meta.parse("$(s)"),i,Sl,true) for i =1:Lx*Ly for s in ["l";"r"]]...])
    ket = FOKet([[FO(Meta.parse("$(s)"),i,Sl,false) for i =1:Lx*Ly for s in ["l";"r"]]...])

    display(fos)
    (expr,redundant,delta) = simplify(fos,bra,ket)
    display(Reduce.Algebra.factor(expr))
    print("\n")
    display(redundant)
    print("\n")
    display(delta)
    print("\n")
end

@testset "Steps 2 by steps contraction for Lx,Ly =2,2; Nv=1 PBC" begin
    Lx,Ly = 2,2
    fos = FOS()

    push!(fos.fos,PFO(FO(Meta.parse("d"),1,Sl,false),:b13))
    push!(fos.fos,PFO(FO(Meta.parse("u"),3,Sl,false),:b13))
    push!(fos.fos,PFO(FO(Meta.parse("d"),3,Sl,false),:b31))
    push!(fos.fos,PFO(FO(Meta.parse("u"),1,Sl,false),:b31))

    push!(fos.fos,PFO(FO(Meta.parse("d"),2,Sl,false),:b24))
    push!(fos.fos,PFO(FO(Meta.parse("u"),4,Sl,false),:b24))
    push!(fos.fos,PFO(FO(Meta.parse("d"),4,Sl,false),:b42))
    push!(fos.fos,PFO(FO(Meta.parse("u"),2,Sl,false),:b42))

    push!(fos.fos,PFO(FO(:u,1,Sl,true),:u1))
    push!(fos.fos,PFO(FO(:d,1,Sl,true),:d1))
    push!(fos.fos,PFO(FO(:σ,1,Sl,true),:s1))
    push!(fos.fos,PFO(FO(:u,2,Sl,true),:u2))
    push!(fos.fos,PFO(FO(:d,2,Sl,true),:d2))
    push!(fos.fos,PFO(FO(:σ,2,Sl,true),:s2))

    push!(fos.fos,PFO(FO(:u,3,Sl,true),:u3))
    push!(fos.fos,PFO(FO(:d,3,Sl,true),:d3))
    push!(fos.fos,PFO(FO(:σ,3,Sl,true),:s3))
    push!(fos.fos,PFO(FO(:u,4,Sl,true),:u4))
    push!(fos.fos,PFO(FO(:d,4,Sl,true),:d4))
    push!(fos.fos,PFO(FO(:σ,4,Sl,true),:s4))

    bra = FOBra([[FO(Meta.parse("$(s)"),i,Sl,true) for i =1:Lx*Ly for s in ["l";"r";"d";"u"]]...])
    ket = FOKet([[FO(Meta.parse("$(s)"),i,Sl,false) for i =1:Lx*Ly for s in ["l";"r";"d";"u"]]...])

    display(fos)
    (expr,redundant,delta) = simplify(fos,bra,ket)
    display(Reduce.Algebra.factor(expr))
    print("\n")
    display(redundant)
    print("\n")
    display(delta)
    print("\n")
end

# @testset  "simplify 2D Lx,Ly = 2,2; Spinless; Nv = 1; OBC" begin
#     Lx,Ly = 2,2
#     fos = FOS()

#     si,so = (1,3)
#     push!(fos.fos,PFO(FO(Meta.parse("r"),si,Sl,false),symbols(Meta.parse("b$(si)$(so)"))))
#     push!(fos.fos,PFO(FO(Meta.parse("l"),so,Sl,false),symbols(Meta.parse("b$(si)$(so)"))))

#     si,so = (2,4)
#     push!(fos.fos,PFO(FO(Meta.parse("r"),si,Sl,false),symbols(Meta.parse("b$(si)$(so)"))))
#     push!(fos.fos,PFO(FO(Meta.parse("l"),so,Sl,false),symbols(Meta.parse("b$(si)$(so)"))))

#     si,so = (1,2)
#     push!(fos.fos,PFO(FO(Meta.parse("d"),si,Sl,false),symbols(Meta.parse("b$(si)$(so)"))))
#     push!(fos.fos,PFO(FO(Meta.parse("u"),so,Sl,false),symbols(Meta.parse("b$(si)$(so)"))))

#     si,so = (3,4)
#     push!(fos.fos,PFO(FO(Meta.parse("d"),si,Sl,false),symbols(Meta.parse("b$(si)$(so)"))))
#     push!(fos.fos,PFO(FO(Meta.parse("u"),so,Sl,false),symbols(Meta.parse("b$(si)$(so)"))))

#     push!(fos.fos,PFO(FO(:σ,2,Sl,true),:s2))
#     push!(fos.fos,PFO(FO(:r,2,Sl,true),:r2))
#     push!(fos.fos,PFO(FO(:u,2,Sl,true),:u2))

#     push!(fos.fos,PFO(FO(:d,1,Sl,true),:d1))
#     push!(fos.fos,PFO(FO(:σ,1,Sl,true),:s1))
#     push!(fos.fos,PFO(FO(:r,1,Sl,true),:r1))

#     push!(fos.fos,PFO(FO(:l,3,Sl,true),:l3))
#     push!(fos.fos,PFO(FO(:d,3,Sl,true),:d3))
#     push!(fos.fos,PFO(FO(:σ,3,Sl,true),:s3))

#     push!(fos.fos,PFO(FO(:u,4,Sl,true),:u4))
#     push!(fos.fos,PFO(FO(:l,4,Sl,true),:l4))
#     push!(fos.fos,PFO(FO(:σ,4,Sl,true),:s4))

#     bra = FOBra([[FO(Meta.parse("$(s)"),i,Sl,true) for i =1:Lx*Ly for s in ["l";"r";"d";"u"]]...])
#     ket = FOKet([[FO(Meta.parse("$(s)"),i,Sl,false) for i =1:Lx*Ly for s in ["l";"r";"d";"u"]]...])

#     display(fos)
#     (expr,redundant,delta) = simplify(fos,bra,ket)
#     display(Reduce.Algebra.factor(expr))
#     print("\n")
#     display(redundant)
#     print("\n")
#     display(delta)
#     print("\n")
# end


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

using Markdown
using InteractiveUtils

# ╔═╡ 954560de-faa3-47e9-94e7-194d45f7a0df
using Manopt, Manifolds, Random

# ╔═╡ 65e3c376-ad9f-11ec-003e-7f033e6865d8
md"""
# How to Record Data during the Iterations

The recording and debug possiblities make it possible to record nearly any data during the iterations.
This tutorial illustrates how to

* record one value during the iterations
* record multiple values during the iterations and access them afterwards
* define an own `RecordAction` to perform individual recordings.

Several predefined recordings exist, for example `RecordCost()` or `RecordGradient()`, depending on the solver used.
For fields of the `Options` this can be directly done using the [`RecordEntry(:field)
For other recordings, for example more advanced computations before storing a value, an own `RecordAction` can be defined.

We illustrate these using the gradient descent from the mean computation tutorial.
"""

# ╔═╡ 4331f82f-9778-4b12-8c2d-f007df8b213b
begin
    Random.seed!(42)
    m = 30
    M = Sphere(m)
    n = 800
    σ = π / 8
    x = zeros(Float64, m + 1)
    x[2] = 1.0
    data = [exp(M, x, random_tangent(M, x, Val(:Gaussian), σ)) for i in 1:n]
end

# ╔═╡ 604b880f-faf5-4d1c-a6b9-cd94e2a8bf46
F(M, y) = sum(1 / (2 * n) * distance.(Ref(M), Ref(y), data) .^ 2)

# ╔═╡ 066267a1-06e1-4eb7-a319-1b05f25b9eac
gradF(M, y) = sum(1 / n * grad_distance.(Ref(M), data, Ref(y)))

# ╔═╡ c822fe70-becd-4b0a-b66f-34366f98dfad
md"""
## Plain examples

For the high level interfaces of the solvers, like [`gradient_descent`](https://manoptjl.org/stable/solvers/gradient_descent.html) we have to set `return_options` to `true` to obtain the whole options structure and not only the resulting minimizer.

Then we can easily use the `record=` option to add recorded values. This kesword accepts `RecordAction`s as well as several Symbols as shortcuts, for example `:Cost` to record the cost or if your options have a field `f` `:f` would record that entry.
"""

# ╔═╡ fbc02960-7f7e-4500-92ac-30df7d058fa9
R = gradient_descent(M, F, gradF, data[1]; record=:Cost, return_options=true)

# ╔═╡ 7552c7ba-6d45-4ed9-856c-b00be28a84a0
md"""
From the returned options, we see that the `Options` are encapsulated (decorated) with
`RecordOptions`.

You can attach different recorders to some operations (`:Start`. `:Stop`, `:Iteration` at time of
 writing), where `:Iteration` is the default, so the following is the same as `get_record(R, :Iteation)`. We get
"""

# ╔═╡ bd4d1bce-7942-498f-a243-bc9090bb3b1c
get_record(R)

# ╔═╡ 0545973a-d77b-4cce-b834-38a9e08f9a17
md"""
To record more than one value, you can pass a array of a mix of symbols and `RecordAction` which formally introduces `RecordGroup`. Such a Group records a tuple of values in every iteration.
"""

# ╔═╡ 68cb3e06-866d-4650-863a-61f965e3320f
R2 = gradient_descent(M, F, gradF, data[1]; record=[:Iteration, :Cost], return_options=true)

# ╔═╡ 4943d34f-0aff-4879-94a6-337ce42b2e36
md"""
Here, the Symbol `:Cost` is mapped to using the `RecordCost` action. The same holds for `:Iteration` and `:Iterate` and any member field of the current `Options`.
To access these you can first extract the group of records (that is where the `:Iteration`s are recorded – note the Pluraö) and then access the `:Cost`
"""

# ╔═╡ 9d17df05-ca92-4a98-b101-c68ba8d3f7b9
get_record_action(R2, :Iteration)

# ╔═╡ edafb0ef-be70-4fed-97b1-eabeae9be5c2
md"""
`:Iteration` is the default here, i.e. something recorded through the iterations – and we can access the recorded data the same way as we specify them in the `record=` keyword, that is, using the indexing operation.
"""

# ╔═╡ 46a79917-e03e-4c99-9b13-c525b1abb585
get_record_action(R2)[:Cost]

# ╔═╡ ef47dee4-3ee9-45c2-9ebc-2598ea56aa75
md"""This can be also done by using a the high level interface `get_record`."""

# ╔═╡ a42b9205-2447-4b64-af23-9261c34ab284
get_record(R2, :Iteration, :Cost)

# ╔═╡ 44091233-12c7-4e2b-81c4-2b2d0aa488fa
md"""
Note that the first symbol again refers to the point where we record (not to the thing we record).
We can also pass a Tuple as second argument to have our own order (not that now the second `:Iteration` refers to the recorded iteratons)
"""

# ╔═╡ 98baba1a-d64d-4852-8f80-85f7248ece98
get_record(R2, :Iteration, (:Iteration, :Cost))

# ╔═╡ de36da12-d1bb-4b91-8d5f-7a2b64f70dd1
md"""
## A more complex example

To illustrate a complicated example let's record
* the iteration number, cost and gradient field, but only every sixth iteration
* the iteration at which we stop

We first generate the problem and the options, to also illustrate the low-level works when not using `gradient_descent`.
"""

# ╔═╡ 5169f6ba-710c-4cef-8b99-9ce19300d634
p = GradientProblem(M, F, gradF)

# ╔═╡ 2015f505-cc26-4c1a-b512-3664fba1c0ed
o = GradientDescentOptions(
    M,
    copy(data[1]);
    stopping_criterion=StopAfterIteration(200) | StopWhenGradientNormLess(10.0^-9),
)

# ╔═╡ f079c326-d6bb-4710-9bd4-1133d74018d6
md"""
We first build a `RecordGroup` to group the three entries we want to record per iteration. We then put this into a `RecordEvery` to only record this every 6th iteration
"""

# ╔═╡ 3a915462-efac-4bc4-94d3-da0f400a2eaf
rI = RecordEvery(
    RecordGroup([
        :Iteration => RecordIteration(),
        :Cost => RecordCost(),
        :Gradient => RecordEntry(similar(data[1]), :gradient),
    ]),
    6,
)

# ╔═╡ cb2b3c30-d930-4de5-b919-39efc5e65fa1
md"""and a small option to record iterations"""

# ╔═╡ 930bb10c-cac0-48e2-99fa-f6faf3d0c519
sI = RecordIteration()

# ╔═╡ 995c264c-793e-4c23-9dfd-3eb7bc640e30
md"""
We now combine both into the `RecordOptions` decorator. It acts completely the same as an `Option` but records something in every iteration additionslly. This is stored in a dictionary of `RecordActions`, where `:Iteration` is the action (here the only every 6th iteration group) and the `sI` which is executed at stop.

Note that the keyword `record=` (in the high level interface `gradient_descent` only would fill the `:Iteration` symbol).
"""

# ╔═╡ 36965fe9-679d-49b8-8f0d-1cb4877cb5ab
r = RecordOptions(o, Dict(:Iteration => rI, :Stop => sI))

# ╔═╡ abc01e8a-9ae5-483c-9cda-87a196d88918
md"""We now call the solver"""

# ╔═╡ 45ee91a6-0377-46e8-bdac-faea18b120e0
res = solve(p, r)

# ╔═╡ 1c2c65f0-ccbd-4566-a8b3-d1449b683707
md"""
And we can check the recorded value at `:Stop` to see how many iterations were performed
"""

# ╔═╡ 8c6e70a2-3ebe-43a2-a96a-0219b64a1528
get_record(res, :Stop)

# ╔═╡ bba94b58-2815-4b83-8ee0-02baf9eefa77
md"""and the other values during the iterations are"""

# ╔═╡ 93af8680-a9d4-4c2d-8423-03eefb0fd752
get_record(res, :Iteration, (:Iteration, :Cost))

# ╔═╡ 5327d3b1-eb73-4725-bfa8-5d39b6993a6a
md"""
## Writing an own `RecordAction`s

Let's investigate where we want to count the number of function evaluations, again just to illustrate, since for the gradient this is just one evaluation per iteration.
We first define a cost, that counts it's own calls.
"""

# ╔═╡ b3a8e43c-6c82-483f-9449-680ee21a0650
begin
    mutable struct MyCost{T}
        data::T
        count::Int
    end
    MyCost(data::T) where {T} = MyCost{T}(data, 0)
    function (c::MyCost)(M, x)
        c.count += 1
        return sum(1 / (2 * length(c.data)) * distance.(Ref(M), Ref(x), c.data) .^ 2)
    end
end

# ╔═╡ 05284884-e05c-4714-bf53-2da071c664f7
md"""
and we define the following RecordAction, which is a Functor, i.e. a struct that is also a function. The function we have to implement is similar to a single solver step in signature, since it might get called every iteration:
"""

# ╔═╡ e0fe662d-edc5-4d1f-9fd6-987f37098cc4
begin
    mutable struct RecordCount <: RecordAction
        recorded_values::Vector{Int}
        RecordCount() = new(Vector{Int}())
    end
    function (r::RecordCount)(p::Problem, ::Options, i)
        if i > 0
            push!(r.recorded_values, p.cost.count)
        elseif i < 0 # reset if negative
            r.recorded_values = Vector{Int}()
        end
    end
end

# ╔═╡ 8eba9b9b-836b-48f3-b4fe-067bab605b87
md"""
Now we can initialize the new cost and call the gradient descent.
Note that this illustrates also the last use case – you can pass symbol-Action pairs into the `record=`array.
"""

# ╔═╡ 26cac0ff-3339-45f9-8021-04adb4d061fa
F2 = MyCost(data)

# ╔═╡ 4acd7072-f442-4433-8d98-77f772dcb35c
R3 = gradient_descent(
    M,
    F2,
    gradF,
    data[1];
    record=[:Iteration, :Count => RecordCount(), :Cost],
    return_options=true,
)

# ╔═╡ c3b78022-1055-4164-aee6-79f0580d72bd
md"""
For `:Cost` we already learned how to access them, the `:Count =>` introduces the following action to obtain the `:Count`. We can again access the whole sets of records
"""

# ╔═╡ c93b9f90-4069-4716-940b-17fcd8db055e
get_record(R3)

# ╔═╡ 8fdf9cbd-80be-4eac-9eff-2ca11f029e7a
md"""
this is equivalent to calling `R[:Iteration]`.
Note that since we introduced `:Count` we can also access a single recorded value using
"""

# ╔═╡ 4396a3c0-b9f6-4dc4-b5ae-f61cd5372671
R3[:Iteration, :Count]

# ╔═╡ 6456503a-6614-43f9-89c1-448e191c7769
md"""
And we see that the cost function is called once per iteration.
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
Manopt = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[compat]
Manifolds = "~0.7.7"
Manopt = "~0.3.20"
"""

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

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.AbstractFFTs]]
deps = ["ChainRulesCore", "LinearAlgebra"]
git-tree-sha1 = "6f1d9bc1c08f9f4a8fa92e3ea3cb50153a1b40d4"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.1.0"

[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.3"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.ArnoldiMethod]]
deps = ["LinearAlgebra", "Random", "StaticArrays"]
git-tree-sha1 = "62e51b39331de8911e4a7ff6f5aaf38a5f4cc0ae"
uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
version = "0.2.0"

[[deps.ArrayInterface]]
deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "6e8fada11bb015ecf9263f64b156f98b546918c7"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "5.0.5"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

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

[[deps.Calculus]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad"
uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
version = "0.5.1"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.14.0"

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

[[deps.ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "12fc73e5e0af68ad3137b886e3f7c1eacfca2640"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.17.1"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "96b0bc6c52df76506efc8a441c6cf1adcb1babc4"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.42.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.CovarianceEstimation]]
deps = ["LinearAlgebra", "Statistics", "StatsBase"]
git-tree-sha1 = "a3e070133acab996660d31dcf479ea42849e368f"
uuid = "587fd27a-f159-11e8-2dae-1979310e6154"
version = "0.2.7"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

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

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.DensityInterface]]
deps = ["InverseFunctions", "Test"]
git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
version = "0.4.0"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Distributions]]
deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
git-tree-sha1 = "c43e992f186abaf9965cc45e372f4693b7754b22"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.52"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.DualNumbers]]
deps = ["Calculus", "NaNMath", "SpecialFunctions"]
git-tree-sha1 = "5837a837389fccf076445fce071c8ddaea35a566"
uuid = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
version = "0.6.8"

[[deps.Einsum]]
deps = ["Compat"]
git-tree-sha1 = "4a6b3eee0161c89700b6c1949feae8b851da5494"
uuid = "b7d42ee7-0b51-5a75-98ca-779d3107e4c0"
version = "0.4.1"

[[deps.FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "246621d23d1f43e3b9c368bf3b72b2331a27c286"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.13.2"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Graphs]]
deps = ["ArnoldiMethod", "Compat", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
git-tree-sha1 = "57c021de207e234108a6f1454003120a1bf350c4"
uuid = "86223c79-3864-5bf0-83f7-82e725a168b6"
version = "1.6.0"

[[deps.HybridArrays]]
deps = ["LinearAlgebra", "Requires", "StaticArrays"]
git-tree-sha1 = "eb6b23460f5544c5d09efae0818b86736cefcd3d"
uuid = "1baab800-613f-4b0a-84e4-9cd3431bfbb9"
version = "0.4.10"

[[deps.HypergeometricFunctions]]
deps = ["DualNumbers", "LinearAlgebra", "SpecialFunctions", "Test"]
git-tree-sha1 = "65e4589030ef3c44d3b90bdc5aac462b4bb05567"
uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
version = "0.3.8"

[[deps.IfElse]]
git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.1"

[[deps.Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

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

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "91b5dcf362c5add98049e6c29ee756910b03051d"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.3"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

[[deps.Kronecker]]
deps = ["LinearAlgebra", "NamedDims", "SparseArrays", "StatsBase"]
git-tree-sha1 = "a51f46415c844dee694cb8b20a3fcbe6dba342c2"
uuid = "2c470bb0-bcc8-11e8-3dad-c9649493f05e"
version = "0.5.1"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "58f25e56b706f95125dcb796f39e1fb01d913a71"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.10"

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

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[deps.Manifolds]]
deps = ["Colors", "Distributions", "Einsum", "Graphs", "HybridArrays", "Kronecker", "LinearAlgebra", "ManifoldsBase", "Markdown", "Random", "RecipesBase", "RecursiveArrayTools", "Requires", "SimpleWeightedGraphs", "SpecialFunctions", "StaticArrays", "Statistics", "StatsBase"]
git-tree-sha1 = "e8fb7c69d6e67d377152fd882a20334535db050f"
uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
version = "0.7.7"

[[deps.ManifoldsBase]]
deps = ["LinearAlgebra", "Markdown"]
git-tree-sha1 = "058ba95cf4a41d4c7b88879f5b961352880ec919"
uuid = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb"
version = "0.12.12"

[[deps.Manopt]]
deps = ["ColorSchemes", "ColorTypes", "Colors", "DataStructures", "Dates", "LinearAlgebra", "ManifoldsBase", "Markdown", "Printf", "Random", "Requires", "SparseArrays", "StaticArrays", "Statistics", "Test"]
git-tree-sha1 = "bcd147170ad2699518155ba9b12ad33f835c5a4d"
uuid = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
version = "0.3.20"

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

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

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

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NaNMath]]
git-tree-sha1 = "737a5957f387b17e74d4ad2f440eb330b39a62c5"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "1.0.0"

[[deps.NamedDims]]
deps = ["AbstractFFTs", "ChainRulesCore", "CovarianceEstimation", "LinearAlgebra", "Pkg", "Requires", "Statistics"]
git-tree-sha1 = "0856b62716585eb90cc1dada226ac9eab5f69aa5"
uuid = "356022a1-0364-5f58-8944-0da4b18d706f"
version = "0.2.47"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "e8185b83b9fc56eb6456200e873ce598ebc7f262"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.7"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "d3538e7f8a790dc8903519090857ef8e1283eecd"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.5"

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

[[deps.QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.2"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

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

[[deps.RecipesBase]]
git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.1"

[[deps.RecursiveArrayTools]]
deps = ["Adapt", "ArrayInterface", "ChainRulesCore", "DocStringExtensions", "FillArrays", "LinearAlgebra", "RecipesBase", "Requires", "StaticArrays", "Statistics", "ZygoteRules"]
git-tree-sha1 = "f5dd036acee4462949cc10c55544cc2bee2545d6"
uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
version = "2.25.1"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[deps.Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

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

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.SimpleTraits]]
deps = ["InteractiveUtils", "MacroTools"]
git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
version = "0.9.4"

[[deps.SimpleWeightedGraphs]]
deps = ["Graphs", "LinearAlgebra", "Markdown", "SparseArrays", "Test"]
git-tree-sha1 = "a6f404cc44d3d3b28c793ec0eb59af709d827e4e"
uuid = "47aef6b3-ad0c-573a-a1e2-d07658019622"
version = "1.2.1"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.1.4"

[[deps.Static]]
deps = ["IfElse"]
git-tree-sha1 = "87e9954dfa33fd145694e42337bdd3d5b07021a6"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.6.0"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "4f6ec5d99a28e1a749559ef7dd518663c5eca3d5"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.4.3"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.StatsAPI]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "c3d8ba7f3fa0625b062b82853a7d5229cb728b6b"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.2.1"

[[deps.StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "8977b17906b0a1cc74ab2e3a05faa16cf08a8291"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.16"

[[deps.StatsFuns]]
deps = ["ChainRulesCore", "HypergeometricFunctions", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "25405d7016a47cf2bd6cd91e66f4de437fd54a07"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.16"

[[deps.SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

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

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.ZygoteRules]]
deps = ["MacroTools"]
git-tree-sha1 = "8c1a8e4dfacb1fd631745552c8db35d0deb09ea0"
uuid = "700de1a5-db45-46bc-99cf-38207098b444"
version = "0.2.2"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─65e3c376-ad9f-11ec-003e-7f033e6865d8
# ╠═954560de-faa3-47e9-94e7-194d45f7a0df
# ╠═4331f82f-9778-4b12-8c2d-f007df8b213b
# ╠═604b880f-faf5-4d1c-a6b9-cd94e2a8bf46
# ╠═066267a1-06e1-4eb7-a319-1b05f25b9eac
# ╟─c822fe70-becd-4b0a-b66f-34366f98dfad
# ╠═fbc02960-7f7e-4500-92ac-30df7d058fa9
# ╟─7552c7ba-6d45-4ed9-856c-b00be28a84a0
# ╠═bd4d1bce-7942-498f-a243-bc9090bb3b1c
# ╟─0545973a-d77b-4cce-b834-38a9e08f9a17
# ╠═68cb3e06-866d-4650-863a-61f965e3320f
# ╟─4943d34f-0aff-4879-94a6-337ce42b2e36
# ╠═9d17df05-ca92-4a98-b101-c68ba8d3f7b9
# ╟─edafb0ef-be70-4fed-97b1-eabeae9be5c2
# ╠═46a79917-e03e-4c99-9b13-c525b1abb585
# ╟─ef47dee4-3ee9-45c2-9ebc-2598ea56aa75
# ╠═a42b9205-2447-4b64-af23-9261c34ab284
# ╟─44091233-12c7-4e2b-81c4-2b2d0aa488fa
# ╠═98baba1a-d64d-4852-8f80-85f7248ece98
# ╟─de36da12-d1bb-4b91-8d5f-7a2b64f70dd1
# ╠═5169f6ba-710c-4cef-8b99-9ce19300d634
# ╠═2015f505-cc26-4c1a-b512-3664fba1c0ed
# ╟─f079c326-d6bb-4710-9bd4-1133d74018d6
# ╟─3a915462-efac-4bc4-94d3-da0f400a2eaf
# ╟─cb2b3c30-d930-4de5-b919-39efc5e65fa1
# ╟─930bb10c-cac0-48e2-99fa-f6faf3d0c519
# ╟─995c264c-793e-4c23-9dfd-3eb7bc640e30
# ╠═36965fe9-679d-49b8-8f0d-1cb4877cb5ab
# ╟─abc01e8a-9ae5-483c-9cda-87a196d88918
# ╠═45ee91a6-0377-46e8-bdac-faea18b120e0
# ╟─1c2c65f0-ccbd-4566-a8b3-d1449b683707
# ╠═8c6e70a2-3ebe-43a2-a96a-0219b64a1528
# ╟─bba94b58-2815-4b83-8ee0-02baf9eefa77
# ╠═93af8680-a9d4-4c2d-8423-03eefb0fd752
# ╟─5327d3b1-eb73-4725-bfa8-5d39b6993a6a
# ╠═b3a8e43c-6c82-483f-9449-680ee21a0650
# ╟─05284884-e05c-4714-bf53-2da071c664f7
# ╠═e0fe662d-edc5-4d1f-9fd6-987f37098cc4
# ╟─8eba9b9b-836b-48f3-b4fe-067bab605b87
# ╠═26cac0ff-3339-45f9-8021-04adb4d061fa
# ╠═4acd7072-f442-4433-8d98-77f772dcb35c
# ╟─c3b78022-1055-4164-aee6-79f0580d72bd
# ╠═c93b9f90-4069-4716-940b-17fcd8db055e
# ╟─8fdf9cbd-80be-4eac-9eff-2ca11f029e7a
# ╠═4396a3c0-b9f6-4dc4-b5ae-f61cd5372671
# ╟─6456503a-6614-43f9-89c1-448e191c7769
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using Dates

using JuliaDB

function filter_vcf_table(vcf_table, chromosome, start_position::Int64, end_position::Int64)

    start_time = now()

    filtered_vcf_table = filter(i -> (i.CHROM == chromosome) && (i.POS > start_position) && (i.POS < end_position), vcf_table)

    return filtered_vcf_table

    end_time = now()

    println("\nTook $(canonicalize(Dates.CompoundPeriod(end_time - start_time))).\n")

end

export filter_vcf_table


text:
using Test

using ReadStat, DataValues # for the Petersen data set

using FamaMacBeth

@testset "general" begin
    y = rand(100,1)
    X = rand(100,3)

    # if all in the same group, we should get OLS coefficients
    g = vec(ones(size(y,1),1))
    @test coef(fmb(X,y,g))==X\y

    g2 = repeat(collect(1:10),inner=10)
    resp = fmb(X,y,g2)
    @test all(nobs(resp).==10)
    @test size(coef(resp),1) == 3
    @test size(coef(resp),2) == 1 # because we take the mean!
end

@testset "petersen_dataset" begin
    # download and prcess data
    file = "https://www.kellogg.northwestern.edu/faculty/petersen/htm/papers/se/test_data.dta"
    loc = download(file)
    data = read_dta(loc)
    firmid = Float64.(get.(data.data[1]))
    year = Float64.(get.(data.data[2]))
    x = Float64.(get.(data.data[3]))
    y = Float64.(get.(data.data[4]))

    # fit
    m = fmb([ones(size(x,1),1) x], y, year)

    # compare with https://www.kellogg.northwestern.edu/faculty/petersen/htm/papers/se/test_data.htm
    @test all(round.(coef(m), digits=4) .== [.0313, 1.0356])
    @test all(round.(stderror(m), digits=4) .== [.0234, .0333])
    @test all(round.(tstat(m), digits=2) .== [1.34, 31.06])
    # @test round(adjr2(m),4)==.2078 # this is what Petersen lists, but cannot reproduce
    @test round(adjr2(m);digits=4)==.2076
end


text:
using HTTP, JSON

data = HTTP.request("GET",
"https://ocw.mit.edu/ans7870/6/6.006/s08/lecturenotes/files/t8.shakespeare.txt")
shakespeare = String(data.body)
shakespeareWords = split(shakespeare)

jsonWords = HTTP.request("GET",
"https://raw.githubusercontent.com/"*
"h-Klok/StatsWithJuliaBook/master/1_chapter/jsonCode.json")
parsedJsonDict = JSON.parse( String(jsonWords.body))

keywords = Array{String}(parsedJsonDict["words"])
numberToShow = parsedJsonDict["numToShow"]
wordCount = Dict([(x,count(w -> lowercase(w) == lowercase(x), shakespeareWords))
                  for x in keywords])

sortedWordCount = sort(collect(wordCount),by=last,rev=true)
sortedWordCount[1:numberToShow]

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
mutable struct ListNode{T}
    data :: T
    prev :: ListNode{T}
    next :: ListNode{T}

    function ListNode(::Type{T}) where T
        return new{T}()
    end

    function ListNode(data::T) where T
        return new{T}(data)
    end
end

Base.show(io::IO, node::ListNode) = print(io, string("ListNode(", node.data, ")"))

struct List{T}
    first :: ListNode{T} # First and last do not represent any real nodes, but used as placeholders to have a reference for beginning and ending of the list
    last  :: ListNode{T}

    function List(::Type{T}) where T
        first = ListNode(T)
        last  = ListNode(T)

        first.next = last
        last.prev  = first

        return new{T}(first, last)
    end
end

Base.show(io::IO, list::List) = print(io, string("List(", eltype(list), ")"))

Base.isempty(list::List)                  = list.first.next === list.last
Base.eltype(::Type{ <: List{T} }) where T = T

function Base.push!(list::List{T}, data::T) where T
    node = ListNode(T)

    node.data = data
    node.prev = list.last.prev
    node.next = list.last

    list.last.prev.next = node
    list.last.prev      = node

    return list
end

function pushnode!(list::List{T}, data::T) where T
    node = ListNode(T)

    node.data = data
    node.prev = list.last.prev
    node.next = list.last

    list.last.prev.next = node
    list.last.prev      = node

    return node
end

function Base.empty!(list::List)
    list.first.next = list.last
    list.last.prev  = list.first
    return list
end

function remove(node::ListNode)
    next = node.next
    prev = node.prev

    next.prev = prev
    prev.next = next

    return node
end

Base.IteratorSize(::Type{ <: List })   = Base.SizeUnknown()
Base.IteratorEltype(::Type{ <: List }) = Base.HasEltype()

Base.iterate(list::List)                  = isempty(list) ? nothing : (list.first.next.data, list.first.next.next)
Base.iterate(list::List, state::ListNode) = list.last === state ? nothing : (state.data, state.next)


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ ebcdbdd0-9e9f-11eb-37cd-9d8e9dd6028f
begin
	using Plots,PlutoUI
	plotlyjs()
end

# ╔═╡ 05405030-12ea-4cb4-b272-0d1552a480c6
f(x) = x / (2x+1)

# ╔═╡ dd39da8c-d2e1-4df2-85e4-f5d6d928bcba
begin
xo_s = @bind xo Scrubbable(-100.0:100.0,default=0.0)
s_s  = @bind s Scrubbable(0.01:0.01:10,default=1.0)
xf_s = @bind xf Scrubbable(-100.0:100.0,default=10.0)
md" $(xo_s) : $(s_s) : $(xf_s)"
end

# ╔═╡ ccdb4441-c1c3-4966-a9c5-bd3e886815b8
begin
r = xo:s:xf
plot(r,f)
end

# ╔═╡ Cell order:
# ╟─ebcdbdd0-9e9f-11eb-37cd-9d8e9dd6028f
# ╟─05405030-12ea-4cb4-b272-0d1552a480c6
# ╟─dd39da8c-d2e1-4df2-85e4-f5d6d928bcba
# ╟─ccdb4441-c1c3-4966-a9c5-bd3e886815b8

max_stars_repo_name: akielbowicz/cbc-analisis-ingenieria

max_issues_repo_name: akielbowicz/cbc-analisis-ingenieria

max_forks_repo_name: akielbowicz/cbc-analisis-ingenieria


text:
ts = DateTime(2018,01,13):Minute(1):DateTime(2018,01,14)
x = randn(length(ts))

daylight(ts,x,loc=Point(42.35,-71.05),tz=tz"EST")