text:
# ### m10.3t.jl

using Pkg, DrWatson

@quickactivate "SR2TuringPluto"
using Turing
using StatisticalRethinking
Turing.setprogress!(false)

delim=';'
df = CSV.read(sr_datadir("chimpanzees.csv"), DataFrame; delim);

# pulled_left, condition, prosoc_left
@model function ppl10_3(y, x₁, x₂)
    α ~ Normal(0, 10)
    βp ~ Normal(0, 10)
    βpC ~ Normal(0, 10)

    logits = α .+ (βp .+ βpC * x₁) .* x₂
    y .~ BinomialLogit.(1, logits)
end;

data = (df.pulled_left, df.condition, df.prosoc_left)
m10_3t = ppl10_3(data...)
nchains = 4; sampler = NUTS(0.65); nsamples=2000
chns10_3t = mapreduce(c -> sample(m10_3t, sampler, nsamples), chainscat, 1:nchains)

# Rethinking result

m_10_03t_result = "
      Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
 a    0.05   0.13      -0.15       0.25  3284    1
 bp   0.62   0.22       0.28       0.98  3032    1
 bpC -0.11   0.26      -0.53       0.29  3184    1
";

chains = sample(m10_3t, NUTS(1000, .65), MCMCThreads(), 1000, 4)

# method for MCMCChains
function pointwise_loglikes(chain::Chains, data, ll_fun)
    samples = Array(Chains(chain, :parameters).value)
    pointwise_loglikes(samples, data, ll_fun)
end

# generic method for arrays
function pointwise_loglikes(samples::Array{Float64,3}, data, ll_fun)
    n_data = length(data)
    n_samples, n_chains = size(samples)[[1,3]]
    pointwise_lls = fill(0.0, n_data, n_samples, n_chains)
    for c in 1:n_chains
        for s in 1:n_samples
            for d in 1:n_data
                pointwise_lls[d,s,c] = ll_fun(samples[s,:,c], data[d])
            end
        end
    end
    return pointwise_lls
end

function compute_loo(psis_output, pointwise_lls)
    dims = size(pointwise_lls)
    lwp = deepcopy(pointwise_lls)
    lwp += psis_output.weights;
    lwpt = reshape(lwp, dims[1], dims[2] * dims[3])';
    loos = reshape(logsumexp(lwpt; dims=1), size(lwpt, 2));
    return sum(loos)
end

# compute the pointwise log likelihoods where indices correspond to [data, sample, chain]
pointwise_lls = pointwise_loglikes(chains, df.pulled_left, (p,d)->logpdf(Normal(p...), d))

# compute the psis object
psis_output = psis(pointwise_lls)

# return loo based on Rob's example
loo = compute_loo(psis_output, pointwise_lls)


text:
using HTTP

include("../src/network.jl")

@testset "network" begin

    r = request(url="https://api.crossref.org/works",
                query="Heisenberg",
                filter=nothing,
                offset=nothing,
                limit=1,
                sample=nothing,
                sort=nothing,
                order=nothing,
                facet=nothing,
                select=nothing,
                cursor=nothing)

    @test r.status == 200

    let err = nothing
        try
            request(url="https://api.crossref.org/works", offset=11000)
        catch err
        end

        @test err isa Exception
        @test sprint(showerror, err) == "Offset; Max: 10000"

        try
            request(url="https://api.crossref.org/works", sample=500)
        catch err
        end

        @test err isa Exception
        @test sprint(showerror, err) == "Sample size; Max: 100"

    end

end

ext: jl

max_stars_repo_name: naustica/Pitaya

max_issues_repo_name: naustica/Pitaya

max_forks_repo_name: naustica/Pitaya


text:
"""
Each [`Interval`](@ref) consists of two `IntervalPoint`s, one of [`RayOrigin`](@ref), [`Intersection`](@ref) or [`Infinity`](@ref).
"""
abstract type IntervalPoint{T<:Real} end
Base.eltype(::IntervalPoint{T}) where {T<:Real} = T

abstract type FinitePoint{T} <: IntervalPoint{T} end
Base.eltype(::FinitePoint{T}) where {T<:Real} = T

"""
    Intersection{T,N} <: IntervalPoint{T}

Represents the point at which an [`Ray`](@ref) hits a [`Surface`](@ref).
This consists of the distance along the ray, the intersection point in world space, the normal in world space, the UV on the surface and the [`OpticalInterface`](@ref) hit.

Has the following accessor methods:
```julia
point(a::Intersection{T,N}) -> SVector{N,T}
normal(a::Intersection{T,N}) -> SVector{N,T}
uv(a::Intersection{T,N}) -> SVector{2,T}
u(a::Intersection{T,N}) -> T
v(a::Intersection{T,N}) -> T
α(a::Intersection{T,N}) -> T
interface(a::Intersection{T,N}) -> OpticalInterface{T}
flippednormal(a::Intersection{T,N}) -> Bool
```
"""
struct Intersection{T,N} <: FinitePoint{T}
    α::T # value of ray parameter at the point of intersection
    point::SVector{N,T}
    normal::SVector{N,T}
    u::T
    v::T
    interface::AllOpticalInterfaces{T} # returns a union of all OpticalInterface subtypes - can't have an abstract type here as it results in allocations
    flippednormal::Bool

    function Intersection(α::T, point::SVector{N,T}, normal::SVector{N,T}, u::T, v::T, interface::AllOpticalInterfaces{T}; flippednormal = false) where {T<:Real,N}
        new{T,N}(α, point, normalize(normal), u, v, interface, flippednormal)
    end

    function Intersection(α::T, point::AbstractVector{T}, normal::AbstractVector{T}, u::T, v::T, interface::AllOpticalInterfaces{T}; flippednormal = false) where {T<:Real}
        @assert length(point) == length(normal)
        N = length(point)
        new{T,N}(α, SVector{N,T}(point), SVector{N,T}(normal), u, v, interface, flippednormal)
    end
end
export Intersection

point(a::Intersection{T,N}) where {T<:Real,N} = a.point
normal(a::Intersection{T,N}) where {T<:Real,N} = a.normal
uv(a::Intersection{T,N}) where {T<:Real,N} = SVector{2,T}(a.u, a.v)
u(a::Intersection{T,N}) where {T<:Real,N} = a.u
v(a::Intersection{T,N}) where {T<:Real,N} = a.v
α(a::Intersection{T,N}) where {T<:Real,N} = a.α
interface(a::Intersection{T,N}) where {T<:Real,N} = a.interface
flippednormal(a::Intersection{T,N}) where {T<:Real,N} = a.flippednormal

function Base.print(io::IO, a::Intersection{T,N}) where {T<:Real,N}
    println(io, Intersection{T,N})
    println(io, "α \t$(α(a))")
    println(io, "Point \t$(point(a))")
    println(io, "Normal \t$(normal(a))")
    println(io, "normal is flipped? \t$(flippednormal(a))")
    println(io, "u \t$(u(a))")
    println(io, "v \t$(v(a))")
    println(io, "interface \t$(interface(a))")
end

"""
    reversenormal(a::Intersection{T,N})

Used by the CSG complement operator (i.e. [`csgdifference`](@ref)) to reverse the inside outside sense of the object.
"""
function reversenormal(a::Intersection{T,N}) where {T<:Real,N}
    return Intersection(α(a), point(a), -normal(a), u(a), v(a), interface(a), flippednormal = !flippednormal(a))
end

"""
    Infinity{T} <: IntervalPoint{T}

Point representing ∞ within an [`Interval`](@ref).

```julia
Infinity(T = Float64)
Infinity{T}()
```
"""
struct Infinity{T} <: IntervalPoint{T}
    Infinity(::Type{T} = Float64) where {T<:Real} = new{T}()
    Infinity{T}() where {T<:Real} = new{T}()
end
"""
    RayOrigin{T} <: IntervalPoint{T}

Point representing 0 within an [`Interval`](@ref), i.e. the start of the ray.

```julia
RayOrigin(T = Float64)
RayOrigin{T}()
```
"""
struct RayOrigin{T} <: FinitePoint{T}
    RayOrigin(::Type{T} = Float64) where {T<:Real} = new{T}()
    RayOrigin{T}() where {T<:Real} = new{T}()
end
export Infinity, RayOrigin

Base.:(==)(::RayOrigin{P}, ::RayOrigin{P}) where {P<:Real} = true
Base.:(==)(::Infinity{P}, ::Infinity{P}) where {P<:Real} = true
Base.:(==)(a::Intersection{P,N}, b::Intersection{P,N}) where {P<:Real,N} = α(a) == α(b)
Base.:(==)(::IntervalPoint{P}, ::IntervalPoint{P}) where {P<:Real} = false

Base.:(<=)(::RayOrigin{P}, ::RayOrigin{P}) where {P<:Real} = true
Base.:(<=)(::Infinity{P}, ::Infinity{P}) where {P<:Real} = true
Base.:(<=)(a::IntervalPoint{P}, b::Infinity{P}) where {P<:Real} = true
Base.:(<=)(a::Intersection{P,N}, b::Intersection{P,N}) where {P<:Real,N} = α(a) <= α(b)
Base.:(<=)(a::RayOrigin{P}, b::Infinity{P}) where {P<:Real} = true
Base.:(<=)(a::RayOrigin{P}, b::IntervalPoint{P}) where {P<:Real} = true
Base.:(<=)(::IntervalPoint{P}, ::IntervalPoint{P}) where {P<:Real} = false

Base.:(>=)(::RayOrigin{P}, ::RayOrigin{P}) where {P<:Real} = true
Base.:(>=)(::Infinity{P}, ::Infinity{P}) where {P<:Real} = true
Base.:(>=)(a::IntervalPoint{P}, b::Infinity{P}) where {P<:Real} = false
Base.:(>=)(a::Intersection{P,N}, b::Intersection{P,N}) where {P<:Real,N} = α(a) >= α(b)
Base.:(>=)(a::RayOrigin{P}, b::Infinity{P}) where {P<:Real} = false
Base.:(>=)(a::RayOrigin{P}, b::IntervalPoint{P}) where {P<:Real} = false
Base.:(>=)(::IntervalPoint{P}, ::IntervalPoint{P}) where {P<:Real} = false

Base.:(<)(::Infinity{P}, ::Infinity{P}) where {P<:Real} = false
Base.:(<)(::IntervalPoint{P}, ::Infinity{P}) where {P<:Real} = true
Base.:(<)(a::Intersection{P,N}, b::Intersection{P,N}) where {P<:Real,N} = α(a) < α(b)

Base.:(>)(::Infinity{P}, ::Infinity{P}) where {P<:Real} = false
Base.:(>)(::IntervalPoint{P}, ::Infinity{P}) where {P<:Real} = false
Base.:(>)(a::Intersection{P,N}, b::Intersection{P,N}) where {P<:Real,N} = α(a) > α(b)

Base.isless(a::IntervalPoint{P}, b::IntervalPoint{P}) where {P<:Real} = α(a) < α(b)

Base.eltype(::Intersection{T,N}) where {T<:Real,N} = T

"""
    isinfinity(a) -> Bool

Returns true if `a` is [`Infinity`](@ref). In performance critical contexts use `a isa Infinity{T}`.

"""
isinfinity(::Infinity) = true
isinfinity(::Any) = false

α(::Infinity{T}) where {T<:AbstractFloat} = typemax(T)
Base.eltype(::Infinity{T}) where {T<:Real} = T

"""
    israyorigin(a) -> Bool

Returns true if `a` is [`RayOrigin`](@ref). In performance critical contexts use `a isa RayOrigin{T}`.

"""
israyorigin(::Any) = false
israyorigin(::RayOrigin) = true

α(::RayOrigin{T}) where {T<:Real} = zero(T)
Base.eltype(::RayOrigin{T}) where {T<:Real} = T


text:
using Images
using TikzGraphs
using LightGraphs
import TikzPictures
using LaTeXStrings
using CGP

function draw_graph(c::Chromosome, filename::String)
    ids = 1:length(c.nodes)
    active = map(i->i.active, c.nodes)
    active_ids = ids[active]
    n_active = length(active_ids)
    g = DiGraph(n_active + c.nin + c.nout)
    names = Array{LaTeXStrings.LaTeXString}(n_active + c.nin + c.nout)

    names[1] = L"i_{red}"
    names[2] = L"i_{green}"
    names[3] = L"i_{blue}"

    for i in 1:n_active
        nid = active_ids[i]
        node = c.nin+i
        for ci in 1:2
            o = c.nodes[nid].connections[ci]
            input_node = o
            if o > c.nin
                input_node = c.nin+findfirst(active_ids .== o)
            end
            add_edge!(g, input_node, node)
        end
        names[node] = latexstring(
            @sprintf("%0.2f", c.nodes[nid].p),
            "f_{", replace(
            split(string(c.nodes[nid].f),".")[3][3:end],
            "_", ""), "}")
    end

    for o in 1:length(c.outputs)
        output_node = n_active + c.nin + o
        i = c.outputs[o]
        input_node = i
        if i > c.nin
            input_node = c.nin+findfirst(active_ids .== i)
        end
        add_edge!(g, input_node, output_node)
        names[output_node] = latexstring("o_{", o,"}")
    end

    t = TikzGraphs.plot(g, Layouts.Layered(), names)

    TikzPictures.save(TikzPictures.PDF(filename), t)
    t
end


text:
# Here are tests which test what interface the solvers require.

################################################################################
# This is to test a scalar-like state variable
# (due to @acroy: https://gist.github.com/acroy/28be4f2384d01f38e577)

import Base: +, -, *, /
using LinearAlgebra

const delta0 = 0.
const V0 = 1.
const g0 = 0.

# define custom type ...
struct CompSol <: AbstractMatrix{ComplexF64}
  rho::Matrix{ComplexF64}
  x::Float64
  p::Float64

  CompSol(r,x,p) = new(copy(r),x,p)
end

Base.getindex(A::CompSol,i...) = A.rho[i...]
Base.setindex!(A,x,i...) = setindex!(A.rho,x,i...)

# ... which has to support the following operations
# to work with odeX
LinearAlgebra.norm(y::CompSol, p::Float64) = maximum([LinearAlgebra.norm(y.rho, p) abs(y.x) abs(y.p)])
LinearAlgebra.norm(y::CompSol) = LinearAlgebra.norm(y::CompSol, 2.0)

