text:
module DAT
export mapCube
import ..Cubes
using ..YAXTools
using Distributed: RemoteChannel, nworkers,pmap,
  @everywhere, workers, remotecall_fetch,
  remote_do, myid, nprocs, RemoteException, remotecall
import ..Cubes: cubechunks, iscompressed, chunkoffset,
  CubeAxis, YAXArray,
  caxes, YAXSlice
import ..Cubes.Axes: AxisDescriptor, axname, ByInference, axsym,
  getOutAxis, getAxis, findAxis
import ..Datasets: Dataset, createdataset
import ...YAXArrays
import ...YAXArrays.workdir
import ProgressMeter: Progress, next!, progress_pmap
using YAXArrayBase
using Dates
global const debugDAT=[false]
#TODO use a logging package
macro debug_print(e)
  debugDAT[1] && return(:(println($e)))
  :()
end

include("registration.jl")

"""
Internal representation of an input cube for DAT operations
"""
mutable struct InputCube{N}
  cube::Any   #The input data cube
  desc::InDims               #The input description given by the user/registration
  axesSmall::Vector{CubeAxis} #List of axes that were actually selected through the desciption
  icolon::Vector{Int}
  colonperm::Union{Vector{Int},Nothing}
  loopinds::Vector{Int}        #Indices of loop axes that this cube does not contain, i.e. broadcasts
  cachesize::Vector{Int}     #Number of elements to keep in cache along each axis
end

function InputCube(c, desc::InDims)
  iaxessmall = findAxis.(desc.axisdesc,Ref(c))
  axlist = caxes(c)
  axesSmall = map(i->axlist[i], iaxessmall)
  colonperm = issorted(iaxessmall) ? nothing : collect(Base.invperm(sortperm(collect(iaxessmall))))
  iaxessmall = collect(iaxessmall)
  any(isequal(nothing),axesSmall) && error("One of the input axes not found in input cubes")
  InputCube{ndims(c)}(c,desc,collect(CubeAxis,axesSmall),iaxessmall,colonperm,Int[],Int[])
end
geticolon(ic::InputCube) = ic.icolon
createworkarrays(T,s,ntr)=[Array{T}(undef,s...) for i=1:ntr]

"""
Internal representation of an output cube for DAT operations
"""
mutable struct OutputCube
  cube::Any #The actual outcube cube, once it is generated
  cube_unpermuted::Any #The unpermuted output cube
  desc::OutDims                 #The description of the output axes as given by users or registration
  axesSmall::Array{CubeAxis}       #The list of output axes determined through the description
  allAxes::Vector{CubeAxis}        #List of all the axes of the cube
  loopinds::Vector{Int}              #Index of the loop axes that are broadcasted for this output cube
  innerchunks
  outtype
end

const InOutCube = Union{InputCube,OutputCube}

getsmallax(c::InOutCube)=c.axesSmall
geticolon(c::OutputCube)=1:length(c.axesSmall)
getAxis(desc,c::InOutCube) = getAxis(desc,c.cube)

function getworkarray(c::InOutCube,ntr)
  wa = createworkarrays(eltype(c.cube),length.(c.axesSmall),ntr)
  map(wa) do w
    wrapWorkArray(c.desc.artype,w,c.axesSmall)
  end
end

function interpretoutchunksizes(desc,axesSmall,incubes)
  if desc.chunksize == :max
    map(ax->axname(ax)=>length(ax),axesSmall)
  elseif desc.chunksize == :input
    map(axesSmall) do ax
      for cc in incubes
        i = findAxis(axname(ax),cc)
        if i !== nothing
          return axname(ax)=>min(length(ax),cubechunks(cc)[i])
        end
      end
      return axname(ax)=>length(ax)
    end
  else
    desc.chunksize
  end
end

getOutAxis(desc::Tuple,inAxes,incubes,pargs,f)=map(i->getOutAxis(i,inAxes,incubes,pargs,f),desc)
function OutputCube(desc::OutDims,inAxes::Vector{CubeAxis},incubes,pargs,f)
  axesSmall     = getOutAxis(desc.axisdesc,inAxes,incubes,pargs,f)
  outtype       = getOuttype(desc.outtype,incubes)
  innerchunks   = interpretoutchunksizes(desc,axesSmall,incubes)
  OutputCube(nothing,nothing,desc,collect(CubeAxis,axesSmall),CubeAxis[],Int[],innerchunks,outtype)
end

"""
Configuration object of a DAT process. This holds all necessary information to perform the calculations.
It contains the following fields:

- `incubes::NTuple{NIN,InputCube}` The input data cubes
- `outcube::NTuple{NOUT,OutputCube}` The output data cubes
allInAxes     :: Vector
LoopAxes      :: Vector
ispar         :: Bool
loopcachesize :: Vector{Int}
max_cache
fu
inplace      :: Bool
include_loopvars:: Bool
ntr
addargs
kwargs
"""
mutable struct DATConfig{NIN,NOUT}
  "The input data cubes"
  incubes       :: NTuple{NIN,InputCube}

  outcubes      :: NTuple{NOUT,OutputCube}
  allInAxes     :: Vector
  LoopAxes      :: Vector
  ispar         :: Bool
  loopcachesize :: Vector{Int}
  max_cache
  fu
  inplace      :: Bool
  include_loopvars:: Bool
  ntr
  addargs
  kwargs
end
function DATConfig(cdata,indims,outdims,inplace,max_cache,fu,ispar,include_loopvars,nthreads,addargs,kwargs)

  isa(indims,InDims) && (indims=(indims,))
  isa(outdims,OutDims) && (outdims=(outdims,))
  length(cdata)==length(indims) || error("Number of input cubes ($(length(cdata))) differs from registration ($(length(indims)))")
  incubes  = ([InputCube(o[1],o[2]) for o in zip(cdata,indims)]...,)
  allInAxes = vcat([ic.axesSmall for ic in incubes]...)
  outcubes = ((map(1:length(outdims),outdims) do i,desc
     OutputCube(desc,allInAxes,cdata,addargs,fu)
  end)...,)

  DATConfig(
    incubes,
    outcubes,
    allInAxes,
    CubeAxis[],                                 # LoopAxes
    ispar,
    Int[],
    max_cache,                                  # max_cache
    fu,                                         # fu                                      # loopcachesize
    inplace,                                    # inplace
    include_loopvars,
    nthreads,
    addargs,                                    # addargs
    kwargs
  )
end

getOuttype(outtype::Int,cdata)=eltype(cdata[outtype])
function getOuttype(outtype::DataType,cdata)
  outtype
end

mapCube(fu::Function,cdata,addargs...;kwargs...)=mapCube(fu,(cdata,),addargs...;kwargs...)

function mapCube(f, in_ds::Dataset, addargs...; indims=InDims(), outdims=OutDims(), inplace=true, kwargs...)
    allars = values(in_ds.cubes)
    allaxes = collect(values(in_ds.axes))
    arnames = keys(in_ds.cubes)
    sarnames = (arnames...,)
    any(ad->findAxis(ad,allaxes)===nothing,indims.axisdesc) && error("One of the Dimensions does not exist in Dataset")

    idar = collect(indims.axisdesc)
    allindims = map(allars) do c
        idshort = filter(idar) do ad
            findAxis(ad,c)!==nothing
        end
        InDims((idshort...,), indims.artype, indims.procfilter)
    end
    isa(outdims, OutDims) || error("Only one output cube currently supported for datasets")
    isempty(addargs) || error("Additional arguments currently not supported for datasets, use kwargs instead")
    if inplace
        fnew = let arnames=collect(arnames), f=f
            function dsfun(xout,xin...; kwargs...)
                incubes = NamedTuple{sarnames, typeof(xin)}(xin)
                f(xout,incubes; kwargs...)
            end
        end
    else
        fnew = let arnames=collect(arnames), f=f
            function dsfun(xin...; kwargs...)
                incubes = NamedTuple{sarnames, typeof(xin)}(xin)
                f(incubes; kwargs...)
            end
        end
    end
    allcubes = collect(values(in_ds.cubes))
    mapCube(fnew,(allcubes...,);indims=allindims,outdims=outdims, inplace=inplace, kwargs...)
end

import Base.mapslices
function mapslices(f,d::Union{YAXArray, Dataset},addargs...;dims,kwargs...)
    isa(dims,String) && (dims=(dims,))
    mapCube(f,d,addargs...;indims = InDims(dims...),outdims = OutDims(ByInference()),inplace=false,kwargs...)
end

"""
    mapCube(fun, cube, addargs...;kwargs)

Map a given function `fun` over slices of the data cube `cube`.

### Keyword arguments

* `max_cache=1e7` maximum size of blocks that are read into memory, defaults to approx 10Mb
* `indims::InDims` List of input cube descriptors of type [`InDims`](@ref) for each input data cube
* `outdims::OutDims` List of output cube descriptors of type [`OutDims`](@ref) for each output cube
* `inplace` does the function write to an output array inplace or return a single value> defaults to `true`
* `ispar` boolean to determine if parallelisation should be applied, defaults to `true` if workers are available.
* `showprog` boolean indicating if a ProgressMeter shall be shown
* `include_loopvars` boolean to indicate if the varoables looped over should be added as function arguments
* `loopchunksize` determines the chunk sizes of variables which are looped over, a dict
* `kwargs` additional keyword arguments passed to the inner function

The first argument is always the function to be applied, the second is the input cube or
a tuple input cubes if needed.
"""
function mapCube(fu::Function,
    cdata::Tuple,addargs...;
    max_cache=YAXArrays.YAXDefaults.max_cache[],
    indims=InDims(),
    outdims=OutDims(),
    inplace=true,
    ispar=nprocs()>1,
    debug=false,
    include_loopvars=false,
    showprog=true,
    nthreads=ispar ? Dict(i=>remotecall_fetch(Threads.nthreads,i) for i in workers()) : [Threads.nthreads()] ,
    loopchunksize = Dict(),
    kwargs...)

  #Translate slices
  if any(i->isa(i,YAXSlice),cdata)
    inew = map(cdata) do d
      isa(d,YAXSlice) ? InDims(axname.(d.sliceaxes[2])...) : InDims()
    end
    cnew = map(i->isa(i,YAXSlice) ? i.c : i, cdata)
    return mapCube(fu, cnew, addargs...;
      indims=inew, outdims=outdims, inplace=inplace, ispar=ispar,
      debug=debug, include_loopvars=include_loopvars, showprog=showprog,
      nthreads=nthreads, kwargs...)
  end
  @debug_print "Generating DATConfig"
  dc=DATConfig(cdata,indims,outdims,inplace,
    max_cache,fu,ispar,include_loopvars,nthreads,addargs,kwargs)
  @debug_print "Analysing Axes"
  analyzeAxes(dc)
  @debug_print "Calculating Cache Sizes"
  getCacheSizes(dc,loopchunksize)
  @debug_print "Generating Output Cube"
  generateOutCubes(dc)
  @debug_print "Running main Loop"
  debug && return(dc)
  runLoop(dc,showprog)
  @debug_print "Finalizing Output Cube"

  if length(dc.outcubes)==1
    return dc.outcubes[1].cube_unpermuted
  else
    return (map(i->i.cube_unpermuted,dc.outcubes)...,)
  end

end

function makeinplace(f)
  (args...;kwargs...)->begin
    first(args).=f(Base.tail(args)...;kwargs...)
    nothing
  end
end

function getchunkoffsets(dc::DATConfig)
  co = zeros(Int,length(dc.LoopAxes))
  lc = dc.loopcachesize
  for ic in dc.incubes
    for (ax,cocur,cs) in zip(caxes(ic.cube),chunkoffset(ic.cube),cubechunks(ic.cube))
      ii = findAxis(ax,dc.LoopAxes)
      if !isa(ii,Nothing) && iszero(co[ii]) && cocur>0 && mod(lc[ii],cs)==0
        co[ii]=cocur
      end
    end
  end
  (co...,)
end

updatears(clist,r,f,caches) = foreach(clist, caches) do ic, ca
  updatear(f,r, ic.cube,geticolon(ic),ic.loopinds, ca )
end
getindsall(indscol, loopinds, rfunc, colfunc = _->Colon()) = getindsall((),1,(sort(indscol)...,),(loopinds...,),rfunc,colfunc)
function getindsall(indsall,inow,indscol,loopinds,rfunc,colfunc)
  if !isempty(indscol) && first(indscol) == inow
    getindsall((indsall..., colfunc(inow)),inow+1,Base.tail(indscol),loopinds,rfunc,colfunc)
  else
    getindsall((indsall..., rfunc(first(loopinds))),inow+1,indscol,Base.tail(loopinds),rfunc,colfunc)
  end
end
getindsall(indsall,inow,::Tuple{},::Tuple{},r,c) = indsall

function updatear(f,r,cube,indscol,loopinds,cache)
  indsall = getindsall(indscol,loopinds,i->r[i])
  l2 = map((i,s)->isa(i,Colon) ? s : length(i),indsall,size(cache))
  if size(cache) != l2
    hinds = map((i,s)->isa(i,Colon) ? (1:s) : 1:length(i),indsall,size(cache))
    if f == :read
      cache[hinds...] = getdata(cube)[indsall...]
    else
      getdata(cube)[indsall...] = cache[hinds...]
    end
  else
    if f == :read
      d = getdata(cube)[indsall...]
      cache[:] = d
    else
      getdata(cube)[indsall...] = cache
    end
  end
