text:
const Dim = Union{Int,Base.Slice}

"""
    BufferType = Union{Vector,BitVector}

Types which can be a buffer.
"""
const BufferType = Union{Vector,BitVector}

"""
    isresizable(A::AbstractArray)

Determines if an array is resizable, if not `resize!(A, args...)` would throw an error.
"""
isresizable(A::AbstractArray) = isresizable(typeof(A))
isresizable(::Type{T}) where {T} = _isresizable(has_parent(T), parent_type(T))
_isresizable(::True, ::Type{T}) where {T} = isresizable(T)
_isresizable(::True, ::Type{<:BufferType}) = true
_isresizable(::False, ::Type) = false
_isresizable(::False, ::Type{<:BufferType}) = true # BufferType is resizable now

"""
    has_resize_buffer(A::AbstractArray)

Determines if an array has `resize_buffer*` methods,
if not `resize!(A, args...)` would resize its parent.
"""
has_resize_buffer(A::AbstractArray) = has_resize_buffer(typeof(A))
has_resize_buffer(::Type{T}) where {T} = parent_type(T) <: BufferType

"""
    to_parentinds(A::AbstractArray, Is::Tuple) -> Is′

Convert the index(s) `Is` of `A` to index(s) `Is′` of `parent(A)`.
"""
to_parentinds(::AbstractArray, Is::Tuple) = Is
"""
    to_parentinds(A::AbstractArray, i::Integer, I) -> (i′, I′)

Convert the index(s) `I` at `d`th dimension of `A` to index(s) `I′` at `d′`th
dimension of `parent(A)`.
"""
to_parentinds(::AbstractArray, d::Integer, I) = d, I
# impl of adjoint and transpose
to_parentinds(A::AdjOrTransAbsVec, (i, J)::Tuple) = (check_dimbounds(A, 1, i); (J,))
to_parentinds(::AdjOrTransAbsMat, (I, J)::Tuple) = (J, I)
to_parentinds(::AdjOrTrans, d::Integer, I) = ifelse(d == 1, 2, 1), I

"""
    copyto_parent!(dst::AbstractArray, src::AbstractArray, dinds...)

Copy the elements of `src` to parent of `dst` at the `dst` indices `dinds`.
This method is called by `resize_buffer!` to copy elements of the parent.
The default method is `copyto!(view(dst, dinds...), src)`.
The default method assumes that
`setindex!(parent(dst), v, i...)` is the same as `setindex!(dst, v, i...)` .
Override this method if they are different.
"""
copyto_parent!(dst::AbstractArray, src::AbstractArray, dinds...) =
    copyto!(view(dst, dinds...), src)

# Base.sizehint!(A, sz)
"""
    sizehint!(A::AbstractArray{T,N}, sz::NTuple{N}) where {T,N}

Suggest that array `A` reserve size for at least `sz`.
This can improve performance.
"""
Base.sizehint!(A::AbstractArray{T,N}, sz::NTuple{N,Any}) where {T,N} =
    sizehint!(A, to_dims(sz))
Base.sizehint!(A::AbstractArray{T,N}, sz::Dims{N}) where {T,N} = sizehint!(A, prod(sz))

"""
    Base.sizehint!(A::AbstractArray, nl::Integer)

Suggest that array `A` reserve capacity for at least `nl` elements.
This can improve performance.
"""
function Base.sizehint!(A::AbstractArray, nl::Integer)
    if isresizable(A)
        sizehint!(parent(A), nl)
        return A
    end
    throw_methoderror(sizehint!, A)
end

"""
    resize_parent!(A::AbstractArray{T,N}, sz::NTuple{N})

The same as `resize_parent!(A, prod(sz)`.
"""
@inline resize_parent!(A::AbstractArray{T,N}, sz::NTuple{N,Any}) where {T,N} =
    resize_parent!(A, prod(to_dims(sz)))

"""
    resize_parent!(A::AbstractArray, nl::Integer)

Resize the parent of A. This method will (and should only) be called by
`resize_buffer!` or `resize_buffer_dim!`, the default implementation is
`resize!(parent(A), nl)`, but for arrays with preserved space, this methods can
be override to keep size of parent.
"""
resize_parent!(A::AbstractArray, nl::Integer) = resize_parent!(A, Int(nl))
function resize_parent!(A::AbstractArray, nl::Int)
    if parent_type(A) <: BufferType
        return resize!(parent(A), nl)
    end
    error("parent_type(A) must be BufferType")
end

"""
    pre_resize!(A::AbstractArray{T,N}, sz::NTuple{N,Any})
    pre_resize!(A::AbstractArray{T,N}, d::Integer, n::Any)

Do something before resize `A` with given arguments (do nothing by default).
This method is called by `resize!` with the same arguments.
"""
pre_resize!
"""
    after_resize!(A::AbstractArray{T,N}, sz::NTuple{N,Any})
    after_resize!(A::AbstractArray{T,N}, d::Integer, n::Any)

Do something after resize `A` with given arguments (do nothing by default).
This methods is called by `resize!` with the same arguments.
"""
after_resize!

@inline pre_resize!(A::AbstractArray{T,N}, ::NTuple{N,Any}) where {T,N} = A
@inline after_resize!(A::AbstractArray{T,N}, ::NTuple{N,Any}) where {T,N} = A

"""
    resize!(A::AbstractArray{T,N}, sz)

Resize `A` to `sz`. `sz` can be a tuple of integer or Colon or iterator.
"""
function Base.resize!(A::AbstractArray{T,N}, dims::NTuple{N,Any}, ::B=False()) where {T,N,B}
    if isresizable(A)
        dims′ = _to_indices(A, dims, B())
        pre_resize!(A, dims′)
        if has_resize_buffer(A)
            resize_buffer!(A, to_indices(A, dims′)...)
        else
            resize!(parent(A), to_parentinds(A, dims′), True())
            setsize!(A, dims′)
        end
        after_resize!(A, dims′)
        return A
    end
    throw_methoderror(resize!, A)
end

"""
    resize_buffer!(A::AbstractArray, nsz...)

Implementation of `resize!(A, nsz)` where `parent(A)` is [`BufferType`](@ref).
"""
resize_buffer!(A::AbstractArray, nsz...) = _resize_buffer!(A, nsz...)
# resize vector will not move any element
function _resize_buffer!(A::AbstractVector{T}, n::Int) where {T}
    checksize(A, (n,))
    resize_parent!(A, n) # resize! buffer by Base.resize!
    setsize!(A, 1, n)
    return A
end
# resize array with Int sz
function _resize_buffer!(A::AbstractArray{T,N}, nsz::Vararg{Int,N}) where {T,N}
    checksize(A, nsz)
    sz = size(A)
    nsz == sz && return A # if sz not change
    if nsz[1:N-1] == sz[1:N-1] # if only last dim changed
        resize_parent!(A, nsz)
        setsize!(A, N, nsz[N])
        return A
    end
    # normal case
    ssz = map(min, sz, nsz)
    # step 1: get reserved elements
    src = if ssz != sz # some elements will be drop
        A[map(Base.OneTo, ssz)...]
    else # all elements will be reserved
        copy(A)
    end
    # step 2: resize A for enough capacity
    resize_parent!(A, nsz)
    setsize!(A, nsz)
    # step 3: copy src to dst region of A
    copyto_parent!(A, src, map(Base.OneTo, ssz)...)
    return A
end
# resize with Colon
function _resize_buffer!(A::AbstractArray{T,N}, dims::Vararg{Dim,N}) where {T,N}
    M = restdim(Base.Slice, dims...)
    # if only last dim changed
    if M == 1
        dims[N] > 0 || error("dimension(s) must be > 0")
        nlen = stride(A, N) * dims[N]
        resize_parent!(A, nlen)
        setsize!(A, N, dims[N])
        return A
    end
    nsz = to_dims(dims)
    checksize(A, nsz)
    sz_tail = tailn(Val(M), size(A)...)
    nsz_tail = tailn(Val(M), nsz...)
    nlen = prod(nsz)
    cpy = similar(A, nlen)
    blocklen = stride(A, N - M + 1)
    ssz_tail = map(min, sz_tail, nsz_tail)
    if ssz_tail == sz_tail
        copyto!(cpy, parent(A))
    else
        strds = _accumulate_rec(identity, *, 1, sz_tail...)
        for (i, sub) in enumerate(CartesianIndices(ssz_tail))
            ind = mapreduce((strd, i) -> strd * (i - 1), +, strds, sub.I) + 1
            copyto!(
                cpy,
                blocklen * (i - 1) + 1,
                parent(A),
                blocklen * (ind - 1) + 1,
                blocklen,
            )
        end
    end
    resize_parent!(A, nlen)
    setsize!(A, dims)
    if ssz_tail == nsz_tail
        copyto!(parent(A), cpy)
    else
        nstrds = _accumulate_rec(identity, *, 1, nsz_tail...)
        for (i, sub) in enumerate(CartesianIndices(ssz_tail))
            ind = mapreduce((strd, i) -> strd * (i - 1), +, nstrds, sub.I) + 1
            copyto!(
                parent(A),
                blocklen * (ind - 1) + 1,
                cpy,
                blocklen * (i - 1) + 1,
                blocklen,
            )
        end
    end
    return A
end
_resize_buffer!(A::AbstractArray{T,N}, ::Vararg{Base.Slice,N}) where {T,N} = A
# resize with inds
function _resize_buffer!(A::AbstractArray{T,N}, inds::Vararg{Any,N}) where {T,N}
    sz = to_dims(inds)
    nlen = prod(sz)
    sinds = map(_to_sind, axes(A), inds)
    src = A[sinds...]
    resize_parent!(A, nlen)
    setsize!(A, inds)
    dinds = map(_to_oneto, sinds)
    copyto_parent!(A, src, dinds...)
    return A
end

_to_sind(ind, n::Integer) = n >= length(ind) ? UnitRange(ind) : UnitRange(ind[1:n])
_to_sind(ind, I) = I
_to_oneto(I) = Base.OneTo(length(I))

@inline pre_resize!(A::AbstractArray, ::Int, ::Any) = A
@inline after_resize!(A::AbstractArray, ::Int, ::Any) = A

"""
    resize!(A::AbstractArray{T,N}, d::Integer, I)

Resize the `d`th dimension to `I`, where `I` can be an integer or a colon or an iterator.
"""
function Base.resize!(A::AbstractArray, d::Integer, I)
    if isresizable(A)
        d′ = Int(d)
        I′ = to_index(A, axes(A, d′), I)
        pre_resize!(A, d′, I′)
        if has_resize_buffer(A)
            resize_buffer_dim!(A, d′, I′)
        else
            resize!(parent(A), to_parentinds(A, d′, I′)...)
            setsize!(A, d′, I′)
        end
        after_resize!(A, d′, I′)
        return A
    end
    throw_methoderror(resize!, A)
end
Base.resize!(A::AbstractArray, ::Integer, ::Colon) = A

"""
    resize_buffer_dim!(A::AbstractArray, d::Int, I)

Implementation of `resize!(A, d, I)` where `parent(A)` is a [`BufferType`](@ref).
"""
resize_buffer_dim!(A::AbstractArray, d::Int, I) = _resize_buffer_dim!(A, d, I)

function _resize_buffer_dim!(A::AbstractArray, d::Int, n::Int)
    N = ndims(A)
    n == size(A, d) && return A
    if d == N
        resize_parent!(A, stride(A, N) * n)
        setsize!(A, N, n)
        return A
    end
    blk_len, blk_num, batch_num = _blkinfo(A, d)
    nlen = blk_len * n * batch_num
    cpy = similar(A, nlen)
    batch_len = blk_len * blk_num
    sbatch_len = blk_len * min(n, blk_num)
    δ = 0
    for j in 1:batch_num
        soffs = batch_len * (j - 1) + 1
        doffs = soffs + δ
        copyto!(cpy, doffs, parent(A), soffs, sbatch_len)
        δ = δ + blk_len * (n - blk_num)
    end
    resize_parent!(A, nlen)
    setsize!(A, d, n)
    copyto!(parent(A), cpy)
    return A
end
function _resize_buffer_dim!(A::AbstractArray, d::Int, I::Base.LogicalIndex)
    @boundscheck check_dimbounds(A, d, I)
    I′ = I.mask
    n = length(I)
    blk_len, blk_num, batch_num = _blkinfo(A, d)
    nlen = blk_len * n * batch_num
    cpy = similar(A, nlen)
    batch_len = blk_len * blk_num
    δ = 0
    for j in 1:batch_num
        for (k, flag) in enumerate(I′)
            if flag
                soffs = batch_len * (j - 1) + blk_len * (k - 1) + 1
                doffs = soffs + δ
                copyto!(cpy, doffs, parent(A), soffs, blk_len)
            else
                δ = δ - blk_len
            end
        end
    end
    resize_parent!(A, nlen)
    setsize!(A, d, I)
    copyto!(parent(A), cpy)
    return A
end
function _resize_buffer_dim!(A::AbstractArray, d::Int, I::AbstractVector)
    @boundscheck check_dimbounds(A, d, I)
    I′ = zeros(Bool, size(A, d))
    I′[I] .= true
    n = length(I)
    blk_len, blk_num, batch_num = _blkinfo(A, d)
    nlen = blk_len * n * batch_num
    cpy = similar(A, nlen)
    batch_len = blk_len * blk_num
    δ = 0
    for j in 1:batch_num
        for (k, flag) in enumerate(I′)
            if flag
                soffs = batch_len * (j - 1) + blk_len * (k - 1) + 1
                doffs = soffs + δ
                copyto!(cpy, doffs, parent(A), soffs, blk_len)
            else
                δ = δ - blk_len
            end
        end
    end
    resize_parent!(A, nlen)
    setsize!(A, d, I)
    copyto!(parent(A), cpy)
    return A
end

# aux methods
restdim(::Type{T}, ::T, rest...) where {T} = restdim(T, rest...)
restdim(::Type{T}, rest...) where {T} = length(rest)

function tailn(::Val{N}, item, items...) where {N}
    M = length(items)
    if N == M
        return items
    elseif N > M
        return item, items...
    elseif N < M
        return tailn(Val(N), items...)
    end
end

to_index(A, ind, I) = Base.to_index(A, I)
to_index(A, ind, I::FunctionIndices.AbstractFunctionIndex) =
    FunctionIndices.to_index(FunctionIndices.indextype(A, I), A, ind, I)

_accumulate_rec(f, op, init, item) = (init, op(init, f(item)))
_accumulate_rec(f, op, init, item, items...) =
    (init, _accumulate_rec(f, op, op(init, f(item)), items...)...)

@inline _to_indices(::AbstractArray, dims, ::True) = dims
@inline _to_indices(A::AbstractArray, dims, ::False) = to_indices(A, dims)

