text:
function fastnorm(u)
  # dest[1] = ...
  n = length(u)
  T = eltype(u)
  s = zero(T)
  @fastmath @inbounds @simd for i in 1:n
      s += u[i]^2
  end
  @fastmath @inbounds return sqrt(s)
end

"""
    $TYPEDSIGNATURES

Calculate the Kernel Embedding MMD 'distance' between sample points (or kernel density estimates).

Notes
- `bw::Vector=[0.001;]` controls the mmd kernel bandwidths.
- Overloading from ApproxManifoldProducts

Related

`AMP.kld`
"""
function mmd( p1::AbstractVector{P1},
              p2::AbstractVector{P2},
              varType::Union{InstanceType{<:InferenceVariable},InstanceType{<:AbstractFactor}};
              bw::AbstractVector{<:Real}=[0.001;] ) where {P1 <: AbstractVector, P2 <: AbstractVector}
  #
  mani = getManifold(varType)
  mmd(p1, p2, mani, bw=bw)
end

# TODO move to AMP?
function mmd( p1::ManifoldKernelDensity,
              p2::ManifoldKernelDensity,
              nodeType::Union{InstanceType{<:InferenceVariable},InstanceType{<:AbstractFactor}};
              bw::AbstractVector{<:Real}=[0.001;])
  #
  mmd(getPoints(p1), getPoints(p2), nodeType, bw=bw)
end

# moved to CalcFactor.jl

function sampleFactor(ccwl::CommonConvWrapper,
                      N::Int  )
  #
  cf = CalcFactor( ccwl.usrfnc!, _getFMdThread(ccwl), 0, length(ccwl.measurement), ccwl.measurement, ccwl.params)
  sampleFactor(cf, N)
end

# part of consolidation, see #927
function sampleFactor!( ccwl::CommonConvWrapper,
                        N::Int,
                        fmd::FactorMetadata=_getFMdThread(ccwl),
                        vnd=nothing )
  #
  # depr warning added before IIF v0.20
  vnd !== nothing ? @warn("sampleFactor! no longer accepts vnd::Vector as meaningful input.") : nothing

  # build a CalcFactor object and get fresh samples.
  cf = CalcFactor( ccwl.usrfnc!, fmd, 0, length(ccwl.measurement), ccwl.measurement, ccwl.params)
  # TODO make this an in-place operation as far possible
  ccwl.measurement = sampleFactor(cf, N)

  nothing
end

sampleFactor(fct::DFGFactor, N::Int=1) = sampleFactor(_getCCW(fct), N)

function sampleFactor(dfg::AbstractDFG,
                      sym::Symbol,
                      N::Int=1 )
  #
  return sampleFactor(getFactor(dfg, sym), N)
end

"""
    $(SIGNATURES)

Update cliq `cliqID` in Bayes (Juction) tree `bt` according to contents of `urt` -- intended use is to update main clique after a upward belief propagation computation has been completed per clique.
"""
function updateFGBT!( fg::AbstractDFG,
                      cliq::TreeClique,
                      IDvals::Dict{Symbol, TreeBelief};
                      dbg::Bool=false,
                      fillcolor::String="",
                      logger=ConsoleLogger()  )
  #
  # if dbg
  #   # TODO find better location for the debug information (this is old code)
  #   cliq.attributes["debug"] = deepcopy(urt.dbgUp)
  # end
  if fillcolor != ""
    setCliqueDrawColor!(cliq, fillcolor)
  end
  for (id,dat) in IDvals
    with_logger(logger) do
      @info "updateFGBT! up -- update $id, infoPerCoord=$(dat.infoPerCoord)"
    end
    updvert = DFG.getVariable(fg, id)
    setValKDE!(updvert, deepcopy(dat), true) ## TODO -- not sure if deepcopy is required
  end
  with_logger(logger) do
    @info "updateFGBT! up -- updated $(getLabel(cliq))"
  end
  nothing
end

"""
    $SIGNATURES

Build a graph given one factor and an ordered vector of `(variables types,nothing).  In addition, init values can be passed instead of nothing.

Notes
- Often used to quickly generate temporary graphs for a variety of local calculations.
- does not yet support split `_` characters in auto-find `lastVar` from `varPattern`.
- Will always add a factor, but will skip adding variable labels that already exist in `dfg`.

DevNotes
- TODO allow pts to be full MKD beliefs, part of replacing old `approxConvCircular`, see #1351
"""
function _buildGraphByFactorAndTypes!(fct::AbstractFactor,
                                      varTypes::Tuple,
                                      pts::Tuple=();
                                      dfg::AbstractDFG = initfg(),
                                      solveKey::Symbol=:default,
                                      newFactor::Bool=true,
                                      destPattern::Regex = r"x\d+",
                                      destPrefix::Symbol = match(r"[a-zA-Z_]+", destPattern.pattern).match |> Symbol,
                                      _allVars::AbstractVector{Symbol} = sortDFG(ls(dfg, destPattern)),
                                      currLabel::Symbol = 0 < length(_allVars) ? _allVars[end] : Symbol(destPrefix, 0),
                                      currNumber::Integer = reverse(match(r"\d+", reverse(string(currLabel))).match) |> x->parse(Int,x),
                                      graphinit::Bool = false,
                                      _blockRecursion::Bool=false  )
  #

  # TODO generalize beyond binary
  len = length(varTypes)
  vars = Symbol[Symbol(destPrefix, s_) for s_ in (currNumber .+ (1:len))]
  for (s_, vTyp) in enumerate(varTypes)
    # add the necessary variables
    exists(dfg, vars[s_]) ? nothing : addVariable!(dfg, vars[s_],  vTyp)
    # set the numerical values if available
    # TODO allow pts to come in as full MKD beliefs, not just one point
    ((0 < length(pts)) && (pts[s_] isa Nothing)) ? nothing : initManual!(dfg,  vars[s_], [pts[s_],], solveKey, bw=ones(getDimension(vTyp)))
  end
  # if newFactor then add the factor on vars, else assume only one existing factor between vars
  _dfgfct = newFactor ? addFactor!(dfg, vars, fct, graphinit=graphinit, _blockRecursion=_blockRecursion) : getFactor(dfg, intersect((ls.(dfg, vars))...)[1] )

  return dfg, _dfgfct
end

"""
    $SIGNATURES

Check if a variable might already be located at the test location, by means of a (default) `refKey=:simulated` PPE stored in the existing variables.

Notes
- Checks, using provided `factor` from `srcLabel` in `fg` to an assumed `dest` variable whcih may or may not yet exist.
- This function was written to aid in building simulation code,
  - it's use in real world usage may have unexpected behaviour -- hence not exported.
- Return `::Tuple{Bool, Vector{Float64}, Symbol}`, eg. already exists `(true, [refVal], :l17)`, or if a refernce variable does not yet `(false, [refVal], :l28)`.
  - Vector contains the PPE reference location of the new variable as calculated.
- Auto `destPrefix` is trying to parse `destRegex` labels like `l\\d+` or `tag\\d+`, won't work with weirder labels e.g. `:l_4_23`.
  - User can overcome weird names by self defining `destPrefix` and `srcNumber`.
  - User can also ignore and replace the generated new label `Symbol(destPrefix, srcNumber)`.
- This function does not add new variables or factors to `fg`, user must do that themselves after.
  - Useful to use in combination with `setPPE!` on new variable.
- At time of writing `accumulateFactorMeans` could only incorporate priors or binary relative factors.
  - internal info, see [`solveFactorParameteric`](@ref),
  - This means at time of writing `factor` must be a binary factor.
- Tip, if simulations are inducing odometry bias, think of using two factors from caller (e.g. simPerfect and simBias).

Example
```julia
# fg has :x5 and :l2 and PPEs :simulated exists in all variables
# user wants to add a factor from :x5 to potential new :l5, but maybe a (simulated) variable, say :l2, is already there.

newFactor = RoME.Pose2Point2BearingRange(Normal(), Normal(20,0.5))
isAlready, simPPE, genLabel = IIF._checkVariableByReference(fg, :x5, r"l\\d+", Point2, newFactor)

# maybe add new variable
if !isAlready
  @info "New variable with simPPE" genLabel simPPE
  newVar = addVariable!(fg, genLabel, Point2)
  addFactor!(fg, [:x5; genLabel], newFactor)

  # also set :simulated PPE for similar future usage
  newPPE = DFG.MeanMaxPPE(:simulated, simPPE, simPPE, simPPE)
  setPPE!(newVar, :simulated, typeof(newPPE), newPPE)   # TODO this API can be improved
else
  @info "Adding simulated loop closure with perfect data association" :x5 genLabel
  addFactor!(fg, [:x5; genLabel], newFactor)
end

# the point is that only the (0,20) values in newFactor are needed, all calculations are abstracted away.
```

See also: [`RoME.generateCanonicalFG_Honeycomb!`](@ref), [`accumulateFactorMeans`](@ref), [`getPPE`](@ref)
"""
function _checkVariableByReference( fg::AbstractDFG,
                                    srcLabel::Symbol,
                                    destRegex::Regex,
                                    destType::Type{<:InferenceVariable},
                                    factor::AbstractRelative;
                                    srcType::Type{<:InferenceVariable} = getVariableType(fg, srcLabel) |> typeof,
                                    doRef::Bool = true,
                                    refKey::Symbol=:simulated,
                                    prior = !doRef ? nothing : DFG._getPriorType(srcType)( MvNormal(getPPE(fg[srcLabel], refKey).suggested, diagm(ones(getDimension(srcType)))) ),
                                    atol::Real = 1e-2,
                                    destPrefix::Symbol = match(r"[a-zA-Z_]+", destRegex.pattern).match |> Symbol,
                                    srcNumber = match(r"\d+", string(srcLabel)).match |> x->parse(Int,x),
                                    overridePPE = doRef ? nothing : zeros(getDimension(destType))  )
  #

  refVal = if overridePPE !== nothing
    overridePPE
  else
    # calculate and add the reference value
    # TODO refactor consolidation to use `_buildGraphByFactorAndTypes!`
    tfg = initfg()
    addVariable!(tfg, :x0, srcType )
    addFactor!(tfg, [:x0], prior )
    addVariable!(tfg, :l0, destType )
    addFactor!( tfg, [:x0; :l0], factor, graphinit=false )

    # calculate where the landmark reference position is
    accumulateFactorMeans(tfg, [:x0f1; :x0l0f1])
  end

  ppe = DFG.MeanMaxPPE(refKey, refVal, refVal, refVal)

  # now check if we already have a landmark at this location
  varLms = ls(fg, destRegex) |> sortDFG
  already = if doRef
    ppeLms = getPPE.(getVariable.(fg, varLms), refKey) .|> x->x.suggested
    errmask = ppeLms .|> ( x -> isapprox(x, refVal; atol=atol) )
    any(errmask)
  else
    false
  end

  if already
    # does exist, ppe, variableLabel
    alrLm = varLms[findfirst(errmask)]
    # @info "Variable on :$refKey does exists at" srcLabel alrLm
    return true, ppe, alrLm
  end

  # Nope does not exist, ppe, generated new variable label only
  return false, ppe, Symbol(destPrefix, srcNumber)
end

function _checkVariableByReference( fg::AbstractDFG,
                                    srcLabel::Symbol,
                                    destRegex::Regex,
                                    destType::Type{<:InferenceVariable},
                                    factor::AbstractPrior;
                                    srcType::Type{<:InferenceVariable} = getVariableType(fg, srcLabel) |> typeof,
                                    doRef::Bool = true,
                                    refKey::Symbol=:simulated,
                                    prior = typeof(factor)( MvNormal(getMeasurementParametric(factor)...) ),
                                    atol::Real = 1e-3,
                                    destPrefix::Symbol = match(r"[a-zA-Z_]+", destRegex.pattern).match |> Symbol,
                                    srcNumber = match(r"\d+", string(srcLabel)).match |> x->parse(Int,x),
                                    overridePPE = doRef ? nothing : zeros(getDimension(destType))  )
  #

  refVal = if overridePPE !== nothing
    overridePPE
  else
    getMeasurementParametric(factor)[1]
  end

  ppe = DFG.MeanMaxPPE(refKey, refVal, refVal, refVal)

  # Nope does not exist, ppe, generated new variable label only
  return false, ppe, Symbol(destPrefix, srcNumber)
end

#


text:
"""
Mutable reference to an array element.
In Julia v1.4 we had some performance regressions problems that could be solved using
a `MutRef`.
"""
mutable struct MutRef{T} <: Ref{T}
    x::Vector{T}
    i::Int
end

@inline Base.getindex(b::MutRef) = b.x[b.i]
@inline Base.setindex!(b::MutRef, x) = (b.x[b.i] = x; b)

ext: jl

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
module CUDA

using GPUCompiler

using GPUArrays

using LLVM
using LLVM.Interop

using Adapt

using Requires

using LinearAlgebra

## source code includes

const root = dirname(@__DIR__)

include("$root/src/pointer.jl")

# core library
include("$root/lib/utils/APIUtils.jl")
include("$root/lib/cudadrv/CUDAdrv.jl")

# essential stuff
include("$root/src/initialization.jl")
include("$root/src/state.jl")

# binary dependencies
include("$root/deps/discovery.jl")
include("$root/deps/compatibility.jl")
include("$root/deps/bindeps.jl")

# device functionality (needs to be loaded first, because of generated functions)
include("$root/src/device/pointer.jl")
include("$root/src/device/array.jl")
include("$root/src/device/cuda.jl")
include("$root/src/device/llvm.jl")
include("$root/src/device/runtime.jl")