+(y1::CompSol, y2::CompSol) = CompSol(y1.rho+y2.rho, y1.x+y2.x, y1.p+y2.p)
-(y1::CompSol, y2::CompSol) = CompSol(y1.rho-y2.rho, y1.x-y2.x, y1.p-y2.p)
*(y1::CompSol, s::Real) = CompSol(y1.rho*s, y1.x*s, y1.p*s)
*(s::Real, y1::CompSol) = y1*s
/(y1::CompSol, s::Real) = CompSol(y1.rho/s, y1.x/s, y1.p/s)

### new for PR #68
Base.abs(y::CompSol) = LinearAlgebra.norm(y, 2.) # TODO not needed anymore once https://github.com/JuliaLang/julia/pull/11043 is in current stable julia
Base.zero(::Type{CompSol}) = CompSol(complex(zeros(2,2)), 0., 0.)
ODE.isoutofdomain(y::CompSol) = any(isnan, vcat(y.rho[:], y.x, y.p))

################################################################################

# define RHSs of differential equations
# delta, V and g are parameters
function rhs(t, y, delta, V, g)
  H = [[-delta/2 V]; [V delta/2]]

  rho_dot = -im*H*y.rho + im*y.rho*H
  x_dot = y.p
  p_dot = -y.x

  return CompSol( rho_dot, x_dot, p_dot)
end

# inital conditons
rho0 = zeros(2,2);
rho0[1,1]=1.;
y0 = CompSol(complex(rho0), 2., 1.)

# solve ODEs
endt = 2*pi;

t,y1 = ODE.ode45((t,y)->rhs(t, y, delta0, V0, g0), y0, [0., endt]) # used as reference
print("Testing interface for scalar-like state... ")
for solver in solvers
    # these only work with some Array-like interface defined:
    if solver in [ODE.ode23s, ODE.ode4s_s, ODE.ode4s_kr]
        continue
    end
    t,y2 = solver((t,y)->rhs(t, y, delta0, V0, g0), y0, range(0., stop=endt, length=500))
    @test LinearAlgebra.norm(y1[end]-y2[end])<0.15
end
println("ok.")

################################################################################
# TODO: test a vector-like state variable, i.e. one which can be indexed.


text:
import Effects: effects
import Effects:_reference_grid
import Effects:typify
import Effects.typify
import Effects:_symequal

using Effects

"""
effects(design::AbstractDict, model::UnfoldModel;typical=mean)

Calculates marginal effects for all term-combinations in `design`.

 Implementation based on Effects Package; likely could repackage in UnfoldEffects; somebody wants to do it? This would make it easier to cross-maintain it to changes/bugfixes in the Effects.jl Package
 `design` is a Dictionary containing those predictors (als keys) with levels (as values), that you want to evaluate. The `typical` refers to the value, that other predictors not in the Dictionary should take on.

 # Example
 ```julia-repl
 julia> f = @formula 0 ~ categoricalA + continuousA + continuousB
 julia> uf = fit(UnfoldModel,(Any=>(f,times)),data,events)
 julia> d = Dict(:categoricalA=>["levelA","levelB"],:continuousB=>[-2,0,2])
 julia> effects(d,uf)
```
 will result in 6 predicted values: A/-2, A/0, A/2, B/-2, B/0, B/2.
"""
Effects.typify(reference_grid,form::Matrix,X;kwargs...) = typify.(Ref(reference_grid),form,Ref(X);kwargs...)

function effects(design::AbstractDict, model::UnfoldModel;typical=mean)
    reference_grid = _reference_grid(design)
    form = Unfold.formula(model) # get formula
    form_typical = typify(reference_grid, form, modelmatrix(model); typical=typical) # replace non-specified fields with "constants"

    eff = yhat(model,form_typical,reference_grid)

    # because coefficients are 2D/3D arry, we have to cast it correctly to one big dataframe
    if isa(eff,Tuple)
        # TimeContinuous Model, we also get back other things like times & fromWhereToWhere a BasisFunction goes
        if typeof(form) <: FormulaTerm
            form = [form]
        end
        # figure out the basisname
        bnames = [[form[ix].rhs.basisfunction.name] for ix in 1:length(form)]
        # repeat it as much as necessary
        bnames = repeat.(bnames,[e.stop+e.step-1 for e in eff[1]])

        result = DataFrame(cast_referenceGrid(reference_grid,eff[3],eff[2] ;basisname=vcat(bnames...)))
        select!(result,Not(:latency)) # remove the latency column if it was added
    else
        # normal mass univariate model
        result = DataFrame(cast_referenceGrid(reference_grid,eff,times(model)[1] ))
    end

return result
end

