text:
# new Sandshark example
# add more julia processes
using Distributed
addprocs(4)

using Caesar, RoME
@everywhere using Caesar, RoME
using Interpolations
using Distributions

using RoMEPlotting, ApproxManifoldProducts
using Gadfly, Fontconfig, Cairo
using DataFrames
using ProgressMeter
using DelimitedFiles

@everywhere setForceEvalDirect!(true)

# const TU = TransformUtils

Gadfly.set_default_plot_size(35cm,25cm)

include(joinpath(@__DIR__,"Plotting.jl"))
include(joinpath(@__DIR__,"SandsharkUtils.jl"))

odonoise = Matrix(Diagonal(10*[0.1;0.1;0.005].^2))

# Step: Selecting a subset for processing and build up a cache of the factors.
epochs = timestamps[50:3:300]

lastepoch = 0
for ep in epochs
  global lastepoch
  if lastepoch != 0
    # @show interp_yaw(ep)
    deltaAng = interp_yaw(ep) - interp_yaw(lastepoch)

    wXi = TU.SE2([interp_x(lastepoch);interp_y(lastepoch);interp_yaw(lastepoch)])
    wXj = TU.SE2([interp_x(ep);interp_y(ep);interp_yaw(ep)])
    iDXj = se2vee(wXi\wXj)
    NAV[ep] = iDXj
    # NAV[ep][1:2] .*= 0.7
    # println("$(iDXj[1]), $(iDXj[2]), $(iDXj[3])")

    odoDict[ep] = Pose2Pose2(MvNormal(NAV[ep], odonoise) )
  end
  rangepts = rangedata[ep][:]
  rangeprob = kde!(rangepts)
  azipts = azidata[ep][:,1]
  aziprob = kde!(azipts)

  # alternative range probability
  rawmf = readdlm("/home/dehann/data/sandshark/full_wombat_2018_07_09/extracted/matchedfilter/raw/$(ep).csv",',')
  range_bss = AliasingScalarSampler(rawmf[:,1], exp.(rawmf[:,2]), SNRfloor=0.6)

  # prep the factor functions
  ppbrDict[ep] = Pose2Point2BearingRange(aziprob, range_bss) # rangeprob
  lastepoch = ep
end

function runEpochs!(fgl, STEP::Int, index::Vector{Int})
    for ep in epochs[(STEP+1):(STEP+10)]
      curvar = Symbol("x$(index[1])")
      addVariable!(fgl, curvar, Pose2)

      # xi -> l1 - nonparametric factor
      if index[1] % 5 == 0
          addFactor!(fgl, [curvar; :l1], ppbrDict[ep], autoinit=true)
      end

      if ep != epochs[1]
        # Odo factor x(i-1) -> xi
        addFactor!(fgl, [Symbol("x$(index[1]-1)"); curvar], odoDict[ep], autoinit=true)
      else
        # Prior to the first pose location (a "GPS" prior)
        initLoc = [interp_x(ep);interp_y(ep);interp_yaw(ep)]
        println("Adding a prior at $curvar, $initLoc")
        addFactor!(fgl, [curvar], PriorPose2( MvNormal(initLoc, Matrix(Diagonal([0.1;0.1;0.05].^2))) ), autoinit=true)
      end
      # Heading partial prior
      addFactor!(fgl, [curvar], RoME.PartialPriorYawPose2(Normal(interp_yaw(ep), deg2rad(3))), autoinit=true)
      index[1] += 1
    end
    nothing
end

function deleteEpochs!(fgl, lbls::Vector{Symbol})
  function delFactsVar!(fgl, lb::Symbol)
    fcts = ls(fgl, lb)
    map(x->deleteFactor!(fgl, x), fcts)
    deleteVariable!(fgl, lb)
    nothing
  end
  map(x->delFactsVar!(fgl, x), lbls)
  nothing
end

## Step: Building the factor graph
fg1 = initfg()
# Add a central beacon with a prior
addVariable!(fg1, :l1, Point2)
# Pinger location is (0.6; -16)
addFactor!(fg1, [:l1], PriorPose2( MvNormal([17; 1.8], Matrix(Diagonal([0.1; 0.1].^2)) ) ), autoinit=true) # [0.6; -16]

# init tree for simpler code later down
tree1 = wipeBuildNewTree!(fg1)

## and some settings

# enable fixed lag operation
# getSolverParams(fg1).isfixedlag = true
# getSolverParams(fg1).qfl = 20

# first solve and initialization
getSolverParams(fg1).drawtree = true
# getSolverParams(fg).showtree = false

mkpath(joinpath(getSolverParams(fg1).logpath, "graphs"))

index1 = Int[0;]
index2 = Int[0;]
storeLast = Dict{Symbol,BallTreeDensity}()

for STEP in 0:10:80
    global fg1, tree1
    global index1, index2
    global storeLast

    runEpochs!(fg1, STEP, index1)
    poses = sortVarNested(ls(fg1, r"x"))

    if STEP-10 >= 0
        @show poses[1], poses[end], poses[end-10]
        @show dellist = poses[1:(end-10-1)]
        deleteEpochs!(fg1, dellist)
        addFactor!(fg1, [poses[end-10];], PriorPose2(storeLast[poses[end-10]]))
    end

    saveDFG(fg1, joinpath(getSolverParams(fg1).logpath, "graphs", "fg_$(STEP)_before.pdf"))
    @time tree1, smt, hist = solveTree!(fg1, tree1)
    saveDFG(fg1, joinpath(getSolverParams(fg1).logpath, "graphs", "fg_$(STEP)_after.pdf"))

    pla = drawPosesLandmarksAndOdo(fg1, ppbrDict, navkeys, X, Y, lblX, lblY)
    plb = plotSandsharkFromDFG(fg1)
    hstack(pla,plb) |> PDF(joinpath(getSolverParams(fg1).logpath, "sandshark-beacon_$STEP.pdf"))

    poses = sortVarNested(ls(fg1, r"x"))
    storeLast[poses[end]] = getKDE(fg1, poses[end])
    @show keys(storeLast)
end

## debugging a plotting

## plot forward convolves

fcts = setdiff(ls(fg1, :l1), lsfPriors(fg1))

pts = approxConv(fg1, fcts[1], :l1)

pts, infd = predictbelief(fg1,:l1,fcts)

pl1 = Gadfly.plot(x=pts[1,:],y=pts[2,:],Geom.hexbin)

union!(PL.layers, pl1.layers)

PL

plotKDE(manikde!(pts, Point2().manifolds))

dfg = fg1
fctsym = :x70l1f1
fct = getFactorType(dfg, fctsym)

#

fctsym

function solveFactorMeasurement(dfg::AbstractDFG,
                                fctsym::Symbol  )
  #

end

fcto = getFactor(dfg, fctsym)
varsyms = fcto._variableOrderSymbols
vars = map(x->getPoints(getKDE(dfg,x)), varsyms)
fcttype = getFactorType(fcto)
zDim = getData(fcto).fnc.zDim

N = size(vars[1])[2]
res = zeros(zDim)
ud = FactorMetadata()
meas = (zeros(zDim, N),)

idx = 1

function makemeas(i, dm, meas)
  meas[1][:,i] = meas[1][:,i] + dm
  return meas
end

ggo = (i, dm) -> fcttype(res,ud,i,makemeas(i, dm, meas),vars...)

ggo(1, [0.0;0.0])

using Optim

optimize((x) -> ggo(1,x), [0.0; 0.0])

Gadfly.plot(z=(x,y)->ggo(1,[x;y]), xmin=[-pi],xmax=[pi],ymin=[-100.0],ymax=[100.0], Geom.contour)

# idea is to find predicted parameters for noise model

#

import RoMEPlotting: plotFactor

function plotFactor(dfg::AbstractDFG, fctsym::Symbol, fct::Pose2Point2Bearing; hdl=[])
  #

  # variables
  vars = ls(dfg, fctsym)

  # the pose
  pose = intersect(vars, ls(dfg, Pose2))[1]
  poin = intersect(vars, ls(dfg, Point2))[1]

  # convolve the yaw angle with bearing rotation model
  pX = marginal(getKDE(dfg, pose), [3])
  pts = approxConvCircular(pX, fct.bearing)

  # draw plots
  measest = manikde!(pts, Sphere1)

  # inverse solve for predicted bearing
  dx = getPoints(getKDE(dfg, poin))[1,:] - getPoints(getKDE(dfg, pose))[1,:]
  dy = getPoints(getKDE(dfg, poin))[2,:] - getPoints(getKDE(dfg, pose))[2,:]
  pred = reshape(atan.(dy,dx), 1,:)

  ppX = manikde!(pred, Sphere1)

  plcl = plotKDECircular( [measest; ppX], logpdf=true, legend=["Meas. Est.";"Predicted"] )

  # plot pose and point by itself
  posepl1 = plotKDE(dfg, pose, dims=[1;2], c=["green"])
  posepl2 = plotKDECircular(marginal(getKDE(dfg, pose), [3]))
  landmpl = plotKDE(dfg, poin, c=["red"])

  tfg = initfg()
  addVariable!(tfg, pose, Pose2)
  addVariable!(tfg, poin, Point2)
  addFactor!(tfg, [pose;poin], fct, autoinit=false)
  manualinit!(tfg, pose, getKDE(dfg,pose))
  manualinit!(tfg, poin, getKDE(dfg,poin))
  plt = drawPosesLandms(tfg, point_size=5pt)

  # plot handles
  push!(hdl, landmpl)
  push!(hdl, posepl1)
  push!(hdl, posepl2)
  push!(hdl, plcl)
  push!(hdl, plt)

  hstack(vstack(landmpl,posepl1, posepl2), vstack(plcl, plt))
end

hdl = []

function plotFactor(dfg::AbstractDFG, fctsym::Sybmol, fct::Pose2Point2BearingRange; hdl = [])
  #

# variables
vars = ls(dfg, fctsym)

# the pose
pose = intersect(vars, ls(dfg, Pose2))[1]
poin = intersect(vars, ls(dfg, Point2))[1]

# plot current pose & point
pl_poin = plotKDE(dfg, poin, levels=5, c=["red";])
pl_pose = plotPose(dfg, pose, c=["black"])

# project landmark

# do range model separately too
pr = Pose2Point2Range(fct.range)

plotFactor(dfg, fctsym, pr, hdl=hdl)

hdl

end

plhist

## development

br = Pose2Point2Bearing(fct.bearing)

plotFactor(dfg, :x70l1f1, br)

pr = Pose2Point2Range(fct.range)

function plotFactor(dfg::G, fctsym::Symbol, fct::Pose2Point2Bearing) where {G <: AbstractDFG}

hdl = []

# variables
vars = ls(dfg, fctsym)

# the pose
pose = intersect(vars, ls(dfg, Pose2))[1]
poin = intersect(vars, ls(dfg, Point2))[1]

poseyaw = plotKDECircular(marginal(getKDE(dfg, pose), [3]), title="yaw angle")
bearpl = plotKDECircular(fct.bearing, title="bearing body frame")

vstack(poseyaw, bearpl)

end

"""
    $SIGNATURES

Build an approximate density `[Y|X,DX,.]=[X|Y,DX][DX|.]` as proposed by the conditional convolution.

Notes
- Assume both are on circular manifold, `manikde!(pts, (:Circular,))`
"""
function approxConvCircular(pX::BallTreeDensity, pDX::BallTreeDensity)
  #

  # building basic factor graph
  tfg = initfg()
  addVariable!(tfg, :s1, Sphere1)
  addVariable!(tfg, :s2, Sphere1)
  addFactor!(tfg, [:s1;:s2], Sphere1Sphere1(pDX), autoinit=false)
  manualinit!(tfg,:s1, pX)

  # solve for outgoing proposal value
  approxConv(tfg,:s1s2f1,:s2)
end

# """
#     $SIGNATURES
#
# Solve for the measurement values, leaving all variables in tact.
# """
# function inverseEvaluateFactor(dfg::G, fct::F) where {G <: AbstractDFG, F <: DFGFactor}
#
#
# end

## overview

PL = drawPosesLandmarksAndOdo(fg1, ppbrDict, navkeys, X, Y, lblX, lblY)

plotSandsharkFromDFG(fg1, scale=1.5)

## slightly random

drawGraph(fg1)

plotLocalProduct(fg1, :x55)

sortVarNested(ls(fg1, r"x"))

pp, parr, partials, lb, infdim = IIF.localProduct(fg1, :x55)

plotKDE([pp;parr], levels=1, legend=["prod"; string.(lb[1:3])], dims=[1;2])

plotPose(fg1, :x59)

## DEBUG beam pattern timestamps

epochs

pll = layerBeamPatternRose(ppbrDict[ep].bearing, wRr=wRr, wTRr=theta[1:2], scale=5.0)

ROT = Float64[]

for (tim, data) in azidata
  global ROT
  ROT = [ROT; data[:,1]]
end

Gadfly.plot(x=ROT, Geom.histogram)

st = sort(collect(keys(azidata)))[150]

Gadfly.plot(x=azidata[st][:,1], Geom.histogram)

Gadfly.plot(x=getPoints(PX)[:], Geom.histogram)