# compiler libraries
include("$root/lib/cupti/CUPTI.jl")
include("$root/lib/nvtx/NVTX.jl")
export CUPTI, NVTX

# compiler implementation
include("$root/src/compiler/gpucompiler.jl")
include("$root/src/compiler/execution.jl")
include("$root/src/compiler/exceptions.jl")
include("$root/src/compiler/reflection.jl")

# array abstraction
include("$root/src/memory.jl")
include("$root/src/array.jl")
include("$root/src/gpuarrays.jl")
include("$root/src/subarray.jl")
include("$root/src/utilities.jl")

# array libraries
include("$root/lib/complex.jl")
include("$root/lib/library_types.jl")
include("$root/lib/cublas/CUBLAS.jl")
include("$root/lib/cusparse/CUSPARSE.jl")
include("$root/lib/cusolver/CUSOLVER.jl")
include("$root/lib/cufft/CUFFT.jl")
include("$root/lib/curand/CURAND.jl")
include("$root/lib/cudnn/CUDNN.jl")
include("$root/lib/cutensor/CUTENSOR.jl")
export CUBLAS, CUSPARSE, CUSOLVER, CUFFT, CURAND, CUDNN, CUTENSOR

# integrations and specialized functionality
include("$root/src/indexing.jl")
include("$root/src/broadcast.jl")
include("$root/src/mapreduce.jl")
include("$root/src/accumulate.jl")
include("$root/src/linalg.jl")
include("$root/src/nnlib.jl")
include("$root/src/iterator.jl")
include("$root/src/statistics.jl")

include("$root/src/deprecated.jl")

end


text:
module SimpleRecoTest

using FactorizationMachines
using Base.Test

# simple recommendation example

 #1 is User1
 #2 is User2
 #3 is Item1,
 #4 is Item2,
 #5 is Item3,
 #6 is Item4
 #7 is the category “18-25″,
 #8 is the category “26-40″,
 #9 is the category “40-60″
 #10 will represent the price feature

#rating     #user   #items     #category    #price

T = [
5           1 0     1 0 0 0    1 0 0        12.5;
5           1 0     0 1 0 0    1 0 0        20;
4           1 0     0 0 1 0    1 0 0        78;
1           0 1     1 0 0 0    0 0 1        12.5;
1           0 1     0 1 0 0    0 0 1        20;
]

info("Testing reading from libsvm format...")

(XFromFile, yFromFile) = read_libsvm("data/small_train.libfm")

println(XFromFile)
println(size(XFromFile))
@test size(XFromFile) == (10, 5)

X = sparse(T[:,2:end])
y = T[:,1]

fm = train(sparse(XFromFile),yFromFile)

(TFromFile, tFromFile) = read_libsvm("data/small_test.libfm")
p = predict(fm,sparse(TFromFile))

info("Predictions: $p")

end


text:
using PyPlot:
    imshow,
    xlabel,
    ylabel,
    contour,
    figure,
    pause,
    hold,
    axis,
    title,
    scatter,
    gcf,
    savefig,
    matplotlib,
    rc,
    tick_params,
    clf,
    cla
function vis2(m::SearchDomain, uav::SimUnit;
                   pause_time=0.3,
                   alpha=0.1,
                   show_mean::Bool=false,
                   show_cov::Bool=false,
                   save_gif::Bool=true
                  )

    frames = Frames(MIME("image/png"), fps=5)

    x_frames = Frames(MIME("image/png"), fps=5)
    es_frames = Frames(MIME("image/png"), fps=5)

    #println("rc_context = ", rc_context())
    #rc("font", family="Times New Roman", size=16)
    #rc("text", usetext=true)

    figure("Simulation")
    plot(m, uav.f, uav.x, alpha=alpha)
    #title("i = 0")
    title("\$t\$ = 0 s")
    push!(frames, gcf())

    # What was the cost to getting this first observation?
    cost_sum = get_cost(uav, m)

    # before doing anything else, we observe
    #  and update filter once
    o = observe(m, uav.x)
    update!(uav, o)

    # This was our first step; steps count number of observations
    step_count = 1

    # plot if need be
    plot(m, uav.f, uav.x, alpha=alpha)
    #title("i = $(step_count)")
    title("\$t\$ = $(step_count) s")
    push!(frames, gcf())

    while !is_complete(uav.f, uav.tc, step_count)
        # act
        #tic()
        a = action(m, uav, o)
        figure("Simulation")
        push!(x_frames, gcf())
        figure("MI")
        push!(es_frames, gcf())
        #ta = toq()
        #println("time = ", ta)
        act!(m, uav.x, a)

        move_target!(m)

        # get cost and update step count
        cost_val = get_cost(uav, m, a)
        cost_sum += cost_val
        step_count += 1

        # observe and update
        o = observe(m, uav.x)
        update!(uav, o)

        # plot if need be
        pause(pause_time)
        figure("Simulation")
        cla()
        plot(m, uav.f, uav.x, alpha=alpha, show_mean=show_mean, show_cov=show_cov)

        #max_b = maximum(uav.f.b)
        #title("i = $(step_count), cost = $(round(cost_val,3))")
        #title("\$t\$ = $(step_count) s, cost = $(round(cost_val,3))")
        title("\$t\$ = $(step_count) s")
        push!(frames, gcf())

        # temporary...
        #cm = covariance(uav.f)      # covariance matrix
        #ev = eigvals( covariance(uav.f) )
        #println("ev[1] = ", round(sqrt(ev[1]),3), ", ev[2] = ", round(sqrt(ev[2]),3))
        #num_out = 0
        #parts = uav.f._b.particles
        #for praxar in parts
        #    if praxar[1] > m.length || praxar[1] < 0.0 || praxar[2] > m.length || praxar[2] < 0.0
        #        num_out += 1
        #    end
        #end
        #p_out = round(100.0 * num_out / length(parts), 2)
        #println("percentage of particles out: ", p_out)
        #ce = cheap_entropy(uav.f, m.length, 50)
        #println("\tcheap_entropy = ", round(ce,3))
        #ce2 = cheap_entropy2(uav.f, m.length, 50)
        #println("\tcheap_entropy2 = ", round(ce2,3))

    end
    write("temp.gif", frames)

    figure("Simulation")
    push!(x_frames, gcf())
    figure("MI")
    push!(es_frames, gcf())

    write("temp1.gif", x_frames)
    write("temp2.gif", es_frames)

    return cost_sum
end


text:
function print_log(logfile::String, printargs; write_log=true, print_to_stdout = true)

    if print_to_stdout
        println(printargs...)
    end

    if write_log
        open(logfile, "a") do lf
            # modify mile content
            for arg in printargs
                write_lf_entries(lf, arg)
            end
            write(lf, "\n")
        end
    end

    return nothing
end

function write_lf_entries(lf, arg)
    if typeof(arg) <: Array
        write(lf,"[")
        for (i,arg_i) in enumerate(arg)
            write_lf_entries(lf, arg_i)
            if i != length(arg)
                write(lf,",")
            end
        end
        write(lf,"]")
    else
        if typeof(arg) <: Number
            write(lf,string(arg))
        elseif typeof(arg) == String
            write(lf, arg)
        elseif typeof(arg) == Char
            write(lf, string(arg))
        else
            println(typeof(arg))
            error("Unknow type for writing to log file")
        end
    end
    return nothing
end

function new_log(logfile::String)
    if isfile(logfile)
        println("Logfile \"",logfile ,"\"allready exist. Continue and overwrite? y/N");
        s = readline()
        if (s == "Y" || s ==  "y")
            rm(logfile)
        else
            error("Manual abort, due to existing logfile!")
        end
    end

    println("Set m1d1d_logfile varible to ", logfile)
    global m1d1d_logfile = logfile

    print_log(m1d1d_logfile, string(Dates.now()))
end

function test_global_logfilename()
    println(m1d1d_logfile)
end


text:
#Functions to be implemented for Dataset sources:
"Return a DiskArray handle to a dataset"
get_var_handle(ds, name) = ds[name]

"Return a list of variable names"
function get_varnames end

"Return a list of dimension names for a given variable"
function get_var_dims end

"Return a dict with the attributes for a given variable"
function get_var_attrs end

"Return a dict with global attributes for the dataset"
function get_global_attrs end

#Functions to be implemented for Dataset sinks
"Initialize and return a handle to a new empty dataset"
function create_empty end

"""
    add_var(ds, T, name, s, dimlist, atts)

Add a new variable to the dataset with element type `T`,
name `name`, size `s` and depending on the dimensions `dimlist`
given by a list of Strings. `atts` is a Dict with attributes.
"""
function add_var end

"""
    allow_parallel_write(ds)

Returns true if different chunks of a dataset can be written to by
2 processes simultaneously
"""
allow_parallel_write(ds) = false

"""
    allow_missings(ds)

Returns true if Union{T,Missing} is an allowed data type for the backend
and if an array containing missings can be written to the array.
"""
allow_missings(ds) = true

#Fallback for writing array
function add_var(ds,x::AbstractArray,name,dimlist,atts;kwargs...)
  a = add_var(ds,eltype(x),name,size(x),dimlist,atts;kwargs...)
  a .= x
  a
end

function create_dataset(T::Type, path, gatts, dimnames, dimvals, dimattrs, vartypes, varnames, vardims, varattrs, varchunks; kwargs...)
  ds = create_empty(T, path, gatts)
  axlengths = Dict{String, Int}()
  for (dname, dval, dattr) in zip(dimnames, dimvals, dimattrs)
    add_var(ds, dval, dname, (dname,), dattr)
    axlengths[dname] = length(dval)
  end
  for (T, vn, vd, va, vc) in zip(vartypes, varnames, vardims, varattrs, varchunks)
    s = getindex.(Ref(axlengths),vd)
    add_var(ds, T, vn, (s...,), vd, va; chunksize = vc, kwargs...)
  end
  ds
end

using DataStructures: OrderedDict
"""
    YAXArrayBase.backendlist

List of symbols defining dataset backends.
"""
backendlist = OrderedDict{Symbol, Any}(
  :array => Array,
)

backendregex = Pair[]

function to_dataset(g::String; driver=:all)
  if driver == :all
    for p in YAXArrayBase.backendregex
      if match(p[1],g) !== nothing
        return to_dataset(p[2],g)
      end
    end
    return to_dataset(g,driver=:zarr)
  else
    return to_dataset(backendlist[driver],g)
  end
end
to_dataset(g; kwargs...) = g
to_dataset(T::Type{<:Any}, g::String) = T(g)


text:
using HTbin:SpaInput, get_za!,
            _loadapexsh, _apxg2q, _apxq2g,
            atm8_chapman,
            hwm_dir, hwm14,
            gtd7

function test_all()
    qlatout, qlonout, glatout, glonout = HTbin.test_apex()
    @test isapprox(qlatout, 16.185024)
    @test isapprox(qlonout, 85.22224)
    @test isapprox(glatout, 25.0)
    @test isapprox(glonout, 11.000001)

    ch1 = HTbin.test_chapman()
    @test isapprox(ch1, 1.154317169142657)

    OUTF, OARR = HTbin.test_iri()
    @test isapprox(OUTF[1], 1.0536364e11)

    u, v = HTbin.test_hwm()
    @test isapprox(u,  -2843.706703186035)
    @test isapprox(v, -34.654366970062256)

    d, temp = HTbin.test_gtd7()
    @test isapprox(d[1], 2.0281334e12)
    @test isapprox(d[2], 0.0)
    @test isapprox(d[3], 3.0222293e17)
    @test isapprox(d[4], 8.107734e16)
    @test isapprox(d[5], 3.6150396e15)
    @test isapprox(d[6], 1.8594203e-5)
    @test isapprox(d[7], 0.0)
    @test isapprox(d[8], 0.0)
    @test isapprox(d[9], 0.0)
    @test isapprox(temp[1], 1027.3185)
    @test isapprox(temp[2], 230.39586)

    zenithout, azimuthout = HTbin.test_spa()
    @test isapprox(zenithout, 9.936110030111635)
    @test isapprox(azimuthout,  202.08736805795803)
end

test_all()


text:
using RigidBodyDynamics
using RigidBodyDynamics: Bounds
using MeshCat
using MeshCatMechanisms
using LCPSim
using JuMP, Gurobi
using Base.Test

@testset "rimless wheel" begin
    urdf = joinpath(@__DIR__, "urdf", "rimless_wheel.urdf")
    mechanism = parse_urdf(Float64, urdf)
    floating_base = findjoint(mechanism, "floating_base")
    position_bounds(floating_base) .= Bounds(-100, 100)
    velocity_bounds(floating_base) .= Bounds(-100, 100)

    vis = Visualizer()
    if !haskey(ENV, "CI")
        open(vis)
        wait(vis)
    end
    mvis = MechanismVisualizer(mechanism, URDFVisuals(urdf), vis[:robot])

    x0 = MechanismState(mechanism)
    set_configuration!(x0, [-1.5, 1.8, -0.05])
    env = LCPSim.parse_contacts(mechanism, urdf);
    controller = LCPSim.passive_controller()
    Δt = 0.005
    N = 1700
    lcp_solver = GurobiSolver(Gurobi.Env(); OutputFlag=0)

    results = LCPSim.simulate(x0, controller, env, Δt, N, lcp_solver)

    for r in results
        set_configuration!(mvis, configuration(r.state))
        sleep(Δt)
    end
end


text:
"""
Function to obtain series of states out of DAE Solution. It receives the solution, the dynamical system
and a tuple containing the symbol name of the Dynamic Injection device and the symbol of the state.

"""