function _blkinfo(A::AbstractArray, i::Integer)
    i <= ndims(A) || throw(ArgumentError("dim must less than ndims(A)"))
    blk_len = 1
    sz = size(A)
    @inbounds for j in 1:(i-1)
        blk_len *= sz[j]
    end
    batch_num = 1
    @inbounds for j in (i+1):ndims(A)
        batch_num *= sz[j]
    end
    return blk_len, @inbounds(sz[i]), batch_num
end

checksize(::Type{Bool}, A, sz::Dims) = (&)(map(>(0), sz)...)
checksize(A, sz::Dims) = checksize(Bool, A, sz) || error("dimension(s) must be > 0")

# checkbounds at dim d
check_dimbounds(A::AbstractArray, d::Integer, I) =
    checkindex(Bool, axes(A, d), I) || throw_dimboundserror(A, d, I)

# Exceptions
struct MethodUndefineError <: Exception
    f::Any
    T::DataType
end
Base.showerror(io::IO, err::MethodUndefineError) =
    print(io, "MethodUndefineError: ", err.f, " is not defined for ", err.T)

throw_methoderror(f, A::AbstractArray{N}) where {N} =
    throw(MethodUndefineError(f, typeof(A)))

struct DimBoundsError <: Exception
    A::AbstractArray
    d::Int
    I::Any
end
function Base.showerror(io::IO, err::DimBoundsError)
    print(io, "DimBoundsError: attempt to access ")
    summary(io, err.A), print(io, " at dimension ", err.d, " with index ")
    if err.I isa Base.LogicalIndex
        print(io, err.I.mask)
    else
        print(io, err.I)
    end
    return nothing
end
throw_dimboundserror(A::AbstractArray{N}, d::Integer, I) where {N} =
    throw(DimBoundsError(A, Int(d), I))


text:
# Introductory demonstration of IterationControl.jl

using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

# Here's a simple iterative mdel that computes Babylonian
# approximations to a square root:

mutable struct SquareRooter
    x::Float64     # input - number to be square rooted
    root::Float64  # current approximation of root
    training_losses::Vector{Float64} # successive approximation differences
    SquareRooter(x) = new(x, 1.0, Float64[])
end

function train!(m::SquareRooter, Δn::Int)
    m.training_losses = Float64[]
    for i in 1:Δn
        next_guess = (m.root + m.x/m.root)/2
        push!(m.training_losses, abs(next_guess - m.root))
        m.root = next_guess
    end
end

loss(m::SquareRooter) = abs(m.root^2 - m.x)
training_losses(m::SquareRooter) = m.training_losses

# And here it in action:

model = SquareRooter(9)
model.root

#-

train!(model, 2) # train for 2 iterations
model.root

#-

train!(model, 1) # train for 1 more iteration
model.root

#  Then we can replace the integer argument `n` in `train!(model, n)`
#  with a number of more sophisticated *controls* by "lifting" the method
# `train!` to the `IterationControl.train!` method defined in this
# package:

using IterationControl
IterationControl.train!(model::SquareRooter, n) =  train!(model, n)

# The lifted `train!` has the same functionality as the original one:

model = SquareRooter(9)
IterationControl.train!(model, 2)

model.root

# But now we can also do this:

IterationControl.train!(model, Step(2), NumberLimit(3), Info(m->m.root));

# Here each control is repeatedly applied until one of them triggers a
# stop. The first control `Step(2)` says "train the model two more
# iterations"; the second says "stop after 3 repetitions" (of the
# sequence of control applications); and the third, "log the value of
# the root to `Info`".

# If `model` admits a method returning a loss (for example, the
# difference between `x` and the square of `root`), then we can lift
# that method to `IterationControl.loss` to enable control using
# loss-based stopping criteria, such as a loss threshold. In the
# demonstation below, we also include a callback:

model = SquareRooter(4)
train!(model, 1)
loss(model)

#-

IterationControl.loss(model::SquareRooter) = loss(model)

losses = Float64[]
callback(model) = push!(losses, loss(model))

IterationControl.train!(model,
                        Step(1),
                        Threshold(0.0001),
                        Callback(callback));
#-

losses

# If training `model` generates user-inspectable "training losses" (one
# per iteration) then similarly lifting the appropriate access function
# to `IterationControl.training_losses` enables Prechelt's
# progress-modified generalization loss stopping criterion, `PQ`.

# `PQ` is the only criterion from the
# [EarlyStopping.jl](https://github.com/ablaom/EarlyStopping.jl) package
# not otherwise enabled when `IterationControl.loss` is overloaded as
# above.

# *Reference.* [Prechelt, Lutz
#  (1998)](https://link.springer.com/chapter/10.1007%2F3-540-49430-8_3):
#  "Early Stopping - But When?", in *Neural Networks: Tricks of the
#  Trade*, ed. G. Orr, Springer.

# The interface just described is sufficient for controlling
# conventional machine learning models with an iteration parameter, as
# this [tree boosting
# example](https://github.com/ablaom/IterationControl.jl/tree/master/examples/tree_booster)
# shows.

using Literate #src
Literate.markdown(@__FILE__, @__DIR__, execute=true) #src
Literate.notebook(@__FILE__, @__DIR__, execute=true) #src

ext: jl


text:
## add lead & lag discharge to main dataframe
function get_leadlag(times::Vector{DateTime}, df::DataFrame, x::String)
   for (i,j) in enumerate(times)
      for (ind,val) in enumerate( filter(row -> j-Dates.Hour(2) .<= row.datetime .<= j+Dates.Hour(2), df)[!,Symbol(x)] )
         bel[i, Symbol(string(x, ind))] += val
      end
   end
end

## Function to impute missing value sin dataframe
function impute_var(df::DataFrame, x::String)
   for i in eachrow(df)
      for j in 1:20
         if i[string(x,j)] == 0
            i[string(x,j)] += i[string(x,j-1)]
         end
      end
   end
end

## Get site input
function get_site_input()
        println("Please enter a site  04168557 or 04165705")
        site = readline()
        if contains(site, r"04168557|04165705") == false
           println("Please input site 04168557 or 04165705")
           site = get_site_input()
        end
return(site)
end

## Get date Time input
function get_dt_input()
        println("Please enter date and time you would like to predict  YYYY-MM-DDTHH:MM:SS")
        dt = readline()
        try
          dt =  DateTime(dt)
        catch e
           println("Please use this format: YYYY-MM-DDTHH:MM:SS")
           dt = get_dt_input()
        end
return(dt)
end

## NWIS URL
function nwis_url(begin_date::DateTime, end_date::Date=Date(begin_date))
   begin_date = Date(begin_date)
   string("https://nwis.waterdata.usgs.gov/nwis/uv?cb_00060=on&cb_00065=on&cb_72255=on&format=rdb&site_no=04165710&period=&begin_date=",begin_date,"&end_date=", end_date)
end

## Get continous data from NWIS Fort Wayne Detroit
function nwis_pull(dt::DateTime)
        url = nwis_url(dt)
        cont_dis = CSV.read(download(url),DataFrame,comment="#", delim="\t")
        deleteat!(cont_dis,1)
        dropmissing!(cont_dis)
return(cont_dis)
end

## Clean Continous Discharge
function clean_dis!(cont_dis::DataFrame=cont_dis)
   rename!(cont_dis, Symbol("240039_00060") => :discharge,
                  Symbol("72269_00065") => :gage_height,
                 Symbol("242971_72255") => :channel_velocity)
   select!(cont_dis, [:datetime, :discharge, :gage_height, :channel_velocity])
   cont_dis.datetime = DateTime.(cont_dis.datetime, dateformat"yyyy-mm-dd HH:MM:SS")
   cont_dis.discharge =  parse.(Int64, cont_dis.discharge)
   for i in [:gage_height, :channel_velocity]
        cont_dis[!,i] =  parse.(Float64, cont_dis[!,i])
   end
return(cont_dis)
end

## Make lead and lag dataframe columns
function make_lead_lag!(bel::DataFrame)
   ## Insert cols for lead and lag of variables
   insertcols!(bel, (Symbol.(string.("discharge", 1:21)) .=> 0)...)
   for i in ["gage_height", "channel_velocity"]
        insertcols!(bel, (Symbol.(string.(i, 1:21)) .=> 0.0)...)
   end
end

## Fill in lead and lag values from continous data
function fill_lead_lag!(bel::DataFrame, cont_dis::DataFrame=cont_dis)
   ## Add lead lag of values from continous discharge to compute current discharge
   for i in ["discharge", "gage_height", "channel_velocity"]
        get_leadlag(bel.datetime, cont_dis, i)
   end
   select!(bel, Not(findall(c -> occursin("21",c), names(bel))) )
end

## Impute Zero or missing values
function impute_z!(bel::DataFrame)
   ## Impute zeros or remove instances with all zeros
   for i in ["discharge", "gage_height", "channel_velocity"]
      impute_var(bel, i)
   end
end

## Take the mean of travel time 1 Hour before
function get_lag(time::DateTime, cont_dis::DataFrame)
        rw = filter(row ->  time-Dates.Hour(2) .<= row.datetime .<= time-Dates.Hour(1) , cont_dis)[!,["datetime","discharge"]]
        round( sum(rw.discharge)/length(rw.discharge), digits=0)
end

max_stars_repo_name: atott/Detroit-Discharge-Prediction

max_issues_repo_name: atott/Detroit-Discharge-Prediction

max_forks_repo_name: atott/Detroit-Discharge-Prediction


text:
mutable struct Mesh
    t
    p
    e
end

function rect_grid2(xmin, xmax, ymin, ymax, nx, ny)
    nt = 2*nx*ny
    np = (nx+1)*(ny+1)

    mesh = Mesh(0, 0, 0)
    mesh.t = zeros(nt, 3)
    mesh.p = zeros(np, 2)

    nxp1 = nx + 1
    nyp1 = ny + 1

    nt  = 0
    for ix = 1:nx
       for iy = 1:ny

          iv  = (ix-1)*nyp1 + iy
          iv1 = iv + nyp1

          nt = nt + 1
          mesh.t[nt, 1] = iv
          mesh.t[nt, 2] = iv1
          mesh.t[nt, 3] = iv1+1

          nt = nt + 1
          mesh.t[nt, 1] = iv
          mesh.t[nt, 2] = iv1+1
          mesh.t[nt, 3] = iv+1
      end
    end

    hx   = (xmax-xmin)/nx
    hy   = (ymax-ymin)/ny
    x    = xmin

    for ix = 1:nx

      i1 = (ix-1)*(ny+1)+1
      i2 = ix*(ny+1)
      mesh.p[i1:i2, 1] = x*ones(nyp1, 1)
      mesh.p[i1:i2, 2] = collect(ymin:hy:ymax)

      x = x + hx
    end

    i1 = nx*(ny+1)+1
    i2 = (nx+1)*(ny+1)
    mesh.p[i1:i2, 1] = xmax*ones(nyp1,1)
    mesh.p[i1:i2, 2] = collect(ymin:hy:ymax)

    mesh.e = ones(2*(nx+ny), 3)

    mesh.e[1:ny, 1] = collect(1:ny)
    mesh.e[1:ny, 2] = collect(2:ny+1)

    mesh.e[ny+1:nx+ny, 1] = collect(ny+1:ny+1:np-1)
    mesh.e[ny+1:nx+ny, 2] = collect(2*(ny+1):ny+1:np)

    mesh.e[nx+ny+1:nx+2*ny, 1] = collect(np-ny:np-1)
    mesh.e[nx+ny+1:nx+2*ny, 2] = collect(np-ny+1:np)

    mesh.e[nx+2*ny+1:2*(nx+ny), 1] = collect(1:ny+1:np-2*ny-1)
    mesh.e[nx+2*ny+1:2*(nx+ny), 2] = collect(ny+2:ny+1:np-ny)
    return mesh
end

function Stiff_Mass(mesh)

    mesh.e[:, 3] .= 1

    N::Int = maximum(maximum(mesh.t[:, 1:3]))

    dirichlet = mesh.e[findall(x->x==1, mesh.e[:, 3]), 1:2]
    dirichlet = unique(sort(dirichlet[:]))
    FreeNodes = setdiff(1:N, dirichlet)

    (Dphi, area, _) = gradbasis(mesh.p, mesh.t)

    NT = size(mesh.t, 1)
    Mt = zeros(NT, 3, 3)
    At = zeros(NT, 3, 3)
    for i = 1:3
        for j = 1:3
            At[:, i, j] = (Dphi[:, 1, i].*Dphi[:, 1, j] + Dphi[:, 2, i].*Dphi[:, 2, j]).*area
            Mt[:, i, j] = area*((i==j)+1)/12
        end
    end

    #%% Assemble the mass matrix in Omega
    M = spzeros(N, N)
    A = spzeros(N, N)
    for i = 1:3
        krow = mesh.t[:, i]
        for j = 1:3
            kcol = mesh.t[:, j]
            M = M + sparse(Int64.(krow), Int64.(kcol), Mt[:, i, j], N, N)
            A = A + sparse(Int64.(krow), Int64.(kcol), At[:, i, j], N, N)
        end
    end
    return A, M, FreeNodes, dirichlet
end
#=
a=rect_grid2(1, 2, 1, 2, 10, 10)
Stiff_Mass(a)
=#

function gradbasis(node, elem)
    NT = size(elem, 1)
    elem = Int64.(elem)

    ve1 = node[elem[:, 3], :]-node[elem[:, 2], :]
    ve2 = node[elem[:, 1], :]-node[elem[:, 3], :]
    ve3 = node[elem[:, 2], :]-node[elem[:, 1], :]
    area = 0.5*(-ve3[:, 1].*ve2[:, 2] + ve3[:, 2].*ve2[:, 1])
    Dlambda = zeros(NT, 2, 3)
    Dlambda[1:NT, :, 3] = [-ve3[:, 2]./(2*area) ve3[:, 1]./(2*area)]
    Dlambda[1:NT, :, 1] = [-ve1[:, 2]./(2*area) ve1[:, 1]./(2*area)]
    Dlambda[1:NT, :, 2] = [-ve2[:, 2]./(2*area) ve2[:, 1]./(2*area)]

    idx = findall(x->x<0, area)
    #idx = (area<0)
    area[idx, :] = -area[idx, :]
    elemSign = ones(NT, 1)
    elemSign[idx] .= -1
    return Dlambda, area, elemSign
end

function f(p)
    f = (sin.(2*pi*p[:, 1]).*sin.(2*pi*p[:, 2]))
    return f
end
function exactu(p, s)
    u = (8*pi^2)^(-s)*sin.(2*pi*p[:, 1]).*sin.(2*pi*p[:, 2])
    return u
end
#s=0.2