mask = rand(Categorical(0.001*ones(1000)),100)
pts = collect(azidata[st][mask,1:1]')
PX = manikde!( pts, (:Circular,))
AMP.plotKDECircular([PX;], rVo=[10.0; 0.0; 0.0])

Gadfly.plot(
  layerBeamPatternRose( PX, wTRr= , wRr=wRr )
)

## Plot just dead reckoning

ep = epochs[1]
initLoc = [interp_x(ep);interp_y(ep);interp_yaw(ep)]

#

XX= Float64[]
YY = Float64[]
TH = Float64[]

PL = []

for i in 1:2:80
  @show ep = epochs[i]

  push!(XX, interp_x(ep))
  push!(YY, interp_y(ep))
  push!(TH, interp_yaw(ep)) # bad at wrap boundary

  mask = rand(Categorical(0.001*ones(1000)),100)
  pts = collect(azidata[ep][mask,1:1]')
  PX = manikde!( TU.wrapRad.(deepcopy(pts)), (:Circular,))
  # PX = manikde!( pts, (:Euclid,))

  pl = plotSandsharkFromDFG(dfg)
  # pll = AMP.plotKDECircular([PX;], rVo=[XX[end]; YY[end]; TH[end]], radix=1.2, scale=0.3)
  # push!(PL, pll)
  # his1 = Gadfly.plot(x=azidata[ep][:,1], Geom.histogram)
  # his2 = Gadfly.plot(x=pts, Geom.histogram)
  # vstack(his1, his2, pll) |> PDF("/tmp/caesar/random/debug/$ep.pdf")
end

pl = Gadfly.plot(
  Gadfly.layer(x=XX, y=YY, Geom.path),
  Gadfly.layer(x=[0.6;], y=[-16.0;], Geom.point, Theme(default_color=colorant"red")),
)

for pll in PL
  union!(pl.layers, pll.layers)
end

pl

## plot raw data in loop

i = 0
for t in timestamps
  @show t
  global i += 1
  pl1 = Gadfly.plot(x=azidata[t][:,1], y=azidata[t][:,2], Geom.point)
  pl2 = Gadfly.plot(x=)
   # |> PNG("/tmp/plots/$i.png")

end

## plot all beam patterns from factor graph

dfg = fg1

plotSandsharkFromDFG(dfg)

ls(dfg, Point2)
ls(dfg, Pose2)

ls(dfg, Pose2Pose2)

## RANGE MODELING

fct_range = getFactorType(fg1, :x45l1f1).range

plotKDE(fct_range)

rawmf = readdlm("/home/dehann/data/sandshark/full_wombat_2018_07_09/extracted/matchedfilter/raw/$(epochs[45]).csv",',')
bss = AliasingScalarSampler(rawmf[:,1], rawmf[:,2], SNRfloor=0.6)

pts = rand(bss,100)
Gadfly.plot(x=rawmf[:,1], y=rawmf[:,2], Geom.path)

Gadfly.plot(x=pts, Geom.histogram)

pts = rand(fct_range, 1000)

#


text:
struct Prior
    N::DiscreteUnivariateDistribution
    v::ContinuousUnivariateDistribution
    theta_penalty::Function

    ## known parameters
    Q         ::Int ## order of spline
    N_max     ::Int ## maximum number of knots
    cand_knots::Vector{Float64} ## candidate knot locations
    outer     ::NTuple{2,Float64} ## outer knot locations
end

Prior() = Prior(Poisson(4), InverseGamma(0.01, 0.01), (theta, v) -> sum(theta.^2)/(2*v), 4, 100, Vector(0.01:0.01:0.99), (0.0, 1.0))

mutable struct Sampler
    loglik::Function
    prior ::Prior
    tuner ::Tuner

    ## state variables
    t      ::Int
    warmup ::Int
    attempt::Vector{Int}
    accept ::Vector{Int}

    ## posterior samples of model parameters
    theta::Array{Float64,2} ## log control points
    v    ::Vector{Float64} ## variance of prior for theta (smoothness penalty parameter)
    N    ::Vector{Int} ## number of internal knots
    knots::Array{Int,2} ## inner knot locations

    knot_ocp::Vector{Bool}  ## knot occupation indicators (same length as prior.cand_knots)

    ## when number and location of knots are fixed then only the weights change
    ## matrices of the basis spline functions evaluated at each event time and its integrals can pre-computed once
    ## and stored for substantial speed-up
    fixed_knots   ::Bool
    precalc_haz   ::Array{Float64,2}
    precalc_cumhaz::Array{Float64,2}
end

mutable struct Param
    N::Int
    theta::Vector{Float64}
    knots::Vector{Int}
    v::Float64
end

"""
    create_sampler(loglik::Function, prior::Prior=Prior())

Create `Sampler` object from log-likelihood `loglik` and prior object `prior`.
The number of iterations of the MCMC and the initial state must be specified either in a subsequent call to `set_initial_state!` or
in the call to `sample!`.
"""
function create_sampler(loglik::Function, prior::Prior=Prior())
    Sampler(loglik, prior, Tuner(0), 0, 0, Vector{Int}(undef, 0), Vector{Int}(undef,0),
            Array{Float64,2}(undef,0,0), Vector{Float64}(undef,0), Vector{Int}(undef,0), Array{Float64,2}(undef,0,0),
            Vector{Bool}(undef,0), false, Array{Float64,2}(undef,0,0), Array{Float64,2}(undef,0,0))
end

"""
    setup_sampler(M::Int, evtime::Vector{Float64}, loglik::Function; nknots=4, N_max=100)

Create `Sampler` object from log-likelihood `loglik` and initialize with default prior and by default 4 randomly selected knots at quantiles of observed event times 'evtime'. The sampler is setup for 'M' iterations and initial state is set to default values.
"""
function setup_sampler(M::Int, evtime::Vector{Float64}, loglik::Function; nknots=4, N_max=100, outer=(0.0, maximum(evtime)))
    cand_knots = quantile(evtime, Vector(0.0:1/(N_max+1):1.0))[2:(end-1)]

    @assert minimum(cand_knots) > outer[1]
    @assert maximum(cand_knots) < outer[2]

    prior = Prior(Poisson(nknots), InverseGamma(0.01, 0.01),
                  (theta, v) -> sum(theta.^2)/(2*v),
                  nknots, N_max, cand_knots, outer)

    s = create_sampler(loglik, prior)

    knots0 = sort(sample(1:length(cand_knots), nknots))

    set_initial_state!(s, M, Param(nknots, zeros(Float64, nknots+4), knots0, 1.0))
    s.tuner = Tuner(Int(ceil(M/2))) ##set_tuner!(s, Tuner(Int(M/2)))
    s
end

"""
    set_initial_state!(s::Sampler, M::Int, init::Param; warmup=Int(floor(M/2)), fixed_knots::Bool=false,
                                time::Vector{Float64}=Vector{Float64}(undef,0), status::Vector{Bool}=Vector{Bool}(undef,0))

Set the initial state of the Markov chain to `init` and allocate arrays holding the posterior samples of `M` total iterations.
When number and location of the knots are fixed (`fixed_knots == true`) then the vector of event times `time` and the censoring
indicators `status`
"""
function set_initial_state!(s::Sampler, M::Int, init::Param; warmup=Int(floor(M/2)), fixed_knots::Bool=false,
                            time::Vector{Float64}=Vector{Float64}(undef,0), status::Vector{Bool}=Vector{Bool}(undef,0))
    ## number of inner knots
    s.N = zeros(Int, M+1)
    s.N[1] = length(init.knots)

    ## candidate knots
    ##s.cand_knots = s.prior.cand_knots ##quantile(data.time[data.status], Array(0.0:1/(N_max+1):1.0))[2:(end-1)]
    n_candk = length(s.prior.cand_knots)

    ## inner knots
    s.knots = Array{Int,2}(undef, M+1, s.prior.N_max)
    s.knots[1,1:s.N[1]] = init.knots

    ## knot occupation indicators
    s.knot_ocp = zeros(Int, n_candk)
    ##ck = sort(sample(1:n_candk, N[1])) ##sort(randperm(n_candk)[1:N[1]])
    s.knot_ocp[init.knots] .= true

    Q = s.prior.Q

    ## log basis spline weights
    s.theta = Array{Float64,2}(undef, M+1, s.prior.N_max + Q)
    @assert length(init.theta) == s.N[1] + Q
    s.theta[1,1:(s.N[1]+Q)] = init.theta

    if fixed_knots
        ## we need the baseline hazard basis function at all failure times, and their
        ## integrals (cumulative hazard basis functions) at every event time
        s.precalc_haz, s.precalc_cumhaz = Spline.eval_basis(time[status], time, init.knots, s.prior.outer, Q)
    end

    ## prior variance of log basis spline weights
    s.v = Vector{Float64}(undef,M+1)
    s.v[1] = init.v

    s.warmup = warmup
    s.attempt = zeros(Int, 5)
    s.accept = zeros(Int, 5)
    s.t = 1
end

function set_tuner!(s::Sampler, tuner::Tuner)
    s.tuner = tuner
end

"""
    sample_prior(prior::Prior)

Sample spline from 'prior' distribution
"""
function sample_prior(prior::Prior)
    outer = prior.outer
    cand_knots = prior.cand_knots

    N = min(rand(prior.N, 1)[1], length(cand_knots))
    K = N + prior.Q
    D = diagm(0 => ones(Float64, K))
    v = rand(InverseGamma(0.01 + K/2, 25), 1)[1]
    knots = sort(sample(cand_knots, N; replace=false))
    theta = rand(MvNormal(D ./ v), 1)[:,1]

    spline(outer, knots, exp.(theta))
end

"""
    sample!(s::Sampler, M::Int, init::Param; warmup=Int(floor(M/2)), tuner=Tuner(warmup))

Perform a total of `M` iterations of the RJ-MCMC starting from initial state `init`. By default half of the total number of iterations are warmup
iterations where the standard deviation of the proposal distribution for the spline weights is tuned.
"""
function sample!(s::Sampler, M::Int, init::Param=extract(s, 1); warmup=Int(floor(M/2)), tuner=Tuner(warmup))
    set_initial_state!(s, M, init; warmup=warmup)
    set_tuner!(s, tuner)

    for t in 1:M
        update!(s, nothing)
    end
end

"""
    extract(s::Sampler, t::Int)

Return posterior sample of parameters at iteration `t`.
"""
function extract(s::Sampler, t::Int)
    Param(s.N[t], s.theta[t,1:(s.N[t]+s.prior.Q)], s.knots[t,1:s.N[t]], s.v[t])
end

"""
    extract_spline(s::Sampler, t::Int)

Returns the posterior spline object of type `Dierckx.Spline1D` of iteration `t`.
"""
function extract_spline(s::Sampler, t::Int)
    p = extract(s, t)
    Spline.spline(s.prior.outer, s.prior.cand_knots[p.knots], exp.(p.theta), s.prior.Q-1)
end

"""
    extract_hazard(s::Sampler, time::Vector{Float64})

Evaluate the posterior hazard function sample of iteration `t` at each element of `time`.
"""
function extract_hazard(s::Sampler, time::Vector{Float64}, t::Int)
    p = extract_spline(s, t)
    Spline.hazard(time, p)
end

"""
    extract_hazard(s::Sampler, time::Vector{Float64}; return_warmup=false)

Evaluate all posterior hazard function samples at each element of `time`. Ignores warmup iterations by default.
"""
function extract_hazard(s::Sampler, time::Vector{Float64}; return_warmup=false)
    start = 1 + s.warmup*(1 - return_warmup)
    res = Array{Float64,2}(undef, s.t-start+1, length(time))
    for t in start:s.t
        res[t-start+1,:] = extract_hazard(s, time, t)
    end
    res
end

"""
    extract_cumhaz(s::Sampler, time::Vector{Float64})

Evaluate the posterior cumulative hazard function sample of iteration `t` at each element of `time`.
"""
function extract_cumhaz(s::Sampler, time::Vector{Float64}, t::Int)
    p = extract_spline(s, t)
    Spline.cumhaz(time, p)
end

"""
    extract_cumhaz(s::Sampler, time::Vector{Float64}; return_warmup=false)

Evaluate all posterior cumulative hazard function samples at each element of `time`. Ignores warmup iterations by default.
"""
function extract_cumhaz(s::Sampler, time::Vector{Float64}; return_warmup=false)
    start = 1 + s.warmup*(1 - return_warmup)
    res = Array{Float64,2}(undef, s.t-start+1, length(time))
    for t in start:s.t
        res[t-start+1,:] = extract_cumhaz(s, time, t)
    end
    res
end

"""Evaluate posterior survival function sample of iteration `t` at each element of `time`."""
extract_survival(s::Sampler, time::Vector{Float64}, t::Int) = exp.(.-extract_cumhaz(s, time, t))

"""Evaluate all posterior survival function samples at each element of `time`. Ignores warmup iterations by default."""
extract_survival(s::Sampler, time::Vector{Float64}; return_warmup=false) = exp.(.-extract_cumhaz(s, time; return_warmup=return_warmup))

"""
    summary(s::Sampler, time::Vector{Float64})

Create `DataFrame` with mean posterior hazard, cumulative hazard and survival functions evaluated at each element of `time` and the corresponding
(pointwise) 95% credible intervals.
"""
function summary(s::Sampler, time::Vector{Float64})
    haz = extract_hazard(s, time)
    cumhaz = extract_cumhaz(s, time)
    surv = exp.(.-cumhaz)

    function msci(x)
        mean(x, dims=1)[1,:], sqrt.(var(x, dims=1))[1,:], [quantile(x[:,k], 0.025) for k in 1:size(x)[2]], [quantile(x[:,k], 0.975) for k in 1:size(x)[2]]
    end

    msci_haz = msci(haz)
    msci_cumhaz = msci(cumhaz)
    msci_surv = msci(surv)

    DataFrame(time=time,
              hazard=msci_haz[1], hazard_sd=msci_haz[2], hazard_lo=msci_haz[3], hazard_hi=msci_haz[4],
              cumhazard=msci_cumhaz[1], cumhazard_sd=msci_cumhaz[2], cumhazard_lo=msci_cumhaz[3], cumhazard_hi=msci_cumhaz[4],
              survival=msci_surv[1], survival_sd=msci_surv[2], survival_lo=msci_surv[3], survival_hi=msci_surv[4])

end

ext: jl


text:
## multi-channel.jl
#
## process multi-channel (mtch) data.
## these functions are not currently used

## Top-level function:
## run qPCR functions on multi-channel fluorescence data and output consolidated results.
## original matrix will become 3-D array, original string or list will be come list.
function process_mtch(
    input_dict_bych::OrderedDict, ## keyed by channels
    arydims2to3::Bool,
    func::Function,
    func_args...
    )

    channels = get_ordered_keys(input_dict_bych)

    result_dict_bych = OrderedDict([
        channel => func(input_dict_bych[channel], func_args...)
        for channel in channels])

    result_keys = get_ordered_keys(result_dict_bych[channels[1]])

    result_consoli = OrderedDict([
        result_key => consolidate_rdc_per_element(result_key, result_dict_bych, arydims2to3)
        for result_key in result_keys])

    return OrderedDict(
        "pre_consoli" => result_dict_bych,
        "post_consoli" => result_consoli)
end

## functions called by `process_mtch`

## function: consolidate for one result element from organization of channel then result element (as in `result_dict_bych`) into organization of result element then channel
function consolidate_rdc_per_element(
    result_key,
    result_dict_bych::OrderedDict,
    arydims2to3::Bool
    )

    channels = get_ordered_keys(result_dict_bych)
    erfc = result_dict_bych[channels[1]][result_key] ## element of interest in the result for the first channel

    if arydims2to3 && isa(erfc, AbstractArray) && ndims(erfc) == 2
        consolidated_ele = cat(3, map(channels) do channel
            result_dict_bych[channel][result_key]
        end...) ## The resulting DataArray is indexed by the linear indice of channel in `channels` instead of the actual channel values.
    else
        consolidated_ele = OrderedDict(map(channels) do channel
            channel => result_dict_bych[channel][result_key]
        end)
    end

    return consolidated_ele
end

#

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
module TestOperator
using MXNet
using Base.Test

using ..Main: rand_dims, reldiff

function test_scalar_op()
  data  = mx.Variable(:data)
  shape = rand_dims()
  info("Operator::scalar_op::dims = $shape")

  data_jl  = 5ones(Float32, shape)
  arr_data = mx.copy(data_jl, mx.cpu())
  arr_grad = mx.zeros(shape)

  test = 2 ./ (4 - ((1+data+1)*2/5) - 0.2)
  exec_test = mx.bind(test, mx.cpu(), [arr_data], args_grad=[arr_grad])
  mx.forward(exec_test)
  out = copy(exec_test.outputs[1])
  jl_out1 = (4 - ((1+data_jl+1)*2/5) - 0.2)
  jl_out = 2 ./ jl_out1
  @test reldiff(copy(out), jl_out) < 1e-6

  out_grad = 2mx.ones(shape)
  jl_grad  = 2copy(out_grad) / 5
  jl_grad  = 2jl_grad ./ (jl_out1 .^ 2)
  mx.backward(exec_test, out_grad)
  @test reldiff(copy(arr_grad), jl_grad) < 1e-6
end

################################################################################
# Run tests
################################################################################
test_scalar_op()

end


text:
using MarkableIntegers
using Test

pos1, pos2 = Int16( 1),  2
neg1, neg2 = Int16(-1), -2

pos1_nomark, pos2_nomark = Unmarked(pos1), Unmarked(pos2)
neg1_nomark, neg2_nomark = Unmarked(neg1), Unmarked(neg2)

pos1_marked, pos2_marked = Marked(pos1), Marked(pos2_nomark)
neg1_marked, neg2_marked = Marked(neg1), Marked(neg2_nomark)

@test ismarked(pos1)        == false
@test ismarked(pos1_nomark) == false
@test ismarked(pos1_marked) == true

@test pos2        == pos2_nomark
@test pos2        == pos2_marked
@test pos2_nomark == pos2_marked

@test pos2        == pos2_nomark
@test pos2        == pos2_marked
@test pos2_nomark == pos2_marked

@test isless(pos1, pos2_nomark)
@test isless(pos1, pos2_marked)
@test isless(pos1_nomark, pos2_nomark)
@test isless(pos1_nomark, pos2_marked)
@test isless(pos1_marked, pos2_nomark)
@test isless(pos1_marked, pos2_marked)

@test isless(neg2_nomark, neg1)
@test isless(neg2_marked, neg1)
@test isless(neg2_nomark, neg1_nomark)
@test isless(neg2_marked, neg1_nomark)
@test isless(neg2_nomark, neg1_marked)
@test isless(neg2_marked, neg1_marked)

@test isequal(pos2       , pos2_nomark)
@test isequal(pos2       , pos2_marked)
@test isequal(pos2_nomark, pos2_marked)
@test isequal(pos2_nomark, pos2)
@test isequal(pos2_marked, pos2)
@test isequal(pos2_marked, pos2_marked)

@test !isless(pos2       , pos2_nomark)
@test !isless(pos2       , pos2_marked)
@test !isless(pos2_nomark, pos2_marked)
@test !isless(pos2_nomark, pos2)
@test !isless(pos2_marked, pos2)
@test !isless(pos2_marked, pos2_marked)

pos3_nomark = pos1_nomark + pos2_nomark
pos3_marked = pos1_nomark + pos2_marked
@test !ismarked(pos3_nomark)
@test ismarked(pos3_marked)
@test pos3_nomark == pos3_marked
@test pos3_nomark == 3

@test pos2_nomark * neg2_nomark == -4
tst = div(pos3_nomark, pos2_marked)
@test ismarked(tst)
@test tst == div(3,2)
tst = fld(pos3_nomark, pos2_marked)
@test ismarked(tst)
@test tst == fld(3,2)
tst = cld(pos3_nomark, pos2_nomark)
@test isunmarked(tst)
@test tst == cld(3,2)
tst = mod(pos3_nomark, pos2_marked)
@test ismarked(tst)
@test tst == mod(3,2)
tst = rem(pos3_nomark, pos2_nomark)
@test isunmarked(tst)
@test tst == rem(3,2)

tst = pos3_nomark | 8
@test typeof(tst) == typeof(pos3_nomark)
@test tst == 3|8


text:
export ANYTHING, anything_of_type, anything_in, anything_containing, Matcher, redescribe

import Base: ==, &, |, !

immutable Matcher
  predicate::Function
  description::String
end

==(matcher1::Matcher, matcher2::Matcher) = error("Can not compare 2 Matchers")
==(value::WeakRef, matcher::Matcher) = matcher == value.value
==(value::Any, matcher::Matcher) = matcher == value
==(matcher::Matcher, value::WeakRef) = matcher.predicate(value.value)
==(matcher::Matcher, value::Any) = matcher.predicate(value)

Base.show(io::IO, matcher::Matcher) = print(io, matcher.description)

(&)(matcher1::Matcher, matcher2::Matcher) = Matcher(x->(matcher1.predicate(x) && matcher2.predicate(x)), "($(matcher1.description) && $(matcher2.description))")

(|)(matcher1::Matcher, matcher2::Matcher) = Matcher(x->(matcher1.predicate(x) || matcher2.predicate(x)), "($(matcher1.description) || $(matcher2.description))")

(!)(matcher::Matcher) = Matcher(x->(!matcher.predicate(x)), "(not $(matcher.description))")

redescribe(matcher::Matcher, description::String) = Matcher(matcher.predicate, description)

anything_of_type(T::Type) = Matcher(v -> isa(v, T), "anything_of_type($T)")
anything_in(value::WeakRef) = anything_in(value.value)
anything_in(value::Any) = Matcher(v -> v in value, "anything_in($value)")
anything_containing(value::WeakRef) = anything_containing(value.value)
anything_containing(value::Any) = Matcher(v -> value in v, "anything_containing($value)")

const ANYTHING = anything_of_type(Any)


text:
_ones_feminine = Dict{Int, NTuple{1, String}}(
    1 => ("одна",),
    2 => ("две",),
    3 => ("три",),
    4 => ("четыре",),
    5 => ("пять",),
    6 => ("шесть",),
    7 => ("семь",),
    8 => ("восемь",),
    9 => ("девять",),
)

_ones = Dict{Int, NTuple{1, String}}(
    1 => ("один",),
    2 => ("два",),
    3 => ("три",),
    4 => ("четыре",),
    5 => ("пять",),
    6 => ("шесть",),
    7 => ("семь",),
    8 => ("восемь",),
    9 => ("девять",),
)

_tens = Dict{Int, NTuple{1, String}}(
    0 => ("десять",),
    1 => ("одиннадцать",),
    2 => ("двенадцать",),
    3 => ("тринадцать",),
    4 => ("четырнадцать",),
    5 => ("пятнадцать",),
    6 => ("шестнадцать",),
    7 => ("семнадцать",),
    8 => ("восемнадцать",),
    9 => ("девятнадцать",),
)

_twenties = Dict{Int, NTuple{1, String}}(
    2 => ("двадцать",),
    3 => ("тридцать",),
    4 => ("сорок",),
    5 => ("пятьдесят",),
    6 => ("шестьдесят",),
    7 => ("семьдесят",),
    8 => ("восемьдесят",),
    9 => ("девяносто",),
)

_hundreds = Dict{Int, NTuple{1, String}}(
    1 => ("сто",),
    2 => ("двести",),
    3 => ("триста",),
    4 => ("четыреста",),
    5 => ("пятьсот",),
    6 => ("шестьсот",),
    7 => ("семьсот",),
    8 => ("восемьсот",),
    9 => ("девятьсот",),
)

_thousands = Dict{Int, NTuple{3, String}}(
    1 => ("тысяча", "тысячи", "тысяч"), # 10^3
    2 => ("миллион", "миллиона", "миллионов"), # 10^6
    3 => ("миллиард", "миллиарда", "миллиардов"), # 10^9
    4 => ("триллион", "триллиона", "триллионов"), # 10^12
    5 => ("квадриллион", "квадриллиона", "квадриллионов"), # 10^15
    6 => ("квинтиллион", "квинтиллиона", "квинтиллионов"), # 10^18
    7 => ("секстиллион", "секстиллиона", "секстиллионов"), # 10^21
    8 => ("септиллион", "септиллиона", "септиллионов"), # 10^24
    9 => ("октиллион", "октиллиона", "октиллионов"), # 10^27
    10 => ("нониллион", "нониллиона", "нониллионов"), # 10^30
)


text:
module YaoGraphs

include("register.jl")

end # module

max_stars_repo_name: QuantumBFS/YaoGraphs

max_issues_repo_name: QuantumBFS/YaoGraphs

max_forks_repo_name: QuantumBFS/YaoGraphs


text:
export ft_setlatencytimer

function ft_setlatencytimer(ft_handle::Culong, timer::Integer)
  @assert timer> @linux? 0:1 # is 1ms valid for linux?
  @assert timer<256
  ft_status = ccall((:FT_SetLatencyTimer, d2xx),
                     Cuint,
                     (Culong, Cuchar),
                     ft_handle, timer)
  checkstatus(ft_status)
  return nothing
end


text:
export get_edge_plot, get_arrow_plot, get_node_plot, get_nlabel_plot, get_elabel_plot

"Get the `EdgePlot` subplot from a `GraphPlot`."
function get_edge_plot(gp::GraphPlot)
    p = gp.plots[1]
    @assert p isa EdgePlot
    return p
end

"Get the scatter plot of the arrow heads from a `GraphPlot`."
function get_arrow_plot(gp::GraphPlot)
    p = gp.plots[2]
    @assert p isa Scatter
    @assert p.marker[] == '➤'
    return p
end

"Get the scatter plot of the nodes from a `GraphPlot`."
function get_node_plot(gp::GraphPlot)
    p = gp.plots[3]
    @assert p isa Scatter
    @assert p[1][] == gp[:node_pos][]
    return p
end

"Get the text plot of the node labels from a `GraphPlot`."
get_nlabel_plot(gp::GraphPlot) = _get_label_plot(gp, gp.nlabels[])
"Get the text plot of the edge labels from a `GraphPlot`."
get_elabel_plot(gp::GraphPlot) = _get_label_plot(gp, gp.elabels[])

function _get_label_plot(gp::GraphPlot, labels)
    ps = filter(p -> p isa Makie.Text && p[1][] == labels, gp.plots)
    if isempty(ps)
        return nothing
    elseif length(ps) == 1
        return ps[1]
    else
        error("Could not determine plot $ps")
    end
end

"""
    getattr(o::Observable, idx, default=nothing)

If observable wraps an AbstractVector or AbstractDict return
the value at idx. If dict has no key idx rerturn default.
Else return the one and only element.
"""
function getattr(o::Observable, idx, default=nothing)
    if o[] isa AbstractVector && !isa(o[], Point)
        return o[][idx]
    elseif o[] isa AbstractDict
        return get(o[], idx, default)
    else
        return o[] === nothing ? default : o[]
    end
end

"""
    Pointf(p::Point{N, T})

Convert Point{N, T} or NTuple{N, T} to Point{N, Float32}.
"""
Pointf(p::Union{Point{N,T}, NTuple{N,T}}) where {N,T} = Point{N, Float32}(p)
Pointf(p::StaticVector{N, T}) where {N,T} = Point{N, Float32}(p)
Pointf(p::Vararg{T,N}) where {N,T} = Point{N, Float32}(p)
Pointf(p::Vector{T}) where {T} = Point{length(p), Float32}(p)

"""
    align_to_dir(align::Tuple{Symbol, Symbol})

Given a tuple of alignment (i.e. `(:left, :bottom)`) return a normalized
2d vector which points in the direction of the offset.
"""
function align_to_dir(align::Tuple{Symbol, Symbol})
    halign, valign = align

    x = 0.0
    if halign === :left
        x = 1.0
    elseif halign === :right
        x = -1.0
    end

    y = 0.0
    if valign === :top
        y = -1.0
    elseif valign === :bottom
        y = 1.0
    end
    norm = x==y==0.0 ? 1 : sqrt(x^2 + y^2)
    return Point2f(x/norm, y/norm)
end

"""
    plot_controlpoints!(ax::Axis, gp::GraphPlot)
    plot_controlpoints!(ax::Axis, path::BezierPath)

Add all the bezier controlpoints of graph plot or a single
path to the axis `ax`.
"""
function plot_controlpoints!(ax::Axis, gp::GraphPlot)
    ep = get_edge_plot(gp)
    ep.plots[1] isa BezierSegments || return
    ep = ep.plots[1]
    paths = ep[:paths][]

    for (i, p) in enumerate(paths)
        p isa Line && continue
        color = getattr(gp.edge_color, i)
        plot_controlpoints!(ax, p; color)
    end
end

function plot_controlpoints!(ax::Axis, p::BezierPath; color=:black)
    for (j, c) in enumerate(p.commands)
        if c isa CurveTo
            segs = [p.commands[j-1].p, c.c1, c.p, c.c2]
            linesegments!(ax, segs; color, linestyle=:dot)
            scatter!(ax, [c.c1, c.c2]; color)
        end
    end
end


text:
fname = path*"/SampleFiles/fdsn.conf"
sac_pz_file = path*"/SampleFiles/SAC/test_sac.pz"
hood_reg = Float64[44.8, 46.0, -122.4, -121.0]
rainier_rad = Float64[46.852886, -121.760374, 0.0, 0.1]
h5file = "sacreq.h5"
seisfile = "sacreq.seis"

printstyled("  FDSN web requests\n", color=:light_green)

safe_isfile(h5file) && safe_rm(h5file)

printstyled("    fdsn_chp\n", color=:light_green)

# fdsn_chp
c1 = "CC.VALT..,PB.B001..BS?,PB.B001..E??,XX.YYY.00.BHZ"
c2 = ["CC.VALT..", "PB.B001..BS?", "PB.B001..E??", "XX.YYY.00.BHZ"]
c3 = ["CC" "VALT" "" ""; "PB" "B001" "" "BS?"; "PB" "B001" "" "E??"; "XX" "YYY" "00" "BHZ"]
try
  @test (minreq(SeisIO.fdsn_chp(c1, 0)))[:,1:4] == SeisIO.fdsn_chp(c2, 0)[:,1:4] == minreq(SeisIO.fdsn_chp(c3, 0))
catch err
  @warn("Inconsistent reslts from fdsn_chp!")
end

# FDSNsta
printstyled("    FDSNsta\n", color=:light_green)
S = FDSNsta("CC.VALT..,PB.B001..BS?,PB.B001..E??")
for i in ("PB.B001.T0.BS1", "PB.B001..EHZ", "CC.VALT..BHZ")
  j = findid(S, i)
  try
    @test j > 0
  catch
    @warn(string("No data from ", i, "; check connection!"))
  end
end

# FDSNsta with MultiStageResp
S = FDSNsta("CC.VALT..,PB.B001..BS?,PB.B001..E??", msr=true)

# With autoname
printstyled("      get_data(\"FDSN\", ..., autoname=true)\n", color=:light_green)
req_f = "2019.001.00.00.00.000.UW.VLL..EHZ.R.mseed"
req_ok = (try
    S = check_get_data("FDSN", "UW.VLL..EHZ", src="IRIS", s="2019-01-01", t=3600, autoname=true)
    true
  catch
    @warn("Station VLL appears to be offline; test skipped.")
    false
  end)
if req_ok
  @test safe_isfile(req_f)
  rm(req_f)

  printstyled("        test match to IRIS filename convention\n", color=:light_green)
  S = check_get_data("IRIS", "UW.VLL..EHZ", s="2019-01-01", t=3600, autoname=true)
  @test safe_isfile(req_f)
end

printstyled("      radius search (rad=)\n", color=:light_green)
rad = Float64[45.373514, -121.695919, 0.0, 0.1]
S = FDSNsta(rad=rainier_rad)
try
  @test S.n > 0 # Test will break if everything around Mt. Rainier is offline
catch
  @warn("Stations around Mt. Rainier appear to be offline.")
end

printstyled("      rectangular search (reg=)\n", color=:light_green)
S = FDSNsta(reg=hood_reg)
try
  @test S.n > 0 # Test will break if everything around Mt. Hood is offline
catch
  @warn("Stations around Mt. Rainier appear to be offline.")
end
printstyled("    get_data\n", color=:light_green)

printstyled("      GeoCSV output\n", color=:light_green)
S = check_get_data("FDSN", "CC.JRO..BHZ,IU.COLA.00.*", src="IRIS", s=-600, t=0, fmt="geocsv", w=true)
S = check_get_data("FDSN", "CC.JRO..BHZ,CC.VALT.*", src="IRIS", s=-300, t=0, fmt="geocsv.slist")

printstyled("      config file for channel spec\n", color=:light_green)

S = SeisData()
check_get_data!(S, "FDSN", fname, src="IRIS", s=-600, t=0, w=true)
# Ensure we got data
if (isempty(S) ? 0 : maximum([length(x) for x in S.x])) == 0
  @warn("Request returned no data")
else
  printstyled("        writesac with channel range on this request\n", color=:light_green)
  if S.n > 3
    writesac(S, chans=1:2)
    for i in 1:2
      @test any([(occursin("writesac", n) && occursin(S.id[i], n)) for n in S.notes[i]])
    end
  end

  printstyled("        bandpass filtering of this request\n", color=:light_green)
  deleteat!(S, findall(S.fs.<25.0))
  filtfilt!(S, fl=0.01, fh=10.0)
  for i in 1:2
    @test any([occursin("processing", n) for n in S.notes[i]])
  end
end

# Ensure station headers are set
ids = ["UW.HOOD..ENE", "CC.VALT..BHZ", "UW.TDH..EHZ", "UW.VLL.EHZ"]
fss = [100.0, 50.0, 100.0, 100.0]
codes = ['N', 'H', 'H', 'H']
for i = 1:4
  j = findid(S, ids[i])
  if j > 0
    @test ≈(S.fs[j], fss[i])
    @test inst_code(S, j) == codes[i]
    break
  else
    @warn(string("No data from ", ids[i], "; check connection!"))
  end
end

# Check that headers get overwritten with SACPZ info when we use read_sacpz
Nc = S.n
read_sacpz!(S, sac_pz_file)
@test S.n > Nc
i = findid("CC.VALT..BHZ", S)
if i > 0
  @test S.misc[i]["OUTPUT UNIT"] == "COUNTS"
end
i = findid("UW.HOOD..ENE", S)
if i > 0
  @test S.misc[i]["INSTTYPE"] == "ES-T-3339=Q330S+-6410"
end

# Check that msr=true works
S = SeisData()
check_get_data!(S, "FDSN", fname, src="IRIS", msr=true, s=-600, t=0)
if isempty(S)
  @warn("Empty request; check connectivity!")
else
  for i in 1:S.n
    @test typeof(S.resp[i]) == MultiStageResp
  end
end

# Try a string array for input
printstyled("      string array for channel spec\n", color=:light_green)
S = SeisData()
check_get_data!(S, "FDSN", ["UW.HOOD..E??", "CC.VALT..???", "UW.XNXNX.99.QQQ"], src="IRIS", s=-600, t=0, opts="szsrecs=true")

# Try a single string
printstyled("      string for channel spec\n", color=:light_green)
S = check_get_data("FDSN", "CC.JRO..BHZ,IU.COLA.00.*", src="IRIS", s=-600, t=0, v=1,
  demean=true,
  detrend=true,
  rr=true,
  taper=true,
  ungap=true,
  unscale=true)

# This should return exactly 4 days of data, which we know IRIS' FDSN server has
printstyled("      multi-day request\n", color=:light_green)
ts = "2018-01-31T00:00:00"
te = "2018-02-02T00:00:00"
S = get_data("FDSN", "CI.ADO..BH?", s=ts, t=te)
id = "CI.ADO..BHE"
i = findid(S, id)
if i == 0
  @warn(string("No data from ", id, "; check connection!"))
elseif isempty(S)
  @warn(string("Request empty; rest of test skipped!"))
else
  # Check that we have two complete days of data with no gaps

  if (length(S.x[i]) / (86400*S.fs[i])) != 2.0
    @warn(string("Partial outage; missing data from ", id, "; check connection!"))
  end

  printstyled("        are data written identically?\n", color=:green)

  # get rid of bad request channels
  k = Int64[]
  append!(k, findall([startswith(id, "XX.FAIL") for id in S.id]))
  append!(k, findall([startswith(id, "XX.FMT") for id in S.id]))
  isempty(k) || deleteat!(S, k)

  # write ASDF first, since this modifies the first sample start time in S
  write_hdf5(h5file, S, add=true, ovr=true, len=Day(2))
  S1 = read_hdf5(h5file, ts, te, msr=false)

  for f in (:id, :loc, :fs, :gain, :resp, :units, :misc)
    @test isequal(getfield(S1,f), getfield(S,f))
  end
  for i in 1:S.n
    (isempty(S.x[i]) || isempty(S1.x[i])) && continue
    x1 = S1.x[i]
    x2 = S.x[i]
    if (any(isnan.(x1)) == false) && (any(isnan.(x2)) == false)
      @test x1 ≈ x2
    end
  end

  # Check that these data can be written and read faithfully in ASDF, SAC, and SeisIO
  writesac(S)
  wseis(seisfile, S)
  S2 = rseis(seisfile)[1]
  if !(isempty(S) || isempty(S2))
    for f in (:id, :loc, :fs, :gain, :resp, :units, :misc)
      @test isequal(getfield(S2,f), getfield(S,f))
    end
    for i in 1:S.n
      x1 = S.x[i]
      x2 = S2.x[i]
      if (any(isnan.(x1)) == false) && (any(isnan.(x2)) == false)
        @test x1 ≈ x2
      end
    end
  end

  # These are the only fields preserved; :loc is preserved to Float32 precision
  d0 = DateTime(ts)
  y0 = Year(d0).value
  j0 = md2j(y0, Month(d0).value, Day(d0).value)
  sac_str = string(y0, ".", lpad(j0, 3, '0'), "*CI.ADO..BH*SAC")
  S1 = read_data("sac", sac_str)
  for f in (:id, :fs, :gain, :t)
    @test getfield(S, f) == getfield(S1, f)
    @test getfield(S, f) == getfield(S2, f)
  end
  for f in (:lat, :lon, :el, :dep, :az, :inc)
    @test isapprox(getfield(S.loc[i], f), getfield(S1.loc[i], f), atol=1.0e-3)
  end
  for i in 1:S.n
    x = S.x[i]
    x1 = S1.x[i]
    x2 = S2.x[i]
    if (any(isnan.(x)) == false) && (any(isnan.(x1)) == false) && (any(isnan.(x2)) == false)
      @test x ≈ x1
      @test x1 ≈ x2
    end
  end

  # clean up
  safe_rm(seisfile)
  safe_rm(h5file)
end

# A bad data format should produce warnings and "error" channels
printstyled("      bad request logging\n", color=:light_green)
redirect_stdout(out) do
  try
    S = SeisData()
    bad_id = "DE.NENA.99.LUFT"

    # this should return a request error channel with ID = XX.FAIL..001
    get_data!(S, "FDSN", bad_id, v=3, si=false, fmt="sac.zip")

    # this should return a format error channel with ID = XX.FMT..001
    get_data!(S, "FDSN", "UW.LON.."; src="IRIS", s=-600, t=0, v=3, fmt="sac.zip")

    # Check that the info appears where it needs to
    @test S.id[1] == "XX.FAIL..001"
    @test any([occursin("request failed", n) for n in S.notes[1]])
    @test haskey(S.misc[1], "msg")

    @test S.id[S.n] == "XX.FMT..001"
    @test any([occursin("unparseable format", n) for n in S.notes[S.n]])
    @test haskey(S.misc[S.n], "raw")
  catch err
    @warn(string("Bad request logging test failed; caught error ", err))
  end
end

# Potsdam test
printstyled("      request from GFZ\n", color=:light_green)
S = check_get_data("FDSN", "GE.BKB..BH?", src="GFZ", s="2011-03-11T06:00:00", t="2011-03-11T06:05:00", v=0, y=false)
if isempty(S)
  @warn(string("No data from GFZ request; check connection!"))
end

# ❄❄❄❄❄❄❄❄❄❄❄❄❄❄❄❄❄❄ (oh, California...)
printstyled("    servers with special headers:\n", color=:light_green)
rubric = [
  "NCEDC" "BK.MOD..BHE"
  "SCEDC" "CI.SDD..BHZ"
]
❄ = size(rubric, 1)

for i = 1:❄
  printstyled("      ", rubric[i,1], ":\n", color=:light_green)
  ds = now()-Day(1)
  ds -= Millisecond(ds)
  s = string(ds)
  t = string(ds+Hour(1))
  try
    printstyled("        station info\n", color=:light_green)
    S = FDSNsta(rubric[i,2], s=s, t=t, msr=true, src=rubric[i,1])
    printstyled("        trace data\n", color=:light_green)
    get_data!(S, "FDSN", rubric[i,2], src=rubric[i,1], s=s, t=t, msr=true)
    if isempty(S)
      printstyled("        No data; check headers & connection!\n", color=:red)
    end
  catch err
    @warn(string("Request errored; error output below.\n\n", err))
  end
end

printstyled("    issues 42 & 43\n", color=:light_green)
redirect_stdout(out) do
  S = get_data("FDSN", "IU.ANMO.00.LHZ", s="2019-02-14", t="2019-02-15", v=3)

  # This should have two "fail" channels; one from each day with no data
  S = get_data("FDSN", "C0.HAYD.00.LHZ", s="2019-01-01", t="2019-01-05", v=3)
  @test S.n == 3
end

# test IRISPH5 mseed
printstyled("    IRISPH5\n", color=:light_green)
s = DateTime(2016,6,23)
t = s + Minute(1) # this is 250 Hz data
id = "YW.1002..DP?"
printstyled("      station info\n", color=:light_green)
S = FDSNsta(id, src="IRISPH5", s=s, t=t, msr=true)
printstyled("      trace data\n", color=:light_green)
printstyled("        mini-SEED\n", color=:light_green)
get_data!(S, "FDSN", id, src="IRISPH5", s=s, t=t, msr=true)
S = get_data("PH5", id, s=s, t=t, opts="reduction=10")
printstyled("        GeoCSV tspair\n", color=:light_green)
S = get_data("FDSN", id, src="IRISPH5", s=s, t=t, msr=true, fmt="geocsv.tspair")
printstyled("        GeoCSV slist\n", color=:light_green)
S = get_data("FDSN", id, src="IRISPH5", s=s, t=t, msr=true, fmt="geocsv.slist")

# these will fail
# need to implement multi-channel version of read_sac_stream
# segy in FDSN-like requests is NYI
printstyled("        SAC (should fail)\n", color=:light_green)
S = get_data("FDSN", id, src="IRISPH5", s=s, t=t, msr=true, fmt="sac")
printstyled("        segy1 (NYI)\n", color=:light_green)
S = get_data("FDSN", id, src="IRISPH5", s=s, t=t, msr=true, fmt="segy1")
printstyled("        segy2 (NYI)\n", color=:light_green)
S = get_data("FDSN", id, src="IRISPH5", s=s, t=t, msr=true, fmt="segy2")

printstyled("        unsupported format\n", color=:light_green)
S = get_data("FDSN", id, src="IRISPH5", s=s, t=t, fmt="nothing")


text:
import HoTRG_lemon.SpinModule:
	initialize!

type Quantum2dSquareLattice{T} <: Quantum2dLattice{T}
	spinmodel::QuantumSpinModel{T}

	tensorT::Array{T,6}
	tensorTtilde::Array{T,6}
	function Quantum2dSquareLattice{T}(qmodel::QuantumSpinModel{T})
		this = new{T}()
 		this.spinmodel = qmodel
 		initialize!(this)
 		return this
 	end
end

Quantum2dSquareLattice{T}(qmodel::QuantumSpinModel{T}) = Quantum2dSquareLattice{T}(qmodel)

#---
# initialize functions for constructor
function initialize!{T}(lattice::Quantum2dSquareLattice{T})
	lattice.tensorT = initializeTensorT2dQ(lattice.spinmodel)
	lattice.tensorTtilde = initializeTensorT2dQ(lattice.spinmodel; tilde = true)
end

#---
# functions
function getCoarserate(lattice::Quantum2dSquareLattice)
	return 2
end

function getCoarserate(lattice::Quantum2dSquareLattice, whichAxis::AbstractString)
	if whichAxis == "z"
		return 2
	else
		return 2
	end
end

#---
# functions for test
function testQuantum2dSquareLattice()
	spinmodel = testQuantumIsingModel()
	return buildLattice("quantum_2d_square", spinmodel)
end


text:
"""
    One-moment bulk microphysics scheme, which includes:

  - terminal velocity of precipitation
  - condensation and evaporation of cloud liquid water and
    deposition and sublimation of cloud ice (relaxation to equilibrium)
  - autoconversion of cloud liquid water into rain and of cloud ice into snow
  - accretion due to collisions between categories of condensed species
  - evaporation and sublimation of hydrometeors
  - melting of snow into rain
"""
module Microphysics_1M

import SpecialFunctions

import Thermodynamics
import CloudMicrophysics
import CLIMAParameters

const SF = SpecialFunctions
const TD = Thermodynamics
const CO = CloudMicrophysics.Common
const CP = CLIMAParameters
const CP_planet = CLIMAParameters.Planet
const CP_micro = CLIMAParameters.Atmos.Microphysics
const APS = CP.AbstractParameterSet

# Additional type hierarchy to dispatch over for some microphysics parameters
"""
    AbstractCloudType

The top-level super-type for cloud liquid water and cloud ice types
"""
abstract type AbstractCloudType end

"""
    AbstractPrecipType

The top-level super-type for precipitation types (rain and snow)
"""
abstract type AbstractPrecipType end

"""
    LiquidType

The type for cloud liquid water condensate
"""
struct LiquidType <: AbstractCloudType end

"""
    IceType

The type for cloud ice condensate
"""
struct IceType <: AbstractCloudType end

"""
    RainType

The type for rain
"""
struct RainType <: AbstractPrecipType end

"""
    SnowType

The type for snow
"""
struct SnowType <: AbstractPrecipType end

E(param_set::APS, ::LiquidType, ::RainType) = CP_micro.E_liq_rai(param_set::APS)
E(param_set::APS, ::LiquidType, ::SnowType) = CP_micro.E_liq_sno(param_set::APS)
E(param_set::APS, ::IceType, ::RainType) = CP_micro.E_ice_rai(param_set::APS)
E(param_set::APS, ::IceType, ::SnowType) = CP_micro.E_ice_sno(param_set::APS)
E(param_set::APS, ::RainType, ::SnowType) = CP_micro.E_rai_sno(param_set::APS)
E(param_set::APS, ::SnowType, ::RainType) = CP_micro.E_rai_sno(param_set::APS)

export τ_relax

export terminal_velocity

export conv_q_vap_to_q_liq_ice

export conv_q_liq_to_q_rai
export conv_q_ice_to_q_sno

export accretion
export accretion_rain_sink
export accretion_snow_rain

export evaporation_sublimation
export snow_melt

"""
    v0_rai(param_set, ρ)

 - `param_set` - abstract set with Earth parameters
 - `ρ` air density

Returns the proportionality coefficient in terminal velocity(r/r0).
"""
function v0_rai(param_set::APS, ρ::FT) where {FT <: Real}

    _ρ_cloud_liq::FT = CP_planet.ρ_cloud_liq(param_set)
    _C_drag::FT = CP_micro.C_drag(param_set)
    _grav::FT = CP_planet.grav(param_set)
    _r0_rai::FT = CP_micro.r0_rai(param_set)

    return sqrt(
        FT(8 / 3) / _C_drag * (_ρ_cloud_liq / ρ - FT(1)) * _grav * _r0_rai,
    )
end

"""
    n0_sno(param_set, q_sno, ρ)

 - `param_set` - abstract set with Earth parameters
 - `q_sno` -  snow specific humidity
 - `ρ` - air density

Returns the intercept parameter of the assumed Marshall-Palmer distribution of
snow particles.
"""
function n0_sno(param_set::APS, q_sno::FT, ρ::FT) where {FT <: Real}

    _ν_sno::FT = CP_micro.ν_sno(param_set)
    _μ_sno::FT = CP_micro.μ_sno(param_set)

    # TODO               this max should be replaced by
    #                    limiting inside a PhasePartition struct for
    #                    precipitation (once it is implemented)
    return _μ_sno * (ρ * max(0, q_sno))^_ν_sno
end

"""
    unpack_params(param_set, micro, ρ, q_)

 - `param_set` - abstract set with Earth parameters
 - `micro` - type for cloud ice, rain or snow
 - `q_` - specific humidity
 - `ρ` - air density

Utility function that unpacks microphysics parameters.
"""
function unpack_params(
    param_set::APS,
    ice::IceType,
    ρ::FT,
    q_ice::FT,
) where {FT <: Real}
    #TODO - make ρ and q_ice optional
    _n0_ice::FT = CP_micro.n0_ice(param_set)
    _r0_ice::FT = CP_micro.r0_ice(param_set)

    _m0_ice::FT = CP_micro.m0_ice(param_set)
    _me_ice::FT = CP_micro.me_ice(param_set)

    _χm_ice::FT = CP_micro.χm_ice(param_set)
    _Δm_ice::FT = CP_micro.Δm_ice(param_set)

    return (_n0_ice, _r0_ice, _m0_ice, _me_ice, _χm_ice, _Δm_ice)
end
function unpack_params(
    param_set::APS,
    rain::RainType,
    ρ::FT,
    q_rai::FT,
) where {FT <: Real}
    #TODO - make q_rai optional
    _n0_rai::FT = CP_micro.n0_rai(param_set)
    _r0_rai::FT = CP_micro.r0_rai(param_set)

    _m0_rai::FT = CP_micro.m0_rai(param_set)
    _me_rai::FT = CP_micro.me_rai(param_set)
    _a0_rai::FT = CP_micro.a0_rai(param_set)
    _ae_rai::FT = CP_micro.ae_rai(param_set)
    _v0_rai::FT = v0_rai(param_set, ρ)
    _ve_rai::FT = CP_micro.ve_rai(param_set)

    _χm_rai::FT = CP_micro.χm_rai(param_set)
    _Δm_rai::FT = CP_micro.Δm_rai(param_set)
    _χa_rai::FT = CP_micro.χa_rai(param_set)
    _Δa_rai::FT = CP_micro.Δa_rai(param_set)
    _χv_rai::FT = CP_micro.χv_rai(param_set)
    _Δv_rai::FT = CP_micro.Δv_rai(param_set)

    return (
        _n0_rai,
        _r0_rai,
        _m0_rai,
        _me_rai,
        _χm_rai,
        _Δm_rai,
        _a0_rai,
        _ae_rai,
        _χa_rai,
        _Δa_rai,
        _v0_rai,
        _ve_rai,
        _χv_rai,
        _Δv_rai,
    )
end
function unpack_params(
    param_set::APS,
    snow::SnowType,
    ρ::FT,
    q_sno::FT,
) where {FT <: Real}

    _n0_sno::FT = n0_sno(param_set, q_sno, ρ)
    _r0_sno::FT = CP_micro.r0_sno(param_set)

    _m0_sno::FT = CP_micro.m0_sno(param_set)
    _me_sno::FT = CP_micro.me_sno(param_set)
    _a0_sno::FT = CP_micro.a0_sno(param_set)
    _ae_sno::FT = CP_micro.ae_sno(param_set)
    _v0_sno::FT = CP_micro.v0_sno(param_set)
    _ve_sno::FT = CP_micro.ve_sno(param_set)

    _χm_sno::FT = CP_micro.χm_sno(param_set)
    _Δm_sno::FT = CP_micro.Δm_sno(param_set)
    _χa_sno::FT = CP_micro.χa_sno(param_set)
    _Δa_sno::FT = CP_micro.Δa_sno(param_set)
    _χv_sno::FT = CP_micro.χv_sno(param_set)
    _Δv_sno::FT = CP_micro.Δv_sno(param_set)

    return (
        _n0_sno,
        _r0_sno,
        _m0_sno,
        _me_sno,
        _χm_sno,
        _Δm_sno,
        _a0_sno,
        _ae_sno,
        _χa_sno,
        _Δa_sno,
        _v0_sno,
        _ve_sno,
        _χv_sno,
        _Δv_sno,
    )
end

"""
    lambda(q, ρ, n0, m0, me, r0, χm, Δm)

 - `q` - specific humidity of rain, ice or snow
 - `ρ` - air density
 - `n0` - size distribution parameter
 - `m0`, `me`, `χm`, `Δm`, `r0` - mass(radius) parameters

Returns the rate parameter of the assumed size distribution of
particles (rain drops, ice crystals, snow crystals).
"""
function lambda(
    q::FT,
    ρ::FT,
    n0::FT,
    m0::FT,
    me::FT,
    r0::FT,
    χm::FT,
    Δm::FT,
) where {FT <: Real}

    λ::FT = FT(0)

    if q > FT(0)
        λ =
            (
                χm * m0 * n0 * SF.gamma(me + Δm + FT(1)) / ρ / q / r0^(me + Δm)
            )^FT(1 / (me + Δm + 1))
    end
    return λ
end

"""
    τ_relax(param_set, liquid)
    τ_relax(param_set, ice)

 - `param_set` - abstract set with Earth parameters
 - `liquid` - a type for cloud liquid water
 - `ice` - a type for cloud ice

Returns the relaxation timescale for condensation and evaporation of
cloud liquid water or the relaxation timescale for sublimation and
deposition of cloud ice.
"""
function τ_relax(param_set::APS, liquid::LiquidType)

    _τ_relax = CP_micro.τ_cond_evap(param_set)
    return _τ_relax
end
function τ_relax(param_set::APS, ice::IceType)

    _τ_relax = CP_micro.τ_sub_dep(param_set)
    return _τ_relax
end

"""
    terminal_velocity(param_set, precip, ρ, q_)

 - `param_set` - abstract set with Earth parameters
 - `precip` - a type for rain or snow
 - `ρ` - air density
 - `q_` - rain or snow specific humidity

Returns the mass weighted average terminal velocity assuming
a Marshall-Palmer (1948) distribution of rain drops and snow crystals.
"""
function terminal_velocity(
    param_set::APS,
    precip::AbstractPrecipType,
    ρ::FT,
    q_::FT,
) where {FT <: Real}
    fall_w = FT(0)
    if q_ > FT(0)

        (_n0, _r0, _m0, _me, _χm, _Δm, _a0, _ae, _χa, _Δa, _v0, _ve, _χv, _Δv) =
            unpack_params(param_set, precip, ρ, q_)

        _λ::FT = lambda(q_, ρ, _n0, _m0, _me, _r0, _χm, _Δm)

        fall_w =
            _χv *
            _v0 *
            (_λ * _r0)^(-_ve - _Δv) *
            SF.gamma(_me + _ve + _Δm + _Δv + FT(1)) /
            SF.gamma(_me + _Δm + FT(1))
    end

    return fall_w
end

"""
    conv_q_vap_to_q_liq_ice(param_set, liquid, q_sat, q)
    conv_q_vap_to_q_liq_ice(param_set, ice, q_sat, q)

 - `param_set` - abstract set with Earth parameters
 - `liquid` - a type for cloud water
 - `ice` - a type for cloud ice
 - `q_sat` - PhasePartition at equilibrium
 - `q` - current PhasePartition

Returns the cloud water tendency due to condensation and evaporation
or cloud ice tendency due to sublimation and vapor deposition.
The tendency is obtained assuming a relaxation to equilibrium with
a constant timescale.
"""
function conv_q_vap_to_q_liq_ice(
    param_set::APS,
    liquid::LiquidType,
    q_sat::TD.PhasePartition{FT},
    q::TD.PhasePartition{FT},
) where {FT <: Real}

    _τ_cond_evap::FT = τ_relax(param_set, liquid)

    return (q_sat.liq - q.liq) / _τ_cond_evap
end
function conv_q_vap_to_q_liq_ice(
    param_set::APS,
    ice::IceType,
    q_sat::TD.PhasePartition{FT},
    q::TD.PhasePartition{FT},
) where {FT <: Real}

    _τ_sub_dep::FT = τ_relax(param_set, ice)

    return (q_sat.ice - q.ice) / _τ_sub_dep
end

"""
    conv_q_liq_to_q_rai(param_set, q_liq)

 - `param_set` - abstract set with Earth parameters
 - `q_liq` - liquid water specific humidity

Returns the q_rai tendency due to collisions between cloud droplets
(autoconversion), parametrized following Kessler (1995).
"""
function conv_q_liq_to_q_rai(param_set::APS, q_liq::FT) where {FT <: Real}

    _τ_acnv_rai::FT = CP_micro.τ_acnv_rai(param_set)
    _q_liq_threshold::FT = CP_micro.q_liq_threshold(param_set)

    return max(0, q_liq - _q_liq_threshold) / _τ_acnv_rai
end

"""
    conv_q_ice_to_q_sno_no_supersat(param_set, q_ice)

 - `param_set` - abstract set with Earth parameters
 - `q_ice` -  cloud ice specific humidity

Returns the q_sno tendency due to autoconversion from ice.
This is a simplified version of a snow autoconversion rate that can be used in
simulations where there is no supersaturation
(for example in TC.jl when using saturation adjustment).
"""
function conv_q_ice_to_q_sno_no_supersat(
    param_set::APS,
    q_ice::FT,
) where {FT <: Real}

    _τ_acnv_sno::FT = CP_micro.τ_acnv_sno(param_set)
    _q_ice_threshold::FT = CP_micro.q_ice_threshold(param_set)

    return max(0, q_ice - _q_ice_threshold) / _τ_acnv_sno
end

"""
    conv_q_ice_to_q_sno(param_set, q, ρ, T)

 - `param_set` - abstract set with Earth parameters
 - `q` - phase partition
 - `ρ` - air density
 - `T` - air temperature

Returns the q_sno tendency due to autoconversion from ice.
Parameterized following Harrington et al. (1996) and Kaul et al. (2015).
"""
function conv_q_ice_to_q_sno(
    param_set::APS,
    q::TD.PhasePartition{FT},
    ρ::FT,
    T::FT,
) where {FT <: Real}
    acnv_rate = FT(0)
    _S::FT = TD.supersaturation(param_set, q, ρ, T, TD.Ice())

    if (q.ice > FT(0) && _S > FT(0))

        _G::FT = CO.G_func(param_set, T, TD.Ice())

        _r_ice_snow::FT = CP_micro.r_ice_snow(param_set)

        (_n0, _r0, _m0, _me, _χm, _Δm) =
            unpack_params(param_set, IceType(), ρ, q.ice)

        _λ::FT = lambda(q.ice, ρ, _n0, _m0, _me, _r0, _χm, _Δm)

        acnv_rate =
            4 * FT(π) * _S * _G * _n0 / ρ *
            exp(-_λ * _r_ice_snow) *
            (
                _r_ice_snow^FT(2) / (_me + _Δm) +
                (_r_ice_snow * _λ + FT(1)) / _λ^FT(2)
            )
    end
    return acnv_rate
end

"""
    accretion(param_set, cloud, precip, q_clo, q_pre, ρ)

 - `param_set` - abstract set with Earth parameters
 - `cloud` - type for cloud water or cloud ice
 - `precip` - type for rain or snow
 - `q_clo` - cloud water or cloud ice specific humidity
 - `q_pre` - rain water or snow specific humidity
 - `ρ` - rain water or snow specific humidity

Returns the source of precipitating water (rain or snow)
due to collisions with cloud water (liquid or ice).
"""
function accretion(
    param_set::APS,
    cloud::AbstractCloudType,
    precip::AbstractPrecipType,
    q_clo::FT,
    q_pre::FT,
    ρ::FT,
) where {FT <: Real}

    accr_rate = FT(0)
    if (q_clo > FT(0) && q_pre > FT(0))

        (_n0, _r0, _m0, _me, _χm, _Δm, _a0, _ae, _χa, _Δa, _v0, _ve, _χv, _Δv) =
            unpack_params(param_set, precip, ρ, q_pre)

        _λ::FT = lambda(q_pre, ρ, _n0, _m0, _me, _r0, _χm, _Δm)
        _E::FT = E(param_set, cloud, precip)

        accr_rate =
            q_clo * _E * _n0 * _a0 * _v0 * _χa * _χv / _λ *
            SF.gamma(_ae + _ve + _Δa + _Δv + FT(1)) /
            (_λ * _r0)^(_ae + _ve + _Δa + _Δv)
    end
    return accr_rate
end

"""
    accretion_rain_sink(param_set, q_ice, q_rai, ρ)

 - `param_set` - abstract set with Earth parameters
 - `q_ice` - cloud ice specific humidity
 - `q_rai` - rain water specific humidity
 - `ρ` - air density

Returns the sink of rain water (partial source of snow) due to collisions
with cloud ice.
"""
function accretion_rain_sink(
    param_set::APS,
    q_ice::FT,
    q_rai::FT,
    ρ::FT,
) where {FT <: Real}

    accr_rate = FT(0)
    if (q_ice > FT(0) && q_rai > FT(0))

        (_n0_ice, _r0_ice, _m0_ice, _me_ice, _χm_ice, _Δm_ice) =
            unpack_params(param_set, IceType(), ρ, q_ice)

        (
            _n0_rai,
            _r0_rai,
            _m0_rai,
            _me_rai,
            _χm_rai,
            _Δm_rai,
            _a0_rai,
            _ae_rai,
            _χa_rai,
            _Δa_rai,
            _v0_rai,
            _ve_rai,
            _χv_rai,
            _Δv_rai,
        ) = unpack_params(param_set, RainType(), ρ, q_rai)

        _E::FT = E(param_set, IceType(), RainType())

        _λ_rai::FT = lambda(
            q_rai,
            ρ,
            _n0_rai,
            _m0_rai,
            _me_rai,
            _r0_rai,
            _χm_rai,
            _Δm_rai,
        )
        _λ_ice::FT = lambda(
            q_ice,
            ρ,
            _n0_ice,
            _m0_ice,
            _me_ice,
            _r0_ice,
            _χm_ice,
            _Δm_ice,
        )

        accr_rate =
            _E / ρ *
            _n0_rai *
            _n0_ice *
            _m0_rai *
            _a0_rai *
            _v0_rai *
            _χm_rai *
            _χa_rai *
            _χv_rai / _λ_ice / _λ_rai * SF.gamma(
                _me_rai +
                _ae_rai +
                _ve_rai +
                _Δm_rai +
                _Δa_rai +
                _Δv_rai +
                FT(1),
            ) /
            (
                _r0_rai * _λ_rai
            )^(_me_rai + _ae_rai + _ve_rai + _Δm_rai + _Δa_rai + _Δv_rai)
    end
    return accr_rate
end

"""
    accretion_snow_rain(param_set, type_i, type_j, q_i, q_j, ρ)

 - `i` - snow for temperatures below freezing
         or rain for temperatures above freezing
 - `j` - rain for temperatures below freezing
         or rain for temperatures above freezing
 - `param_set` - abstract set with Earth parameters
 - `type_i`, `type_j` - a type for snow or rain
 - `q_` - specific humidity of snow or rain
 - `ρ` - air density

Returns the accretion rate between rain and snow.
Collisions between rain and snow result in
snow at temperatures below freezing and in rain at temperatures above freezing.
"""
function accretion_snow_rain(
    param_set::APS,
    type_i::AbstractPrecipType,
    type_j::AbstractPrecipType,
    q_i::FT,
    q_j::FT,
    ρ::FT,
) where {FT <: Real}

    accr_rate = FT(0)
    if (q_i > FT(0) && q_j > FT(0))

        (
            _n0_i,
            _r0_i,
            _m0_i,
            _me_i,
            _χm_i,
            _Δm_i,
            _a0_i,
            _ae_i,
            _χa_i,
            _Δa_i,
            _v0_i,
            _ve_i,
            _χv_i,
            _Δv_i,
        ) = unpack_params(param_set, type_i, ρ, q_i)
        (
            _n0_j,
            _r0_j,
            _m0_j,
            _me_j,
            _χm_j,
            _Δm_j,
            _a0_j,
            _ae_j,
            _χa_j,
            _Δa_j,
            _v0_j,
            _ve_j,
            _χv_j,
            _Δv_j,
        ) = unpack_params(param_set, type_j, ρ, q_j)

        _E_ij::FT = E(param_set, type_i, type_j)

        _λ_i::FT = lambda(q_i, ρ, _n0_i, _m0_i, _me_i, _r0_i, _χm_i, _Δm_i)
        _λ_j::FT = lambda(q_j, ρ, _n0_j, _m0_j, _me_j, _r0_j, _χm_j, _Δm_j)

        _v_ti = terminal_velocity(param_set, type_i, ρ, q_i)
        _v_tj = terminal_velocity(param_set, type_j, ρ, q_j)

        accr_rate =
            FT(π) / ρ *
            _n0_i *
            _n0_j *
            _m0_j *
            _χm_j *
            _E_ij *
            abs(_v_ti - _v_tj) / _r0_j^(_me_j + _Δm_j) * (
                FT(2) * SF.gamma(_me_j + _Δm_j + FT(1)) / _λ_i^FT(3) /
                _λ_j^(_me_j + _Δm_j + FT(1)) +
                FT(2) * SF.gamma(_me_j + _Δm_j + FT(2)) / _λ_i^FT(2) /
                _λ_j^(_me_j + _Δm_j + FT(2)) +
                SF.gamma(_me_j + _Δm_j + FT(3)) / _λ_i /
                _λ_j^(_me_j + _Δm_j + FT(3))
            )
    end
    return accr_rate
end

"""
    evaporation_sublimation(param_set, rain, q, q_rai, ρ, T)
    evaporation_sublimation(param_set, snow, q, q_sno, ρ, T)

 - `param_set` - abstract set with Earth parameters
 - `rain` - a type for rain
 - `snow` - a type for snow
 - `q` - phase partition
 - `q_rai` - rain specific humidity
 - `q_sno` - snow specific humidity
 - `ρ` - air density
 - `T` - air temperature

Returns the tendency due to rain evaporation or snow sublimation.
"""
function evaporation_sublimation(
    param_set::APS,
    rain::RainType,
    q::TD.PhasePartition{FT},
    q_rai::FT,
    ρ::FT,
    T::FT,
) where {FT <: Real}
    evap_subl_rate = FT(0)
    _S::FT = TD.supersaturation(param_set, q, ρ, T, TD.Liquid())

    if (q_rai > FT(0) && _S < FT(0))

        _a_vent::FT = CP_micro.a_vent_rai(param_set)
        _b_vent::FT = CP_micro.b_vent_rai(param_set)
        _ν_air::FT = CP_micro.ν_air(param_set)
        _D_vapor::FT = CP_micro.D_vapor(param_set)

        _G::FT = CO.G_func(param_set, T, TD.Liquid())

        (_n0, _r0, _m0, _me, _χm, _Δm, _a0, _ae, _χa, _Δa, _v0, _ve, _χv, _Δv) =
            unpack_params(param_set, rain, ρ, q_rai)

        _λ::FT = lambda(q_rai, ρ, _n0, _m0, _me, _r0, _χm, _Δm)

        evap_subl_rate =
            4 * FT(π) * _n0 / ρ * _S * _G / _λ^FT(2) * (
                _a_vent +
                _b_vent * (_ν_air / _D_vapor)^FT(1 / 3) /
                (_r0 * _λ)^((_ve + _Δv) / FT(2)) *
                (FT(2) * _v0 * _χv / _ν_air / _λ)^FT(1 / 2) *
                SF.gamma((_ve + _Δv + FT(5)) / FT(2))
            )
    end
    # only evaporation is considered for rain
    return min(0, evap_subl_rate)
end
function evaporation_sublimation(
    param_set::APS,
    snow::SnowType,
    q::TD.PhasePartition{FT},
    q_sno::FT,
    ρ::FT,
    T::FT,
) where {FT <: Real}
    evap_subl_rate = FT(0)
    if q_sno > FT(0)

        _a_vent::FT = CP_micro.a_vent_sno(param_set)
        _b_vent::FT = CP_micro.b_vent_sno(param_set)
        _ν_air::FT = CP_micro.ν_air(param_set)
        _D_vapor::FT = CP_micro.D_vapor(param_set)

        _S::FT = TD.supersaturation(param_set, q, ρ, T, TD.Ice())
        _G::FT = CO.G_func(param_set, T, TD.Ice())

        (_n0, _r0, _m0, _me, _χm, _Δm, _a0, _ae, _χa, _Δa, _v0, _ve, _χv, _Δv) =
            unpack_params(param_set, snow, ρ, q_sno)
        _λ::FT = lambda(q_sno, ρ, _n0, _m0, _me, _r0, _χm, _Δm)

        evap_subl_rate =
            4 * FT(π) * _n0 / ρ * _S * _G / _λ^FT(2) * (
                _a_vent +
                _b_vent * (_ν_air / _D_vapor)^FT(1 / 3) /
                (_r0 * _λ)^((_ve + _Δv) / FT(2)) *
                (FT(2) * _v0 * _χv / _ν_air / _λ)^FT(1 / 2) *
                SF.gamma((_ve + _Δv + FT(5)) / FT(2))
            )
    end
    return evap_subl_rate
end

"""
    snow_melt(param_set, q_sno, ρ, T)

 - `param_set` - abstract set with Earth parameters
 - `q_sno` - snow water specific humidity
 - `ρ` - air density
 - `T` - air temperature

Returns the tendency due to snow melt.
"""
function snow_melt(param_set::APS, q_sno::FT, ρ::FT, T::FT) where {FT <: Real}

    snow_melt_rate = FT(0)
    _T_freeze = CP_planet.T_freeze(param_set)

    if (q_sno > FT(0) && T > _T_freeze)

        _a_vent::FT = CP_micro.a_vent_sno(param_set)
        _b_vent::FT = CP_micro.b_vent_sno(param_set)
        _ν_air::FT = CP_micro.ν_air(param_set)
        _D_vapor::FT = CP_micro.D_vapor(param_set)
        _K_therm::FT = CP_micro.K_therm(param_set)

        L = TD.latent_heat_fusion(param_set, T)

        (_n0, _r0, _m0, _me, _χm, _Δm, _a0, _ae, _χa, _Δa, _v0, _ve, _χv, _Δv) =
            unpack_params(param_set, SnowType(), ρ, q_sno)
        _λ::FT = lambda(q_sno, ρ, _n0, _m0, _me, _r0, _χm, _Δm)

        snow_melt_rate =
            4 * FT(π) * _n0 / ρ * _K_therm / L * (T - _T_freeze) / _λ^FT(2) * (
                _a_vent +
                _b_vent * (_ν_air / _D_vapor)^FT(1 / 3) /
                (_r0 * _λ)^((_ve + _Δv) / FT(2)) *
                (FT(2) * _v0 * _χv / _ν_air / _λ)^FT(1 / 2) *
                SF.gamma((_ve + _Δv + FT(5)) / FT(2))
            )
    end
    return snow_melt_rate
end

end #module Microphysics_1M.jl


text:
using SemanticModels.Parsers
import Base: ==

mutable struct Edge{S,T,R}
    func::S
    args::T
    ret::R
end

==(e::Edge,E::Edge) = begin (@show e.func==E.func && @show e.args == E.args && @show e.ret == E.ret) end

mutable struct CallArg{S,T,U,V}
    f::S
    arg::T
    T::U
    val::V
end
write(io, ca::CallArg) = write(io, "$(repr(ca.f))($(repr(f.arg))::$(repr(f.T))=$(repr(f.val)))")

mutable struct RetArg{S,T,U,V}
    f::S
    arg::T
    T::U
    val::V
end

function Edges(snapshots)
    stack = Any[]
    fs = Any[]
    for snap in snapshots
        #@show (x->x.f).(stack)
        if !isa(snap, RetArg)
            push!(stack, snap)
        end
        if isa(snap, RetArg)
            if length(stack) <= 0
                @show snap
                @show fs
                error("stack shouldn't have been empty here")
            end
            callargs = pop!(stack)
            framename = first(callargs).f
            t = Edge(framename, (x->x.T).(callargs), snap.T)
            push!(fs, t)
        end
    end
    return fs
end

annotate(x::Any) = x

function annotate(x::Expr)
    # Grabing the values and types of the inputs
    if head(x) == :function
        b = bodyblock(x)
        argl= argslist(x)
        funcname = string(argl[1])
        if length(argl) > 1
            argl= argslist(x)[2:end]
        else
            argl = []
        end
        sargl = string.(argl)
        v = Expr(:vect)
        for (i, t) in enumerate(zip(sargl,argl))
            name, val = t
            T = :(typeof($val))
            q = :(CallArg($funcname, $name, $T, $val))
            push!(v.args, q)
        end
        insert!(b, 1, :(store($v)))

    end
    # Grabbing the return types
    if head(x) == :return
        g = gensym("ann")
        s = string(x)
        r = x.args[1]

        return quote
            $g = $r
            T = typeof($g)
            store(RetArg("ret", $s, T, $g))
            return $g
        end
     end
    return x
end

function wrap(expr::Expr)
    b = expr.args
    g = gensym("wrap_storage")
    insert!(b, 1, :(store(args) = begin push!($g, args) end ))
    insert!(b, 1, :($g = Any[]))
    push!(b, :(edgelist=$g))
    # push!(b, :(Edges($g)))
    expr
end

"""    typegraph(expr::Expr)

annotate a code expression so that when you eval it, you get the typegraph.
used in the macro @typegraph.

Note: Does not yet support docstrings, kwargs, or varargs.
"""
function typegraph(expr::Expr)
    return wrap(postwalk(annotate, expr))
end

"""    @typegraph(expr::Expr)

extract a typegraph from an expression by annotation and execution.

Note: Does not yet support docstrings, kwargs, or varargs.
"""
macro typegraph(expr)
    return typegraph(expr)
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

# TODO: deprecated
export
    ComplementGraphView, ComplementDiGraphView,
    complementgraph

struct ComplementGraphView{T<:UndirectedGraph} <: GraphView
    graph::T
end

struct ComplementDiGraphView{T<:DirectedGraph} <: DiGraphView
    graph::T
end

# TODO: use traits
ComplementGraph = Union{ComplementGraphView,ComplementDiGraphView}

complementgraph(graph::UndirectedGraph) = ComplementGraphView(graph)
complementgraph(graph::DirectedGraph) = ComplementDiGraphView(graph)

function getedge(view::ComplementGraphView, idx)
    throw(ErrorException("Edge indexing not supported for this view."))
end

function getedge(view::ComplementGraphView, u, v)
    (v in neighbors(view.graph, u)) || throw(KeyError(v))
    return edgetype(view.graph)(u, v)
end

hasedge(view::ComplementGraphView, u, v) = !hasedge(view.graph, u, v)

function edgesiter(view::ComplementGraphView)
    throw(ErrorException("Edge indexing not supported for this view."))
end

function neighbors(view::ComplementGraphView, idx)
    throw(ErrorException("Edge indexing not supported for this view."))
end

function adjacencies(view::ComplementGraphView, idx)
    ns = nodeset(view.graph)
    pop!(ns, idx)
    return setdiff(ns, adjacencies(view.graph, idx))
end

function incidences(view::ComplementGraphView, idx)
    throw(ErrorException("Edge indexing not supported for this view."))
end

function incidentedges(view::ComplementGraphView, idx)
    throw(ErrorException("Edge indexing not supported for this view."))
end

function edgecount(view::ComplementGraphView)
    n = nodecount(view.graph)
    return n * (n - 1) / 2 - edgecount(view.graph)
end


text:
using Pigeon
using PyCall

function load_path_msg(fname, msg=convert(PyCall.PyObject, Pigeon.osprey.msg.path()))
    open(fname) do f
        Pigeon.osprey.msg.path[:deserialize](msg, read(f))
        msg[:header][:frame_id] = String(msg[:header][:frame_id])
        convert(Pigeon.osprey.msg.path, msg)
    end
end

function load_all_path_msgs(dir=joinpath(@__DIR__, "path"))
    Dict(fname => load_path_msg(joinpath(dir, fname)) for fname in readdir(dir) if endswith(fname, ".msg"))
end

@eval(Pigeon, const TEST_PATHS = $(load_all_path_msgs()))

max_stars_repo_name: Michelle-NYX/CollisionAvoidance

max_issues_repo_name: Michelle-NYX/CollisionAvoidance

max_forks_repo_name: Michelle-NYX/CollisionAvoidance


text:
""" Compute the Mahalanobis Distance to center (MDC)
    in a dataset 'data' for a given perturbation of indices 'indpert'
    compared to a reference of indices 'indref'."""
function distance_mahalanobis_center(data, indpert, indref)
    setPert = Matrix(data[indpert,:])
    setRef = Matrix(data[indref,:])

    mdCenter = dropdims(mean(setRef, dims = 1), dims = 1)
    mdCov = cov(setRef)

    pertCenter = dropdims(mean(setPert, dims = 1), dims = 1)

    MD = mahalanobis(pertCenter, mdCenter, mdCov)

    return(MD)
end

""" Permute labels and compute the Mahalanobis Distance to center (MDC)
    in a dataset 'data' for a given perturbation of indices 'indpert'
    compared to a reference of indices 'indref', to create an empirical distribution."""
function shuffled_distance_mahalanobis_center(data, indpert, indref; nb_rep = 250)
    setPert = data[indpert,:]
    setRef = data[indref,:]
    set = Matrix(vcat(setRef, setPert))

    function iterShufMD()
        nset = size(set, 1)
        npert = size(setPert, 1)
        nref = size(setRef, 1)
        shuffSet = set[sample(1:nset, nset; replace = false),:]
        # Take random subsets of corresponding sizes
        shuffSetPert = shuffSet[1:npert,:]
        shuffSetRef = shuffSet[(npert+1):(npert+nref),:]

        # Compute Mahalanobis Distance

        mdCenter = dropdims(mean(shuffSetRef, dims = 1), dims = 1)
        mdCov = cov(shuffSetRef)

        pertCenter = dropdims(mean(shuffSetPert, dims = 1), dims = 1)

        MD = mahalanobis(pertCenter, mdCenter, mdCov)
        return(MD)
    end

    return(map(x -> iterShufMD(), 1:nb_rep))
end

""" Compute the median Mahalanobis Distance (MD)
    in a dataset 'data' for a given perturbation of indices 'indpert'
    compared to a reference of indices 'indref'."""
function distance_mahalanobis_median(data, indpert, indref)
    setPert = data[indpert,:]
    setRef = Matrix(data[indref,:])

    mdCenter = dropdims(mean(setRef, dims = 1), dims = 1)
    mdCov = cov(setRef)

    MD = median(map(x -> mahalanobis(x, mdCenter, mdCov), eachrow(setPert)))
    return(MD)
end

""" Permute labels and compute the median Mahalanobis Distance (RMD)
    in a dataset 'data' for a given perturbation of indices 'indpert'
    compared to a reference of indices 'indref', to create an empirical distribution."""
function shuffled_distance_mahalanobis_median(data, indpert, indref; nb_rep = 250)
    setPert = data[indpert,:]
    setRef = data[indref,:]
    set = Matrix(vcat(setRef, setPert))

    function iterShufMD()
        nset = size(set, 1)
        npert = size(setPert, 1)
        nref = size(setRef, 1)
        shuffSet = set[sample(1:nset, nset; replace = false),:]
        # Take random subsets of corresponding sizes
        shuffSetPert = shuffSet[1:npert,:]
        shuffSetRef = shuffSet[(npert+1):(npert+nref),:]

        # Compute Mahalanobis Distance

        mdCenter = dropdims(mean(shuffSetRef, dims = 1), dims = 1)
        mdCov = cov(shuffSetRef)

        MD = median(map(x -> mahalanobis(x, mdCenter, mdCov), eachrow(DataFrame(shuffSetPert))))
        return(MD)
    end

    return(map(x -> iterShufMD(), 1:nb_rep))
end

""" Compute the median Robust Mahalanobis Distance (RMD)
    in a dataset 'data' for a given perturbation of indices 'indpert'
    compared to a reference of indices 'indref'.
    See https://e-archivo.uc3m.es/bitstream/handle/10016/24613/ws201710.pdf
    This function calls R using RCall and a seed is set by default to ensure
    the results are reproducible. If you don't want that, for instance if you
    also use RCall and rely on another seed, set 'r_seed' to false."""
function distance_robust_mahalanobis_median(data, indpert, indref; r_seed = true)
    setPert = data[indpert,:]
    setRef = data[indref,:]

    # Ensure that we have enough points to compute distance
    if ((size(setPert)[1] < 2*size(data, 2))|(size(setRef)[1] < 2*size(data, 2)))
        return(missing)
    end
    # NB: having less points than twice the number of dimensions leads to singularity

    # Compute Minimum Covariance Determinant and corresponding Robust Mahalanobis Distance
    @rput setRef
    @rput r_seed

    R"""
    if (!require("robustbase")) install.packages("robustbase",
                                repos = "https://cloud.r-project.org")
    library(robustbase)

    if (r_seed){set.seed(777)}
    mcd <- covMcd(setRef)
    mcdCenter <- mcd$center
    mcdCov <- mcd$cov
    """
    @rget mcdCenter
    @rget mcdCov

    RMD = median(map(x -> mahalanobis(x, mcdCenter, mcdCov), eachrow(setPert)))
    return(RMD)
end

""" Permute labels and compute the median Robust Mahalanobis Distance (RMD)
    in a dataset 'data' for a given perturbation of indices 'indpert'
    compared to a reference of indices 'indref', to create an empirical distribution.
    This function calls R using RCall and a seed is set by default to ensure
    the results are reproducible. If you don't want that, for instance if you
    also use RCall and rely on another seed, set 'r_seed' to false."""
function shuffled_distance_robust_mahalanobis_median(data, indpert, indref; nb_rep = 250, r_seed = true)
    setPert = data[indpert,:]
    setRef = data[indref,:]
    set = vcat(setRef, setPert)

    nset = size(set, 1)
    npert = size(setPert, 1)
    nref = size(setRef, 1)

    # Ensure that we have enough points to compute distance
    if ((size(setPert)[1] < 2*size(data, 2))|(size(setRef)[1] < 2*size(data, 2)))
        return(repeat([missing], nb_rep))
    end
    # NB: having less points than twice the number of dimensions leads to singularity

    function iterShufRMD()
        shuffSet = set[sample(1:nset, nset; replace = false),:]
        # Take random subsets of corresponding sizes
        shuffSetPert = shuffSet[1:npert,:]
        shuffSetRef = shuffSet[(npert+1):(npert+nref),:]

        # Compute Minimum Covariance Determinant and corresponding Robust Mahalanobis Distance
        @rput shuffSetRef
        @rput r_seed

        R"""
        if (!require("robustbase")) install.packages("robustbase",
                                    repos = "https://cloud.r-project.org")
        library(robustbase)

        if (r_seed){set.seed(3895)}
        mcd <- covMcd(shuffSetRef)
        mcdCenter <- mcd$center
        mcdCov <- mcd$cov
        """
        @rget mcdCenter
        @rget mcdCov

        RMD = median(map(x -> mahalanobis(x, mcdCenter, mcdCov), eachrow(shuffSetPert)))
        return(RMD)
    end

    return(map(x -> iterShufRMD(), 1:nb_rep))
end

""" Compute the Robust Hellinger Distance (RHD)
    in a dataset `data` for a given perturbation of indices `indpert`
    compared to a reference of indices `indref`.
    This function calls R using RCall and a seed is set by default to ensure
    the results are reproducible. If you don't want that, for instance if you
    also use RCall and rely on another seed, set 'r_seed' to false."""
function distance_robust_hellinger(data, indpert, indref; r_seed = true)
    setPert = data[indpert,:]
    setRef = data[indref,:]

    # Ensure that we have enough points to compute distance
    if ((size(setPert)[1] < 2*size(data, 2))|(size(setRef)[1] < 2*size(data, 2)))
        return(missing)
    end
    # NB: having less points than twice the number of dimensions leads to singularity

    # Compute Minimum Covariance Determinant and corresponding Robust Hellinger Distance
    @rput setRef
    @rput setPert
    @rput r_seed

    R"""
    if (!require("robustbase")) install.packages("robustbase",
                                repos = "https://cloud.r-project.org")
    library(robustbase)

    if (r_seed){set.seed(777)}
    mcd1 <- covMcd(setRef)
    mcdCenter1 <- mcd1$center
    mcdCov1 <- mcd1$cov

    # We set the seed twice to always
    # find the same estimators given
    # the same sample
    if (r_seed){set.seed(777)}
    mcd2 <- covMcd(setPert)
    mcdCenter2 <- mcd2$center
    mcdCov2 <- mcd2$cov
    """
    @rget mcdCenter1
    @rget mcdCov1
    @rget mcdCenter2
    @rget mcdCov2

    RHD = hellinger(mcdCenter1, mcdCov1, mcdCenter2, mcdCov2)
    return(RHD)
end

""" Permute labels and compute the Robust Hellinger Distance (RHD)
    in a dataset `data` for a given perturbation of indices `indpert`
    compared to a reference of indices `indref`, to create an empirical distribution.
    This function calls R using RCall and a seed is set by default to ensure
    the results are reproducible. If you don't want that, for instance if you
    also use RCall and rely on another seed, set 'r_seed' to false."""
function shuffled_distance_robust_hellinger(data, indpert, indref; nb_rep = 250, r_seed = true)
    setPert = data[indpert,:]
    setRef = data[indref,:]
    set = vcat(setRef, setPert)

    nset = size(set, 1)
    npert = size(setPert, 1)
    nref = size(setRef, 1)

    # Ensure that we have enough points to compute distance
    if ((size(setPert)[1] < 2*size(data, 2))|(size(setRef)[1] < 2*size(data, 2)))
        return(repeat([missing], nb_rep))
    end
    # NB: having less points than twice the number of dimensions leads to singularity

    function iterShufRHD()
        shuffSet = set[sample(1:nset, nset; replace = false),:]
        # Take random subsets of corresponding sizes
        shuffSetPert = shuffSet[1:npert,:]
        shuffSetRef = shuffSet[(npert+1):(npert+nref),:]

        # Compute Minimum Covariance Determinant and corresponding Robust Mahalanobis Distance
        @rput shuffSetRef
        @rput shuffSetPert
        @rput r_seed

        R"""
        if (!require("robustbase")) install.packages("robustbase",
                                    repos = "https://cloud.r-project.org")
        library(robustbase)

        if (r_seed){set.seed(777)}
        mcd <- covMcd(shuffSetRef)
        mcdCenter1 <- mcd$center
        mcdCov1 <- mcd$cov

        # We set the seed twice to always
        # find the same estimators given
        # the same sample
        if (r_seed){set.seed(777)}
        mcd <- covMcd(shuffSetPert)
        mcdCenter2 <- mcd$center
        mcdCov2 <- mcd$cov
        """
        @rget mcdCenter1
        @rget mcdCov1
        @rget mcdCenter2
        @rget mcdCov2

        RHD = hellinger(mcdCenter1, mcdCov1, mcdCenter2, mcdCov2)
        return(RHD)
    end

    return(map(x -> iterShufRHD(), 1:nb_rep))
end

""" Compute the Robust Morphological Perturbation Value (RMPV)
    for a given Experiment `e`, for all levels of a column `s`,
    compared to rows matching a given filter `f` or where `s`
    is equal to `ref`.
    The RMPV quantifies the significance of changes between all
    conditions (levels in `s`) and a reference condition (defined
    by the filter `f`).
    In brief, the distance of type `dist` between points of each
    perturbation and points of the reference is computed and its
    statistical significance is defined using a permutation test
    in which the perturbation and reference labels are shuffled
    `nb_rep` times.
    If `process_pool` is a pool of worker processes, they will
    be used for parallel computation in the permutation test.
    This function might call R using RCall and a seed is set by default
    to ensure the results are reproducible. If you don't want that, for
    instance if you also use RCall and rely on another seed, set 'r_seed'
    to false.
    This returns a DataFrame with three columns:
    * `Condition`: the levels in `s`
    * `Distance`: the distance between a condition and the
    reference
    * `RMPV`: the RMPV (empirical p-value corrected for multiple
    testing)
    """
function robust_morphological_perturbation_value end

# If filter provided for reference
function robust_morphological_perturbation_value(e::AbstractExperiment,
                                                 s::Symbol,
                                                 f::AbstractFilter;
                                                 nb_rep::Int64 = 250,
                                                 dist::Symbol = :RobHellinger,
                                                 process_pool = nothing,
                                                 r_seed = true)
    if dist == :RobHellinger
        selected_distance = (x...; kw...)->distance_robust_hellinger(x...;
                                                                     r_seed = r_seed, kw...)
        shuffled_distance = (x...; kw...)->shuffled_distance_robust_hellinger(x...;
                                                                     r_seed = r_seed, kw...)
    elseif dist == :RobMedMahalanobis
        selected_distance = (x...; kw...)->distance_robust_mahalanobis_median(x...;
                                                                     r_seed = r_seed, kw...)
        shuffled_distance = (x...; kw...)->shuffled_distance_robust_mahalanobis_median(x...;
                                                                     r_seed = r_seed, kw...)
    elseif dist == :MedMahalanobis
        selected_distance = distance_mahalanobis_median
        shuffled_distance = shuffled_distance_mahalanobis_median
    elseif dist == :CenterMahalanobis
        selected_distance = distance_mahalanobis_center
        shuffled_distance = shuffled_distance_mahalanobis_center
    else
        throw(DomainError(dist, "Invalid `dist` argument. "*
                                "Only :RobHellinger "*
                                "and :RobMedMahalanobis "*
                                "and :MedMahalanobis "*
                                "and :CenterMahalanobis "*
                                "are supported"))
    end

    # All conditions considered
    cnd_levels = levels(e.data[e.selected_entries,s])

    # Actual observed distances
    allRD = map(x -> selected_distance(e.data[:,e.selected_features],
                                        filter_entries(e, Filter(x, s)),
                                        filter_entries(e, f)),
                cnd_levels)

    # Shuffled distances
    if isnothing(process_pool)
        allShuffRD = map(x -> shuffled_distance(e.data[:,e.selected_features],
                                                filter_entries(e, Filter(x, s)),
                                                filter_entries(e, f),
                                                nb_rep = nb_rep),
                         cnd_levels)
    else
        sendto(workers(), e=e,
                          s=s,
                          f=f,
                          nb_rep=nb_rep)
        allShuffRD = pmap(x -> shuffled_distance(e.data[:,e.selected_features],
                                            filter_entries(e, Filter(x, s)),
                                            filter_entries(e, f),
                                            nb_rep = nb_rep),
                     process_pool,
                     cnd_levels)
    end

    # Missing values might need to be handled explicitely
    @assert !any(ismissing.(allRD))

    # Compute the Robust Morphological Perturbation Value
    plateRMPV = DataFrame()
    plateRMPV.RMPV = adjust([mean(obs .< sim) for (obs, sim)
                in zip(allRD, allShuffRD)], BenjaminiHochberg())
    plateRMPV.Distance = allRD
    plateRMPV.Condition = cnd_levels

    return(plateRMPV)
end

# If reference value provided for reference
function robust_morphological_perturbation_value(e::AbstractExperiment,
                                                 s::Symbol,
                                                 ref;
                                                 nb_rep::Int64 = 250,
                                                 dist::Symbol = :RobHellinger,
                                                 process_pool = nothing,
                                                 r_seed = true)
    ref_filter = Filter(ref, s)
    return(robust_morphological_perturbation_value(e,
                                                   s,
                                                   ref_filter;
                                                   nb_rep=nb_rep,
                                                   dist=dist,
                                                   process_pool=process_pool,
                                                   r_seed=r_seed))
end


text:
using Printf, LinearAlgebra, Plots, ParallelStencil, ParallelStencil.FiniteDifferences3D

# Initialize Globals
if !@isdefined USE_GPU
    const USE_GPU = false
end

if !@isdefined BENCHMARK
    const BENCHMARK = true
end

if !@isdefined VISUALIZE
    const VISUALIZE = false
end

if !@isdefined STEADY
    const STEADY = false
end

# Initialize `ParallelStencil.jl`
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 3);
else
    @init_parallel_stencil(Threads, Float64, 3);
end

# Macro for shorter notation
macro inner(A) esc(:( $A[2:end-1,2:end-1,2:end-1] )) end

@parallel function compute_ResH!(H, Hold, ResH, _dt, D_dx², D_dy², D_dz²)
    @all(ResH) = -(@inn(H) - @inn(Hold)) * _dt + (@d2_xi(H)*D_dx² + @d2_yi(H)*D_dy² + @d2_zi(H)*D_dz²)
    return
end

@parallel function compute_ResH_steady!(H, ResH, D_dx², D_dy², D_dz²)
    @all(ResH) = @d2_xi(H)*D_dx² + @d2_yi(H)*D_dy² + @d2_zi(H)*D_dz²
    return
end

@parallel function compute_dHdt!(ResH, dHdt, dHdt2, damp)
    @all(dHdt2) = @all(ResH) + damp * @all(dHdt)
    return
end

@parallel function compute_H!(H2, H, dHdt, dτ)
    @inn(H2) = @inn(H) + dτ * @all(dHdt)
    return
end

"""
    diffusion_3D(res[,tol])

Runs the 3D dualtime diffusion simulation

# Arguments
- `res::Int`: xyz-resolutions of the simulation
- `tol::Float64`(optional): tolerance for convergence - default: `1e-8`

# Returns (Depending on globals)
- If BENCHMARK is true ->
    `T_eff`: Effective memory throughput [GB/s]
- If STEADY is true and BENCHMARK is false ->
    `ittot`, `xc` and `H`: Total number of iterations needed to reach steadystate, the global x-coordinate Array and the global solution array
- If STEADY and BENCHMARK are false ->
    `xc` and `H`: The global x-coordinate Array and the global solution array
"""
@views function diffusion_3D(res::Int; tol=1e-8)
    # Physics
    lx, ly, lz = 10.0, 10.0, 10.0 # domain size
    D     = 1.0                   # diffusion coefficient
    ttot  = 1.0                   # total simulation time
    dt    = 0.2                   # physical time step

    # Numerics
    nx, ny, nz = res, res, res
    itMax      = 1e4              # max number of iterations
    nout       = 10               # tol check

    # Derived numerics
    dx, dy, dz = lx/nx, ly/ny, lz/nz            # cell sizes
    dτ         = (1.0/(dx^2/D/6.1) + 1.0/dt)^-1 # iterative timestep
    damp       = 1-29/nx                        # damping (this is a tuning parameter, dependent on e.g. grid resolution)
    xc         = LinRange(dx/2, lx-dx/2, nx)    #
    yc         = LinRange(dy/2, ly-dy/2, ny)
    zc         = LinRange(dz/2, lz-dz/2, nz)
    _dt        = 1.0/dt
    D_dx²      = D/dx^2
    D_dy²      = D/dy^2
    D_dz²      = D/dz^2

    # Array allocation
    H     = @zeros(nx  , ny  , nz  )
    H2    = @zeros(nx  , ny  , nz  )
    Hᵗ    = @zeros(nx  , ny  , nz  )
    dHdt  = @zeros(nx-2, ny-2, nz-2)
    dHdt2 = @zeros(nx-2, ny-2, nz-2)
    ResH  = @zeros(nx-2, ny-2, nz-2)

    # Initial condition
    @inner(H) .= @inner(Data.Array([2.0 * exp(-(xc[ix] - 0.5*lx)^2 -(yc[iy] - 0.5*ly)^2 -(zc[iz] - 0.5*lz)^2) for ix=1:nx, iy=1:ny, iz=1:nz]))
    H2        .= H
    Hᵗ        .= H

    # Benchmark variables
    @static if BENCHMARK
        warmup = 10         # Number of warmup-iterations
        t_tic  = 0.0        # Holds start time
        itMax  = 1e4/nx     # scale max number of iterations down my grid size
    end

    # Run to Steady State variables
    @static if STEADY
        damp = 1-4/nx
    end

    # Animation object for visualization
    @static if VISUALIZE
        gr()
        ENV["GKSwstype"] = "nul"
        anim = Animation()
        H_v = zeros(nx, ny, nz)
        init_extrema = extrema(H)
    end

    # initial physical time loop conditions
    t_tic = 0.0
    t = 0.0
    it_t = 0
    ittot = 0

    @static if !STEADY
        println("Running with dual-time")

        # Physical time loop
        while t<ttot

            # initial pseudo-transient loop conditions
            it_τ = 0
            err = 2*tol

            # Pseudo-transient iteration
            while err>tol && it_τ<itMax

                # Step 1
                @parallel compute_ResH!(H, Hᵗ, ResH, _dt, D_dx², D_dy², D_dz²)
                # Step 2
                @parallel compute_dHdt!(ResH, dHdt, dHdt2, damp)
                dHdt, dHdt2 = dHdt2, dHdt # Pointer Swap
                # Step 3
                @parallel compute_H!(H2, H, dHdt, dτ)
                H, H2 = H2, H # Pointer Swap

                # Calculate L2 norm
                if (it_τ % nout == 0)
                    err = norm(ResH)/sqrt(length(ResH))
                end

                # Start the clock after warmup-iterations
                @static if BENCHMARK
                    if (it_t == 0 && it_τ == warmup+1)
                        t_tic = Base.time();
                    end
                end

                # Render a slice of the 3D diffusion-map and save as an animation frame
                @static if VISUALIZE
                    if (it_τ % nout == 0)
                        H_v .= H
                        opts = (aspect_ratio=1, xlims=extrema(xc), ylims=extrema(yc), clims=init_extrema, xlabel="lx", ylabel="ly", title="3D Diffusion (it_t=$it_t, it_τ=$it_τ)")
                        heatmap(xc, yc, H_v[:,:,nz÷2]; opts...)
                        frame(anim)
                    end
                end

                it_τ += 1
            end

            # update physical time step
            ittot += it_τ; it_t += 1; t += dt
            Hᵗ .= H
            if isnan(err) error("NaN") end
        end
    else
        println("Running in pseudo-time to steady-state")

        # initial pseudo-transient loop conditions
        it_τ = 0
        err = 2*tol

        # Pseudo-transient iteration
        while err>tol && it_τ<itMax
            # Step 1
            @parallel compute_ResH_steady!(H, ResH, D_dx², D_dy², D_dz²)
            # Step 2
            @parallel compute_dHdt!(ResH, dHdt, dHdt2, damp)
            dHdt, dHdt2 = dHdt2, dHdt # Pointer Swap
            # Step 3
            @parallel compute_H!(H2, H, dHdt, dτ)
            H, H2 = H2, H # Pointer Swap

            # Calculate L2 norm (error)
            if (it_τ % nout == 0)
                err = norm(ResH)/sqrt(length(ResH))
            end
            it_τ += 1
        end

        # update physical time step
        ittot += it_τ; it_t += 1
        if isnan(err) error("NaN") end
    end

    # Calculate performance
    @static if BENCHMARK
        t_toc = Base.time() - t_tic                            # Stop the clock
        reads = length(Hᵗ)                                     # Read Only Memory Access: Hᵗ
        updates = length(H) + length(dHdt)                     # Update Memory access: H and dHdt
        A_eff = 1e-9 * (2 * updates + reads) * sizeof(Float64) # Effective main memory access per iteration [GB]
        t_it  = t_toc/(ittot-warmup)                           # Execution time per iteration [s]
        T_eff = A_eff/t_it                                     # Effective memory throughput [GB/s]
        println("time(sec)=$t_toc T_eff=$T_eff")
    end
    @printf("Ttime steps = %d, nx = %d, iterations tot = %d \n", it_t, nx, ittot)

    # Save the animation as a GIF
    @static if VISUALIZE
        gif(anim, "diffusion3D_xpu.gif", fps = 15)
    end

    # Return T_eff if BENCHMARK and Result otherwise
    @static if BENCHMARK
        return T_eff
    elseif STEADY
        return ittot, xc, Array(H)
    else
        return xc, Array(H)
    end
end