function get_state_series(sim::Simulation, ref::Tuple{String, Symbol})
    global_state_index = get_global_index(sim.simulation_inputs)
    ix = get(global_state_index[ref[1]], ref[2], nothing)
    return sim.solution.t, [value[ix] for value in sim.solution.u]
end

"""
Function to obtain the voltage magnitude series out of the DAE Solution. It receives the solution, the dynamical system and the bus number.

"""
function get_voltagemag_series(sim::Simulation, bus_number::Int64)
    n_buses = get_bus_count(sim.simulation_inputs)
    bus_ix = get(get_lookup(sim.simulation_inputs), bus_number, nothing)
    if isnothing(bus_ix)
        error("Bus number $(bus_number) not found.")
    else
        return sim.solution.t,
        [sqrt(value[bus_ix]^2 + value[bus_ix + n_buses]^2) for value in sim.solution.u]
    end
end

"""
Function to print initial states. It receives the vector of initial states and the dynamical system.
"""
function print_init_states(sim::Simulation)
    for (ix, val_sys) in get_global_index(sim.simulation_inputs)
        println("Differential States")
        println(ix)
        println("====================")
        for (k, val) in val_sys
            print(k, " ", sim.x0_init[val], "\n")
        end
        println("====================")
    end

    # println("Algebraic States") # TODO: Print Buses Voltages
    return
end

"""
Function to print initial states. It receives the vector of initial states and the dynamical system.
"""
function print_device_states(sim::Simulation)
    bus_size = get_bus_count(sim.simulation_inputs)
    system = get_system(sim.simulation_inputs)
    println("Voltage Variables")
    println("====================")
    buses_sorted =
        sort(collect(PSY.get_components(PSY.Bus, system)); by = x -> PSY.get_number(x))
    for bus in buses_sorted
        name = PSY.get_name(bus)
        println(name)
        println("====================")
        bus_n = PSY.get_number(bus)
        bus_ix = get_lookup(sim.simulation_inputs)[bus_n]
        V_R = sim.x0_init[bus_ix]
        V_I = sim.x0_init[bus_ix + bus_size]
        Vm = sqrt(V_R^2 + V_I^2)
        θ = angle(V_R + V_I * 1im)
        print("Vm ", round(Vm, digits = 4), "\n")
        print("θ ", round(θ, digits = 4), "\n")
        println("====================")
    end
    println("====================")
    for device in PSY.get_components(PSY.DynamicInjection, system)
        states = PSY.get_states(device)
        name = PSY.get_name(device)
        println("Differential States")
        println(name)
        println("====================")
        global_index = get_global_index(sim.simulation_inputs)[name]
        for s in states
            print(s, " ", round(sim.x0_init[global_index[s]], digits = 4), "\n")
        end
        println("====================")
    end
    dyn_branches = PSY.get_components(PSY.DynamicBranch, system)
    if !isempty(dyn_branches)
        println("====================")
        println("Line Current States")
        for br in dyn_branches
            states = PSY.get_states(br)
            name = PSY.get_name(br)
            printed_name = "Line " * name
            println("====================")
            println(printed_name)
            global_index = get_global_index(sim.simulation_inputs)[name]
            x0_br = Dict{Symbol, Float64}()
            for (i, s) in enumerate(states)
                print(s, " ", round(sim.x0_init[global_index[s]], digits = 5), "\n")
            end
            println("====================")
        end
    end
    return
end


text:
using NCDatasets
using Test

filename = tempname()

ds = Dataset(filename,"c")
x = collect(1:10)
defVar(ds, "x", x, ("x",))
defDim(ds, "Time", Inf)
sync(ds)
defVar(ds, "Time", Float64, ("Time",))

defVar(ds, "a", Float64, ("x", "Time"))
defVar(ds, "u", Float64, ("x", "Time"))
defVar(ds, "v", Float64, ("x", "Time"))
defVar(ds, "w", Float64, ("x", "Time"))

for i in 1:10
    ds["Time"][i] = i
    ds["a"][:,i] = 1
    @test_throws NCDatasets.NetCDFError ds["u"][:,i] = collect(1:9)
    @test_throws NCDatasets.NetCDFError ds["v"][:,i] = collect(1:11)
    @test_throws NCDatasets.NetCDFError ds["w"][:,i] = reshape(collect(1:20), 10, 2)

    # ignore singleton dimension
    ds["w"][:,i] = reshape(collect(1:10), 1, 1, 10, 1)
end

ds["w"][:,:] = ones(10,10)

# w should grow along the unlimited dimension
ds["w"][:,:] = ones(10,15)
@test size(ds["w"]) == (10,15)

# w cannot grow along a fixed dimension
@test_throws NCDatasets.NetCDFError ds["w"][:,:] = ones(11,15)

# NetCDF: Index exceeds dimension bound
@test_throws NCDatasets.NetCDFError ds["u"][100,100]
close(ds)


text:
data = [
    8.9 14.0 4.3 19.9 2.1 28.0 3.6 1.3 4.3;
    13.5 9.3 4.1 17.5 4.5 26.6 5.7 2.1 4.0;
    18.0 9.9 3.3 19.5 5.7 28.1 4.8 2.4 6.5;
    13.9 10.0 4.7 25.8 2.2 24.0 6.2 1.6 2.9;
    9.5 13.6 3.6 23.4 2.5 22.4 4.2 1.8 3.7;
    13.1 10.1 3.1 23.8 2.3 25.6 2.8 2.4 4.9;
    17.4 5.7 4.7 20.6 4.3 24.3 4.7 3.4 3.3;
    11.4 12.5 4.1 18.8 3.4 18.6 5.2 1.5 3.8
]

N = size(data)[1]
prototypes = copy(data)
clusters = [[i] for i = 1:N]
C = length(clusters)

d²(x⃗, y⃗) = sum((x⃗ - y⃗).^2)
d′(x⃗ᵢ, x⃗ⱼ, nᵢ, nⱼ) = ((nᵢ * nⱼ) / (nᵢ + nⱼ)) * sum((x⃗ᵢ - x⃗ⱼ).^2)

T = 2 # C

for t = 1:T
    D = [i<j ? d′(prototypes[i, :], prototypes[j, :], length(clusters[i]), length(clusters[j])) : (i == j ? Inf : Inf) for i = 1:C, j = 1:C]
    _, idx = findmin(D)
    i, j = idx[1], idx[2]

    append!(clusters[i], clusters[j])
    deleteat!(clusters, j)
    C = length(clusters)

    prototypes = Matrix{Float64}(undef, C, 9)
    for c = 1:C
        vectors = [data[k, :] for k in clusters[c]]
        mean = sum(vectors) / length(vectors)
        prototypes[c, :] = mean
    end

    variances = zeros(C)
    for r = 1:C
        for k in clusters[r]
            variances[r] += d²(data[k, :], prototypes[r, :])
        end
    end
    total_variance = sum(variances)
    println("total variance = $total_variance")
end


text:
using Gen
import Random

@gen function model()
    if ({:z} ~ bernoulli(0.5))
        m1 = ({:m1} ~ gamma(1, 1))
        m2 = ({:m2} ~ gamma(1, 1))
    else
        m = ({:m} ~ gamma(1, 1))
        (m1, m2) = (m, m)
    end
    {:y1} ~ normal(m1, 0.1)
    {:y2} ~ normal(m2, 0.1)
end

@gen function mean_random_walk_proposal(trace)
    if trace[:z]
        {:m1} ~ normal(trace[:m1], 0.1)
        {:m2} ~ normal(trace[:m2], 0.1)
    else
        {:m} ~ normal(trace[:m], 0.1)
    end
end

@gen function split_merge_proposal(trace)
    if trace[:z]
        # currently two means, switch to one
    else
        # currently one mean, switch to two
        {:u} ~ uniform_continuous(0, 1)
    end
end

function merge_mean(m1, m2)
    m = sqrt(m1 * m2)
    u = m1 / (m1 + m2)
    (m, u)
end

function split_mean(m, u)
    m1 = m * sqrt((u / (1 - u)))
    m2 = m * sqrt(((1 - u) / u))
    (m1, m2)
end

@involution function split_merge_involution(
        model_args, proposal_args, proposal_retval)
    if @read_discrete_from_model(:z)
        # currently two means, switch to one
        @write_discrete_to_model(:z, false)
        m1 = @read_continuous_from_model(:m1)
        m2 = @read_continuous_from_model(:m2)
        (m, u) = merge_mean(m1, m2)
        @write_continuous_to_model(:m, m)
        @write_continuous_to_proposal(:u, u)
    else
        # currently one mean, switch to two
        @write_discrete_to_model(:z, true)
        m = @read_continuous_from_model(:m)
        u = @read_continuous_from_proposal(:u)
        (m1, m2) = split_mean(m, u)
        @write_continuous_to_model(:m1, m1)
        @write_continuous_to_model(:m2, m2)
    end
end

function do_inference_simple(y1, y2)
    trace, = generate(model, (), choicemap((:y1, y1), (:y2, y2), (:z, false), (:m, 1.2)))
    zs = Bool[]
    m = Float64[]
    m1 = Float64[]
    m2 = Float64[]
    for iter=1:100
        trace, = mh(trace, select(:z))
        trace, = mh(trace, mean_random_walk_proposal, ())
        push!(zs, trace[:z])
        push!(m, trace[:z] ? NaN : trace[:m])
        push!(m1, trace[:z] ? trace[:m1] : NaN)
        push!(m2, trace[:z] ? trace[:m2] : NaN)
    end
    (zs, m, m1, m2)
end

function do_inference_rjmcmc(y1, y2)
    trace, = generate(model, (), choicemap((:y1, y1), (:y2, y2), (:z, false), (:m, 1.2)))
    zs = Bool[]
    m = Float64[]
    m1 = Float64[]
    m2 = Float64[]
    for iter=1:100
        trace, = mh(
            trace, split_merge_proposal, (), split_merge_involution; check=true)
        trace, = mh(trace, mean_random_walk_proposal, ())
        push!(zs, trace[:z])
        push!(m, trace[:z] ? NaN : trace[:m])
        push!(m1, trace[:z] ? trace[:m1] : NaN)
        push!(m2, trace[:z] ? trace[:m2] : NaN)
    end
    (zs, m, m1, m2)
end

y1, y2 = (1.0, 1.3)
(zs, m, m1, m2) = do_inference_rjmcmc(y1, y2)
println(zs)

function plots()
    Random.seed!(2)

    figure(figsize=(6, 3))

    subplot(2, 2, 1)
    y1, y2 = (1.0, 1.3)
    (zs, m, m1, m2) = do_inference_rjmcmc(y1, y2)
    plot(m, label="m")
    plot(m1, label="m1")
    plot(m2, label="m2")
    title("Involution MH (RJMCMC)")
    legend(loc="lower right")
    gca().set_ylim(0.5, 1.5)

    subplot(2, 2, 3)
    plot(zs, label="z", color="black")
    xlabel("\\# MCMC moves")
    legend(loc="center right")
    yticks(ticks=[0, 1], labels=["F", "T"])
    gca().set_ylim(-0.1, 1.1)

    subplot(2, 2, 2)
    y1, y2 = (1.0, 1.3)
    (zs, m, m1, m2) = do_inference_simple(y1, y2)
    plot(m, label="m")
    plot(m1, label="m1")
    plot(m2, label="m2")
    title("Selection MH")
    legend(loc="lower right")
    gca().set_ylim(0.5, 1.5)

    subplot(2, 2, 4)
    plot(zs, label="z", color="black")
    xlabel("\\# MCMC moves")
    legend(loc="center right")
    yticks(ticks=[0, 1], labels=["F", "T"])
    gca().set_ylim(-0.1, 1.1)

    tight_layout()
    savefig("rjmcmc.png")
end

using PyPlot
plots()


text:
"""
Estimate the period of `x`, from an initial guess `L`. The result will be the
peak of the continuous cross-correlation between `0.5L` and `1.5L`, i.e. the
period is not constrained to be an integer.

The algorithm is to first compute a discrete cross-correlation and oversample by
4x to account for inter-sample peaks. We use the peak of the oversampled
cross-correlation as a starting point for a numerical optimization of the
average-squared-distance function in the frequency domain.
"""
function getperiod(x::AbstractVector, L; max_ppm=500_000)
    maxD = L * max_ppm ÷ 1_000_000

    # create two signals that are offset by L (our expected period length), and
    # zero-pad by to do a linear (rather than circular) convolution. We set up
    # the padding to pre-shift so the cross-correlation starts start at their
    # left-most lag (`-maxD`), and the `maxD`th sample will actually be L lag)
    # we also add extra padding to account for time aliasing from the fractional
    # delay. Otherwise it can wrap around. 512 samples is enough for a sinc
    # function to decay more than 60dB

    # note this could be problematic if the signal is temporally compact, e.g.
    # an impulse train. Innacuracy in `L` could chop off an impulse from one or
    # the other of these.
    x1 = [zeros(maxD); x[L+1:end]; zeros(512)]
    x2 = [x[1:(end-L)]; zeros(maxD); zeros(512)]

    m = length(x1)
    X1 = rfft(x1)
    X2 = rfft(x2)
    M = length(X1)

    # do a discrete cross-correlation and upsample. The maximum here should be
    # a pretty good estimate of the true peak. Should be valid from 1:8L, which
    # represents lags from -L:L (upsampled by 4)
    # not sure why `xc` wasn't being type-inferrred`
    xc = resample(irfft(X1 .* conj.(X2), m)[1:2maxD+1], 4)::Vector{Float64}
    D_est1 = (findmax(xc)[2]-1)/4-maxD

    # note if `m` were guaranteed to be even the upper bound here would be π
    ω = range(0, (M-1)*2π/m, length=M)

    # error function for  frequency-domain delay estimation. This can be
    # interpreted as a frequency-domain ASDF (average-squared distance function)
    # as a function of the continuous delay
    err(D) = sum(abs2(X1[i] - X2[i] * exp(-im*ω[i]*(D+maxD)))
                 for i in 1:M)

    # opt = optimize(D->err(first(D)), [D_est1], NewtonTrustRegion(), autodiff=:forward)
    opt = optimize(err, D_est1-0.5, D_est1+0.5)
    if !converged(opt)
        throw(ErrorException("Delay estimation failed to converge"))
    end

    # return the delay and the normalized minimum error
    norm = 2sqrt(sum(abs2, X1)*sum(abs2, X2))
    L+minimizer(opt)[1], minimum(opt)/norm