function solve()

    for ind = 1:5
        mesh = rect_grid2(0, 1, 0, 1, 2^(ind+1), 2^(ind+1))

        (A, M, FreeNodes, dirichlet) = Stiff_Mass(mesh)

        N = size(A, 1)

        rhs = M*f(mesh.p)

        global u = zeros(N, 1)

        h = N^(-1/2)

        k = 1/log(1/h)

        Nplus  = ceil(pi^2/(4*s*k^2))
        Nminus = ceil(pi^2/(4*(1-s)*k^2))

        for ell=-Nminus:1:Nplus
            y = k*ell;
            u[FreeNodes] = u[FreeNodes] + exp((1-s)*y) * ( (exp(y)*M[FreeNodes, FreeNodes]+A[FreeNodes, FreeNodes]) \ rhs[FreeNodes])
        end

        u[FreeNodes] = (sin(s*pi)/pi) * k * u[FreeNodes]

        uD = zeros(N, 1)
        dirichlet = Int.(dirichlet)
        uD[dirichlet] = exactu(mesh.p[dirichlet, :], s)
        v = uD
        rhs = -A*uD
        v[FreeNodes] = A[FreeNodes, FreeNodes] \ rhs[FreeNodes]

        u = u+v

        #errL2[ind, 1] = sqrt((u-exactu(mesh.p, s))'*M*(u-exactu(mesh.p, s)))

        #hmax[ind, 1] = h
    end
    return u
end
#=
m = rect_grid2(0, 1, 0, 1, 2^(5+1), 2^(5+1))

using Plots
gr()
#pyplot()
#plotlyjs()
#pyplot()
plot3d(m.p[:, 1], m.p[:, 2], u)
=#


text:
include("../include.jl")

#nlp_raw = CUTEstModel("KISSING") # 180 ITS
#nlp_raw = CUTEstModel("KISSING2") # 151 ITS
nlp_raw = CUTEstModel("KISSING2","-param","m=100", "n=8")
#nlp_raw = CUTEstModel("CHANDHEU")
if true
using Ipopt
solver = IpoptSolver(print_level=5, max_iter=3000, bound_relax_factor=0.0, nlp_scaling_method="none")
#,mehrotra_algorithm="yes") #, tol_dual_abs=1e-6)
#solver = IpoptSolver(print_level=5, tol=1e-8)
mp = NLPModels.NLPtoMPB(nlp_raw, solver)
MathProgBase.optimize!(mp)
@show norm(mp.inner.mult_g, Inf)
#y = MathProgBase.getdual(mp)
solver = MathProgBase.getrawsolver(mp)
#finalize(nlp_raw)
end

begin
nlp = Class_CUTEst(nlp_raw)
## FEASIBLE (probably)
timer = class_advanced_timer()
start_advanced_timer(timer)
#include("include.jl")
#intial_it = initial_point_satisfy_bounds(nlp, my_par)
start_advanced_timer(timer, "INIT")
intial_it = init(nlp, my_par, timer)
pause_advanced_timer(timer, "INIT")

#intial_it = initial_point_generic(nlp, my_par, nlp_raw.meta.x0)

@assert(is_feasible(intial_it, my_par.ls.comp_feas))
iter, status, hist, t, err = one_phase_IPM(intial_it, my_par, timer);

pause_advanced_timer(timer)

print_timer_stats(timer)

end
finalize(nlp_raw)

#
# aggressive steps do max LP step

if false
include("include.jl")
x = nlp_raw.meta.x0;
m = nlp_raw.meta.ncon;
@time for i = 1:20 obj(nlp_raw, x) end;
@time for i = 1:20 grad(nlp_raw, x) end;
@time for i = 1:20 cons(nlp_raw, x) end;
@time for i = 1:20 J = jac(nlp_raw, x) end;
@time for i = 1:20 p = jtprod(nlp_raw, x, randn(m))  end;
@time for i = 1:20  jac_coord(nlp_raw, x) end;
nlp = Class_CUTEst(nlp_raw);
@time for i = 1:20 eval_jac(nlp, zeros(10000)) end;
end


text:
module ReadBasis

using BuildBasis,TypesBasis

export readgaussparams, retreadbasis

"""
	use: readgaussparams(basisfile,format="nwchem")
input:
    basisfile - name of file containing basis parameters.
    format (optional) - the format of the basis file.

returns: nbasisfunc - number of basis functions.
         basisfunc - an Array{Any} of basis function parameters

Description: Read the exponent and coefficients for a gaussian basis set. Returns an array of
strings that need to be later processed using TypesBasis.jl and Basis.jl.

Supports: mod NWChem EMSL format, i.e., but ignores species and orbital info by
commenting them out in the file with #.

""" function readgaussparams(filename;format="modnwchem")
            ele = geteledict();
            basisread = String[];
	    if format == "modnwchem"
	    	  file = open(filename,"r");
		  for line in readlines(file)
           	       if occursin("#",line)
               		   continue;
           	       elseif isempty(line)
               		    continue;
           	       else
			    push!(basisread,line);
           	       end
       		  end
		close(file);

                nbasisread = length(basisread);
                i = 1;
                nbasisfunc = 0;
                basisfunc = Array{Array{String,1},1}();
                flag = false;
                while flag != true
                    line = basisread[i]
                    sline = split(line);
                    if haskey(ele,sline[1])
                        symbol = sline[1];
                        orbital = sline[2];
                        nprims = parse(Int,sline[3]);
                        nbasisfunc += 1;
                    else
                        error("Basis file needs start line having element oribtal number-primitives, example, H S 3")
                    end
                    params = Array{String,1}();
                    for ii=1:nprims
                        i += 1
                        push!(params,basisread[i]);
                    end
                    push!(basisfunc,params);
                    i += 1;
                    if i > nbasisread
                        flag = true;
                    end
                end # while flag

	        return nbasisfunc,basisfunc
             else
		error("Format $format is not supported yet.")
	     end
end #readgaussparams

"""
     use: retgaussbasis(myfile, format="nwchem", btype="STO-3G", atomicZ=1)
Input:
Returns: TypesBasis.BasisFunctions

Description: constructs TypeBasis.Gaussian.

""" function retreadbasis(filename;format="modnwchem",btype="6-31G")

    if isequal(btype,"6-31")
        #Parse file
        nbasisfunc,basisfuncstring = readgaussparams(filename,format=format);

        #Construct datatype
        basisfuncarry = Array{Gaussian,1}(undef,nbasisfunc);
        for (ibfs,bfs) in enumerate(basisfuncstring)
            nprims = length(bfs);
            #Unpack Array{String,1} to two Array{Float64,1}
            exponents = Array{Float64,1}(undef,nprims);
            coefficients = Array{Float64,1}(undef,nprims);
            for i=1:nprims
                strexpo,strcoef = split(bfs[i]);
                exponents[i] = parse(Float64,strexpo);
                coefficients[i] = parse(Float64,strcoef);
            end
            basisfuncarry[ibfs] = Gaussian(btype,nprims,coefficients,exponents)
        end #enumerate(basisfuncstring)
        return BasisFunctions(nbasisfunc,basisfuncarry)
    else
        error("The basis set type $btype is not implemented")
    end

end #retreadbasis

geteledict() = Dict("H" => 1, "He" => 2, "Li" => 3);

end #ReadBasis


text:
#=
For this example you need to checkout this package:
Pkg.clone("https://github.com/dpsanders/BilliardModels.jl")
Pkg.checkout("BilliardModels", "time_step")
=#
using GLAbstraction, MeshIO, Colors
using GLVisualize, GeometryTypes, Reactive, ColorTypes
if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(linspace(0f0,1f0,360))
end
const interactive_example = true

using BilliardModels

# create the billiard table
const table = Sinai_billiard(0.1)
const max_particles = 8_000

# function that steps through the simulation
function BilliardModels.step!(particles, table, _)
    for particle in particles
         BilliardModels.step!(particle, table, 0.01)
     end
     particles
end

# convert a particle to a point
function to_points(data, particles)
    @inbounds for (i,p) in enumerate(particles)
        data[i] = to_points(p)
    end
    data
end
to_points(p::BilliardModels.Particle) = Point3f0(p.x.x*2pi, p.x.y*2pi, atan2(p.v.x, p.v.y))

# color lookup table
const colorramp = map(RGBA{Float32}, colormap("RdBu", 100))
function to_color(p::BilliardModels.Particle)
    l = (atan2(p.v.x, p.v.y) + pi) / 2pi
    colorramp[round(Int, clamp(l, 0, 1) * (length(colorramp)-1))+1]
end
function to_color(data, particles)
    @inbounds for (i,p) in enumerate(particles)
        data[i] = to_color(p)
    end
    data
end

cubecamera(window)

x0              = Vector2D(0.3, 0.1)
particles       = [BilliardModels.Particle(x0, Vector2D(1.0, 0.001*i)) for i=1:max_particles]
colors          = RGBA{Float32}[RGBA{Float32}(1., 0.1, clamp(0.001*i, 0.0, 1.0), 1.0) for i=1:max_particles]
particle_stream = const_lift(BilliardModels.step!, particles, table, timesignal)
v0              = map(to_points, particles)
vc0             = map(to_color, particles)
colors          = const_lift(to_color, vc0, particle_stream)
pointstream     = const_lift(to_points, v0, particle_stream)
primitive       = Circle(Point2f0(0), 0.05f0)

# we know that the particles will only be in this range
boundingbox = AABB{Float32}(Vec3f0(-pi), Vec3f0(2pi))
particles = visualize(
    (primitive, pointstream),
    color=colors, # set color array. This is per particle
    billboard=true, # set billboard to true, making the particles always face the camera
    boundingbox=Signal(boundingbox) # set boundingbox, to avoid bb re-calculation when particles update( is expensive)
)

# visualize the boundingbox
boundingbox = visualize(boundingbox, :lines)
# _view them (add them to the windows render list)
_view(particles, window, camera=:perspective)
_view(boundingbox, window, camera=:perspective)

if !isdefined(:runtests)
    renderloop(window)
end


text:
export ObservablesProblem

"""
    ObservablesProblem

The problem representing the estimation of several observables.
"""
struct ObservablesProblem{B,SM} <: OperatorEstimationProblem where {SM}
    HilbSpace::B            # 0
    ObservablesTransposed::Vector{SM}
    Names::Vector{Symbol}
end

"""
    ObservablesProblem([T=STD_REAL_PREC], obs1, obs2, ....)

Constructs an ObservablesProblem from the observables provided.
Userd to compute observables.
"""
ObservablesProblem(args...; kwargs...) = ObservablesProblem(STD_REAL_PREC, args...; kwargs...)

function ObservablesProblem(T::Type{<:Number}, obs::Any...; operator=true)
    if length(obs) == 1
        obs = first(obs)
    end
    if obs isa DataOperator || obs isa AbstractOperator
        obs = [obs]
    end

    names          = Vector{Symbol}()
    matrices_trans = Vector{Any}()

    if operator
        b = nothing
        for (i, el) = enumerate(obs)
            if el isa Tuple
                name = first(el)
                op = last(el)
            else
                name = "obs_$i"
                op = el
            end
            push!(names, Symbol(name))
            push!(matrices_trans, to_linear_operator(op))
            b = basis(op)
        end

        T = typeof(first(matrices_trans))
        if all(T.==typeof.(matrices_trans))
            matrices_trans = Vector{T}(matrices_trans)
        end

        return ObservablesProblem(b, matrices_trans, names)
    else
        b = nothing
        matrices_trans = Vector{SparseMatrixCSC}()
        for (i, el) = enumerate(obs)
            if el isa Tuple
                name = first(el)
                op = last(el)
            else
                name = "obs_$i"
                op = el
            end
            push!(names, Symbol(name))
            if op isa QuantumOptics.DataOperator
                if op isa DenseOperator
                    op = SparseOperator(op)
                end
                push!(matrices_trans, Complex{T}.(op.data))
            elseif isa(op, QuantumLattices.GraphOperator)
                push!(matrices_trans, Complex{T}.(SparseOperator(op).data))
            else
                push!(matrices_trans, Complex{T}.(last(el)))
            end
            b = basis(op)
        end
        T = typeof(first(matrices_trans))
        return ObservablesProblem(b, T.(matrices_trans), names)
    end
end

basis(prob::ObservablesProblem) = prob.HilbSpace

state(T::Type{<:Number}, prob::ObservablesProblem, net::MatrixNet) =
    DiagonalStateWrapper(state(T, basis(prob), net))

#state(T::Type{<:Number}, prob::ObservablesProblem{<:Any, <:AbsLinearOperator}, net) =
#    DiagonalStateWrapper(state_lut(T, basis(prob), net))

Base.show(io::IO, p::ObservablesProblem) = print(io,
    "ObservablesProblem on space : $(basis(p)) for the observables:"*
    prod(["\n\t$obj" for obj=p.Names]))


text:
module TestTypes
using Test
using EcologicalNetworks

@test typeof(UnipartiteNetwork(rand(Bool, (5, 5)))) <: UnipartiteNetwork
@test typeof(UnipartiteNetwork(rand(Bool, (5, 5)))) <: AbstractUnipartiteNetwork
@test typeof(UnipartiteNetwork(rand(Bool, (5, 5)))) <: DeterministicNetwork
@test typeof(UnipartiteNetwork(rand(Bool, (2, 2)), ["a", "b"])) <: DeterministicNetwork

@test typeof(BipartiteNetwork(rand(Bool, (3, 5)))) <: BipartiteNetwork
@test typeof(BipartiteNetwork(rand(Bool, (2, 3)), [:a, :b], [:c, :d, :e])) <: BipartiteNetwork
@test typeof(BipartiteNetwork(rand(Bool, (3, 5)))) <: AbstractBipartiteNetwork
@test typeof(BipartiteNetwork(rand(Bool, (3, 5)))) <: DeterministicNetwork

@test typeof(BipartiteProbabilisticNetwork(rand(7, 3))) <: BipartiteProbabilisticNetwork
@test typeof(BipartiteProbabilisticNetwork(rand(3, 4))) <: AbstractBipartiteNetwork
@test typeof(BipartiteProbabilisticNetwork(rand(2, 3))) <: ProbabilisticNetwork
@test typeof(BipartiteProbabilisticNetwork(rand(2, 3), ["a", "b"], ["c", "d", "e"])) <: ProbabilisticNetwork

@test typeof(UnipartiteProbabilisticNetwork(rand(3, 3))) <: UnipartiteProbabilisticNetwork
@test typeof(UnipartiteProbabilisticNetwork(rand(3, 3))) <: AbstractUnipartiteNetwork
@test typeof(UnipartiteProbabilisticNetwork(rand(3, 3))) <: ProbabilisticNetwork
@test typeof(UnipartiteProbabilisticNetwork(rand(3, 3), ["a", "b", "c"])) <: ProbabilisticNetwork

@test typeof(UnipartiteQuantitativeNetwork(rand(Int64, (5, 5)))) <: UnipartiteQuantitativeNetwork
@test typeof(UnipartiteQuantitativeNetwork(rand(Float64, (5, 5)))) <: QuantitativeNetwork
@test typeof(UnipartiteQuantitativeNetwork(rand(Float64, (2, 2)), [:a, :b])) <: QuantitativeNetwork
@test typeof(UnipartiteQuantitativeNetwork(rand(Float64, (2, 2)), ["a", "b"])) <: QuantitativeNetwork