end
updateinars(dc,r,incaches)=updatears(dc.incubes,r,:read,incaches)
writeoutars(dc,r, outcaches)=updatears(dc.outcubes,r,:write,outcaches)

function loopworker(dcchan::RemoteChannel, ranchan, reschan)
    dc = try
      take!(dcchan)
    catch e
      println("Error serializing DATConfig, make sure all required package are loaded on all workers. ")
      put!(reschan,e)
      return
    end
    loopworker(dc, ranchan, reschan)
    return
end
function loopworker(dc::DATConfig, ranchan, reschan)
  incaches, outcaches, args = try
    getallargs(dc)
  catch e
    println("Error during initialization: ", e)
    put!(reschan, e)
  end
  try
    loopworker(dc,ranchan, reschan, incaches, outcaches, args)
  catch e
    if e isa InvalidStateException
      return nothing
    elseif e isa RemoteException && e.captured.ex isa InvalidStateException
      return nothing
    else
      println("Error during running loop: ", e)
      put!(reschan, e)
    end
  end
end

function loopworker(dc,ranchan, reschan, incaches, outcaches, args)
  while true
    r = take!(ranchan)
    updateinars(dc,r,incaches)
    innerLoop(r,args...)
    writeoutars(dc,r,outcaches)
    @async put!(reschan,r)
  end
end

function moduleloadedeverywhere()
  try
    isloaded = map(workers()) do w
      #We try calling a function defined inside this module, thi will error when YAXArrays is not loaded on the remote workers
      remotecall(()->true,w)
    end
    fetch.(isloaded)
  catch e
      return false
  end
  return true
end

function runLoop(dc::DATConfig,showprog)
  allRanges=distributeLoopRanges((dc.loopcachesize...,),(map(length,dc.LoopAxes)...,),getchunkoffsets(dc))
  chnlfunc() = Channel{Union{eltype(allRanges),Nothing}}(length(allRanges))
  outchanfunc() = Channel(length(allRanges))
  inchan, outchan, dcpass = if dc.ispar
    #Test if YAXArrays is loaded on all workers:
    moduleloadedeverywhere() || error("YAXArrays is not loaded on all workers. Please run `@everywhere using YAXArrays` to fix.")
    dcpass = RemoteChannel(()->Channel{DATConfig}(nworkers()))
    for i=1:nworkers()
      put!(dcpass, dc)
    end
    RemoteChannel(chnlfunc), RemoteChannel(outchanfunc), dcpass
  else
    chnlfunc(), outchanfunc(), dc
  end
  #Now distribute the jobs
  for r in allRanges
    put!(inchan,r)
  end
  #And start the workers
  if dc.ispar
    for p in workers()
      remote_do(YAXArrays.DAT.loopworker, p, dcpass, inchan, outchan)
    end
  else
    @async loopworker(dc, inchan, outchan)
  end
  showprog && (pm=Progress(length(allRanges)))
  for i in 1:length(allRanges)
    r = take!(outchan)
    if isa(r,Exception)
      throw(r)
    end
    showprog && next!(pm)
  end
  close(inchan)
  dc.outcubes
end

abstract type AxValCreator end
struct NoLoopAxes <:AxValCreator end
struct AllLoopAxes{S,V} <:AxValCreator
  loopsyms::S
  loopaxvals::V
end
AllLoopAxes(a) = AllLoopAxes(map(axsym,a), map(i->i.values,a))
getlaxvals(::NoLoopAxes,cI,offscur) = ()
getlaxvals(a::AllLoopAxes,cI,offscur) = (NamedTuple{a.loopsyms}(map((ax,i,of)->(i+of,ax[i+of]),a.loopaxvals,cI.I, offscur)),)

function getallargs(dc::DATConfig)
  incache, outcache = getCubeCache(dc)
  filters = map(ic->ic.desc.procfilter,dc.incubes)
  inworkar, outworkar = generateworkarrays(dc)
  axvals = if dc.include_loopvars
    lax = (dc.LoopAxes...,)
    AllLoopAxes(lax)
  else
    NoLoopAxes()
  end
  adda = dc.addargs
  kwa = dc.kwargs
  fu = if !dc.inplace
    makeinplace(dc.fu)
  else
    dc.fu
  end
  inarsbc = map(dc.incubes, incache) do ic, cache
    allax = getindsall(geticolon(ic), ic.loopinds, i->true)
    if ic.colonperm === nothing
      pa = PickAxisArray(cache,allax)
    else

      pa = PickAxisArray(cache,allax,ic.colonperm)
    end
  end
  outarsbc = map(dc.outcubes, outcache) do oc, cache
    allax = getindsall(1:length(oc.axesSmall), oc.loopinds, i->true)
    pa = PickAxisArray(cache,allax)
    pa
  end
  incache, outcache, (fu,inarsbc,outarsbc,filters,inworkar,outworkar,
  axvals,adda,kwa)
end

function getbackend(oc,ispar,max_cache)
  eltype=Union{oc.outtype,Missing}
  outsize=sizeof(eltype)*(length(oc.allAxes)>0 ? prod(map(length,oc.allAxes)) : 1)
  rt = oc.desc.backend
  if rt == :auto
    if ispar[] || outsize>max_cache
      rt = :zarr
    else
      rt = :array
    end
  end
  b = YAXArrayBase.backendlist[Symbol(rt)]
  if !allow_parallel_write(b)
    ispar[] = false
  end
  eltype,b
end

function generateOutCube(::Type{T},eltype,oc::OutputCube,loopcachesize,co;kwargs...) where T
  cs_inner = (map(length,oc.axesSmall)...,)
  cs = (cs_inner..., loopcachesize...)
  for (i,cc) in enumerate(oc.innerchunks)
    if cc !== nothing
      cs = Base.setindex(cs,cc,i)
    end
  end
  cube1, cube2 = createdataset(T, oc.allAxes; chunksize=cs, chunkoffset=co, kwargs...)
  oc.cube=cube1
  oc.cube_unpermuted = cube2
end
function generateOutCube(::Type{T},eltype,oc::OutputCube,loopcachesize,co;kwargs...) where T<:Array
  newsize=map(length,oc.allAxes)
  outar=Array{eltype}(undef,newsize...)
  map!(_->_zero(eltype),outar,1:length(outar))
  oc.cube = YAXArray(oc.allAxes,outar)
  oc.cube_unpermuted = oc.cube
end
_zero(T) = zero(T)
_zero(T::Type{<:AbstractString}) = convert(T,"")

function generateOutCubes(dc::DATConfig)
  co = getchunkoffsets(dc)
  foreach(dc.outcubes) do c
    co2 = (zeros(Int,length(c.axesSmall))...,co...)
    generateOutCube(c,Ref(dc.ispar),dc.max_cache,dc.loopcachesize,co2)
  end
end
function generateOutCube(oc::OutputCube,ispar::Ref{Bool},max_cache,loopcachesize,co)
  eltype,cubetype = getbackend(oc,ispar,max_cache)
  generateOutCube(cubetype,eltype,oc,loopcachesize,co;oc.desc.backendargs...)
end

function getCubeCache(dc::DATConfig)
  outcaches = map(i->allocatecachebuf(i,dc.loopcachesize),dc.outcubes)
  incaches = map(i->allocatecachebuf(i,dc.loopcachesize),dc.incubes)
  incaches, outcaches
end

function allocatecachebuf(ic::Union{InputCube,OutputCube},loopcachesize)
  s = size(ic.cube)
  indsall = getindsall(geticolon(ic), ic.loopinds, i->loopcachesize[i], i->s[i])
  zeros(eltype(ic.cube),indsall...)
end

function init_DATworkers()
  freshworkermodule()
end

function analyzeAxes(dc::DATConfig{NIN,NOUT}) where {NIN,NOUT}

  for cube in dc.incubes
    for a in caxes(cube.cube)
      in(a,cube.axesSmall) || in(a,dc.LoopAxes) || push!(dc.LoopAxes,a)
    end
  end
  length(dc.LoopAxes)==length(unique(map(typeof,dc.LoopAxes))) || error("Make sure that cube axes of different cubes match")
  for cube=dc.incubes
    myAxes = caxes(cube.cube)
    for (il,loopax) in enumerate(dc.LoopAxes)
      in(typeof(loopax),map(typeof,myAxes)) && push!(cube.loopinds,il)
    end
  end
  #Add output broadcast axes
  for outcube=dc.outcubes
    LoopAxesAdd=CubeAxis[]
    for (il,loopax) in enumerate(dc.LoopAxes)
      push!(outcube.loopinds,il)
      push!(LoopAxesAdd,loopax)
    end
    outcube.allAxes=CubeAxis[outcube.axesSmall;LoopAxesAdd]
    dold = outcube.innerchunks
    newchunks = Union{Int,Nothing}[nothing for _ in 1:length(outcube.allAxes)]
    for (k,v) in dold
      ii = findAxis(k,outcube.allAxes)
      if ii !== nothing
        newchunks[ii] = v
      end
    end
    outcube.innerchunks = newchunks
  end
  #And resolve names in chunk size dicts
  return dc
end

mysizeof(x)=sizeof(x)
mysizeof(x::Type{String})=1

"""
Function that compares two cache miss specifiers by their importance
"""
function cmpcachmisses(x1,x2)
  #First give preference to compressed misses
  if xor(x1.iscompressed,x2.iscompressed)
    return x1.iscompressed
  #Now compare the size of the miss multiplied with the inner size
  else
    return x1.cs * x1.innerleap > x2.cs * x2.innerleap
  end
end

function getCacheSizes(dc::DATConfig,loopchunksizes)

  inAxlengths      = Vector{Int}[Int.(length.(cube.axesSmall)) for cube in dc.incubes]
  inblocksizes     = map((x,T)->isempty(x) ? mysizeof(eltype(T.cube)) : prod(x)*mysizeof(eltype(T.cube)),inAxlengths,dc.incubes)
  inblocksize,imax = findmax(inblocksizes)
  outblocksizes    = map(C->length(C.axesSmall)>0 ? sizeof(C.outtype)*prod(map(Int∘length,C.axesSmall)) : 1,dc.outcubes)
  outblocksize     = length(outblocksizes) > 0 ? findmax(outblocksizes)[1] : 1
  #Now add cache miss information for each input cube to every loop axis
  cmisses= NamedTuple{(:iloopax,:cs, :iscompressed, :innerleap,:preventpar),Tuple{Int64,Int64,Bool,Int64,Bool}}[]
  userchunks = Dict{Int,Int}()
  for (k,v) in loopchunksizes
    ii = findAxis(k,dc.LoopAxes)
    if ii !== nothing
      userchunks[ii] = v
    end
  end
  foreach(dc.LoopAxes,1:length(dc.LoopAxes)) do lax,ilax
    haskey(userchunks,ilax) && return nothing
    for ic in dc.incubes
      ii = findAxis(lax,ic.cube)
      if !isa(ii,Nothing)
        inax = isempty(ic.axesSmall) ? 1 : prod(map(length,ic.axesSmall))
        push!(cmisses,(iloopax = ilax,cs = cubechunks(ic.cube)[ii],iscompressed = iscompressed(ic.cube), innerleap=inax,preventpar=false))
      end
    end
    for oc in dc.outcubes
      cs = oc.innerchunks
      ii = findAxis(lax,oc.allAxes)
      if !isa(ii,Nothing) && !isa(cs[ii],Nothing)
        innerleap = isempty(oc.axesSmall) ? 1 : prod(map(length,oc.axesSmall))
        push!(cmisses,(iloopax = ilax,cs = cs[ii],iscompressed = iscompressed(oc.cube), innerleap=innerleap, preventpar=true))
      end
    end
  end
  sort!(cmisses,lt=cmpcachmisses)
  loopcachesize,nopar    = getLoopCacheSize(max(inblocksize,outblocksize),map(length,dc.LoopAxes),dc.max_cache, cmisses, userchunks)
  for cube in dc.incubes
    cube.cachesize = map(length,cube.axesSmall)
    for (cs,loopAx) in zip(loopcachesize,dc.LoopAxes)
      in(typeof(loopAx),map(typeof,caxes(cube.cube))) && push!(cube.cachesize,cs)
    end
  end
  nopar && (dc.ispar = false)
  dc.loopcachesize=loopcachesize
  return dc
end