end


text:
# The 2D model:
# - interpolate data onto flow band
# - pass it to 1D model
# - extrapolate data back to 2D

################
# Forward model 2D
################

# """
# The forward model needs to implement the function:

#     fwdm{P<:APara, M<:APara, N<:APara}(dem, b, fq, po::P, mp::M, nu::N) -> (h, us)

# Input:

# - 2D fields (or a scalar):
#   - dem: DEM of the glacier, includes a mask of the glacierized cells
#   - bdot: mass balance
#   - dhdt: dh/dt
#   - fq: fraction of deformational ice flux
#   - A: ice flow parameter
# - pp: misc. physical parameters: physical constants
# - mp: misc. model parameters: somewhat physical things, but quite ad-hoc
# - np: misc. numerical parameters: e.g. discretisation size, etc.

# Output:
# - h: ice thickness
# - us: surface flow speed

# Units: use SI (m,kg,s)
# """
# function fwdm end

# NOTE: do not access gb.fsl, gb.btilde, gb.temp, etc here!  Use
# values passed in as function arguments

"""
    fwdm(gl::Glacier, pp::Phys, pm::MPara, pn::Num)
    fwdm(gb, pp::Phys, pm::MPara, pn::Num,
              btilde = gb.bdot-gb.dhdt,
              fsl = gb.fsl,
              temp = gb.temp,
              )

The 2D forward model function does the following steps:
- interpolate onto 1D (this can be done separately with `make_bands!(gl,pp,pm,pn)`)
- run 1D model
- extrapolate onto 2D

Input:
- `fwdm(gl::Glacier, pp::Phys, pm::MPara, pn::Num)`
- `fwdm(gb, pp::Phys, pm::MPara, pn::Num[, btilde, fsl, temp])` does
  no 2D->1D and uses the 1D-fields if provided.  Otherwise uses
  gb.bdot, etc.

Parameters (to be fitted):
- A: ice flow factor (scalar)
- fq: ice sliding ratio (field)

Return:

    hs1d, taus1d, ivs1d, hs2d, taus2d, ivs2d, taus1d_local, vol_ratio, gb

Units: use SI (m,kg,s) except for flow velocity and SMB (.../a)

TODO: add method which takes the map_onto_band variables
"""
function fwdm(gl::Glacier, pp::Phys, pm::MPara, pn::Num)
    # map onto bands
    gb, pm = make_bands(gl, pp, pm, pn)
    return fwdm(gb, pp, pm, pn, gb.bdot-gb.dhdt, gb.fsl, gb.temp)
end

function fwdm(gb, pp::Phys, pm::MPara, pn::Num,
              btilde, fsl, temp)

    gl = gb.gl
    # call fwdm1d
    hs1d_n, taus1d_local, taus1d_n, ivs1d_n, qtot1d_n = fwdm1d(gb, btilde, fsl, temp, pp, pm, pn)
    # extrapolate back to 2D
    hs2d, taus2d, ivs2d, ivs2d_band, ivs2d_band_mask, vol, vol_above =
        extrapolate(hs1d_n, taus1d_n, ivs1d_n, btilde, fsl, gl, gb, pp, pm)
    vol_ratio = check_vol_1D_2D(gb, hs1d_n, vol)
    if pn.verbose && !(0.9<vol_ratio<1.1)
        warn("Volume of 1D and 2D model is more than 10% different.  Ratio 1D/2D=$vol_ratio")
    end
    return FWDSol(hs1d_n, ivs1d_n, taus1d_n, taus1d_local,
                  hs2d, ivs2d, ivs2d_band, ivs2d_band_mask, taus2d,
                  vol, vol_above, vol_ratio,
                  gb, pp, pm, pn, btilde, fsl, temp, Dict{Symbol,Any}())
end

################
# 2D->1D
################
#
# The 2D data gets binned into elevation bands.

"""

Holds elevation bands information for a glacier, which are geometric
properties and 1D fields (their prior expectation value in a Bayesian
context).  Band indices start at the top of the glacier.

Also holds entities to allow extrapolation from bands back to 2D.

Note that gb.bands are the bands starting from the top and their location corresponds to the band-start.
gb.ele gives the band midpoints, gb.bands_ is like gb.bands but includes also the end-point (terminus).

"""
@with_kw struct Bands{Gl<:Glacier}
    gl::Gl
    dem::Gridded{F} # possibly smoothed surface DEM [m] (if pm.window_dem_smooth>0) otherwise original
    bands::SRangeL{Float64} # elevation bands starts.  Runs from top (divide) to bottom (terminus)
    bands_::SRangeL{Float64} # elevation bands starts and end
    bandi::Vector{Vector{Int}} # indices into gl.dem.v for each band
    ele::SRangeL{Float64} # mean elevation of bands, i.e. bands mid-points
    ws::Vector{F}       # width
    ls::Vector{F}       # length
    malphas::Vector{F}  # mean slope of a band
    x::Vector{F}        # distance from divide to beginning of bands (i.e. to nodes)
    xmid::Vector{F}     # distance from divide to middle of bands (i.e. cell-centers)

    # Smoothed 1D fields which represent the expectation value of
    # their prior distribution.
    # TODO: figure out how to band the sigmas of gl.bdot, etc.
    bdot::Vector{F} # [m ice/a]
    dhdt::Vector{F} # [m ice/a]
    fsl::Vector{F} # []
    temp::Vector{F} # [C]
    # original fields (not smoothed)
    fields1d_original::Dict{Symbol,Vector{F}}

    # 1D->2D extrapolation:
    # weights used in extrapolation (TODO put elsewhere?)
    weights::Vector{Vector{F}}
    # precualculated boxcar filter matrix
    boxcar_M::Tuple{SparseMatrixCSC{F,Int},F}
    boxcar_M_iv::Tuple{SparseMatrixCSC{F,Int},F,Vector{Dict{Int,Vector{VAWTools.Line}}},
                       Matrix{F}, Matrix{F}}
    # 2D surface slope calculated with a possibly smoothed dem (the dem in above field)
    alpha2d::Matrix{F}
end
# function Base.show(io::IO, gb::Bands)
#     print(io, "Elevation bands discretisation of:\n  ")
#     Base.show(io, gb.gl)
# end
Base.length(gb::Bands) = length(gb.bands)

"Return a 1D field, including the non-field btilde"
function getfield1d(gb, name)
    if name==:btilde
        return apparent_mass_balance(gb.bdot, gb.dhdt)
    else
        return getfield(gb, name)
    end
end

"""
$(SIGNATURES)

Calculate bands.  This is one of the core functions!

In:
- glacier
- pp
- MPara: model parameters

Out:
- Bands
"""
function make_bands(gl::Glacier, pp::Phys, pm::MPara, pn::Num)
    @unpack dem, glaciermask = gl
    @unpack bandsize, window_width_smooth, window_dem_smooth,
            slope_min, slope_max, iv_window_frac, window = pm
    @unpack calc_boxcar_M = pn

    @assert bandsize>0 "Need a pm.bandsize>0"

    bands, bandi, malphas, areas, lengths, widths, x, xmid, dem_smoothed, alpha2d_smoothed =
        make_1Dglacier(dem, -bandsize, glaciermask;
                       window_dem_smooth=window_dem_smooth,
                       window_width_smooth=window_width_smooth,
                       alpha_min=slope_min, alpha_max=slope_max,
                       FILL=FILL,
                       min_bands=3,
                       min_bin_number_ends=5,
                       verbose=pn.verbose
                       )

    ele = bands+step(bands)/2 # mean elevation of each band
    pm = MPara(pm, bandsize=-step(bands))

    # h-extrapolation weights `sin(alpha)^(-n/(n+2))`
    weights = extrapolation_weigths(bands, bandi, widths, lengths,
                                       pm.alpha_min, alpha2d_smoothed, pp.n)

    # thickness smoothing operator
    dx = step(dem.x)
    if calc_boxcar_M
        boxcar_M = boxcar_matrix(F, round(Int,window/dx), gl.glaciermask, (!).(gl.glaciermask))
        window_boxcar_M = window
        # IV smoothing operator for :masscons method
        boxcar_M_iv, boundaries, ux, uy = VAWTools.get_iv_boxcar_M(F, dem_smoothed, glaciermask, bands, bandi, lengths, iv_window_frac)
        window_boxcar_M_iv = iv_window_frac
    else # set to dummies which will trigger the use of the straight
         # boxcar in extrapolate
        boxcar_M = speye(0)
        window_boxcar_M = -9999.0
        boxcar_M_iv = speye(0)
        window_boxcar_M_iv = -9999.0
        boundaries = Dict{Int,Vector{VAWTools.Line}}[]
        ux = zeros(F,0,0)
        uy = zeros(F,0,0)
    end

    # 1d fields: see what is in gl and do the appropriate thing.
    # Steps:
    # - if 2D: interpolate to 1D & smooth
    # - if 1D: interpolate to elevation bands & smooth
    # - if scalar: make constant for all elevation bands

    fields1d = Dict{Symbol,Vector{F}}()
    fields1d_original = Dict{Symbol,Vector{F}}()
    for name in [:iv, :bdot, :dhdt, :fsl, :temp]
        glfld = getfield(gl, name)
        if isa(glfld, SomeData{Gridded{F}}) # 2D
            gr = glfld.vals
            bandi_ = VAWTools.bandi_for_other_grid(bands, bandi, dem_smoothed, gr, gr.v.!=FILL)
            vals = map_onto_bands(bandi_, gr.v, mean, FILL)
            vals = [(isnan(v) || v==FILL) ? glfld.vals_prior : v for v in vals]
            for v in vals
                if name!=:iv && (isnan(v) || v==FILL)
                    @show name, v, glfld.vals_prior
                end
            end
            fields1d_original[name] = vals
            # smooth: This does not conserve mass!  TODO is this good?
            fields1d[name] = smooth_vector(ele, vals, glfld.smooth_window)
        elseif isa(glfld, SomeData{Gridded1d{F}}) # 1D
            # interpolate onto elevation bands
            if glfld.vals.x==ele # no interpolation
                fields1d_original[name] = glfld.vals.v
                # smooth: This does not conserve mass!  TODO is this good?
                fields1d[name] = smooth_vector(ele, glfld.vals.v, glfld.smooth_window)
            else
                # linear interpolation
                f = piecewiselinear(glfld.vals.x, glfld.vals.v)
                fields1d_original[name] = f.(ele)
                # smooth: This does not conserve mass!  TODO is this good?
                fields1d[name] = smooth_vector(glfld.vals.x, glfld.vals.v,
                                               glfld.smooth_window, ele)
            end
        elseif isa(glfld, SomeData{F}) # 0D
            nb = length(bands)
            fields1d_original[name] = zeros(F,nb) + glfld.vals
            fields1d[name] = zeros(F,nb) + glfld.vals
        elseif isa(glfld, SomeData{Traj{F}})
            fields1d_original[name] = glfld.vals.v
            fields1d[name] = F[] # TODO
        else
            if name==:iv && isa(glfld, SomeData{Void})
                fields1d_original[name] = F[]
                fields1d[name] = F[]
            else
                error("Need to provide a value for field $name in Glacier datastruct.")
            end
        end
        if name!=:iv
            any(isnan.(fields1d_original[name])) && error("1D field $name has NaN: $(fields1d_original[name])")
            any(isnan.(fields1d[name])) && error("Smoothed 1D field $name has NaN: $(fields1d_original[name])")
        end
    end

    return Bands{typeof(gl)}(gl, dem_smoothed,
                             bands, first(bands):step(bands):last(bands)+step(bands), bandi,
                             ele, widths, lengths, malphas, x, xmid,
                             fields1d[:bdot], fields1d[:dhdt], fields1d[:fsl], fields1d[:temp],
                             fields1d_original,
                             weights, (boxcar_M, window_boxcar_M),
                             (boxcar_M_iv, window_boxcar_M_iv, boundaries, ux, uy),
                             alpha2d_smoothed), pm::MPara
end

VAWTools.map_onto_bands(gb::Bands, field, fn=mean) = map_onto_bands(gb.bandi, field, fn, FILL)
VAWTools.bandi_for_other_grid(gb::Bands, othergrid::Gridded, othermask=trues(size(othergrid.v)) ) =
    VAWTools.bandi_for_other_grid(gb.bands, gb.bandi, gb.dem_smoothed, othergrid, othermask)

"""
    map_iv_onto_bands(gl::Glacier, gb::Bands)

Maps the IV in `gl.iv` (which is usually on a different grid to
`gl.dem`) onto the elevation bands.

Does not take mass-flux into account.

Not in use.
"""
function map_iv_onto_bands(gl::Glacier, gb::Bands)
    bandi = VAWTools.bandi_for_other_grid(gb, gl.iv.vals, gl.ivmask)
    return map_onto_bands(bandi, gl.iv.vals.v, mean, FILL)
end