function cast_referenceGrid(r,eff,times;basisname=nothing)
    nchan = size(eff, 2) # correct
    neff = size(r,1) # how many effects requested
    neffCol = size(r,2) # how many predictors
    ncols = size(eff,1) ÷ neff # typically ntimes

    #@show size(basisname)
    #@show size(times)
    @show nchan
    @show ncols
    @show neff
    @show neffCol
    @show r
    # replicate
    # for each predictor in r (reference grid), we need this at the bottom
    coefs_rep = Array{Float64}(undef,nchan,ncols,neff,neffCol)
    for k = 1:neffCol
        # repeat it for nchan + ncols
        coefs_rep[:,:,:,k] = permutedims(repeat(r[:,k], outer = [1, nchan, ncols]), [2, 3, 1])
    end

    # often the "times" vector
    if length(times)  == neff*ncols
        # in case we have timeexpanded, times is already in long format and doesnt need to be repeated for each coefficient
        colnames_basis_rep = permutedims(repeat(times, 1, nchan, 1), [2 1 3])
    else
        colnames_basis_rep = permutedims(repeat(times, 1, nchan, neff), [2 1 3])
    end

    # for multiple channels
    chan_rep = repeat(1:nchan, 1, ncols, neff)

    # for mass univariate there is no basisname
    if isnothing(basisname)
        basisname = fill(nothing,ncols)
        basisname_rep = permutedims(repeat(basisname,1,nchan,neff),[2,1,3])
    else
        basisname_rep = permutedims(repeat(basisname,1,nchan,1),[2,1,3])
    end

    # repeat it!

    #eff is already 2D, but we have to flip it.
    @show (size(eff'), size(colnames_basis_rep),size(chan_rep), size(basisname_rep))
    #@assert size(colnames_basis_rep) == size(chan_rep) == size(basisname_rep)
    result = Dict(   :yhat => linearize(eff'),
            :time => linearize(colnames_basis_rep),
            :channel => linearize(chan_rep),
            :basisname =>linearize(basisname_rep))

    for k = 1:neffCol
            push!(result,Symbol(names(r)[k]) =>linearize(coefs_rep[:,:,:,k]))
    end

    return result
end

_symequal(t1::AbstractTerm,t2::Unfold.TimeExpandedTerm) = _symequal(t1,t2.term)


text:
"""
peelAndBound(
    model::SequencingModel,
    maxWidth::U;

    orderingType::OrderingOption=edge,
    bestKnownSolution::Union{Array{U},Nothing}=nothing,
    bestKnownValue::Union{T,Nothing}=nothing,
    peelSetting::PeelSetting=peelf,
    rrbSetting::Union{RRBSetting,Nothing}=nothing,
    rrbPacket=nothing,
    timeLimit::Union{Nothing, W}=nothing,
    timeElapsed::X=0.0,
    loggingOn::Bool=true,
    debugOn::Bool=false,
    fileName::Union{Nothing,String}=nothing
)where{T<:Real,U<:Int, W<:Real, X<:Real}

Solve a sequencing problem to optimality using a peel and bound scheme
Returns bestKnownSolution (Vector{Int}), bestKnownValue (Real), bound (Real)

# Arguments
- `model::SequencingModel`: The sequencing problem
- `maxWidth::Int`: The max width allowed in a decision diagram layer
- `orderingType::OrderingOption`:Which heuristic to use to sort the edges, defaults to 'edge'
- `bestKnownSolution::Union{Vector{Int},Nothing}`: An optional parameter that gives the DD a starting solution
- `bestKnownValue::Union{T,Nothing}`: An optional parameter, The value of the best known feasible solution
- `peelSetting::PeelSetting`:Which heuristic to use to decide where to start the peel from, defaults to 'peelf'
- `rrbSetting::Union{RRBSetting,Nothing}`: An optional parameter that chooses an rough relaxed bounding method
- `rrbPacket`: Optional data to be passed to the rrbFunction
- `timeLimit::Union{Nothing, W}`: An optional parameter that adds a time limit
- `timeElapsed::X`: An optional parameter denoting how much time passed already
- `loggingOn::Bool`: An optional parameter that can be used to turn off logging to the console and local file
- `debugOn::Bool`: An optional parameter that can be used to turn on debug statements in other functions
- `fileName::Union{Nothing,String}`: The file to write results to
"""
function peelAndBound(
    model::SequencingModel,
    maxWidth::U;
    orderingType::OrderingOption=edge,
    bestKnownSolution::Union{Array{U},Nothing}=nothing,
    bestKnownValue::Union{T,Nothing}=nothing,
    peelSetting::PeelSetting=peelf,
    rrbSetting::Union{RRBSetting,Nothing}=nothing,
    rrbPacket=nothing,
    timeLimit::Union{Nothing, W}=nothing,
    timeElapsed::X=0.0,
    loggingOn::Bool=true,
    debugOn::Bool=false,
    fileName::Union{Nothing,String}=nothing
)where{T<:Real,U<:Int, W<:Real, X<:Real}
    startTime = time()
    bestBound = nothing
    timedOut() = isTimedOut(timeLimit,startTime, timeElapsed)
    cleanRestrictedDD(node) = restrictedSequencingDD(model, maxWidth,rootNode = node, debugOn=debugOn,bestKnownSolution=bestKnownSolution, bestKnownValue=bestKnownValue, rrbSetting=rrbSetting, rrbPacket=rrbPacket)
    cleanRelaxedDD(dd) = relaxedSequencingDD(model, maxWidth,splitOrder, dd = dd, debugOn=debugOn,bestKnownSolution=bestKnownSolution, bestKnownValue=bestKnownValue, rrbSetting=rrbSetting, rrbPacket=rrbPacket, timeLimit=timeLimit, timeElapsed=timeElapsed, startTime=startTime)
    splitOrder = findSplitOrdering(model,orderingType)
    cleanPeel!(relaxedDD, relaxedSolution, frontierNode, fNodeIndex) = peelDD!(model, maxWidth, relaxedDD, relaxedSolution, frontierNode, fNodeIndex, splitOrder,bestKnownValue,rrbSetting,rrbPacket,loggingOn,debugOn)
    cleanPeelAndRemove!(dd, frontierNode, fNodeIndex, bestKnownValue) = peelAndRemove!(model,dd,frontierNode,fNodeIndex,bestKnownValue,rrbSetting,rrbPacket,loggingOn,debugOn)

    #get initial bounds
    if isnothing(bestKnownSolution)
        bestKnownSolution, bestKnownValue, exact = restrictedSequencingDD(model,maxWidth,debugOn=debugOn)
        if exact
            if loggingOn
                finishPeelAndBound(bestKnownSolution, bestKnownValue,startTime=startTime, timeElapsed=timeElapsed, fileName=fileName)
            end
            return bestKnownSolution, bestKnownValue, bestKnownValue
        end
    end
    if timedOut()
        return bestKnownSolution, bestKnownValue, 0
    end
    #improve initial bounds using rrb
    if !isnothing(rrbSetting)
        bestKnownSolution, bestKnownValue, exact = cleanRestrictedDD(nothing)
        #if solution is found finish
        if exact
            if loggingOn
                finishPeelAndBound(bestKnownSolution, bestKnownValue,startTime=startTime, timeElapsed=timeElapsed, fileName=fileName)
            end
            return bestKnownSolution, bestKnownValue, bestKnownValue
        end
        if timedOut()
            return bestKnownSolution, bestKnownValue, 0
        end
    end
    if loggingOn
        logSequencingSolution(bestKnownSolution, bestKnownValue,firstLine="Initial Solution:", fileName=fileName)
        logRuntime(startTime, timeElapsed, fileName=fileName)
    end
    #generate initial relaxed DD
    relaxedDD, relaxedBound, relaxedSolution, bestKnownSolution, bestKnownValue, exact = cleanRelaxedDD(nothing)

    if exact || relaxedBound == bestKnownValue
        if loggingOn
            finishPeelAndBound(bestKnownSolution, bestKnownValue,startTime=startTime, timeElapsed=timeElapsed, fileName=fileName)
        end
        return bestKnownSolution, bestKnownValue, bestKnownValue
    elseif timedOut()
        return bestKnownSolution, bestKnownValue, 0
    else
        bestBound = relaxedBound
        if loggingOn
            logOptimalityGap(bestBound, bestKnownValue, fileName=fileName)
            logRuntime(startTime, timeElapsed, fileName=fileName)
        end
    end

    #initialize queue
    ddType = Vector{Vector{RelaxedSequencingNode}}
    valueType = typeof(bestKnownValue)
    queue = Vector{Tuple{ddType, valueType, Vector{Int}, Bool}}()
    push!(queue, (relaxedDD, bestBound, relaxedSolution, false))

    if loggingOn
        logQueueLength(length(queue),fileName=fileName)
    end
    #process queue
    while !isempty(queue)
        currentValue = bestKnownValue
        currentBlock = pop!(queue)
        if debugOn
            if isempty(queue)
                println("In Queue: ",length(queue)+1, "  From: ",currentBlock[2], " - ",currentBlock[2], "\n")
            else
                println("In Queue: ",length(queue)+1, "  From: ",currentBlock[2], " - ",queue[1][2], "\n")
            end
        end

        if peelSetting == peelf
            frontierNode, fNodeIndex = getFrontierNode(currentBlock[1], currentBlock[3])
        else
            frontierNode, fNodeIndex = getLastExactNode(currentBlock[1], currentBlock[3])
        end
        filterArcs!(frontierNode.arcsOut,model;bestKnownValue=bestKnownValue, rrbSetting=rrbSetting, rrbPacket=rrbPacket)

        bestKnownSolution, bestKnownValue, restrictedExact = cleanRestrictedDD(frontierNode)
        #time limit block
        if timedOut()
            return bestKnownSolution, bestKnownValue, bestBound
        end

        if restrictedExact
            #delete the node
            ddInfo, bestPeeledSolution, bestPeeledValue = cleanPeelAndRemove!(currentBlock[1], frontierNode, fNodeIndex, bestKnownValue)

            if !ddInfo[4] && isBetterSolutionValue(model,ddInfo[2], bestKnownValue)
                addDDToQueue!(queue, ddInfo, model)
            elseif isBetterSolutionValue(model, bestPeeledValue, bestKnownValue)
                bestKnownValue = bestPeeledValue
                bestKnownSolution = bestPeeledSolution
            end
        else
            #peel the dd
            ddInfo, newDDInfo, bestPeeledSolution, bestPeeledValue = cleanPeel!(currentBlock[1], currentBlock[3], frontierNode, fNodeIndex)

            if isBetterSolutionValue(model, bestPeeledValue, bestKnownValue)
                bestKnownValue = bestPeeledValue
                bestKnownSolution = bestPeeledSolution
            end

            #add the DDs to the queue
            if !ddInfo[4] && isBetterSolutionValue(model,ddInfo[2], bestKnownValue)
                addDDToQueue!(queue, ddInfo, model)
            end
            if !newDDInfo[4] && isBetterSolutionValue(model,newDDInfo[2], bestKnownValue)
                addDDToQueue!(queue, newDDInfo, model)
            end
        end

        while !isempty(queue) && isBetterSolutionValue(model, bestKnownValue, queue[1][2])
            deleteat!(queue, 1)
        end
        if !isempty(queue)
            if loggingOn && isBetterSolutionValue(model, bestKnownValue, currentValue)
                logSequencingSolution(bestKnownSolution, bestKnownValue,firstLine="Improved Solution:", fileName=fileName)
                logOptimalityGap(last(queue)[2], bestKnownValue, fileName=fileName)
                logRuntime(startTime, timeElapsed, fileName=fileName)
                logQueueLength(length(queue),fileName=fileName)
            elseif loggingOn && isBetterSolutionValue(model, bestBound, last(queue)[2])
                logOptimalityGap(last(queue)[2], bestKnownValue, fileName=fileName)
                logRuntime(startTime, timeElapsed, fileName=fileName)
                logQueueLength(length(queue),fileName=fileName)
            end
            bestBound = last(queue)[2]
        end

        #time limit block
        if timedOut()
            return bestKnownSolution, bestKnownValue, bestBound
        end
    end

    if loggingOn
        finishPeelAndBound(bestKnownSolution, bestKnownValue, startTime=startTime, timeElapsed=timeElapsed, fileName=fileName)
    end
    return bestKnownSolution, bestKnownValue, bestKnownValue
end


text:
export SingularPoint,
    LogPointSingularity,
    AlgebraicPointSingularity,
    DiagonalSingularity

"The supertype of all kinds of singularities of an integrand."
abstract type Singularity end

"The integrand has no (known) singularity."
struct NoSingularity <: Singularity end

"Supertype of all kinds of point singularities."
abstract type PointSingularity <: Singularity end

point(s::PointSingularity) = s.point

"An unspecified singularity occurring at a specific point."
struct SingularPoint{T} <: PointSingularity
    point   ::  T
end

"A logarithmic singularity at a specific point."
struct LogPointSingularity{T} <: PointSingularity
    point   ::  T
end

"An algebraic singularity at a specific point."
struct AlgebraicPointSingularity{O,T} <: PointSingularity
    point   ::  T
    order   ::  O
end

order(s::AlgebraicPointSingularity) = s.order

"Supertype of singularities along a curve."
abstract type CurveSingularity <: Singularity end

"For 2D integrands, a singularity along the line `x=y`."
struct DiagonalSingularity <: CurveSingularity
end

"""
A list of breakpoints of the integrand, for example because the integrand
is piecewise smooth.
"""
struct Waypoints{P} <: Singularity
    points  ::  P
end


text:
module StatsModels

using Tables
using StatsBase
using ShiftedArrays: lag, lead
using DataStructures
using DataAPI: levels

using SparseArrays
using LinearAlgebra

using Tables: ColumnTable

export
    #re-export from StatsBase:
    StatisticalModel,
    RegressionModel,

    @formula,
    ModelFrame,
    ModelMatrix,

    AbstractContrasts,
    EffectsCoding,
    DummyCoding,
    HelmertCoding,
    ContrastsCoding,

    coefnames,
    dropterm,
    setcontrasts!,

    AbstractTerm,
    ConstantTerm,
    Term,
    ContinuousTerm,
    CategoricalTerm,
    InteractionTerm,
    FormulaTerm,
    InterceptTerm,
    FunctionTerm,
    LagTerm,
    MatrixTerm,

    lag, lead, # Rexported from ShiftedArrays

    term,
    terms,
    drop_term,
    schema,
    concrete_term,
    apply_schema,
    width,
    modelcols,
    modelmatrix,
    response

include("traits.jl")
include("contrasts.jl")
include("terms.jl")
include("schema.jl")
include("temporal_terms.jl")
include("formula.jl")
include("modelframe.jl")
include("statsmodel.jl")

end # module StatsModels


text:
@testset "798.smallest-rotation-with-highest-score.jl" begin
    @test best_rotation([2, 3, 1, 4, 0]) == 3
    @test best_rotation([1, 3, 0, 2, 4]) == 0
end


text:
function is_function(ex)
    if isexpr(ex, :(=)) || isexpr(ex, :->)
        isexpr(ex.args[1], :call) && return true
        isexpr(ex.args[1], :where) && return true
    elseif isexpr(ex, :function)
        return true
    end
    return false
end


text:
using Distributed
using StatsBase: mad, iqr, zscore
using Printf
import Glob
import Dierckx
import SparseArrays: spdiagm

include("./blocks.jl")

"""Helper to make a group path similar to `dset`, with new base"""
_match_dset_path(dset, newgroup) = join([newgroup; split(dset, '/')[2:end]], '/')

# E.g. `counts/1` when passed `stack/1`
_count_dset(dset) = _match_dset_path(dset, "counts")
_excl_dset(dset) = _match_dset_path(dset, "excluded")
# Used if we want to track the standard deviaition
# _stddev_dset(dset) = _match_dset_path(dset, "stddev")
# _stddev_raw_dset(dset) = _match_dset_path(dset, "stddev_raw")

# I'll be passing the included, clean geos, and want to record the bad missing ones
# is_excluded(total, included) = isnothing.(indexin(total, included))  # Too much space :(

function proc_pixel_linear(
    unw_stack_file,
    in_dset,
    valid_igram_indices,
    unmasked_idx_set,
    outfile,
    outdset,
    geolist,
    intlist,
    alpha,
    L1 = false;
    prune_outliers = true,
    sigma = 4,
    row_slice = nothing,
    col_slice = nothing,
)
    # in_buf = zeros(Float32, (length(col_slice), length(row_slice), length(valid_igram_indices)))
    in_buf = zeros(Float32, (length(row_slice), length(col_slice), length(valid_igram_indices)))
    try
        # in_buf .= h5read(unw_stack_file, in_dset, (col_slice, row_slice, :))[:, :, valid_igram_indices]
        in_buf .= h5read(unw_stack_file, in_dset, (row_slice, col_slice, :))[:, :, valid_igram_indices]
    catch e
        println(e, " ", "fail ", row_slice, " ", col_slice)
        return
    end
    # out_buf = Array{Float32, 2}(undef, (length(row_slice), length(col_slice)))
    # out_buf = Array{Float32, 2}(undef, (length(col_slice), length(row_slice)))
    # out_buf = zeros(Float32, (length(col_slice), length(row_slice)))
    out_buf = zeros(Float32, (length(row_slice), length(col_slice)))
    # println(row_slice, col_slice, length(row_slice), length(col_slice))

    for (i,row) in enumerate(row_slice), (j,col) in enumerate(col_slice)
    # seems to be transposed
    # for (i,row) in enumerate(col_slice), (j,col) in enumerate(row_slice)
        !((row, col) in unmasked_idx_set) && continue

        unw_pixel_raw = in_buf[i, j, :]
        # unw_pixel_raw = h5read(unw_stack_file, in_dset, (row, col, :))[1, 1, valid_igram_indices]
        if any(isnan.(unw_pixel_raw))
            println("$row, $col: nan")
            return
        end

        # Also load correlations for cutoff
        # cor_pixel = h5read(CC_FILENAME, "stack", (row, col, :))[1, 1, valid_igram_indices]
        # cor_thresh = 0.05
        # cor_pixel, cor_thresh = nothing, 0.0

        linear = true
        soln_phase, igram_count, geo_clean, unw_clean = calc_soln(
            unw_pixel_raw,
            geolist,
            intlist,
            alpha,
            linear;
            L1 = L1,
            prune_outliers = prune_outliers,
            sigma = sigma,
        )

        # println(size(soln_phase), size(out_buf))
        # println(soln_phase, " $i $j")
        out_buf[i, j] = soln_phase[1] * P2MM
    end
    # println(row_slice, " ", col_slice)
    # println(out_buf)

    dist_outfile = string(Distributed.myid()) * outfile
    h5open(dist_outfile, "r+") do f
        f[outdset][row_slice, col_slice] = out_buf
        # f[outdset][col_slice, row_slice] = out_buf
        # f[outdset][row, col] = P2MM * soln_phase
        # f[_count_dset(outdset)][row, col] = igram_count
        # f[_excl_dset(outdset)][row, col] = encode_missing(geolist, geo_clean)
        # f[_stddev_dset(outdset)][row, col] = std(unw_clean) .* abs(PHASE_TO_CM)
        # f[_stddev_raw_dset(outdset)][row, col] = std(unw_pixel_raw) .* abs(PHASE_TO_CM)
    end
    return
end

function proc_pixel_daily(
    unw_stack_file,
    in_dset,
    valid_igram_indices,
    unmasked_idx_set,
    outfile,
    outdset,
    geolist,
    intlist,
    alpha,
    L1 = false;
    prune_outliers = true,
    sigma = 4,
    row_slice = nothing,
    col_slice = nothing,
)
    in_buf = zeros(Float32, (length(row_slice), length(col_slice), length(valid_igram_indices)))
    try
        # in_buf .= h5read(unw_stack_file, in_dset, (col_slice, row_slice, :))[:, :, valid_igram_indices]
        in_buf .= h5read(unw_stack_file, in_dset, (row_slice, col_slice, :))[:, :, valid_igram_indices]
    catch e
        println(e, " ", "fail ", row_slice, " ", col_slice)
        return
    end
    out_buf = zeros(Float32, (length(row_slice), length(col_slice), length(geolist)))

    for (i,row) in enumerate(row_slice), (j,col) in enumerate(col_slice)
        !((row, col) in unmasked_idx_set) && continue

        unw_pixel_raw = in_buf[i, j, :]
        if any(isnan.(unw_pixel_raw))
            println("$row, $col: nan")
            return
        end
        linear = false
        soln_velos, igram_count, geo_clean, unw_clean = calc_soln(
            unw_pixel_raw,
            geolist,
            intlist,
            alpha,
            linear;
            L1 = L1,
            prune_outliers = prune_outliers,
            sigma = sigma,
        )
        timeseries_cm = _unreg_to_cm(soln_velos, geo_clean, geolist)
        out_buf[i, j, :] = timeseries_cm
    end

    dist_outfile = string(Distributed.myid()) * outfile
    h5open(dist_outfile, "r+") do f
        f[outdset][row_slice, col_slice, :] = out_buf
        # f[outdset][row, col, :] = phi_arr
        # f[_count_dset(outdset)][row, col] = igram_count
        # f[_excl_dset(outdset)][row, col] = encode_missing(geolist, geo_clean)
        # f[_stddev_dset(outdset)][row, col] = std(unw_clean) .* abs(PHASE_TO_CM)
        # f[_stddev_raw_dset(outdset)][row, col] = std(unw_pixel_raw) .* abs(PHASE_TO_CM)
    end
    return
end

function _unreg_to_cm(soln_velos, geolist_short, geolist_full)
    # First find the phase for only the clean dates we solved for
    timediffs = day_diffs(geolist_short)
    phi_clean = integrate_velocities(soln_velos, timediffs)

    # Then linearly interpolate between these for the removed phases
    return PHASE_TO_CM * interpolate_phase(geolist_short, phi_clean, geolist_full)
end

function calc_soln(
    unw_pixel,
    geolist,
    intlist,
    alpha,
    constant_velocity;
    L1 = true,
    cor_pixel = nothing,
    cor_thresh = 0.0,
    prune_outliers = true,
    sigma = 4,
)::Tuple{Array{Float32,1},Int64,Array,Array}
    geo_clean, intlist_clean, unw_clean = geolist, intlist, unw_pixel
    if prune_outliers
        geo_clean, intlist_clean, unw_clean =
            remove_outliers(geo_clean, intlist_clean, unw_clean, mean_sigma_cutoff = sigma)
    end
    if cor_thresh > 0 && !isnothing(cor_pixel)
        intlist_clean, unw_clean = prune_cor(
            intlist_clean,
            unw_clean,
            cor_pixel = cor_pixel,
            cor_thresh = cor_thresh,
        )
    end

    igram_count = length(unw_clean)

    B = prepB(geo_clean, intlist_clean, constant_velocity, alpha)
    # Last, pad with zeros if doing Tikh. regularization
    unw_final = alpha > 0 ? augment_zeros(B, unw_clean) : unw_clean

    if L1
        soln = invert_pixel_l1(unw_final, B)
    elseif constant_velocity
        # TODO: Do i ever want the backslash for this, which weights dt like 1/dt^2 ?
        soln = [sum(unw_final) / sum(B)]
    else
        soln = B \ unw_final
    end
    soln_phase = Float32.(soln)
    # end
    return soln_phase, igram_count, geo_clean, unw_clean
end

"""Run sbas on multiple processes"""
function run_sbas(
    unw_stack_file::String,
    dset::String,
    outfile::String,
    outdset::String,
    geolist,
    intlist,
    valid_igram_indices,
    constant_velocity::Bool,
    alpha::Real,
    L1::Bool = false,
    prune_outliers = true,
    sigma = 4,
)

    L1 ? println("Using L1 penalty for fitting") :
    println("Using least squares for fitting")
    prune_outliers ? println("Pruning .geo dates for outliers") :
    println("Not pruning outliers.")
    alpha > 0 ? println("Regularizing solution with alpha = $alpha") :
    println("No regularization")

    # TODO: do I ever really care about these to change as variable?
    # rho, alpha, abstol = 1.0, 1.6, 1e-3
    nrows, ncols, _ = size(unw_stack_file, dset)

    # Choose the correct processing function and size based on if we choose constant
    # TODO: is using multiple dispatch here better in any way (clarity/speed)?
    if constant_velocity
        println("Using constant velocity for inversion solution")
        proc_func = proc_pixel_linear
        outsize = (nrows, ncols)
        chunksize = (500, 500)
    else
        println("Not using linear: Finding solution for each day")
        proc_func = proc_pixel_daily
        outsize = (nrows, ncols, length(geolist))
        chunksize = (10, 10, length(geolist))
    end

    println("Making new writing file for each of $(length(Distributed.workers()))")
    @time @sync @distributed for id in Distributed.workers()
        h5open(string(id) * outfile, "w") do fout
            d_create(fout, outdset, datatype(Float32), dataspace(outsize)) #, "chunk", chunksize)
            d_create(fout, _count_dset(outdset), datatype(Int32), dataspace(nrows, ncols))
            d_create(fout, _excl_dset(outdset), datatype(Int64), dataspace(nrows, ncols))
            # d_create(fout, _stddev_dset(outdset), datatype(Float64), dataspace(nrows, ncols))
            # d_create(fout, _stddev_raw_dset(outdset), datatype(Float64), dataspace(nrows, ncols))
        end
    end

    # Get blocks to load and process
    unmasked_idx_set = Set(get_unmasked_idxs(geolist))  # Set of tuples
    # windows returns (row_slice, col_slice)
    w = collect(windows(unw_stack_file, dset))
    # print(w[1:10])

    @time @sync @distributed for (row_slice, col_slice) in w
    # @time @sync @distributed for (row_slice, col_slice) in w[1:3]
        # @time @sync @distributed for (row, col) in collect(Iterators.product(300:400, 300:400))
        proc_func(
            unw_stack_file,
            dset,
            valid_igram_indices,
            unmasked_idx_set,
            outfile,
            outdset,
            geolist,
            intlist,
            alpha,
            L1;
            prune_outliers = prune_outliers,
            sigma = sigma,
            row_slice = row_slice,
            col_slice = col_slice,
        )
    end
    println("Merging files into $outfile")
    @time merge_partial_files(outfile, outdset, _count_dset(outdset), _excl_dset(outdset))
    # @time merge_partial_files(outfile, outdset)
    # _stddev_dset(outdset), _stddev_raw_dset(outdset))

    # _save_std_attrs(outfile, outdset)
    return outfile, outdset
end

# Need the .I so we can use to load from h5
get_unmasked_idxs(geolist, do_permute = false) = [
    cart_idx.I for cart_idx in findall(.!Sario.load_mask(geolist, do_permute = do_permute))
]

function _save_std_attrs(outfile, outdset)
    attr_dict1 = Dict(
        "description" => "original unwrapped (shifted) std. dev of pixels",
        "units" => "cm",
    )
    h5writeattr(outfile, _stddev_raw_dset(outdset), attr_dict1)
    attr_dict2 =
        Dict("description" => "Final, outlier-removed std. dev of pixels", "units" => "cm")
    h5writeattr(outfile, _stddev_dset(outdset), attr_dict2)
end

"""Linearly interpolate the values between dates that we didn't solve for"""
function interpolate_phase(geo_clean, phis, geolist_full)
    return interp1d(_get_day_nums(geo_clean), phis, _get_day_nums(geolist_full))
end

function interp1d(x, y, xv, k = 1)
    sp = Dierckx.Spline1D(x, y, k = k)
    return sp(xv)
end

function prepB(geolist, intlist, constant_velocity = false, alpha = 0)
    if constant_velocity
        return Float32.(reshape(temporal_baseline(intlist), :, 1))
    end
    # Prepare A and B matrix used for each pixel inversion
    # A = build_A_matrix(geolist, intlist)
    B = build_B_matrix(geolist, intlist)

    if alpha > 0
        B = augment_B(B, alpha)
    end
    return B
end

# In case we make each pixel inversion more complicated (extra penalty functions, etc.)
# TODO: does a Val type for the case with no extra zeros help here?
function invert_pixel(pixel::Array{Float32,1}, pB::Array{Float32,2}, extra_zeros = 0)
    if extra_zeros == 0
        return pB * pixel
    else
        return pB * vcat(pixel, zeros(extra_zeros))
    end
end

# Defined in optimize.jl
function invert_pixel_l1(
    pixel::AbstractArray{<:AbstractFloat},
    B::AbstractArray{<:AbstractFloat};
    rho = 1.0,
    alpha = 1.6,
    lu_tuple = nothing,
    abstol = 1e-3,
)
    return Float32.(huber_fit(
        Float64.(B),
        Float64.(pixel),
        rho,
        alpha,
        lu_tuple = lu_tuple,
        abstol = abstol,
    ))
end

"""Takes the list of igram dates and builds the SBAS B (velocity coeff) matrix

Each row corresponds to an igram, each column to a .geo
Values will be t_k+1 - t_k for columns after the -1 in A,
up to and including the +1 entry
"""
function build_B_matrix(geolist::Array{Date,1}, intlist::Array{Igram,1})

    timediffs = day_diffs(geolist)
    # We take the first .geo to be time 0, leave out of matrix
    # Match on date (not time) to find indices
    geolist2 = @view geolist[2:end]

    M = length(intlist)  # Number of igrams,
    N = length(geolist2)

    B = zeros(Float32, (M, N))

    for i = 1:M
        # row = A[i, :]
        early_date, late_date = intlist[i]
        idx_early = findfirst(isequal(early_date), geolist2)
        if isnothing(idx_early)  # The first SLC will not be in the matrix
            start_idx = 1
        else
            start_idx = idx_early + 1
        end
        idx_late = findfirst(isequal(late_date), geolist2)
        end_idx = idx_late

        # Now fill in the time diffs in the index range
        B[i, start_idx:end_idx] .= @view timediffs[start_idx:end_idx]
    end

    return B
end

"""
    build_A_matrix(geolist::Array{Date, 1}, intlist::Array{Igram, 1}) -> Array{Int8, 2}

Takes the list of igram dates and builds the SBAS A matrix

Returns the incident-like matrix from the SBAS paper: A*phi = dphi
    Each row corresponds to an igram, each column to a .geo
    value will be -1 on the early (slave) igrams, +1 on later (master)
"""
function build_A_matrix(geolist::Array{Date,1}, intlist::Array{Igram,1})
    # We take the first .geo to be time 0, leave out of matrix
    # Match on date (not time) to find indices
    geolist = geolist[2:end]

    M = length(intlist)  # Number of igrams,
    N = length(geolist)
    A = zeros(Int8, M, N)
    for i = 1:M
        early_date, late_date = intlist[i]
        idx_early = findfirst(isequal(early_date), geolist)
        if ~isnothing(idx_early)  # The first SLC will not be in the matrix
            A[i, idx_early] = -1
        end
        idx_late = findfirst(isequal(late_date), geolist)
        A[i, idx_late] = 1

    end
    return A
end

_D(n) = spdiagm(0 => -ones(n - 1), 1 => ones(n - 1))[1:end-1, :]

"""For Tikhonov regularization, pad the B matrix with alpha*D"""
function augment_B(B::Array{Float32,2}, alpha::Real)
    extra = alpha * I
    n = size(B, 2)
    return Float64.(vcat(B, alpha * _D(n)))  # Note: as of 9/29/19, qr has bug for sparse float32 matrices
end

augment_zeros(B, unw) = vcat(unw, zeros(size(B, 1) - length(unw)))

# function augment_B(B::Array{Float32, 2}, alpha::Float32)
#     extra = alpha*I
#     return Float32.(vcat(B, ))
# end

# TODO: this can't work for huge stacks
function augment_matrices(B::Array{Float32,2}, unw_stack::Array{Float32,3}, alpha::Real)
    B = Float32.(vcat(B, alpha * I))
    # Now make num rows match
    nrows, ncols, nlayers = size(unw_stack)
    zeros_shape = (nrows, ncols, size(B, 1) - nlayers)
    unw_stack = Float32.(cat(unw_stack, zeros(zeros_shape), dims = 3))
    return B, unw_stack
end

partial_files(outfile) = Glob.glob("[0-9]*" * outfile)

function merge_partial_files(outfile, dsets...)
    f1 = partial_files(outfile)[1]
    for dset in dsets
        tmp = zeros(eltype(f1, dset), size(f1, dset))

        for f in partial_files(outfile)
            cur = h5read(f, dset)
            tmp += cur
        end
        h5open(outfile, "cw") do fout
            fout[dset] = tmp
        end
    end
    for f in partial_files(outfile)
        rm(f)
    end
end

"""Finds the number of days between successive .geo files"""
function day_diffs(geolist::Array{Date,1})
    [difference.value for difference in diff(geolist)]
end

"""Takes velocity solution output and finds phases

Arguments:
    velocities come from invert_sbas
    timediffs are the days between each SAR acquisitions
        length will be 1 less than num SAR acquisitions
"""
function integrate_velocities(
    velocities::AbstractArray{<:AbstractFloat,1},
    timediffs::Array{Int,1},
)
    phi_diffs = velocities .* timediffs
    phi_arr = cumsum(coalesce.(phi_diffs, 0))
    pushfirst!(phi_arr, 0)
    return phi_arr
end

### Outlier functions ######
# 3 Reasons for removing igrams:
# 1. remove all from .geo dates with huge averages (whole day is garbage)
# 2. remove very low correlation
# 3. remove longer igrams w/ longer baseline than is possible
#   e.g. if we have 2.5 cm/year velocity, anything longer than ~1 year
#   will be all noise- we can't reliably sense such quickly moving ground
#
# TODO: figure out which should go in separate file for cleanliness
function remove_outliers(geolist, intlist, unw_pixel; mean_sigma_cutoff = 4)
    # TODO: verify this third criteria?

    # @show mean_sigma_cutoff
    # @show size(intlist), size(geolist), std(unw_pixel)
    # 1. find outliers in this pixels' values and remove them
    geo_clean, intlist_clean, unw_clean =
        peel_nsigma(geolist, intlist, unw_pixel, nsigma = mean_sigma_cutoff)
    # @show size(intlist_clean), size(geo_clean), std(unw_clean)
    return geo_clean, intlist_clean, unw_clean
end

# 2. low correlation cleaning
function prune_cor(intlist, unw_pixel, cor_pixel, cor_thresh = 0.0)
    low_cor_igrams = intlist[cor_pixel.<cor_thresh]
    intlist_clean, unw_clean = remove_igrams(intlist, unw_pixel, low_cor_igrams)
    return intlist_clean, unw_clean
end

function remove_igrams(
    intlist,
    unw_vals,
    bad_items::Union{Date,AbstractArray{Date},AbstractArray{Igram}}...,
)
    good_idxs = trues(size(intlist))
    for item in bad_items
        good_idxs .&= _good_idxs(item, intlist)
    end
    return remove_by_idx(good_idxs, intlist, unw_vals)
end

"""Pass directly the indexes of good ones igrams (removes all BUT the good_idxs)"""
function remove_by_idx(good_idxs::AbstractArray, intlist, unw_vals)
    unw_clean = unw_vals[good_idxs]
    intlist_clean = intlist[good_idxs]
    return intlist_clean, unw_clean
end

# If we pass the B, also remove rows for bad idxs
function remove_by_idx(good_idxs::AbstractArray, intlist, unw_vals, B)
    B_clean = B[good_idxs, :]
    unw_clean = unw_vals[good_idxs]
    intlist_clean = intlist[good_idxs]
    return intlist_clean, unw_clean, B_clean
end

function _good_idxs(bad_date_arr::AbstractArray{Date}, intlist)
    return .!reduce(
        (x, y) -> x .| y,
        (b in intlist for b in bad_date_arr),
        init = falses(size(intlist)),
    )
end

function _good_idxs(bad_igram_arr::AbstractArray{Igram}, intlist)
    return [!(ig in bad_igram_arr) for ig in intlist]
end

_good_idxs(bad_date::Date, intlist) = .!(bad_date in intlist)

""" 4* the sigma valud as calculated using MAD
Used for robust variance est. (as a cutoff for outliers)
See https://en.wikipedia.org/wiki/Robust_measures_of_scale#IQR_and_MAD
Normalize makes it equal to stddev for normally dist. data"""
mednsigma(arr, n = 4) = n * mad(arr, normalize = true)

modzscore(arr) = (arr .- median(arr)) ./ mad(arr, normalize = true)

# min_spread added in case `mad` produces a very low number,
# so if you want to avoid removing lots, set based on your data
function two_way_cutoff(arr, nsigma, min_spread = 0)
    # min_spread = abs(2 * median(arr))
    # spread = 1.5 * iqr(arr)
    spread = max(min_spread, mednsigma(arr, nsigma))
    #@show spread
    # spread = max(min_spread, std(arr)*nsigma)
    #@show std(arr)*nsigma

    # Make we dont cut off very low var points
    low = min(0, median(arr) - spread)
    high = median(arr) + spread
    return (low, high)
end

function two_way_outliers(arr, nsigma, min_spread = 0)
    # return abs.(modzscore(arr)) .> nsigma
    low, high = two_way_cutoff(arr, nsigma, min_spread)
    return (arr .< low) .| (arr .> high)
end

"""Simpler means by day: just abs. value of phases, either mean or median"""
mean_abs_val(geolist, intlist, unw_vals) =
    [mean(abs.(vals_by_date(d, intlist, unw_vals))) for d in geolist];
median_abs_val(geolist, intlist, unw_vals) =
    [median(abs.(vals_by_date(d, intlist, unw_vals))) for d in geolist];
# valsq(geolist, intlist, unw_vals) = [mean(vals_by_date(d, intlist, unw_vals) .^ 2) for d in geolist];

"""Get the values (unw, cc, etc.) of one date"""
vals_by_date(date::Date, intlist, vals) = vals[date in intlist]

"""Get the values (unw, cc, etc.) grouped by each date"""
vals_by_date(date_arr::Array{Date}, intlist, vals) = [vals[d in intlist] for d in date_arr]

"""Remove all igrams corresponding to the dates with means falling outside an `nsigma` interval"""
function peel_nsigma(geo, int, val; nsigma = 4)
    dates_to_remove = nsigma_days(geo, int, val, nsigma)
    int2, val2 = remove_igrams(int, val, dates_to_remove)
    geo2 = [g for g in geo if !(g in dates_to_remove)]
    return geo2, int2, val2
end

"""Return the days of `geo` which are more than nsigma away from mean"""
function nsigma_days(geo, int, val, nsigma = 4, min_spread = 2)
    means = mean_abs_val(geo, int, val)
    out_idxs = two_way_outliers(means, nsigma, min_spread)

    # FOR PRINTING ONLY
    # low, high = two_way_cutoff(means, nsigma, min_spread)
    # println("Using cutoff around $(median(means)), spread $(max(min_spread, mednsigma(means, nsigma))) : ($low, $high) ")
    # println("Number removed: $(sum(out_idxs))")
    # @show std(means), mad(means, normalize=true)
    return geo[out_idxs]
end

# TODO: this is gonna fail for more than 63 dates ( log2(typemax(Int64)) = 63)...
function encode_missing(total, included)
    out = Int64(0)
    for (i, t) in enumerate(total)
        if !(t in included)
            out += 2^(i - 1)
        end
    end
    return out
end

function decode_excluded(excl_int::Int64, geos::Array{Date,1})
    excl_bitstring = reverse(bitstring(excl_int)[end-length(geos)+1:end])
    return [geos[i] for (i, g) in enumerate(geos) if excl_bitstring[i] == '1']
end

function decode_excluded(excl_int_map::AbstractArray{Int64,2}, geos::Array{Date,1})
    out = Array{Array{Date,1},2}(undef, size(excl_int_map))
    for idx in eachindex(excl_int_map)
        out[idx] = decode_excluded(excl_int_map[idx], geos)
    end
    return out
end


text:
export HierarchicalMatrix, isfactored, condest, blockrank

condest(A)=1
blocksize(A) = (1,1)
blockrank(A)=rank(A)

##
# Represent a binary hierarchical matrix
# [ x x \ /
#   x x / \
#   \ / x x
#   / \ x x ]
# where the diagonal blocks are of one type, and
# the off-diagonal blocks are of a type with low-rank structure.
#
# To have simple recursive construction, we use the ordering:
# [ D₁   L₄
#      ⋯      L₂
#   L₃   D₂
#          D₃   L₆
#      L₁     ⋯
#          L₅   D₄ ]
# The diagonal blocks are ordered along the diagonal,
# and the off-diagonal blocks are ordered from: bottom left,
# top right, then followed recursively by top left and bottom right.
##

typealias AbstractHierarchicalMatrix{S,V,T,HS,HV} AbstractHierarchicalArray{Tuple{S,V},T,Tuple{HS,HV},2}

type HierarchicalMatrix{S,V,T,HS,HV} <: AbstractHierarchicalMatrix{S,V,T,HS,HV}
    diagonaldata::HS
    offdiagonaldata::HV

    A::Matrix{T} # Cache of matrix for pivot computation
    factorization::PivotLDU{T,Matrix{T}} # Cache of factorization of A for pivot computation
    factored::Bool

    function HierarchicalMatrix(diagonaldata::HS,offdiagonaldata::HV)
        H = new()
        H.diagonaldata = diagonaldata
        H.offdiagonaldata = offdiagonaldata

        r = mapreduce(rank,+,offdiagonaldata)
        H.A = eye(T,r,r)
        H.factored = false

        H
    end
end

HierarchicalMatrix{S1,S2,V1,V2}(diagonaldata::Tuple{S1,S2},offdiagonaldata::Tuple{V1,V2}) = HierarchicalMatrix{promote_type(S1,S2),promote_type(V1,V2),promote_type(eltype(S1),eltype(S2),eltype(V1),eltype(V2)),Tuple{S1,S2},Tuple{V1,V2}}(diagonaldata,offdiagonaldata)
HierarchicalMatrix{S1,S2,V1,V2,V3,V4,T1,T2,HS1,HS2,HV1,HV2}(diagonaldata::Tuple{HierarchicalMatrix{S1,V1,T1,HS1,HV1},HierarchicalMatrix{S2,V2,T2,HS2,HV2}},offdiagonaldata::Tuple{V3,V4}) = HierarchicalMatrix{promote_type(S1,S2),promote_type(V1,V2,V3,V4),promote_type(eltype(S1),eltype(S2),eltype(V1),eltype(V2),eltype(V3),eltype(V4),T1,T2),Tuple{HierarchicalMatrix{S1,V1,T1,HS1,HV1},HierarchicalMatrix{S2,V2,T2,HS2,HV2}},Tuple{V3,V4}}(diagonaldata,offdiagonaldata)
HierarchicalMatrix{S,S1,V,V1,V2,T,HS,HV}(diagonaldata::Tuple{S1,HierarchicalMatrix{S,V,T,HS,HV}},offdiagonaldata::Tuple{V1,V2}) = HierarchicalMatrix{promote_type(S,S1),promote_type(V,V1,V2),promote_type(eltype(S),eltype(S1),eltype(V),eltype(V1),eltype(V2),T),Tuple{S1,HierarchicalMatrix{S,V,T,HS,HV}},Tuple{V1,V2}}(diagonaldata,offdiagonaldata)
HierarchicalMatrix{S,S1,V,V1,V2,T,HS,HV}(diagonaldata::Tuple{HierarchicalMatrix{S,V,T,HS,HV},S1},offdiagonaldata::Tuple{V1,V2}) = HierarchicalMatrix{promote_type(S,S1),promote_type(V,V1,V2),promote_type(eltype(S),eltype(S1),eltype(V),eltype(V1),eltype(V2),T),Tuple{HierarchicalMatrix{S,V,T,HS,HV},S1},Tuple{V1,V2}}(diagonaldata,offdiagonaldata)

HierarchicalMatrix(diagonaldata::Vector,offdiagonaldata::Vector)=HierarchicalMatrix(diagonaldata,offdiagonaldata,round(Int,log2(length(diagonaldata))))

function HierarchicalMatrix(diagonaldata::Vector,offdiagonaldata::Vector,n::Int)
    @assert length(diagonaldata) == 2^n
    @assert length(offdiagonaldata) == 2^(n+1)-2

    if n == 1
        return HierarchicalMatrix(tuple(diagonaldata...),tuple(offdiagonaldata...))
    elseif n ≥ 2
        dldp2 = div(length(offdiagonaldata)+2,2)
        return HierarchicalMatrix((HierarchicalMatrix(diagonaldata[1:2^(n-1)],offdiagonaldata[3:dldp2],n-1),
                                   HierarchicalMatrix(diagonaldata[1+2^(n-1):end],offdiagonaldata[dldp2+1:end],n-1)),
                                   (offdiagonaldata[1],offdiagonaldata[2]))
    end
end

Base.similar(H::HierarchicalMatrix) = HierarchicalMatrix(map(similar,diagonaldata(H)),map(similar,offdiagonaldata(H)))
Base.similar{SS,V,T}(H::HierarchicalMatrix{SS,V,T}, S) = HierarchicalMatrix(map(A->similar(A,S),diagonaldata(H)),map(A->similar(A,S),offdiagonaldata(H)))

diagonaldata(H::HierarchicalMatrix) = H.diagonaldata
offdiagonaldata(H::HierarchicalMatrix) = H.offdiagonaldata

nlevels(H::HierarchicalMatrix) = 1+mapreduce(nlevels,max,diagonaldata(H))

partition(H::HierarchicalMatrix) = diagonaldata(H),offdiagonaldata(H)

collectoffdiagonaldata{S,V,T}(H::HierarchicalMatrix{S,V,T,NTuple{2,S}}) = collect(offdiagonaldata(H))
function collectoffdiagonaldata{S,V,T,HS}(H::HierarchicalMatrix{S,V,T,HS})
    data = collect(offdiagonaldata(H))
    append!(data,mapreduce(collectoffdiagonaldata,vcat,diagonaldata(H)))
    data
end

collectdiagonaldata{S,V,T}(H::HierarchicalMatrix{S,V,T,NTuple{2,S}}) = collect(diagonaldata(H))
function collectdiagonaldata{S,V,T,HS}(H::HierarchicalMatrix{S,V,T,HS})
    data = S[]
    append!(data,mapreduce(collectdiagonaldata,vcat,diagonaldata(H)))
    data
end

isfactored(H::HierarchicalMatrix) = H.factored

function condest(H::HierarchicalMatrix)
    !isfactored(H) && factorize!(H)
    return cond(H.A)*mapreduce(condest,+,diagonaldata(H))
end

Base.convert{S,V,T,HS,HV}(::Type{HierarchicalMatrix{S,V,T,HS,HV}},M::HierarchicalMatrix) = HierarchicalMatrix(convert(Vector{S},collectdiagonaldata(M)),convert(Vector{V},collectoffdiagonaldata(M)))
Base.convert{T}(::Type{Matrix{T}},M::HierarchicalMatrix) = full(M)
Base.promote_rule{S,V,T,HS,HV,SS,VV,TT,HSS,HVV}(::Type{HierarchicalMatrix{S,V,T,HS,HV}},::Type{HierarchicalMatrix{SS,VV,TT,HSS,HVV}})=HierarchicalMatrix{promote_type(S,SS),promote_type(V,VV),promote_type(T,TT),promote_type(HS,HSS),promote_type(HV,HVV)}
Base.promote_rule{T,SS,VV,TT,HSS,HVV}(::Type{Matrix{T}},::Type{HierarchicalMatrix{SS,VV,TT,HSS,HVV}})=Matrix{promote_type(T,TT)}
Base.promote_rule{T,SS,VV,TT,HSS,HVV}(::Type{LowRankMatrix{T}},::Type{HierarchicalMatrix{SS,VV,TT,HSS,HVV}})=Matrix{promote_type(T,TT)}

Base.transpose(H::HierarchicalMatrix) = HierarchicalMatrix(map(transpose,diagonaldata(H)),map(transpose,reverse(offdiagonaldata(H))))
Base.ctranspose(H::HierarchicalMatrix) = HierarchicalMatrix(map(ctranspose,diagonaldata(H)),map(ctranspose,reverse(offdiagonaldata(H))))

function Base.size{S<:AbstractMatrix,V<:AbstractMatrix}(H::HierarchicalMatrix{S,V})
    H11,H22 = diagonaldata(H)
    H21,H12 = offdiagonaldata(H)
    m1,n2 = size(H12)
    m2,n1 = size(H21)
    @assert (m1,n1) == size(H11)
    @assert (m2,n2) == size(H22)
    m1+m2,n1+n2
end

function Base.getindex{S<:AbstractMatrix,V<:AbstractMatrix}(H::HierarchicalMatrix{S,V},i::Int,j::Int)
    H11,H22 = diagonaldata(H)
    H21,H12 = offdiagonaldata(H)
    m1,n2 = size(H12)
    m2,n1 = size(H21)
    if 1 ≤ i ≤ m1
        if 1 ≤ j ≤ n1
            return getindex(H11,i,j)
        elseif n1 < j ≤ n1+n2
            return getindex(H12,i,j-n1)
        else
            throw(BoundsError())
        end
    elseif m1 < i ≤ m1+m2
        if 1 ≤ j ≤ n1
            return getindex(H21,i-m1,j)
        elseif n1 < j ≤ n1+n2
            return getindex(H22,i-m1,j-n1)
        else
            throw(BoundsError())
        end
    else
        throw(BoundsError())
    end
end
Base.getindex{S<:AbstractMatrix,V<:AbstractMatrix}(H::HierarchicalMatrix{S,V},i::Int,jr::Range) = eltype(H)[H[i,j] for j=jr].'
Base.getindex{S<:AbstractMatrix,V<:AbstractMatrix}(H::HierarchicalMatrix{S,V},ir::Range,j::Int) = eltype(H)[H[i,j] for i=ir]
Base.getindex{S<:AbstractMatrix,V<:AbstractMatrix}(H::HierarchicalMatrix{S,V},ir::Range,jr::Range) = eltype(H)[H[i,j] for i=ir,j=jr]
Base.full{S<:AbstractMatrix,V<:AbstractMatrix}(H::HierarchicalMatrix{S,V})=H[1:size(H,1),1:size(H,2)]

Base.copy(H::HierarchicalMatrix) = HierarchicalMatrix(map(copy,diagonaldata(H)),map(copy,offdiagonaldata(H)))
Base.copy!(H::HierarchicalMatrix,J::HierarchicalMatrix) = (map(copy!,diagonaldata(H),diagonaldata(J));map(copy!,offdiagonaldata(H),offdiagonaldata(J));H)

Base.rank(H::HierarchicalMatrix) = rank(full(H))
Base.cond(H::HierarchicalMatrix) = cond(full(H))

function blockrank(H::HierarchicalMatrix)
    m,n = blocksize(H)
    A = Array{Int64}(m,n)
    (m1,n1),(m2,n2) = map(blocksize,diagonaldata(H))
    r1,r2 = map(rank,offdiagonaldata(H))
    for j=1:n1,i=1:m2
        A[i+m1,j] = r1
    end
    for j=1:n2,i=1:m1
        A[i,j+n1] = r2
    end
    A11,A22 = map(blockrank,diagonaldata(H))
    for j=1:n1,i=1:m1
        A[i,j] = A11[i,j]
    end
    for j=1:n2,i=1:m2
        A[i+m1,j+n1] = A22[i,j]
    end
    A
end

for op in (:+,:-,:.+,:.-)
    @eval begin
        $op(H::HierarchicalMatrix) = HierarchicalMatrix(map($op,diagonaldata(H)),map($op,offdiagonaldata(H)))
        $op(H::HierarchicalMatrix,J::HierarchicalMatrix) = HierarchicalMatrix(map($op,diagonaldata(H),diagonaldata(J)),map($op,offdiagonaldata(H),offdiagonaldata(J)))
        $op(H::HierarchicalMatrix,L::LowRankMatrix) = $op(promote(H,L)...)
        $op(L::LowRankMatrix,H::HierarchicalMatrix) = $op(promote(L,H)...)
        $op(H::HierarchicalMatrix,A::Matrix) = $op(promote(H,A)...)
        $op(A::Matrix,H::HierarchicalMatrix) = $op(promote(A,H)...)
    end
end

function *(H::HierarchicalMatrix, x::HierarchicalVector)
    T = promote_type(eltype(H),eltype(x))
    A_mul_B!(similar(x,T),H,x)
end

function Base.A_mul_B!(b::AbstractVector,D::Base.LinAlg.Diagonal,x::AbstractVector)
    d = D.diag
    for i=1:min(length(b),length(x))
        @inbounds b[i] += d[i]*x[i]
    end
    b
end

function Base.A_mul_B!(b::HierarchicalVector,H::HierarchicalMatrix,h::HierarchicalVector)
    (H11,H22),(H21,H12) = partition(H)
    h1,h2 = partition(h)
    b1,b2 = partition(b)
    A_mul_B!(b1,H12,h2)
    A_mul_B!(b1,H11,h1)
    A_mul_B!(b2,H21,h1)
    A_mul_B!(b2,H22,h2)
    b
end

function *(H::HierarchicalMatrix,b::Vector)
    (H11,H22),(H21,H12) = partition(H)
    m1,m2 = size(H12,1),size(H21,1)
    (b1,b2) = (b[1:m1],b[1+m1:m1+m2])
    HierarchicalVector((H11*b1⊕H12*b2,H21*b1⊕H22*b2))
end

\(H::HierarchicalMatrix,b::AbstractVecOrMat) = full(H)\b


text:
function ZeroValue()
    value = Ref{API.pmix_value_t}()
    ccall(:memset, Ptr{Cvoid}, (Ptr{Cvoid}, Cint, Csize_t), value, 0, sizeof(API.pmix_value_t))
    return value[]
end

function UndefValue()
    value = ZeroValue()
    value = @set value.type = API.PMIX_UNDEF
    return value
end

function destruct(value)
    r_value = Ref{API.pmix_value_t}(value)
    @check API.pmix_value_destruct(r_value)
    return nothing
end

# function unload(value)
#     r_value = Ref{API.pmix_value_t}(value)
#     data = Ref{Ptr{Cvoid}}()
#     sz = Ref{Csize_t}()
#     @check API.pmix_value_unload(r_value, data, sz)
#     return (data[], sz[])
# end

function load(value, data, type)
    r_value = Ref{API.pmix_value_t}(value)
    API.pmix_value_load(r_value, data, type)
    return r_value[]
end

function Value(data, type)
    value = ZeroValue()
    value = load(value, Ref(data), type)
    return value
end

# function xfer()
# end

macro get_number(T, PMIX_T, value)
    quote
        let value = $(esc(value))
            if value.type != $API.$(PMIX_T)
                throw($PMIxException($API.PMIX_ERR_TYPE_MISMATCH))
            end
            data = Ref(value.data)
            GC.@preserve data begin
                ptr = Base.unsafe_convert(Ptr{Cvoid}, data)
                result = unsafe_load(reinterpret(Ptr{$T}, ptr))
            end
            result
        end
    end
end

import Base: convert
function convert(::Type{Int32}, value::API.pmix_value_t)
    @assert Cint == Int32
    @assert Cint == API.pid_t
    if !(value.type == API.PMIX_INT || value.type ==  API.PMIX_INT32 || value.type == API.PMIX_PID)
        throw(PMIxException(API.PMIX_ERR_TYPE_MISMATCH))
    end
    data = Ref(value.data)
    GC.@preserve data begin
        ptr = Base.unsafe_convert(Ptr{Cvoid}, data)
        return unsafe_load(reinterpret(Ptr{Int32}, ptr))
    end
end

function convert(::Type{UInt32}, value::API.pmix_value_t)
    @assert Cuint == UInt32
    if !(value.type == API.PMIX_UINT || value.type ==  API.PMIX_UINT32)
        throw(PMIxException(API.PMIX_ERR_TYPE_MISMATCH))
    end
    data = Ref(value.data)
    GC.@preserve data begin
        ptr = Base.unsafe_convert(Ptr{Cvoid}, data)
        return unsafe_load(reinterpret(Ptr{UInt32}, ptr))
    end
end

function convert(::Type{UInt64}, value::API.pmix_value_t)
    @assert Csize_t == UInt64
    if !(value.type == API.PMIX_SIZE || value.type ==  API.PMIX_UINT64)
        throw(PMIxException(API.PMIX_ERR_TYPE_MISMATCH))
    end
    data = Ref(value.data)
    GC.@preserve data begin
        ptr = Base.unsafe_convert(Ptr{Cvoid}, data)
        return unsafe_load(reinterpret(Ptr{UInt64}, ptr))
    end
end

convert(::Type{Bool}, value::API.pmix_value_t) = @get_number(Bool, PMIX_BOOL, value)
convert(::Type{Int8}, value::API.pmix_value_t) = @get_number(Int8, PMIX_INT8, value)
convert(::Type{Int16}, value::API.pmix_value_t) = @get_number(Int16, PMIX_INT16, value)
convert(::Type{Int64}, value::API.pmix_value_t) = @get_number(Int64, PMIX_INT64, value)
convert(::Type{UInt8}, value::API.pmix_value_t) = @get_number(UInt8, PMIX_UINT8, value)
convert(::Type{UInt16}, value::API.pmix_value_t) = @get_number(UInt16, PMIX_UINT16, value)
convert(::Type{Float32}, value::API.pmix_value_t) = @get_number(Float32, PMIX_FLOAT, value)
convert(::Type{Float64}, value::API.pmix_value_t) = @get_number(Float64, PMIX_DOUBLE, value)

function get_number(value::API.pmix_value_t)
    if value.type == API.PMIX_SIZE
        return @get_number(Csize_t, PMIX_SIZE, value)
    elseif value.type == API.PMIX_PID
        return @get_number(API.pid_t, PMIX_PID, value)
    elseif value.type == API.PMIX_INT
        return @get_number(Cint, PMIX_INT, value)
    elseif value.type == API.PMIX_INT8
        return @get_number(Int8, PMIX_INT8, value)
    elseif value.type == API.PMIX_INT16
        return @get_number(Int16, PMIX_INT16, value)
    elseif value.type == API.PMIX_INT32
        return @get_number(Int32, PMIX_INT32, value)
    elseif value.type == API.PMIX_INT64
        return @get_number(Int64, PMIX_INT64, value)
    elseif value.type == API.PMIX_UINT
        return @get_number(Cuint, PMIX_UINT, value)
    elseif value.type == API.PMIX_UINT8
        return @get_number(UInt8, PMIX_UINT8, value)
    elseif value.type == API.PMIX_UINT16
        return @get_number(UInt16, PMIX_UINT16, value)
    elseif value.type == API.PMIX_UINT32
        return @get_number(UInt32, PMIX_UINT32, value)
    elseif value.type == API.PMIX_UINT64
        return @get_number(UInt64, PMIX_UINT64, value)
    elseif value.type == API.PMIX_FLOAT
        return @get_number(Float32, PMIX_FLOAT, value)
    elseif value.type == API.PMIX_DOUBLE
        return @get_number(Float64, PMIX_DOUBLE, value)
    end
    error("Unkown numeric typed $(value.type)")
end

max_issues_repo_name: JuliaParallel/PMIx

max_forks_repo_name: JuliaParallel/PMIx


text:
# Interface -------

function watchers(owner, repo; auth = AnonymousAuth(), options...)
    watchers(auth, owner, repo; options...)
end

function watchers(auth::Authorization, owner, repo; headers = Dict(), options...)
    authenticate_headers(headers, auth)
    r = get(URI(API_ENDPOINT; path = "/repos/$owner/$repo/subscribers");
            headers = headers,
            options...)

    handle_error(r)

    data = JSON.parse(r.data)
    map!( u -> User(u), data)
end

function watched(user; auth = AnonymousAuth(), options...)
    watched(auth, user; options...)
end

function watched(auth::Authorization, user; headers = Dict(), options...)
    authenticate_headers(headers, auth)
    r = get(URI(API_ENDPOINT; path = "/users/$user/subscriptions"); headers = headers,
                                                                    options...)
    handle_error(r)

    data = JSON.parse(r.data)
    map!( r -> Repo(r), data)
end

function watching(owner, repo; auth = AnonymousAuth(), options...)
    watching(auth, owner, repo; options...)
end

function watching(auth::Authorization, owner, repo; headers = Dict(), options...)
    authenticate_headers(headers, auth)
    r = get(URI(API_ENDPOINT; path = "/repos/$owner/$repo/subscription");
            headers = headers,
            options...)

    r.status == 200 && return true
    r.status == 404 && return false

    handle_error(r)  # 404 is not an error in this case

    return false  # who knows... assume no
end

function watch(owner, repo; auth = AnonymousAuth(), options...)
    watch(auth, owner, repo; options...)
end

function watch(auth::Authorization, owner, repo; headers = Dict(),
                                                 query = Dict(),
                                                 subscribed = nothing,
                                                 ignored = nothing,
                                                 options...)
    authenticate_headers(headers, auth)

    subscribed != nothing && (query["subscribed"] = subscribed)
    ignored != nothing && (query["ignored"] = ignored)

    println("QUERY: ", query)

    r = put(URI(API_ENDPOINT; path = "/repos/$owner/$repo/subscription"); headers = headers,
                                                                          query = query,
                                                                          options...)
    handle_error(r)
end

function unwatch(owner, repo; auth = AnonymousAuth(), options...)
    unwatch(auth, owner, repo; options...)
end

function unwatch(auth::Authorization, owner, repo; headers = Dict(), options...)
    authenticate_headers(headers, auth)

    r = delete(URI(API_ENDPOINT; path = "/repos/$owner/$repo/subscription"); headers = headers,
                                                                             options...)
    handle_error(r)
end


text:
"""
    RProp_params
"""
struct RProp_params{T<:Real}
    eta_plus::T
    eta_minus::T
    Delta_min::T
    Delta_max::T
    Delta_0::T
    function RProp_params()
        return new{Float64}(1.2,0.5,1000.0,0.001,1.0)
    end
end

"""
    train_with_RProp(X::Array{Float64,2}, y::Array{Float64,1}, cov_func_parameters::GaussianKernelHyperparameters, MaxIter::Int, noise::Float64, params::RProp_params)
Trains kriging hyperparameters with RProp.
"""
function calibrate_by_ML_with_Rprop(X::AbstractArray{<:Real,2}, y::AbstractArray{R,1}, cov_func_parameters::GaussianKernelHyperparameters, MaxIter::Integer, noise::Real, params::RProp_params) where R<:Real
    # An implementation of this paper https://pdfs.semanticscholar.org/aa65/042ae494455a14811927eb0574871d276454.pdf
    iter = 0
    w_0 = cov_func_parameters.w_0
    w_i = cov_func_parameters.w_i
    threshold = params.Delta_min * params.eta_plus + 10*eps()
    old_0_sign = 1
    old_i_signs = ones(length(w_i))
    while iter < MaxIter
        marginal_likelihood, K, invK = marginal_likelihood_gaussian_derivatives(X, y, GaussianKernelHyperparameters(w_0, w_i), noise)
        Delta0 = params.Delta_0
        Deltas = Array{typeof(w_0),1}(undef, length(w_i)) .= params.Delta_0
        new_0_sign = sign(marginal_likelihood[1])
        new_i_signs = sign.(marginal_likelihood[2:length(marginal_likelihood)])
        w_0 = w_0 + new_0_sign*Delta0
        w_i = w_i .+ new_i_signs .* Deltas
        if iter > 0
            Delta0 = new_0_sign * old_0_sign > 0 ? Delta0*params.eta_plus : Delta0*params.eta_minus
            Delta0 = min(max(Delta0, params.Delta_min ), params.Delta_max )
            for i in 1:length(w_i)
                Deltas[i] = new_i_signs[i]*old_i_signs[i] > 0 ? Deltas[i]*params.eta_plus : Deltas[i] *params.eta_minus
                Deltas[i] = min(max(Deltas[i], params.Delta_min ), params.Delta_max )
            end
        end
        if all(Deltas .< threshold)
            return GaussianKernelHyperparameters(w_0, w_i)
        end
        old_0_sign  = new_0_sign
        old_i_signs = new_i_signs
        iter = iter + 1
    end
    return GaussianKernelHyperparameters(w_0, w_i)
end

"""
    sample(twister::MersenneTwister, dim::Integer, batch_size::Integer)
This does sampling with or without replacement.
"""
function sample(twister::MersenneTwister, dim::Integer, batch_size::Integer)
    return randperm(twister, dim)[1:batch_size]
end

"""
    calibrate_by_ML_with_SGD(X::Array{Float64,2}, y::Array{Float64,1}, cov_func_parameters::GaussianKernelHyperparameters, steps::Int, batch_size::Int, step_multiple::Float64, noise::Float64, seed::Int = 1988)
This trains a kriging model by using maximum likelihood with stochastic gradient descent.
"""
function calibrate_by_ML_with_SGD(X::AbstractArray{<:Real,2}, y::AbstractArray{R,1}, cov_func_parameters::GaussianKernelHyperparameters, steps::Integer, batch_size::Integer,
                                  step_multiple::Real = 0.05, noise::Real = 0.00001, twister::MersenneTwister = MersenneTwister(1988)) where R<:Real
    ow_0 = cov_func_parameters.w_0
    ow_i = cov_func_parameters.w_i
    ndims = length(ow_i)
    nobs = length(y)
    for s in 1:steps
        samples = sample(twister, nobs, batch_size)
        XSample = X[samples, :]
        ySample = y[samples]
        marginal_likelihood, K, invK = marginal_likelihood_gaussian_derivatives(XSample, ySample, GaussianKernelHyperparameters(ow_0, ow_i), noise)
        normalised_grad =  sign.(marginal_likelihood) .* ( abs.(marginal_likelihood) ./ (abs.(marginal_likelihood) .+ maximum(abs.(marginal_likelihood))) )
        ow_0 = ow_0 .* (1 .+ normalised_grad[1] .* step_multiple)
        ow_i = ow_i .* (1 .+ normalised_grad[2:(ndims+1)] .* step_multiple)
    end
    return GaussianKernelHyperparameters(ow_0, ow_i)
end


text:
@testset begin "AWESOME OCIM"
    AO_path = AO.download_and_unpack()
    @test isdir(AO_path)
end


text:
import Statistics
import ChainRulesCore
import ChainRulesCore: rrule, Tangent, ZeroTangent, @opt_out
import Yota: is_chainrules_primitive

loss_simple(W, b, x) = sum(W * x .+ b)
loss_double_broadcast(W, b, x) = sum(sin.(W * x) .+ b)
loss_kw_mean(W, b, x) = Statistics.mean(W * x .+ b; dims=1)[1]

function rrule(::typeof(Broadcast.broadcasted), ::typeof(sin), x)
    sin_pullback(dy) = (ZeroTangent(), ZeroTangent(), cos.(x) .* dy)
    return sin.(x), sin_pullback
end

chain_foo(x::Number) = :ok

rrule(::typeof(chain_foo), ::Number) = ZeroTangent(), dy -> ZeroTangent()
rrule(::typeof(chain_foo), ::Real) = ZeroTangent(), dy -> ZeroTangent()
rrule(::typeof(chain_foo), ::Float64) = ZeroTangent(), dy -> ZeroTangent()
@opt_out rrule(::typeof(chain_foo), ::Real)

update_chainrules_primitives!()

# TODO: uncomment when we have a general rrule() or @drule for primitive broadcasting
# @testset "grad: simple" begin
#     args3 = (rand(4, 3), rand(4), rand(3))
#     args4 = (rand(4, 3), rand(4), rand(3), rand(4))
#     @test gradcheck((W, b, x) -> sum(W * x .+ b), args3...)
#     @test gradcheck((W, b, x) -> sum(tanh.(W * x .+ b)), args3...)
#     @test gradcheck((W, b, x, y) -> sum(abs2.(tanh.(W * x .+ b) .- y)), args4...)
# end

@testset "grad: no_rrule" begin
    @test is_chainrules_primitive(Tuple{typeof(chain_foo), ComplexF64}) == true  # should hit ::Number
    @test is_chainrules_primitive(Tuple{typeof(chain_foo), Float32}) == false    # should hit ::Real
    @test is_chainrules_primitive(Tuple{typeof(chain_foo), Float64}) == true
end

@testset "grad: kw" begin
    args = (rand(3, 4), rand(3), rand(4))
    @test gradcheck(loss_simple, args...)
    @test gradcheck(loss_double_broadcast, args...)

    val, g = grad(loss_kw_mean, args...)
    @test val == loss_kw_mean(args...)
    @test gradcheck(loss_kw_mean, args...)

    @test gradcheck(x -> sum(sum(x, dims=1)), rand(2, 3))
end

@testset "grad: getindex" begin
    x = rand(3, 4, 5)
    x1 = zero(x); x1[1] = 1
    @test grad(x -> x[1], x)[2][2] == x1
    x2 = zero(x); x2[1, 2, 1] = 1
    @test grad(x -> x[1, 2, 1], x)[2][2] == x2
    x3 = zero(x); x3[:, 1, :] .= 1
    @test grad(x -> sum(x[:, 1, :]), x)[2][2] == x3
end

@testset "grad: cat" begin
    @test gradcheck((a, b) -> sum(vcat(a, b)), ones(2, 3), ones(3, 3))
    @test gradcheck((a, b) -> sum(hcat(a, b)), ones(2, 3), ones(2, 4))

    arrs = Any[i * ones(3, 4, i, 6) for i=1:4]
    @test gradcheck((a, b, c, d) -> sum(cat(a, b, c, d; dims=3)), arrs...)
end

@testset "grad: iterate" begin
    # iterate over tuple, e.g. for x in (1.0, 2.0, 3.0)
    x = (1.0, 2.0, 3.0)
    CT = Tangent{typeof(x)}
    @test grad(x -> iterate(x)[1], x)[2][2] == CT(1.0, ZeroTangent(), ZeroTangent())
    @test grad(x -> iterate(x, 2)[1], x)[2][2] == CT(ZeroTangent(), 1.0, ZeroTangent())
    @test grad(x -> iterate(x, 3)[1], x)[2][2] == CT(ZeroTangent(), ZeroTangent(), 1.0)

    # iterate over array, e.g. for x in [1.0, 2.0, 3.0]
    x = [1.0, 2.0, 3.0]
    # TODO (uncomment when scatter_add is fixed)
    # @test grad(x -> iterate(x)[1], x)[2][1] == [1.0, 0, 0]
    # @test grad(x -> iterate(x, 2)[1], x)[2][1] == [0, 1, 0]
    # @test grad(x -> iterate(x, 3)[1], x)[2][1] == [0, 0, 1.0]

    x = (1:3)
    @test grad(x -> iterate(x)[1], x)[2][2] == ZeroTangent()
    @test grad(x -> iterate(x, 1)[1], x)[2][2] == ZeroTangent()

end

sum_bcast(x, y) = sum(x .+ y)

@testset "special bcast" begin
    for args in [
        (rand(3, 4), rand(3)),
        (rand(3, 4), rand(3, 1)),
        (rand(3, 4), rand(1, 4)),
        (rand(3), rand(3, 4)),
        (rand(3, 1), rand(3, 4)),
        (rand(1, 4), rand(3, 4)),
        # (rand(3, 4), rand()),
        # (rand(), rand(3, 4)),
    ]
        @test gradcheck(sum_bcast, args...)
        val, g = grad(sum_bcast, args...)
        for i=1:length(args)
            @test size(g[i + 1]) == size(args[i])
        end
    end
end

sum_prod_bcast(x, y) = sum(x .* y)

@testset "special bcast 2" begin
    for args in [
        (rand(3, 4), rand(3)),
        (rand(3, 4), rand(3, 1)),
        (rand(3, 4), rand(1, 4)),
        (rand(3), rand(3, 4)),
        (rand(3, 1), rand(3, 4)),
        (rand(1, 4), rand(3, 4)),
        # (rand(3, 4), rand()),
        # (rand(), rand(3, 4)),
    ]
        @test gradcheck(sum_prod_bcast, args...)
        val, g = grad(sum_prod_bcast, args...)
        for i=1:length(args)
            @test size(g[i + 1]) == size(args[i])
        end
    end
end

@testset "grad: literal_pow" begin
    @test gradcheck(x -> sum(x .^ 2), rand(3, 4))
end

@testset "grad: transpose" begin
    a = rand(5)
    b = rand(2, 3)
    @test grad(x -> sum(transpose(x)), a)[2][2] == ones(size(a))
    @test grad(x -> sum(adjoint(x)), a)[2][2] == ones(size(a))
    @test grad(x -> sum(transpose(x)), b)[2][2] == ones(size(b))
    @test grad(x -> sum(adjoint(x)), b)[2][2] == ones(size(b))
end

mutable struct Linear{T}
    W::AbstractArray{T,2}
    b::AbstractArray{T}
end

forward(m::Linear, X) = m.W * X .+ m.b

loss(m::Linear, X) = sum(forward(m, X))

mutable struct Point x; y end
constructor_loss(a) = (p = Point(a, a); p.x + p.y)

@testset "grad: structs" begin
    args = Any[Linear(rand(3,4), rand(3)), rand(4,5)]
    tape = gradtape(loss, args...)

    @test tape[tape.result].val[1] == loss(args...)

    play!(tape)
    @test tape[tape.result].val[1] == loss(args...)

    val, g = grad(loss, args...)
    @test val == loss(args...)

    _, g = grad(constructor_loss, 4.0)
    @test g[2][1] == 2
end

mutable struct Line p1; p2 end

function add_points(x)
    p = Point(x, x)
    l = Line(p, p)
    return 2*l.p1.x + 5*l.p2.y
end

@testset "grad: structs/new" begin
    _, tape = trace(add_points, rand())

    @test grad(add_points, rand())[2][2] == 7
end

make_t(x) = tuple(x, x)

function make_t_loss(x)
    a, b = make_t(x)
    return a + b
end

@testset "grad: tuple unpack" begin
    _, g = grad(make_t_loss, 3.0)
    @test g[2] == 2.0
end

@testset "grad: seed" begin
    val, g = grad(x -> 2x, [1.0, 2.0, 3.0]; seed=ones(3))
    @test val == [2.0, 4.0, 6.0]
    @test g == (ZeroTangent(), [2.0, 2.0, 2.0])
end

@testset "grad: compiled" begin
    args = Any[Linear(rand(3,4), rand(3)), rand(4,5)]
    tape = gradtape(loss, args...)

    # generate new args
    args = Any[Linear(rand(3,4), rand(3)), rand(4,5)]

    # compiled vs. non-compiled
    res1 = play!(tape, loss, args...)
    grad_loss = compile(tape)
    res2 = grad_loss(loss, args...)
    @test res1 == res2

    # fresh vs. cached version
    res3 = grad(loss, args...)
    res4 = grad(loss, args...)
    @test res3 == res2
    @test res4 == res2
end


text:
module TestTSMLextra
using Test
using TSML
using TSMLextra
using TSMLextra.System

include("test_readerwriter.jl")

if LIB_CRT_AVAILABLE
    include("test_caret.jl")
else
    @info "Skipping CARET tests."
end

if LIB_SKL_AVAILABLE
    include("test_scikitlearn.jl")
    include("test_skpreprocessing.jl")
    include("test_crossvalidator.jl")
else
    @info "Skipping scikit-learn tests."
end

end

max_stars_repo_name: ppalmes/TSMLextra

max_issues_repo_name: ppalmes/TSMLextra

max_forks_repo_name: ppalmes/TSMLextra


text:
"""
    compute_adc_sta(model)

Compute the ADC for each compartment in the short diffusion time regime.
"""
function compute_adc_sta(model::Model{T}, gradient::ScalarGradient) where {T}
    ncompartment = length(model.mesh.points)

    adcs = zeros(T, ncompartment)
    for icmpt = 1:ncompartment
        points = model.mesh.points[icmpt]
        elements = model.mesh.elements[icmpt]
        facets = reduce(hcat, model.mesh.facets[icmpt, :])

        volumes, = get_mesh_volumes(points, elements)
        v = sum(volumes)
        _, areas, _, normals = get_mesh_surface(points, facets)
        d = gradient.dir
        D₀ = d' * model.D[icmpt] * d

        # Project surface areas onto plane orthogonal to gradient direction
        # (result is 1×1 matrix, extract scalar component)
        SAu = ((d'normals) .^ 2*areas)[1]

        t = diffusion_time_sta(gradient.profile)
        adcs[icmpt] = D₀ * (1 - 4 / 3 * √(D₀ / π) * SAu / v * √t)
    end

    all(≥(0), adcs) ||
        @warn "Obtained negative STA ADC, STA approximation does not hold" filter(
            <(0),
            adcs,
        )

    adcs
end


text:
#
# Recursive Fourier Propagation for BoostFractor
#
# V: 2019-08-06
#
# Stefan Knirck
#

export dancer, dance_intro, cheerleader

"""
Propagates the fields through the system

* `amin`:           Mimum (local) amplitude of a field, in order to be propagated
* `nmax`:           Maximum number of beam iteration steps, directly equivalent to how many boundaries a beam 'sees' at most
* `bdry`:           SetupBoundaries-Objekt, containing all relevant geometrical information (disk positions, epsilon, etc)
* `f`:              Frequency
* `prop`:           Propagator Function to use. Standard is propagator()
* `reflect`:        If nothing (standar value), the axion-induced signal is computed.
                    If set, this field defines a beam, for which the reflected beam will be calculated
* `Xset`, `Yset`:   Explicitly set the coordinate system for the fields
* `returnsum`:      If false, the out-propagating contributions after each iteration will be returned, without summing.
* `immediatesum`:   If false, the out-propagating contributions will be saved and summed up at the end.
"""
function dancer(amin, nmax, sbdry::SetupBoundaries, coords::CoordinateSystem; f=10.0e9, prop=propagator, emit=nothing, reflect=nothing, diskR=0.1, returnsum=true, immediatesum=true)

    # Make dancer swallow the same SetupBoundaries object as cheerleader and transformer
    bdry = deepcopy(sbdry)
    bdry.eps = bdry.eps[2:end]
    bdry.distance = bdry.distance[2:end]
    bdry.relative_tilt_x = bdry.relative_tilt_x[2:end]
    bdry.relative_tilt_y = bdry.relative_tilt_y[2:end]
    bdry.relative_surfaces = bdry.relative_surfaces[2:end,:,:]
    append!(bdry.r, 0.0)

    rightmoving = 1
    leftmoving = 2

    fields = Array{Complex{Float64}}(zeros(length(bdry.distance), 2, length(coords.X), length(coords.Y)))
    #fields = SharedArray(Complex{Float64},  length(bdry.distance), 2, length(coords.X), length(coords.Y))
    #                       number of regions -----^                ^              ^
    #                       number of propagation directions -------^              ^
    #                       dimensions of the fields at each position -------------^

    # Pre-allocate memory
    if immediatesum
        Eout = Array{Complex{Float64}}(zeros(length(coords.X), length(coords.Y),1))
    else
        Eout = Array{Complex{Float64}}(zeros(length(coords.X), length(coords.Y),nmax+1))
    end

    # TODO: propagation through and emission from last bdry to the right
    if reflect === nothing
        if emit === nothing
            fields = dance_intro(bdry,coords;diskR=diskR)
        else
            fields = emit
        end
        # Eout =
    else
        fields[length(bdry.distance), leftmoving, :, :] = reflect
        # Eout =
    end

    lambda = wavelength(f)

    n = 0   # Iteration Count
    a = 1.0 # Maximum Field Amplitude in the iteration

    # The measure of music (Takt) .... da capo al fine
    while n <= nmax && a >= amin
        # The "legs" of the dancer dancing in parallel
        # @sync is important to make sure that the outer loop waits for the inner one to finish...
        #@sync @parallel
        for i in 1:1:length(bdry.r)
            # The beats

            # Move left leg (rightmoving)
            if i > 1
                # In place propagate left leg:
                fields[i-1,rightmoving,:,:] = prop(fields[i-1,rightmoving,:,:], bdry.distance[i-1], diskR, bdry.eps[i-1], bdry.relative_tilt_x[i-1], bdry.relative_tilt_y[i-1], bdry.relative_surfaces[i-1,:,:], lambda, coords)
            end

            # Move right leg (leftmoving)
            if i < length(bdry.r)
                fields[i,leftmoving,:,:] = prop(fields[i,leftmoving,:,:], bdry.distance[i], diskR, bdry.eps[i], bdry.relative_tilt_x[i], bdry.relative_tilt_y[i], bdry.relative_surfaces[i,:,:], lambda, coords)
            end

            # Reflection & Transmission
            if i == 1 # Leftmost case
                fields[i, leftmoving,:,:] .*= -bdry.r[i]
                # The right-moving does not exist and about the transmitted one we dont care (yet)
            elseif i == length(bdry.r) # Rightmost case
                # The rightmoving is transmitted to Eout (yeah! that'll be our result!)
                # old version without pre-allocation
                #Eout                         = cat(3,Eout, (1+bdry.r[i])*fields[i-1, rightmoving,:,:])
                # with pre-allocated array:
                if immediatesum
                    Eout[:,:,1]              .+= (1+bdry.r[i]).*fields[i-1, rightmoving,:,:]
                else
                    Eout[:,:,n+1]              = (1+bdry.r[i]).*fields[i-1, rightmoving,:,:]
                end
                # And reflected as well.
                fields[i-1, rightmoving,:,:] .*= bdry.r[i]
            else # Standard Case

                # field(i-1) = transmit(field(i)) + reflect(field(i-1))
                # field(i) = transmit(field(i-1)) + reflect(field(i))
                # Basically the above, but as much in-place as possible

                # Always the field on the left gets a copy
                FieldArrivedLeft = copy(fields[i-1, rightmoving,:,:])
                # Then we in-place update it
                fields[i-1, rightmoving,:,:] .*= bdry.r[i]
                fields[i-1, rightmoving,:,:] .+= (1-bdry.r[i]).*fields[i, leftmoving,:,:]
                # Then we in-place update the other one with the previously made backup
                fields[i, leftmoving,:,:] .*= -bdry.r[i]
                fields[i, leftmoving,:,:] .+= (1+bdry.r[i]).*FieldArrivedLeft
            end
        end

        # Now the directions have changed
        lmv = copy(leftmoving)
        leftmoving = rightmoving
        rightmoving = lmv

        # a is the maximum field that still occurs in the system
        if amin != 0
            a = maximum(abs.(fields))
        end
        # n is the iteration count
        n += 1
    end

    # Dance Outro
    # The summation at the end is to give the opportunity to intelligently sum
    # since different rays might have very different orders of magnitude
    # The other reason is, that we are able to return the different components
    # such that we can get immediate access to less accurate results
    # TODO: See if Julia is really using the most intelligent summing strategy
    if returnsum
        return sum(reverse(Eout, dims = 3), dims = 3)[:,:,1]
    else
        return Eout
    end
end

"""
    dance_intro(bdry::SetupBoundaries, X, Y; bfield=nothing, velocity_x=0, f=10e9,diskR=0.1)

Initialize EM fields. Can include velocity effects.
"""
function dance_intro(bdry::SetupBoundaries, coords::CoordinateSystem; bfield=nothing, velocity_x=0, f=10e9,diskR=0.1)
    # Initialize the variable we want to return
    fields_initial = Array{Complex{Float64}}(zeros(length(bdry.distance), 2, length(coords.X), length(coords.Y)))

    # Inaccuracies of the emitted fields, BField and Velocity Effects ###################
    if bfield === nothing
        # Make sure that there is only emission on the disk surfaces
        bfield = [x^2 + y^2 > diskR^2 ? 0.0 : 1.0 for z in ones(length(bdry.distance)+1), x in coords.X, y in coords.Y]
    end
    if velocity_x != 0
        lambda = wavelength(f)
        Ma_PerMeter = 2pi/lambda # k = 2pi/lambda (c/f = lambda)
        bfield .*= [exp(-1im*Ma_PerMeter*(-velocity_x)*x) for z in ones(length(bdry.distance)+1), x in coords.X, y in coords.Y]
    end
    ####################################################################################

    # Iterate over the gaps and initialize the emissions from them #####################
    # This implements Theoretical Foundations (arxiv: 1612.07057) (3.3)
    for n in 1:length(bdry.distance)
        ax_rightmoving = 0
        if n == 1
            if bdry.r[1] == 1
                ax_rightmoving = -1.0
            else #TODO
                ax_rightmoving = 0
            end
        else
            # Right-Moving in that gap
            eps_i = bdry.eps[n-1]
            eps_m = bdry.eps[n]
            denominator = eps_i * sqrt(eps_m) + eps_m * sqrt(eps_i)
            #ax_i = sqrt(eps_i) * (1 - eps_m/eps_i) / denominator + 0im
            ax_rightmoving = sqrt(eps_m) * (1 - eps_i/eps_m) / denominator +0.0im
        end

        # Left-Moving in that gap
        eps_i = bdry.eps[n]
        eps_m = (n == length(bdry.distance)) ? 1 : bdry.eps[n+1] # Rightmost epsilon is 1.
        denominator = eps_i * sqrt(eps_m) + eps_m * sqrt(eps_i)
        ax_leftmoving = sqrt(eps_i) * (1 - eps_m/eps_i) / denominator +0.0im
        #ax_m = sqrt(eps_m) * (1 - eps_i/eps_m) / denominator +0.0im

        # Fill the initial array
        fields_initial[n,1,:,:] = [1.0*ax_rightmoving + 0.0im for x in coords.X, y in coords.Y] .* bfield[n,:,:]
        fields_initial[n,2,:,:] = [1.0*ax_leftmoving + 0.0im for x in coords.X, y in coords.Y] .* bfield[n+1,:,:]
    end

    return fields_initial
end

"""
    cheerleader(amin, nmax, bdry::SetupBoundaries; f=10.0e9, prop=propagator, emit=nothing, reflect=nothing, Xset=X, Yset=Y, diskR=0.1, returnboth=false)

New Recursive Fourier Propagation implementation.

# Arguments:
- `amin`: Mimum (local) amplitude of a field, in order to be propagated
- `nmax`: Maximum number of beam iteration steps, directly equivalent to how many boundaries a beam 'sees' at most
- `bdry::SetupBoundaries`: Properties of dielectric boundaries
- `f::Float64` ```> 0```: Frequency of EM radiation
- `prop`: Propagator Function to use. Standard is propagator().
- `emit`:  Explicitly set the axion-induced fields emitted from each boundary (to the left and to the right).
               If ``nothing`` fields are initialized according to uniform,
               homogeneous external B-field with zero axion velocity.
- `reflect`: If `nothing` (standar value), the axion-induced signal is computed.
             If set, this field defines a beam, for which the reflected beam will be calculated
- `Xset` and `Yset`: Explicitly set the coordinate system for the fields
- `diskR`: Radius of dielectric disk
- `returnboth::Bool`: If `true` cheerleader returns fields leaving on left and right.
                      If `false` only returns fields leaving on right.

See [`dancer`](@ref) for old version.
"""
function cheerleader(amin, nmax, bdry::SetupBoundaries, coords::CoordinateSystem; f=10.0e9, prop=propagator, emit=nothing, reflect=nothing, diskR=0.1, returnboth=false)

    # Before speed of light was 3e8 here, but this means an error at the permil level, i.e. order ~20MHz at 20GHz,
    # if fixing lambda to 1.5 cm, one gets a shift of roughly 10MHz
    lambda = wavelength(f)

    # Pre-allocate memory
    # Note that fields[0,:,:] contains the fields leaving the system on the left
    # and fields[length(bdry.distance)+1,:,:] the fields leaving on the right
    #fields = OffsetArray(::Array{Complex{Float64}}}, 0:length(bdry.distance)+1, 1:length(X), 1:length(Y))
    fields = Array{Complex{Float64}}(zeros(     length(bdry.distance)+2, length(coords.X), length(coords.Y)))
    # number of regions + 2 outporpagating --^                                ^
    # dimensions of the fields at each position ------------------------------^

    # In a next step this could/should be generalized in the SetupBoundaries structure..
    reflectivities_leftmoving =  -bdry.r
    reflectivities_rightmoving = bdry.r
    transmissivities_leftmoving = 1 .- bdry.r
    transmissivities_rightmoving = 1 .+ bdry.r

    ### Indexing of arrays, note the specific indexing of the regions, different from dancer()!.
    #Boundaries:     [ 1,      2,      3,      4 ]
    #          ####### |       |       |       |   #######
    #          # <-- # |       |       |       |   # --> #
    #          ####### |       |       |       |   #######
    #Regions:    [  1  ,   2   ,   3   ,   4   ,   5  ]

    # emit in this function expects for each region left- and right-propagating fields.
    # Note that this is different from dancer() since here we do not take the mirror
    # explicit, such that also a transmissivity can be calculated

    # TODO: propagation through and emission from last bdry to the right
    if reflect === nothing
        if emit === nothing
            emit = Array{Complex{Float64}}(zeros(length(bdry.distance), 2, length(coords.X), length(coords.Y)))
            # we may reuse the dance_intro-function in the standard case
            bdrycpy = deepcopy(bdry)
            bdrycpy.distance = bdrycpy.distance[2:end]
            bdrycpy.eps = bdrycpy.eps[2:end]
            emit[2:end,:,:,:] = dance_intro(bdrycpy,coords,diskR=diskR)
        end
    else
        emit = Array{Complex{Float64}}(zeros(length(bdry.distance), 2, length(coords.X), length(coords.Y)))
        # This should be right like this, it is explicitly taken care of at the bottom...
        emit[length(bdry.distance), 2, :, :] = reflect
    end

    # Initialize fields ........................................

    # Add everything which is emitted to the right
    fields = copy(emit[:,1,:,:])

    # Push forward the fields rightmoving...
    for i in 2:1:length(bdry.distance)-1
        # Propagate to the right
        fields[i,:,:] = prop(fields[i,:,:], bdry.distance[i], diskR, bdry.eps[i], bdry.relative_tilt_x[i], bdry.relative_tilt_y[i], bdry.relative_surfaces[i,:,:], lambda, coords)
        # Reflect and Transmit
        fields[i+1,:,:] .+= transmissivities_rightmoving[i].*fields[i,:,:]
        fields[i,:,:]   .*= reflectivities_rightmoving[i]
    end

    # Add everything which is emitted to the left (except last gap)
    fields[1:length(bdry.distance)-1,:,:] .+= emit[1:length(bdry.distance)-1,2,:,:]
    #fields .+= emit[:,2,:,:]

    # The last gap is not supposed to contain left-propagating stuff,
    # so we have to treat it explicitly (not elegant, but ok)
    let i = length(bdry.distance)-1
        # Transmit
        fields[i,:,:]   .+= transmissivities_leftmoving[i].*emit[i+1,2,:,:]
        # Reflect
        fields[i+1,:,:] .+= reflectivities_leftmoving[i].*emit[i+1,2,:,:]
    end

    # Main iterations ...............................................
    n = 0   # Iteration Count
    a = 1.0 # Maximum Field Amplitude in the iteration

    while n <= nmax && a >= amin
        # Iterate over regions

        # The first region explicitly
        let i = 1
            # Do NOT Propagate to the fields in the current gap to the right
            # Region 1 always contains left-moving fields.
            # Propagate to the fields in the next gap to the left
            fields[i+1,:,:] = prop(fields[i+1,:,:], bdry.distance[i+1], diskR, bdry.eps[i+1], bdry.relative_tilt_x[i+1], bdry.relative_tilt_y[i+1], bdry.relative_surfaces[i+1,:,:], lambda, coords)
            # Reflect and Transmit
            # No backup copy needed since nothing is coming from the i = 1 region
            # Transmit
            fields[i,:,:]   .+= transmissivities_leftmoving[i].*fields[i+1,:,:]
            # Reflect
            fields[i+1,:,:] .*= reflectivities_leftmoving[i]
        end

        # Push forward...
        for i in 2:1:length(bdry.distance)-2
            # Propagate to the fields in the current gap to the right
            fields[i,:,:] = prop(fields[i,:,:], bdry.distance[i], diskR, bdry.eps[i], bdry.relative_tilt_x[i], bdry.relative_tilt_y[i], bdry.relative_surfaces[i,:,:], lambda, coords)
            # Propagate to the fields in the next gap to the left
            fields[i+1,:,:] = prop(fields[i+1,:,:], bdry.distance[i+1], diskR, bdry.eps[i+1], bdry.relative_tilt_x[i+1], bdry.relative_tilt_y[i+1], bdry.relative_surfaces[i+1,:,:], lambda, coords)
            # Reflect and Transmit
            fields_next_gap_copy = copy(fields[i+1,:,:])
            fields[i+1,:,:] .*= reflectivities_leftmoving[i]
            fields[i+1,:,:] .+= transmissivities_rightmoving[i].*fields[i,:,:]
            fields[i,:,:]   .*= reflectivities_rightmoving[i]
            fields[i,:,:]   .+= transmissivities_leftmoving[i].*fields_next_gap_copy
        end

        let i = length(bdry.distance)-1
            # Propagate to the fields in the current gap to the right
            fields[i,:,:] = prop(fields[i,:,:], bdry.distance[i], diskR, bdry.eps[i], bdry.relative_tilt_x[i], bdry.relative_tilt_y[i], bdry.relative_surfaces[i,:,:], lambda, coords)
            # DO NOT Propagate to the fields in the next gap to the left
            # Since it only contains fields propagating to the right
            # Reflect and Transmit
            fields[i+1,:,:] .+= transmissivities_rightmoving[i].*fields[i,:,:]
            fields[i,:,:]   .*= reflectivities_rightmoving[i]
        end

        # The last region always contains right-moving fields and is automatically added up correctly

        # Check convergence
        if amin != 0
        	a = maximum(abs.(fields))
		end
        # n is the iteration count
        n += 1
    end

    if returnboth
        return fields[1,:,:], fields[length(bdry.distance),:,:]
    else
        return fields[length(bdry.distance),:,:]
    end
end