text:
# Define some new primitives: conv4 and pool

"""

`conv4(w,x;kwargs...)` executes convolutions or cross-correlations
using filters specified with `w` over tensor `x`.  Currently 4 or 5
dimensional KnetArrays with Float32 or Float64 entries are supported.

If `w` has dimensions (W1,W2,...,I,O) and `x` has dimensions
(X1,X2,...,I,N), the result y will have dimensions (Y1,Y2,...,O,N)
where

    Yi=1+floor(Xi-Wi+2*padding[i]/stride[i])

Here I is the number of input channels, O is the number of output
channels, N is the number of instances, and Wi,Xi,Yi are spatial
dimensions.  Padding and stride are keyword arguments that can be
specified as a single number (in which case they apply to all
dimensions), or an array/tuple with entries for each spatial
dimension.

Here is a description of all available keyword arguments:

* padding: the number of extra zeros implicitly concatenated at the start and at the end of each dimension. Default=0.
* stride: the number of elements to slide to reach the next filtering window. Default=1.
* upscale: upscale factor for each dimension. Default=1.
* mode: 0 for convolution and 1 for cross-correlation.  Default=0.
* alpha: can be used to scale the result. Default=1.
* algo: specifies which convolution algorithm shoud be used to compute the results. Default=0. See the CUDNN User Guide for details.
* workSpace: data pointer to GPU memory to a workspace needed to able to execute the specified algorithm. Default=C_NULL.
* workSpaceSizeInBytes: the size in bytes of the provided workSpace. Default=0.
* handle: handle to a previously created cuDNN context. Default=Knet allocated context.

"""
function conv4{T}(w::KnetArray{T},x::KnetArray{T};
                  handle=cudnnhandle, alpha=one(T), beta=zero(T),
                  algo=0, workSpace=C_NULL, workSpaceSizeInBytes=0, o...)
    y = similar(x, cdims(w,x;o...))
    @cuda(cudnn, cudnnConvolutionForward,
          (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,   UInt32,Cptr,     Csize_t,             Ptr{T},Cptr,Ptr{T}),
          handle,Ref(alpha),TD(x),x,FD(w),w,CD(w,x;o...),algo,workSpace,workSpaceSizeInBytes,Ref(beta),TD(y),y)
    return y
end

function conv4x{T}(w::KnetArray{T},x::KnetArray{T},dy::KnetArray{T};
                   handle=cudnnhandle, alpha=one(T), beta=zero(T),
                   algo=0, workSpace=C_NULL, workSpaceSizeInBytes=0, o...)
    dx = similar(x)
    @cuda(cudnn,cudnnConvolutionBackwardData,
          (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,     UInt32,Cptr,     Csize_t,             Ptr{T},Cptr,Ptr{T}),
          handle,Ref(alpha),FD(w),w,TD(dy),dy,CD(w,x;o...),algo,workSpace,workSpaceSizeInBytes,Ref(beta),TD(dx),dx)
    return dx
end

function conv4w{T}(w::KnetArray{T},x::KnetArray{T},dy::KnetArray{T};
                   handle=cudnnhandle, alpha=one(T), beta=zero(T),
                   algo=0, workSpace=C_NULL, workSpaceSizeInBytes=0, o...)
    dw = similar(w)
    @cuda(cudnn,cudnnConvolutionBackwardFilter,
          (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,     UInt32,Cptr,     Csize_t,             Ptr{T},Cptr,Ptr{T}),
          handle,Ref(alpha),TD(x),x,TD(dy),dy,CD(w,x;o...),algo,workSpace,workSpaceSizeInBytes,Ref(beta),FD(dw),dw)
    return dw
end

@primitive conv4(w,x; o...),dy  conv4w(w,x,dy;o...)  conv4x(w,x,dy;o...)
@zerograd  conv4x(w,x,dy;o...)
@zerograd  conv4w(w,x,dy;o...)

"""

`pool(x;kwargs...)` computes pooling of input values (i.e., the
maximum or average of several adjacent values) to produce an output
with smaller height and/or width.  Currently 4 or 5 dimensional
KnetArrays with Float32 or Float64 entries are supported.

If `x` has dimensions (X1,X2,...,I,N), the result y will have
dimensions (Y1,Y2,...,I,N) where

   Yi=1+ceil(Xi-window[i]+2*padding[i])/stride[i]

Here I is the number of input channels, N is the number of instances,
and Xi,Yi are spatial dimensions.  Window, padding and stride are
keyword arguments that can be specified as a single number (in which
case they apply to all dimensions), or an array/tuple with entries for
each spatial dimension.

Here is a description of all available keyword arguments:

* window: the pooling window size for each dimension. Default=2.
* padding: the number of extra zeros implicitly concatenated at the start and at the end of each dimension. Default=0.
* stride: the number of elements to slide to reach the next pooling window. Default=same as window.
* mode: 0 for max, 1 for average including padded values, 2 for average excluding padded values.  Default=0.
* maxpoolingNanOpt: Nan numbers are not propagated if 0, they are propagated if 1. Default=0.
* alpha: can be used to scale the result. Default=1.
* handle: Handle to a previously created cuDNN context. Default=Knet allocated context.

"""
function pool{T}(x::KnetArray{T}; handle=cudnnhandle, alpha=one(T), beta=zero(T), o...)
    y = similar(x, pdims(x; o...))
    @cuda(cudnn, cudnnPoolingForward,
          (Cptr, Cptr,      Ptr{T},    Cptr,Ptr{T},Ptr{T},   Cptr,Ptr{T}),
          handle,PD(x;o...),Ref(alpha),TD(x),x,    Ref(beta),TD(y),y)
    return y
end

function poolx{T}(x::KnetArray{T},y::KnetArray{T},dy::KnetArray{T};
                  handle=cudnnhandle, alpha=one(T), beta=zero(T), o...)
    dx = similar(x)
    @cuda(cudnn,cudnnPoolingBackward,
          (Cptr,Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Ptr{T},Cptr,Ptr{T}),
          handle,PD(x;o...),Ref(alpha),TD(y),y,TD(dy),dy,TD(x),x,Ref(beta),TD(dx),dx)
    return dx
end

@primitive pool(x;o...),dy,y  poolx(x,y,dy;o...)
@zerograd  poolx(x,y,dy;o...)

# cudnn descriptors

type TD; ptr
    function TD(a::KnetArray)
        d = Cptr[0]
        @cuda(cudnn,cudnnCreateTensorDescriptor,(Ptr{Cptr},),d)
        n = ndims(a)
        sz = [Cint(size(a,n-i+1)) for i=1:n]
        st = [Cint(stride(a,n-i+1)) for i=1:n]
        @cuda(cudnn,cudnnSetTensorNdDescriptor,
              (Cptr,UInt32,Cint,Ptr{Cint},Ptr{Cint}),
              d[1], DT(a), n, sz, st)
        td = new(d[1])
        finalizer(td, x->@cuda(cudnn,cudnnDestroyTensorDescriptor,(Cptr,),x.ptr))
        return td
    end
end

type FD; ptr
    function FD(a::KnetArray)
        d = Cptr[0]
        @cuda(cudnn,cudnnCreateFilterDescriptor,(Ptr{Cptr},),d)
        n = ndims(a)
        sz = [Cint(size(a,n-i+1)) for i=1:n]
        @cuda(cudnn,cudnnSetFilterNdDescriptor,
              (Cptr,UInt32,UInt32,Cint,Ptr{Cint}),
              d[1], DT(a), 0,     n,   sz)
        fd = new(d[1])
        finalizer(fd, x->@cuda(cudnn,cudnnDestroyFilterDescriptor,(Cptr,),x.ptr))
        return fd
    end
end

type CD; ptr
    function CD(w::KnetArray,x::KnetArray; padding=0, stride=1, upscale=1, mode=0)
        d = Cptr[0]
        @cuda(cudnn,cudnnCreateConvolutionDescriptor,(Ptr{Cptr},),d)
        nd = ndims(x)-2
        @cuda(cudnn,cudnnSetConvolutionNdDescriptor,
              (Cptr,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},UInt32,UInt32),
              d[1],nd,cdsize(padding,nd),cdsize(stride,nd),cdsize(upscale,nd),mode,DT(x))
        cd = new(d[1])
        finalizer(cd, x->@cuda(cudnn,cudnnDestroyConvolutionDescriptor,(Cptr,),x.ptr))
        return cd
    end
end

type PD; ptr
    function PD(x::KnetArray; window=2, padding=0, stride=window, mode=0, maxpoolingNanOpt=0)
        d = Cptr[0]
        @cuda(cudnn,cudnnCreatePoolingDescriptor,(Ptr{Cptr},),d)
        nd = ndims(x)-2
        @cuda(cudnn,cudnnSetPoolingNdDescriptor,
              (Cptr,UInt32,UInt32,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint}),
              d[1],mode,maxpoolingNanOpt,nd,cdsize(window,nd),cdsize(padding,nd),cdsize(stride,nd))
        pd = new(d[1])
        finalizer(pd, x->@cuda(cudnn,cudnnDestroyPoolingDescriptor,(Cptr,),x.ptr))
        return pd
    end
end

import Base: unsafe_convert
unsafe_convert(::Type{Cptr}, td::TD)=td.ptr
unsafe_convert(::Type{Cptr}, fd::FD)=fd.ptr
unsafe_convert(::Type{Cptr}, cd::CD)=cd.ptr
unsafe_convert(::Type{Cptr}, pd::PD)=pd.ptr

function cdsize(w, nd)
    if isa(w,Integer)
        fill(Cint(w),nd)
    elseif length(w)==nd
        [ Cint(w[nd-i+1]) for i=1:nd ]
    else
        throw(DimensionMismatch())
    end
end

DT(::KnetArray{Float32})=UInt32(0)
DT(::KnetArray{Float64})=UInt32(1)
DT(::KnetArray{Float16})=UInt32(2)

function cdims{T,N}(w::KnetArray{T,N},x::KnetArray{T,N}; padding=0, stride=1, o...)
    ntuple(N) do i
        if i < N-1
            pi = (if isa(padding,Number); padding; else padding[i]; end)
            si = (if isa(stride,Number); stride; else stride[i]; end)
            1 + div(size(x,i) - size(w,i) + 2*pi, si)
        elseif i == N-1
            size(w,N)
        else # i == N
            size(x,N)
        end
    end
end

function pdims{T,N}(x::KnetArray{T,N}; window=2, padding=0, stride=window, o...)
    ntuple(N) do i
        if i < N-1
            wi = (if isa(window,Number); window; else window[i]; end)
            pi = (if isa(padding,Number); padding; else padding[i]; end)
            si = (if isa(stride,Number); stride; else stride[i]; end)
            1 + ceil(Int, (size(x,i) + 2*pi - wi) / si)
        else
            size(x,i)
        end
    end
end

"""

mat(x) reshapes x into a two-dimensional matrix.  For 1-D inputs mat
returns `reshape(x, (length(x),1))`.  For inputs with more than two
dimensions of size (X1,X2,...,XD), mat returns

    reshape(x, (X1*X2*...*X[D-1],XD))

"""
function mat(x)
    if ndims(x) > 2
        xn = size(x,ndims(x))
        reshape(x, (div(length(x),xn),xn))
    elseif ndims(x)==2
        x
    elseif ndims(x)==1
        reshape(x, (length(x),1))
    else
        throw(MethodError(mat,x))
    end
end


text:
export PIBasis

# ---------------------- Implementation of the PIBasisSpec

"""
`struct PIBasisSpec`
"""
struct PIBasisSpec
   orders::Vector{Int}     # order (length) of ith basis function
   iAA2iA::Matrix{Int}     # where in A can we find the ith basis function
end

==(B1::PIBasisSpec, B2::PIBasisSpec) = _allfieldsequal(B1, B2)

Base.length(spec::PIBasisSpec) = length(spec.orders)

maxcorrorder(spec::PIBasisSpec) = size(spec.iAA2iA, 2)

function _get_pibfcn(spec0, Aspec, vv)
   vv1 = vv[2:end]
   vvnz = vv1[findall(vv1 .!= 0)]
   return (spec0[vv[1]], Aspec[vvnz])
end

function _get_pibfcn(Aspec, vv)
   vvnz = vv[findall(vv .!= 0)]
   return Aspec[vvnz]
end

# TODO: maybe instead of property == nothing, there should be a
#       generic property with no symmetry attached to it.

function PIBasisSpec( basis1p::OneParticleBasis,
                      symgrp::SymmetryGroup,
                      Bsel::DownsetBasisSelector;
                      property = nothing,
                      filterfun = _->true,
                      init1pbasis = true )

   # we initialize the 1p-basis here; to prevent this it must be manually
   # avoided by passing in init1pbasis = false
   if init1pbasis
      init1pspec!(basis1p, Bsel)
   end

   # get the basis spec of the one-particle basis
   #  Aspec[i] described the basis function that will get written into A[i]
   Aspec = get_spec(basis1p)

   # we assume that `Aspec` is sorted by degree, but best to double-check this
   # since the notion of degree used to construct `Aspec` might be different
   # from the one used to construct AAspec.
   if !issorted(Aspec; by = b -> level(b, Bsel, basis1p))
      error("""PIBasisSpec : AAspec construction failed because Aspec is not
               sorted by degree. This could e.g. happen if an incompatible
               notion of degree was used to construct the 1-p basis spec.""")
   end
   # An AA basis function is given by a tuple 𝒗 = vv. Each index 𝒗ᵢ = vv[i]
   # corresponds to the basis function Aspec[𝒗ᵢ] and the tuple
   # 𝒗 = (𝒗₁, ...) to a product basis function
   #   ∏ A_{vₐ}
   tup2b = vv -> _get_pibfcn(Aspec, vv)

   #  degree or level of a basis function ↦ is it admissible?
   admissible = bb -> (level(bb, Bsel, basis1p) <= maxlevel(bb, Bsel, basis1p))

   if property != nothing
      filter1 = bb -> filterfun(bb) && filter(bb, Bsel, basis1p) && filter(property, symgrp, bb)
   else
      filter1 = bb -> filterfun(bb) && filter(bb, Bsel, basis1p)
   end

   # we can now construct the basis specification; the `ordered = true`
   # keyword signifies that this is a permutation-invariant basis
   maxord = maxorder(Bsel)
   AAspec = gensparse(; NU = maxorder(Bsel),
                        tup2b = tup2b,
                        admissible = admissible,
                        ordered = true,
                        maxvv = [length(Aspec) for _=1:maxord],
                        filter = filter1)

   return PIBasisSpec(AAspec)