"Calculate optimal Cache size to DAT operation"
function getLoopCacheSize(preblocksize,loopaxlengths,max_cache,cmisses,userchunks)
  totcachesize=max_cache
  loopcachesize = ones(Int,length(loopaxlengths))
  #Go through user definitions
  incfac=totcachesize/preblocksize/prod(loopcachesize)
  incfac<1 && error("The requested slices do not fit into the specified cache. Please consider increasing max_cache")

  # Go through list of cache misses first and decide
  imiss = 1
  while imiss<=length(cmisses)
    il = cmisses[imiss].iloopax
    s = min(cmisses[imiss].cs,loopaxlengths[il])/loopcachesize[il]
    #Check if cache size is already set for this axis
    if loopcachesize[il]==1
      if s<incfac
        loopcachesize[il]=min(cmisses[imiss].cs,loopaxlengths[il])
        incfac=totcachesize/preblocksize/prod(loopcachesize)
      else
        ii=floor(Int,incfac)
        while ii>1 && rem(s,ii)!=0
          ii=ii-1
        end
        loopcachesize[il]=ii
        break
      end
    end
    imiss+=1
  end
  #Now second run to read multiple blocks at once
  if incfac >= 2
    for i=1:length(loopcachesize)
      haskey(userchunks,i) && continue
      imul = min(floor(Int,incfac),loopaxlengths[i]÷loopcachesize[i])
      loopcachesize[i] = loopcachesize[i] * imul
      incfac = incfac/imul
      incfac < 2 && break
    end
  end
  if imiss<length(cmisses)+1
    @warn "There are still cache misses"
    cmisses[imiss].iscompressed && @warn "There are compressed caches misses, you may want to use a different cube chunking"
  else
    #TODO continue increasing cache sizes on by one...
  end
  nopar = any(i->i.preventpar,cmisses[imiss:end])
  return loopcachesize, nopar
end

function distributeLoopRanges(block_size::NTuple{N,Int},loopR::NTuple{N,Int},co) where N
    allranges = map(block_size,loopR,co) do bs,lr,cocur
      collect(filter(!isempty,[max(1,i):min(i+bs-1,lr) for i in (1-cocur):bs:lr]))
    end
    Iterators.product(allranges...)
end

function generateworkarrays(dc::DATConfig)
  inwork = map(i->getworkarray(i,dc.ntr[myid()]),dc.incubes)
  outwork = map(i->getworkarray(i,dc.ntr[myid()]),dc.outcubes)
  inwork, outwork
end

macro innercode()
  esc(quote
    ithr = Threads.threadid()
    #Pick the correct array according to thread
    myinwork = map(i->i[ithr],inwork)
    myoutwork = map(i->i[ithr],outwork)
    #Copy data into work arrays
    foreach(myinwork,xinBC) do iw,x
      iw.=view(x,cI.I...)
    end
    #Apply filters
    mvs = map(docheck,filters,myinwork)
    if any(mvs)
      # Set all outputs to missing
      foreach(ow->fill!(ow,missing),myoutwork)
    else
      #Compute loop axis values if necessary
      laxval = getlaxvals(axvalcreator,cI,offscur)
      #Finally call the function
      f(myoutwork...,myinwork...,laxval...,addargs...;kwargs...)
    end
    #Copy data into output array
    foreach((iw,x)->view(x,cI.I...).=iw,myoutwork,xoutBC)
  end)
end

using DataStructures: OrderedDict
using Base.Cartesian
@noinline function innerLoop(loopRanges,f,xinBC,xoutBC,filters,
  inwork,outwork,axvalcreator,addargs,kwargs)
  offscur = map(i->(first(i)-1), loopRanges)
  if length(inwork[1])==1
    for cI in CartesianIndices(map(i->1:length(i),loopRanges))
      @innercode
    end
  else
    Threads.@threads for cI in CartesianIndices(map(i->1:length(i),loopRanges))
      @innercode
    end
  end
end

"Calculate an axis permutation that brings the wanted dimensions to the front"
function getFrontPerm(dc,dims)
  ax=caxes(dc)
  N=length(ax)
  perm=Int[i for i=1:length(ax)];
  iold=Int[]
  for i=1:length(dims) push!(iold,findAxis(dims[i],ax)) end
  iold2=sort(iold,rev=true)
  for i=1:length(iold) splice!(perm,iold2[i]) end
  perm=Int[iold;perm]
  return ntuple(i->perm[i],N)
end

include("dciterators.jl")
include("tablestats.jl")
include("CubeIO.jl")
end


text:
type VolatilityBumpInstrumentJacobianSwaption
  startIndex::Int
  endIndex::Int
end

type VolatilityBumpInstrumentJacobianCap
  startIndex::Int
  endIndex::Int
  strike::Float64
end

type VolatilityBumpInstrumentJacobian
  bumps::VegaBumpCollection
  swaptions::Vector{VolatilityBumpInstrumentJacobianSwaption}
  caps::Vector{VolatilityBumpInstrumentJacobianCap}
  computed::BitArray{1}
  allComputed::Bool
  derivatives::Vector{Vector{Float64}}
  onePercentBumps::Vector{Vector{Float64}}
  bumpMatrix::Matrix{Float64}
end

function VolatilityBumpInstrumentJacobian(bumps::VegaBumpCollection,
                                          swaptions::Vector{VolatilityBumpInstrumentJacobianSwaption},
                                          caps::Vector{VolatilityBumpInstrumentJacobianCap})

  swaptionsPlusCaps = length(swaptions) + length(caps)
  computed = falses(swaptionsPlusCaps)
  derivatives = [zeros(number_of_bumps(bumps)) for _ = 1:swaptionsPlusCaps]
  bumpMatrix = zeros(swaptionsPlusCaps, number_of_bumps(bumps))

  onePercentBumps = deepcopy(derivatives)
  allComputed = false

  return VolatilityBumpInstrumentJacobian(bumps, swaptions, caps, computed, allComputed, derivatives, onePercentBumps, bumpMatrix)
end

get_input_bumps(voljacobian::VolatilityBumpInstrumentJacobian) = voljacobian.bumps

function derivatives_volatility!(voljacobian::VolatilityBumpInstrumentJacobian, j::Int)
  j <= length(voljacobian.swaptions) + length(voljacobian.caps) || error("too high index passed to derivatives_volatility")

  if voljacobian.computed[j]
    return voljacobian.derivatives[j]
  end

  resize!(voljacobian.derivatives[j], number_of_bumps(voljacobian.bumps))
  resize!(voljacobian.onePercentBumps[j], number_of_bumps(voljacobian.bumps))

  sizesq = 0.0
  voljacobian.computed[j] = true

  initj = j

  if j <= length(voljacobian.swaptions)
    # it's a swaption
    thisPseudo = SwaptionPseudoDerivative(associated_model(voljacobian.bumps), voljacobian.swaptions[j].startIndex, voljacobian.swaptions[j].endIndex)

    for k = 1:number_of_bumps(voljacobian.bumps)
      v = 0.0
      for i = voljacobian.bumps.allBumps[k].stepBegin:voljacobian.bumps.allBumps[k].stepEnd
        fullDerivative = thisPseudo.volatilityDerivatives[i]
        for f = voljacobian.bumps.allBumps[k].factorBegin:voljacobian.bumps.allBumps[k].factorEnd, r = voljacobian.bumps.allBumps[k].rateBegin:voljacobian.bumps.allBumps[k].rateEnd
          v += fullDerivative[r, f]
        end
      end

      voljacobian.derivatives[j][k] = v
      sizesq += v * v
    end
  else
    # it's a cap
    j -= length(voljacobian.swaptions) # need to get back to index 1

    thisPseudo = CapPseudoDerivative(associated_model(voljacobian.bumps), voljacobian.caps[j].strike, voljacobian.caps[j].startIndex, voljacobian.caps[j].endIndex, 1.0)

    for k = 1:number_of_bumps(voljacobian.bumps)
      v = 0.0
      for i = voljacobian.bumps.allBumps[k].stepBegin:voljacobian.bumps.allBumps[k].stepEnd
        fullDerivative = thisPseudo.volatilityDerivatives[i]
        for f = voljacobian.bumps.allBumps[k].factorBegin:voljacobian.bumps.allBumps[k].factorEnd, r = voljacobian.bumps.allBumps[k].rateBegin:voljacobian.bumps.allBumps[k].rateEnd
          v += fullDerivative[r, f]
        end
      end

      sizesq += v * v
      voljacobian.derivatives[initj][k] = v
    end
  end

  for k = 1:number_of_bumps(voljacobian.bumps)
    voljacobian.bumpMatrix[initj, k] = voljacobian.onePercentBumps[initj][k] = 0.01 * voljacobian.derivatives[initj][k] / sizesq
  end

  return voljacobian.derivatives[initj]
end

function get_all_one_percent_bumps!(voljacobian::VolatilityBumpInstrumentJacobian)
  if ~voljacobian.allComputed
    for i = 1:length(voljacobian.swaptions) + length(voljacobian.caps)
      derivatives_volatility!(voljacobian, i)
    end
  end

  voljacobian.allComputed = true

  return voljacobian.bumpMatrix
end

type OrthogonalizedBumpFinder
  derivativesProducer::VolatilityBumpInstrumentJacobian
  multiplierCutoff::Float64
  tolerance::Float64
end

OrthogonalizedBumpFinder(bumps::VegaBumpCollection,
                        swaptions::Vector{VolatilityBumpInstrumentJacobianSwaption},
                        caps::Vector{VolatilityBumpInstrumentJacobianCap},
                        multiplierCutoff::Float64,
                        tolerance::Float64) = OrthogonalizedBumpFinder(VolatilityBumpInstrumentJacobian(bumps, swaptions, caps), multiplierCutoff, tolerance)

function get_vega_bumps!(obf::OrthogonalizedBumpFinder, theBumps::Vector{Vector{Matrix{Float64}}})
 # todo

  projector = OrthogonalProjection(get_all_one_percent_bumps!(obf.derivativesProducer), obf.multiplierCutoff, obf.tolerance)
  numberRestrictedBumps = projector.numberValidVectors
  marketmodel = associated_model(get_input_bumps(obf.derivativesProducer))
  evolution = marketmodel.evolution

  numberSteps = number_of_steps(evolution)
  numberRates = evolution.numberOfRates
  factors = marketmodel.numberOfFactors

  resize!(theBumps, numberSteps)

  modelMatrix = zeros(numberRates, factors)

  for i in eachindex(theBumps)
    theBumps[i] = Matrix[copy(modelMatrix) for i = 1:numberRestrictedBumps]
  end

  # for i = 1:numberSteps, j = 1:numberRestrictedBumps
  #   theBumps[i][j] = copy(modelMatrix)
  # end

  bumpClusters = get_input_bumps(obf.derivativesProducer).allBumps

  bumpIndex = 1
  for inst in eachindex(projector.validVectors)
    if projector.validVectors[inst]
      for cluster in eachindex(bumpClusters)
        magnitude = get_vector(projector, inst)[cluster]
        for s = bumpClusters[cluster].stepBegin:bumpClusters[cluster].stepEnd, r = bumpClusters[cluster].rateBegin:bumpClusters[cluster].rateEnd, factor = bumpClusters[cluster].factorBegin:bumpClusters[cluster].factorEnd
          theBumps[s][bumpIndex][r, factor] = magnitude
        end
      end
      bumpIndex += 1
    end
  end
  return theBumps
end

ext: jl


text:
module DivConformingFESpacesTests

using Test
using Gridap.Geometry
using Gridap.ReferenceFEs
using Gridap.FESpaces
using Gridap.CellData
using Gridap.TensorValues
using Gridap.Fields
using Gridap.Io

function test_div_v_q_equiv(U,V,P,Q,Ω)
  v=get_fe_basis(V)
  u=get_trial_fe_basis(U)

  q=get_fe_basis(Q)
  p=get_trial_fe_basis(P)

  dΩ=Measure(Ω,1)
  dΩᵣ=Measure(Ω,1,ReferenceDomain())

  a1(p,v)=∫(divergence(v)*p)dΩ
  a2(p,v)=∫(DIV(v)*p)dΩᵣ

  tol=1.0e-12
  assem = SparseMatrixAssembler(P,V)
  data = collect_cell_matrix(P,V,a1(p,v))
  A1 = assemble_matrix(assem,data)
  data = collect_cell_matrix(P,V,a2(p,v))
  A2 = assemble_matrix(assem,data)
  @test norm(A1-A2) < tol

  a3(u,q)=∫(q*divergence(u))dΩ
  a4(u,q)=∫(q*DIV(u))dΩᵣ
  assem = SparseMatrixAssembler(U,Q)
  data = collect_cell_matrix(U,Q,a3(u,q))
  A3 = assemble_matrix(assem,data)
  data = collect_cell_matrix(U,Q,a4(u,q))
  A4 = assemble_matrix(assem,data)
  @test norm(A3-A4) < tol

  uh=FEFunction(U,rand(num_free_dofs(U)))
  l1(q)=∫(q*divergence(uh))dΩ
  l2(q)=∫(q*DIV(uh))dΩᵣ
  v1=assemble_vector(l1,Q)
  v2=assemble_vector(l2,Q)
  @test norm(v1-v2) < tol
end