"""
    extrapolation_weigths(bands, bandi, ws, ls,
                               alpha_min, alpha2d, n)

Pre-calculates the weights used in the extrapolation due to local
slope.
"""
function extrapolation_weigths(bands, bandi, ws, ls,
                               alpha_min, alpha2d, n)
    weights = Vector{F}[]
    for (i,band) in enumerate(bands)  # @inbounds does not help
        w, l = ws[i], ls[i]
        inds = bandi[i]
        weight = similar(inds, F)
        for (j,jj) in enumerate(inds)
            alpha_ = max(alpha_min, alpha2d[jj])
            weight[j] = sin(alpha_)^( -n/(n+2))
        end
        push!(weights, weight)
    end
    return weights
end

# """
# Values given on bands.
# """
# struct BandValues{N<:Number}
#     gb::Bands
#     v::Vector{N}
# end
# BandValues{T}(gb::Bands, field::Gridded{T}) = BandValues{T}(gb, map_onto_bands(b, field, mean, FILL))

################
# 1D->2D
################
#

"""
$(SIGNATURES)

    -> hs2d, taus2d, ivs2d, ivs2d_band, ivs2d_band_mask

Extrapolation of the output of the 1D model to 2D.

TODO:
- use directional/an-isotropic smoothing
"""
function extrapolate(hs_n, taus_n, ivs_n, btilde, fsl, gl::Glacier, gb::Bands,
                               pp::Phys, pm::MPara)
    @unpack n, rho, rhosw, g, r = pp
    @unpack alpha_min, window, cap_at_floatation,
            dist_land_maxdist, dist_exp,
            iv_extrapolation_scaling, iv_tau_exp, iv_h_exp, iv_dist_exp1,
            iv_dist_exp2, iv_nye,
            iv_flux_dir_window, iv_window_frac, iv_scale_factor = pm
    @unpack dist_land, dist_land_maxdist_calc, glaciermask = gl
    @unpack weights, ws, ls, bands, bandi, boxcar_M, boxcar_M_iv = gb
    @assert dist_land_maxdist_calc>=dist_land_maxdist "Calculated distance in gl.dist_land is smaller than pm.dist_land_maxdist."

    dx = step(gl.dem.x)
    hs2d = zeros(gl.dem.v)
    taus2d = zeros(gl.dem.v)  # effective driving stress
    ivs2d = zeros(gl.dem.v)

    hs = on_cells(hs_n)
    taus = on_cells(taus_n)
    ivs = on_cells(ivs_n)

    ## Local ice thickness weighted proportional to
    # - (distance to land)^dist_exp
    # - (sin α)^(-n/(n+2))
    for (i,band) in enumerate(bands)  # @inbounds does not help
        h, w, l = hs[i], ws[i], ls[i]
        vol = w*l*h
        inds = bandi[i]
        weight = weights[i] # pre-calculated weight due to surface slope
        vol2d = 0.0
        for (j,jj) in enumerate(inds)
            hs2d[jj] = weight[j] * min(dist_land[jj], dist_land_maxdist)^dist_exp
            vol2d += hs2d[jj]*dx^2
        end
        # scale to conserve total volume of band (note this is non-linear in h)
        for jj in inds
            hs2d[jj] *= vol/vol2d
        end
    end

    # Smooth.  Does not quite conserve mass (because of edge effects),
    # also increases thickness at the edge!
    #
    # TODO this takes a lot of time, even with the sparse multiply;
    # in particular for larger smoothing windows.
    # Consider: for the fitting, one could only smooth at the radar points.
    if boxcar_M[2]==window
        hs2d = apply_boxcar_matrix(boxcar_M[1], hs2d)
    else
        # hs2d = boxcar(hs2d, round(Int,window/dx), ones(size(hs2d)), !gl.glaciermask)
        hs2d = boxcar(hs2d, round(Int,window/dx), gl.glaciermask, .!gl.glaciermask)
    end

    # Remove bits above flotation (does not conserve mass) and calculate
    # volume above flotation.
    dA = step(gl.dem.x) * step(gl.dem.y)
    vol, vol_above = 0.0, 0.0
    fac = 1/((rhosw-rho)/rhosw)
    for j=1:length(hs2d)
        !(glaciermask[j]) && continue
        zs = gl.dem.v[j]
        h = hs2d[j]
        hmax = zs*fac
        if cap_at_floatation && h>hmax
            hs2d[j] = hmax
            h = hmax
        end
        vol += h*dA
        # thickness such that just grounded
        hmin = max(0.0, - (zs - h) * rhosw/rho)
        vol_above += (h-hmin)*dA
    end

    ## Tau
    # Clarke et al 2013 assume that tau is the same over an elevation
    # band, we do this too:
    map_back_to_2D!(taus2d, bandi, taus)
    ## Ice velocity
    for (i,band) in enumerate(bands)  # @inbounds does not help
        # First make a trial u
        inds = bandi[i]
        for jj in inds
            d = min(dist_land[jj], dist_land_maxdist)

            len_scale = min(ws[i]/2,hs[i])
            d2 = min(dist_land[jj], min(len_scale, dist_land_maxdist) )
            # u-trial made with extrapolation, which can happen in a few different ways.
            # It is a trial-solution in the sense that it has not the correct scale
            # but the desired cross-valley/closest-to-outcrop shape.
            ivs2d[jj] =
                 !iv_nye * (
                    taus[i]^iv_tau_exp * # surface slope term on band
                    hs2d[jj]^iv_h_exp *  # ice thickness term
                    d^iv_dist_exp1 ) + # distance to land term 1
                  iv_nye * (
                      1 - ( (len_scale-d2)/len_scale )^iv_dist_exp2 ) # Nye 1965 inspired, Eq 6
        end

        # Scale the trial such that it "resembles" 1D field.
        if iv_extrapolation_scaling==:huss
            # The makes the 2D mean equal to the 1D flow.  Not mass conserving.
            mean_ = mean(ivs2d[inds])
            for jj in inds
                ivs2d[jj] = ivs[i]*ivs2d[jj]/mean_
            end
        elseif iv_extrapolation_scaling==:masscons
            # The makes the 2D flow such that the flux from one
            # elevation band to the next is correct.
            #
            # Whole step is performed below
        elseif iv_extrapolation_scaling==:none
            # Just scale by a constant.  Probably bad.
            for jj in inds
                ivs2d[jj] *= iv_scale_factor # a scale factor.  Only needed if no scaling is applied later.
            end
        else
            error("Unsupported value of iv_extrapolation_scaling = $iv_extrapolation_scaling")
        end
    end
    # Smooth IV
    if iv_extrapolation_scaling==:huss || iv_extrapolation_scaling==:none
        if boxcar_M[2]==window
            ivs2d = apply_boxcar_matrix(boxcar_M[1], ivs2d)
        else
            ivs2d = boxcar(ivs2d, round(Int,2*window/dx), gl.glaciermask, .!gl.glaciermask)
        end
    elseif iv_extrapolation_scaling==:masscons
        # This is both scaling and smoothing in one
        scale_surface_iv = 1./((1-r)*fsl + r)
        if boxcar_M_iv[2]==iv_window_frac # use the pre-calculated smoothing matrix
            boundaries = boxcar_M_iv[3]
            bM = boxcar_M_iv[1]
            ux, uy = boxcar_M_iv[4], boxcar_M_iv[5]
            ivs2d, ivs2d_band, ivs2d_band_mask = VAWTools.calc_u(flux_n(gb, btilde)*year2sec,
                                                     boundaries, ivs2d, hs2d,
                                                     ux, uy, dx, glaciermask, bands, bandi, ls,
                                                     iv_flux_dir_window,
                                                     bM, scale_surface_iv)
        else
            ux,uy = VAWTools.get_ux_uy(gb.dem, gl.glaciermask)
            binmat = VAWTools.bins2matrix(gl.dem, bands, bandi)
            boundaries = VAWTools.calc_boundaries(bands,bandi,binmat)
            ivs2d, ivs2d_band, ivs2d_band_mask = VAWTools.calc_u(flux_n(gb, btilde)*year2sec, boundaries, ivs2d, hs2d,
                                                                 ux, uy, dx, glaciermask, bands, bandi, ls,
                                                                 iv_flux_dir_window,
                                                                 iv_window_frac, scale_surface_iv)
        end
    else
        error("Unsupported value of iv_extrapolation_scaling = $iv_extrapolation_scaling")
    end

    return hs2d, taus2d, ivs2d, ivs2d_band, ivs2d_band_mask, vol, vol_above
end

"""
Forward model solution

```
@with_kw struct FWDSol
    # Fields from the 1D calculation
    hs1d::Vector{F}
    ivs1d::Vector{F}
    taus1d::Vector{F}
    taus1d_local::Vector{F}

    # Fields from the 2D extrapolation
    hs2d::Matrix{F}
    ivs2d::Matrix{F} # the smoothed ivs2d_band
    ivs2d_band::Matrix{F} # the calculated, mass conserving IV across bands
    ivs2d_band_mask::BitMatrix
    taus2d::Matrix{F}

    vol::Float64 # total ice volume (calculated from 2D field)
    vol_above::Float64 # ice volume above flotation
    vol_ratio::Float64 # ration between ice volume calculated from 1D and 2D field.
                       # Gives an indication of how much volume as not conserved in the extrapolation.

    # model input
    gb::Bands
    pp::Phys
    pm::MPara
    pn::Num
    btilde::Vector{F}
    fsl::Vector{F}
    temp::Vector{F}

    # to add various other stuff
    misc::Dict{Symbol,Any}
end
```
"""
@with_kw struct FWDSol
    # Fields from the 1D calculation
    hs1d::Vector{F}
    ivs1d::Vector{F}
    taus1d::Vector{F}
    taus1d_local::Vector{F}

    # Fields from the 2D extrapolation
    hs2d::Matrix{F}
    ivs2d::Matrix{F} # the smoothed ivs2d_band
    ivs2d_band::Matrix{F} # the calculated, mass conserving IV across bands
    ivs2d_band_mask::BitMatrix
    taus2d::Matrix{F}

    vol::Float64 # total ice volume (calculated from 2D field)
    vol_above::Float64 # ice volume above flotation
    vol_ratio::Float64 # ration between ice volume calculated from 1D and 2D field.
                       # Gives an indication of how much volume as not conserved in the extrapolation.

    # model input
    gb::Bands
    pp::Phys
    pm::MPara
    pn::Num
    btilde::Vector{F}
    fsl::Vector{F}
    temp::Vector{F}

    # to add various other stuff
    misc::Dict{Symbol,Any}
end

volume(gb::Bands, hs1d_n) = sum(on_cells(hs1d_n).*gb.ws.*gb.ls)
volume(hs2d, dx) = sum(hs2d)*dx^2
volume(fwdsol) = volume(fwdsol.hs2d, step(fwdsol.gb.gl.dem.x))

##########
# Checks
##########
"Check volume of 1D vs 2D.  Returns mass1D/mass2D"
check_vol_1D_2D(gb, hs1d_n, vol) = volume(gb, hs1d_n)/vol


text:
using ITensors,
      Test

@testset "QN" begin

  @testset "QNVal Basics" begin
    qv = ITensors.QNVal()
    @test !isactive(qv)

    qv = ITensors.QNVal("Sz",0)
    @test ITensors.name(qv) == ITensors.SmallString("Sz")
    @test val(qv) == 0
    @test modulus(qv) == 1
    @test isactive(qv)

    qv = ITensors.QNVal("A",1,2)
    @test ITensors.name(qv) == ITensors.SmallString("A")
    @test val(qv) == 1
    @test modulus(qv) == 2
    @test !isfermionic(qv)

    qv = ITensors.QNVal("Nf",1,-1)
    @test ITensors.name(qv) == ITensors.SmallString("Nf")
    @test val(qv) == 1
    @test modulus(qv) == -1
    @test isfermionic(qv)
  end

  @testset "QN Basics" begin
    q = QN()
    @test length(sprint(show,q)) > 1

    q = QN(("Sz",1))
    @test length(sprint(show,q)) > 1
    @test isactive(q[1])
    @test val(q,"Sz") == 1

    q = QN("Sz",1)
    @test length(sprint(show,q)) > 1
    @test isactive(q[1])
    @test val(q,"Sz") == 1

    q = QN("P",1,2)
    @test length(sprint(show,q)) > 1
    @test isactive(q[1])
    @test val(q,"P") == 1
    @test modulus(q,"P") == 2

    q = QN(("A",1),("B",2))
    @test isactive(q[1])
    @test isactive(q[2])
    @test val(q,"A") == 1
    @test val(q,"B") == 2
    @test modulus(q,"A") == 1
    @test modulus(q,"B") == 1

    q = QN(("B",2),("A",1))
    @test val(q,"A") == 1
    @test val(q,"B") == 2
  end

  @testset "Comparison" begin
    @test QN() == QN()
    @test QN("A",1) == QN("A",1)
    @test QN(("A",1),("B",3)) == QN(("A",1),("B",3))
    @test QN(("A",1),("B",3)) == QN(("B",3),("A",1))

    # Zero value and missing sector treated the same:
    @test QN(("A",0),("B",3)) == QN("B",3)
    @test QN(("B",3),("A",0)) == QN("B",3)
  end

  @testset "Arithmetic" begin
    @test QN("Sz",1) + QN() == QN("Sz",1)
    @test QN("Sz",1) + QN("Sz",2) == QN("Sz",3)
    @test QN("Sz",1) + QN("Sz",-2) == QN("Sz",-1)

    @test QN(("A",1),("Sz",0)) + QN(("A",0),("Sz",1)) == QN(("A",1),("Sz",1))

    @test QN("P",0,2) + QN("P",1,2) == QN("P",1,2)
    @test QN("P",1,2) + QN("P",1,2) == QN("P",0,2)

    # Arithmetic involving mixed-label QNs
    @test QN()-QN("Sz",2) == QN("Sz",-2)
    @test QN("Sz",2)-QN() == QN("Sz",2)
    @test QN()-QN(("Sz",2),("N",1)) == QN(("Sz",-2),("N",-1))
    @test QN("N",1)-QN("Sz",2) == QN(("N",1),("Sz",-2))
  end

  @testset "Ordering" begin
    z = QN()
    qa = QN(("Sz",1),("Nf",1))
    qb = QN(("Sz",0),("Nf",2))
    qc = QN(("Sz",1),("Nf",2))
    qd = QN(("Sz",1),("Nf",2))
    qe = QN(("Sz",-1),("Nf",-2))

    @test !(z < z)
    @test !(qa < z)
    @test (z < qa)
    @test (z < qb)
    @test !(qb < z)
    @test (z < qc)
    @test !(qc < z)
    @test (z < qd)
    @test !(qd < z)
    @test !(z < qe)
    @test (qe < z)

    @test !(qa > qb)
    @test qb > qa
    @test !(qb == qa)
    @test (qb < qc)
    @test !(qc < qb)
    @test !(qc == qb)
    @test (qc == qd)
    @test !(qc < qd)
    @test !(qd < qc)
  end

  @testset "Hashing" begin
    @test hash(QN(("Sz",0))) == hash(QN())
    @test hash(QN("Sz",0)) == hash(QN("N",0))
    @test hash(QN(("Sz",1),("N",2))) == hash(QN(("N",2),("Sz",1)))
  end