end

function PIBasisSpec(AAspec)
   orders = zeros(Int, length(AAspec))
   iAA2iA = zeros(Int, (length(AAspec), length(AAspec[1])))
   for (iAA, vv) in enumerate(AAspec)
      # we use reverse because gensparse constructs the indices in
      # ascending order, but we want descending here.
      # (I don't remember why though)
      iAA2iA[iAA, :] .= reverse(vv)
      orders[iAA] = length( findall( vv .!= 0 ) )
   end
   return PIBasisSpec(orders, iAA2iA)
end

get_spec(AAspec::PIBasisSpec, i::Integer) = AAspec.iAA2iA[i, 1:AAspec.orders[i]]

# ------------------ PISpec sparsification

"""
returns a new `PIBasisSpec` constructed from the old one, but keeping only
the basis indices `Ikeep`. This maintains the order of the basis functions.
"""
sparsify(spec::PIBasisSpec, Ikeep::AbstractVector{<: Integer}) =
      PIBasisSpec(spec.orders[Ikeep], spec.iAA2iA[Ikeep, :])

function _fix_A_indices!(spec::PIBasisSpec, new_inds::AbstractVector{<: Integer})
   for iAA = 1:size(spec.iAA2iA, 1)
      for α = 1:spec.orders[iAA]
         vα = spec.iAA2iA[iAA, α]
         new_vα = new_inds[vα]
         @assert new_vα > 0
         spec.iAA2iA[iAA, α] = new_vα
      end
   end
   return nothing
end

# --------------------------------- PIBasis implementation

"""
`mutable struct PIBasis:` implementation of a permutation-invariant
basis based on the density projection trick.

The standard constructor is
```
PIBasis(basis1p, N, D, maxdeg)
```
* `basis1p` : a one-particle basis
* `N` : maximum interaction order
* `D` : an abstract degee specification, e.g., SparsePSHDegree
* `maxdeg` : the maximum polynomial degree as measured by `D`
"""
mutable struct PIBasis{BOP, REAL, TB, TA} <: ACEBasis
   basis1p::BOP             # a one-particle basis
   spec::PIBasisSpec
   real::REAL     # could be `real` or `identity` to keep AA complex
   # evaluator    # classic vs graph
   B_pool::VectorPool{TB}
   dAA_pool::VectorPool{TA}
end

cutoff(basis::PIBasis) = cutoff(basis.basis1p)

==(B1::PIBasis, B2::PIBasis) =
      ( (B1.basis1p == B2.basis1p) &&
        (B1.spec == B2.spec) &&
        (B1.real == B2.real) )

valtype(basis::PIBasis) = basis.real( valtype(basis.basis1p) )

valtype(basis::PIBasis, cfg::AbstractConfiguration) =
      basis.real( valtype(basis.basis1p, cfg) )

gradtype(basis::PIBasis, cfgorX) =
      basis.real( gradtype(basis.basis1p, cfgorX) )

Base.length(basis::PIBasis) = length(basis.spec)

# default symmetry group
PIBasis(basis1p, Bsel::AbstractBasisSelector; kwargs...) =
   PIBasis(basis1p, O3(), Bsel; kwargs...)

PIBasis(basis1p, symgrp, Bsel::AbstractBasisSelector;
        isreal = false, kwargs...) =
   PIBasis(basis1p,
           PIBasisSpec(basis1p, symgrp, Bsel; kwargs...),
           isreal ? Base.real : Base.identity )

function PIBasis(basis1p::OneParticleBasis, spec::PIBasisSpec, real)
   VT1 = valtype(basis1p)
   VT = real(VT1)  # default valtype
   B_pool = VectorPool{VT}()
   dAA_pool = VectorPool{VT1}()
   return PIBasis(basis1p, spec, real, B_pool, dAA_pool)
end

get_spec(pibasis::PIBasis) =
   [ get_spec(pibasis, i) for i = 1:length(pibasis) ]

get_spec(pibasis::PIBasis, i::Integer) =
      get_spec.( Ref(pibasis.basis1p), get_spec(pibasis.spec, i) )

setreal(basis::PIBasis, isreal::Bool) =
   PIBasis(basis.basis1p, basis.spec, isreal)

maxcorrorder(basis::PIBasis) = maxcorrorder(basis.spec)

# ------------------ sparsification

function sparsify!(basis::PIBasis, Ikeep::AbstractVector{<: Integer})
   basis.spec = sparsify(basis.spec, Ikeep)
   return basis
end

"""
This should allow the 1p basis to sparsify itself, then feed back to the
pibasis what the correct indices are.
"""
function clean_1pbasis!(basis::PIBasis)
   spec = get_spec(basis)
   B1p = basis.basis1p
   spec1p = NamedTuple[]
   for bb in spec
      append!(spec1p, bb)
   end
   identity.(unique!(spec1p))
   # sparsify the product 1p basis
   _, new_inds = sparsify!(basis.basis1p, spec1p)
   # now fix the indexing of the PIBasis specification
   _fix_A_indices!(basis.spec, new_inds)
   return basis
end

# syms = symbols(B1p)
# rgs = Dict{Symbol, Any}([sym => [] for sym in syms]...)
# for bb in spec, b in bb, sym in keys(b)
#    push!(rgs[sym], getproperty(b, sym))
# end
# for sym in syms
#    rgs[sym] = identity.(unique(rgs[sym]))
# end

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

_scaling_absvalue(x::Number) = abs(x)
_scaling_absvalue(x::Symbol) = 0

# TODO: this is a hack; cf. #68
function scaling(pibasis::PIBasis, p)
   ww = zeros(Float64, length(pibasis))
   bspec = get_spec(pibasis)
   for i = 1:length(pibasis)
      for b in bspec[i]
         # TODO: revisit how this should be implemented for a general basis
         ww[i] += sum(x -> _scaling_absvalue(x)^p, b)  #  abs.(values(b)).^p
      end
   end
   return ww
end

# function scaling(pibasis::PIBasis, p)
#    ww = zeros(Float64, length(pibasis))
#    for iz0 = 1:numz(pibasis)
#       wwin = @view ww[pibasis.inner[iz0].AAindices]
#       for i = 1:length(pibasis.inner[iz0])
#          bspec = get_basis_spec(pibasis, iz0, i)
#          wwin[i] = scaling(bspec, p)
#       end
#    end
#    return ww
# end

# graphevaluator(basis::PIBasis) =
#    PIBasis(basis.basis1p, zlist(basis), basis.inner, DAGEvaluator())
#
# standardevaluator(basis::PIBasis) =
#    PIBasis(basis.basis1p, zlist(basis), basis.inner, StandardEvaluator())

# -------------------------------------------------
# FIO codes

write_dict(basis::PIBasis) =
   Dict(  "__id__" => "ACE_PIBasis",
         "basis1p" => write_dict(basis.basis1p),
            "spec" => write_dict(basis.spec),
            "real" => basis.real == Base.real ? true : false )

read_dict(::Val{:ACE_PIBasis}, D::Dict) =
   PIBasis( read_dict(D["basis1p"]),
            read_dict(D["spec"]),
            D["real"] ? Base.real : Base.identity )

write_dict(spec::PIBasisSpec) =
   Dict( "__id__" => "ACE_PIBasisSpec",
         "orders" => spec.orders,
         "iAA2iA" => write_dict(spec.iAA2iA) )

read_dict(::Val{:ACE_PIBasisSpec}, D::Dict) =
   PIBasisSpec( D["orders"], read_dict(D["iAA2iA"]) )

# -------------------------------------------------
# Evaluation codes

function evaluate!(AA, basis::PIBasis, config::AbstractConfiguration)
   A = acquire_B!(basis.basis1p, config)   #  THIS ALLOCATES!!!!
   evaluate!(A, basis.basis1p, config)
   fill!(AA, 1)
   for iAA = 1:length(basis)
      aa = one(eltype(A))
      for t = 1:basis.spec.orders[iAA]
         aa *= A[ basis.spec.iAA2iA[ iAA, t ] ]
      end
      AA[iAA] = basis.real(aa)
   end
   release_B!(basis.basis1p, A)
   return AA
end

# -------------------------------------------------
# gradients

ACE.evaluate_d(basis::PIBasis, cfg::AbstractConfiguration, args...) =
         ACE.evaluate_ed(basis::PIBasis, cfg::AbstractConfiguration, args...)[2]

function evaluate_ed!(AA, dAA, basis::PIBasis,
                      cfg::AbstractConfiguration, args...)
   A = acquire_B!(basis.basis1p, cfg)
   dA = acquire_dB!(basis.basis1p, cfg)   # TODO: THIS WILL ALLOCATE!!!!!
   evaluate_ed!(A, dA, basis.basis1p, cfg, args...)
   evaluate_ed!(AA, dAA, basis, A, dA)
   release_dB!(basis.basis1p, dA)
   release_B!(basis.basis1p, A)
   return AA, dAA
end

function _AA_local_adjoints!(dAAdA, A, iAA2iA, iAA, ord, _real)
   if ord == 1
      return _AA_local_adjoints_1!(dAAdA, A, iAA2iA, iAA, ord, _real)
   elseif ord == 2
      return _AA_local_adjoints_2!(dAAdA, A, iAA2iA, iAA, ord, _real)
   # elseif ord == 3
   #    return _AA_local_adjoints_3!(dAAdA, A, iAA2iA, iAA, ord, _real)
   # elseif ord == 4
   #    return _AA_local_adjoints_4!(dAAdA, A, iAA2iA, iAA, ord, _real)
   else
      return _AA_local_adjoints_x!(dAAdA, A, iAA2iA, iAA, ord, _real)
   end
end

function _AA_local_adjoints_1!(dAAdA, A, iAA2iA, iAA, ord, _real)
   @inbounds dAAdA[1] = 1
   @inbounds A1 = A[iAA2iA[iAA, 1]]
   return _real(A1)
end

function _AA_local_adjoints_2!(dAAdA, A, iAA2iA, iAA, ord, _real)
   @inbounds A1 = A[iAA2iA[iAA, 1]]
   @inbounds A2 = A[iAA2iA[iAA, 2]]
   @inbounds dAAdA[1] = A2
   @inbounds dAAdA[2] = A1
   return _real(A1 * A2)
end

# function _AA_local_adjoints_3!(dAAdA, A, iAA2iA, iAA, ord, _real)
#    A1 = A[iAA2iA[iAA, 1]]
#    A2 = A[iAA2iA[iAA, 2]]
#    A3 = A[iAA2iA[iAA, 3]]
#    A12 = A1 * A2
#    dAAdA[1] = A2 * A3
#    dAAdA[2] = A1 * A3
#    dAAdA[3] = A12
#    return _real(A12 * A3)
# end

# function _AA_local_adjoints_4!(dAAdA, A, iAA2iA, iAA, ord, _real)
#    @inbounds A1 = A[iAA2iA[iAA, 1]]
#    @inbounds A2 = A[iAA2iA[iAA, 2]]
#    @inbounds A3 = A[iAA2iA[iAA, 3]]
#    @inbounds A4 = A[iAA2iA[iAA, 4]]
#    A12 = A1 * A2
#    A34 = A3 * A4
#    @inbounds dAAdA[1] = A2 * A34
#    @inbounds dAAdA[2] = A1 * A34
#    @inbounds dAAdA[3] = A12 * A4
#    @inbounds dAAdA[4] = A12 * A3
#    return _real(A12 * A34)
# end

function _AA_local_adjoints_x!(dAAdA, A, iAA2iA, iAA, ord, _real)
   @assert length(dAAdA) >= ord
   @assert ord >= 2
   # TODO - optimize a bit more? can move one operation out of the loop
   # Forward pass:
   @inbounds A1 = A[iAA2iA[iAA, 1]]
   @inbounds A2 = A[iAA2iA[iAA, 2]]
   @inbounds dAAdA[1] = 1
   @inbounds dAAdA[2] = A1
   @inbounds AAfwd = A1 * A2
   @inbounds for a = 3:ord-1
      dAAdA[a] = AAfwd
      AAfwd *= A[iAA2iA[iAA, a]]
   end
   @inbounds dAAdA[ord] = AAfwd
   @inbounds Aend = A[iAA2iA[iAA, ord]]
   aa = _real(AAfwd * Aend)
   # backward pass
   @inbounds AAbwd = Aend
   @inbounds for a = ord-1:-1:3
      dAAdA[a] *= AAbwd
      AAbwd *= A[iAA2iA[iAA, a]]
   end
   dAAdA[2] *= AAbwd
   AAbwd *= A2
   dAAdA[1] *= AAbwd

   return aa
end

_acquire_dAAdA!(basis::PIBasis) = acquire!(basis.dAA_pool, maxcorrorder(basis))

function evaluate_ed!(AA, dAA, basis::PIBasis,
                      A::AbstractVector, dA::AbstractMatrix)
   orders = basis.spec.orders
   iAA2iA = basis.spec.iAA2iA
   dAAdA = _acquire_dAAdA!(basis)

   # Must treat the constants separately. This is not so elegant and could
   # maybe be improved?
   if orders[1] == 0  # SHOULD BE THE ONLY ONE with ord=0!!
      iAAinit = 2
      AA[1] = 1.0
      dAA[1, :] .= Ref(zero(eltype(dAA)))
   else
      iAAinit = 1
   end

   for iAA = iAAinit:length(basis)
      ord = orders[iAA]

      # ----- compute the local adjoints dAA / dA
      # dAAdA[a] ← ∏_{t ≂̸ a} A_{v_t}
      AA[iAA] = _AA_local_adjoints!(dAAdA, A, iAA2iA, iAA, orders[iAA], basis.real)

      # ----- now convert them into dAA / dX
      for j = 1:size(dA, 2)
         dAA[iAA, j] = sum(dAAdA[a] * dA[iAA2iA[iAA, a], j]
                           for a = 1:ord) |> basis.real
      end
   end

   return AA, dAA
end

ext: jl


text:
#!/usr/bin/env julia