@test typeof(BipartiteQuantitativeNetwork(rand(Float64, (3, 5)))) <: BipartiteQuantitativeNetwork
@test typeof(BipartiteQuantitativeNetwork(rand(Float64, (5, 6)))) <: AbstractBipartiteNetwork
@test typeof(BipartiteQuantitativeNetwork(rand(Int64, (5, 2)))) <: QuantitativeNetwork
@test typeof(BipartiteQuantitativeNetwork(rand(Float64, (3, 2)), ["a", "b", "c"], ["A", "B"])) <: BipartiteQuantitativeNetwork

@test_throws ArgumentError UnipartiteNetwork(rand(Bool, (5, 3)))
@test_throws ArgumentError UnipartiteNetwork(rand(Bool, (5, 5)), ["a", "b", "c"])
@test_throws ArgumentError UnipartiteNetwork(rand(Bool, (2, 2)), ["a", "b", "c"])
@test_throws ArgumentError BipartiteNetwork(rand(Bool, (3, 2)), ["a", "b", "c"], ["c", "d", "e"])
@test_throws ArgumentError BipartiteNetwork(rand(Bool, (3, 2)), ["a", "b", "c"], ["c", "e"])
@test_throws MethodError BipartiteNetwork(rand(Bool, (3, 2)), ["a", "b", "c"], [:e, :f])

@test typeof(UnipartiteProbabilisticNetwork(rand(3, 3), [:a, :b, :c])) <: UnipartiteProbabilisticNetwork
@test_throws ArgumentError UnipartiteProbabilisticNetwork(rand(3, 3).*2.0, [:a, :b, :c])

end


text:
function check_flux(flux, aquifer::UnconfinedAquifer, index::Int)
    # Check if cell is dry
    if aquifer.head[index] <= aquifer.bottom[index]
        # If cell is dry, no negative flux is allowed
        return max(0, boundary.flux[index])
    else
        return flux
    end
end

# Do nothing for a confined aquifer: aquifer can always provide flux
check_flux(flux, aquifer::ConfinedAquifer, index::Int) = flux

struct River{T} <: AquiferBoundaryCondition
    stage::Vector{T}
    infiltration_conductance::Vector{T}
    exfiltration_conductance::Vector{T}
    bottom::Vector{T}
    flux::Vector{T}
    index::Vector{Int}
end

function flux!(Q, river::River, aquifer)
    for (i, index) in enumerate(river.index)
        ϕ = aquifer.head[index]
        stage = river.stage[i]
        if stage > ϕ
            cond = river.infiltration_conductance[i]
            Δϕ = min(stage - river.bottom[i], stage - ϕ)
        else
            cond = river.exfiltration_conductance[i]
            Δϕ = stage - ϕ
        end
        river.flux[i] = check_flux(cond * Δϕ, aquifer, index)
        Q[index] += river.flux[i]
    end
    return Q
end

struct Drainage{T} <: AquiferBoundaryCondition
    elevation::Vector{T}
    conductance::Vector{T}
    flux::Vector{T}
    index::Vector{Int}
end

function flux!(Q, drainage::Drainage, aquifer)
    for (i, index) in enumerate(drainage.index)
        cond = drainage.conductance[i]
        Δϕ = min(0, drainage.elevation[i] - aquifer.head[index])
        drainage.flux[i] = check_flux(cond * Δϕ, aquifer, index)
        Q[index] += drainage.flux[i]
    end
    return Q
end

struct HeadBoundary{T} <: AquiferBoundaryCondition
    head::Vector{T}
    conductance::Vector{T}
    flux::Vector{T}
    index::Vector{Int}
end

function flux!(Q, headboundary::HeadBoundary, aquifer)
    for (i, index) in enumerate(headboundary.index)
        cond = headboundary.conductance[i]
        Δϕ = headboundary.head[i] - aquifer.head[index]
        headboundary.flux[i] = check_flux(cond * Δϕ, aquifer, index)
        Q[index] += headboundary.flux[i]
    end
    return Q
end

struct Recharge{T} <: AquiferBoundaryCondition
    rate::Vector{T}
    flux::Vector{T}
    index::Vector{Int}
end

function flux!(Q, recharge::Recharge, aquifer)
    for (i, index) in enumerate(recharge.index)
        recharge.flux[i] = check_flux(recharge.rate[i] * aquifer.area[index], aquifer, index)
        Q[index] += recharge.flux[i]
    end
    return Q
end

struct Well{T} <: AquiferBoundaryCondition
    volumetric_rate::Vector{T}
    flux::Vector{T}
    index::Vector{Int}
end

function flux!(Q, well::Well, aquifer)
    for (i, index) in enumerate(well.index)
        well.flux[i] = check_flux(well.volumetric_rate[i], aquifer, index)
        Q[index] += well.flux[i]
    end
    return Q
end


text:
# Two different plot functions depending on the need of having a continuous or discontinuous mesh.
function to_plot_dg_mesh(grid::Grid)
  UnstructuredGrid(grid) |> to_plot_dg_mesh
end

function to_plot_dg_mesh(grid::UnstructuredGrid)
    to_simplex_grid(grid) |> to_dg_mesh
end

function to_plot_mesh(grid::Grid)
  UnstructuredGrid(grid) |> to_plot_mesh
end

function to_plot_mesh(grid::UnstructuredGrid)
    to_simplex_grid(grid) |> to_mesh
end

# For the moment, we simplexify all grids (missing support for quads or general polytopes).
function to_simplex_grid(grid::UnstructuredGrid)
  reffes = get_reffes(grid)
  polys = map(get_polytope,reffes)
  all(is_simplex,polys) ? grid : simplexify(grid)
end

# Create continuous GeometryBasics mesh equivalent of our grid.
function to_mesh(grid::UnstructuredGrid)
  reffes = get_reffes(grid)
  @assert all(reffe->get_order(reffe)==1,reffes)
  @assert all(reffe->is_simplex(get_polytope(reffe)),reffes)
  xs = get_node_coordinates(grid)
  Tp = eltype(eltype(xs))
  Dp = length(eltype(xs))
  ps = collect(reinterpret(GeometryBasics.Point{Dp,Tp},xs))
  cns = get_cell_node_ids(grid)
  Tc = eltype(eltype(cns))
  Dc = num_cell_dims(grid)
  fs = collect(lazy_map(GeometryBasics.NgonFace{Dc+1,Tc},cns))
  GeometryBasics.Mesh(GeometryBasics.connect(ps,fs))
end

# Create discontinuous grid with the corresponding handling of nodal or cell fields.
function to_point(x::VectorValue{D,T}) where {D,T}
  GeometryBasics.Point{D,T}(Tuple(x))
end

function to_dg_points(grid::UnstructuredGrid)
  node_x = get_node_coordinates(grid)
  coords = to_dg_node_values(grid, node_x)
  [to_point(x) for x in coords]
end

function to_dg_mesh(grid::UnstructuredGrid)
  ps = to_dg_points(grid)
  Dc = num_cell_dims(grid)
  Tc = Int32
  cns = Vector{Vector{Tc}}(undef, num_cells(grid))
  i = 1
  for cell in 1:num_cells(grid)
    cn = zeros(Tc, Dc+1)
    for lnode in 1:(Dc+1)
      cn[lnode] = i
      i += one(Tc)
    end
    cns[cell] = cn
  end
  fs = lazy_map(GeometryBasics.NgonFace{Dc+1,Tc}, cns) |> collect
  GeometryBasics.connect(ps, fs) |> GeometryBasics.Mesh
end

function to_dg_node_values(grid::Grid, node_value::Vector)
  cell_nodes = get_cell_node_ids(grid)
  i = 0
  cache = array_cache(cell_nodes)
  for cell in 1:length(cell_nodes)
    nodes = getindex!(cache, cell_nodes, cell)
    for node in nodes
      i += 1
    end
  end
  T = eltype(node_value)
  values = zeros(T,i)
  i = 0
  for cell in 1:length(cell_nodes)
    nodes = getindex!(cache, cell_nodes, cell)
    for node in nodes
      i += 1
      values[i] = node_value[node]
    end
  end
  values
end

function to_dg_cell_values(grid::Grid, cell_value::Vector)
  cell_nodes = get_cell_node_ids(grid)
  i = 0
  cache = array_cache(cell_nodes)
  for cell in 1:length(cell_nodes)
    nodes = getindex!(cache, cell_nodes, cell)
    for node in nodes
      i += 1
    end
  end
  T = eltype(cell_value)
  values = zeros(T,i)
  i = 0
  for cell in 1:length(cell_nodes)
    nodes = getindex!(cache, cell_nodes, cell)
    for node in nodes
      i += 1
      values[i] = cell_value[cell]
    end
  end
  values
end

# Obtain edge and vertex skeletons:
function to_lowdim_grid(grid::Grid, ::Val{D}) where D
  topo = GridTopology(grid)
  labels = FaceLabeling(topo)
  model = DiscreteModel(grid, topo, labels)
  Grid(ReferenceFE{D}, model)
end

to_lowdim_grid(grid::Grid{D}, ::Val{D}) where D = grid

to_face_grid(grid::Grid)   = to_lowdim_grid(grid, Val(2))
to_edge_grid(grid::Grid)   = to_lowdim_grid(grid, Val(1))
to_vertex_grid(grid::Grid) = to_lowdim_grid(grid, Val(0))

function to_face_grid_with_map(grid::Grid)
  function _fixnodeids!(face_to_nodes,face_to_m,face_to_n)
    for face in 1:length(face_to_nodes)
      m = face_to_m[face]
      n = face_to_n[face]
      p = face_to_nodes.ptrs[face]
      if m⋅n < 0
        a = face_to_nodes.data[p]
        face_to_nodes.data[p] = face_to_nodes.data[p+1]
        face_to_nodes.data[p+1] = a
      end
    end
  end
  D = 2
  topo = GridTopology(grid)
  labels = FaceLabeling(topo)
  model = DiscreteModel(grid, topo, labels)
  nfaces = num_faces(model,D)
  Γ = BoundaryTriangulation(model,IdentityVector(nfaces))
  n_Γ = get_normal_vector(Γ)
  x_Γ = CellPoint(Fill(VectorValue(0.0,0.0),nfaces),Γ,ReferenceDomain())
  face_to_n = collect(n_Γ(x_Γ))
  face_to_ϕ = get_cell_map(Γ)
  face_to_Jt = lazy_map(∇,face_to_ϕ)
  Jt = GenericCellField(face_to_Jt,Γ,ReferenceDomain())
  m_Γ = Operation(a->VectorValue(a[1],a[3],a[5])×VectorValue(a[2],a[4],a[6]))(Jt)
  face_to_m = collect(m_Γ(x_Γ))
  g = UnstructuredGrid(Grid(ReferenceFE{D}, model))
  face_to_nodes = get_cell_node_ids(g)
  _fixnodeids!(face_to_nodes,face_to_m,face_to_n)
  face_to_cells = get_faces(topo,D,num_cell_dims(grid))
  face_to_cell = lazy_map(getindex,face_to_cells,Fill(1,nfaces))
  g,face_to_cell
end

# Obtain grid and cellfield from triangulation:
function to_grid(trian::Triangulation)
  vds = visualization_data(trian,"")
  first(vds).grid
end

function to_grid(trian::Triangulation, uh::Any)
  vds = visualization_data(trian, "", cellfields=[""=>uh])
  grid = first(vds).grid
  nodaldata = first(first(vds).nodaldata)[2]
  scalarnodaldata = map(to_scalar, nodaldata)
  grid, scalarnodaldata
end

to_scalar(x) = x
to_scalar(x::VectorValue) = norm(x)


text:
using Pkg.TOML, StaticArrays, LinearAlgebra

@testset "MagneticMoments.jl" begin
	params = TOML.parsefile("MagneticMoments.toml")

	magneticMoment = MPISimulations.initialize(MPISimulations.Langevin,"Langevin",params,Float64)
	for T in [Float64,Float32,]
		magneticMomentT = MPISimulations.Langevin(elementType=T)
		@test magneticMomentT.msat ≈ magneticMoment.msat
		@test magneticMomentT.beta ≈ magneticMoment.beta
	end

	t = randn()
	r = rand(SVector{3,Float64})
	H = rand(SVector{3,Float64})
	H = zero(SVector{3,Float64})
	@test MPISimulations.meanMagneticMoment(magneticMoment,H,r,t) ≈ zero(SVector{3,Float64})
	H = SVector{3,Float64}(1e5,0,0)
	@test MPISimulations.meanMagneticMoment(magneticMoment,H,r,t) ≈ SVector{3,Float64}(magneticMoment.msat,0,0)
	H = SVector{3,Float64}(0,1e5,0)
	@test MPISimulations.meanMagneticMoment(magneticMoment,H,r,t) ≈ SVector{3,Float64}(0,magneticMoment.msat,0)
	H = SVector{3,Float64}(0,0,1e5)
	@test MPISimulations.meanMagneticMoment(magneticMoment,H,r,t) ≈ SVector{3,Float64}(0,0,magneticMoment.msat)
	H = randn(SVector{3,Float64})
	x = magneticMoment.beta*norm(H)
	@test MPISimulations.meanMagneticMoment(magneticMoment,H,r,t) ≈ magneticMoment.msat*(coth(x) - 1/x)*normalize(H)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Methods operating on different special matrix types

# Interconversion between special matrix types

# conversions from Diagonal to other special matrix types
Bidiagonal(A::Diagonal) = Bidiagonal(A.diag, fill!(similar(A.diag, length(A.diag)-1), 0), :U)
SymTridiagonal(A::Diagonal) = SymTridiagonal(A.diag, fill!(similar(A.diag, length(A.diag)-1), 0))
Tridiagonal(A::Diagonal) = Tridiagonal(fill!(similar(A.diag, length(A.diag)-1), 0), A.diag,
                                       fill!(similar(A.diag, length(A.diag)-1), 0))

# conversions from Bidiagonal to other special matrix types
Diagonal(A::Bidiagonal) = Diagonal(A.dv)
SymTridiagonal(A::Bidiagonal) =
    iszero(A.ev) ? SymTridiagonal(A.dv, A.ev) :
        throw(ArgumentError("matrix cannot be represented as SymTridiagonal"))
Tridiagonal(A::Bidiagonal) =
    Tridiagonal(A.uplo == 'U' ? fill!(similar(A.ev), 0) : A.ev, A.dv,
                A.uplo == 'U' ? A.ev : fill!(similar(A.ev), 0))

# conversions from SymTridiagonal to other special matrix types
Diagonal(A::SymTridiagonal) = Diagonal(A.dv)
Bidiagonal(A::SymTridiagonal) =
    iszero(A.ev) ? Bidiagonal(A.dv, A.ev, :U) :
        throw(ArgumentError("matrix cannot be represented as Bidiagonal"))
Tridiagonal(A::SymTridiagonal) =
    Tridiagonal(copy(A.ev), A.dv, A.ev)