domain =(0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

order = 1

u(x) = x

reffe = ReferenceFE(raviart_thomas,order)

V = TestFESpace(model,reffe,dirichlet_tags = [1,6])
test_single_field_fe_space(V)
U = TrialFESpace(V,u)

reffe = ReferenceFE(lagrangian,Float64,order)
Q = TestFESpace(model,reffe,conformity=:L2)
P = TrialFESpace(Q)

uh = interpolate(u,U)

e = u - uh

Ω = Triangulation(model)
dΩ = Measure(Ω,2*order)

el2 = sqrt(sum( ∫( e⋅e )*dΩ ))
@test el2 < 1.0e-10

test_div_v_q_equiv(U,V,P,Q,Ω)

#using Gridap.Visualization
#
#writevtk(Ω,"trian",nsubcells=10,cellfields=["uh"=>uh])

order = 1

reffe = ReferenceFE(TET,raviart_thomas,order)

domain =(0,1,0,1,0,1)
partition = (3,3,3)
model = simplexify(CartesianDiscreteModel(domain,partition))

labels = get_face_labeling(model)
dir_tags = Array{Integer}(undef,0)

V = FESpace(model,reffe,conformity=DivConformity())
U = TrialFESpace(V,u)

reffe = ReferenceFE(lagrangian,Float64,order)
Q = TestFESpace(model,reffe,conformity=:L2)
P = TrialFESpace(Q)

v(x) = VectorValue(-0.5*x[1]+1.0,-0.5*x[2],-0.5*x[3])
vh = interpolate(v,V)

e = v - vh

Ω = Triangulation(model)
dΩ = Measure(Ω,2*order)

el2 = sqrt(sum( ∫( e⋅e )*dΩ ))
@test el2 < 1.0e-10

test_div_v_q_equiv(U,V,P,Q,Ω)

#using Gridap.Visualization
#
#writevtk(trian,"test",order=3,cellfields=["vh"=>vh])

# Tests on manifold

# Create domain
domain = (0,1,0,1,0,1)
cells  = (2,2,2)
model  = CartesianDiscreteModel(domain,cells)

# Restrict model to cube surface (using new BoundaryDiscreteModel)
labels = get_face_labeling(model)
bgface_to_mask = get_face_mask(labels,"boundary",2)
Γface_to_bgface = findall(bgface_to_mask)
Dc2Dp3model = BoundaryDiscreteModel(Polytope{2},model,Γface_to_bgface)

order  = 0
degree = 1

reffe_rt = ReferenceFE(raviart_thomas,Float64,order)
V  = FESpace(Dc2Dp3model, reffe_rt ; conformity=:HDiv)
U = TrialFESpace(V,u)
reffe = ReferenceFE(lagrangian,Float64,order)
Q = TestFESpace(Dc2Dp3model,reffe,conformity=:L2)
P = TrialFESpace(Q)
uh = FEFunction(V,rand(num_free_dofs(V)))
vh = interpolate_everywhere(uh,V)

Ω = Triangulation(Dc2Dp3model)
dΩ = Measure(Ω,2*order)
e=sqrt(sum(∫((uh-vh)⋅(uh-vh))dΩ))
@test e < 1.0e-12

test_div_v_q_equiv(U,V,P,Q,Ω)

end # module


text:
abstract AbstractDataExporter
export AbstractDataExporter

write_data(fp::FEProblem, exporter::AbstractDataExporter, n_print::Int) = error("Not implemented")

@lazymod VTKExportMod "vtkexport_xml.jl"


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 12a6dc11-0856-43b7-8460-09fb2c28e27b
begin
    import Pkg
    Pkg.activate(mktempdir())
	Pkg.add("PlutoUI")
    using Plots, PlutoUI
end

# ╔═╡ 09a781eb-a30f-468d-a067-dfa48c4c0068
md"# Trigonometric dance"

# ╔═╡ 9c31c3fd-354d-4a98-802a-61962e16ae9d
md"""
Triangles and circles might not to appear to get along. Triangles are full of straight lines and sharp angles. Circles bendy and are smooth. But Trigonometry is the music that allows them to dance together, despite their apparent differences.
"""

# ╔═╡ 0865f3d4-9e25-11eb-1fd0-0104f43e7534
side_opposite(angle_in_rad) = sin(angle_in_rad)

# ╔═╡ bde04203-706d-4ed9-9915-2584783246be
md"""
parametric version of cos of radian
"""

# ╔═╡ 0d8618b9-be99-40a8-8151-95bbe8abca21
side_adjacent(angle_in_rad) = cos(angle_in_rad)

# ╔═╡ 3edef3a8-a1fb-42bf-9459-1db16b303930
@bind n Slider(0.000 : 0.001: 2π)

# ╔═╡ c868ec19-5172-4223-b876-2645e1eac057
n

# ╔═╡ fb57b26f-5f5e-481c-8cab-627e7de2b468
plot(side_adjacent, side_opposite, 0, n, aspect_ratio=1,xlims=(-1.50, 1.50), ylims=(-1.50,1.50))

# ╔═╡ 7905c9b2-299c-49df-986d-399c8e5993e5
cos(n)

# ╔═╡ a9f208f0-15e2-4ab4-adda-d9cc49300cb5
sin(n)

# ╔═╡ 3a0d9ec7-6a8c-47d5-940c-fbd2c5003a3c
plot!([0,cos(n)],[0,0])

# ╔═╡ 9fcf0f7a-5f11-4497-93f5-11bdd694f4c1
plot!([cos(n),cos(n)],[0,sin(n)])

# ╔═╡ 8c330e74-f42c-4857-8e29-483b05f9d861
plot!([0,cos(n)],[0,sin(n)])

# ╔═╡ Cell order:
# ╟─09a781eb-a30f-468d-a067-dfa48c4c0068
# ╠═9c31c3fd-354d-4a98-802a-61962e16ae9d
# ╠═12a6dc11-0856-43b7-8460-09fb2c28e27b
# ╠═c868ec19-5172-4223-b876-2645e1eac057
# ╠═0865f3d4-9e25-11eb-1fd0-0104f43e7534
# ╟─bde04203-706d-4ed9-9915-2584783246be
# ╠═0d8618b9-be99-40a8-8151-95bbe8abca21
# ╠═fb57b26f-5f5e-481c-8cab-627e7de2b468
# ╠═7905c9b2-299c-49df-986d-399c8e5993e5
# ╠═a9f208f0-15e2-4ab4-adda-d9cc49300cb5
# ╠═3a0d9ec7-6a8c-47d5-940c-fbd2c5003a3c
# ╠═9fcf0f7a-5f11-4497-93f5-11bdd694f4c1
# ╠═3edef3a8-a1fb-42bf-9459-1db16b303930
# ╠═8c330e74-f42c-4857-8e29-483b05f9d861


text:
using Base.Meta

import Base: complex, conj, real, imag, max, min, abs, round, sign, floor, hypot
import Base: &, |, $, .>, .>=, .<, .<=, !, .==, .!=, ^, .^, /, ./
import Base: +, .+, -, .-, *, .*, /, ./, %, .%, <<, .<<, >>, .>>, ^, .^
import Base: sin, cos, tan, sinh, cosh, tanh, asin, acos, atan, cbrt, erf, erfc,
             exp, expm1, factorial, lgamma, log, log10, log1p, sqrt, gamma, log2,
             atan2

export sigmoid

for (op,fn) in ((:+,:af_add), (:.+,:af_add), (:-,:af_sub), (:.-,:af_sub), (:.*,:af_mul),
                (:./,:af_div), (:%, :af_mod), (:.%, :af_mod), (:<<, :af_bitshiftl),
                (:.<<, :af_bitshiftl),(:>>, :af_bitshiftr), (:.>>, :af_bitshiftr),
                (:^, :af_pow), (:.^, :af_pow ))

    @compat @eval function $op{T,S}(a::AFArray{T}, b::AFArray{S}; batched = true)
        ptr = new_ptr()
        eval($(quot(fn)))(ptr, a, b, batched)
        AFArray{af_promote(T,S)}(ptr[])
    end

end

@compat Base.:-(a::AFArray) = 0 - a
@compat Base.:+(a::AFArray{Bool}, v::Bool) = +(a,Int(v))
@compat Base.:+(v::Bool, a::AFArray{Bool}) = +(a,v)
@compat Base.:-(a::AFArray{Bool}, v::Bool) = -(a,Int(v))
@compat Base.:-(v::Bool, a::AFArray{Bool}) = -(a,v)

@compat ^{T<:Real}(a::AFArray{T}, v::Integer) = ^(a, Real(v))
@compat ^{T<:Real}(a::AFArray{Complex{T}}, v::Integer) = ^(a, Real(v))
@compat .^{T<:Real}(v::Irrational{:e}, a::AFArray{T}) = ^(Real(e), a)
@compat .^{T<:Real}(v::Irrational{:e}, a::AFArray{Complex{T}}) = ^(Real(e), a)

for (op,fn) in ((:+, :af_add), (:.+, :af_add), (:-, :af_sub), (:.-, :af_sub), (:*, :af_mul),
                (:.*, :af_mul), (:./, :af_div), (:%, :af_mod), (:.%, :af_mod),
                (:^, :af_pow), (:.^, :af_pow))

    @compat @eval function $op{T<:Real,S<:Real}(a::AFArray{T}, v::S)
        b = constant((af_promote(T,S))(v), size(a)...)
        ptr = new_ptr()
        eval($(quot(fn)))(ptr, a, b, true)
        AFArray{af_promote(T,S)}(ptr[])
    end

    @compat @eval function $op{T<:Real,S<:Real}(a::AFArray{T}, v::Complex{S})
        b = constant(Complex{af_promote(T,S)}(v), size(a)...)
        ptr = new_ptr()
        eval($(quot(fn)))(ptr, a, b, true)
        AFArray{Complex{af_promote(T,S)}}(ptr[])
    end

    @compat @eval function $op{T<:Real,S<:Real}(v::S, a::AFArray{T})
        b = constant((af_promote(T,S))(v), size(a)...)
        ptr = new_ptr()
        eval($(quot(fn)))(ptr, b, a, true)
        AFArray{af_promote(T,S)}(ptr[])
    end

    @compat @eval function $op{T<:Real,S<:Real}(v::Complex{S}, a::AFArray{T})
        b = constant(Complex{af_promote(T,S)}(v), size(a)...)
        ptr = new_ptr()
        eval($(quot(fn)))(ptr, b, a, true)
        AFArray{Complex{af_promote(T,S)}}(ptr[])
    end

    @compat @eval function $op{T<:Real,S<:Real}(a::AFArray{Complex{T}}, v::S)
        b = constant((af_promote(T,S))(v), size(a)...)
        ptr = new_ptr()
        eval($(quot(fn)))(ptr, a, b, true)
        AFArray{Complex{af_promote(T,S)}}(ptr[])
    end

    @compat @eval function $op{T<:Real,S<:Real}(a::AFArray{Complex{T}}, v::Complex{S})
        b = constant(Complex{af_promote(T,S)}(v), size(a)...)
        ptr = new_ptr()
        eval($(quot(fn)))(ptr, a, b, true)
        AFArray{Complex{af_promote(T,S)}}(ptr[])
    end

    @compat @eval function $op{T<:Real,S<:Real}(v::S, a::AFArray{Complex{T}})
        b = constant((af_promote(T,S))(v), size(a)...)
        ptr = new_ptr()
        eval($(quot(fn)))(ptr, b, a, true)
        AFArray{Complex{af_promote(T,S)}}(ptr[])
    end

    @compat @eval function $op{T<:Real,S<:Real}(v::Complex{S}, a::AFArray{Complex{T}})
        b = constant(Complex{af_promote(T,S)}(v), size(a)...)
        ptr = new_ptr()
        eval($(quot(fn)))(ptr, b, a, true)
        AFArray{Complex{af_promote(T,S)}}(ptr[])
    end
end
/{T,S<:Real}(a::AFArray{T}, v::S) = ./(a, v)
/{T,S<:Real}(a::AFArray{T}, v::Complex{S}) = ./(a, v)

for (op,fn) in ((:sin, :af_sin), (:cos, :af_cos), (:tan, :af_tan), (:asin, :af_asin),
                (:acos, :af_acos), (:atan, :af_atan), (:sinh, :af_sinh), (:cosh, :af_cosh),
                (:tanh, :af_tanh), (:asinh, :af_asinh), (:acosh, :af_acosh), (:atanh, :af_atanh),
                (:cbrt, :af_cbrt), (:erf, :af_erf), (:erfc, :af_erfc), (:exp, :af_exp),
                (:expm1, :af_expm1), (:factorial, :af_factorial), (:lgamma, :af_lgamma),
                (:log, :af_log), (:log10, :af_log10), (:log1p, :af_log1p), (:sqrt, :af_sqrt),
                (:gamma, :af_tgamma), (:log2, :af_log2))
    @compat @eval function $op(a::AFArray)
        out = new_ptr()
        eval($(quot(fn)))(out, a)
        AFArray{backend_eltype(out[])}(out[])
    end

end

for (op,fn) in ((:atan2, :af_atan2),)
    @eval @compat function $op{T,S}(a::AFArray{T}, b::AFArray{S}; batched = true)
        ptr = new_ptr()
        eval($(quot(fn)))(ptr, a, b, batched)
        AFArray{af_promote(T,S)}(ptr[])
    end

    @eval @compat function $op{T<:Real,S<:Real}(a::AFArray{T}, v::S)
        b = constant((af_promote(T,S))(v), size(a)...)
        ptr = new_ptr()
        eval($(quot(fn)))(ptr, a, b, true)
        AFArray{af_promote(T,S)}(ptr[])
    end

    @eval @compat function $op{T<:Real,S<:Real}(v::S, a::AFArray{T})
        b = constant((af_promote(T,S))(v), size(a)...)
        ptr = new_ptr()
        eval($(quot(fn)))(ptr, b, a, true)
        AFArray{af_promote(T,S)}(ptr[])
    end
end

function sigmoid(a::AFArray)
    out = new_ptr()
    af_sigmoid(out, a)
    AFArray{backend_eltype(out[])}(out[])
end

# Complex functions

function complex{T<:Real}(a::AFArray{T})
    out = new_ptr()
    af_cplx(out, a)
    AFArray{Complex{T}}(out[])
end

function complex{T<:Real}(a::AFArray{T}, b::AFArray{T}; batch = true)
    out = new_ptr()
    af_cplx2(out, a, b, batch)
    AFArray{Complex{T}}(out[])
end

function conj{T<:Complex}(a::AFArray{T})
    out = new_ptr()
    af_conjg(out, a)
    AFArray{T}(out[])
end

function real{T<:Real}(a::AFArray{Complex{T}})
    out = new_ptr()
    af_real(out, a)
    AFArray{T}(out[])
end

function real{T<:Real}(a::AFArray{T})
    out = new_ptr()
    af_real(out, a)
    AFArray{T}(out[])
end

function imag{T<:Real}(a::AFArray{Complex{T}})
    out = new_ptr()
    af_imag(out, a)
    AFArray{T}(out[])
end

function imag{T}(a::AFArray{T})
    out = new_ptr()
    af_imag(out, a)
    AFArray{T}(out[])
end

# Numeric functions

for (op,fn) in ((:max, :af_maxof), (:min, :af_minof))

    @eval function ($op){T,S}(a::AFArray{T}, b::AFArray{S}; batched = true)
        out = new_ptr()
        eval($fn)(out, a, b, true)
        AFArray{af_promote(T,S)}(out[])
    end

    @eval function ($op){T,S<:Real}(a::AFArray{T}, b::S)
        out = new_ptr()
        tmp = constant(b, size(a))
        eval($fn)(out, a, tmp, true)
        AFArray{af_promote(T,S)}(out[])
    end
    @eval ($op){T,S<:Real}(b::S, a::AFArray{T}) = max(a, b)

    @eval function ($op){T<:Complex,S<:Complex}(a::AFArray{T}, b::S)
        out = new_ptr()
        tmp = constant(b, size(a))
        eval($fn)(out, a, tmp, true)
        AFArray{af_promote(T,S)}(out[])
    end
    @eval ($op){T<:Complex,S<:Complex}(b::S, a::AFArray{T}) = max(a, b)

end

for (op,fn) in ((:abs, :af_abs), (:arg, :af_arg),
                (:ceil, :af_ceil), (:sign, :af_sign),
                (:floor, :af_floor), (:round, :af_round),
                (:trunc, :af_trunc))

    @eval function ($op){T}(a::AFArray{T})
        out = new_ptr()
        eval($fn)(out, a)
        AFArray{T}(out[])
    end

end

function hypot{T}(a::AFArray{T}, b::AFArray{T}; batched = true)
    out = new_ptr()
    af_hypot(out, a, b, batched)
    AFArray{T}(out[])
end

# Logical Operations

for (op,fn) in ((:&, :af_bitand),(:|, :af_bitor), (:$, :af_bitxor),
                (:.==, :af_eq), (:.!=, :af_neq), (:.>, :af_gt),
                (:.>=, :af_ge), (:.<, :af_lt), (:.<=, :af_le))

    @eval function ($op)(a::AFArray, b::AFArray; batched = true)
        out = new_ptr()
        eval($fn)(out, a, b, batched)
        AFArray{backend_eltype(out[])}(out[])
    end

    @eval function ($op)(a::AFArray, b::Real; batched = true)
        out = new_ptr()
        tmp = constant(b, size(a))
        eval($fn)(out, a, tmp, batched)
        AFArray{backend_eltype(out[])}(out[])
    end
    @eval ($op)(b::Real, a::AFArray) = ($op)(a, b)

end

function !{T}(a::AFArray{T})
    out = new_ptr()
    af_not(out, a)
    AFArray{T}(out[])
end


text:
"""
MCS
"""

export ftmcs, ftmcs!
using DD

function ftmcs(f::AbstractFaultTreeNode)
    top, = bdd(f)
    return ftmcs(top)
end

function ftmcs(f::AbstractDDNode{Tv,Ti}) where {Tv,Ti}
    cache = Dict{AbstractDDNode{Tv,Ti},Vector{Vector{Symbol}}}()
    return ftmcs!(f, cache)
end

function ftmcs!(f::DDVariable{Tv,Ti,2}, cache::Dict{AbstractDDNode{Tv,Ti},Vector{Vector{Symbol}}}) where {Tv,Ti}
    _minimalset(_ftmcs!(f, cache))
end

function _ftmcs!(f::DDVariable{Tv,Ti,2}, cache::Dict{AbstractDDNode{Tv,Ti},Vector{Vector{Symbol}}}) where {Tv,Ti}
    get(cache, f) do
        res = [_ftmcs!(x, cache) for x = f.nodes]
        res2 = [push!(copy(x), f.label) for x = res[2]]
        res1 = [x for x = res[1]]
        cache[f] = vcat(res1, res2)
    end
end

function _ftmcs!(f::DDValue{Tv,Ti}, cache::Dict{AbstractDDNode{Tv,Ti},Vector{Vector{Symbol}}}) where {Tv,Ti}
    (f.val == Tv(1)) ? Vector{Symbol}[Symbol[]] : Vector{Symbol}[]
end

function _remove!(x, xs)
    for y = xs
        isempty(setdiff(x, y)) && return false
    end
    return true
end

function _minimalset(args)
    xs = sort(args, by=x->length(x), rev=true)
    result = Vector{Symbol}[]
    while !isempty(xs)
        x = pop!(xs)
        push!(result, sort(x))
        xs = [y for y = xs if !issubset(x, y)]
    end
    sort(result)
end


text:
#
# Correctness Tests
#

using Test, DelimitedFiles, .MathConstants
using GeoEfficiency
const G = GeoEfficiency		# to qualify non exported objects

include("Helper.jl")

"""# UnExported

	@to_string ex

convert the expresion `ex` to a string.
"""
macro to_string(ex)
	:($(G.to_string(ex)))
end

const SourceFiles = [
					"Helper",
					"Error",
					"Input_Console",
					"Physics_Model",
					"Input_Batch",
    				"Calculations",
					"Output_Interface"
					]

@testset "GeoEfficiency" begin
	@test G.srcType === G.srcUnknown  		# the initial program condition
	@test typeofSrc() === G.srcUnknown  	# the initial program condition
	@test setSrcToPoint() === false      	# not defined, set to not point
	println("\n")

	@testset "$SourceFile" for SourceFile = SourceFiles
		@debug("Begin test", SourceFile)
		include("test_$SourceFile.jl")
		println("\n")
	end #testset_SourceFile

end #testset_GeoEfficiency

max_stars_repo_name: DrKrar/GeoEfficiency

max_issues_repo_name: DrKrar/GeoEfficiency

max_forks_repo_name: DrKrar/GeoEfficiency


text:
struct PolarPoint{T} <: AbstractSet{T}
    a::Vector{T}
end
dimension(h::PolarPoint) = length(h.a)
space_variables(::PolarPoint) = nothing

function scaling_function(h::PolarPoint)
    @assert dimension(h) == 2
    return (x, y) -> begin
        return h.a[1] * x + h.a[2] * y
    end
end

function _print_gauge_function(h::PolarPoint; digits=6)
    DynamicPolynomials.@polyvar x[1:2]
    print(" ")
    a = h.a
    if digits !== nothing
        a = round.(a, digits=6)
    end
    println(x' * a)
end


text:
using TrustPilotCalc
using Dates, JuliaDB

# load sample dataset and generate yearmon list
example_data = loadtable("../../TrustPlotCalc.jl-old/data/trustpilot_working.csv"; type_detect_rows = 500)
yearmon_list = unique(map(x-> round(x.review_created_at_utc, Dates.Month), example_data))
sort!(yearmon_list; rev = true)
pop!(yearmon_list) # removing last entry to avoid empty arrays

# task channels
c_output = Channel(length(yearmon_list))

function run_through(base_dataset, max_yearmon)
    while true
        result = playBack(base_dataset; max_yearmon = max_yearmon)
        put!(c_output, result)
    end
end

for i in 1:length(yearmon_list)
    @async run_through(example_data, yearmon_list[i])
end

[x for x in c_output]

# [print(playBack(example_data; max_yearmon = x).score, " | ", x, "\n") for x in yearmon_list]


text:
using DynamicGrids, Test

init = [1.0 4.0 7.0;
        2.0 5.0 8.0;
        3.0 6.0 9.0]

mask = Bool[0 1 0;
            0 1 1;
            1 1 0]

struct DoNothingRule <: Rule end
DynamicGrids.applyrule(::DoNothingRule, data, state, args...) = state

rules = Ruleset(DoNothingRule(); init=init, mask=mask)

output = ArrayOutput(init, 8)
sim!(output, rules; tspan=(1, 8))

@test output[1] == [1.0 4.0 7.0;
                    2.0 5.0 8.0;
                    3.0 6.0 9.0]

@test output[2] == [0.0 4.0 0.0;
                    0.0 5.0 8.0;
                    3.0 6.0 0.0]


text:
# distributions
@rule typeof(+)(:in1, Marginalisation) (m_out::UnivariateNormalDistributionsFamily, m_in2::UnivariateNormalDistributionsFamily) = begin
    min2, vin2 = mean_var(m_in2)
    mout, vout = mean_var(m_out)
    return NormalMeanVariance(mout - min2, vout + vin2)
end

@rule typeof(+)(:in1, Marginalisation) (m_out::MultivariateNormalDistributionsFamily, m_in2::MultivariateNormalDistributionsFamily) = begin
    min2, vin2 = mean_cov(m_in2)
    mout, vout = mean_cov(m_out)
    return MvNormalMeanCovariance(mout - min2, vout + vin2)
end

# PointMass
@rule typeof(+)(:in1, Marginalisation) (m_out::PointMass, m_in2::PointMass) = PointMass(mean(m_out) - mean(m_in2))

@rule typeof(+)(:in1, Marginalisation) (m_out::PointMass, m_in2::NormalMeanPrecision) = begin
    min2, win2 = mean_precision(m_in2)
    return NormalMeanPrecision(mean(m_out) - min2, win2)
end

@rule typeof(+)(:in1, Marginalisation) (m_out::NormalMeanPrecision, m_in2::PointMass) = begin
    mout, wout = mean_precision(m_out)
    return NormalMeanPrecision(mout - mean(m_in2), wout)
end

@rule typeof(+)(:in1, Marginalisation) (m_out::PointMass, m_in2::MvNormalMeanPrecision) = begin
    min2, win2 = mean_precision(m_in2)
    return MvNormalMeanPrecision(mean(m_out) - min2, win2)
end

@rule typeof(+)(:in1, Marginalisation) (m_out::MvNormalMeanPrecision, m_in2::PointMass) = begin
    mout, wout = mean_precision(m_out)
    return MvNormalMeanPrecision(mout - mean(m_in2), wout)
end

@rule typeof(+)(:in1, Marginalisation) (m_out::PointMass, m_in2::UnivariateNormalDistributionsFamily) = begin
    min2, vin2 = mean_var(m_in2)
    return NormalMeanVariance(mean(m_out) - min2, vin2)
end

@rule typeof(+)(:in1, Marginalisation) (m_out::UnivariateNormalDistributionsFamily, m_in2::PointMass) = begin
    mout, vout = mean_var(m_out)
    return NormalMeanVariance(mout - mean(m_in2), vout)
end

@rule typeof(+)(:in1, Marginalisation) (m_out::PointMass, m_in2::MultivariateNormalDistributionsFamily) = begin
    min2, vin2 = mean_cov(m_in2)
    return MvNormalMeanCovariance(mean(m_out) - min2, vin2)
end

@rule typeof(+)(:in1, Marginalisation) (m_out::MultivariateNormalDistributionsFamily, m_in2::PointMass) = begin
    mout, vout = mean_cov(m_out)
    return MvNormalMeanCovariance(mout - mean(m_in2), vout)
end

# specialized
@rule typeof(+)(:in1, Marginalisation) (m_out::MvNormalWeightedMeanPrecision{T1}, m_in2::MvNormalWeightedMeanPrecision{T2}) where { T1 <: LinearAlgebra.BlasFloat, T2 <: LinearAlgebra.BlasFloat } = begin

    min2, vin2 = mean_cov(m_in2)
    vout = cov(m_out)
    T = promote_type(T1, T2)
    BLAS.gemv!('N', -one(T), vout, weightedmean(m_out), one(T), min2)
    vin2 .+= vout
    return MvNormalMeanCovariance(min2, vin2)

end

@rule typeof(+)(:in1, Marginalisation) (m_out::MvNormalWeightedMeanPrecision{T1}, m_in2::PointMass) where { T1 } = begin

    ξout, wout = weightedmean_precision(m_out)
    ξin1 = wout*mean(m_in2)
    ξin1 .-= ξout
    T = promote_type(T1, eltype(m_in2))
    ξin1 .*= -one(T)
    return MvNormalWeightedMeanPrecision(ξin1, wout)

end

@rule typeof(+)(:in1, Marginalisation) (m_out::PointMass, m_in2::MvNormalWeightedMeanPrecision) = begin

    ξin2, win2 = weightedmean_precision(m_in2)
    ξin1 = win2*mean(m_out)
    ξin1 .-= ξin2
    return MvNormalWeightedMeanPrecision(ξin1, win2)

end


text:
"""
read_varint reads a variable integer from a stream
"""
 function read_varint(s::Base.GenericIOBuffer{Array{UInt8,1}})
     i = UInt8[]
     readbytes!(s, i, 1)
     if i == [0xfd]
         # 0xfd means the next two bytes are the number
         readbytes!(s, i, 2)
         return reinterpret(Int16, i)[1]
     elseif i == [0xfe]
         # 0xfe means the next four bytes are the number
         readbytes!(s, i, 4)
         return reinterpret(Int32, i)[1]
     elseif i == [0xff]
         # 0xff means the next eight bytes are the number
         readbytes!(s, i, 8)
         return reinterpret(Int64, i)[1]
     else
         # anything else is just the integer
         return reinterpret(Int8, i)[1]
     end
 end

"""
Encodes an integer as a varint
"""
 function encode_varint(n::Integer)
    if n < 0xfd
        return [UInt8(n)]
    elseif n < 0x10000
        return prepend!(int2bytes(n, 2, true), [0xfd])
    elseif n < 0x100000000
        return prepend!(int2bytes(n, 4, true), [0xfe])
    elseif n < 0x10000000000000000
        return prepend!(int2bytes(n, 8, true), [0xff])
    else
        error("Integer, ", n, " is too large")
    end
 end

"""
Double sha256 function
"""
function hash256(x::Array{UInt8, 1})
    return sha256(sha256(x))
end


text:
include("default_formulas.jl")

struct ReckonerInstance{R, T}
    aup::Function
    weight::Function
    challenge_window::Function
    skill::Function
    rating::Function
    eff_challenge::Function
    win_chances::Function
    rank_interval::Function
    display_rank::Function

    function ReckonerInstance{R, T}(;
                aup::Function = default_aup,
                weight::Function = default_weight,
                challenge_window::Function = default_challenge_window,
                skill::Function = default_skill,
                rating::Function = default_rating,
                eff_challenge::Function = default_eff_challenge,
                win_chances::Function = default_win_chances,
                rank_interval::Function = default_rank_interval,
                display_rank::Function = default_display_rank
                )::ReckonerInstance where {R, T}

        new(aup, weight, challenge_window, skill, rating, eff_challenge, win_chances, rank_interval, display_rank)
    end
end

const reckoner_defaults = ReckonerInstance{DefaultMatch, DefaultMatches}()

# Redefine elemental functions in terms of a specific ReckonerInstance

function aup(curr::R, inst::ReckonerInstance{R,T} = reckoner_defaults)::T where {R, T}
    inst.aup(curr)
end

function weight(curr::R, prev::R, inst::ReckonerInstance{R,T} = reckoner_defaults)::Float64 where {R, T}
    inst.weight(curr, prev)
end

function challenge_window(curr::R, prev::R, inst::ReckonerInstance{R,T} = reckoner_defaults)::Float64 where {R, T}
    inst.challenge_window(curr, prev)
end

function skill(wins::Vector{Bool}, weights::Vector{<:Real}, challenge_windows::Vector{<:Real}, inst::ReckonerInstance{R,T} = reckoner_defaults)::Beta{Float64} where {R, T}
    inst.skill(wins, weights, challenge_windows)
end

function rating(wins::Vector{Bool}, weights::Vector{<:Real}, win_chances::Vector{<:Real}, inst::ReckonerInstance{R,T} = reckoner_defaults)::Beta{Float64} where {R, T}
    inst.skill(wins, weights, win_chances)
end

function eff_challenge(benchmarks::Vector{Beta{Float64}}, teams::Vector{<:Integer}, inst::ReckonerInstance{R,T} = reckoner_defaults)::Vector{Beta{Float64}} where {R, T}
    inst.eff_challenge(benchmarks, teams)
end

function win_chances(local_skills::Vector{Beta{Float64}}, teams::Vector{<:Integer}, inst::ReckonerInstance{R,T} = reckoner_defaults)::Vector{Beta{Float64}} where {R, T}
    inst.win_chances(local_skills, teams)
end

function display_rank(win_chance::Real, inst::ReckonerInstance{R,T} = reckoner_defaults)::Float64 where {R, T}
    inst.display_rank(win_chance)
end

# Define variations of the elemental functions with overlapping interfaces

function weights(curr::R, prev::Union{R, T}, inst::ReckonerInstance{R,T} = reckoner_defaults)::Vector{Float64} where {R, T}
    new_prev = Tables.rows(prev)
    inst.weight.((curr,), new_prev)
end

function weights(curr::Vector{R}, prev::Union{Vector{T}, Vector{R}}, inst::ReckonerInstance{R,T} = reckoner_defaults)::Vector{Vector{Float64}} where {R, T}
    weights.(curr, prev, (inst,))
end

function challenge_windows(curr::R, prev::Union{R, T}, inst::ReckonerInstance{R,T} = reckoner_defaults)::Vector{Float64} where {R, T}
    new_prev = Tables.rows(prev)
    inst.challenge_window.((curr,), new_prev)
end

function challenge_windows(curr::Vector{R}, prev::Union{Vector{T}, Vector{R}}, inst::ReckonerInstance{R,T} = reckoner_defaults)::Vector{Vector{Float64}} where {R, T}
    challenge_windows.(curr, prev, (inst,))
end

function skill(curr::R, prev::T, inst::ReckonerInstance{R,T} = reckoner_defaults)::Beta{Float64} where {R, T}
    calc_weights::Vector{Float64} = weights(curr, prev, inst)

    calc_windows::Vector{Float64} = challenge_windows(curr, prev, inst)

    inst.skill(win(prev), calc_weights, calc_windows)
end

function skills(curr::Vector{R}, prev::Vector{T}, inst::ReckonerInstance{R,T} = reckoner_defaults)::Vector{Beta{Float64}} where {R, T}
    skill.(curr, prev, (inst,))
end

function rating(curr::R, prev::T, inst::ReckonerInstance{R,T} = reckoner_defaults)::Beta{Float64} where {R, T}
    calc_weights::Vector{Float64} = weights(curr, prev, inst)

    inst.rating(win(prev), calc_weights, win_chance(prev))
end

function ratings(curr::Vector{R}, prev::Vector{T}, inst::ReckonerInstance{R,T} = reckoner_defaults)::Vector{Beta{Float64}} where {R, T}
    rating.(curr, prev, (inst,))
end

function eff_challenge(curr::Vector{R}, prev::Vector{T}, inst::ReckonerInstance{R,T} = reckoner_defaults)::Vector{Beta{Float64}} where {R,T}
    benches::Vector{Beta{Float64}} = ratings(curr, prev, inst)

    inst.eff_challenge(benches, team_id.(curr))
end

function player_win_chances(local_skills::Vector{Beta{Float64}}, teams::Vector{<:Integer}, inst::ReckonerInstance{R,T} = reckoner_defaults)::Vector{Float64} where {R, T}
    team_chances::Vector{Float64} = mean.(inst.win_chances(local_skills, teams))

    [team_chances[i] for i in teams]
end

function player_win_chances(curr::Vector{R}, prev::Vector{T}, inst::ReckonerInstance{R,T} = reckoner_defaults)::Vector{Float64} where {R, T}
    local_skills::Vector{Beta{Float64}} = skills(curr, prev, inst)

    player_win_chances(local_skills, team_id.(curr), inst)
end


text:
using DocumenterTools

GitHub_username = "Shoichiro-Tsutsui"
repository_name = "LearnTravisCI"

Travis.genkeys(
    user = GitHub_username,
    repo = "git@github.com:" * GitHub_username * "/" * repository_name *".git")

max_stars_repo_name: Shoichiro-Tsutsui/LearnTravisCI

max_stars_repo_licenses: FTL RSA-MD

max_issues_repo_name: Shoichiro-Tsutsui/LearnTravisCI

max_issues_repo_licenses: FTL RSA-MD

max_forks_repo_name: Shoichiro-Tsutsui/LearnTravisCI

max_forks_repo_licenses: FTL RSA-MD


text:
@testset "auxiliary functions" begin
    test = [reshape((1:18) .+ i,3,6) for i in 1:8];
    test_fl = LadderDGA.flatten_2D(test)
    res = true
    for i in 1:8
        !(test[i] == reshape(test_fl[i,:],3,6)) && (res = false)
    end
    @test res

    @testset "reduce_range" begin
        @test all(LadderDGA.reduce_range(1:3, 1.0) .== 2:2)
        @test all(LadderDGA.reduce_range(1:4, 1.0) .== 2:3)
        @test all(LadderDGA.reduce_range(1:5, 1.0) .== 3:3)
        @test all(LadderDGA.reduce_range(1:6, 1.0) .== 3:4)
        @test all(LadderDGA.reduce_range(1:10, 0.2) .== 2:9)
        @test all(LadderDGA.reduce_range(1:20, 0.1) .== 2:19)
        @test all(LadderDGA.reduce_range(1:21, 0.1) .== 2:20)
        @test all(LadderDGA.reduce_range(1:20, 0.5) .== 6:15)
        @test all(LadderDGA.reduce_range(1:21, 0.5) .== 6:16)
    end
end


text:
using ShallowWaters
using Test

@testset "No Forcing" begin
    Prog = run_model(Fx0=0)
    @test all(Prog.u .== 0.0f0)
    @test all(Prog.v .== 0.0f0)
    @test all(Prog.η .== 0.0f0)
end

@testset "Boundary Conditions" begin
    Prog = run_model(bc="periodic")
    @test all(abs.(Prog.u) .< 10)    # velocity shouldn't exceed +-10m/s

    Prog = run_model(bc="nonperiodic")
    @test all(abs.(Prog.u) .< 10)

    Prog = run_model(α=0)     # free-slip
    @test all(abs.(Prog.u) .< 10)
end

@testset "Linear dynamics" begin
    P = run_model(dynamics="linear")
    @test all(abs.(P.u) .< 10)   # velocity shouldn't exceed +-10m/s
end

@testset "Continuity Forcing" begin
    Prog = run_model(surface_relax=true)
    @test all(abs.(Prog.u) .< 10)    # velocity shouldn't exceed +-10m/s
    @test all(Prog.u .!= 0.0f0)

    Prog = run_model(surface_forcing=true)
    @test all(abs.(Prog.u) .< 10)
    @test all(Prog.u .!= 0.0f0)
end

@testset "Mixed Precision" begin
    Prog = run_model(Float32,Tprog=Float64)
    @test all(abs.(Prog.u) .< 10)    # velocity shouldn't exceed +-10m/s
    @test all(Prog.u .!= 0.0f0)

    Prog = run_model(Float16,Tprog=Float32,Ndays=50)
    @test all(abs.(Prog.u) .< 10)    # velocity shouldn't exceed +-10m/s
    @test all(Prog.u .!= 0.0f0)
end

@testset "RK3 tests" begin
    Prog = run_model(time_scheme="RK",RKo=3,cfl=0.6)
    @test all(abs.(Prog.u) .< 10)    # velocity shouldn't exceed +-10m/s
    @test all(Prog.u .!= 0.0f0)
end

@testset "SSPRK2 tests" begin
    Ndays = 500
    Prog = run_model(time_scheme="SSPRK2",RKs=2,cfl=0.7;Ndays)
    @test all(abs.(Prog.u) .< 10)    # velocity shouldn't exceed +-10m/s
    @test all(Prog.u .!= 0.0f0)

    Prog = run_model(time_scheme="SSPRK2",RKs=3,cfl=1.2;Ndays)
    @test all(abs.(Prog.u) .< 10)    # velocity shouldn't exceed +-10m/s
    @test all(Prog.u .!= 0.0f0)

    Prog = run_model(time_scheme="SSPRK2",RKs=4,cfl=1.3;Ndays)
    @test all(abs.(Prog.u) .< 10)    # velocity shouldn't exceed +-10m/s
    @test all(Prog.u .!= 0.0f0)
end

@testset "SSPRK3 tests" begin
    Ndays = 500
    Prog = run_model(time_scheme="SSPRK3",RKn=2,cfl=1.2;Ndays)
    @test all(abs.(Prog.u) .< 10)    # velocity shouldn't exceed +-10m/s
    @test all(Prog.u .!= 0.0f0)

    Prog = run_model(time_scheme="SSPRK3",RKn=3,cfl=2.5;Ndays)
    @test all(abs.(Prog.u) .< 10)    # velocity shouldn't exceed +-10m/s
    @test all(Prog.u .!= 0.0f0)
end

@testset "Scaling" begin
    P64 = run_model(nx=100,scale=64);
    P1 = run_model(nx=100,scale=1);

    @test P64.u == P1.u
    @test P64.v == P1.v
    @test P64.η == P1.η
    @test P64.sst == P1.sst
end


text:
using UNet
using Flux
using Zygote
using MAT
using Statistics
using LinearAlgebra
using Random
using BSON
using JLD
using PyPlot

#=
# load variables
mf = matopen("/home/naxos2-raid1/maxiao/Projects/pTXDL/Data_raw.mat")
Xtrain_AP_1 = read(mf, "Xtrain_AP_1")
Xtrain_AP_2 = read(mf, "Xtrain_AP_2")
Ytrain_AP_1 = read(mf, "Ytrain_AP_1_rot")
Ytrain_AP_2 = read(mf, "Ytrain_AP_2_rot")
close(mf)

#=
fig, ax = plt.subplots()
ax.imshow(abs.(Xtrain_AP_1[:,:,50,1,1]))
display(fig)
ax.imshow(abs.(Ytrain_AP_1[:,:,50,1,1]))
display(fig)
=#
Xtrain_AP_1_subj1 = Xtrain_AP_1[:,:,:,:,1]./maximum(maximum(maximum(Xtrain_AP_1[:,:,:,1,1])))
Xtrain_AP_1_subj2 = Xtrain_AP_1[:,:,:,:,2]./maximum(maximum(maximum(Xtrain_AP_1[:,:,:,1,2])))
Xtrain_AP_1_subj3 = Xtrain_AP_1[:,:,:,:,3]./maximum(maximum(maximum(Xtrain_AP_1[:,:,:,1,3])))
Xtrain_AP_2_subj1 = Xtrain_AP_2[:,:,:,:,1]./maximum(maximum(maximum(Xtrain_AP_1[:,:,:,1,1])))
Xtrain_AP_2_subj2 = Xtrain_AP_2[:,:,:,:,2]./maximum(maximum(maximum(Xtrain_AP_1[:,:,:,1,2])))
Xtrain_AP_2_subj3 = Xtrain_AP_2[:,:,:,:,3]./maximum(maximum(maximum(Xtrain_AP_1[:,:,:,1,3])))
Xtrain_AP_1 = nothing
Xtrain_AP_2 = nothing

Ytrain_AP_1_subj1 = Ytrain_AP_1[:,:,:,:,1]./maximum(maximum(maximum(Ytrain_AP_1[:,:,:,1,1])))
Ytrain_AP_1_subj2 = Ytrain_AP_1[:,:,:,:,2]./maximum(maximum(maximum(Ytrain_AP_1[:,:,:,1,2])))
Ytrain_AP_1_subj3 = Ytrain_AP_1[:,:,:,:,3]./maximum(maximum(maximum(Ytrain_AP_1[:,:,:,1,3])))
Ytrain_AP_2_subj1 = Ytrain_AP_2[:,:,:,:,1]./maximum(maximum(maximum(Ytrain_AP_1[:,:,:,1,1])))
Ytrain_AP_2_subj2 = Ytrain_AP_2[:,:,:,:,2]./maximum(maximum(maximum(Ytrain_AP_1[:,:,:,1,2])))
Ytrain_AP_2_subj3 = Ytrain_AP_2[:,:,:,:,3]./maximum(maximum(maximum(Ytrain_AP_1[:,:,:,1,3])))
Ytrain_AP_1 = nothing
Ytrain_AP_2 = nothing

#Xtrain = cat(Xtrain_AP_1_subj1,Xtrain_AP_2_subj1,Xtrain_AP_1_subj2,Xtrain_AP_2_subj2,dims=4)
Xtrain = cat(Xtrain_AP_1_subj1,Xtrain_AP_2_subj1,dims=4)
Sall = size(Xtrain)
Xtrain = reshape(Xtrain,Sall[1],Sall[2],1,Sall[3]*Sall[4])

#Ytrain = cat(Ytrain_AP_1_subj1,Ytrain_AP_2_subj1,Ytrain_AP_1_subj2,Ytrain_AP_2_subj2,dims=4)
Ytrain = cat(Ytrain_AP_1_subj1,Ytrain_AP_2_subj1,dims=4)
Sall = size(Ytrain)
Ytrain = reshape(Ytrain,Sall[1],Sall[2],1,Sall[3]*Sall[4])

XValid = cat(Xtrain_AP_1_subj3,Xtrain_AP_2_subj3,dims=4)
Sall = size(XValid)
XValid = reshape(XValid,Sall[1],Sall[2],1,Sall[3]*Sall[4])

YValid = cat(Ytrain_AP_1_subj3,Ytrain_AP_2_subj3,dims=4)
Sall = size(YValid)
YValid = reshape(YValid,Sall[1],Sall[2],1,Sall[3]*Sall[4])

Xtrain_AP_1_subj1 = nothing;
Xtrain_AP_1_subj2 = nothing;
Xtrain_AP_1_subj3 = nothing;
Xtrain_AP_2_subj1 = nothing;
Xtrain_AP_2_subj2 = nothing;
Xtrain_AP_2_subj3 = nothing;
Ytrain_AP_1_subj1 = nothing;
Ytrain_AP_1_subj2 = nothing;
Ytrain_AP_1_subj3 = nothing;
Ytrain_AP_2_subj1 = nothing;
Ytrain_AP_2_subj2 = nothing;
Ytrain_AP_2_subj3 = nothing;

idx_shuffle = shuffle(collect(1:size(Xtrain,4)))
Xtrain = Xtrain[:,:,:,idx_shuffle]
Ytrain = Ytrain[:,:,:,idx_shuffle]

JLD.save("XYtrain_subj1.jld","Xtrain",Xtrain,"Ytrain",Ytrain,"XValid",XValid,"YValid",YValid)
=#

Xtrain = JLD.load("XYtrain_subj1.jld","Xtrain")
Ytrain = JLD.load("XYtrain_subj1.jld","Ytrain")

train_loader = Flux.Data.DataLoader(Xtrain,Ytrain, batchsize=25,shuffle=true)

u = Unet(1)

loss(x,y) = Flux.mse(u(x),y)
validationloss() = Flux.mse(u(XValid),YValid)

cb = Flux.throttle(() -> @show(validationloss()), 10)
#@time Flux.train!(loss, θ, dataset, opt, cb = cb)
opt = ADAM(0.001)

function my_custom_train!(loss, ps, data, opt)
    # training_loss is declared local so it will be available for logging outside the gradient calculation.
    local training_loss
    ps = Params(ps)
    for d in data
      if d isa AbstractArray{<:Number}
        gs = Flux.gradient(ps) do
          training_loss = loss(d)
          return training_loss
        end
      else
        gs = Flux.gradient(ps) do
          training_loss = loss(d...)
          return training_loss
        end
      end
      # Insert whatever code you want here that needs training_loss, e.g. logging.
      # logging_callback(training_loss)
      # Insert what ever code you want here that needs gradient.
      # E.g. logging with TensorBoardLogger.jl as histogram so you can see if it is becoming huge.
      #println("training_loss = ",training_loss,"; gradient[1] = ",gs[ps[1]],"; W1 = ",ps.order.data[1][1],"; b1 = ",ps.order.data[2][1])
      println("training_loss = ",training_loss)
      Flux.update!(opt, ps, gs)
      # Here you might like to check validation set accuracy, and break out to do early stopping.
    end
  end

for epoch in 1:1
#@time Flux.train!(loss, Flux.params(u), train_loader, opt, cb = cb)
@time my_custom_train!(loss, Flux.params(u), train_loader, opt)
end


text:
export rv2kep,
        get_hyperbolic_tof,
        get_hyperbolic_anomaly,
        get_soa,
        get_soi,
        getv,
        getvhyp,
        mjd2jd,
        jd2mjd,
        get_approach_phase_angle,
        get_HORIZONS_ephem,
        get_sse

"""
Transform a Cartesian ECI state vector to a set of
classical/Keplerian elements.

x_kep - [a,e,i,RAAN,AOP,nu], where:
    a    = semimajor axis [km]
    e    = eccentricity
    i    = inclination [rad]
    RAAN = right ascension of the ascending node [rad]
    AOP  = argument of perigee [rad]
    nu   = true anomaly [rad]

Reference: GTDS math spec

**Undefined quantities return as 0.0**
"""
function rv2kep(x_rv::AbstractVector{T}, mu = muEarth) where {T<:AbstractFloat}
    r = norm(x_rv[1:3])
    v = norm(x_rv[4:6])
    h = cross(x_rv[1:3],x_rv[4:6])
    n = cross([0;0.;1],h)

    sma = mu * r/(2*mu - r*v^2)
    evec = eccvec(x_rv)
    ecc = norm(evec)
    inc = angle_between(h, [0.0; 0.0; 1.0])

    if norm(n) < eps() # equatorial
        RAAN = 0.
    else
        RAAN = angle_between([1.0; 0.0; 0.0], n)
        if n[2] < 0.
            RAAN = 2pi - RAAN
        end
    end

    # circular
    if ecc < eps()
        AOP = 0.
        nu = 0.
    else
        if norm(n) < eps() # equatorial
            AOP = 0.0
        else
            AOP = angle_between(n, evec)
            if evec[3]<0.
                AOP = 2pi - AOP
            end
        end
        nu = angle_between(evec, x_rv[1:3])
        if dot(x_rv[1:3],x_rv[4:6]) < 0.
            nu = 2pi - nu
        end
    end
    return [sma;ecc;inc;RAAN;AOP;nu]
end

function get_hyperbolic_anomaly(r, a_hyp, e_hyp)
    return acosh((1. - r/a_hyp)/e_hyp)
end

"""
Return the TOF in seconds given the SMA and eccentricity of a
hyperbola and r0 and r1
"""
function get_hyperbolic_tof(r0, r1, a_hyp, e_hyp)
    H0 = get_hyperbolic_anomaly(r0,a_hyp,e_hyp)
    H1 = get_hyperbolic_anomaly(r1,a_hyp,e_hyp)

    return sqrt(-a_hyp^3/muEarth) *
            (e_hyp * sinh(H1) - H1 - (e_hyp * sinh(H0) - H0))
end

"""
Calculate the Sphere of Influence radius for a body, given it's
properties and the properties of the larger body
"""
function get_soi(m_smaller, M_larger, R_CB)
    return R_CB*(m_smaller/M_larger)^(2/5.)
end

"""
Calculate the Sphere of Activity radius for a body, given it's
properties and the properties of the larger body
"""
function get_soa( m_smaller, M_larger, R_CB )
    return R_CB*(m_smaller/(3.0 * M_larger))^(1/3)
end

"Calculate elliptical orbit velocity"
function getv(mu, r, a)
    return sqrt(mu*((2/r) - 1/a))
end

"Calculate hyperbolic orbit velocity"
function getvhyp(mu, r, vInf)
    return sqrt(vInf^2 + 2*mu/r)
end

function getperiod(mu, a)
    return 2pi*sqrt(a^3/mu)
end

"Convert Julian date to Modified Julian date"
function jd2mjd(jd::Float64)
    return jd - 2400000.5
end

"Convert Modified Julian date to Julian date"
function mjd2jd(mjd::Float64)
    return mjd + 2400000.5
end

"Calculate the approach phase angle (Sun - x_target - x)"
function get_approach_phase_angle(x::Vector{T}, x_target::Vector{T}) where {T<:AbstractFloat}
    v_rel = (x[4:6] - x_target[4:6])./norm(x[4:6] - x_target[4:6])
    return acos( dot(v_rel, x_target[1:3]./norm(x_target[1:3]) ) )
end

"Calculate the Sun-Spacecraft-Earth angle"
function get_sse(r_SC_S::Vector{T}, r_E_S::Vector{T}) where {T<:AbstractFloat}
    r_S_SC = -r_SC_S
    r_E_SC = r_E_S - r_SC_S
    return acos( dot(r_S_SC,r_E_SC)/(norm(r_S_SC)*norm(r_E_SC)) )
end

"Parse a HORIZONS ephemeris file"
function get_HORIZONS_ephem(f::IOStream)
    ephem = Vector{Vector{Float64}}(undef,0)
    regex_jdate = r"\d{7}.\d{9}"
    #regex_datetime = r"\d{4}-[a-zA-Z]{3}-\d{2}" # not used
    #regex_time = r"\d{2}:\d{2}:\d{2}.\d+"       # not used
    regex_state = r"[-+]?\d+?.\d+[eE][+-]\d+.*$" #matches from the first match to the end of the line
    while !eof(f)
        x = readline(f)
        if contains(x,"\$\$SOE") # start of ephem
            x = readline(f)
            while !contains(x,"\$\$EOE") #end of ephemeris
                eph_current = float64(split(match(regex_state,x).match,',')[1:6])' #extra comma is off the end
                jdate_current = float64(match(regex_jdate,x).match)
                push!(ephem, vcat(jdate_current,eph_current))
                x= readline(f)
            end
        end
    end
    return hcat(ephem...)'
end


text:
#
# Minimize the acceleration of a composite Bézier curve on the Sphere S2 with interpolation
#
# This example appeared in Sec. 5.2, second example, of
#
# > Bergmann, R. and Gousenbourger, P.-Y.: _A variational model for data fitting on manifolds
# > by minimizing the acceleration of a Bézier curve_.
# > Frontiers in Applied Mathematics and Statistics, 2018.
# > doi: [10.3389/fams.2018.00059](https://dx.doi.org/10.3389/fams.2018.00059),
# > arXiv: [1807.10090](https://arxiv.org/abs/1807.10090)
#
using Manopt, Manifolds, Colors, ColorSchemes
#
# Load Manopt and required packages
#
using Manopt, Manifolds, Colors, ColorSchemes
import Printf.@sprintf
import ColorSchemes.viridis

results_folder = joinpath(@__DIR__, "Minimize_Acceleration_Bezier", "video")
!isdir(results_folder) && mkdir(results_folder)
final_folder = joinpath(@__DIR__, "Minimize_Acceleration_Bezier")
!isdir(final_folder) && mkdir(final_folder)
experiment_name = "Bezier_Approximation_video"

render_detail = 2
asy_export = true
asy_export_summary = true
render_video = true
λRange = collect(range(10.0, 0.0; length=1001))[1:(end - 1)] #exclude zero.
colors = RGBA.(get.(Ref(viridis), range(0.0, 1.0; length=length(λRange))))

curve_samples = [range(0, 3; length=101)...] # sample curve for the gradient
curve_samples_plot = [range(0, 3; length=201)...] # sample curve for asy exports

camera_position = (-1.0, -0.7, 0.3)
curveColor = RGBA{Float64}(colorant"#AAAAAA")
M = Sphere(2)
B = artificial_S2_composite_bezier_curve()
degs = get_bezier_degrees(M, B)
cP = de_casteljau(M, B, curve_samples_plot)
d = get_bezier_junctions(M, B)
pB = get_bezier_points(M, B, :differentiable)
results = [similar(pB) for i in eachindex(λRange)]
resulting_curves = [similar(cP) for i in eachindex(λRange)]
N = PowerManifold(M, NestedPowerRepresentation(), length(pB))
x0 = pB
for i in eachindex(λRange)
    λ = λRange[i]
    F(pB) = cost_L2_acceleration_bezier(M, pB, degs, curve_samples, λ, d)
    ∇F(pB) = ∇L2_acceleration_bezier(M, pB, degs, curve_samples, λ, d)
    results[i] = gradient_descent(
        N,
        F,
        ∇F,
        x0;
        stepsize=ArmijoLinesearch(1.0, ExponentialRetraction(), 0.5, 0.001),
        stopping_criterion=StopWhenAny(
            StopWhenChangeLess(5 * 10.0^(-6)), StopAfterIteration(15000)
        ),
        debug=[
            :Stop,
            "$(λ) ",
            :Iteration,
            " | ",
            :Cost,
            " | ",
            DebugGradientNorm(),
            " | ",
            DebugStepsize(),
            " | ",
            :Change,
            "\n",
            200,
        ],
    )
    B_opt = get_bezier_segments(M, results[i], degs, :differentiable)
    resulting_curves[i] = de_casteljau(M, B_opt, curve_samples_plot)
end
if asy_export
    print("Exporting...\n")
    s = joinpath(results_folder, experiment_name)
    for i in eachindex(λRange)
        print("$(λRange[i])..")
        B_opt = get_bezier_segments(M, results[i], degs, :differentiable)
        asymptote_export_S2_signals(
            s * "-$(@sprintf "%04.0f" i)-result.asy";
            curves=[resulting_curves[i], cP],
            points=[get_bezier_junctions(M, B_opt), get_bezier_inner_points(M, B_opt)],
            colors=Dict(
                :curves => [colors[i], curveColor], :points => [colors[i], colors[i]]
            ),
            camera_position=camera_position,
            line_widths=[1.0, 0.5],
            dot_size=2.0,
        )
        render_asymptote(s * "-$(@sprintf "%04.0f" i)-result.asy"; render=render_detail)
    end
end
if asy_export_summary
    s = joinpath(final_folder, experiment_name)
    asymptote_export_S2_signals(
        s * "-Summary-result.asy";
        curves=[cP, resulting_curves...],
        colors=Dict(:curves => [curveColor, colors...]),
        camera_position=camera_position,
        line_widths=[0.75, [1.5 for i in eachindex(λRange)]...],
    )
    render_asymptote(s * "-Summary-result.asy"; render=render_detail)
end
if render_video
    s = joinpath(results_folder, experiment_name)
    s2 = joinpath(final_folder, experiment_name)
    cmd = `ffmpeg -i $(s)-%04d-result.png -r 25 -framerate 60 -c:v libx264 -crf 20 -pix_fmt yuv420p -y $(s2)-movie.mp4`
    run(cmd)
end


text:
export
    SolverOptions

abstract type AbstractSolverOptions{T<:Real} end
abstract type DirectSolverOptions{T} <: AbstractSolverOptions{T} end

"""
    set_options!(opts::AbstractSolverOptions; kwargs...)
    set_options!(solver::AbstractSolver; opts...)

Set solver options via keyword arguments, supporting either solvers or solver option
types directly. Will set any and all options that match the provided arguments. For
example, `set_options!(solver::ALTROSolver, constraint_tolerance=1e-4)` will set the
constraint tolerance option in the ALTRO, Augmented Lagrangian, and Project Newton
solvers.

The only exeption is the `verbose` setting, which always accepts a boolean, while
ALTRO and Augmented Lagrangian solvers accept integers 0-2, with 1 providing output
for the outer AL iterations but not the iLQR iterations.
"""
function set_options!(opts::OPT; d...) where OPT <: AbstractSolverOptions
    for (key,val) in pairs(d)
        if hasfield(OPT, key)
            val = convert(fieldtype(OPT, key), val)
            setfield!(opts, key, val)
        end
    end
end

@with_kw mutable struct SolverOptions{T} <: AbstractSolverOptions{T}
    # Optimality Tolerances
    constraint_tolerance::T = 1e-6
    cost_tolerance::T = 1e-4
    cost_tolerance_intermediate::T = 1e-4
    gradient_tolerance::T = 10.0
    gradient_tolerance_intermediate::T = 1.0

    # iLQR
    iterations_inner::Int = 300
    dJ_counter_limit::Int = 10
    square_root::Bool = false
    line_search_lower_bound::T = 1e-8
    line_search_upper_bound::T = 10.0
    iterations_linesearch::Int = 20
    max_cost_value::T = 1.0e8
    max_state_value::T = 1.0e8
    max_control_value::T = 1.0e8
    static_bp::Bool = true
	save_S::Bool = false

    # Backward pass regularization
    bp_reg::Bool = false
    bp_reg_initial::T = 0.0
    bp_reg_increase_factor::T = 1.6
    bp_reg_max::T = 1.0e8
    bp_reg_min::T = 1.0e-8
    bp_reg_type::Symbol = :control
    bp_reg_fp::T = 10.0

    # Augmented Lagrangian
    penalty_initial::T = 1.0
    penalty_scaling::T = 10.0
    active_set_tolerance_al::T = 1e-3
    dual_max::T = 1e8
    penalty_max::T = 1e8
    iterations_outer::Int = 30
    kickout_max_penalty::Bool = false
    reset_duals::Bool = true
    reset_penalties::Bool = true

    # Projected Newton
    force_pn::Bool = false
    verbose_pn::Bool = false
    n_steps::Int = 2
    solve_type::Symbol = :feasible
    projected_newton_tolerance::T = 1e-3
    active_set_tolerance_pn::T = 1e-3
    multiplier_projection::Bool = true
    ρ_chol::T = 1e-2     # cholesky factorization regularization
    ρ_primal::T = 1.0e-8 # primal regularization
    ρ_dual::T = 1.0e-8   # regularization for multiplier projection
    r_threshold::T = 1.1

    # General options
    projected_newton::Bool = true
    reuse_jacobians::Bool = false
    iterations::Int = 1000   # max number of iterations
    show_summary::Bool = true
    verbose::Int = 0
end

function Base.copy(opts::SolverOptions)
    SolverOptions([getfield(opts, fname) for fname in fieldnames(SolverOptions)]...)
end

@with_kw mutable struct SolverStats{T}
    # Iteration counts
    iterations::Int = 0
    iterations_outer::Int = 0
    iterations_pn::Int = 0

    # Iteration stats
    iteration::Vector{Int} = Int[]
    iteration_outer::Vector{Int} = Int[]
    iteration_pn::Vector{Bool} = Bool[]   # is the iteration a pn iteration
    cost::Vector{T} = Float64[]
    dJ::Vector{T} = Float64[]
    c_max::Vector{T} = Float64[]
    gradient::Vector{T} = Float64[]
    penalty_max::Vector{T} = Float64[]

    # iLQR
    dJ_zero_counter::Int = 0

    # Other
    tstart::Float64 = time()
    tsolve::Float64 = Inf
    to::TimerOutput = TimerOutput()
    status::TerminationStatus = UNSOLVED
    is_reset::Bool = false
    "Which solver is the top-level solver and responsible for resetting and trimming."
    parent::Symbol
end

function reset!(stats::SolverStats, N::Int, parent::Symbol)
    if parent == stats.parent
        stats.is_reset = false  # force a reset
        reset!(stats, N)
    end
    return nothing
end

function reset!(stats::SolverStats, N::Int=0)
    stats.is_reset && return nothing
    stats.iterations = 0
    stats.iterations_outer = 0
    stats.iterations_pn = 0
    function reset!(v::AbstractVector{T}, N::Int) where T
        resize!(v, N)
        v .= zero(T)
    end
    reset!(stats.iteration, N)
    reset!(stats.iteration_outer, N)
    reset!(stats.iteration_pn, N)
    reset!(stats.cost, N)
    reset!(stats.dJ, N)
    reset!(stats.c_max, N)
    reset!(stats.gradient, N)
    reset!(stats.penalty_max, N)
    stats.tstart = time_ns()
    stats.tsolve = Inf
    # reset_timer!(stats.to)
    stats.dJ_zero_counter = 0
    stats.is_reset = true
    stats.status = UNSOLVED
    return nothing
end

function record_iteration!(stats::SolverStats;
        cost=NaN, dJ=NaN, c_max=NaN, gradient=NaN, penalty_max=NaN,
        is_pn::Bool=false, is_outer::Bool=false
    )
    # don't increment the iteration for an outer loop
    if is_outer
        stats.iterations_outer += 1
    else
        stats.iterations += 1
    end
    i = stats.iterations
    function record!(vec, val, i)
        if isnan(val)
            if vec[i] != 0
                val = vec[i]
            elseif i > 1
                val = vec[i-1]
            else
                val = vec[i]
            end
        end
        vec[i] = val
    end
    stats.iteration[i] = i
    stats.iteration_outer[i] = stats.iterations_outer
    stats.iteration_pn[i] = is_pn
    record!(stats.cost, cost, i)
    record!(stats.dJ, dJ, i)
    record!(stats.c_max, c_max, i)
    record!(stats.gradient, gradient, i)
    record!(stats.penalty_max, penalty_max, i)
    is_pn && (stats.iterations_pn += 1)
    return stats.iterations
end

function trim!(stats::SolverStats, parent::Symbol)
    if parent == stats.parent
        trim!(stats)
    end
end
function trim!(stats::SolverStats)
    N = stats.iterations
    resize!(stats.iteration, N)
    resize!(stats.iteration_outer, N)
    resize!(stats.iteration_pn, N)
    resize!(stats.cost, N)
    resize!(stats.dJ, N)
    resize!(stats.c_max, N)
    resize!(stats.gradient, N)
    resize!(stats.penalty_max, N)
    stats.is_reset = false
    return N
end

function Dict(stats::SolverStats)
    fields = (:iteration, :iteration_outer, :iteration_pn,
        :cost, :dJ, :c_max, :gradient, :penalty_max)
    pairs = [field=>getfield(stats,field) for field in fields]
    Dict(pairs...)
end


text:
"""
    poisson_eigenvalues(N, L, dim, ::Periodic)

Return the eigenvalues satisfying the discrete form of Poisson's equation
with periodic boundary conditions along the dimension `dim` with `N` grid
points and domain extent `L`.
"""
function poisson_eigenvalues(N, L, dim, ::Periodic)
    inds = reshape(1:N, reshaped_size(N, dim)...)
    return @. (2sin((inds - 1) * π / N) / (L / N))^2
end

"""
    poisson_eigenvalues(N, L, dim, ::Bounded)

Return the eigenvalues satisfying the discrete form of Poisson's equation
with staggered Neumann boundary conditions along the dimension `dim` with
`N` grid points and domain extent `L`.
"""
function poisson_eigenvalues(N, L, dim, ::Bounded)
    inds = reshape(1:N, reshaped_size(N, dim)...)
    return @. (2sin((inds - 1) * π / 2N) / (L / N))^2
end

"""
    poisson_eigenvalues(N, L, dim, ::Flat)

Return N-element array of `0.0` reshaped to three-dimensions.
This is also the first `poisson_eigenvalue` for `Bounded` and `Periodic` directions.
"""
poisson_eigenvalues(N, L, dim, ::Flat) = reshape(zeros(N), reshaped_size(N, dim)...)