using Luxor, Colors
width, height = 1000, 1000
fname = "/tmp/subpath-clip-image.pdf"
Drawing(width, height, fname)
origin()
background("ivory")
sethue("black")
image = readpng("../figures/sector-test.png")
w = image.width
h = image.height

# all these graphics will make up the clipping mask:
# reversepath makes holes in the path
b = false
for i in 3:15
    randomhue()
    poly(ngon(0, 0, 50 + (i * 30), i, pi/2, vertices=true), :path, reversepath=b)
    if i < 15
        newsubpath()
    end
    b  = !b
end

# clip them
clippreserve()

# draw image 'through' clipping mask
placeimage(image, -w/2, -h/2)

clipreset()

finish()
println("finished test: output in $(fname)")


text:
const _pkg_root = dirname(dirname(@__FILE__))
const _pkg_assets = joinpath(_pkg_root,"assets")

!isdir(_pkg_assets) && mkdir(_pkg_assets)

deps = [
    "https://cdn.jsdelivr.net/gh/piever/InteractResources@0.4.0/bulma/main.min.css" => "bulma.min.css",
    "https://cdn.jsdelivr.net/gh/piever/InteractResources@0.4.0/bulma/main_confined.min.css" => "bulma_confined.min.css",
]

for (dep, name) in deps
    download(dep, joinpath(_pkg_assets, name))
end


text:
# This file defines the functions pertaining to the ManpowerSimulation type.

# The functions of the ManpowerSimulation type require SimJulia,
#   ResumableFunctions, and Distributions.

# The functions of the ManpowerSimulation type require all types.
requiredTypes = [ "recruitment",
    "state",
    "transition",
    "retirement",
    "attrition",
    "simulationReport",
    "manpowerSimulation" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

# This function tests if the simulation has been properly initialised.
export isInitialised
function isInitialised( mpSim::ManpowerSimulation )

    return mpSim.isInitialised

end  # isInitialised( mpSim )

# This function sets the ID key in the simulation to the given attribute. This
#   does NOT update the databases.
function setKey( mpSim::ManpowerSimulation, id::Union{Symbol, String} = "id" )

    mpSim.idKey = String( id )

end  # setKey!( mpSim, id )

include( joinpath( funcPath, "parFileProcessing.jl" ) )
include( joinpath( funcPath, "sqLiteToSim.jl" ) )
include( joinpath( funcPath, "configToSQLite.jl" ) )

# This function sets the cap on the number of personnel in the simulation. Note
#   that this cap can be (temporarily) violated if the initial manpower force is
#   larger than the personnel cap. If this function receives a value ⩽ 0, there
#   will not be a personnel cap.
export setPersonnelCap
function setPersonnelCap( mpSim::ManpowerSimulation, cap::T ) where T <: Integer

    mpSim.personnelTarget = cap > 0 ? cap : 0

end  # setPersonnelCap( mpSim, cap )

# This function sets the start date of the simulation.
export setSimStartDate
function setSimStartDate( mpSim::ManpowerSimulation, startDate::Date )::Void

    mpSim.simStartDate = startDate
    return

end

# This function sets the time between two database commits.
export setDatabaseCommitTime
function setDatabaseCommitTime( mpSim::ManpowerSimulation, commitFreq::T ) where T <: Real

    if commitFreq <= 0.0
        warn( "Time between two database commits must be > 0.0." )
        return
    end  # if commitFreq <= 0.0

    mpSim.commitFrequency = commitFreq

end  # setDatabaseCommitTime( mpSim, commitFreq )

# This function adds a personnel attribute to the manpower simulation. If an
#   attribute with the same name already exists in the simulation, the atrtibute
#   will be overwritten and a warning wil be issued.
# This function only works if the simulation hasn't started yet.
export addAttribute!
function addAttribute!( mpSim::ManpowerSimulation,
    attr::PersonnelAttribute )::Void

    if now( mpSim ) > 0
        return
    end  # if now( mpSim ) > 0

    isInit = !isempty( attr.values )
    listToCheck = isInit ? mpSim.initAttrList : mpSim.otherAttrList
    indexOfAttr = findfirst( tmpAttr -> tmpAttr.name == attr.name, listToCheck )

    if indexOfAttr > 0
        warn( "Attribute with name '$(attr.name)' already exists, overwriting the attribute." )
        listToCheck[ indexOfAttr ] = attr
    else
        push!( listToCheck, attr )
    end  # if indexOfAttr > 0

    return

end  # addAttribute!( mpSim, attr )

# This function clears the list of personnel attributes.
# This function only works if the simulation is in a virgin state.
export clearAttributes!
function clearAttributes!( mpSim::ManpowerSimulation )

    if now( mpSim ) == 0
        empty!( mpSim.initAttrList )
        empty!( mpSim.otherAttrList )
    end  # if now( mpSim ) == 0

    return

end  # clearAttributes!( mpSim )

# This function adds a possible personnel state.
# This function does nothing if the simulation is already running.
export addState!
function addState!( mpSim::ManpowerSimulation, state::State,
    isInitial::Bool = false )

    if now( mpSim ) == 0
        if isInitial
            mpSim.initStateList[ state ] = Vector{Transition}()
        else
            mpSim.otherStateList[ state ] = Vector{Transition}()
        end  # if isinitial

        mpSim.stateList[ state.name ] = state
    end  # if now( mpSim ) == 0

    return

end  # addState!( mpSim )

export removeState!
function removeState!( mpSim::ManpowerSimulation, stateName::String )::Void

    # Nothing to remove if state doesn't exist.
    if !haskey( mpSim.stateList, stateName )
        return
    end  # if !haskey( mpSim.stateList, stateName )

    # Clear the state from all compound states it's a part of.
    for compStateName in keys( mpSim.compoundStateList )
        compState = mpSim.compoundStateList[ compStateName ]
        stateInd = findfirst( stateName .== compState.stateList )

        if stateInd != 0
            deleteat!( compState.stateList, stateInd )
        end  # if stateInd != 0
    end  # for compStateName in keys( mpSim.compoundStateList )

    state = mpSim.stateList[ stateName ]
    isInitial = state.isInitial
    delete!( isInitial ? mpSim.initStateList : mpSim.otherStateList, state )
    delete!( mpSim.stateList, stateName )
    return

end  # removeState!( mpSim, stateName )

# This function clears the list of possible personnel states.
# This function does nothing if the simulation is already running.
export clearStates!
function clearStates!( mpSim::ManpowerSimulation )

    if now( mpSim ) == 0
        empty!( mpSim.compoundStateList )
        empty!( mpSim.initStateList )
        empty!( mpSim.otherStateList )
        empty!( mpSim.stateList )
        empty!( mpSim.transList )
    end  # if now( mpSim ) == 0

    return

end  # clearStates!( mpSim, state, isInitial )

export addCompoundState!,
       removeCompoundState!,
       clearCompoundStates!
"""
```
addCompoundState!( mpSim::ManpowerSimulation,
                   compState::CompoundState )
```
This function adds the compound state `compState` to the manpower simulation
`mpSim`. Ensure that all the states comprising the compound exist in the
simulation as non-existent states will be removed from the compound state's
composition. If a compound state with the same name already exists in the
simulation, it will be overwritten.

This function returns `nothing`.
"""
function addCompoundState!( mpSim::ManpowerSimulation,
    compState::CompoundState )::Void

    # Filter out the states that aren't defined.
    stateInd = map( stateName -> !haskey( mpSim.stateList, stateName ),
        compState.stateList )
    deleteat!( compState.stateList, stateInd )

    if haskey( mpSim.compoundStateList, compState.name )
        warn( "Compound state with name '$(compState.name)' already exists. Overwriting." )
    end  # if haskey( mpSim.compoundStateList, compState.name )

    mpSim.compoundStateList[ compState.name ] = compState
    return

end  # addCompoundState!( mpSim, compState )

"""
```
addCompoundState!( mpSim::ManpowerSimulation,
                   compName::String,
                   target::T,
                   stateList::String... )
    where T <: Integer
```
This function adds the compound state with name `compName`, personnel target
`target`, and composed of the states in `stateList` to the manpower simulation
`mpSim`. Ensure that all the states comprising the compound exist in the
simulation as non-existent states will be removed from the compound state's
composition. If a compound state with the same name already exists in the
simulation, it will be overwritten.

This function returns `nothing`.
"""
function addCompoundState!( mpSim::ManpowerSimulation, compName::String,
    target::T, stateList::String... )::Void where T <: Integer

    newCompState = CompoundState( compName )
    setStateTarget!( newCompState, target )
    addStateToCompound!( newCompState, stateList... )
    addCompoundState!( mpSim, newCompState )
    return

end  # addCompoundState!( mpSim, compName, target, stateList )

"""
```
removeCompoundState!( mpSim::ManpowerSimulation,
                      stateList::String... )
```
This function removes the compound states with names in `stateList` from the
manpower simulation `mpSim`.

This function returns `nothing`.
"""
function removeCompoundState!( mpSim::ManpowerSimulation,
    stateList::String... )::Void

    foreach( stateList ) do
        delete!( mpSim.compoundStateList, stateName )
        delete!( mpSim.compoundStates, stateName )
    end  # foreach( stateList ) do

    return

end  # removeCompoundState!( mpSim, stateList )

"""
```
clearCompoundStates!( mpSim::ManpowerSimulation )
```
This function removes all compound states from the manpower simulation `mpSim`.

This function returns `nothing`.
"""
function clearCompoundStates!( mpSim::ManpowerSimulation )::Void

    empty!( mpSim.compoundStates )
    empty!( mpSim.compoundStateList )
    return

end  # clearCompoundStates!( mpSim )

export setTransTypePriority!
function setTransTypePriority!( mpSim::ManpowerSimulation, transName::String,
    prio::Int )::Void

    mpSim.transList[ transName ] = prio
    return

end  # setTransTypePriority!( mpSim, transName, prio )

# This function adds a personnel state transition.
# This function does nothing if the simulation is already running.
export addTransition!
function addTransition!( mpSim::ManpowerSimulation, trans::Transition )

    if now( mpSim ) == 0
        if trans.startState.isInitial
            push!( mpSim.initStateList[ trans.startState ], trans )
        else
            push!( mpSim.otherStateList[ trans.startState ], trans )
        end  # if isStartInit

        if !haskey( mpSim.transList, trans.name )
            mpSim.transList[ trans.name ] = 0
        end  # if !haskey( mpSim.transList, trans.name )
    end  # if now( mpSim ) == 0

    return

end  # addTransition!( mpSim, state )

# This function clears the the personnel state transitions.
# This function does nothing if the simulation is already running.
export clearTransitions!
function clearTransitions!( mpSim::ManpowerSimulation )

    if now( mpSim ) == 0
        foreach( state -> empty!( mpSim.initStateList[ state ] ),
            keys( mpSim.initStateList ) )
        foreach( state -> empty!( mpSim.otherStateList[ state ] ),
            keys( mpSim.otherStateList ) )
        empty!( mpSim.transList )
    end  # if now( mpSim ) == 0

    return

end  # clearTransitions!( mpSim )

function processTransPriorities( mpSim::ManpowerSimulation )::Void

    transNames = collect( keys( mpSim.transList ) )
    prios = map( transName -> mpSim.transList[ transName ], transNames )
    prios[ prios .<= 0 ] = maximum( prios ) + 1
    orderIndex = sortperm( prios )
    transNames = transNames[ orderIndex ]
    foreach( ii -> mpSim.transList[ transNames[ ii ] ] = ii,
        eachindex( transNames ) )
    return

end  # processTransPriorities( mpSim )

function assignTransPriorities( mpSim::ManpowerSimulation )::Void

    stateList = orderStates( mpSim )
    nTypes = length( mpSim.transList )
    transPrio = 1

    for ii in eachindex( stateList )
        state = mpSim.stateList[ stateList[ ii ] ]
        isInitial = haskey( mpSim.initStateList, state )
        transList = isInitial ? mpSim.initStateList[ state ] :
            mpSim.otherStateList[ state ]
        transPrios = similar( transList, Int )

        for jj in eachindex( transList )
            trans = transList[ jj ]
            transPrios[ jj ] = trans.hasPriority ? 0 :
                mpSim.transList[ trans.name ]
        end  # for jj in eachindex( transList )

        orderIndex = sortperm( transPrios )

        for jj in eachindex( transList )
            kk = orderIndex[ jj ]
            trans = transList[ kk ]
            trans.transPriority = trans.hasPriority ? 0 : transPrio

            if !trans.hasPriority
                transPrio += 1
            end  # if !trans.hasPriority
        end  # for jj in eachindex( transList )
    end  # for state in stateList

    return

end

# This function adds a recruitment scheme to the simulation.
# This function does nothing if the simulation is already running. XXX (desirable??)
export addRecruitmentScheme!
function addRecruitmentScheme!( mpSim::ManpowerSimulation,
    recScheme::Recruitment )

    if now( mpSim ) == 0
        targetState = recScheme.recState

        # Don't do anything if the recruitment state isn't defined.
        if ( targetState != "" ) && !haskey( mpSim.stateList, targetState )
        else
            # Switch the recruitment state (if any) to an initial state if it
            #   isn't an initial state already.
            if ( targetState != "" ) &&
                !haskey( mpSim.initStateList, mpSim.stateList[ targetState ] )
                state = mpSim.stateList[ targetState ]
                removeState!( mpSim, targetState )
                state.isInitial = true
                addState!( mpSim, state, true )
            end  # if ( targetState != "" ) && ...

            push!( mpSim.recruitmentSchemes, recScheme )
        end  # if ( targetState != "" ) &&
    end  # if now( mpSim ) == 0

    return

end  # addRecruitmentScheme!( mpSim, recScheme )

# This function clears the recruitment schemes from the simulation.
# This function should NOT be called while running a simulation!
export clearRecruitmentSchemes!
function clearRecruitmentSchemes!( mpSim::ManpowerSimulation )

    if now( mpSim ) == 0
        empty!( mpSim.recruitmentSchemes )
    end  # if now( mpSim ) == 0

    return

end  # clearRecruitmentSchemes!( mpSim )

# This function sets the retirement scheme of the simulation.
export setRetirement
function setRetirement( mpSim::ManpowerSimulation,
    retScheme::Union{Void, Retirement} = nothing )

    mpSim.retirementScheme = retScheme

end  # setRetirement( mpSim, retScheme )

# This function sets the attrition scheme.
export setAttrition
function setAttrition( mpSim::ManpowerSimulation,
    attrScheme::Union{Void, Attrition} = nothing )

    mpSim.defaultAttritionScheme = attrScheme

end  # setAttrition( mpSim, attrScheme )

# This function sets the length of the simulation.
export setSimulationLength
function setSimulationLength( mpSim::ManpowerSimulation, simLength::T ) where T <: Real

    if simLength <= 0.0
        warn( "Simulation length must be > 0.0" )
        return
    end  # if simLength <= 0

    mpSim.simLength = simLength

end  # setSimulationLength( mpSim, simLength )

export setPhasePriority
function setPhasePriority( mpSim::ManpowerSimulation, phase::Symbol,
    priority::T ) where T <: Integer

    if ( phase ∉ [ :recruitment, :retirement, :attrition, :attrCheck ] ) &&
        ( phase ∉ mpSim.workingDbase.attrs )
        warn( "Unknown simulation phase, not setting priority." )
        return
    end

    mpSim.phasePriorities[ phase ] = priority

end  # setPhasePriority( mpSim, phase, priority )

export resetSimulation
function resetSimulation( mpSim::ManpowerSimulation )

    mpSim.sim = Simulation()
    mpSim.simTimeElapsed = Dates.Millisecond( 0 )

    # First, drop the tables with the same name if they exist.
    # XXX The personnel table must be dropped last due to the FOREIGN KEY
    #   constraints in the other tables.
    SQLite.drop!( mpSim.simDB, mpSim.historyDBname, ifexists = true )
    SQLite.drop!( mpSim.simDB, mpSim.transitionDBname, ifexists = true )
    SQLite.drop!( mpSim.simDB, mpSim.personnelDBname, ifexists = true )

    # Then, create the personnel and history databases.
    # XXX Other defined attributes need to be introduced here as well.
    command = "CREATE TABLE $(mpSim.personnelDBname)(
        $(mpSim.idKey) varchar(16) NOT NULL PRIMARY KEY,
        timeEntered float,
        timeExited float,
        ageAtRecruitment float,
        expectedRetirementTime float,
        expectedAttritionTime float,
        attritionScheme varchar(64)"

    if !( isempty( mpSim.initAttrList ) && isempty( mpSim.otherAttrList ) )
        command *= ",\n" * join( map( attr -> "'" * attr.name * "' varcar(64)",
            vcat( mpSim.initAttrList, mpSim.otherAttrList ) ), ",\n" )
    end  # if !( isempty( mpSim.initAttrList ) &&

    command *= ",\nstatus varchar(16) )"
    SQLite.execute!( mpSim.simDB, command )

    command = "CREATE TABLE $(mpSim.historyDBname)(
        $(mpSim.idKey) varchar(16),
        attribute varchar(255),
        timeIndex float,
        numValue float,
        strValue varchar(255),
        FOREIGN KEY ($(mpSim.idKey)) REFERENCES $(mpSim.personnelDBname)($(mpSim.idKey))
    )"
    SQLite.execute!( mpSim.simDB, command )

    command = "CREATE TABLE $(mpSim.transitionDBname)(
        $(mpSim.idKey) varchar(16),
        timeIndex float,
        transition varchar(255),
        startState varchar(255),
        endState varchar(255),
        FOREIGN KEY ($(mpSim.idKey)) REFERENCES $(mpSim.personnelDBname)($(mpSim.idKey))
    )"
    SQLite.execute!( mpSim.simDB, command )

    mpSim.personnelSize = 0
    mpSim.resultSize = 0

    for state in keys( merge( mpSim.initStateList, mpSim.otherStateList ) )
        empty!( state.inStateSince )
        empty!( state.isLockedForTransition )
    end  # for state in keys( merge(

    # And wipe all existing simulation reports.
    # empty!( mpSim.simReports )

    mpSim.isVirgin = true