# conversions from Tridiagonal to other special matrix types
Diagonal(A::Tridiagonal) = Diagonal(A.d)
Bidiagonal(A::Tridiagonal) =
    iszero(A.dl) ? Bidiagonal(A.d, A.du, :U) :
    iszero(A.du) ? Bidiagonal(A.d, A.dl, :L) :
        throw(ArgumentError("matrix cannot be represented as Bidiagonal"))

# conversions from AbstractTriangular to special matrix types
Bidiagonal(A::AbstractTriangular) =
    isbanded(A, 0, 1) ? Bidiagonal(diag(A, 0), diag(A,  1), :U) : # is upper bidiagonal
    isbanded(A, -1, 0) ? Bidiagonal(diag(A, 0), diag(A, -1), :L) : # is lower bidiagonal
        throw(ArgumentError("matrix cannot be represented as Bidiagonal"))

const ConvertibleSpecialMatrix = Union{Diagonal,Bidiagonal,SymTridiagonal,Tridiagonal,AbstractTriangular}
const PossibleTriangularMatrix = Union{Diagonal, Bidiagonal, AbstractTriangular}

convert(T::Type{<:Diagonal},       m::ConvertibleSpecialMatrix) = m isa T ? m :
    isdiag(m) ? T(m) : throw(ArgumentError("matrix cannot be represented as Diagonal"))
convert(T::Type{<:SymTridiagonal}, m::ConvertibleSpecialMatrix) = m isa T ? m :
    issymmetric(m) && isbanded(m, -1, 1) ? T(m) : throw(ArgumentError("matrix cannot be represented as SymTridiagonal"))
convert(T::Type{<:Tridiagonal},    m::ConvertibleSpecialMatrix) = m isa T ? m :
    isbanded(m, -1, 1) ? T(m) : throw(ArgumentError("matrix cannot be represented as Tridiagonal"))

convert(T::Type{<:LowerTriangular}, m::Union{LowerTriangular,UnitLowerTriangular}) = m isa T ? m : T(m)
convert(T::Type{<:UpperTriangular}, m::Union{UpperTriangular,UnitUpperTriangular}) = m isa T ? m : T(m)

convert(T::Type{<:LowerTriangular}, m::PossibleTriangularMatrix) = m isa T ? m :
    istril(m) ? T(m) : throw(ArgumentError("matrix cannot be represented as LowerTriangular"))
convert(T::Type{<:UpperTriangular}, m::PossibleTriangularMatrix) = m isa T ? m :
    istriu(m) ? T(m) : throw(ArgumentError("matrix cannot be represented as UpperTriangular"))

# Constructs two method definitions taking into account (assumed) commutativity
# e.g. @commutative f(x::S, y::T) where {S,T} = x+y is the same is defining
#     f(x::S, y::T) where {S,T} = x+y
#     f(y::T, x::S) where {S,T} = f(x, y)
macro commutative(myexpr)
    @assert myexpr.head===:(=) || myexpr.head===:function # Make sure it is a function definition
    y = copy(myexpr.args[1].args[2:end])
    reverse!(y)
    reversed_call = Expr(:(=), Expr(:call,myexpr.args[1].args[1],y...), myexpr.args[1])
    esc(Expr(:block, myexpr, reversed_call))
end

for op in (:+, :-)
    for (matrixtype, uplo, converttype) in ((:UpperTriangular, 'U', :UpperTriangular),
                                            (:UnitUpperTriangular, 'U', :UpperTriangular),
                                            (:LowerTriangular, 'L', :LowerTriangular),
                                            (:UnitLowerTriangular, 'L', :LowerTriangular))
        @eval begin
            function ($op)(A::$matrixtype, B::Bidiagonal)
                if B.uplo == $uplo
                    ($op)(A, convert($converttype, B))
                else
                    ($op).(A, B)
                end
            end

            function ($op)(A::Bidiagonal, B::$matrixtype)
                if A.uplo == $uplo
                    ($op)(convert($converttype, A), B)
                else
                    ($op).(A, B)
                end
            end
        end
    end
end

# specialized +/- for structured matrices. If these are removed, it falls
# back to broadcasting which has ~2-10x speed regressions.
# For the other structure matrix pairs, broadcasting works well.

# For structured matrix types with different non-zero diagonals the underlying
# representations must be promoted to the same type.
# For example, in Diagonal + Bidiagonal only the main diagonal is touched so
# the off diagonal could be a different type after the operation resulting in
# an error. See issue #28994

function (+)(A::Bidiagonal, B::Diagonal)
    newdv = A.dv + B.diag
    Bidiagonal(newdv, typeof(newdv)(A.ev), A.uplo)
end

function (-)(A::Bidiagonal, B::Diagonal)
    newdv = A.dv - B.diag
    Bidiagonal(newdv, typeof(newdv)(A.ev), A.uplo)
end

function (+)(A::Diagonal, B::Bidiagonal)
    newdv = A.diag + B.dv
    Bidiagonal(newdv, typeof(newdv)(B.ev), B.uplo)
end

function (-)(A::Diagonal, B::Bidiagonal)
    newdv = A.diag-B.dv
    Bidiagonal(newdv, typeof(newdv)(-B.ev), B.uplo)
end

function (+)(A::Diagonal, B::SymTridiagonal)
    newdv = A.diag+B.dv
    SymTridiagonal(A.diag+B.dv, typeof(newdv)(B.ev))
end

function (-)(A::Diagonal, B::SymTridiagonal)
    newdv = A.diag-B.dv
    SymTridiagonal(newdv, typeof(newdv)(-B.ev))
end

function (+)(A::SymTridiagonal, B::Diagonal)
    newdv = A.dv+B.diag
    SymTridiagonal(newdv, typeof(newdv)(A.ev))
end

function (-)(A::SymTridiagonal, B::Diagonal)
    newdv = A.dv-B.diag
    SymTridiagonal(newdv, typeof(newdv)(A.ev))
end

# this set doesn't have the aforementioned problem

+(A::Tridiagonal, B::SymTridiagonal) = Tridiagonal(A.dl+B.ev, A.d+B.dv, A.du+B.ev)
-(A::Tridiagonal, B::SymTridiagonal) = Tridiagonal(A.dl-B.ev, A.d-B.dv, A.du-B.ev)
+(A::SymTridiagonal, B::Tridiagonal) = Tridiagonal(A.ev+B.dl, A.dv+B.d, A.ev+B.du)
-(A::SymTridiagonal, B::Tridiagonal) = Tridiagonal(A.ev-B.dl, A.dv-B.d, A.ev-B.du)

function (+)(A::Diagonal, B::Tridiagonal)
    newdv = A.diag+B.d
    Tridiagonal(typeof(newdv)(B.dl), newdv, typeof(newdv)(B.du))
end

function (-)(A::Diagonal, B::Tridiagonal)
    newdv = A.diag-B.d
    Tridiagonal(typeof(newdv)(-B.dl), newdv, typeof(newdv)(-B.du))
end

function (+)(A::Tridiagonal, B::Diagonal)
    newdv = A.d+B.diag
    Tridiagonal(typeof(newdv)(A.dl), newdv, typeof(newdv)(A.du))
end

function (-)(A::Tridiagonal, B::Diagonal)
    newdv = A.d-B.diag
    Tridiagonal(typeof(newdv)(A.dl), newdv, typeof(newdv)(A.du))
end

function (+)(A::Bidiagonal, B::Tridiagonal)
    newdv = A.dv+B.d
    Tridiagonal((A.uplo == 'U' ? (typeof(newdv)(B.dl), newdv, A.ev+B.du) : (A.ev+B.dl, newdv, typeof(newdv)(B.du)))...)
end

function (-)(A::Bidiagonal, B::Tridiagonal)
    newdv = A.dv-B.d
    Tridiagonal((A.uplo == 'U' ? (typeof(newdv)(-B.dl), newdv, A.ev-B.du) : (A.ev-B.dl, newdv, typeof(newdv)(-B.du)))...)
end

function (+)(A::Tridiagonal, B::Bidiagonal)
    newdv = A.d+B.dv
    Tridiagonal((B.uplo == 'U' ? (typeof(newdv)(A.dl), newdv, A.du+B.ev) : (A.dl+B.ev, newdv, typeof(newdv)(A.du)))...)
end

function (-)(A::Tridiagonal, B::Bidiagonal)
    newdv = A.d-B.dv
    Tridiagonal((B.uplo == 'U' ? (typeof(newdv)(A.dl), newdv, A.du-B.ev) : (A.dl-B.ev, newdv, typeof(newdv)(A.du)))...)
end

function (+)(A::Bidiagonal, B::SymTridiagonal)
    newdv = A.dv+B.dv
    Tridiagonal((A.uplo == 'U' ? (typeof(newdv)(B.ev), A.dv+B.dv, A.ev+B.ev) : (A.ev+B.ev, A.dv+B.dv, typeof(newdv)(B.ev)))...)
end

function (-)(A::Bidiagonal, B::SymTridiagonal)
    newdv = A.dv-B.dv
    Tridiagonal((A.uplo == 'U' ? (typeof(newdv)(-B.ev), newdv, A.ev-B.ev) : (A.ev-B.ev, newdv, typeof(newdv)(-B.ev)))...)
end

function (+)(A::SymTridiagonal, B::Bidiagonal)
    newdv = A.dv+B.dv
    Tridiagonal((B.uplo == 'U' ? (typeof(newdv)(A.ev), newdv, A.ev+B.ev) : (A.ev+B.ev, newdv, typeof(newdv)(A.ev)))...)
end

function (-)(A::SymTridiagonal, B::Bidiagonal)
    newdv = A.dv-B.dv
    Tridiagonal((B.uplo == 'U' ? (typeof(newdv)(A.ev), newdv, A.ev-B.ev) : (A.ev-B.ev, newdv, typeof(newdv)(A.ev)))...)
end

# fixing uniform scaling problems from #28994
# {<:Number} is required due to the test case from PR #27289 where eltype is a matrix.

function (+)(A::Tridiagonal{<:Number}, B::UniformScaling)
    newd = A.d .+ B.λ
    Tridiagonal(typeof(newd)(A.dl), newd, typeof(newd)(A.du))
end

function (+)(A::SymTridiagonal{<:Number}, B::UniformScaling)
    newdv = A.dv .+ B.λ
    SymTridiagonal(newdv, typeof(newdv)(A.ev))
end

function (+)(A::Bidiagonal{<:Number}, B::UniformScaling)
    newdv = A.dv .+ B.λ
    Bidiagonal(newdv, typeof(newdv)(A.ev), A.uplo)
end

function (+)(A::Diagonal{<:Number}, B::UniformScaling)
    Diagonal(A.diag .+ B.λ)
end

function (+)(A::UniformScaling, B::Tridiagonal{<:Number})
    newd = A.λ .+ B.d
    Tridiagonal(typeof(newd)(B.dl), newd, typeof(newd)(B.du))
end

function (+)(A::UniformScaling, B::SymTridiagonal{<:Number})
    newdv = A.λ .+ B.dv
    SymTridiagonal(newdv, typeof(newdv)(B.ev))
end

function (+)(A::UniformScaling, B::Bidiagonal{<:Number})
    newdv = A.λ .+ B.dv
    Bidiagonal(newdv, typeof(newdv)(B.ev), B.uplo)
end

function (+)(A::UniformScaling, B::Diagonal{<:Number})
    Diagonal(A.λ .+ B.diag)
end

function (-)(A::UniformScaling, B::Tridiagonal{<:Number})
    newd = A.λ .- B.d
    Tridiagonal(typeof(newd)(-B.dl), newd, typeof(newd)(-B.du))
end

function (-)(A::UniformScaling, B::SymTridiagonal{<:Number})
    newdv = A.λ .- B.dv
    SymTridiagonal(newdv, typeof(newdv)(-B.ev))
end

function (-)(A::UniformScaling, B::Bidiagonal{<:Number})
    newdv = A.λ .- B.dv
    Bidiagonal(newdv, typeof(newdv)(-B.ev), B.uplo)
end

function (-)(A::UniformScaling, B::Diagonal{<:Number})
    Diagonal(A.λ .- B.diag)
end

rmul!(A::AbstractTriangular, adjB::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) =
    (B = adjB.parent; rmul!(full!(A), adjoint(B)))
*(A::AbstractTriangular, adjB::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) =
    (B = adjB.parent; *(copyto!(similar(parent(A)), A), adjoint(B)))

# fill[stored]! methods
fillstored!(A::Diagonal, x) = (fill!(A.diag, x); A)
fillstored!(A::Bidiagonal, x) = (fill!(A.dv, x); fill!(A.ev, x); A)
fillstored!(A::Tridiagonal, x) = (fill!(A.dl, x); fill!(A.d, x); fill!(A.du, x); A)
fillstored!(A::SymTridiagonal, x) = (fill!(A.dv, x); fill!(A.ev, x); A)

_small_enough(A::Union{Diagonal, Bidiagonal}) = size(A, 1) <= 1
_small_enough(A::Tridiagonal) = size(A, 1) <= 2
_small_enough(A::SymTridiagonal) = size(A, 1) <= 2