end

nothing


text:
# this user supplied file creates a dictionary of arguments
# if a key is repeated, the last use of the key is used
# it is a little bit dangerous letting the user run arbitrary code
# as part of the solver
# now that this file is read inside a function, it is better encapsulated

arg_dict = Dict{Any, Any}(
# specify physcs, SBP operator, and solve type (homotopy continuation)
"physics" => "Euler",
"operator_type" => "SBPDiagonalE",
"dimensions" => 2,
"run_type" => 40,  # use homotopy to solve the nonlinear problem
"jac_method" => 2,
"jac_type" => 2,
"order" => 1,
"use_DG" => true,
"use_lps" => true,  # use LPS stabilization
"Flux_name" => "HLLFlux",  # HLL Riemann solver used for interface flux

"real_time" => false,
"volume_integral_type" => 2,
"Volume_flux_name" => "IRFlux",  # entropy stable volume integrals

# set up initial and boundary conditions
"IC_name" => "ICFreeStream",
"numBC" => 2,
"BC1" => [36],
"BC1_name" => "FreeStreamBC",
"BC2" => [5],
"BC2_name" => "noPenetrationESBC",  # entropy stable Euler wall BC

# problem-specific parameters
"aoa" => 1.25,
"Ma" => 0.8,

# mesh and geometry
"smb_name" => "../meshes/airfoil40_mesh1_curve_.smb",
"dmg_name" => ".null",  # if using Simmetrix tools, can be airfoil40.smd
"addShockCapturing" => true,  # enable shock capturing
"shock_capturing_name" => "Volume",  # volume dissipation term
"shock_sensor_name" => "SensorPP",  # Persons and Peraire's shock sensor
"homotopy_addBoundaryIntegrals" => true,  # include boundary dissipation in
                                          # homotopy function
# solve options
"res_abstol" => 1e-10,
"res_reltol" => 1e-10,
"itermax" => 500,
"output_freq" => 1,
"do_postproc" => false,
)


text:
## Distributions to JSON/Packed types

packDistribution(dtr::Categorical) = PackedCategorical(; p=dtr.p )
packDistribution(dtr::Uniform) = PackedUniform(; a=dtr.a, b=dtr.b )
packDistribution(dtr::Normal) = PackedNormal(; mu=dtr.μ, sigma=dtr.σ )
packDistribution(dtr::ZeroMeanDiagNormal) = PackedZeroMeanDiagNormal(; diag=dtr.Σ.diag )
packDistribution(dtr::ZeroMeanFullNormal) = PackedZeroMeanFullNormal(; cov=dtr.Σ.mat[:] )
packDistribution(dtr::DiagNormal) = PackedDiagNormal(; mu=dtr.μ, diag=dtr.Σ.diag )
packDistribution(dtr::FullNormal) = PackedFullNormal(; mu=dtr.μ, cov=dtr.Σ.mat[:] )
packDistribution(dtr::Rayleigh) = PackedRayleigh(; sigma=dtr.σ )

packDistribution(dtr::AliasingScalarSampler) = PackedAliasingScalarSampler(; domain=dtr.domain, weights=dtr.weights.values )

# A more specialized constructor to help serialization processes reading Any[Any[1,2,3,..]] rather than floats.
function PackedHeatmapGridDensity(
    _type::String,
    data::AbstractVector, # {Any}
    domain::AbstractVector, # {Any}
    hint_callback::String,
    bw_factor::Float64,
    N::Int64
  )
  #
  # TODO data might not be type Float64, should store and recover as performance enhancement (if user specified different element type)
  data_ = Vector{Vector{Float64}}(undef, length(data))
  for (i,dat) in enumerate(data)
    dat_ = replace(dat, nothing=>0)
    data_[i] = float.(dat_)
  end
  domain_ = tuple(float.(domain[1]), float.(domain[2]))

  PackedHeatmapGridDensity(_type, data_, domain_, hint_callback, bw_factor, N)
end

function packDistribution( obj::HeatmapGridDensity )
  #
  data_ = obj.data
  @cast data[j][i] := data_[i,j]
  # str = convert(SamplableBelief, obj.densityFnc)
  N = Npts(obj.densityFnc)
  # TODO misses the hint...
  PackedHeatmapGridDensity( "IncrementalInference.PackedHeatmapGridDensity",
                            data,
                            obj.domain,
                            "",
                            obj.bw_factor,
                            N )
end

packDistribution(dtr::LevelSetGridNormal) = PackedLevelSetGridNormal( "IncrementalInference.PackedLevelSetGridNormal",
                                                                      dtr.level,
                                                                      dtr.sigma,
                                                                      dtr.sigma_scale,
                                                                      convert(PackedHeatmapGridDensity, dtr.heatmap) )
#

function parchDistribution(hgd::HeatmapGridDensity)
  @assert 2 <= size(hgd.data,1) "parchDistribution of HeatmapGridDensity can only be done when `.data` is larger than 2x1"

  data = Matrix{eltype(hgd.data)}(undef, 2,2)
  data[1,1] = hgd.data[1,1]
  # data[2,2] = hgd.data[2,2] # disable since data might be a single column in unusual cases
  data[2,1] = size(hgd.data,1)
  data[1,2] = size(hgd.data,2)

  domain = hgd.domain
  hint_callback = hgd.hint_callback
  bw_factor = hgd.bw_factor
  densityFnc = parchDistribution(hgd.densityFnc)

  HeatmapGridDensity(
    data,
    domain,
    hint_callback,
    bw_factor,
    densityFnc
  )
end

## Unpack JSON/Packed to Distribution types

unpackDistribution(dtr::PackedCategorical) = Categorical( dtr.p ./ sum(dtr.p) )
unpackDistribution(dtr::PackedUniform) = Uniform(dtr.a, dtr.b )
unpackDistribution(dtr::PackedNormal) = Normal( dtr.mu, dtr.sigma )
unpackDistribution(dtr::PackedZeroMeanDiagNormal) = MvNormal( LinearAlgebra.Diagonal(map(abs2, sqrt.(dtr.diag))) ) # sqrt.(dtr.diag)
unpackDistribution(dtr::PackedZeroMeanFullNormal) = MvNormal( reshape(dtr.cov, length(dtr.mu), :) )
unpackDistribution(dtr::PackedDiagNormal) = MvNormal( dtr.mu, sqrt.(dtr.diag) )
unpackDistribution(dtr::PackedFullNormal) = MvNormal( dtr.mu, reshape(dtr.cov, length(dtr.mu), :) )
unpackDistribution(dtr::PackedRayleigh) = Rayleigh( dtr.sigma )

unpackDistribution(dtr::PackedAliasingScalarSampler) = AliasingScalarSampler( dtr.domain, dtr.weights ./ sum(dtr.weights) )

function unpackDistribution(obj::PackedHeatmapGridDensity)
  #
  # do intermediate conversions
  data_ = obj.data
  data__ = map(x->collect(x), data_)
  @cast data[i,j] := data__[j][i]
  _data__ = collect(data)
  # densFnc = convert(SamplableBelief, obj.densityFnc)
  # build the final object, misses the hint...
  HeatmapGridDensity( _data__,
                      obj.domain,
                      obj.hint_callback == "" ? nothing : nothing,
                      obj.bw_factor;
                      N=obj.N )
end

unpackDistribution(dtr::PackedLevelSetGridNormal) = LevelSetGridNormal( dtr.level,
                                                                        dtr.sigma,
                                                                        dtr.sigma_scale,
                                                                        convert(HeatmapGridDensity, dtr.heatmap) )
#

# ## strip field from NamedTuple

# function _delete( nt::Union{<:NamedTuple, <:Dict{K,T}},
#                   key::K=:_type ) where {K,T}
#   #
#   kys = keys(nt)
#   # rm index
#   ridx = findfirst(k->k==key, kys)
#   # keep indices
#   idxs = setdiff(1:length(nt), ridx)
#   # to Dict
#   dict = OrderedDict{K,Any}()
#   for id in idxs
#     ky = kys[id]
#     dict[ky] = nt[ky]
#   end
#   #

#   NamedTuple{Tuple(keys(dict))}(values(dict))
# end

## ===========================================================================================
## Converts still necessary?
## ===========================================================================================

# NOTE part of new effort to overhaul the SamplableBelief serialization approach
convert(::Type{<:PackedSamplableBelief}, obj::StringThemSamplableBeliefs) = packDistribution(obj)
convert(::Type{<:SamplableBelief}, obj::PackedSamplableBelief) = unpackDistribution(obj)

function convert(::Type{<:PackedSamplableBelief}, nt::NamedTuple)
  distrType = DFG.getTypeFromSerializationModule(nt._type)
  distrType(;nt...)
end

##===================================================================================

# FIXME ON FIRE, must deprecate nested JSON written fields in all serialization
# TODO is string necessary, because unpacking templated e.g. PackedType{T} has problems, see DFG #668
convert(::Type{String}, dtr::StringThemSamplableBeliefs) = JSON2.write(packDistribution(dtr))

function convert(::Type{<:SamplableBelief}, str_obj::AbstractString)
  #

  # go from stringified to generic packed (no type info)
  _pck = JSON2.read(str_obj)
  # NOTE, get the packed type from strong assumption that field `_type` exists in the
  T = DFG.getTypeFromSerializationModule( _pck[:_type] )
  # unpack again to described packedType
  pckT = JSON2.read(str_obj, T)

  # unpack to regular <:SamplableBelief
  return unpackDistribution(pckT)
end

#


text:
module test_pointattributelist
using Triangulate
function test()

    x=[1.9  1.2  1.5  1.6  1.3  1.9  1.2  1.4  1.1  1.0]
    y=[1.5  1.4  1.7  1.2  1.3  1.9  2.0  1.9  1.5  1.0]
    z=[1.7  1.0  1.1  1.6  2.0  1.8  1.7  1.8  1.7  1.5]

    pointset = Triangulate.TriangulateIO(; pointlist=[x; y])
    pointset.pointattributelist = z
    triout, vorout=Triangulate.triangulate("Q", pointset);

    test1= triout.trianglelist== [10   5  7   2  2   9  2  4  8  8  6  8  3;
                                  5  10  9   9  3   7  5  1  3  9  3  6  5;
                                  2   4  8  10  9  10  3  3  6  3  1  7  4]
    test2=    triout.pointattributelist==z

    triout, vorout=Triangulate.triangulate("Qa0.05", pointset);

    test3 = isapprox(triout.pointattributelist,[1.7  1.0  1.1  1.6  2.0  1.8  1.7  1.8  1.7  1.5  1.6  1.72149  1.55  1.42347  1.13548  1.50625],atol=1.0e-5)

    test4= triout.trianglelist== [13   5   7   2  12   9  14   9  16   2   3  12   6  8   9  14   5   7   1  15   3   8  14;
                                  2  13  12   9   9  11  15   7  14  13  16   3  16  6   2   3   4   9  14   2  12  16  16;
                                  10   4   8  10  15  10   5  11   1   5   8   8   1  7  15  15  14  12   4   5  15   6   3]
    test1&test2&test3&test4
end
end


text:
module ResizingTools

using ArrayInterface
using ArrayInterface: has_parent, parent_type
using LinearAlgebra: AdjOrTrans, AdjOrTransAbsVec, AdjOrTransAbsMat
using Static
using FunctionIndices

export Size, set!, to_dims, getsize, setsize!
export resize_buffer!, resize_buffer_dim!
export SimpleRDArray

# tools for define and access size
include("size.jl")
# core, defined sizehint! and resize! and related methods
include("resize.jl")
# some abstract type with some pre-defined methods
include("abstract.jl")
# some implementation of resizable array
include("example/simplerdarray.jl")

end # module


text:
module SecretSanta

import Combinatorics
import Dates
import GLPK
import JSON
import JuMP
import MathOptInterface
import Random
import SMTPClient