end  # resetSimulation( mpSim )

# This function generates a random initial population for the simulation. If the
#   cap is < 0, the initial population satisfies the cap in the simulation.
# If the simulation has already started, this will reset the simulation!
export populate
function populate( mpSim::ManpowerSimulation, cap::T = 0 ) where T <: Integer
    # Throw an error if a populate to cap is requested for a simulation without
    #   personnel cap.
    if ( cap < 0 ) && ( mpSim.personnelTarget == 0 )
        error( "Simulation has no cap, cannot populate to cap." )
    end  # if ( cap < 0 ) ...

    # (Re)initialise the simulation and the population.
    resetSimulation( mpSim )

    # If the initial population is zero, stop.
    if cap == 0
        return
    end  # if cap == 0

    # Set the actual number of personnel to seed.
    tmpCap = cap < 0 ? mpSim.personnelTarget : cap
    isRetirement = mpSim.retirementScheme !== nothing
    timeDist = isRetirement ?
        Uniform( 0.0, mpSim.retirementScheme.maxCareerLength ) : nothing

    # Create the required number of persons.
    addPersonnel!( mpSim.workingDbase,
        map( ii -> "Init" * string( ii ), 1:tmpCap ) )
    addPersonnel!( mpSim.simResult,
        map( ii -> "Init" * string( ii ), 1:tmpCap ) )

    # Fill in the necessary attributes.
    for ii in 1:tmpCap
        person = mpSim.workingDbase[ ii ]
        result = mpSim.simResult[ ii ]
        person[ :status ] = :active
        result[ :history ] = Dict{Symbol, History}()
        result[ :history ][ :status ] = History( :status )

        # Generate time in system.
        timeInSystem = isRetirement ? rand( timeDist ) : 1.0

        # Compute moment of entry in system and create retirement process if
        #   necessary.
        person[ :timeEntered ] = -timeInSystem
        result[ :history ][ :status ][ person[ :timeEntered ] ] = :active

        if isRetirement
            person[ :processRetirement ] = @process retireProcess( mpSim.sim,
                person, result, mpSim )
        end  # if isRetirement
    end  # for ii in 1:tmpCap

    # Database has entries.
    mpSim.isVirgin = false
end

# This function (re)initialises the simulation.
export initialise
function initialise( mpSim::ManpowerSimulation,
    id::Union{Symbol, String} = "id" )

    setKey( mpSim, id )
    # setPhasePriority( mpSim, :recruitment, 1 )
    # setPhasePriority( mpSim, :retirement, 2 )
    # setPhasePriority( mpSim, :attrition, 3 )
    resetSimulation( mpSim )

    # The simulation is correctly initialised.
    mpSim.isInitialised = true

end  # initialise( mpSim, id, cap, recSchemes, retScheme, initPop )

# This function returns the current time of the manpower simulation.
function Dates.now( mpSim::ManpowerSimulation )

    return now( mpSim.sim )

end

# This file holds the database commit process.
include( joinpath( funcPath, "dbManagement.jl" ) )

export configureSimFromGraph
function configureSimFromGraph( mpSim::ManpowerSimulation, graphFile::String,
    showPlot::Bool = true,
    exportFileName::String = "networkExport.graphml" )::Void

    if !ispath( graphFile )
        error( "File '$graphFile' does not exist." )
    end  # if !ispath(  )

    resetSimulation( mpSim )
    empty!( mpSim.initAttrList )
    empty!( mpSim.otherAttrList )
    empty!( mpSim.initStateList )
    empty!( mpSim.otherStateList )

    # Get the graph
    xmlGraph = readxml( graphFile )
    fileRoot = root( xmlGraph )
    fileElements = elements( fileRoot )
    graphIndex = findfirst( el -> nodename( el ) == "graph", fileElements )
    graphRoot = fileElements[ graphIndex ]

    # Find the indices of the nodes (States) and edges (Transitions) in the
    #   graph element.
    graphElements = elements( graphRoot )
    nodeIndices = find( el -> nodename( el ) == "node", graphElements )
    edgeIndices = find( el -> nodename( el ) == "edge", graphElements )

    # Read all states.
    systemNetwork = MetaDiGraph( SimpleDiGraph() )
    nodeDict = Dict{String, String}()
    nNodes = 0

    for ii in nodeIndices
        nodeRoot = graphElements[ ii ]
        nodeID = nodeRoot[ "id" ]
        shapeIndex = findfirst( el -> el[ "key" ] == "d6", elements( nodeRoot ) )
        shapeNode = elements( nodeRoot )[ shapeIndex ]
        shapeNode = elements( shapeNode )[ 1 ]
        shapeElements = elements( shapeNode )
        labelIndex = findfirst( el -> nodename( el ) == "NodeLabel",
            shapeElements )
        nodeLabel = "State $ii"

        if labelIndex == 0
            warn( "Node has no label, setting it to $nodeLabel." )
        else
            nodeLabel = nodecontent( shapeElements[ labelIndex ] )
        end  # if labelIndex == 0

        if lowercase( nodeLabel ) ∉ [ "in", "external" ]
            nNodes += 1
            nodeDict[ nodeID ] = nodeLabel
            add_vertex!( systemNetwork )
            set_prop!( systemNetwork, nNodes, :state, nodeLabel )
        end  # if lowercase( nodeLabel ) != "in"
    end  # for ii in nodeIndices

    stateList = map( node -> get_prop( systemNetwork, node, :state ),
        vertices( systemNetwork ) )

    # Read all transitions.
    for ii in edgeIndices
        edgeRoot = graphElements[ ii ]
        edgeID = edgeRoot[ "id" ]
        sourceName = get( nodeDict, edgeRoot[ "source" ], "in" )
        targetName = get( nodeDict, edgeRoot[ "target" ], "in" )
        shapeIndex = findfirst( el -> el[ "key" ] == "d10",
            elements( edgeRoot ) )
        shapeNode = elements( edgeRoot )[ shapeIndex ]
        shapeNode = elements( shapeNode )[ 1 ]
        shapeElements = elements( shapeNode )
        labelIndex = findfirst( el -> nodename( el ) == "EdgeLabel",
            shapeElements )

        if "in" ∉ [ sourceName, targetName ]
            edgeLabel = "$(sourceName)_to_$(targetName)"

            if labelIndex == 0
                warn( "Edge has no label. Setting it to $edgeLabel" )
            else
                edgeLabel = nodecontent( shapeElements[ labelIndex ] )
            end  # if labelIndex == 0

            sourceVertex = findfirst( node ->
                get_prop( systemNetwork, node, :state ) == sourceName,
                vertices( systemNetwork ) )
            targetVertex = findfirst( node ->
                get_prop( systemNetwork, node, :state ) == targetName,
                vertices( systemNetwork ) )
            add_edge!( systemNetwork, sourceVertex, targetVertex )
            set_prop!( systemNetwork, sourceVertex, targetVertex, :trans,
                edgeLabel )
        end  # if "in" ∉ [ sourceName, targetName ]
    end

    stateDict = Dict{String, State}()

    # Add states to simulation.
    for stateNode in vertices( systemNetwork )
        isInitial = isempty( inneighbors( systemNetwork, stateNode ) )
        stateName = get_prop( systemNetwork, stateNode, :state )
        newState = State( stateName, isInitial)
        stateDict[ stateName ] = newState
        addState!( mpSim, newState, isInitial )
    end  # for stateNode in vertices( systemNetwork )

    # Add transitions to simulation.
    for transEdge in edges( systemNetwork )
        transName = get_prop( systemNetwork, transEdge, :trans )
        sourceName = stateList[ src( transEdge ) ]
        targetName = stateList[ dst( transEdge ) ]
        newTrans = Transition( transName, stateDict[ sourceName ],
            stateDict[ targetName ] )
        addTransition!( mpSim, newTrans )
    end  # for transEdge in edges( systemNetwork )

    if showPlot
        plotTransitionMap( mpSim, stateList..., fileName = exportFileName )
    end  # if showPlot

    return

end  # configureSimFromGraph( mpSim::ManpowerSimulation, graphFile::String )

# This file holds the functions to upload an initial population snapshot.
include( joinpath( funcPath, "snapshot.jl" ) )

# This function runs the manpower simulation if it has been properly
#   initialised.
# function SimJulia.run( mpSim::ManpowerSimulation, toTime::T = 0.0 ) where T <: Real
function SimJulia.run( mpSim::ManpowerSimulation )

    if !mpSim.isInitialised
        error( "Simulation not properly initialised. Cannot run." )
    end  # if !mpSim.isInitialised

    # Set up the recruitment processes and read initial population snapshot.
    if mpSim.isVirgin
        # We cannot write this with a map statement. The @process macro messes
        #   that up.
        for ii in eachindex( mpSim.recruitmentSchemes )
            @process recruitProcess( mpSim.sim, ii, mpSim )
        end  # for ii in eachindex( mpSim.recruitmentSchemes )

        mpSim.isWellDefined = true
        mpSim.isVirgin = false
        readSnapshot( mpSim )
    end  # if mpSim.isVirgin

    toTime = 0.0
    oldSimTime = now( mpSim )
    mpSim.attrExecTimeElapsed = Dates.Millisecond( 0 )

    # Process the compound states read from the catalogue.
    processCompoundStates( mpSim )

    # Start the database commits.
    SQLite.execute!( mpSim.simDB, "BEGIN TRANSACTION" )
    @process dbCommitProcess( mpSim.sim,
        toTime == 0.0 ? mpSim.simLength : toTime, mpSim )
    initiateTransitionProcesses( mpSim )
    @process retireProcess( mpSim.sim, mpSim )
    @process checkAttritionProcess( mpSim.sim, mpSim )
    startTime = now()

    if toTime > 0.0
        run( mpSim.sim, toTime )
    else
        run( mpSim.sim )
    end  # if toTime > 0.0

    # Final commit.
    SQLite.execute!( mpSim.simDB, "COMMIT" )
    saveSimConfigToDatabase( mpSim )
    mpSim.simTimeElapsed += now() - startTime
    println( "Attrition execution process took $(mpSim.attrExecTimeElapsed.value / 1000) seconds." )

    return
    # Wipe the simulation reports if the simulation time has advanced.
    # if oldSimTime < now( mpSim )
    #     empty!( mpSim.simReports )
    # end  # if oldSimTime < now( mpSim )

end  # run( mpSim, toTime )