function fill!(A::Union{Diagonal,Bidiagonal,Tridiagonal,SymTridiagonal}, x)
    xT = convert(eltype(A), x)
    (iszero(xT) || _small_enough(A)) && return fillstored!(A, xT)
    throw(ArgumentError("array of type $(typeof(A)) and size $(size(A)) can
    not be filled with $x, since some of its entries are constrained."))
end

one(A::Diagonal{T}) where T = Diagonal(fill!(similar(A.diag, typeof(one(T))), one(T)))
one(A::Bidiagonal{T}) where T = Bidiagonal(fill!(similar(A.dv, typeof(one(T))), one(T)), fill!(similar(A.ev, typeof(one(T))), zero(one(T))), A.uplo)
one(A::Tridiagonal{T}) where T = Tridiagonal(fill!(similar(A.du, typeof(one(T))), zero(one(T))), fill!(similar(A.d, typeof(one(T))), one(T)), fill!(similar(A.dl, typeof(one(T))), zero(one(T))))
one(A::SymTridiagonal{T}) where T = SymTridiagonal(fill!(similar(A.dv, typeof(one(T))), one(T)), fill!(similar(A.ev, typeof(one(T))), zero(one(T))))
# equals and approx equals methods for structured matrices
# SymTridiagonal == Tridiagonal is already defined in tridiag.jl

# SymTridiagonal and Bidiagonal have the same field names
==(A::Diagonal, B::Union{SymTridiagonal, Bidiagonal}) = iszero(B.ev) && A.diag == B.dv
==(B::Bidiagonal, A::Diagonal) = A == B

==(A::Diagonal, B::Tridiagonal) = iszero(B.dl) && iszero(B.du) && A.diag == B.d
==(B::Tridiagonal, A::Diagonal) = A == B

function ==(A::Bidiagonal, B::Tridiagonal)
    if A.uplo == 'U'
        return iszero(B.dl) && A.dv == B.d && A.ev == B.du
    else
        return iszero(B.du) && A.dv == B.d && A.ev == B.dl
    end
end
==(B::Tridiagonal, A::Bidiagonal) = A == B

==(A::Bidiagonal, B::SymTridiagonal) = iszero(B.ev) && iszero(A.ev) && A.dv == B.dv
==(B::SymTridiagonal, A::Bidiagonal) = A == B


text:
#=
utilities for arrays
=#

"""
$(SIGNATURES)

The dimension of the real vectorization of a real or complex vector of length
`len::Int`.
"""
vec_length(::Type{<:Real}, len::Int) = len

vec_length(::Type{<:Complex}, len::Int) = len + len

"""
$(SIGNATURES)

Copy a vector in-place.
"""
vec_copyto!(
    v1::AbstractVecOrMat{T},
    v2::AbstractVecOrMat{T},
    ) where T = copyto!(v1, v2)

"""
$(SIGNATURES)

Copy a complex vector to a real vector in-place.
"""
function vec_copyto!(
    rvec::AbstractVecOrMat{T},
    cvec::AbstractVecOrMat{Complex{T}},
    ) where T
    @assert length(rvec) == 2 * length(cvec)
    k = 1
    @inbounds for i in eachindex(cvec)
        ci = cvec[i]
        rvec[k] = real(ci)
        rvec[k + 1] = imag(ci)
        k += 2
    end
    return rvec
end

"""
$(SIGNATURES)

Copy a real vector to a complex vector in-place.
"""
function vec_copyto!(
    cvec::AbstractVecOrMat{Complex{T}},
    rvec::AbstractVecOrMat{T},
    ) where T
    @assert length(rvec) == 2 * length(cvec)
    k = 1
    @inbounds for i in eachindex(cvec)
        cvec[i] = Complex(rvec[k], rvec[k + 1])
        k += 2
    end
    return cvec
end

# symmetric/Hermitian matrices and svec rescalings

"""
$(SIGNATURES)

The dimension of the vectorized triangle of a real symmetric matrix with side
dimension `side::Int`.
"""
svec_length(side::Int) = div(side * (side + 1), 2)

"""
$(SIGNATURES)

The dimension of the real vectorized triangle of a real symmetric or complex
Hermitian matrix with side dimension `side::Int`.
"""
svec_length(::Type{<:Real}, side::Int) = svec_length(side)

svec_length(::Type{<:Complex}, side::Int) = side^2

"""
$(SIGNATURES)

The side dimension of a real symmetric matrix with vectorized triangle length
`len::Int`.
"""
function svec_side(len::Int)
    side = div(isqrt(1 + 8 * len), 2)
    @assert side * (side + 1) == 2 * len
    return side
end

"""
$(SIGNATURES)

The side dimension of a real symmetric or complex Hermitian matrix with real
vectorized triangle length  `len::Int`.
"""
svec_side(::Type{<:Real}, len::Int) = svec_side(len)

function svec_side(::Type{<:Complex}, len::Int)
    side = isqrt(len)
    @assert side^2 == len
    return side
end

"""
$(SIGNATURES)

The index in the vectorized triangle of a symmetric matrix for element
(`row::Int`, `col::Int`).
"""
function svec_idx(row::Int, col::Int)
    if row < col
        (row, col) = (col, row)
    end
    return div((row - 1) * row, 2) + col
end

"""
$(SIGNATURES)

The indices corresponding to block `block::Int` in a vector of blocks with equal
length `incr::Int`.
"""
block_idxs(incr::Int, block::Int) = (incr * (block - 1) .+ (1:incr))

"""
$(SIGNATURES)

Rescale the elements corresponding to off-diagonals in `arr::AbstractVecOrMat`,
with scaling `scal::Real` and default block increment `incr::Int = 1`.
"""
function scale_svec!(
    arr::AbstractVecOrMat,
    scal::Real;
    incr::Int = 1,
    )
    @assert incr > 0
    n = size(arr, 1)
    (d, r) = divrem(n, incr)
    @assert iszero(r)
    side = svec_side(d)
    k = 1
    for i in 1:side
        for j in 1:(i - 1)
            @inbounds @views @. arr[k:(k + incr - 1), :] *= scal
            k += incr
        end
        k += incr
    end
    @assert k == 1 + n
    return arr
end

"""
$(SIGNATURES)

Copy a real symmetric matrix upper triangle to a svec-scaled vector in-place.
"""
function smat_to_svec!(
    vec::AbstractVector{T},
    mat::AbstractMatrix{T},
    rt2::Real,
    ) where T
    k = 1
    m = size(mat, 1)
    @assert m == size(mat, 2)
    for j in 1:m, i in 1:j
        @inbounds if i == j
            vec[k] = mat[i, j]
        else
            vec[k] = mat[i, j] * rt2
        end
        k += 1
    end
    @assert k == length(vec) + 1
    return vec
end

"""
$(SIGNATURES)

Copy a complex Hermitian matrix upper triangle to a svec-scaled real vector
in-place.
"""
function smat_to_svec!(
    vec::AbstractVector{T},
    mat::AbstractMatrix{Complex{T}},
    rt2::Real,
    ) where T
    k = 1
    m = size(mat, 1)
    @assert m == size(mat, 2)
    for j in 1:m, i in 1:j
        @inbounds if i == j
            vec[k] = real(mat[i, j])
            k += 1
        else
            ck = mat[i, j] * rt2
            vec[k] = real(ck)
            k += 1
            vec[k] = -imag(ck)
            k += 1
        end
    end
    @assert k == length(vec) + 1
    return vec
end

"""
$(SIGNATURES)

Copy a svec-scaled vector to a real symmetric matrix upper triangle in-place.
"""
function svec_to_smat!(
    mat::AbstractMatrix{T},
    vec::AbstractVector{T},
    rt2::Real,
    ) where T
    k = 1
    m = size(mat, 1)
    @assert m == size(mat, 2)
    for j in 1:m, i in 1:j
        @inbounds if i == j
            mat[i, j] = vec[k]
        else
            mat[i, j] = vec[k] / rt2
        end
        k += 1
    end
    @assert k == length(vec) + 1
    return mat
end

"""
$(SIGNATURES)

Copy a svec-scaled real vector to a complex Hermitian matrix upper triangle
in-place.
"""
function svec_to_smat!(
    mat::AbstractMatrix{Complex{T}},
    vec::AbstractVector{T},
    rt2::Real,
    ) where T
    k = 1
    m = size(mat, 1)
    @assert m == size(mat, 2)
    @inbounds for j in 1:m, i in 1:j
        if i == j
            mat[i, j] = vec[k]
            k += 1
        else
            mat[i, j] = Complex(vec[k], -vec[k + 1]) / rt2
            k += 2
        end
    end
    @assert k == length(vec) + 1
    return mat
end

# Kronecker utilities

# compute the real symmetric Kronecker product of a matrix in-place
function symm_kron!(
    skr::AbstractMatrix{T},
    mat::AbstractMatrix{T},
    rt2::T,
    ) where {T <: Real}
    side = size(mat, 1)

    col_idx = 1
    @inbounds for l in 1:side
        for k in 1:(l - 1)
            row_idx = 1
            for j in 1:side
                for i in 1:(j - 1)
                    skr[row_idx, col_idx] =
                        mat[i, k] * mat[j, l] + mat[i, l] * mat[j, k]
                    row_idx += 1
                end
                skr[row_idx, col_idx] = rt2 * mat[j, k] * mat[j, l]
                row_idx += 1
                (row_idx > col_idx) && break
            end
            col_idx += 1
        end

        row_idx = 1
        for j in 1:side
            for i in 1:(j - 1)
                skr[row_idx, col_idx] = rt2 * mat[i, l] * mat[j, l]
                row_idx += 1
            end
            skr[row_idx, col_idx] = abs2(mat[j, l])
            row_idx += 1
            (row_idx > col_idx) && break
        end
        col_idx += 1
    end

    return skr
end

# compute the complex Hermitian Kronecker product of a matrix in-place
function symm_kron!(
    skr::AbstractMatrix{T},
    mat::AbstractMatrix{Complex{T}},
    rt2::T,
    ) where {T <: Real}
    side = size(mat, 1)

    col_idx = 1
    @inbounds for l in 1:side
        for k in 1:(l - 1)
            row_idx = 1
            for j in 1:side
                for i in 1:(j - 1)
                    a = mat[i, k] * mat[l, j]
                    b = mat[j, k] * mat[l, i]
                    spectral_kron_element!(skr, row_idx, col_idx, a, b)
                    row_idx += 2
                end
                c = rt2 * mat[j, k] * mat[l, j]
                skr[row_idx, col_idx] = real(c)
                skr[row_idx, col_idx + 1] = imag(c)
                row_idx += 1
                (row_idx > col_idx) && break
            end
            col_idx += 2
        end

        row_idx = 1
        for j in 1:side
            for i in 1:(j - 1)
                c = rt2 * mat[i, l] * mat[l, j]
                skr[row_idx, col_idx] = real(c)
                skr[row_idx + 1, col_idx] = -imag(c)
                row_idx += 2
            end
            skr[row_idx, col_idx] = abs2(mat[j, l])
            row_idx += 1
            (row_idx > col_idx) && break
        end
        col_idx += 1
    end

    return skr
end

# compute an element of the real spectral Kronecker in-place
function spectral_kron_element!(
    skr::AbstractMatrix{T},
    i::Int,
    j::Int,
    a::T,
    b::T,
    ) where {T <: Real}
    @inbounds skr[i, j] = a + b
    return skr
end

# compute an element of the complex spectral Kronecker in-place
function spectral_kron_element!(
    skr::AbstractMatrix{T},
    i::Int,
    j::Int,
    a::Complex{T},
    b::Complex{T},
    ) where {T <: Real}
    apb = a + b
    amb = a - b
    @inbounds begin
        skr[i, j] = real(apb)
        skr[i + 1, j] = -imag(amb)
        skr[i, j + 1] = imag(apb)
        skr[i + 1, j + 1] = real(amb)
    end
    return skr
end

# compute a real symmetric Kronecker-like outer product of a real or complex
# matrix of eigenvectors and a real symmetric matrix
function eig_dot_kron!(
    skr::AbstractMatrix{T},
    inner::Matrix{T},
    vecs::Matrix{R},
    temp1::Matrix{R},
    temp2::Matrix{R},
    temp3::Matrix{R},
    V::Matrix{R},
    rt2::T,
    ) where {T <: Real, R <: RealOrComplex{T}}
    @assert issymmetric(inner) # must be symmetric (wrapper is less efficient)
    rt2i = inv(rt2)
    d = size(inner, 1)
    copyto!(V, vecs') # allows fast column slices
    V_views = [view(V, :, i) for i in 1:size(inner, 1)]
    scals = (R <: Complex{T} ? (rt2i, rt2i * im) : (rt2i,)) # real and imag parts

    col_idx = 1
    @inbounds for (j, V_j) in enumerate(V_views)
        for i in 1:(j - 1), scal in scals
            mul!(temp3, V_j, V_views[i]', scal, false)
            @. temp2 = inner * (temp3 + temp3')
            mul!(temp1, Hermitian(temp2, :U), V)
            mul!(temp2, V', temp1)
            @views smat_to_svec!(skr[:, col_idx], temp2, rt2)
            col_idx += 1
        end

        mul!(temp2, V_j, V_j')
        temp2 .*= inner
        mul!(temp1, Hermitian(temp2, :U), V)
        mul!(temp2, V', temp1)
        @views smat_to_svec!(skr[:, col_idx], temp2, rt2)
        col_idx += 1
    end

    return skr
end


text:
# this file downloads files which are used during testing the package
# if they are already present and their checksum matches, they are not downloaded again

using BinaryProvider

const testdatadir = @__DIR__

REPO_URL = "https://github.com/yeesian/ArchGDALDatasets/blob/master/"

# remote files with SHA-2 256 hash
remotefiles = [
    ("data/A.tif", "8e8654c56dcb6cfd4b446f811cf3ca84713dbc4768107c1d84a614abeec77898"),
    ("data/point.geojson", "8744593479054a67c784322e0c198bfa880c9388b39a2ddd4c56726944711bd9"),
    ("data/utmsmall.tif", "f40dae6e8b5e18f3648e9f095e22a0d7027014bb463418d32f732c3756d8c54f"),
    ("gdalworkshop/world.tif", "b376dc8af62f9894b5050a6a9273ac0763ae2990b556910d35d4a8f4753278bb"),
    ("ospy/data1/sites.dbf", "7df95edea06c46418287ae3430887f44f9116b29715783f7d1a11b2b931d6e7d"),
    ("ospy/data1/sites.prj", "81fb1a246728609a446b25b0df9ede41c3e7b6a133ce78f10edbd2647fc38ce1"),
    ("ospy/data1/sites.sbn", "198d9d695f3e7a0a0ac0ebfd6afbe044b78db3e685fffd241a32396e8b341ed3"),
    ("ospy/data1/sites.sbx", "49bbe1942b899d52cf1d1b01ea10bd481ec40bdc4c94ff866aece5e81f2261f6"),
    ("ospy/data1/sites.shp", "69af5a6184053f0b71f266dc54c944f1ec02013fb66dbb33412d8b1976d5ea2b"),
    ("ospy/data1/sites.shx", "1f3da459ccb151958743171e41e6a01810b2a007305d55666e01d680da7bbf08"),
    ("ospy/data2/ut_counties.dbf", "94607dfbcb1b547bef9956e5704520003e20d3c1b7d9c38ae6c275a68a359c35"),
    ("ospy/data2/ut_counties.prj", "f10dda97619145fe1e0c6b067675171c5b7f3861755a2daecb18f9a688df4620"),
    ("ospy/data2/ut_counties.shp", "e05104edc9c2120cccad9bd4e316fc3e937f706e8c3ad79118784fcafe78f695"),
    ("ospy/data2/ut_counties.shx", "20d9acbf261bfbcb8a751e6945aeb1589809c9bd5361f5714577d9513bb32bc6"),
    ("ospy/data2/ut_counties.txt", "06585b736091f5bbc62eb040918b1693b2716f550ab306026732e1dfa6cd49a7"),
    ("ospy/data3/cache_towns.dbf", "2344b5195e1a7cbc141f38d6f3214f04c0d43058309b162e877fca755cd1d9fa"),
    ("ospy/data3/cache_towns.sbn", "217e938eb0bec1cdccf26d87e5127d395d68b5d660bc1ecc1d7ec7b3f052f4e3"),
    ("ospy/data3/cache_towns.sbx", "e027b3f67bbb60fc9cf67ab6f430b286fd8a1eaa6c344edaa7da4327485ee9f2"),
    ("ospy/data3/cache_towns.shp", "635998f789d349d80368cb105e7e0d61f95cc6eecd36b34bf005d8c7e966fedb"),
    ("ospy/data3/cache_towns.shx", "0cafc504b829a3da2c0363074f775266f9e1f6aaaf1e066b8a613d5862f313b7"),
    ("ospy/data3/sites.dbf", "7df95edea06c46418287ae3430887f44f9116b29715783f7d1a11b2b931d6e7d"),
    ("ospy/data3/sites.sbn", "198d9d695f3e7a0a0ac0ebfd6afbe044b78db3e685fffd241a32396e8b341ed3"),
    ("ospy/data3/sites.sbx", "49bbe1942b899d52cf1d1b01ea10bd481ec40bdc4c94ff866aece5e81f2261f6"),
    ("ospy/data3/sites.shp", "69af5a6184053f0b71f266dc54c944f1ec02013fb66dbb33412d8b1976d5ea2b"),
    ("ospy/data3/sites.shx", "1f3da459ccb151958743171e41e6a01810b2a007305d55666e01d680da7bbf08"),
    ("ospy/data4/aster.img", "2423205bdf820b1c2a3f03862664d84ea4b5b899c57ed33afd8962664e80a298"),
    ("ospy/data4/aster.rrd", "18e038aabe8fd92b0d12cd4f324bb2e0368343e20cc41e5411a6d038108a25cf"),
    ("ospy/data4/sites.dbf", "7df95edea06c46418287ae3430887f44f9116b29715783f7d1a11b2b931d6e7d"),
    ("ospy/data4/sites.sbn", "198d9d695f3e7a0a0ac0ebfd6afbe044b78db3e685fffd241a32396e8b341ed3"),
    ("ospy/data4/sites.sbx", "49bbe1942b899d52cf1d1b01ea10bd481ec40bdc4c94ff866aece5e81f2261f6"),
    ("ospy/data4/sites.shp", "69af5a6184053f0b71f266dc54c944f1ec02013fb66dbb33412d8b1976d5ea2b"),
    ("ospy/data4/sites.shx", "1f3da459ccb151958743171e41e6a01810b2a007305d55666e01d680da7bbf08"),
    ("ospy/data5/aster.img", "2423205bdf820b1c2a3f03862664d84ea4b5b899c57ed33afd8962664e80a298"),
    ("ospy/data5/aster.rrd", "18e038aabe8fd92b0d12cd4f324bb2e0368343e20cc41e5411a6d038108a25cf"),
    ("ospy/data5/doq1.img", "70b8e641c52367107654962e81977be65402aa3c46736a07cb512ce960203bb7"),
    ("ospy/data5/doq1.rrd", "f9f2fe57d789977090ec0c31e465052161886e79a4c4e10805b5e7ab28c06177"),
    ("ospy/data5/doq2.img", "1e1d744f17e6a3b97dd9b7d8705133c72ff162613bae43ad94417c54e6aced5d"),
    ("ospy/data5/doq2.rrd", "8274dad00b27e008e5ada62afb1025b0e6e2ef2d2ff2642487ecaee64befd914"),
    ("pyrasterio/example.tif", "e4f50bf1f92d8b045c41c6ec101c60e88d663213db9ab90983a66219f28e1b8e"),
    ("pyrasterio/example2.tif", "08145d14e37d68dc586e0ddc74ad5047577ab86c52f2999cd34c0d851a82214c"),
    ("pyrasterio/example3.tif", "e37f6e54d3fb565905a53e9c4a54cc2dd4937df92e35a0aba40470ac166e8688"),
    ("pyrasterio/float_nan.tif", "ae2b8df3ff521d358c1cdd06a3d721ebb028fb19ce89f7e7e83f5efe66d39e1d"),
    ("pyrasterio/float.tif", "758905d08f0f577571bdc083b50c638872c4e1c6465e75076c4ecdd1599e3b24"),
    ("pyrasterio/RGB.byte.tif", "8a9079098409b516d0c3d1965b6ebf64ab69fec12e0c2495d6b200c2ea9dd12a"),
    ("pyrasterio/shade.tif", "a763eb1f223845eeaf6309a11d03acb221b7d54cda526a639fb71a75ba5ec49b")
]

for (f, sha) in remotefiles
    # create the directories if they don't exist
    currdir = dirname(f)
    isdir(currdir) || mkpath(currdir)
    # download the file if it is not there or if it has a different checksum
    currfile = normpath(joinpath(testdatadir, f))
    url = REPO_URL * f * "?raw=true"
    download_verify(url, sha, currfile; force=true)
end


text:
include("quadrature.jl")
# forward recurrence
I1 = log(3.0/2.0)
for n in 1:100
    global I1 = 1.0/n - 2.0*I1
end
println(I)

# backward recurrence
I1 = 1.0/200
for n in 200:-1:101
    global I1 = 1.0/(2.0*n) - I1/2
end
println(I1)

# check result
f(x) = x^100/(x+2)
q = ctr(f)
println(q)


text:
abstract type GradientMode end
abstract type SemiGradient <: GradientMode end
abstract type FullGradient <: GradientMode end

export GradientMode, SemiGradient, FullGradient

include("accumulating.jl")
include("trueonline.jl")
include("tidbd.jl")

export AbstractAccumulatingTraceOptimizer, AccumulatingTraceOptimizer, Parl2AccumulatingTraceOptimizer
export TIDBDOptimizer, AutoTIDBDOptimizer
export AbstractTrueOnlineTraceOptimizer, TrueOnlineTraceOptimizer

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
module Self

export
    indexeswhere

"""
    indexeswhere(A, is)

Find indexes where `A` is `is`.
"""
function indexeswhere(A::AbstractArray{T}, is::T) where T
    IndexStyle(A) == IndexLinear() || throw(ArgumentError("Must be `IndexLinear`"))
    (1:length(A))[(A .== is)[:]]
end

end  # module


text:
struct Gaussian_WG_HNE_DL{K,H} <: UnknownAcvTimeSeriesModel{3,Float64}
    α::Float64
    ωₚ::Float64
    κ::Float64
    ϕₘ::Float64
    σ::Float64

    norm::Float64
    halfκ⁻²::Float64
    cosϕₘ::Float64
    sinϕₘ::Float64
    cos2ϕₘ::Float64
    sin2ϕₘ::Float64
    invexp2σ²::Float64
    invexphalfσ²::Float64
    function Gaussian_WG_HNE_DL{K,H}(α,ωₚ,κ,ϕₘ,σ) where {K,H}
        α > 0 || throw(ArgumentError("α must be > 0."))
        ωₚ > 0 || throw(ArgumentError("ωₚ must be > 0."))
        κ > 0 || throw(ArgumentError("κ must be > 0."))
        σ > 0 || throw(ArgumentError("Gaussian_WG_HNE_DL requires σ > 0"))
        new{K,H}(α,ωₚ,κ,ϕₘ,σ,
        α/(κ*sqrt(2π)),inv(2κ^2),
        cos(ϕₘ),sin(ϕₘ),cos(2ϕₘ),sin(2ϕₘ),exp(-2.0 * σ^2),exp(-0.5 * σ^2)
        )
    end
    function Gaussian_WG_HNE_DL{K,H}(x::AbstractVector{Float64}) where {K,H}
        @boundscheck checkparameterlength(x,Gaussian_WG_HNE_DL{K,H})
        @inbounds Gaussian_WG_HNE_DL{K,H}(x[1], x[2], x[3], x[4], x[5])
    end
end
Gaussian_WG_HNE_DL(x::AbstractVector{Float64}) = Gaussian_WG_HNE_DL(ones(5)...)
Gaussian_WG_HNE_DL(α,ωₚ,κ,ϕₘ,σ) = error("Gaussian_WG_HNE_DL process requires the ammount of aliasing specified as a type parameter. Use Gaussian_WG_HNE_DL{K,H}() where K ∈ N₀.")

WhittleLikelihoodInference.npars(::Type{Gaussian_WG_HNE_DL{K,H}}) where {K,H} = 5
WhittleLikelihoodInference.nalias(::Gaussian_WG_HNE_DL{K,H}) where {K,H} = K

WhittleLikelihoodInference.lowerbounds(::Type{Gaussian_WG_HNE_DL{K,H}}) where {K,H} = [0,0,0,-Inf,0]
WhittleLikelihoodInference.upperbounds(::Type{Gaussian_WG_HNE_DL{K,H}}) where {K,H} = [Inf,Inf,Inf,Inf,Inf]

@propagate_inbounds @fastmath function WhittleLikelihoodInference.add_sdf!(out, model::Gaussian_WG_HNE_DL{K,H}, ω) where {K,H}
    @boundscheck checkbounds(out,1:6)
    @inbounds begin
        signω = sign(ω)
        ω = abs(ω)
        if ω > 1e-10 # note this is due to the depth limiting, and results in a spectra that for some impractical parameters will be meaningfully discontinuous near 0.
            sdf = model.norm*exp(-(model.ωₚ-abs(ω))^2*model.halfκ⁻²)

            cos2part = model.cos2ϕₘ * model.invexp2σ²
            imexppart = sdf * 1.0im * signω * model.invexphalfσ²

            tanhkh = tanh(approx_dispersion(ω, H))
            inv_tanhkh = inv(tanhkh)
            inv_tanhkh² = inv_tanhkh^2

            out[1] += sdf
            out[2] += imexppart * model.cosϕₘ * inv_tanhkh
            out[3] += imexppart * model.sinϕₘ * inv_tanhkh
            out[4] += sdf * (0.5 + 0.5cos2part) * inv_tanhkh²
            out[5] += sdf * model.sin2ϕₘ *0.5*model.invexp2σ² * inv_tanhkh²
            out[6] += sdf * (0.5 - 0.5cos2part) * inv_tanhkh²
        end
    end
    return nothing
end

@propagate_inbounds @fastmath function WhittleLikelihoodInference.grad_add_sdf!(out, model::Gaussian_WG_HNE_DL{K,H}, ω::Real) where {K,H}
    @boundscheck checkbounds(out,1:6,1:5)
    @inbounds begin
        signω = sign(ω)
        ω = abs(ω)
        if ω > 1e-10 # note this is due to the depth limiting, and results in a spectra that for some impractical parameters will be meaningfully discontinuous near 0.
            σ = model.σ

            sdf = model.norm*exp(-(model.ωₚ-ω)^2 * model.halfκ⁻²)
            ∂S∂α = sdf/model.α
            ∂S∂ωₚ = 2(ω-model.ωₚ) * model.halfκ⁻² * sdf
            ∂S∂κ = sdf * ((ω-model.ωₚ)^2/model.κ^3-1/model.κ)

            tanhkh = tanh(approx_dispersion(ω, H))
            inv_tanhkh = inv(tanhkh)
            inv_tanhkh² = inv_tanhkh^2

            cos2part = model.cos2ϕₘ * model.invexp2σ²
            imexppart = 1.0im * signω * model.invexphalfσ²
            dirxx = 0.5 + 0.5cos2part
            diryy = 0.5 - 0.5cos2part
            dirxz = imexppart * model.cosϕₘ
            diryz = imexppart * model.sinϕₘ
            diryx = model.sin2ϕₘ *0.5 * model.invexp2σ²
            ## α
            out[1, 1] +=∂S∂α
            out[2, 1] +=∂S∂α * dirxz * inv_tanhkh
            out[3, 1] +=∂S∂α * diryz * inv_tanhkh
            out[4, 1] +=∂S∂α * dirxx * inv_tanhkh²
            out[5, 1] +=∂S∂α * diryx * inv_tanhkh²
            out[6, 1] +=∂S∂α * diryy * inv_tanhkh²
            ## ωₚ
            out[1, 2] +=∂S∂ωₚ
            out[2, 2] +=∂S∂ωₚ * dirxz * inv_tanhkh
            out[3, 2] +=∂S∂ωₚ * diryz * inv_tanhkh
            out[4, 2] +=∂S∂ωₚ * dirxx * inv_tanhkh²
            out[5, 2] +=∂S∂ωₚ * diryx * inv_tanhkh²
            out[6, 2] +=∂S∂ωₚ * diryy * inv_tanhkh²
            ## γ
            out[1, 3] +=∂S∂κ
            out[2, 3] +=∂S∂κ * dirxz * inv_tanhkh
            out[3, 3] +=∂S∂κ * diryz * inv_tanhkh
            out[4, 3] +=∂S∂κ * dirxx * inv_tanhkh²
            out[5, 3] +=∂S∂κ * diryx * inv_tanhkh²
            out[6, 3] +=∂S∂κ * diryy * inv_tanhkh²
            ## ϕₘ
            horTemp = sdf * model.sin2ϕₘ * model.invexp2σ²
            # out[1, 4] += 0
            out[2, 4] +=-1im * sdf * model.sinϕₘ * model.invexphalfσ² * signω * inv_tanhkh
            out[3, 4] +=1im * sdf * model.cosϕₘ * model.invexphalfσ² * signω * inv_tanhkh
            out[4, 4] +=-horTemp * inv_tanhkh²
            out[5, 4] +=sdf * model.cos2ϕₘ * model.invexp2σ² * inv_tanhkh²
            out[6, 4] +=horTemp  * inv_tanhkh²
            ## σ
            sigtemp = 2.0σ * sdf * cos2part * inv_tanhkh²
            # out[1, 5] += 0
            out[2, 5] +=-σ * sdf * dirxz * inv_tanhkh
            out[3, 5] +=-σ * sdf * diryz * inv_tanhkh
            out[4, 5] +=-sigtemp
            out[5, 5] +=-4.0σ * sdf * diryx  * inv_tanhkh²
            out[6, 5] += sigtemp
        end
    end
    return nothing
end


text:
baremodule GitDataManager

module Types
struct CLI{API <: Function} <: Function
    api::API
end
end

module API
function commit end
function upload end
function empty end
end

"""
    GitDataManager.commit(dir; <keyword arguments>)

Commit files in the Git repository `dir`.

# Keyword Arguments
- `message = "Update @HOSTNAME"`: Commit message.
"""
const commit = Types.CLI(API.commit)

"""
    GitDataManager.upload(dir; <keyword arguments>)

Commit files in the Git repository `dir` and push it to remote.

# Keyword Arguments
- `merge = true`: Merge remote head if required.
- `message = "Update @HOSTNAME"`: Commit message.
"""
const upload = Types.CLI(API.upload)

"""
    GitDataManager.empty(dir)

Remove all tracked files in `dir`.
"""
const empty = Types.CLI(API.empty)

module Internal

using ..GitDataManager: API, GitDataManager, Types
using Sockets: gethostname

include("internal.jl")

end  # module Internal

end  # baremodule GitDataManager


text:
mutable struct Model <: TwoDotsModel
    real_time::Bool
end

Model(;real_time=false) = Model(real_time)

function decide!(model::Model, game, gui)
    options = search_all(game)
    n_dots = length.(options)
    _,idx = findmax(n_dots)
    choice = options[idx]
    select_choice!(model, game, gui, choice)
    model.real_time ? sleep(1) : nothing
    click_submit!(game, gui)
    return nothing
end

function select_choice!(model, game, gui, choice)
    for dot in choice
        click_dot!(game, dot, gui)
        model.real_time ? sleep(.5) : nothing
    end
end


text:
using Statistics
using DataAPI

abstract type FitResult end
"""
A measured set of correlated k-ratios.
"""
struct BasicFitResult <: FitResult
    label::Label
    kratios::UncertainValues
end

function findlabel(ffr::FitResult, cxr::CharXRay)
    lbls = labels(ffr.kratios)
    fa = findall(lbl -> (lbl isa CharXRayLabel) && (cxr in lbl.xrays), lbls)
    @assert length(fa) ≠ 0 "No label found for $cxr."
    @assert length(fa) == 1 "Multiple labels found for $cxr."
    return lbls[fa[1]]
end

Base.show(io::IO, ffr::FitResult) = print(io, "$(ffr.label)")
NeXLUncertainties.NeXLUncertainties.value(ffr::FitResult, label::ReferenceLabel) =
    NeXLUncertainties.value(ffr.kratios, label)
NeXLUncertainties.σ(ffr::FitResult, label::ReferenceLabel) = σ(ffr.kratios, label)
NeXLUncertainties.uncertainty(ffr::FitResult, label::ReferenceLabel) =
    uncertainty(ffr.kratios, label)
NeXLUncertainties.getindex(ffr::FitResult, label::ReferenceLabel) =
    getindex(ffr.kratios, label)
Base.keys(ffr::FitResult) = keys(ffr.kratios)
NeXLUncertainties.labels(ffr::FitResult) = labels(ffr.kratios)
function kratio(cxr::CharXRay, ffr::FitResult)
    lbl = findlabel(ffr, cxr)
    return uv(NeXLUncertainties.value(ffr, lbl), σ(ffr, lbl))
end
unknown(ffr::FitResult)::Label = ffr.label
Base.values(ffrs::Vector{<:FitResult}, lbl::ReferenceLabel) =
    [NeXLUncertainties.value(ffr.kratios, lbl, 0.0) for ffr in ffrs]
σs(ffrs::Vector{<:FitResult}, lbl::ReferenceLabel) =
    [σ(ffr.kratios, lbl, 0.0) for ffr in ffrs]

"""
    kratios(ffr::FitResult)::Vector{KRatio}

The k-ratios associated with each `CharXRayLabel` as a vector 'KRatio' objects.
"""
function kratios(ffr::FitResult)::Vector{KRatio}
    lbls = filter(labels(ffr.kratios)) do l
        (l isa CharXRayLabel) && haskey(properties(l)[:Composition], element(l))
    end
    return map(lbls) do lbl
        props = properties(lbl)
        KRatio(lbl.xrays, properties(unknown(ffr)), props, props[:Composition], ffr.kratios[lbl])
    end
end

"""
    kratio(ffr::FitResult, cxr::CharXRay)::KRatio

Extract the k-ratio associated with the specified CharXRay (zero if not measured)
"""
function kratio(ffr::FitResult, cxr::CharXRay)
    lbls = labels(ffr.kratios)
    i = findfirst(lbl->(lbl isa CharXRayLabel) && (cxr in lbl.xrays), lbls)
    @assert !isnothing(i) "There is no k-ratio associated with $cxr."
    lbl, props = lbls[i], properties(lbls[i])
    KRatio(lbl.xrays, unknown(ffr).spectrum.properties, props, props[:Composition], ffr.kratios[lbl])
  end

"""
    extractStandards(ffr::FitResult, elm::Element, mat::Material)::Vector{KRatio}

Extract a `Vector{KRatio}` for `elm::Element` from a `ffr::FilterFitResult` measured from `mat::Material`.
"""
function extractStandards(ffr::FitResult, elm::Element, mat::Material)::Vector{KRatio}
    @assert haskey(mat, elm) "The element $elm is not contained within $mat."
    lbls = filter(labels(ffr.kratios)) do l
        (l isa CharXRayLabel) && (elm == element(l))
    end
    krs = map(lbls) do lbl
        props = properties(lbl)
        kr = KRatio(lbl.xrays, properties(unknown(ffr)), props, props[:Composition], ffr.kratios[lbl])
        kr.unkProps[:Composition] = mat
        kr
    end
    length(krs) == 0 && @info  "There are no k-ratios associated with the $elm in $ffr."
    return krs
end

"""
    extractStandards(ffr::FitResult, cxrs::AbstractVector{CharXRay}, mat::Material)::Vector{KRatio}

Extract a `KRatio` for the `CharXRay` from a `FilterFitResult` associated with
the `Material`.
"""
function extractStandards(ffr::FitResult, cxrs::AbstractVector{CharXRay}, mat::Material)::Vector{KRatio}
    lbls = filter(labels(ffr.kratios)) do lbl
        (lbl isa CharXRayLabel) && any(cxr->cxr in xrays(lbl), cxrs)
    end
    krs = map(lbls) do lbl
        props = properties(lbl)
        kr = KRatio(lbl.xrays, properties(unknown(ffr)), props, props[:Composition], ffr.kratios[lbl])
        kr.unkProps[:Composition] = mat
        kr
    end
    length(krs) == 0 && @info "There are no k-ratios associated with any of $cxrs in $ffr."
    return krs
end

"""
    heterogeneity(ffrs::Vector{FilterFitResult}, lbl::ReferenceLabel)

Computes the ratio of the standard deviation of the measured values over the mean calculated uncertainty
from the fit.  A value near 1 means the sample appears homogeneous and a value greater than 1 means the sample
appears heterogeneous.
"""
heterogeneity(ffrs::Vector{<:FitResult}, lbl::ReferenceLabel) =
    std(values(ffrs, lbl)) / mean(σs(ffrs, lbl))

"""
    NeXLUncertainties.asa(::Type{DataFrame}, ffrs::AbstractVector{<:FitResult}; charOnly = true, withUnc = false, pivot = false)

Return generic `FitResult` as a DataFrame.
"""
function NeXLUncertainties.asa(
    ::Type{DataFrame},
    ffrs::AbstractVector{<:FitResult};
    charOnly = true,
    withUnc = false,
    pivot = false,
)::DataFrame
    lbls = sort(collect(Set(mapreduce(labels, append!, ffrs))))
    lbls = charOnly ? filter(lbl -> lbl isa CharXRayLabel, lbls) : lbls
    if pivot
        res = DataFrame(ROI = lbls)
        for ffr in ffrs
            vals = [value(ffr.kratios, lbl, missing) for lbl in lbls]
            insertcols!(res, ncol(res) + 1, Symbol(repr(ffr.label)) => vals)
            if withUnc
                unc = [σ(ffr.kratios, lbl, missing) for lbl in lbls]
                insertcols!(res, ncol(res) + 1, Symbol('Δ' * repr(ffr.label)) => unc)
            end
        end
    else
        rowLbls = [repr(ffr.label) for ffr in ffrs]
        res = DataFrame(Symbol("Spectra") => rowLbls)
        for lbl in lbls
            vals = [value(ffr.kratios, lbl, missing) for ffr in ffrs]
            insertcols!(res, ncol(res) + 1, Symbol(repr(lbl)) => vals)
            if withUnc
                unc = [σ(ffr.kratios, lbl, missing) for ffr in ffrs]
                insertcols!(res, ncol(res) + 1, Symbol('Δ' * repr(lbl)) => unc)
            end
        end
    end
    return res
end

function NeXLUncertainties.asa(::Type{DataFrame}, fr::FitResult; withUnc = false)
    rois = labels(fr)
    res =  DataFrame(
        ROI = rois,
        k = map(l -> value(fr.kratios, l), rois),
    )
    if withUnc
        res[:, "σ(k)"] = map(l -> σ(fr.kratios, l), rois)
    end
    return res
end

function DataAPI.describe(ffrs::Vector{<:FitResult})::DataFrame
    df = asa(DataFrame, ffrs)[:, 2:end]
    desc = describe(df, :mean, :std, :min, :q25, :median, :q75, :max)
    lbls = filter(
        lbl -> lbl isa CharXRayLabel,
        sort(collect(Set(mapreduce(labels, append!, ffrs)))),
    )
    insertcols!(desc, 4, :hetero => heterogeneity.(Ref(ffrs), lbls))
    return desc
end

"""
    FilterFitResult

Represents the result of fitting a FilteredUnknownW to a FilteredUnknown.

Struct elements

    label::UnknownLabel  # Identifies the unknown
    kratios::UncertainValues # Labeled with ReferenceLabel objects
    roi::UnitRange{Int} # Range of channels fit
    raw::Vector{Float64} # Raw spectrum data
    residual::Vector{Float64} # Residual spectrum
"""
struct FilterFitResult <: FitResult
    label::UnknownLabel
    kratios::UncertainValues
    roi::UnitRange{Int}
    raw::Vector{Float64}
    residual::Vector{Float64}
    peakback::Dict{ReferenceLabel,NTuple{2,Float64}}
end

"""
    residual(ffr::FilterFitResult)::Spectrum

A Spectrum containing the histogram representing the unknown spectrum
minus the fitted characteristic peaks shapes times the best fit coefficient.
"""
function residual(ffr::FilterFitResult)::Spectrum
    props = copy(properties(ffr.label.spectrum))
    props[:Name] = "Residual[$(props[:Name])]"
    return Spectrum(ffr.label.spectrum.energy, ffr.residual, props)
end

"""
    spectrum(ffr::FilterFitResult)::Spectrum

Returns the original unknown spectrum.
"""
spectrum(ffr::FilterFitResult)::Spectrum = ffr.label.spectrum

"""
    characteristiccounts(ffr::FiterFitResult, strip)

Number of spectrum counts that were accounted for by the fitted elements with the `strip` Element(s) removed.
"""
characteristiccounts(ffr::FilterFitResult, strip) =
    sum(element(ref) in strip ? 0.0 : v[1] - v[2] for (ref, v) in ffr.peakback) # sum(ffr.raw[ffr.roi]-ffr.residual[ffr.roi])

"""
    peaktobackground(ffr::FilterFitResult, backwidth::Float64=10.0)::Float64

The peak-to-background ratio as determined from the raw and residual spectra integrated over
the fit region-of-interest and scaled to `backwidth` eV of continuum (nominally 10 eV).
"""
function peaktobackground(
    ffr::FilterFitResult,
    klabel::ReferenceLabel,
    backwidth::Float64 = 10.0,
)::Float64
    unk = spectrum(unknown(ffr))
    peak, back = ffr.peakback[klabel]
    return (peak * (energy(klabel.roi.stop, unk) - energy(klabel.roi.start, unk))) /
           (backwidth * back)
end

"""
    NeXLUncertainties.asa(::Type{DataFrame}, ffr::FilterFitResult; charOnly::Bool=true, material=nothing)::DataFrame

Tabulate details about each region-of-interest in the 'FilterFitResult' in a 'DataFrame'.
  * If charOnly then only display characteristic X-ray data (not escapes etc.)
  * If `material` is a Material then the computed k-ratio will also be tabulated along with kmeas/kcalc.
"""
function NeXLUncertainties.asa(
    ::Type{DataFrame},
    ffr::FilterFitResult;
    charOnly::Bool = true,
    material = nothing,
)::DataFrame
    lbls, klbl, std, kr, dkr, roi1, roi2, peak, back, ptob = UnknownLabel[],
    ReferenceLabel[],
    String[],
    Float64[],
    Float64[],
    Int[],
    Int[],
    Float64[],
    Float64[],
    Float64[]
    kcalc, ratio = Float64[], Float64[]
    for lbl in NeXLUncertainties.sortedlabels(ffr.kratios)
        if (!charOnly) || (lbl isa CharXRayLabel)
            stdspec, unkspec = properties(lbl), ffr.label.spectrum
            push!(lbls, ffr.label)
            push!(std, stdspec[:Name])
            push!(klbl, lbl)
            push!(roi1, lbl.roi.start)
            push!(roi2, lbl.roi.stop)
            push!(kr, NeXLUncertainties.value(ffr.kratios, lbl))
            push!(dkr, σ(ffr.kratios, lbl))
            pb = ffr.peakback[lbl]
            push!(peak, pb[1])
            push!(back, pb[2])
            push!(ptob, peaktobackground(ffr, lbl))
            if material isa Material
                kc = NaN64
                if lbl isa CharXRayLabel
                    try
                        zcu = zafcorrection(
                            XPP,
                            ReedFluorescence,
                            NullCoating,
                            material,
                            lbl.xrays,
                            unkspec[:BeamEnergy],
                        )
                        zcs = zafcorrection(
                            XPP,
                            ReedFluorescence,
                            NullCoating,
                            stdspec[:Composition],
                            lbl.xrays,
                            stdspec[:BeamEnergy],
                        )
                        kc = k(zcu, zcs, unkspec[:TakeOffAngle], stdspec[:TakeOffAngle])
                    catch
                        kc = NaN64
                    end
                end
                push!(kcalc, kc)
                push!(ratio, NeXLUncertainties.value(ffr.kratios, lbl) / kc)
            end
        end
    end
    return material isa Material ?
           DataFrame(
        Spectrum = lbls,
        Feature = klbl,
        Reference = std,
        Start = roi1,
        Stop = roi2,
        K = kr,
        dK = dkr,
        Peak = peak,
        Back = back,
        PtoB = ptob,
        Kcalc = kcalc,
        Ratio = ratio,
    ) :
           DataFrame(
        Spectrum = lbls,
        Feature = klbl,
        Reference = std,
        Start = roi1,
        Stop = roi2,
        K = kr,
        dK = dkr,
        Peak = peak,
        Back = back,
        PtoB = ptob,
    )
end

"""
    filteredresidual(fit::FilterFitResult, unk::FilteredUnknown, ffs::AbstractVector{FilteredReference})::Vector{Float64}

Computes the difference between the best fit and the unknown filtered spectral data.
"""
function filteredresidual(
    fit::FilterFitResult,
    unk::FilteredUnknown,
    ffs::AbstractVector{FilteredReference},
)::Vector{Float64}
    return unk.filtered - mapreduce(+, ffs) do ff
        (value(fit, ff.label) * (ff.scale / unk.scale)) * extract(ff, unk.ffroi)
    end
end

"""
    NeXLMatrixCorrection.estimatecoating(fr::FitResult, substrate::Material, coating::Material, cxr::CharXRay, mc::Type{<:MatrixCorrection}=XPP)::Film

Estimate the mass-thickness of `coating` on bulk `substrate` using the X-ray `cxr` for which there is KRatio in `fr`.  The result is
assigned to the properties of `fr` so that it can be account for in the matrix correction process.

This method is intended for use on standards where the substrate composition is known a priori.
"""
function NeXLMatrixCorrection.estimatecoating(fr::FitResult, substrate::Material, coating::Material, cxr::CharXRay, mc::Type{<:MatrixCorrection}=XPP)::Film
    coating[:Density] = get(coating, :Density, 1.0)
    properties(fr.label.spectrum)[:Coating] = estimatecoating(substrate, coating, kratio(fr, cxr), mc)
end