struct SecretSantaModel
    model::JuMP.Model
    data::Dict{String, Any}
    constraints::Dict{Symbol, Any} # JuMP constraint references.
    variables::Dict{Symbol, Any} # JuMP variable references.
    solution::Dict{String, Any} # Solution reference.
end

function SecretSantaModel(data::Dict{String, Any})
    # Create the set of the participants.
    P = data["participants"]

    # Create the node set.
    N = [x["email"] for x in P]

    # Create arcs for a complete bipartite graph.
    exclude = vcat([[[x["email"], y] for y in x["exclude"]] for x in P]...)
    A = collect(Combinatorics.combinations(N, 2)) # Collect arcs from i to j.
    A = vcat([reverse(a) for a in A], A) # Concatenate arcs from j to i.
    A = filter(x -> !(x in exclude), A) # Remove excluded arcs.
    A = [(a[1], a[2]) for a in A] # Convert to array of tuples.

    # Shuffle the array of arcs to induce random solutions.
    A = Random.shuffle(A)

    # Create the JuMP model.
    model = JuMP.Model(GLPK.Optimizer)
    variables = Dict{Symbol, Any}(:x => nothing)
    constraints = Dict{Symbol, Any}(:out_flow => nothing, :in_flow => nothing)
    constraints[:out_flow] = Dict{String, JuMP.ConstraintRef}()
    constraints[:in_flow] = Dict{String, JuMP.ConstraintRef}()

    # Create variables corresponding to arc selection.
    variables[:x] = JuMP.@variable(model, [a in A], binary=true, base_name="x")

    for i in N
        out_arcs = collect(filter(x -> (x[1] == i), A))
        out_vars = Array{JuMP.VariableRef}([variables[:x][a] for a in out_arcs])
        constraints[:out_flow][i] = JuMP.@constraint(model, sum(out_vars) == 1)

        in_arcs = collect(filter(x -> (x[2] == i), A))
        in_vars = Array{JuMP.VariableRef}([variables[:x][a] for a in in_arcs])
        constraints[:in_flow][i] = JuMP.@constraint(model, sum(in_vars) == 1)
    end

    solution = Dict{String, Any}()
    ssm = SecretSantaModel(model, data, constraints, variables, solution)
    return ssm # Return the SecretSantaModel instance.
end

function build_model(input_path::String)
    data = JSON.parsefile(input_path)
    return SecretSantaModel(data)
end

function solve_model(ssm::SecretSantaModel)
    JuMP.optimize!(ssm.model)

    if JuMP.termination_status(ssm.model) == MathOptInterface.OPTIMAL
        A = ssm.variables[:x].axes[1]
        return filter(a -> isapprox(JuMP.value(ssm.variables[:x][a]), 1.0), A)
    else
        error("Secret Santa assignment is not possible. Adjust participants.")
    end
end

function send_email(ssm::SecretSantaModel, sender::Dict{String,Any}, recipient::Dict{String,Any}, test::Bool=true)
    # Prepare the subject of the email.
    subject = ssm.data["email"]["subject"]
    recipient_name = recipient["name"]
    subject = replace(subject, "{recipient}" => recipient_name)

    # Prepare the body of the email.
    message = ssm.data["email"]["message"]
    sender_name = sender["name"]
    sender_email = sender["email"]
    recipient_email = recipient["email"]
    message = replace(message, "{sender}" => sender_name)
    message = replace(message, "{recipient}" => recipient_name)
    message = replace(message, "{recipient_email}" => recipient_email)

    time_now = Dates.now(Dates.UTC)
    datetime = Dates.format(time_now, "e, dd u yyyy HH:MM:SS")

    body = "Date: $(datetime) +0000\n" *
           "From: Santa Claus <$(ssm.data["email"]["username"])>\n" *
           "To: $(sender_email)\n" * "Subject: $(subject)\n" * "\n" *
           message * "\n"

    body_io = IOBuffer(body)

    # Prepare email sending options.
    opt = SMTPClient.SendOptions(isSSL=true,
                                 username=ssm.data["email"]["username"],
                                 passwd=ssm.data["email"]["password"])

    # Prepare the email.
    server = ssm.data["email"]["smtp_server"]
    port = string(ssm.data["email"]["smtp_port"])
    url = "smtps://$(server):$(port)"
    rcpt = ["$(sender_email)"]
    from = "$(ssm.data["email"]["username"])"

    if !test
        # Send the email.
        resp = SMTPClient.send(url, rcpt, from, body_io, opt)
    else
        println("------------------------------------------------------------")
        println("Message to $(sender_name) ($(sender_email))")
        println("Subject: $(subject)")
        println("$(message)")
        println("------------------------------------------------------------")
    end
end

function send_matchings(ssm::SecretSantaModel, solution::Array{Tuple{String, String}, 1}, test::Bool=true)
    participants = ssm.data["participants"]

    for matching in solution
        sender = findfirst(x -> x["email"] == matching[1], participants)
        recipient = findfirst(x -> x["email"] == matching[2], participants)
        send_email(ssm, participants[sender], participants[recipient], test)
    end
end

function run(input_path::String; test::Bool=true)
    ssm = build_model(input_path)
    solution = solve_model(ssm)
    send_matchings(ssm, solution, test)
end

end


text:
function move_forward!(bobby, Δx, traj, obj_traj)
    θ = bobby.pose[3]
    x = Δx[1]*cos(θ)
    y = Δx[2]*sin(θ)
    pose = [bobby.pose[1]+x, bobby.pose[2]+y, θ]
    push!(traj, pose)
    bobby.pose = pose
    if !isnothing(bobby.holding)
        obpose = [bobby.holding.pose[1]+x, bobby.holding.pose[2]+y, θ]
        bobby.holding.pose = obpose
        push!(obj_traj, [obpose...,bobby.holding.id])
    else
        push!(obj_traj, nothing)
    end
end

function translate!(bobby, pose, traj, obj_traj; tol=1e-2)
    while(abs(norm(pose[1:2]-bobby.pose[1:2])) > tol)
        u = bobby.kp*abs.(pose[1:2]-bobby.pose[1:2]) + bobby.kv*(bobby.velocity[1:2])
        move_forward!(bobby, u, traj, obj_traj)
    end
    return traj
end

function rotate!(bobby, Δθ, traj, obj_traj)
    θ = bobby.pose[3] + Δθ
    pose = [bobby.pose[1], bobby.pose[2], θ]
    push!(traj, pose)
    bobby.pose = pose
    if !isnothing(bobby.holding)
        R = [cos(θ) -sin(θ);
             sin(θ)  cos(θ)]
        Δxy = [0.6, 0.]
        oxy = R*Δxy
        obpose = [bobby.pose[1]+oxy[1], bobby.pose[2]+oxy[2], θ]

        bobby.holding.pose = obpose
        push!(obj_traj, [obpose...,bobby.holding.id])
    else
        push!(obj_traj, nothing)
    end
end

function turn!(bobby, angle, traj, obj_traj; tol=1e-2)
    while(abs(norm(angle-bobby.pose[3])) > tol)
        u = bobby.kp*(angle-bobby.pose[3]) + bobby.kv*(bobby.velocity[3])
        rotate!(bobby, u, traj, obj_traj)
    end
    return traj
end

#RTR motion controller
function go_to!(bobby, pose, traj, obj_traj)
    θ = atan(pose[2]-bobby.pose[2], pose[1]-bobby.pose[1])
    turn!(bobby, θ, traj, obj_traj)
    translate!(bobby, pose, traj, obj_traj)
    turn!(bobby, pose[3], traj, obj_traj)
end


text:
using ModelObjectsLH, ModelParams, Random, Test

mdl = ModelParams;

function param_info_test()
    @testset "ParamInfo" begin
        pName = :p1;
        startIdx = 10;
        valueV = [1.0 2.0 3.0; 4.1 5.2 6.3];
        lbV = fill(0.5, size(valueV));
        ubV = fill(9.3, size(valueV));
        pInfo = ParamInfo(pName, startIdx, lbV, ubV);
        pInfo2 = ParamInfo(pName, startIdx, lbV, ubV);

        @test mdl.n_values(pInfo) == length(lbV);
        @test mdl.indices(pInfo)[1] == startIdx;
        @test length(mdl.indices(pInfo)) == mdl.n_values(pInfo);
        @test isapprox(pInfo, pInfo2);
        rv = mdl.random_values(pInfo, MersenneTwister(121));
        @test size(rv) == size(lbV);

        pInfo = ParamInfo(:p1, 2, 0.1, 9.3);
        rv = mdl.random_values(pInfo, MersenneTwister(121));
        @test rv isa Float64;
        # set_value!(pInfo, 3.3);
        # @test pInfo.valueV == 3.3;
        # set_value!(pInfo, [3.5]);
        # @test pInfo.valueV == 3.5;

        # pInfo = ParamInfo(:p1, 10, [2.1, 2.2], [0.0, 0.0], [5.0, 3.8]);
        # set_value!(pInfo, [3.0 2.9]);
        # @test pInfo.valueV == [3.0, 2.9];
    end
end

function value_vector_test()
    @testset "ValueVector" begin
        n = 5;
        vv = mdl.make_test_value_vector(n);
        nValues = mdl.n_values(vv);
        # valueV = collect(LinRange(0.5, 0.6, nValues));
        # mdl.set_values!(vv, valueV);
        # @test valueV == mdl.get_values(vv);
        rValueV = mdl.random_guess(vv, MersenneTwister(123));
        @test size(rValueV) == (nValues, );
    end
end

# function guess_test(isCalibrated)
#     @testset "Guess" begin
#         n = 4;
#         pvecV = mdl.make_test_pvector_collection(n);
#         g = make_guess(pvecV, isCalibrated);
#         valueV = get_values(g);
#         set_values!(g, valueV);
#         value2V = get_values(g);
#         @test isapprox(valueV, value2V);
#     end
# end

@testset "Guess" begin
    param_info_test();
    value_vector_test();
    # for isCalibrated in (true, false)
    #     guess_test(isCalibrated);
    # end
end

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


text:
using LinearAlgebra
using SparseArrays

export warmstart!, initialstate!
export gettau, getdelta

"""
Create ScSTO problem with uniform switching cost
    scstoproblem(x0, dyna, d_dyna, uvec, [constr, d_constr, swc, ngrid, t0, tf, Q, E, tau0ws, name])
"""
function scstoproblem(
    x0::Vector{Float64},                        # initial state [nx]
    dynam::Function,                            # Nonlinear Dynamics
    ddynam::Function,                           # Nonlinear Dynamics Derivative
    U::Matrix{Float64};                         # inputs sequence [nu x N]
    t0::Float64 = 0.0,                          # initial time
    tf::Float64 = 1.0,                          # desired final time
    Q::Maybe{Matrix{Float64}} = nothing,        # state cost matrix [nx x nx]
    E::Maybe{Matrix{Float64}} = nothing,        # final state cost matrix [nx x nx]
    swc::Maybe{Float64} = nothing,              # switching cost
    constr::Maybe{Function} = nothing,          # Constraints
    dconstr::Maybe{Function} = nothing,         # Constraints Derivative
    pconstr::Maybe{Function} = nothing,         # Constraints projection
    ncon::Int = 0,                              # number of constraints
    ngrid::Int64 = 100,                         # Number of Linearization points in the fixed grid
    tau0ws::Maybe{Vector{Float64}} = nothing,   # switching times guess [N]
    name::String = "Generic",                   # problem name
)
    @assert ngrid > 1
    @assert tf > t0
    if swc < 0.0
        @warn "negative swc"
    end
    if ncon > 0
        @assert constr !== nothing
        @assert dconstr !== nothing
        @assert pconstr !== nothing
    else
        ncon = 0
    end
    nx = length(x0)                             # state dimension
    N = size(U, 2)                              # number of switching intervals
    if swc === nothing
        swc = 0.0                               # default switching cost
    end
    if Q === nothing
        Q = Diagonal(ones(Float64, nx))         # default cost matrix
    else
        @assert size(Q,1) == nx
        @assert size(Q,2) == nx
    end
    if E === nothing
        E = zeros(Float64, nx, nx)              # default final cost matrix
    else
        @assert size(E,1) == nx
        @assert size(E,2) == nx
    end
    if tau0ws === nothing
        tau0ws = collect(range(t0, stop = tf, length = N + 1))
        tau0ws = tau0ws[2:N]                    # default switching times guess
    else
        @assert length(tau0ws) == N - 1
    end

    # time vectors
    nvec = ngrid + N - 1
    # augmented system state and matrices
    x0 = [x0; 1]
    nx = nx + 1
    spz = Matrix{Float64}(undef, 1, 1)
    spz[1, 1] = 0.0
    spz = sparse(spz)  # Sparse scalar
    Q = Matrix(blockdiag(sparse(Q), spz))
    E = Matrix(blockdiag(sparse(E), spz))

    # problem data
    data = ScSTOdata(
        x0,
        nx,
        N,
        t0,
        tf,
        tf - t0,
        swc,
        Q,
        E,
        U,
        dynam,
        ddynam,
        constr,
        dconstr,
        pconstr,
        ngrid,
        nvec,
    )

    # discretization
    tgrid = collect(range(t0, stop = tf, length = ngrid)) # fixed time grid
    tvec = Vector{Float64}(undef, nvec)         # Complete time grid
    tauIdx = Vector{Int}(undef, N + 1)          # Indeces of switching times in the complete grid
    tgridIdx = Vector{Int}(undef, ngrid)        # Indeces of grid times in the complete grid
    dvec = Vector{Float64}(undef, nvec - 1)     # Intervals of comlete time grid
    # states on complete grid
    X = Matrix{Float64}(undef, nx, nvec)
    X[:, 1] = x0
    # matrices of linearized dynamics
    A = Array{Float64}(undef, nx, nx, nvec - 1)
    exm = Array{Float64}(undef, nx, nx, nvec - 1)
    Phi = Array{Float64}(undef, nx, nx, N + 1, N + 1)
    M = Array{Float64}(undef, nx, nx, nvec - 1)
    S = Array{Float64}(undef, nx, nx, nvec)
    C = Array{Float64}(undef, nx, nx, N)
    delta_old_C = Vector{Float64}(undef, N) # memory, to avoid recomputing
    delta_old_S = Vector{Float64}(undef, N) # memory, to avoid recomputing
    jttv = Vector{Float64}(undef, N-1)
    # evaluator
    eval = ScSTOeval(
        tf,
        tgrid,
        tvec,
        tauIdx,
        tgridIdx,
        dvec,
        A,
        X,
        exm,
        Phi,
        M,
        S,
        C,
        jttv,
        delta_old_C,
        delta_old_S,
    )

    # initial switching intervals
    delta0ws = tau2delta(tau0ws, t0, tf)
    # metadata
    meta = OptiModelMeta(N, ncon, x0=delta0ws, minimize=true, name=name)

    # counters
    ndyna = 0
    nobjf = 0
    ngrad = 0
    nobjg = 0
    nprox = 0
    ncons = 0
    ncjtv = 0
    nproj = 0
    # traces
    objf = Vector{Float64}(undef, 0)
    delta = Matrix{Float64}(undef, N, 0)
    # reporter
    repo = ScSTOrepo(ndyna, nobjf,ngrad, nobjg, nprox, ncons, ncjtv, nproj, objf, delta)

    return ScSTOptiModel(meta, data, eval, repo)