# This function runs the simulation from file.
export runSimFromFile
function runSimFromFile( fName::String )

    if !ispath( fName )
        warn( "File does not exist." )
        return
    elseif !endswith( fName, ".xlsx" )
        warn( "File is not a .xlsx file" )
        return
    end  # if !ispath( fName )

    # Initialise simulation.
    mpSim = ManpowerSimulation( fName )
    println( "Simulation initialised." )

    isRunRequested = true

    XLSX.openxlsx( fName ) do xf
        # Make network plot if requested.
        if XLSX.hassheet( xf, "State Map" )
            sheet = xf[ "State Map" ]

            if sheet[ "B3" ] == "YES"
                println( "Creating network plot." )
                tStart = now()
                plotTransitionMap( mpSim, sheet )
                tEnd = now()
                timeElapsed = (tEnd - tStart).value / 1000
                println( "Network plot time: $timeElapsed seconds." )
            end  # if sheet[ "B3" ] == "YES"
        end  # if XLSX.hassheet( xf, "State Map" )

        # Run only if flag is okay.
        sheet = xf[ "General" ]
        isRunRequested = sheet[ "B14" ] == "YES"
    end  # XLSX.openxlsx( fName ) do xf

    if !isRunRequested
        println( "No simulation run requested." )
        return
    end  # if !isRunRequested

    println( "Running simulation." )
    tStart = now()
    run( mpSim )
    tEnd = now()
    timeElapsed = (tEnd - tStart).value / 1000
    println( "Simulation time: $timeElapsed seconds." )

    # Generate plots.
    println( "Generating plots. This can take a while..." )
    tStart = now()
    showPlotsFromFile( mpSim, fName )
    tEnd = now()
    timeElapsed = (tEnd - tStart).value / 1000
    println( "Plot generation time: $timeElapsed seconds." )

end  # runSimFromFile( fName )

# These are the functions that process simulation results.
include( joinpath( funcPath, "simProcessing.jl" ) )

function Base.show( io::IO, mpSim::ManpowerSimulation )

    if mpSim.dbName == ""
        print( io, "Simulation database kept in memory." )
    else
        print( io, "Simulation database file: \"$(mpSim.dbName)\"" )
    end  # if mpSim.dbName == ""

    print( io, "\nSimulation name: \"$(mpSim.simName)\"" )
    print( io, "\nInitialization state: " )
    print( io, isInitialised( mpSim ) ? "OK" : "not properly initialised" )

    if mpSim.personnelTarget > 0
        print( io, "\nPersonnel cap: $(mpSim.personnelTarget)" )
    end  # if mpSim.personnelTarget > 0

    if isempty( mpSim.initAttrList )
        print( io, "\nNo initialised attributes" )
    else
        print( io, "\nInitialised attributes" )
        foreach( attr -> print( io, "\n$attr" ), mpSim.initAttrList )
    end  # if isempty( mpSim.initAttrList )

    if isempty( mpSim.otherAttrList )
        print( io, "\nNo other attributes" )
    else
        print( io, "\nOther attributes" )
        foreach( attr -> print( io, "\n$attr" ), mpSim.otherAttrList )
    end  # if isempty( mpSim.otherAttrList )

    if isempty( mpSim.initStateList )
        print( io, "\nNo initial personnel states" )
    else
        print( io, "\nInitial personnel states & transitions" )

        for state in keys( mpSim.initStateList )
            print( io, "\n$state" )

            if isempty( mpSim.initStateList[ state ] )
                print( io, "\n    No transitions from this state" )
            else
                foreach( trans -> print( io, "\n$trans" ),
                    mpSim.initStateList[ state ] )
            end  # if isempty( mpSim.initStateList[ state ] )
        end  # for state in  keys( mpSim.initStateList )
    end  # if isempty( mpSim.initStateList )

    if isempty( mpSim.otherStateList )
        print( io, "\nNo other personnel states" )
    else
        print( io, "\nOther personnel states & transitions" )

        for state in keys( mpSim.otherStateList )
            print( io, "\n$state" )

            if isempty( mpSim.otherStateList[ state ] )
                print( io, "\n    No transitions from this state" )
            else
                foreach( trans -> print( io, "\n$trans" ),
                    mpSim.otherStateList[ state ] )
            end  # if isempty( mpSim.otherStateList[ state ] )
        end  # for state in  keys( mpSim.otherStateList )
    end  # if isempty( mpSim.otherStateList )

    if !isempty( mpSim.recruitmentSchemes )
        print( io, "\nRecruitment schemes" )
        foreach( recScheme -> print( io, "\n$recScheme" ),
            mpSim.recruitmentSchemes )
    end  # if !isempty( mpSim.recruitmentSchemes )

    if isa( mpSim.defaultAttritionScheme, Attrition )
        print( io, "\nDefault attrition scheme: $(mpSim.defaultAttritionScheme)" )
    end  #if isa( mpSim.defaultAttritionScheme, Attrition )

    if isa( mpSim.retirementScheme, Retirement )
        print( io, "\nRetirement scheme: $(mpSim.retirementScheme)" )
    end  # if isa( mpSim.retirementScheme, Retirement )

    print( io, "\nSimulation length: $(mpSim.simLength)" )

    if isInitialised( mpSim )
        print( io, "\nCurrent simulation time: $(now( mpSim ))" )

        if now( mpSim ) > 0.0
            print( io, "\nSimulation run time: $(mpSim.simTimeElapsed.value / 1000) seconds" )
        end  # if now( mpSim ) > 0.0
    end  # if isInitialised( mpSim )

    print( io, "\nID key: $(string( mpSim.idKey ))" )

    print( io, "\nSimulation phase priorities: " )
    print( io, join( map( phase -> "$phase: $(mpSim.phasePriorities[ phase ])",
        keys( mpSim.phasePriorities ) ), "; " ) )

    if !isInitialised( mpSim )
        return
    end  # if !isInitialised( mpSim )

    print( io, "\nCurrent number of active personnel in database: ",
        "$(mpSim.personnelSize)" )
    print( io, "\nTotal number of personnel in simulation: ",
        "$(mpSim.resultSize)" )

end  # show( io, mpSim )


text:
# This definitely needs cleaning
using OrdinaryDiffEq, DiffEqDevTools, Test, Random
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear
probArr = Vector{ODEProblem}(undef, 2)
probArr[1] = prob_ode_linear

probArr[2] = prob_ode_2Dlinear
Random.seed!(100)
## Convergence Testing
dts = 1 .//2 .^(8:-1:4)
dts1 = 1 .//2 .^(9:-1:5)
testTol = 0.2

for i = 1:2
  global dts
  prob = probArr[i]
  sim = test_convergence(dts,prob,Euler())
  @test sim.𝒪est[:final] ≈ 1 atol=testTol
  sim2 = test_convergence(dts,prob,Heun())
  @test sim2.𝒪est[:l∞] ≈ 2 atol=testTol
  sim2 = test_convergence(dts,prob,Ralston())
  @test sim2.𝒪est[:l∞] ≈ 2 atol=testTol
  sim2 = test_convergence(dts,prob,Midpoint())
  @test sim2.𝒪est[:l∞] ≈ 2 atol=testTol
  sim3 = test_convergence(dts,prob,RK4())
  @test sim3.𝒪est[:l∞] ≈ 4 atol=testTol
  sim4 = test_convergence(dts,prob,BS3())
  @test sim4.𝒪est[:l2] ≈ 3 atol=testTol
  alg = CarpenterKennedy2N54()
  @test test_convergence(dts,prob,alg).𝒪est[:l∞] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim5 = test_convergence(dts, prob, AB3())
  @test sim5.𝒪est[:l2] ≈ 3 atol=testTol
  sim6 = test_convergence(dts,prob,ABM32())
  @test sim6.𝒪est[:l2] ≈ 3 atol=testTol
  sim7 = test_convergence(dts, prob, AB4())
  @test sim7.𝒪est[:l2] ≈ 4 atol=testTol
  sim8 = test_convergence(dts1,prob,ABM43())  #using dts1 due to floating point error in convergence test
  @test sim8.𝒪est[:l2] ≈ 4 atol=testTol
  sim9 = test_convergence(dts,prob,AB5())
  @test sim9.𝒪est[:l2] ≈ 5 atol=testTol
  sim10 = test_convergence(dts,prob,ABM54())
  @test sim10.𝒪est[:l2] ≈ 5 atol=testTol
  sim101 = test_convergence(dts,prob,VCAB3())
  @test sim101.𝒪est[:l2] ≈ 3 atol=testTol
  sim102 = test_convergence(dts,prob,VCAB4())
  @test sim102.𝒪est[:l2] ≈ 4 atol=testTol
  sim103 = test_convergence(dts,prob,VCAB5())
  @test sim103.𝒪est[:l2] ≈ 5 atol=testTol
  sim104 = test_convergence(dts,prob,VCABM3())
  @test sim104.𝒪est[:l2] ≈ 3 atol=testTol
  sim105 = test_convergence(dts,prob,VCABM4())
  @test sim105.𝒪est[:l2] ≈ 4 atol=testTol
  sim106 = test_convergence(dts,prob,VCABM5())
  @test sim106.𝒪est[:l2] ≈ 5 atol=testTol
  sim160 = test_convergence(dts,prob,Anas5(w=2))
  @test sim160.𝒪est[:l2] ≈ 4 atol=2*testTol

  println("Stiff Solvers")

  dts = 1 .//2 .^(9:-1:5)

  sim11 = test_convergence(dts,prob,ImplicitEuler(extrapolant = :linear))
  @test sim11.𝒪est[:final] ≈ 1 atol=testTol

  sim112 = test_convergence(dts,prob,ImplicitEuler(nlsolve = NLFunctional()))
  @test sim112.𝒪est[:final] ≈ 1 atol=testTol

  sim12 = test_convergence(dts,prob,
          GenericImplicitEuler(nlsolve=OrdinaryDiffEq.NLSOLVEJL_SETUP(autodiff=true)))
  @test sim12.𝒪est[:final] ≈ 1 atol=testTol

  sim122 = test_convergence(dts,prob,
           GenericImplicitEuler(nlsolve=OrdinaryDiffEq.NLSOLVEJL_SETUP(autodiff=false)))

  sim13 = test_convergence(dts,prob,ImplicitMidpoint())
  @test sim13.𝒪est[:final] ≈ 2 atol=testTol

  sim132 = test_convergence(dts,prob,ImplicitMidpoint(nlsolve = NLFunctional()))
  @test sim132.𝒪est[:final] ≈ 2 atol=testTol

  sim13 = test_convergence(dts,prob,Trapezoid())
  @test sim13.𝒪est[:final] ≈ 2 atol=testTol

  sim133 = test_convergence(dts,prob,Trapezoid(nlsolve = NLFunctional()))
  @test sim133.𝒪est[:final] ≈ 2 atol=testTol

  sim14 = test_convergence(dts,prob,
          GenericTrapezoid(nlsolve=OrdinaryDiffEq.NLSOLVEJL_SETUP(autodiff=true)))
  @test sim14.𝒪est[:final] ≈ 2 atol=testTol

  sim142 = test_convergence(dts,prob,
           GenericTrapezoid(nlsolve=OrdinaryDiffEq.NLSOLVEJL_SETUP(autodiff=false)))
  @test sim142.𝒪est[:final] ≈ 2 atol=testTol

  sim14 = test_convergence(dts,prob,TRBDF2())
  @test sim14.𝒪est[:final] ≈ 2 atol=testTol

  sim152 = test_convergence(dts,prob,TRBDF2(autodiff=false))
  @test sim152.𝒪est[:final] ≈ 2 atol=testTol+0.1

  sim15 = test_convergence(dts,prob,SDIRK2())
  @test sim15.𝒪est[:final] ≈ 2 atol=testTol

  sim152 = test_convergence(dts,prob,SSPSDIRK2())
  @test sim152.𝒪est[:final] ≈ 2 atol=testTol

  sim16 = test_convergence(dts,prob,Kvaerno3())
  @test sim16.𝒪est[:final] ≈ 3 atol=testTol

  sim162 = test_convergence(dts,prob,Kvaerno3(nlsolve = NLFunctional()))
  @test sim162.𝒪est[:final] ≈ 3 atol=testTol

  sim17 = test_convergence(dts,prob,KenCarp3())
  @test sim17.𝒪est[:final] ≈ 3 atol=testTol

  #####################################
  # BDF
  #####################################

  sim = test_convergence(dts,prob,ABDF2(extrapolant = :linear))
  @test sim.𝒪est[:final] ≈ 2 atol=testTol
  @test sim.𝒪est[:l2] ≈ 2 atol=testTol
  @test sim.𝒪est[:l∞] ≈ 2 atol=testTol

  sim = test_convergence(dts,prob,ABDF2(nlsolve = NLFunctional()))
  @test sim.𝒪est[:final] ≈ 2 atol=testTol
  @test sim.𝒪est[:l2] ≈ 2 atol=testTol
  @test sim.𝒪est[:l∞] ≈ 2 atol=testTol

  # QBDF
  sim = test_convergence(dts,prob,QBDF1())
  @test sim.𝒪est[:final] ≈ 1 atol=testTol
  @test sim.𝒪est[:l2] ≈ 1 atol=testTol
  @test sim.𝒪est[:l∞] ≈ 1 atol=testTol

  sim = test_convergence(dts,prob,QBDF2())
  @test sim.𝒪est[:final] ≈ 2 atol=testTol
  @test sim.𝒪est[:l2] ≈ 2 atol=testTol
  @test sim.𝒪est[:l∞] ≈ 2 atol=testTol

  # QNDF
  sim = test_convergence(dts,prob,QNDF1())
  @test sim.𝒪est[:final] ≈ 1 atol=testTol
  @test sim.𝒪est[:l2] ≈ 1 atol=testTol
  @test sim.𝒪est[:l∞] ≈ 1 atol=testTol

  sim = test_convergence(dts,prob,QNDF2())
  @test sim.𝒪est[:final] ≈ 2 atol=testTol
  @test sim.𝒪est[:l2] ≈ 2 atol=testTol
  @test sim.𝒪est[:l∞] ≈ 2 atol=testTol

  sim = test_convergence(dts,prob,QNDF2(nlsolve = NLFunctional()))
  @test sim.𝒪est[:final] ≈ 2 atol=testTol
  @test sim.𝒪est[:l2] ≈ 2 atol=testTol
  @test sim.𝒪est[:l∞] ≈ 2 atol=testTol

  dts = 1 .//2 .^(7:-1:4)
  println("Higher Order")

  sim18 = test_convergence(dts,prob,Cash4())
  @test sim18.𝒪est[:final] ≈ 4 atol=testTol

  sim19 = test_convergence(dts,prob,Hairer4())
  @test sim19.𝒪est[:final] ≈ 4 atol=testTol

  sim110 = test_convergence(dts,prob,Hairer42())
  @test sim110.𝒪est[:final] ≈ 4 atol=testTol

  sim111 = test_convergence(dts,prob,Kvaerno4())
  @test sim111.𝒪est[:final] ≈ 4 atol=testTol

  sim112 = test_convergence(dts,prob,KenCarp4())
  @test sim112.𝒪est[:final] ≈ 4 atol=testTol

  sim113 = test_convergence(dts,prob,Kvaerno5())
  @test sim113.𝒪est[:final] ≈ 5 atol=testTol

  sim114 = test_convergence(dts,prob,KenCarp5())
  @test sim114.𝒪est[:final] ≈ 5 atol=testTol

  sim115 = test_convergence(dts,prob,KenCarp5(nlsolve = NLFunctional()))
  @test_broken sim115.𝒪est[:final] ≈ 5 atol=testTol
end


text:
import ProximalBase
import CoordinateDescent

reload("ProximalBase")
reload("CoordinateDescent")

n = 3000
p = 5000
s = 100

srand(123)
X = randn(n, p)
Y = X[:,1:s] * (randn(s) .* (1. .+ rand(s))) + 6. * randn(n)

stdX = std(X, 1)[:]

options = CoordinateDescent.ScaledLassoOptions(;optTol=1e-3, maxIter=50)
x = ProximalBase.SparseIterate(p)
λ = sqrt(2. * log(p) / n)
@time CoordinateDescent.scaledLasso!(x, X, Y, λ, stdX, options)

@show σinit = CoordinateDescent.findInitSigma(X, Y, 30)
options = CoordinateDescent.ScaledLassoOptions(;optTol=1e-2, maxIter=10, σinit=σinit)
x = ProximalBase.SparseIterate(p)
λ = sqrt(2. * log(p) / n)
@time CoordinateDescent.scaledLasso!(x, X, Y, λ, stdX, options)

λ = 0.001

options = CoordinateDescent.CDOptions(;warmStart=true)
x = ProximalBase.SparseIterate(p)
f = CoordinateDescent.CDLeastSquaresLoss(Y,X)
g = ProximalBase.ProxL1(λ)
@time CoordinateDescent.coordinateDescent!(x, f, g, options)

options = CoordinateDescent.CDOptions(;warmStart=false)
x = ProximalBase.SparseIterate(p)
f = CoordinateDescent.CDLeastSquaresLoss(Y,X)
g = ProximalBase.ProxL1(λ)
@time CoordinateDescent.coordinateDescent!(x, f, g, options)

options = CoordinateDescent.CDOptions(;warmStart=true)
x = ProximalBase.SparseIterate(p)
f = CoordinateDescent.CDLeastSquaresLoss(Y,X)
g = ProximalBase.ProxL1(λ, stdX)
@time CoordinateDescent.coordinateDescent!(x, f, g, options)

options = CoordinateDescent.CDOptions(;warmStart=false, numSteps=100)
x = ProximalBase.SparseIterate(p)
f = CoordinateDescent.CDLeastSquaresLoss(Y,X)
g = ProximalBase.ProxL1(λ, stdX)
@time CoordinateDescent.coordinateDescent!(x, f, g, options)

options = CoordinateDescent.CDOptions(;warmStart=false, numSteps=100)
x = ProximalBase.SparseIterate(p)
f = CoordinateDescent.CDLeastSquaresLoss(Y,X)
g = ProximalBase.ProxL1(λ, stdX)
@time CoordinateDescent.coordinateDescent!(x, f, g, options)


text:
using StatsBase, Distributions

@testset "moments" begin
    @testset "two moments" begin
        m = @inferred Moments(1,0.5)
        @test n_moments(typeof(m)) == 2
        @test n_moments(m) == 2
        #@test n_moments(typeof(M)) == 2
        @test (@inferred mean(m)) == m[1] == 1
        @test (@inferred var(m)) == m[2] == 0.5
        @test (@inferred std(m)) == sqrt(0.5)
        @test_throws Exception skewness(m)
        @test_throws Exception kurtosis(m)
        @test (@inferred convert(AbstractArray, m)) == [1.0, 0.5]
    end;
    @testset "zero moments" begin
        m = @inferred  Moments()
        @test (@inferred n_moments(m)) == 0
        @test_throws Exception mean(m)
        @test_throws Exception m[1]
    end;
end;

@testset "moments of Distribution" begin
    function testMoments(d)
        m = @inferred moments(d, Val(4))
        @test [m[i] for i in 1:n_moments(m)] == [mean(d), var(d), skewness(d), kurtosis(d)]
        m = @inferred moments(d, Val(3))
        @test [m[i] for i in 1:n_moments(m)] == [mean(d), var(d), skewness(d)]
        m = @inferred moments(d, Val(2))
        @test [m[i] for i in 1:n_moments(m)] == [mean(d), var(d)]
        m = @inferred moments(d, Val(1))
        @test [m[i] for i in 1:n_moments(m)] == [mean(d)]
        m = @inferred moments(d, Val(0))
        @test [m[i] for i in 1:(@inferred n_moments(m))] == []
        @test_throws Exception moments(d, Val(-1))
        @test_throws Exception moments(d, Val(:bla))
    end
    @testset "Normal" begin
        d = @inferred Normal(2,5)
        testMoments(d)
    end;
    @testset "Lognormal" begin
        testMoments((@inferred LogNormal(2,5)))
    end;
    @testset "Logitnormal" begin
        # mean not defined fro LogitNormal
        @test_throws MethodError testMoments((@inferred LogitNormal(2,5)))
    end;
end;

@testset "QuantilePoint" begin
    qp1 = @inferred QuantilePoint(2,0.25);
    qp2 = @inferred QuantilePoint(7.5,0.95);
    @testset "isless" begin
        #println("isless")
        # p and q both increasing/decreasing
        @test (@inferred(QuantilePoint(2,0.4) < QuantilePoint(3,0.5))) == true
        @test (@inferred(QuantilePoint(1,0.6) < QuantilePoint(1,0.5))) == false
        # same q may have different p as same samples can be repeated
        @test (@inferred(QuantilePoint(2,0.4) < QuantilePoint(2,0.5))) == false
        @test (@inferred(QuantilePoint(2,0.6) < QuantilePoint(2,0.5))) == false
        # q decreasing but p increasing or vise-versa
        @test_throws ErrorException QuantilePoint(1,0.6) < QuantilePoint(2,0.5)
        @test_throws ErrorException QuantilePoint(2,0.4) < QuantilePoint(1,0.5)
        # for same p equal, also q must be equal
        @test_throws ErrorException  QuantilePoint(1,0.5) < QuantilePoint(2,0.5)
        @test_throws ErrorException  QuantilePoint(2,0.5) < QuantilePoint(1,0.5)
    end;
    @testset "modifying Quantilepoint" begin
        qp3 = QuantilePoint(qp1, q = 3);
        @test qp3.q == 3
        @test qp3.p == qp1.p == 0.25
    end;
    @testset "macros" begin
        @test @qp(0.4,0.7) == QuantilePoint(0.4,0.7)
        @test @qp_ll(0.7) == QuantilePoint(0.7,0.025)
        @test @qp_l(0.7) == QuantilePoint(0.7,0.05)
        @test @qp_m(0.7) == QuantilePoint(0.7,0.5)
        @test @qp_u(0.7) == QuantilePoint(0.7,0.95)
        @test @qp_uu(0.7) == QuantilePoint(0.7,0.975)
        @test @qs_cf90(0.2,0.7) == Set([@qp_l(0.2),@qp_u(0.7)])
        @test @qs_cf95(0.2,0.7) == Set([@qp_ll(0.2),@qp_uu(0.7)])
    end;
end;

@testset "global fit functions" begin
    @test_throws ErrorException fit(Distribution, Moments())
    qpl = @qp_m(3)
    qpu = @qp_u(5)
    @test_throws ErrorException fit(Distribution, qpl, qpu)
    @testset "fit to stats and quantilepoint of unknown distribution" begin
        @test_throws MethodError fit(Distribution, 3.0, qpu)
        @test_throws MethodError fit(Distribution, 3.0, qpu, Val(:mean))
        @test_throws MethodError fit(Distribution, 3.0, qpu, Val(:mode))
        @test_throws ErrorException fit(Distribution, 3.0, qpu, Val(:median))
        @test_throws ErrorException fit(Distribution, 3.0, qpu, Val(:bla))
    end;
end;


text:
using TestTUI
using Test

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


text:
struct SubSpace{DS,IT,DD,RR} <: Space{DD,RR}
    space::DS
    indexes::IT
    SubSpace{DS,IT,DD,RR}(sp::DS,ind::IT) where {DS,IT,DD,RR} = new(sp,ind)
end

SubSpace(sp::Space,kr) =
    SubSpace{typeof(sp),typeof(kr),domaintype(sp),rangetype(sp)}(sp,kr)

SubSpace(sp::SubSpace,kr) = SubSpace(sp.space,reindex(sp,sp.indexes,to_indexes(kr))[1])

domain(DS::SubSpace) = domain(DS.space)
dimension(sp::SubSpace) = length(sp.indexes)

|(sp::Space,kr::AbstractRange) = SubSpace(sp,kr)

function |(f::Fun,kr::AbstractInfUnitRange)
    @assert dimension(space(f)) == ∞
    Fun(space(f)|kr,f.coefficients[kr[1]:end])
end

block(sp::SubSpace, k::Integer) = block(sp.space,reindex(sp,(sp.indexes,),(k,))[1])

function blocklengths(sp::SubSpace{DS,UnitRange{Int}}) where DS
    N = first(sp.indexes)
    M = last(sp.indexes)
    B1=block(sp.space,N)
    B2=block(sp.space,M)
    # if the blocks are equal, we have only one bvlock
    B1 == B2 && return [zeros(Int,B1.n[1]-1);length(sp.indexes)]

    [zeros(Int,B1.n[1]-1);
         blockstop(sp.space,B1)-N+1;blocklengths(sp.space)[B1.n[1]+1:B2.n[1]-1];
        M-blockstart(sp.space,B2)+1]
end

function blocklengths(sp::SubSpace{DS,<:AbstractInfUnitRange{Int}}) where DS
    N = first(sp.indexes)
    B1=block(sp.space,N)

    Vcat([zeros(Int,B1.n[1]-1); blockstop(sp.space,B1)-N+1],
            blocklengths(sp.space)[B1.n[1]+1:∞])
end

blocklengths(sp::SubSpace{DS,Block{1,T}}) where {DS, T} =
    [blocklengths(sp.space)[Int(sp.indexes)]]

blocklengths(sp::SubSpace) = error("Not implemented for non-unitrange subspaces")

## Block reindexing for SubSpace
reindex(sp::SubSpace, b::Tuple{Block{1}}, ks::Tuple{Any}) = (blockstart(sp.space,b[1]).+ks[1].-1,)
reindex(sp::SubSpace, br::Tuple{BlockRange1}, ks::Tuple{Block{1}}) = (Block(Int.(br[1])[first(ks[1].n)]),)
reindex(sp::SubSpace, br::Tuple{BlockRange1}, ks::Tuple{Any}) = (blockstart(sp.space,first(Int.(br[1]))).+ks[1].-1,)

# blocks stay the same with unit range indices
reindex(sp::SubSpace, br::Tuple{AbstractVector{Int}}, ks::Tuple{Block{1}}) =
    reindex(sp, br, (blockrange(sp,first(ks)),))
reindex(sp::SubSpace, br::Tuple{AbstractVector{Int}}, ks::Tuple{BlockRange1}) =
    reindex(sp, br, (blockrange(sp,first(ks)),))

## Block
blocklengths(sp::SubSpace{DS,Block}) where {DS} = [blocklengths(sp.space)[sp.indexes.n[1]]]
dimension(sp::SubSpace{DS,Block}) where {DS} = blocklengths(sp.space)[sp.indexes.n[1]]

blocklengths(sp::SubSpace{DS,BlockRange1}) where {DS} = blocklengths(sp.space)[Int.(sp.indexes)]
dimension(sp::SubSpace{DS,BlockRange1}) where {DS} = sum(blocklengths(sp.space)[Int.(sp.indexes)])

block(_,B::Block) = B

##

spacescompatible(S1::SubSpace{DS,IT,DD,RR},S2::SubSpace{DS,IT,DD,RR}) where {DS,IT,DD,RR} =
    spacescompatible(S1.space,S2.space) && S1.indexes == S2.indexes

==(S1::SubSpace{DS,IT,DD,RR},S2::SubSpace{DS,IT,DD,RR}) where {DS,IT,DD,RR} =
    S1.space == S2.space && S1.indexes == S2.indexes

canonicalspace(a::SubSpace) = a.space

setdomain(DS::SubSpace,d::Domain) = SubSpace(setdomain(DS.space,d),DS.indexes)

Conversion(a::SubSpace,b::Space) = ConcreteConversion(a,b)

function Conversion(a::S,b::SubSpace{S,<:AbstractInfUnitRange{Int}}) where S<:Space
    @assert first(b.indexes) == 1
    ConversionWrapper(SpaceOperator(Operator(I,a),a,b))
end

bandwidths(C::ConcreteConversion{<:SubSpace{<:Any,<:AbstractInfUnitRange{Int}}}) =
    first(domainspace(C).indexes)-1,1-first(domainspace(C).indexes)

getindex(C::ConcreteConversion{<:SubSpace}, k::Integer,j::Integer) =
    domainspace(C).indexes[j]==k ? one(eltype(C)) : zero(eltype(C))