end

"""
    warmstart!( prob, tau )
"""
function warmstart!(p::ScSTOptiModel, tau::Vector{Float64})
    @assert length(tau) == p.data.N - 1
    delta = getdelta(p, tau)
    p.meta = OptiModelMeta(
        p.data.N,
        p.meta.ncon,
        x0 = delta,
        y0 = p.meta.y0,
        minimize = true,
        name = p.meta.name,
    )
    p.repo.objf = Vector{Float64}(undef, 0)
    p.repo.delta = Matrix{Float64}(undef, p.data.N, 0)
end

"""
    initialstate!( prob, x0 )
"""
function initialstate!(p::ScSTOptiModel, x0::Vector{Float64})
    @assert length(x0) == p.data.nx - 1
    p.data.x0 = [x0; 1]
end

"getdelta(prob, tau)"
function getdelta(p::ScSTOptiModel, tau::Vector{Float64})
    return tau2delta(tau, p.data.t0, p.data.tf)
end

"gettau(prob, delta)"
function gettau(p::ScSTOptiModel, delta::Vector{Float64})
    return delta2tau(delta, p.data.t0)
end


text:
# Copyright 2020 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using FirstOrderLp

# Built-in Julia packages.
import GZip
import LinearAlgebra
import Random
import Test: @test, @testset, @test_throws, @test_logs
import SparseArrays

const diagm = LinearAlgebra.diagm
const norm = LinearAlgebra.norm
const sparse = SparseArrays.sparse
const RestartScheme = FirstOrderLp.RestartScheme

Random.seed!(123)

include("utilities.jl")
include("shared_test_qp_problems.jl")
include("test_saddle_point.jl")
include("test_qp_processing.jl")
include("test_qp_io.jl")
include("test_sparse_linalg.jl")
include("test_mirror_prox.jl")
include("test_primal_dual_hybrid_gradient.jl")
include("test_iteration_stats.jl")
include("test_trust_region_utils.jl")
include("test_termination.jl")


text:
"""
    DiskTensors

Disk based rank 1,2, and 4 tensors. Defines functions such as +, -, dot,...
with buffered I/O. Uses HDF5 format for storing tensors. One tensor per file.

## structs
    DiskVector -> see ?Fermi.DiskTensors.DiskVector
    DiskMatrix -> see ?Fermi.DiskTensors.DiskMatrix
    DiskFourTensor -> see ?Fermi.DiskTensors.DiskFourTensor

## methods
    squeeze -> see ?Fermi.DiskTensors.squeeze
    eltype -> returns element type of a DiskFourTensor
    ranger -> converts input to a UnitRange, if possible
    blockfill! -> fills a DiskTensor with the value entered
"""
module DiskTensors

#include("DiskVectors.jl")
#include("DiskMatrices.jl")
include("DiskFourTensors.jl")

"""
	ranger
converts inp::Int64 -> UnitRange(inp:inp) for use in indexing DiskTensors
"""
function ranger(inp::Union{UnitRange{Int64},Int64,Colon})
    if typeof(inp) == Int64
        return UnitRange(inp:inp)
    else
        return inp
    end
end
export squeeze
export eltype
export ranger
#export DiskVector
export blockfill!
#export printdv
#export dvwrite!
#export dvread
#export dvdot
#export Base.getindex
#include("DiskVectorOps.jl")

#export DiskMatrix
#export dmread
#export dmwrite
#export dmdot

export DiskFourTensor
    #export dmelmult
end


text:
module BaryTrees

    using Plots, LinearAlgebra

    import Base.∈, Base.∉, Plots.Shape, Base.insert!

    export BaryTree, Point, insert!, ∈, ∉, Draw, Draw!, Query, BaryCentric,
           Cartesian2BaryCentric, BaryCentric2Cartesian

    const Point = Tuple{Float64,Float64}

    Point(x::Float64,y::Float64) = Point((x,y))

    struct BaryCentric
        """
            Represents the 3 cartesian (2D) points of a Barycentric
            coordinate system
        """
        r₁::Point
        r₂::Point
        r₃::Point
    end

    function BaryCentric2Cartesian(λ₁,λ₂,λ₃,p::BaryCentric)
        return Point(
            λ₁*p.r₁[1]+λ₂*p.r₂[1]+λ₃*p.r₃[1],
            λ₁*p.r₁[2]+λ₂*p.r₂[2]+λ₃*p.r₃[2],
        )
    end

    function Cartesian2BaryCentric(p::Point,b::BaryCentric)
        x,y = p
        x1,x2,x3 = b.r₁[1],b.r₂[1],b.r₃[1]
        y1,y2,y3 = b.r₁[2],b.r₂[2],b.r₃[2]

        l1 = ( (y2-y3)*(x-x3)+(x3-x2)*(y-y3) ) / ( (y2-y3)*(x1-x3) + (x3-x2)*(y1-y3) )
        l2 = ( (y3-y1)*(x-x3) + (x1-x3)*(y-y3) ) / ( (y2-y3)*(x1-x3) + (x3-x2)*(y1-y3) )
        return (
            l1,l2,1.0-l1-l2,
        )
    end

    function Medians(b::BaryCentric)::Vector{Point}
        """
            Returns the three mid points of each triangular edge of
            the triangle defining b in cartesian units
        """
        return [
            BaryCentric2Cartesian(1/2,0,1/2,b)
            BaryCentric2Cartesian(0,1/2,1/2,b)
            BaryCentric2Cartesian(1/2,1/2,0,b)
        ]
    end

    struct AxisAlignedTriangle
        """
            wrapper around BaryCentric
        """
        b::BaryCentric
    end

    Medians(t::AxisAlignedTriangle)::Vector{Point} = Medians(t.b)

    mutable struct BaryTree
        root::Bool
        boundary::AxisAlignedTriangle
        point::Union{Nothing,Point}
        c::Union{Nothing,BaryTree}
        t::Union{Nothing,BaryTree}
        l::Union{Nothing,BaryTree}
        r::Union{Nothing,BaryTree}
    end

    function ∈(p::Point,t::AxisAlignedTriangle)::Bool
        """
            Test if a cartesian point is inside a triangle
                using BaryCentric Coordinates of the triangle
        """
        b = Cartesian2BaryCentric(p,t.b)
        if (b[1] >= 0 && b[2] >= 0 && b[1]+b[2] < 1)
            return true
        else
            return false
        end
    end

    function points(a::AxisAlignedTriangle)::Vector{Point}
        """
            Short hand for the corners of an AxisAlignedTriangle
        """
        return  Point.([
                a.b.r₁,
                a.b.r₂,
                a.b.r₃
            ])
    end

    function intersects(a::AxisAlignedTriangle,b::AxisAlignedTriangle)::Bool
        """
            Tests if any of the corners of triangle a are in triangle b
        """
        for p in points(a)
            if p ∈ b
                return true
            end
        end
        return false
    end

    ∉(p::Point,t::AxisAlignedTriangle)::Bool = !∈(p,t)

    function subdivide(t::AxisAlignedTriangle)::Vector{AxisAlignedTriangle}
        """
            Subdivides a triangle (as in the Sierpinski construction) and returns
            them each initialised with their own coordinate systems
        """
        m1,m2,m3 = Medians(t)
        return AxisAlignedTriangle.([
            BaryCentric(t.b.r₁,m3,m1),
            BaryCentric(m3,m2,m1),
            BaryCentric(m3,t.b.r₂,m2),
            BaryCentric(m1,m2,t.b.r₃)
        ])

    end

    ∈(p::Point,q::BaryTree)::Bool = ∈(p,q.boundary)
    ∉(p::Point,q::BaryTree)::Bool = ∉(p,q.boundary)
    isempty(q::BaryTree)::Bool = q.point == nothing

    BaryTree(a::AxisAlignedTriangle)::BaryTree = BaryTree(false,a,nothing,nothing,nothing,nothing,nothing)
    BaryTree(r::Bool,a::AxisAlignedTriangle)::BaryTree = BaryTree(r,a,nothing,nothing,nothing,nothing,nothing)
    BaryTree(a::BaryCentric)::BaryTree = BaryTree(AxisAlignedTriangle(a))

    function subdivide!(q::BaryTree)::Nothing
        """
            Subdivide this tree node as in the Sierpinski construction
            and store it in t,c,l,r child nodes
        """
        trigs = subdivide(q.boundary)
        q.t = BaryTree(trigs[1])
        q.c = BaryTree(trigs[2])
        q.l = BaryTree(trigs[3])
        q.r = BaryTree(trigs[4])
        nothing
    end

    function insert!(p::Point,q::BaryTree)::Bool
        """
            Trial an insert of point p (cartesian) in q's
            triangular boundary
        """
        if (p ∉ q)
            return false
        end

        if isempty(q) && q.t == nothing
            subdivide!(q)
            if p ∈ q.t
                q.t.point = p
            end

            if p ∈ q.c
                q.c.point = p
            end
            if p ∈ q.l
                q.l.point = p
            end
            if p ∈ q.r
                q.r.point = p
            end
            return true
        end

        if (q.t == nothing)
            subdivide!(q)
        end

        if insert!(p,q.t)
            return true
        end
        if insert!(p,q.c)
            return true
        end
        if insert!(p,q.l)
            return true
        end
        if insert!(p,q.r)
            return true
        end

        return false
    end

    Shape(a::AxisAlignedTriangle) = Shape(
        Shape([a.b.r₁[1],a.b.r₂[1],a.b.r₃[1]],[a.b.r₁[2],a.b.r₂[2],a.b.r₃[2]])
    )

    function Draw!(a::AxisAlignedTriangle)
        plot!(Shape(a),label="",fillalpha=0.0)
    end

    function Draw(Q::BaryTree,points=true)
        p = plot(aspect_ratio=:equal,label="")
        Draw!(Q.boundary)
        for q in [Q.t,Q.c,Q.l,Q.r]
            q != nothing ? Draw!(q,points) : nothing
        end
        return p
    end

    function Draw!(Q::BaryTree,points=true)
        Draw!(Q.boundary)
        for q in [Q.t,Q.c,Q.l,Q.r]
            q != nothing ? Draw!(q) : nothing
        end
    end

    function query(a::AxisAlignedTriangle, Q::BaryTree)
        if ~intersects(a,Q.boundary)
            return Vector{Point}([])
        end

        if Q.point == nothing
            return Vector{Point}([])
        end

        result = Vector{Point}()

        if (Q.point ∈ a)
            push!(result,Q.point)
        end

        for q in [Q.t,Q.c,Q.l,Q.r]
            if q != nothing
                for p in query(a,q)
                    push!(result,p)
                end
            end
        end

        return result
    end

    function size(Q::BaryTree)
        s = 1
        for q in [Q.t,Q.c,Q.l,Q.r]
            if q != nothing
                s += size(q)
            end
        end
        return s
    end
end

max_stars_repo_name: Jerboa-app/BaryTrees

max_issues_repo_name: Jerboa-app/BaryTrees

max_forks_repo_name: Jerboa-app/BaryTrees


text:
@testset "6.2 Hyperbolic cosine" begin
    include("6.2.1 (c+d x)^m (a+b cosh)^n.jl")
    include("6.2.2 (e x)^m (a+b x^n)^p cosh.jl")
    include("6.2.3 (e x)^m (a+b cosh(c+d x^n))^p.jl")
    include("6.2.4 (d+e x)^m cosh(a+b x+c x^2)^n.jl")
    include("6.2.5 Hyperbolic cosine functions.jl")
    include("6.2.7 hyper^m (a+b cosh^n)^p.jl")
end

max_stars_repo_path: test/6 Hyperbolic functions/6.2 Hyperbolic cosine/.jl


text:
using PyPlot

function plot_fit(progr)
    m = length(progr)
    b = 0
    for j in eachindex(progr)
        n = size(progr[j],1)
        for i = 1:n
            subplot(m,n,b + i)
            plot(progr[j][i,:]', ".-")
        end
        b += n
    end
    tight_layout()
end

export plot_fit