# avoid ambiguity
for OP in (:leftendpoint,:rightendpoint)
    @eval getindex(E::ConcreteEvaluation{SubSpace{DS,IT,DD,RR},typeof($OP)},k::Integer) where {IT,DS,DD,RR}=
        Evaluation(E.space.space,E.x,E.order)[E.space.indexes[k]]
end
getindex(E::ConcreteEvaluation{SubSpace{DS,IT,DD,RR}},k::Integer) where {IT,DS,DD,RR}=
    Evaluation(E.space.space,E.x,E.order)[E.space.indexes[k]]
getindex(E::ConcreteEvaluation{SubSpace{DS,IT,DD,RR}},kr::AbstractRange) where {IT,DS,DD,RR} =
    Evaluation(E.space.space,E.x,E.order)[E.space.indexes[kr]]

function conversion_rule(a::SubSpace,b::SubSpace)
     if a==b
        a
    else
        NoSpace()
    end
end
# return the space that has banded Conversion to the other
function conversion_rule(a::SubSpace,b::Space)
    if a.space==b
        a  # we can write droping coefficients as a banded operator
    else
        NoSpace()
    end
end

function union_rule(a::SubSpace,b::SubSpace)
     if a == b
        a
    else
        NoSpace()
    end
end
# return the space that has banded Conversion to the other
function union_rule(a::SubSpace,b::Space)
    if a.space==b
        b  # we can write droping coefficients as a banded operator
    else
        NoSpace()
    end
end

function subspace_coefficients(v::AbstractVector,sp::SubSpace,dropsp::SubSpace)
    if sp == dropsp
        v
    else
        coefficients(v,sp,canonicalspace(sp),dropsp)
    end
end

function subspace_coefficients(v::AbstractVector,sp::Space,dropsp::SubSpace)
    n=length(v)
    if sp == dropsp.space
        ret = Array{eltype(v)}(undef, 0)
        for k in dropsp.indexes
            if k > n
                return ret
            end
            push!(ret,v[k])
        end
        ret
    else
        coefficients(v,sp,canonicalspace(dropsp),dropsp)
    end
end

function subspace_coefficients(v::AbstractVector,dropsp::SubSpace,sp::Space)
    if sp==dropsp.space
        isempty(v) && return v
        ret = zeros(eltype(v),dropsp.indexes[length(v)])
        for k = eachindex(v)
            ret[dropsp.indexes[k]] = v[k]
        end
        ret
    else
        coefficients(v,dropsp,canonicalspace(dropsp),sp)
    end
end

coefficients(v::AbstractVector,sp::SubSpace,dropsp::SubSpace) = subspace_coefficients(v,sp,dropsp)

## points

## transform
function transform(sp::SubSpace, vals::AbstractVector)
    ret = transform(sp.space,vals)
    coefficients(ret,sp.space,sp)
end

itransform(sp::SubSpace,cfs::AbstractVector) =
    itransform(sp.space,coefficients(cfs,sp,sp.space))

for Tran in (:plan_transform, :plan_transform!, :plan_itransform, :plan_itransform!)
    @eval $Tran(sp::SubSpace) = $Tran(sp, dimension(sp))
end

points(sp::SubSpace, n) = points(sp.space, n)
points(sp::SubSpace) = points(sp, dimension(sp))

coefficients(v::AbstractVector,::SubSpace{DS,IT,Segment{Vec{2,TT}}},::TensorSpace{SV,DD}) where {DS,IT,TT,SV,DD<:Domain2d} =
    error("Not callable, only defined for ambiguity errors.")
coefficients(v::AbstractVector,::SubSpace{DS,IT,D},::TensorSpace{SV,DD}) where {DS,IT,D,SV,DD<:Domain2d} =
    error("Not callable, only defined for ambiguity errors.")

for TYP in (:SumSpace,:PiecewiseSpace,:TensorSpace,:ConstantSpace,:Space) # Resolve conflict
    @eval begin
        coefficients(v::AbstractVector,sp::$TYP,dropsp::SubSpace) = subspace_coefficients(v,sp,dropsp)
        coefficients(v::AbstractVector,dropsp::SubSpace,sp::$TYP) = subspace_coefficients(v,dropsp,sp)
    end
end

## ProductFUn

# values{S<:SubSpace,V<:SubSpace}(f::ProductFun{S,V})=values(ProductFun(f,space(f,1).space,space(f,2).space))
# values{S<:SubSpace}(f::ProductFun{S})=values(ProductFun(f,space(f,1).space,space(f,2)))
#
#
# function coefficients{n,DS,TT,DD}(f::ProductFun{SubSpace{n,1,DS,DD,RR}},ox::Space,oy::Space)
#     T=eltype(f)
#     m=size(f,1)
#     A=[pad!(coefficients(fx,ox),m+n) for fx in f.coefficients]
#     B=hcat(A...)::Array{T,2}
#     for k=1:size(B,1)
#         ccfs=coefficients(B[k,:],space(f,2),oy)
#         if length(ccfs)>size(B,2)
#             B=pad(B,size(B,1),length(ccfs))
#         end
#         B[k,1:length(ccfs)]=ccfs
#         #B[k,length(ccfs):1:end]=zero(T)
#     end
#
#     B
# end


text:
using CSV
using DelimitedFiles

function loadtestdata(testCase::String, dataPath::String=@__DIR__)
    #println(pwd()) # to debug github action
    data = readdlm(joinpath(dataPath, "data/$(testCase)_data.csv"), ',', Float64, '\n')
    data = dropdims(data,dims=1) # convert to vector
    evts = CSV.read(joinpath(dataPath, "data/$(testCase)_events.csv"), DataFrame)
    return data, evts
end


text:
using LinearAlgebra, Random, Test
using SafeTestsets

@time @safetestset "metrics" begin include("metrics.jl") end


text:
module ModuleInfo

using CodeTracking: pkgfiles, whereis
using StructArrays: StructArray
using InlineTest
using Markdown

const Maybe{T} = Union{Nothing, T}

function moduleinfo(m::Module)
    return moduleinfo!(createinfostore(), m)
end

function moduleinfo(ms::AbstractVector{Module})
    info = createinfostore()
    foreach(m -> moduleinfo!(info, m), ms)
    return info
end

export moduleinfo

"""
    createinfostore()

Create a new database to hold module information.
"""
function createinfostore()
    return Dict(
        :packages => StructArray(
            package_id = String[],
            name = String[],
            basedir = String[],
        ),
        :modules => StructArray(
            module_id = String[],
            parent = Maybe{String}[],
            package_id = Maybe{String}[],  # foreign key
            instance = Module[],
        ),
        :sourcefiles => StructArray(
            package_id = String[],  # foreign key
            file = String[],
        ),
        :symbols => StructArray(
            symbol_id = String[],
            name = String[],
            public = Bool[],
            kind = String[],
            instance = Any[],
            module_id = String[],  # foreign key
        ),
        :methods => StructArray(
            method_id = String[],
            symbol_id = String[],  # foreign key
            line = Int[],
            file = String[],
            signature = Any[],
        ),
        :docstrings => StructArray(
            symbol_id = String[],  # foreign key
            docstring = String[],
        )
    )
end

function moduleinfo!(db, m::Module)
    info = pkgfiles(m)
    pkgid = string(info.id.uuid)

    # If `m` is a top-level module
    if parentmodulerec(m) == m
        # Add package entry
        push!(db[:packages], (
            package_id = pkgid,
            name = info.id.name,
            basedir = info.basedir
        ))

    # Add source files of package
        for file in info.files
            if isfile(joinpath(pkgdir(m), file))
                push!(db[:sourcefiles], (package_id = pkgid, file = file))
            end
        end
    end

    # Add module
    modulename = join(fullname(m), ".")
    push!(db[:modules], (
        module_id = modulename,
        parent = parentmodule(m) === m ? nothing : join(fullname(parentmodule(m)), "."),
        package_id = pkgid,
        instance = m,
    ))

    gathermodulesymbols!(db, m)

    return db
end

getmoduleid(m::Module) = join(fullname(m), ".")
getsymbolid(T) = join((getmoduleid(parentmodule(T)), string(nameof(T))), ".")

function gathermodulesymbols!(db, m::Module)
    moduleid = getmoduleid(m)
    childmodules = Module[]

    for subm in submodules(m)
        moduleinfo!(db, subm)
    end

    for symbol in names(m, all=true)
        if (!isdefined(m, symbol) || symbol === :eval ||
                symbol == :include || !isfrommodule(m, symbol) || startswith(string(symbol), "#"))
            continue
        end

        instance = getfield(m, symbol)
        kind = symbolkind(m, symbol)
        symbol_id = m === instance ? moduleid : "$moduleid.$symbol"

        push!(db[:symbols], (;
            symbol_id = symbol_id,
            name = string(symbol),
            public = Base.isexported(m, symbol),
            kind = kind,
            instance,
            module_id = moduleid,
        ))

        gathermethods!(db, m, symbol)
        gatherdocstring!(db, m, symbol)
    end

    for m in childmodules
        moduleinfo!(db, m)
    end
end

getsymbolid(m, symbol) = join(fullname(m), ".") * "." * string(symbol)

function gatherdocstring!(db, m, symbol)
    # TODO: store method metadata
    docstrs, _ = getdocstrings(m, symbol)
    isempty(docstrs) && return
    docstring = join(docstrs, "\n\n---\n\n")

    instance = getfield(m, symbol)
    symbol_id = m === instance ? getmoduleid(m) : getsymbolid(m, symbol)

    push!(db[:docstrings], (; symbol_id, docstring))
end

function parentmodulerec(x)
    try
        parentmodule(x) == x ? x : parentmodulerec(parentmodule(x))
    catch e
        e isa MethodError && return true
        rethrow()
    end
end

parentmodulerec(e::E) where {E<:Enum} = parentmodulerec(E)

function submodules(m::Module)
    return map(s -> getfield(m, s), filter(names(m, all=true)) do s
        isdefined(m, s) || return false
        subm = getfield(m, s)
        return subm isa Module && subm !== m && parentmodule(subm) === m
    end)
end

isfrommodule(m::Module, symbol) = isfrommodule(m, symbol, submodules(m))
function isfrommodule(m::Module, symbol, submodules)
    if any(isfrommodule(sm, symbol) for sm in submodules)
        return false
    end
    isdefined(m, symbol) || return false
    x = getfield(m, symbol)
    try
        return parentmodule(x) === m
    catch
        isconst(m, symbol) && return true
        return false
    end
end

function gathermethods!(db, m::Module, symbol::Symbol)
    modulename = join(fullname(m), ".")
    symbol_id = "$modulename.$symbol"

    x = getfield(m, symbol)
    for (i, method) in enumerate(methods(x))
        file, line = whereis(method)
        sig = join(split(string(method), ")")[1:end-1]) * ")"
        push!(db[:methods], (
            method_id = "$(symbol_id)_$i",
            symbol_id = symbol_id,
            line = line,
            file = file,
            signature = sig
        ))
    end
    return db
end

function symbolkind(m::Module, symbol::Symbol)
    x = getfield(m, symbol)
    x isa DataType && return (isconcretetype(x) ? "struct" : "abstract type")
    x isa UnionAll && return "struct"
    x isa Function && return "function"
    x isa Module && return "module"
    isconst(m, symbol) && return "const"
    return "unknown"
end

struct DocString
    m::Module
    symbol::Symbol
    docstring::String
    data
end

function getdocstrings(m, sym::Symbol)
    multidoc = getmultidoc(m, sym)
    isnothing(multidoc) && return String[], []
    docstrings = [__plain_text(multidoc.docs[sig]) for sig in multidoc.order]
    metadata = [multidoc.docs[sig].data for sig in multidoc.order]
    return docstrings, metadata
end

function getmultidoc(m::Module, symbol::Symbol)
    binding = Base.Docs.Binding(m, symbol)
    try
        docs = getproperty(m, Docs.META)
        return get(docs, binding, nothing)
    catch e
        e isa UndefVarError && return nothing
        rethrow()
    end
end

function __plain_text(d::Base.Docs.DocStr)
    if d.object isa Markdown.MD
        return Markdown.plain(d.object)
    else
        buf = IOBuffer()
        for part in d.text
            print(buf, part)
        end
        return String(take!(buf))
    end
end

module A
    h(z) = 4

    module B
        f(x::Float16) = 1
        g(y) = 3
        export f
    end

    f(a::Int) = 2
    using .B

    export f, B
end
using .A

@testset "isfrommodule" begin
    @test !isfrommodule(A, :f)
    @test isfrommodule(B, :f)
    @test isfrommodule(B, :g)
    @test !isfrommodule(A, :g)
    @test isfrommodule(A, :h)
    @test !isfrommodule(B, :h)
end

end  # module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayBackendHealthPool <: SwaggerModel
    backendAddressPool::Any # spec type: Union{ Nothing, ApplicationGatewayBackendAddressPool } # spec name: backendAddressPool
    backendHttpSettingsCollection::Any # spec type: Union{ Nothing, Vector{ApplicationGatewayBackendHealthHttpSettings} } # spec name: backendHttpSettingsCollection

    function ApplicationGatewayBackendHealthPool(;backendAddressPool=nothing, backendHttpSettingsCollection=nothing)
        o = new()
        validate_property(ApplicationGatewayBackendHealthPool, Symbol("backendAddressPool"), backendAddressPool)
        setfield!(o, Symbol("backendAddressPool"), backendAddressPool)
        validate_property(ApplicationGatewayBackendHealthPool, Symbol("backendHttpSettingsCollection"), backendHttpSettingsCollection)
        setfield!(o, Symbol("backendHttpSettingsCollection"), backendHttpSettingsCollection)
        o
    end
end # type ApplicationGatewayBackendHealthPool

const _property_map_ApplicationGatewayBackendHealthPool = Dict{Symbol,Symbol}(Symbol("backendAddressPool")=>Symbol("backendAddressPool"), Symbol("backendHttpSettingsCollection")=>Symbol("backendHttpSettingsCollection"))
const _property_types_ApplicationGatewayBackendHealthPool = Dict{Symbol,String}(Symbol("backendAddressPool")=>"ApplicationGatewayBackendAddressPool", Symbol("backendHttpSettingsCollection")=>"Vector{ApplicationGatewayBackendHealthHttpSettings}")
Base.propertynames(::Type{ ApplicationGatewayBackendHealthPool }) = collect(keys(_property_map_ApplicationGatewayBackendHealthPool))
Swagger.property_type(::Type{ ApplicationGatewayBackendHealthPool }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayBackendHealthPool[name]))}
Swagger.field_name(::Type{ ApplicationGatewayBackendHealthPool }, property_name::Symbol) =  _property_map_ApplicationGatewayBackendHealthPool[property_name]

function check_required(o::ApplicationGatewayBackendHealthPool)
    true
end

function validate_property(::Type{ ApplicationGatewayBackendHealthPool }, name::Symbol, val)
end


text:
import Base.show
export show

function show{F}(io::IO, dr::DimRedux{F})
    k,n = size(dr.Ξ)
    print(io, "$(typeof(dr)): dimension reduction map over the field $F from $n to $k dimensions")
end

function show{F,DR}(io::IO, sk::Sketch{F,DR})
    k,n = size(sk.X)
    m,k = size(sk.Y)
    s,s = size(sk.Z)
    print(io, "$(typeof(sk)): sketch for an $m×$n matrix over the field $F using a $DR dimension reduction map")
end


text:
module GrowthMaps

# Use the README as the module docs
@doc read(joinpath(dirname(@__DIR__), "README.md"), String) GrowthMaps

using ConstructionBase,
      Dates,
      Rasters,
      InteractModels,
      LsqFit,
      Plots,
      Reexport,
      Unitful,
      UnitfulRecipes

@reexport using ModelParameters

using Rasters.LookupArrays
using Unitful: °C, K, R, Units
using Base: tail

export mapgrowth, fit, fit!, manualfit!, mapfit!

export RateModel

export GrowthModel, SchoolfieldIntrinsicGrowth

export StressModel, LowerStress, UpperStress

export Layer

include("models.jl")
include("framework.jl")
include("fit.jl")

end # module


text:
using DistributedArrays

function dualcdForADMM2(
x::DistributedArrays.DArray{Float64,2,Array{Float64,2}},
y::DistributedArrays.DArray{Float64,1,Array{Float64,1}},
v::Array{Float64,1},
rho,
C,
shrinking,
iterations,
eps)
    dualcdForADMM(x,y,v,rho=rho,C=C,shrinking=shrinking,iterations=iterations,eps=eps)
end

function dualcdForADMM(
x::DistributedArrays.DArray{Float64,2,Array{Float64,2}},
y::DistributedArrays.DArray{Float64,1,Array{Float64,1}},
v::Array{Float64,1};
rho=1.0,
C=2.0,
shrinking=false,
iterations=1000,
eps=0.001)
    x = localpart(x)
    y = localpart(y)

    s,d = size(x) # s: number of data points
    alpha = zeros(s)

    # compute initial w:
    # TODO: 1/rho missing
    w = v + (1.0/rho)*At_mul_B(x, y.*alpha)

    # Precompute (Q+D)_ii:
    QD = Array{Float64,1}(s)
    for i in 1:s
        xi = vec(x[i,:])
        QD[i] = (0.5/C)
        QD[i] += vecdot(xi,xi) # doesnt have to be multiplied by y[i]*y[i] since this product is always 1
    end

    # the active set as a set of booleans.
    # if an entry is true, then the corresponding point
    # is active
    active = trues(s)
    num_active = s

    PGmax_old = Inf
    PGmin_old = -Inf

    t=0 # number of outer iterations

    while true
        t+=1

        PGmax_new = -Inf
        PGmin_new = Inf

        #TODO: permutate data in active set

        for i in 1:s
            if shrinking == false || active[i] == true
                xi = vec(x[i,:])
                G = y[i]*vecdot(w,xi) - 1.0 + (0.5/C)*alpha[i]

                # project gradient into feasible set:
                PG = 0
                if alpha[i] == 0 # exact zero check is fine here
                    if G > PGmax_old
                        active[i] = false
                        num_active-=1
                        continue
                    elseif G < 0
                        PG = G
                        #PG = min(0,G)
                    end
                else
                    PG = G
                end

                PGmax_new = max(PGmax_new, PG)
                PGmin_new = min(PGmin_new, PG)

                # perform update of alpha and w:
                if abs(PG) > 1e-6
                    alpha_old = alpha[i]
                    alpha[i] = max(0, alpha[i] - G/QD[i])
                    # w update:
                    w = w + (alpha[i] - alpha_old)*y[i]*xi
                end

            end
        end # end inner for loop

        if t >= iterations
            return w
        end

        if PGmax_new - PGmin_new < eps
            if num_active == s
                return w
            else # reset active set to see if we have a solution for the full set of data points
                num_active = s
                PGmax_old = Inf
                PGmin_old = -Inf
                continue # this missing might have caused problems before
            end
        end

        PGmax_old = PGmax_new
        PGmin_old = PGmin_new

        if PGmax_old <= 0
            PGmax_old = Inf
        end

        if PGmin_old >= 0
            PGmin_old = -Inf
        end

    end # end outer while loop
end


text:
# BioSequences.jl
# =======
#
# A julia package for the representation and manipulation of biological sequences.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

module BioSequences

export
    NucleicAcid,
    DNA,
    RNA,
    DNA_A,
    DNA_C,
    DNA_G,
    DNA_T,
    DNA_M,
    DNA_R,
    DNA_W,
    DNA_S,
    DNA_Y,
    DNA_K,
    DNA_V,
    DNA_H,
    DNA_D,
    DNA_B,
    DNA_N,
    DNA_Gap,
    ACGT,
    ACGTN,
    RNA_A,
    RNA_C,
    RNA_G,
    RNA_U,
    RNA_M,
    RNA_R,
    RNA_W,
    RNA_S,
    RNA_Y,
    RNA_K,
    RNA_V,
    RNA_H,
    RNA_D,
    RNA_B,
    RNA_N,
    RNA_Gap,
    ACGU,
    ACGUN,
    isGC,
    iscompatible,
    isambiguous,
    iscertain,
    isgap,
    ispurine,
    ispyrimidine,
    Sequence,
    BioSequence,
    DNASequence,
    RNASequence,
    AminoAcidSequence,
    CharSequence,
    NucleicAcidSequence,
    SeqRecord,
    seqname,
    hasseqname,
    sequence,
    hassequence,
    metadata,
    @dna_str,
    @rna_str,
    @aa_str,
    @char_str,
    @biore_str,
    @prosite_str,
    @kmer_str,
    matched,
    captured,
    alphabet,
    gap,
    complement,
    complement!,
    reverse_complement,
    reverse_complement!,
    ungap,
    ungap!,
    mismatches,
    ispalindromic,
    hasambiguity,
    isrepetitive,
    ambiguous_positions,
    gc_content,
    SequenceGenerator,
    randdnaseq,
    randrnaseq,
    randaaseq,
    canonical,
    neighbors,
    eachkmer,
    each,
    Composition,
    composition,
    NucleicAcidCounts,
    Kmer,
    DNAKmer,
    RNAKmer,
    DNACodon,
    RNACodon,
    translate,
    ncbi_trans_table,
    AminoAcid,
    AA_A,
    AA_R,
    AA_N,
    AA_D,
    AA_C,
    AA_Q,
    AA_E,
    AA_G,
    AA_H,
    AA_I,
    AA_L,
    AA_K,
    AA_M,
    AA_F,
    AA_P,
    AA_S,
    AA_T,
    AA_W,
    AA_Y,
    AA_V,
    AA_O,
    AA_U,
    AA_B,
    AA_J,
    AA_Z,
    AA_X,
    AA_Term,
    AA_Gap,
    FASTA,
    FASTQ,
    TwoBit,
    ABIF,

    # Alphabets
    Alphabet,
    DNAAlphabet,
    RNAAlphabet,
    NucAlphs,
    AminoAcidAlphabet,
    CharAlphabet,

    # search
    ExactSearchQuery,
    ApproximateSearchQuery,
    approxsearch,
    approxsearchindex,
    approxrsearch,
    approxrsearchindex,
    PFM,
    PWM,
    maxscore,
    scoreat,

    ReferenceSequence,
    Demultiplexer,
    demultiplex,
    seqmatrix,
    majorityvote,
    tryread!,
    isfilled,
    MissingFieldException,
    MinHashSketch,
    minhash,
    Site,
    Certain,
    Ambiguous,
    Gap,
    Match,
    Mismatch,
    count_pairwise

import Automa
import Automa.RegExp: @re_str

using BioCore
using BioSymbols
using BufferedStreams: BufferedStreams, BufferedInputStream, BufferedOutputStream
import Combinatorics
import IndexableBitVectors
import IntervalTrees: IntervalValue
import Twiddle: enumerate_nibbles,
    nibble_mask,
    count_zero_nibbles,
    count_nonzero_nibbles,
    count_zero_bitpairs,
    count_nonzero_bitpairs
using Random

BioSymbols.gap(::Type{Char}) = '-'

include("alphabet.jl")
include("bitindex.jl")
include("sequence.jl")
include("bioseq/bioseq.jl")
include("bioseq/hash.jl")
include("bioseq/randseq.jl")
include("kmer.jl")
include("nmask.jl")
include("refseq.jl")
include("eachkmer.jl")
include("composition.jl")
include("geneticcode.jl")
include("demultiplexer.jl")

# Parsing file types
include("fasta/fasta.jl")
include("fastq/fastq.jl")
include("twobit/twobit.jl")
include("abif/abif.jl")

include("search/exact.jl")
include("search/approx.jl")
include("search/re.jl")
include("search/pwm.jl")

include("minhash.jl")

end  # module BioSequences


text:
"""
Gridap, grid-based approximation of PDEs in the Julia programming language

This module provides rich set of tools for the numerical solution of PDE, mainly based
on finite element methods.

The module is structured in the following sub-modules:

- [`Gridap.Helpers`](@ref)
- [`Gridap.Inference`](@ref)
- [`Gridap.Io`](@ref)
- [`Gridap.Algebra`](@ref)
- [`Gridap.TensorValues`](@ref)
- [`Gridap.Arrays`](@ref)
- [`Gridap.Fields`](@ref)
- [`Gridap.Polynomials`](@ref)
- [`Gridap.Integration`](@ref)
- [`Gridap.ReferenceFEs`](@ref)
- [`Gridap.CellData`](@ref)
- [`Gridap.Geometry`](@ref)
- [`Gridap.FESpaces`](@ref)
- [`Gridap.MultiField`](@ref)
- [`Gridap.Visualization`](@ref)

The exported names are:
$(EXPORTS)
"""
module Gridap

using DocStringExtensions

include("Helpers/Helpers.jl")

include("Inference/Inference.jl")

include("Io/Io.jl")

include("Algebra/Algebra.jl")

include("Arrays/Arrays.jl")

include("TensorValues/TensorValues.jl")

include("Fields/Fields.jl")

include("Polynomials/Polynomials.jl")

include("Integration/Integration.jl")

include("ReferenceFEs/ReferenceFEs.jl")

include("CellData/CellData.jl")

include("Geometry/Geometry.jl")

include("FESpaces/FESpaces.jl")

include("MultiField/MultiField.jl")

include("Visualization/Visualization.jl")

include("Exports.jl")

end # module


text:
"""
$(SIGNATURES)
Propagate dynamics with a line search (in-place)
"""
function forwardpass!(prob::Problem, solver::iLQRSolver, ΔV::Array, J_prev::Float64)
    # Pull out values from results
    X = prob.X; U = prob.U; X̄ = solver.X̄; Ū = solver.Ū

    J = Inf
    alpha = 1.0
    iter = 0
    z = -1.
    expected = 0.

    logger = current_logger()
    to = solver.stats[:timer]
    # @logmsg InnerIters :iter value=0
    # @logmsg InnerIters :cost value=J_prev
    while (z ≤ solver.opts.line_search_lower_bound || z > solver.opts.line_search_upper_bound) && J >= J_prev

        # Check that maximum number of line search decrements has not occured
        if iter > solver.opts.iterations_linesearch
            # set trajectories to original trajectory
            copyto!(X̄,X)
            copyto!(Ū,U)

            J = cost(prob.obj, X̄, Ū, get_dt_traj(prob,Ū))

            z = 0.
            alpha = 0.0
            expected = 0.

            # @logmsg InnerLoop "Max iterations (forward pass)"
            regularization_update!(solver,:increase) # increase regularization
            solver.ρ[1] += solver.opts.bp_reg_fp
            break
        end

        # Otherwise, rollout a new trajectory for current alpha
        @timeit to "rollout" flag = rollout!(prob,solver,alpha)

        # Check if rollout completed
        if ~flag
            # Reduce step size if rollout returns non-finite values (NaN or Inf)
            # @logmsg InnerIters "Non-finite values in rollout"
            iter += 1
            alpha /= 2.0
            continue
        end

        # Calcuate cost
        @timeit to "cost" J = cost(prob.obj, X̄, Ū, get_dt_traj(prob,Ū))   # Unconstrained cost

        expected = -alpha*(ΔV[1] + alpha*ΔV[2])
        if expected > 0
            z  = (J_prev - J)/expected
        else
            @logmsg InnerIters "Non-positive expected decrease"
            z = -1
        end

        iter += 1
        alpha /= 2.0

        # Log messages
        # @logmsg InnerIters :iter value=iter
        # @logmsg InnerIters :α value=2*alpha
        # @logmsg InnerIters :cost value=J
        # @logmsg InnerIters :z value=z

    end  # forward pass loop

    # @logmsg InnerLoop :cost value=J
    # @logmsg InnerLoop :dJ value=J_prev-J
    @logmsg InnerLoop :expected value=expected
    @logmsg InnerLoop :z value=z
    @logmsg InnerLoop :α value=2*alpha
    @logmsg InnerLoop :ρ value=solver.ρ[1]

    if J > J_prev
        error("Error: Cost increased during Forward Pass")
    end

    return J
end

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
# This file inclues utilities.

"""
    $TYPEDSIGNATURES

Updates the diagonal elements of `mat` so that the `mat` is a diffusive matric is a zero row sum matrix.
"""
makediffusive!(mat::AbstractMatrix) = (mat[diagind(mat)] -= sum(mat, dims = 2); mat)

"""
    $TYPEDSIGNATURES

Scales `mat` by making its diagonale elements 1.
"""
normalize!(mat::AbstractMatrix) = (mat ./= diag(mat); mat)