text:
__precompile__()
module Decomp

using StaticArrays: MVector, SVector
export eigen!, eigen, inphdecomp, propdecomp, inphmodel, propmodel, mixmodel

 function eigen!(e11::T,e22::T,e33::T,e12::T,e13::T,e23::T,eig::MVector{3,T},eigv1::MVector{3,T},eigv2::MVector{3,T},eigv3::MVector{3,T}) where T<:AbstractFloat

  p1 = muladd(e12, e12, muladd(e13, e13, e23*e23))
  if p1 == 0.0
    # E is diagonal.
    # eig[1] = e11
    # eig[2] = e22
    # eig[3] = e33

    if abs(e11) >= abs(e22)
      if abs(e11) >= abs(e33)
        eig[1]=e11

        eigv1[1] = oneunit(T)
        eigv1[2] = zero(T)
        eigv1[3] = zero(T)

        if abs(e22) >= abs(e33)
          eig[2] = e22
          eig[3] = e33

          eigv2[1] = zero(T)
          eigv2[2] = oneunit(T)
          eigv2[3] = zero(T)
        else
          eig[3] = e22
          eig[2] = e33

          eigv2[1] = zero(T)
          eigv2[2] = zero(T)
          eigv2[3] = oneunit(T)
        end
      else
        eig[1] = e33
        eig[2] = e11
        eig[3] = e22

        eigv1[1] = zero(T)
        eigv1[2] = zero(T)
        eigv1[3] = oneunit(T)

        eigv2[1] = oneunit(T)
        eigv2[2] = zero(T)
        eigv2[3] = zero(T)
      end
    else
      if abs(e22) >= abs(e33)
        eig[1] = e22

        eigv1[1] = zero(T)
        eigv1[2] = oneunit(T)
        eigv1[3] = zero(T)

      elseif abs(e33) >= abs(e11)
        eig[2] = e33
        eig[3] = e11

        eigv2[1] = zero(T)
        eigv2[2] = zero(T)
        eigv2[3] = oneunit(T)

      else
        eig[3] = e33
        eig[2] = e11

        eigv2[1] = oneunit(T)
        eigv2[2] = zero(T)
        eigv2[3] = zero(T)
      end
    end

    eigv3[1] = eigv1[2]*eigv2[3] - (eigv1[3]*eigv2[2])
    eigv3[2] = eigv1[3]*eigv2[1] - (eigv1[1]*eigv2[3])
    eigv3[3] = eigv1[1]*eigv2[2] - (eigv1[2]*eigv2[1])
  else
    q = (e11 + e22 + e33)/3
    p2 = (e11-q)^2 + (e22-q)^2 + (e33-q)^2 + 2*p1
    p = sqrt(p2/6)
    r = ((e11-q)*(e22-q)*(e33-q) - (e11-q)*(e23^2) - (e12^2)*(e33-q) + 2*(e12*e13*e23) - (e13^2)*(e22-q))/(2*p*p*p)

    # In exact arithmetic for a symmetric matrix  -1 <= r <= 1
    # but computation error can leave it slightly outside this range.
    if r <= -1
      ϕ = T(π/3)
    elseif r >= 1
      ϕ = zero(T)
    else
      ϕ = acos(r)/3
    end

    # the eigenvalues satisfy eig[3] <= eig[2] <= eig[1]
    eig[1] = q + 2*p*cos(ϕ)
    eig[3] = q + 2*p*cos(ϕ+(2*π/3))
    eig[2] = 3*q - eig[1] - eig[3]     # since trace(E) = eig[1] + eig[2] + eig[3] = 3q

    #sort!(eig,by=abs,rev=true)
    if abs(eig[1]) < abs(eig[3])
      big = eig[3]
      if abs(eig[1]) < abs(eig[2])
        eig[3] = eig[1]
        eig[1] = big
      else
        eig[3] = eig[2]
        eig[2] = eig[1]
        eig[1] = big
      end
    else
      if abs(eig[2]) < abs(eig[3])
        medium = eig[3]
        eig[3] = eig[2]
        eig[2] = medium
      end
    end

    bla = ((e22 - eig[1])*(e33 - eig[1]) - e23*e23)
    if bla != 0
      eigv1[1] = oneunit(T)
      eigv1[2] = (e23*e13 - (e33-eig[1])*e12)/bla
      eigv1[3] = (-e13 -e23*eigv1[2])/(e33-eig[1])
      aux = sqrt(1 + eigv1[2]^2 + eigv1[3]^2)
      eigv1[1] = 1/aux
      eigv1[2] = eigv1[2]/aux
      eigv1[3] = eigv1[3]/aux
    else
      bla = ((e11 - eig[1])*(e22 - eig[1]) - e12*e12)
      eigv1[3] = oneunit(T)
      eigv1[1] = (e23*e12 - (e22-eig[1])*e13)/bla
      eigv1[2] = (-e23 -e12*eigv1[1])/(e22-eig[1])
      aux = sqrt(1 + eigv1[2]^2 + eigv1[1]^2)
      eigv1[1] = eigv1[1]/aux
      eigv1[2] = eigv1[2]/aux
      eigv1[3] = 1/aux
    end
    bla = ((e22 - eig[2])*(e33 - eig[2]) - e23*e23)
    if bla != 0
      eigv2[1] = oneunit(T)
      eigv2[2] = (e23*e13 - (e33-eig[2])*e12)/bla
      eigv2[3] = (-e13 -e23*eigv2[2])/(e33-eig[2])
      aux = sqrt(1 + eigv2[2]^2 + eigv2[3]^2)
      eigv2[1] = 1/aux
      eigv2[2] = eigv2[2]/aux
      eigv2[3] = eigv2[3]/aux
    else
      bla = ((e11 - eig[2])*(e22 - eig[2]) - e12*e12)
      eigv2[3] = oneunit(T)
      eigv2[1] = (e23*e12 - (e22-eig[2])*e13)/bla
      eigv2[2] = (-e23 -e12*eigv2[1])/(e22-eig[2])
      aux = sqrt(1 + eigv2[2]^2 + eigv2[1]^2)
      eigv2[1] = eigv2[1]/aux
      eigv2[2] = eigv2[2]/aux
      eigv2[3] = 1.0/aux
    end

    eigv3[1] = eigv1[2]*eigv2[3] - (eigv1[3]*eigv2[2])
    eigv3[2] = eigv1[3]*eigv2[1] - (eigv1[1]*eigv2[3])
    eigv3[3] = eigv1[1]*eigv2[2] - (eigv1[2]*eigv2[1])

  end
  return 0
end

function eigen(e11::T,e22::T,e33::T,e12::T,e13::T,e23::T) where T<:AbstractFloat
  eig = MVector{3,T}()
  eigv1 = MVector{3,T}()
  eigv2 = MVector{3,T}()
  eigv3 = MVector{3,T}()
  eigen!(e11,e22,e33,e12,e13,e23,eig,eigv1,eigv2,eigv3)
  return SVector(eig), SVector(eigv1), SVector(eigv2), SVector(eigv3)
end

function eigen(e11::Real,e22::Real,e33::Real,e12::Real,e13::Real,e23::Real)
  return eigen(Float64(e11),Float64(e22),Float64(e33),Float64(e12),Float64(e13),Float64(e23))
end

function eigen(A::AbstractArray{<:Real,2})
  size(A) == (3,3) || throw(ArgumentError("Input must be a 3x3 Matrix"))
  (A[1,2] == A[2,1] && A[1,3] == A[3,1] && A[2,3] == A[3,2]) || throw(ArgumentError("Input must be symmetric matrix"))
  return eigen(A[1,1],A[2,2],A[3,3],A[1,2],A[1,3],A[2,3])
end

 function inphdecomp!(t11::T,t22::T,t33::T,t12::T,t13::T,t23::T,
                e11::T,e22::T,e33::T,e12::T,e13::T,e23::T,
                eig::MVector{3,T},eigv1::MVector{3,T},eigv2::MVector{3,T},eigv3::MVector{3,T}) where T<:AbstractFloat

  eigen!(e11,e22,e33,e12,e13,e23,eig, eigv1, eigv2, eigv3)

  m1 = muladd(t11, eigv1[1]*eigv1[1], muladd(t22, eigv1[2]*eigv1[2], muladd(t33, eigv1[3]*eigv1[3],
       2*muladd(t12, eigv1[1]*eigv1[2], muladd(t13, eigv1[1]*eigv1[3], t23*eigv1[2]*eigv1[3])))))
  m2 = muladd(t11, eigv2[1]*eigv2[1], muladd(t22, eigv2[2]*eigv2[2], muladd(t33, eigv2[3]*eigv2[3],
       2*muladd(t12, eigv2[1]*eigv2[2], muladd(t13, eigv2[1]*eigv2[3], t23*eigv2[2]*eigv2[3])))))
  m3 = t11 + t22 + t33 - m1 - m2

  lam2 = muladd(eig[1], eig[1], muladd(eig[2], eig[2], eig[3]*eig[3]))

  a = 1 - (3*(eig[1]^2)/lam2)
  b = 1 - (3*(eig[2]^2)/lam2)
  c = b*eig[1] - a*eig[2]

  alpha1 = (b*m1 - a*m2)/c
  alpha0 = (m1 - eig[1]*alpha1)/a
  alpha2 = -3*alpha0/lam2

  tm11 = muladd(m1, eigv1[1]*eigv1[1], muladd(m2, eigv2[1]*eigv2[1], m3*eigv3[1]*eigv3[1]))
  tm22 = muladd(m1, eigv1[2]*eigv1[2], muladd(m2, eigv2[2]*eigv2[2], m3*eigv3[2]*eigv3[2]))
  tm33 = muladd(m1, eigv1[3]*eigv1[3], muladd(m2, eigv2[3]*eigv2[3], m3*eigv3[3]*eigv3[3]))
  tm12 = muladd(m1, eigv1[1]*eigv1[2], muladd(m2, eigv2[1]*eigv2[2], m3*eigv3[1]*eigv3[2]))
  tm13 = muladd(m1, eigv1[1]*eigv1[3], muladd(m2, eigv2[1]*eigv2[3], m3*eigv3[1]*eigv3[3]))
  tm23 = muladd(m1, eigv1[2]*eigv1[3], muladd(m2, eigv2[2]*eigv2[3], m3*eigv3[2]*eigv3[3]))

  return tm11, tm22, tm33, tm12, tm13, tm23, alpha0, alpha1, alpha2
end

function inphdecomp(t11::T,t22::T,t33::T,t12::T,t13::T,t23::T,
                e11::T,e22::T,e33::T,e12::T,e13::T,e23::T) where T<:AbstractFloat

  eig = MVector{3,T}()
  eigv1 = MVector{3,T}()
  eigv2 = MVector{3,T}()
  eigv3 = MVector{3,T}()

  return inphdecomp!(t11,t22,t33,t12,t13,t23,e11,e22,e33,e12,e13,e23,eig,eigv1,eigv2,eigv3)
end

inphdecomp(var...) = inphdecomp((Float64(v) for v in var)...)

 function propdecomp(t11::T,t22::T,t33::T,t12::T,t13::T,t23::T,
                    e11::T,e22::T,e33::T,e12::T,e13::T,e23::T) where T<:AbstractFloat

  mode2 = (e11^2 + 2*e12^2 + 2*e13^2 + e22^2 + 2*e23^2 + e33^2)
  alphat = (e11*t11 + 2*e12*t12 + 2*e13*t13 + e22*t22 + 2*e23*t23 + e33*t33)/mode2

  m12 = alphat*e12
  m13 = alphat*e13
  m23 = alphat*e23
  m11 = alphat*e11
  m22 = alphat*e22
  m33 = alphat*e33

  return m11, m22, m33, m12, m13, m23, alphat
end

propdecomp(var...) = propdecomp((Float64(v) for v in var)...)

 function inphmodel(t11::Vector{T},t22::Vector{T},t33::Vector{T},t12::Vector{T},t13::Vector{T},t23::Vector{T}, #input
                   e11::Vector{T},e22::Vector{T},e33::Vector{T},e12::Vector{T},e13::Vector{T},e23::Vector{T}, #input
                   tm11::Vector{T},tm22::Vector{T},tm33::Vector{T},tm12::Vector{T},tm13::Vector{T},tm23::Vector{T}, #output
                   alpha0::Vector{T},alpha1::Vector{T},alpha2::Vector{T},ratio::Vector{T}) where T<:Real

  eigg  = [MVector{3,T}() for i in 1:Threads.nthreads()]
  eigv1 = [MVector{3,T}() for i in 1:Threads.nthreads()]
  eigv2 = [MVector{3,T}() for i in 1:Threads.nthreads()]
  eigv3 = [MVector{3,T}() for i in 1:Threads.nthreads()]

  Threads.@threads for i in 1:length(e11)
    j = Threads.threadid()
    @inbounds tm11[i], tm22[i], tm33[i], tm12[i], tm13[i], tm23[i], alpha0[i], alpha1[i], alpha2[i] = inphdecomp!(t11[i],t22[i],t33[i],t12[i],t13[i],t23[i],
                            e11[i],e22[i],e33[i],e12[i],e13[i],e23[i],
                            eigg[j],eigv1[j],eigv2[j],eigv3[j])

    @inbounds ratio[i] = (tm11[i]^2 + 2*tm12[i]^2 + 2*tm13[i]^2 + tm22[i]^2 + 2*tm23[i]^2 + tm33[i]^2)/(t11[i]^2 + 2*t12[i]^2 + 2*t13[i]^2 + t22[i]^2 + 2*t23[i]^2 + t33[i]^2)
  end
  return 0
end

 function propmodel(t11::A,t22::A,t33::A,t12::A,t13::A,t23::A, #input
                   e11::A,e22::A,e33::A,e12::A,e13::A,e23::A, #input
                   tm11::A,tm22::A,tm33::A,tm12::A,tm13::A,tm23::A, #output
                   alpha::A,ratio::A) where {A<:Vector{<:Real}}

  Threads.@threads for i in 1:length(e11)
    @inbounds tm11[i], tm22[i], tm33[i], tm12[i], tm13[i], tm23[i], alpha[i] = propdecomp(t11[i],t22[i],t33[i],t12[i],t13[i],t23[i],e11[i],e22[i],e33[i],e12[i],e13[i],e23[i])

    @inbounds ratio[i] = (tm11[i]^2 + 2*tm12[i]^2 + 2*tm13[i]^2 + tm22[i]^2 + 2*tm23[i]^2 + tm33[i]^2)/(t11[i]^2 + 2*t12[i]^2 + 2*t13[i]^2 + t22[i]^2 + 2*t23[i]^2 + t33[i]^2)
  end
  return 0
end

 function mixmodel(t11::A,t22::A,t33::A,t12::A,t13::A,t23::A, #input
                  e11::A,e22::A,e33::A,e12::A,e13::A,e23::A, #input
                  c11::A,c22::A,c33::A,c12::A,c13::A,c23::A, #input
                  tm11::A,tm22::A,tm33::A,tm12::A,tm13::A,tm23::A, #output
                  ratio::A) where {A<:Vector{<:Real}}

  Threads.@threads for i in 1:length(e11)
    @inbounds tm11[i] = e11[i] + c11[i]
    @inbounds tm22[i] = e22[i] + c22[i]
    @inbounds tm33[i] = e33[i] + c33[i]
    @inbounds tm12[i] = e12[i] + c12[i]
    @inbounds tm13[i] = e13[i] + c13[i]
    @inbounds tm23[i] = e23[i] + c23[i]

    @inbounds ratio[i] = (tm11[i]^2 + 2*tm12[i]^2 + 2*tm13[i]^2 + tm22[i]^2 + 2*tm23[i]^2 + tm33[i]^2)/(t11[i]^2 + 2*t12[i]^2 + 2*t13[i]^2 + t22[i]^2 + 2*t23[i]^2 + t33[i]^2)
  end
  return 0
end

end # module


text:
module CUSOLVER

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, cuComplex, cuDoubleComplex, libraryPropertyType, cudaDataType
using ..CUDA: libcusolver, libcusolvermg, @allowscalar, assertscalar, unsafe_free!, @retry_reclaim, @context!

using ..CUBLAS: cublasFillMode_t, cublasOperation_t, cublasSideMode_t, cublasDiagType_t
using ..CUSPARSE: cusparseMatDescr_t

using CEnum

using Memoization

using DataStructures

# core library
include("libcusolver_common.jl")
include("error.jl")
include("libcusolver.jl")

# low-level wrappers
include("base.jl")
include("sparse.jl")
include("dense.jl")
include("multigpu.jl")

# high-level integrations
include("linalg.jl")

# cache for created, but unused handles
const idle_dense_handles = HandleCache{CuContext,cusolverDnHandle_t}()
const idle_sparse_handles = HandleCache{CuContext,cusolverSpHandle_t}()

function dense_handle()
    state = CUDA.active_state()
    handle, stream = get!(task_local_storage(), (:CUSOLVER, :dense, state.context)) do
        new_handle = pop!(idle_dense_handles, state.context) do
            cusolverDnCreate()
        end

        finalizer(current_task()) do task
            push!(idle_dense_handles, state.context, new_handle) do
                @context! skip_destroyed=true state.context cusolverDnDestroy()
            end
        end

        cusolverDnSetStream(new_handle, state.stream)

        new_handle, state.stream
    end::Tuple{cusolverDnHandle_t,CuStream}

    if stream != state.stream
        cusolverDnSetStream(handle, state.stream)
        task_local_storage((:CUSOLVER, :dense, state.context), (handle, state.stream))
    end

    return handle
end

function sparse_handle()
    state = CUDA.active_state()
    handle, stream = get!(task_local_storage(), (:CUSOLVER, :sparse, state.context)) do
        new_handle = pop!(idle_sparse_handles, state.context) do
            cusolverSpCreate()
        end

        finalizer(current_task()) do task
            push!(idle_sparse_handles, state.context, new_handle) do
                @context! skip_destroyed=true state.context cusolverSpDestroy(new_handle)
            end
        end

        cusolverSpSetStream(new_handle, state.stream)

        new_handle, state.stream
    end::Tuple{cusolverSpHandle_t,CuStream}

    if stream != state.stream
        cusolverSpSetStream(handle, state.stream)
        task_local_storage((:CUSOLVER, :sparse, state.context), (handle, state.stream))
    end

    return handle
end

function mg_handle()
    ctx = context()
    get!(task_local_storage(), (:CUSOLVER, :mg, ctx)) do
        # we can't reuse cusolverMg handles because they can only be assigned devices once
        handle = cusolverMgCreate()

        finalizer(current_task()) do task
            tls = task.storage
            if haskey(tls, (:CUSOLVER, :mg, ctx)) && CUDA.isvalid(ctx)
                # look-up the handle again, because it might have been destroyed already
                # (e.g. in `devices!`)
                handle = tls[(:CUSOLVER, :mg, ctx)]
                cusolverMgDestroy(handle)
            end
        end

        # select devices
        cusolverMgDeviceSelect(handle, ndevices(), devices())

        handle
    end::cusolverMgHandle_t
end

# module-local version of CUDA's devices/ndevices to support flexible device selection
# TODO: make this task-local
const __devices = Cint[0]
devices() = __devices
ndevices() = length(__devices)

# NOTE: this invalidates any existing cusolverMg handle
function devices!(devs::Vector{CuDevice})
    resize!(__devices, length(devs))
    __devices .= deviceid.(devs)

    # we can't select different devices _after_ having initialized the handle,
    # so just destroy it and wait for initialization to kick in again
    ctx = context()
    if haskey(task_local_storage(), (:CUSOLVER, :mg, ctx))
        handle = task_local_storage((:CUSOLVER, :mg, ctx))
        cusolverMgDestroy(handle)
        delete!(task_local_storage(), (:CUSOLVER, :mg, ctx))
    end

    return
end

end

ext: jl


text:
using Glob
if split(pwd(),'/')[end]!="test"
  try
     cd("test")
  catch
     @error("Could not find test directory from " * pwd() )
  end
end

using Pkg
if haskey(Pkg.installed(),"NBInclude")
   nbinclude_is_installed = true
   using NBInclude
else
   nbinclude_is_installed = false
   macro nbinclude(name)
         return :( @warn("NBInclude not installed, not parsing " * $name * "."); false )
   end
end

if haskey(Pkg.installed(),"Weave")
   weave_is_installed = true
   using Weave
else
   weave_is_installed = false
   @warn("Weave not installed, not parsing .jmd files.")
end

test_files = glob("test*.jl")
println("# test_files=",test_files)
for fn_test in test_files
    println("# Processing test file ",fn_test,"...")
    fn_ex_jl = replace(fn_test,r"test(\d)\." => s"../ex\1.")
    fn_ex_jmd = replace(fn_ex_jl,r"\.jl$" => ".jmd")
    fn_ex_nb = replace(fn_ex_jl,r"\.jl$" => ".ipynb")
    if nbinclude_is_installed && isfile(fn_ex_nb)
       println("# Parsing ",fn_ex_nb, "..." )
       @nbinclude(fn_ex_nb)
    elseif weave_is_installed && isfile(fn_ex_jmd)
       println("# Tangling ", fn_ex_jmd, " into ", fn_ex_jl, "..." )
       tangle(fn_ex_jmd)
       include(fn_ex_jl)
    elseif isfile(fn_ex_jl)
       include(fn_ex_jl)
    else
       println("# Can not parse either ", fn_ex_nb, " or ", fn_ex_jmd, " or ", fn_ex_jl )
       false
    end
    println("# Running test file ",fn_test,"...")
    include(fn_test)
end


text:
# initialisation functions for GeMM

"""
    createpop(cellsize)

Create a new, random individual and replicates it a certain number of times
(depending on metabolic variables) to create a new population of organisms.
Returns an array of individuals.
"""
function createpop(cellsize::Float64)
    # start with an archetypical individual
    archetype = createind(true)
    population = [archetype]
    # calculate the population size
    popsize = 0
    if occursin("metabolic", setting("popsize")) || occursin("single", setting("popsize"))
        # population size determined by adult size
        popsize = round(setting("fertility") * archetype.traits["repsize"] ^ (-1 / 4) *
                        exp(-act / (boltz * 298.0)))
    elseif occursin("bodysize", setting("popsize"))
        # population size up to 25% of the maximum possible in this cell
        quarterpopsize = Integer(ceil(cellsize / archetype.traits["repsize"]) / 4)
        popsize = rand(1:quarterpopsize) + 1
    elseif occursin("minimal", setting("popsize")) || popsize == 0
        popsize = 2 #Takes two to tangle ;-)
    else
        simlog("Invalid value for `popsize`: $(setting("popsize"))", 'e')
    end
    # clone & mutate the archetype N times
    locivar = rand()
    for i in 2:popsize
        #TODO avoid `deepcopy` here
        ind = deepcopy(archetype) # XXX this is potentially very slow!
        ind.id = rand(UInt64)
        varyalleles!(ind.genome, locivar)
        ind.traits = gettraitdict(ind.genome, setting("traitnames"))
        if !(setting("indsize") == "adult" || setting("indsize") == "seed")
            #XXX sizes uniformally distributed?
            ind.size = ind.traits["seedsize"] + rand() * (ind.traits["repsize"] - ind.traits["seedsize"])
        end
        push!(population, ind)
    end
    population
end

"""
    createind(marked=false)

Create an individual organism of a new species with a random genome.
"""
function createind(marked::Bool = false)
    id = rand(UInt64)
    lineage = randstring(6)
    ngenes = rand(1:setting("maxloci")) * length(setting("traitnames"))
    setting("degpleiotropy") == 0 && (ngenes = length(setting("traitnames")))
    ngenes < 1 && (simlog("0 genes would be created", 'e')) #shouldn't happen
    traits = createtraits()
    genes = creategenes(ngenes, traits)
    if setting("linkage") == "none"
        nchrms = length(genes)
    elseif setting("linkage") == "full"
        nchrms = 1
    else
        linkage = rand(1:length(genes))
        nchrms = Integer(round(ngenes/linkage))
    end
    chromosomes = createchrms(nchrms, genes, lineage)
    varyalleles!(chromosomes, rand())
    traitdict = gettraitdict(chromosomes, setting("traitnames"))
    if setting("indsize") == "adult"
        indsize = traitdict["repsize"]
    elseif setting("indsize") == "seed"
        indsize = traitdict["seedsize"]
    else
        indsize = traitdict["seedsize"] + rand() * (traitdict["repsize"] - traitdict["seedsize"])
    end
    Individual(lineage, chromosomes, traitdict, marked, 1.0, 1.0, indsize, hermaphrodite, 0, id)
end

"""
    genesis(cellsize)

Create a new community, composed of random new species populations, for a patch.
Returns an array of individuals.
"""
function genesis(cellsize::Float64)
    community = Individual[]
    totalmass = 0.0
    while true
        population = createpop(cellsize)
        popsize = length(population)
        # Check the cell capacity
        popmass = sum(map(x -> x.size, population))
        if totalmass + popmass > cellsize * setting("overfill")
            # stop loop if cell is full
            if totalmass >= cellsize * 0.9 || occursin("single", setting("popsize"))
                #make sure the cell is full enough
                @simlog("Cell is now $(round((totalmass/cellsize)*100))% full.", 'd')
                break
            else
                continue
            end
        end
        totalmass += popmass
        append!(community, population)
        occursin("single", setting("popsize")) && break
    end
    @simlog("Patch initialized with $(length(community)) individuals.", 'd')
    community
end

"""
    createpatch(patchentry)

Create a new patch and set its parameters by parsing one line of the map file.
"""
function createpatch(patchentry::Array{String,1})
    size(patchentry,1) < 3 && simlog("please check your map file for incomplete or faulty entries. \n
Each line must contain patch information with at least \n
\t - a unique integer ID, \n
\t - an integer x coordinate, \n
\t - an integer y coordinate, \n
separated by a whitespace character (<ID> <x> <y>).", 'e')
    # create the basic patch
    id = parse(Int, patchentry[1])
    xcord = parse(Int, patchentry[2])
    ycord = parse(Int, patchentry[3])
    defaultcapacity = setting("cellsize")
    # XXX the 'global' here is a hack so that I can use eval() later on
    # (eval() always works on the global scope)
    global newpatch = Patch(id, (xcord, ycord), defaultcapacity)
    @simlog("Creating patch $id at $xcord/$ycord", 'd')
    # parse other parameter options
    for p in patchentry[4:end]
        varval = split(p, '=')
        var = varval[1]
        if !(var in map(string, fieldnames(Patch)))
            simlog("Unrecognized patch parameter $var.", 'w')
            continue
        elseif length(varval) < 2
            # if no value is specified, assume 'true'
            # (unless the variable is prepended with an exclamation mark, then take 'false')
            if var[1] == '!'
                var = var[2:end]
                val = false
            else
                val = true
            end
        else
            val = Meta.parse(varval[2])
        end
        # check for correct type and modify the new patch
        vartype = typeof(eval(Meta.parse("newpatch."*var)))
        if !isa(val, vartype)
            try
                val = convert(vartype, val)
            catch
                simlog("Invalid patch parameter type $var: $val", 'w')
                continue
            end
        end
        eval(Meta.parse("newpatch."*string(var)*" = $val"))
    end
    newpatch
end

"""
    createworld(maptable)

Use a parsed map file (as returned by `readmapfile`) to create the world.
Initialises each patch with its parameters and a new community, then returns
an array of patches.
"""
function createworld(maptable::Array{Array{String,1},1})
    simlog("Creating world...")
    world = Array{Patch}(undef, length(maptable))
    for entry in eachindex(maptable)
        newpatch = createpatch(maptable[entry])
        if newpatch.initpop
            if setting("mode") == "zosterops"
                append!(newpatch.community, zgenesis(newpatch))
            elseif setting("indsize") != "seed"
                append!(newpatch.community, genesis(newpatch.capacity))
            else
                append!(newpatch.seedbank, genesis(newpatch.capacity))
            end
        end
        world[entry] = newpatch
    end
    global newpatch = nothing # remove variable used in `createpatch()`
    (setting("mode") == "zosterops") && findneighbours!(world)
    world
end

"""
    updateworld(world, maptable, cellsize)

Reinitialise the world from another parsed map file. Works analogously to
`createworld`. Intended for use in scenarios where the model world changes
during a run (e.g. through global warming or island ontogeny).
If the ID of an old patch matches the ID of a new patch, its community and
seed bank are transferred to the new patch.
"""
function updateworld(oldworld::Array{Patch,1}, maptable::Array{Array{String,1},1})
    simlog("Updating world...")
    newworld = Array{Patch}(undef, length(maptable))
    for entry in eachindex(maptable)
        newpatch = createpatch(maptable[entry])
        updates = findfirst(p -> p.id==newpatch.id, oldworld)
        if !isnothing(updates) # did this patch exist in the old world?
            newpatch.community = oldworld[updates].community
            newpatch.seedbank = oldworld[updates].seedbank
        end
        newworld[entry] = newpatch
    end
    global newpatch = nothing # remove variable used in `createpatch()`
    (setting("mode") == "zosterops") && findneighbours!(newworld)
    return newworld
end

max_stars_repo_name: CCTB-Ecomods/gemm

max_issues_repo_name: CCTB-Ecomods/gemm

max_forks_repo_name: CCTB-Ecomods/gemm


text:
monattick = ["Jan", "Feb", "Mär", "Apr", "Mai", "Jun", "Jul", "Aug", "Sep", "Okt", "Nov", "Dez" ]

## für Bundesländer
Bundesländer = ["Schleswig-Holstein",
                "Hamburg",
                "Niedersachsen",
                "Bremen",
                "Nordrhein-Westfalen",
                "Hessen",
                "Rheinland-Pfalz",
                "Baden-Württemberg",
                "Bayern",
                "Saarland",
                "Berlin",
                "Brandenburg",
                "Mecklenburg-Vorpommern",
                "Sachsen",
                "Sachsen-Anhalt",
                "Thüringen"]

# https://www.orte-in-deutschland.de/amtlicher-gemeindeschluessel-ags.html
"""
    AGS

Are used in vacc data.
"""
AGS = Dict(
    17 => "Bundesressorts",
    1 => "Schleswig-Holstein",
    2 => "Hamburg",
    3 => "Niedersachsen",
    4 => "Bremen",
    5 => "Nordrhein-Westfalen",
    6 => "Hessen",
    7 => "Rheinland-Pfalz",
    8 => "Baden-Württemberg",
    9 => "Bayern",
    10 => "Saarland",
    11 => "Berlin",
    12 => "Brandenburg",
    13 => "Mecklenburg-Vorpommern",
    14 => "Sachsen",
    15 => "Sachsen-Anhalt",
    16 => "Thüringen",
)

# Write your package code here.
function __init__()
# begin
    register(DataDep("Sterbestatistik 2016-2021","""
    Autor: Statistisches Bundesamt
    Dataset: Sterbefälle - Fallzahlen nach Tagen, Wochen, Monaten, Altersgruppen, Geschlecht und Bundesländern für Deutschland 2016 - 2021
    Website: https://www.destatis.de/DE/Themen/Gesellschaft-Umwelt/Bevoelkerung/Sterbefaelle-Lebenserwartung/Tabellen/sonderauswertung-sterbefaelle.html

    Diese Sonderauswertung enthält Sterbefallzahlen nach Tagen, Wochen und Monaten seit dem 1. Januar 2016.
    Die Auswertung für die Jahre 2016 bis 2020 basiert dabei auf den endgültigen plausibilisierten Daten dieser Berichtsjahre.
    Daten ab dem 1. Januar 2021 sind vorläufig – hierbei handelt es sich zunächst um eine reine Fallzahlauszählung der eingegangenen Sterbefallmeldungen aus den Standesämtern – ohne die übliche statistische Aufbereitung.
    Die Daten wurden nicht plausibilisiert und es wurde keine Vollständigkeitskontrolle durchgeführt.
    Um dennoch möglichst genaue, schnelle und vergleichbare Daten bereitzustellen, wird ein Schätzverfahren zur Hochrechnung unvollständiger Sterbefallmeldungen für die jeweils aktuellsten Daten eingesetzt.
    Im Vergleich zu den endgültigen Daten liegt zudem nur ein begrenzter Merkmalsumfang vor (Rohdaten).
    """,
                     "https://www.destatis.de/DE/Themen/Gesellschaft-Umwelt/Bevoelkerung/Sterbefaelle-Lebenserwartung/Tabellen/sonderauswertung-sterbefaelle.xlsx?__blob=publicationFile"
                     ))

    register(DataDep("Sterbestatistik 2000-2015","""
    Autor: Statistisches Bundesamt
    Dataset: Sterbefälle - Fallzahlen nach Tagen, Wochen, Monaten, Altersgruppen, Geschlecht und Bundesländern für Deutschland 2000 - 2015
    Website: https://www.destatis.de/DE/Themen/Gesellschaft-Umwelt/Bevoelkerung/Sterbefaelle-Lebenserwartung/Tabellen/sonderauswertung-sterbefaelle-endgueltige-daten.html;jsessionid=8C5CBDB9DF2D7285D72634AAF9834457.live711?nn=209016

    Diese Sonderauswertung enthält Sterbefallzahlen nach Tagen, Kalenderwochen, Monaten, Altersgruppen, Geschlecht und Bundesländern der Jahre 2000 bis 2015. Für diese Sonderauswertung stand teilweise das Originalmaterial nicht mehr zur Verfügung. Deshalb gibt es in den Jahren 2000, 2002, 2003, 2005, 2006 und 2008 geringfügige Abweichungen gegenüber den an anderer Stelle veröffentlichten endgültigen Ergebnissen. Die Aussagekraft der hier präsentierten Daten wird dadurch nicht eingeschränkt.
    """,
                     "https://www.destatis.de/DE/Themen/Gesellschaft-Umwelt/Bevoelkerung/Sterbefaelle-Lebenserwartung/Tabellen/sonderauswertung-sterbefaelle-endgueltige-daten.xlsx?__blob=publicationFile"
                     ))

    register(DataDep("COVID-19-Impfungen","""
    Autor: Robert-Koch-Institut
    Dataset: Robert Koch-Institut (2021): COVID-19-Impfungen in Deutschland, Berlin: Zenodo. DOI:10.5281/zenodo.5126652.
    Website: https://github.com/robert-koch-institut/COVID-19-Impfungen_in_Deutschland

Die Impfdaten bilden einen tagesaktuellen Stand (8:30 Uhr) aller an das RKI gemeldeten Impfungen in Deutschland ab. Im Dateinamen repräsentiert die Sequenz "JJJJ-MM-TT" das Erstellungsdatum der Datei und gleichzeitig das Datum des enthaltenen Datenstands. "JJJJ" steht dabei für das Jahr, "MM" für den Monat und "TT" für den Tag der Erstellung bzw. des enthaltenen Datenstands. Die "Aktuell_Deutschland_Bundeslaender_COVID-19-Impfungen.csv" ist identisch mit dem neusten Datenstand des Archivs.
Variablen

Die Impfdaten differenzieren nach verschiedenen Merkmalen einer Impfgruppe. Pro Eintrag bzw. Zeile ist eine eineindeutige Impfgruppe abgebildet. Eine Impfgruppe umfasst in der Regel keine Einzelfälle. Jedoch ist es möglich, dass in einer Impfgruppe nur ein Fall enthalten ist. Eine Impfgruppe wird grundlegend durch folgende Merkmale charakterisiert (in den Klammern finden sich die Variablen dieser Merkmale):

    Ort der Impfung (BundeslandId_Impfort)
    Impfung (Impfstoff, Impfserie)

Zusätzlich werden folgende Variablen angegeben:

    Datum der Impfung (Impfdatum)
    Anzahl der Impfungen in der Gruppe (Anzahl)

Eine Impfgruppe nimmt eine eineindeutige Ausprägung hinsichtlich der Anzahl der Impfungen in einem Bundesland, des Impfstoffes und der Impfserie an. Für jede Impfgruppe wird die tägliche Anzahl neuer Impfungen ausgewiesen, sofern diese größer null sind. Für jedes Datum ist angegeben, wie viele Personen, differenziert nach den oben aufgeführten Variablen, geimpft wurden.
Variablenausprägungen

Die Impfdaten enthalten die in der folgenden Tabelle abgebildeten Variablen und deren Ausprägungen:

Impfdatum JJJJ-MM-TT
Datum der Impfungen

BundeslandId_Impfort
01 bis 16 : Bundesland ID
17 : Bundesressorts
Identifikationsnummer des Bundeslandes basierend auf dem Amtlichen Gemeindeschlüssel (AGS).
Impfungen des Bundesressorts werden separat ausgewiesen, da die Impfstellen des Bundes ohne exakte Angabe des Impfortes melden

Impfstoff
AstraZeneca: AstraZeneca, Moderna: Moderna, Comirnaty: BioNTech/Pfizer, Janssen: Janssen‑Cilag/Johnson & Johnson
Verabreichter Impfstoff

Impfserie
1: Erstimpfung
2: Zweitimpfung
3: Auffrischungsimpfung
Angabe zur Erst-, Zweit- oder Auffrischungsimpfung

Anzahl 	Natürliche Zahl 	≥1
Anzahl der Impfungen in der Impfgruppe
    """,
                     "https://github.com/robert-koch-institut/COVID-19-Impfungen_in_Deutschland/raw/master/Aktuell_Deutschland_Bundeslaender_COVID-19-Impfungen.csv"
                     ))

    register(ManualDataDep("Bevölkerung Deutschland", """
    Autor: Statistisches Bundesamt
    Data: Bevölkerung: Deutschland, Stichtag, Altersjahre, Nationalität/Geschlecht/Familienstand
    Code: 12411-0006

    Website: https://www-genesis.destatis.de/genesis//online?operation=table&code=12411-0006&bypass=true&levelindex=1&levelid=1639608813890#abreadcrumb
    """))

    register(ManualDataDep("Bevölkerung Deutschland Bundesländer", """
    Autor: Statistisches Bundesamt
    Data: Bevölkerung: Bundesländer, Stichtag, Geschlecht, Altersjahre
    Code: 12411-0013

    Website: https://www-genesis.destatis.de/genesis//online?operation=table&code=12411-0013&bypass=true&levelindex=1&levelid=1639608813890#abreadcrumb
    """))
end

max_stars_repo_name: gkappler/DeutschlandDaten

max_issues_repo_name: gkappler/DeutschlandDaten

max_forks_repo_name: gkappler/DeutschlandDaten


text:
"""
DIJKSTRA
--------
    compute shortest paths using Dijkstra's algorithm.
    d = dijkstra(A,u) computes the shortest path from vertex u to all nodes
    reachable from vertex u using Dijkstra's algorithm for the problem.
    The graph is given by the weighted sparse matrix A, where A(i,j) is
    the distance between vertex i and j.  In the output vector d,
    the entry d[v] is the minimum distance between vertex u and vertex v.
    A vertex w unreachable from u has d(w)=Inf.
    pred is the predecessor tree to generate the actual shorest paths.
    In the predecessor tree pred[v] is the vertex preceeding v in the
    shortest path and pred[u]=0. Any unreachable vertex has pred[w]=0 as well.
    If your network is unweighted, then use bfs instead.

Functions
---------
- (d,pred) = dijkstra(A::MatrixNetwork,u::Int64)
- (d,pred) = dijkstra{F}(A::SparseMatrixCSC{F,Int64},u::Int64)

Example
-------
~~~
# Find the minimum travel time between Los Angeles (LAX) and
# Rochester Minnesota (RST).

(A,xy,labels) = load_matrix_network_metadata("airports")
A = -A; # fix funny encoding of airport data
lax = 247; rst = 355
(d,pred) = dijkstra(A,lax)
@printf("Minimum time: %d",d[rst]); #Print the path
@printf("Path:")
u = rst
while(u != lax)
    @printf("%s <-- ", labels[u])
    u = pred[u];
    if (u == lax)
        @printf("%s", labels[lax])
    end
end
~~~
"""

function dijkstra(A::MatrixNetwork,u::Int64)
    (rp,ci,ai) = (A.rp, A.ci, A.vals)
    return dijkstra_internal(rp,ci,ai,u)
end

function dijkstra{F}(A::SparseMatrixCSC{F,Int64},u::Int64)
    (rp,ci,ai) = sparse_to_csr(A)
    return dijkstra_internal(rp,ci,ai,u)
end

function dijkstra_internal{F}(rp::Vector{Int64},ci::Vector{Int64},ai::Vector{F},u::Int64)

    if any(ai.<0)
        error("dijkstra''s algorithm cannot handle negative edge weights")
    end

    n = length(rp) - 1
    d = Inf*ones(Float64,n)
    T = zeros(Int64,n)
    L = zeros(Int64,n)
    pred = zeros(Int64,length(rp)-1)

    n = 1
    T[n] = u
    L[u] = n # oops, n is now the size of the heap

    # enter the main dijkstra loop
    d[u] = 0
    while n > 0
        v = T[1]
        ntop = T[n]
        T[1] = ntop
        L[ntop] = 1
        n = n - 1 # pop the head off the heap
        k = 1
        kt = ntop # move element T[1] down the heap
        while true
            i = 2*k
            if i > n
                break
            end      # end of heap
            if i == n
                it = T[i]        # only one child, so skip
            else               # pick the smallest child
                lc = T[i]
                rc = T[i+1]
                it = lc
                if d[rc] < d[lc]
                    i = i+1
                    it = rc
                end # right child is smaller
            end
            if d[kt] < d[it]
                break   # at correct place, so end
            else
                T[k] = it
                L[it] = k
                T[i] = kt
                L[kt] = i
                k=i # swap
            end
        end       # end heap down
        # for each vertex adjacent to v, relax it
        for ei = rp[v]:rp[v+1]-1       # ei is the edge index
            w = ci[ei]
            ew = ai[ei]          # w is the target, ew is the edge weight
            # relax edge (v,w,ew)
            if d[w] > d[v] + ew
                d[w] = d[v] + ew
                pred[w] = v
                # check if w is in the heap
                k = L[w]
                onlyup=false
                if k == 0
                    # element not in heap, only move the element up the heap
                    n = n+1
                    T[n] = w
                    L[w] = n
                    k = n
                    kt = w
                    onlyup = true
                else
                    kt = T[k]
                end
                # update the heap, move the element down in the heap
                while !onlyup
                    i = 2 * k
                    if i > n
                        break
                    end          # end of heap
                    if i == n
                        it = T[i]    # only one child, so skip
                    else            # pick the smallest child
                        lc = T[i]
                        rc = T[i+1]
                        it = lc
                        if d[rc] < d[lc]
                            i = i+1
                            it = rc
                        end # right child is smaller
                    end
                    if d[kt] < d[it]
                        break      # at correct place, so end
                    else
                        T[k] = it
                        L[it] = k
                        T[i] = kt
                        L[kt] = i
                        k = i # swap
                    end
                end
                # move the element up the heap
                j = k
                tj = T[j]
                while j > 1                      # j==1 => element at top of heap
                    j2 = round(Int,floor(j/2))
                    tj2 = T[j2]    # parent element
                    if d[tj2] < d[tj]
                        break      # parent is smaller, so done
                    else                         # parent is larger, so swap
                        T[j2] = tj
                        L[tj] = j2
                        T[j] = tj2
                        L[tj2] = j
                        j=j2
                    end
                end
            end
        end
    end
    return (d,pred)
end


text:
function iterate!(job::BasicMCJob, ::Type{ARS}, ::Type{Univariate})
  if job.tuner.verbose
    job.sstate.tune.proposed += 1
  end

  job.sstate.pstate.value = job.pstate.value+job.sampler.jumpscale*randn()

  job.parameter.logtarget!(job.sstate.pstate)

  job.sstate.logproposal = job.sampler.logproposal(job.sstate.pstate.value)

  job.sstate.weight = job.sstate.pstate.logtarget-job.sampler.proposalscale-job.sstate.logproposal

  if (job.sstate.weight > log(rand()))
    job.pstate.value = job.sstate.pstate.value

    job.pstate.logtarget = job.sstate.pstate.logtarget
    if in(:loglikelihood, job.outopts[:monitor]) && job.parameter.loglikelihood! != nothing
      job.pstate.loglikelihood = job.sstate.pstate.loglikelihood
    end
    if in(:logprior, job.outopts[:monitor]) && job.parameter.logprior! != nothing
      job.pstate.logprior = job.sstate.pstate.logprior
    end

    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = true
      end
    end
    if job.tuner.verbose
      job.sstate.tune.accepted += 1
    end
  else
    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = false
      end
    end
  end

  if job.tuner.verbose
    if job.sstate.tune.totproposed <= job.range.burnin && mod(job.sstate.tune.proposed, job.tuner.period) == 0
      rate!(job.sstate.tune)
      println(
        "Burnin iteration ",
        job.fmt_iter(job.sstate.tune.totproposed),
        " of ",
        job.range.burnin,
        ": ",
        job.fmt_perc(100*_job.sstate.tune.rate),
        " % acceptance rate"
      )
      reset_burnin!(job.sstate.tune)
    end
  end
end

function iterate!(job::BasicMCJob, ::Type{ARS}, ::Type{Multivariate})
  if job.tuner.verbose
    job.sstate.tune.proposed += 1
  end

  job.sstate.pstate.value[:] = job.pstate.value+job.sampler.jumpscale*randn(job.pstate.size)

  job.parameter.logtarget!(job.sstate.pstate)

  job.sstate.logproposal = job.sampler.logproposal(job.sstate.pstate.value)

  job.sstate.weight = job.sstate.pstate.logtarget-job.sampler.proposalscale-job.sstate.logproposal

  if (job.sstate.weight > log(rand()))
    job.pstate.value = copy(job.sstate.pstate.value)

    job.pstate.logtarget = job.sstate.pstate.logtarget
    if in(:loglikelihood, job.outopts[:monitor]) && job.parameter.loglikelihood! != nothing
      job.pstate.loglikelihood = job.sstate.pstate.loglikelihood
    end
    if in(:logprior, job.outopts[:monitor]) && job.parameter.logprior! != nothing
      job.pstate.logprior = job.sstate.pstate.logprior
    end

    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = true
      end
    end
    if job.tuner.verbose
      job.sstate.tune.accepted += 1
    end
  else
    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = false
      end
    end
  end

  if job.tuner.verbose
    if job.sstate.tune.totproposed <= job.range.burnin && mod(job.sstate.tune.proposed, job.tuner.period) == 0
      rate!(job.sstate.tune)
      println(
        "Burnin iteration ",
        job.fmt_iter(job.sstate.tune.totproposed),
        " of ",
        job.range.burnin,
        ": ",
        job.fmt_perc(100*_job.sstate.tune.rate),
        " % acceptance rate"
      )
      reset_burnin!(job.sstate.tune)
    end
  end
end


text:
using MinimalRides
using Plots
using DataFrames
using Dates
using Geodesy
using Underscores

import MinimalRides: make_subset, results, DistanceSmoothing

########################################
# Data loading
########################################

RAW_DATA_URL_UPDATED = ("https://dl.dropboxusercontent.com" *
    "/sh/hibzl6fkzukltk9/AABSMicBJlwMnlmA3ljt1uY5a" *
    "/data_samples-json2.zip")
CACHE_DATA = "/media/win/Data/data_samples-json2.zip"

rides = MinimalRides.load(RAW_DATA_URL_UPDATED, CACHE_DATA);
rides_ecef = MinimalRides.Ride{ECEF}.(rides)

########################################
# Data sampling
########################################

DAYS = Date.(["2020-05-21"])
TYPES = [:freight]

subset = make_subset(rides, DAYS, TYPES)
subset_ecef = make_subset(rides_ecef, DAYS, TYPES)

# Sanity check
@assert map(x -> x.info.id, subset) == map(x -> x.info.id, subset_ecef)

########################################
# Main calculations
########################################

df_ecef = results(subset_ecef; smoothalg = DistanceSmoothing(0.5), search_radius = 2.5 * 1.2) |> DataFrame
pairs_df_ecef = @_ sort(df_ecef, :op, rev = true) |> filter((_.op > 0.75) & (_.len > 400), __)

for i in axes(pairs_df_ecef, 1)
    plot(map(x -> (x.coord.lat, x.coord.lon), subset[pairs_df_ecef[i, :i]].route),
        title = "Rides ($(pairs_df_ecef[i, :i]), $(pairs_df_ecef[i, :j]))",
        label = "Ride $(pairs_df_ecef[i, :i])")
    d = plot!(map(x -> (x.coord.lat, x.coord.lon), subset[pairs_df_ecef[i, :j]].route),
        label = "Ride $(pairs_df_ecef[i, :j])")
    display(d)
end

########################################
# Sanity check
########################################

df = results(subset; smoothalg = DistanceSmoothing(0.5), search_radius = 2.5 * 1.2) |> DataFrame

pairs_df = @_ sort(df, :op, rev = true) |> filter((_.op > 0.75) & (_.len > 400), __)

println("LatLon algorithm: ")
println(pairs_df)
println("ECEF algorithm: ")
println(pairs_df_ecef)

########################################
# Benchmarking
########################################

@benchmark(results($subset_ecef))
# BenchmarkTools.Trial:
#   memory estimate:  3.65 MiB
#   allocs estimate:  1721
#   --------------
#   minimum time:     41.444 ms (0.00% GC)
#   median time:      43.135 ms (0.00% GC)
#   mean time:        43.243 ms (0.00% GC)
#   maximum time:     50.773 ms (0.00% GC)
#   --------------
#   samples:          116
#   evals/sample:     1


text:
using MKL, MKL_jll, Test
import LinearAlgebra

@show LinearAlgebra.BLAS.get_config()

if VERSION > MKL.JULIA_VER_NEEDED
    @test LinearAlgebra.BLAS.get_config().loaded_libs[1].libname == libmkl_rt
else
    @test LinearAlgebra.BLAS.vendor() == :mkl
end

@test LinearAlgebra.peakflops() > 0

include(joinpath(Sys.STDLIB, "LinearAlgebra", "test", "runtests.jl"))


text:
__precompile__(true)
module SortingLab

# export sort, sort!
# package code goes here

#using InternedStrings, StatsBase, CategoricalArrays

# import InternedStrings
import StatsBase: countmap
import CategoricalArrays: CategoricalArray
import Base.Threads: @threads

export fsortperm, radixsort, radixsort!, fsort, fsort!

include("fsortperm.jl")
include("sorttwo!.jl")
include("uint_hist.jl")
include("uint_mapping.jl")
include("radixsort_string.jl")
include("fsortperm_string.jl")
include("fsortperm_Integer.jl")
include("fsort_CategoricalArrays.jl")
include("fsort.jl")
include("fsort-missing.jl")
# include("radixsort-StrFs.jl")
# include("../benchmarks/is_parallel_hist_faster_YES.jl")
end # module


text:
# phi = (galpha, alpha) -> linefunc(galpha, alpha, func, x, d, xtmp, g)

# Display flags are represented as a bitfield
# (not exported, but can use via OptimizeMod.ITER, for example)
const one64 = convert(UInt64, 1)
const FINAL       = one64
const ITER        = one64 << 1
const PARAMETERS  = one64 << 2
const GRADIENT    = one64 << 3
const SEARCHDIR   = one64 << 4
const ALPHA       = one64 << 5
const BETA        = one64 << 6
# const ALPHAGUESS  = one64 << 7 TODO: not needed
const BRACKET     = one64 << 8
const LINESEARCH  = one64 << 9
const UPDATE      = one64 << 10
const SECANT2     = one64 << 11
const BISECT      = one64 << 12
const BARRIERCOEF = one64 << 13
display_nextbit = 14

# There are some modifications and/or extensions from what's in the
# paper (these may or may not be extensions of the cg_descent code
# that can be downloaded from Hager's site; his code has undergone
# numerous revisions since publication of the paper):
#   cgdescent: the termination condition employs a "unit-correct"
#     expression rather than a condition on gradient
#     components---whether this is a good or bad idea will require
#     additional experience, but preliminary evidence seems to suggest
#     that it makes "reasonable" choices over a wider range of problem
#     types.
#   linesearch: the Wolfe conditions are checked only after alpha is
#     generated either by quadratic interpolation or secant
#     interpolation, not when alpha is generated by bisection or
#     expansion. This increases the likelihood that alpha will be a
#     good approximation of the minimum.
#   linesearch: In step I2, we multiply by psi2 only if the convexity
#     test failed, not if the function-value test failed. This
#     prevents one from going uphill further when you already know
#     you're already higher than the point at alpha=0.
#   both: checks for Inf/NaN function values
#   both: support maximum value of alpha (equivalently, c). This
#     facilitates using these routines for constrained minimization
#     when you can calculate the distance along the path to the
#     disallowed region. (When you can't easily calculate that
#     distance, it can still be handled by returning Inf/NaN for
#     exterior points. It's just more efficient if you know the
#     maximum, because you don't have to test values that won't
#     work.) The maximum should be specified as the largest value for
#     which a finite value will be returned.  See, e.g., limits_box
#     below.  The default value for alphamax is Inf. See alphamaxfunc
#     for cgdescent and alphamax for linesearch_hz.

const DEFAULTDELTA = 0.1
const DEFAULTSIGMA = 0.9

function hz_linesearch!{T}(df,
                           x::Array{T},
                           s::Array,
                           xtmp::Array,
                           g::Array,
                           lsr::LineSearchResults{T},
                           c::Real,
                           mayterminate::Bool,
                           delta::Real = DEFAULTDELTA,
                           sigma::Real = DEFAULTSIGMA,
                           alphamax::Real = convert(T,Inf),
                           rho::Real = convert(T,5),
                           epsilon::Real = convert(T,1e-6),
                           gamma::Real = convert(T,0.66),
                           linesearchmax::Integer = 50,
                           psi3::Real = convert(T,0.1),
                           iterfinitemax::Integer = ceil(Integer, -log2(eps(T))),
                           display::Integer = 0)
    if display & LINESEARCH > 0
        println("New linesearch")
    end
    f_calls = 0
    g_calls = 0

    phi0 = lsr.value[1]
    dphi0 = lsr.slope[1]
    (isfinite(phi0) && isfinite(dphi0)) || error("Initial value and slope must be finite")
    philim = phi0 + epsilon * abs(phi0)
    @assert c > 0
    @assert isfinite(c) && c <= alphamax
    phic, dphic, f_up, g_up = linefunc!(df, x, s, c, xtmp, g, true)
    f_calls += f_up
    g_calls += g_up
    iterfinite = 1
    while !(isfinite(phic) && isfinite(dphic)) && iterfinite < iterfinitemax
        mayterminate = false
        lsr.nfailures += 1
        iterfinite += 1
        c *= psi3
        phic, dphic, f_up, g_up = linefunc!(df, x, s, c, xtmp, g, true)
        f_calls += f_up
        g_calls += g_up
    end
    if !(isfinite(phic) && isfinite(dphic))
        println("Warning: failed to achieve finite new evaluation point, using alpha=0")
        return zero(T), f_calls, g_calls # phi0
    end
    push!(lsr, c, phic, dphic)
    # If c was generated by quadratic interpolation, check whether it
    # satisfies the Wolfe conditions
    if mayterminate &&
          satisfies_wolfe(c, phic, dphic, phi0, dphi0, philim, delta, sigma)
        if display & LINESEARCH > 0
            println("Wolfe condition satisfied on point alpha = ", c)
        end
        return c, f_calls, g_calls # phic
    end
    # Initial bracketing step (HZ, stages B0-B3)
    isbracketed = false
    ia = 1
    ib = 2
    @assert length(lsr) == 2
    iter = 1
    cold = -one(T)
    while !isbracketed && iter < linesearchmax
        if display & BRACKET > 0
            println("bracketing: ia = ", ia,
                    ", ib = ", ib,
                    ", c = ", c,
                    ", phic = ", phic,
                    ", dphic = ", dphic)
        end
        if dphic >= 0
            # We've reached the upward slope, so we have b; examine
            # previous values to find a
            ib = length(lsr)
            for i = (ib - 1):-1:1
                if lsr.value[i] <= philim
                    ia = i
                    break
                end
            end
            isbracketed = true
        elseif lsr.value[end] > philim
            # The value is higher, but the slope is downward, so we must
            # have crested over the peak. Use bisection.
            ib = length(lsr)
            ia = ib - 1
            if c != lsr.alpha[ib] || lsr.slope[ib] >= 0
                error("c = ", c, ", lsr = ", lsr)
            end
            # ia, ib = bisect(phi, lsr, ia, ib, philim) # TODO: Pass options
            ia, ib, f_up, g_up = bisect!(df, x, s, xtmp, g, lsr, ia, ib, philim, display)
            f_calls += f_up
            g_calls += g_up
            isbracketed = true
        else
            # We'll still going downhill, expand the interval and try again
            cold = c
            c *= rho
            if c > alphamax
                c = (alphamax + cold)/2
                if display & BRACKET > 0
                    println("bracket: exceeding alphamax, bisecting: alphamax = ", alphamax,
                            ", cold = ", cold, ", new c = ", c)
                end
                if c == cold || nextfloat(c) >= alphamax
                    return cold, f_calls, g_calls
                end
            end
            phic, dphic, f_up, g_up = linefunc!(df, x, s, c, xtmp, g, true)
            f_calls += f_up
            g_calls += g_up
            iterfinite = 1
            while !(isfinite(phic) && isfinite(dphic)) && c > nextfloat(cold) && iterfinite < iterfinitemax
                alphamax = c
                lsr.nfailures += 1
                iterfinite += 1
                if display & BRACKET > 0
                    println("bracket: non-finite value, bisection")
                end
                c = (cold + c) / 2
                phic, dphic, f_up, g_up = linefunc!(df, x, s, c, xtmp, g, true)
                f_calls += f_up
                g_calls += g_up
            end
            if !(isfinite(phic) && isfinite(dphic))
                return cold, f_calls, g_calls
            elseif dphic < 0 && c == alphamax
                # We're on the edge of the allowed region, and the
                # value is still decreasing. This can be due to
                # roundoff error in barrier penalties, a barrier
                # coefficient being so small that being eps() away
                # from it still doesn't turn the slope upward, or
                # mistakes in the user's function.
                if iterfinite >= iterfinitemax
                    println("Warning: failed to expand interval to bracket with finite values. If this happens frequently, check your function and gradient.")
                    println("c = ", c,
                            ", alphamax = ", alphamax,
                            ", phic = ", phic,
                            ", dphic = ", dphic)
                end
                return c, f_calls, g_calls
            end
            push!(lsr, c, phic, dphic)
        end
        iter += 1
    end
    while iter < linesearchmax
        a = lsr.alpha[ia]
        b = lsr.alpha[ib]
        @assert b > a
        if display & LINESEARCH > 0
            println("linesearch: ia = ", ia,
                    ", ib = ", ib,
                    ", a = ", a,
                    ", b = ", b,
                    ", phi(a) = ", lsr.value[ia],
                    ", phi(b) = ", lsr.value[ib])
        end
        if b - a <= eps(b)
            return a, f_calls, g_calls # lsr.value[ia]
        end
        iswolfe, iA, iB, f_up, g_up = secant2!(df, x, s, xtmp, g, lsr, ia, ib, philim, delta, sigma, display)
        f_calls += f_up
        g_calls += g_up
        if iswolfe
            return lsr.alpha[iA], f_calls, g_calls # lsr.value[iA]
        end
        A = lsr.alpha[iA]
        B = lsr.alpha[iB]
        @assert B > A
        if B - A < gamma * (b - a)
            if display & LINESEARCH > 0
                println("Linesearch: secant succeeded")
            end
            if nextfloat(lsr.value[ia]) >= lsr.value[ib] && nextfloat(lsr.value[iA]) >= lsr.value[iB]
                # It's so flat, secant didn't do anything useful, time to quit
                if display & LINESEARCH > 0
                    println("Linesearch: secant suggests it's flat")
                end
                return A, f_calls, g_calls
            end
            ia = iA
            ib = iB
        else
            # Secant is converging too slowly, use bisection
            if display & LINESEARCH > 0
                println("Linesearch: secant failed, using bisection")
            end
            c = (A + B) / convert(T, 2)
            # phic = phi(gphi, c) # TODO: Replace
            phic, dphic, f_up, g_up = linefunc!(df, x, s, c, xtmp, g, true)
            f_calls += f_up
            g_calls += g_up
            @assert isfinite(phic) && isfinite(dphic)
            push!(lsr, c, phic, dphic)
            # ia, ib = update(phi, lsr, iA, iB, length(lsr), philim) # TODO: Pass options
            ia, ib, f_up, g_up = update!(df, x, s, xtmp, g, lsr, iA, iB, length(lsr), philim, display)
            f_calls += f_up
            g_calls += g_up
        end
        iter += 1
    end
    error("Linesearch failed to converge")
end

# Check Wolfe & approximate Wolfe
function satisfies_wolfe{T<:Number}(c::T,
                                    phic::Real,
                                    dphic::Real,
                                    phi0::Real,
                                    dphi0::Real,
                                    philim::Real,
                                    delta::Real,
                                    sigma::Real)
    wolfe1 = delta * dphi0 >= (phic - phi0) / c &&
               dphic >= sigma * dphi0
    wolfe2 = (2.0 * delta - 1.0) * dphi0 >= dphic >= sigma * dphi0 &&
               phic <= philim
    return wolfe1 || wolfe2
end

# HZ, stages S1-S4
function secant(a::Real, b::Real, dphia::Real, dphib::Real)
    return (a * dphib - b * dphia) / (dphib - dphia)
end
function secant(lsr::LineSearchResults, ia::Integer, ib::Integer)
    return secant(lsr.alpha[ia], lsr.alpha[ib], lsr.slope[ia], lsr.slope[ib])
end
# phi
function secant2!{T}(df,
                     x::Array,
                     s::Array,
                     xtmp::Array,
                     g::Array,
                     lsr::LineSearchResults{T},
                     ia::Integer,
                     ib::Integer,
                     philim::Real,
                     delta::Real = DEFAULTDELTA,
                     sigma::Real = DEFAULTSIGMA,
                     display::Integer = 0)
    f_calls = 0
    g_calls = 0

    phi0 = lsr.value[1]
    dphi0 = lsr.slope[1]
    a = lsr.alpha[ia]
    b = lsr.alpha[ib]
    dphia = lsr.slope[ia]
    dphib = lsr.slope[ib]
    if !(dphia < 0 && dphib >= 0)
        error(string("Search direction is not a direction of descent; ",
                     "this error may indicate that user-provided derivatives are inaccurate. ",
                      @sprintf "(dphia = %f; dphib = %f)" dphia dphib))
    end
    c = secant(a, b, dphia, dphib)
    if display & SECANT2 > 0
        println("secant2: a = ", a, ", b = ", b, ", c = ", c)
    end
    @assert isfinite(c)
    # phic = phi(tmpc, c) # Replace
    phic, dphic, f_up, g_up = linefunc!(df, x, s, c, xtmp, g, true)
    f_calls += f_up
    g_calls += g_up
    @assert isfinite(phic) && isfinite(dphic)
    push!(lsr, c, phic, dphic)
    ic = length(lsr)
    if satisfies_wolfe(c, phic, dphic, phi0, dphi0, philim, delta, sigma)
        if display & SECANT2 > 0
            println("secant2: first c satisfied Wolfe conditions")
        end
        return true, ic, ic, f_calls, g_calls
    end
    # iA, iB = update(phi, lsr, ia, ib, ic, philim)
    iA, iB, f_up, g_up = update!(df, x, s, xtmp, g, lsr, ia, ib, ic, philim, display)
    f_calls += f_up
    g_calls += g_up
    if display & SECANT2 > 0
        println("secant2: iA = ", iA, ", iB = ", iB, ", ic = ", ic)
    end
    a = lsr.alpha[iA]
    b = lsr.alpha[iB]
    doupdate = false
    if iB == ic
        # we updated b, make sure we also update a
        c = secant(lsr, ib, iB)
    elseif iA == ic
        # we updated a, do it for b too
        c = secant(lsr, ia, iA)
    end
    if a <= c <= b
        if display & SECANT2 > 0
            println("secant2: second c = ", c)
        end
        # phic = phi(tmpc, c) # TODO: Replace
        phic, dphic, f_up, g_up = linefunc!(df, x, s, c, xtmp, g, true)
        f_calls += f_up
        g_calls += g_up
        @assert isfinite(phic) && isfinite(dphic)
        push!(lsr, c, phic, dphic)
        ic = length(lsr)
        # Check arguments here
        if satisfies_wolfe(c, phic, dphic, phi0, dphi0, philim, delta, sigma)
            if display & SECANT2 > 0
                println("secant2: second c satisfied Wolfe conditions")
            end
            return true, ic, ic, f_calls, g_calls
        end
        iA, iB, f_up, g_up = update!(df, x, s, xtmp, g, lsr, iA, iB, ic, philim, display)
        f_calls += f_up
        g_calls += g_up
    end
    if display & SECANT2 > 0
        println("secant2 output: a = ", lsr.alpha[iA], ", b = ", lsr.alpha[iB])
    end
    return false, iA, iB, f_calls, g_calls
end

# HZ, stages U0-U3
# Given a third point, pick the best two that retain the bracket
# around the minimum (as defined by HZ, eq. 29)
# b will be the upper bound, and a the lower bound
function update!(df,
                 x::Array,
                 s::Array,
                 xtmp::Array,
                 g::Array,
                 lsr::LineSearchResults,
                 ia::Integer,
                 ib::Integer,
                 ic::Integer,
                 philim::Real,
                 display::Integer = 0)
    a = lsr.alpha[ia]
    b = lsr.alpha[ib]
    # Debugging (HZ, eq. 4.4):
    @assert lsr.slope[ia] < 0
    @assert lsr.value[ia] <= philim
    @assert lsr.slope[ib] >= 0
    @assert b > a
    c = lsr.alpha[ic]
    phic = lsr.value[ic]
    dphic = lsr.slope[ic]
    if display & UPDATE > 0
        println("update: ia = ", ia,
                ", a = ", a,
                ", ib = ", ib,
                ", b = ", b,
                ", c = ", c,
                ", phic = ", phic,
                ", dphic = ", dphic)
    end
    if c < a || c > b
        return ia, ib, 0, 0  # it's out of the bracketing interval
    end
    if dphic >= 0
        return ia, ic, 0, 0  # replace b with a closer point
    end
    # We know dphic < 0. However, phi may not be monotonic between a
    # and c, so check that the value is also smaller than phi0.  (It's
    # more dangerous to replace a than b, since we're leaving the
    # secure environment of alpha=0; that's why we didn't check this
    # above.)
    if phic <= philim
        return ic, ib, 0, 0  # replace a
    end
    # phic is bigger than phi0, which implies that the minimum
    # lies between a and c. Find it via bisection.
    return bisect!(df, x, s, xtmp, g, lsr, ia, ic, philim, display)
end

# HZ, stage U3 (with theta=0.5)
function bisect!{T}(df,
                    x::Array,
                    s::Array,
                    xtmp::Array,
                    g::Array,
                    lsr::LineSearchResults{T},
                    ia::Integer,
                    ib::Integer,
                    philim::Real,
                    display::Integer = 0)
    f_calls = 0
    g_calls = 0
    gphi = convert(T, NaN)
    a = lsr.alpha[ia]
    b = lsr.alpha[ib]
    # Debugging (HZ, conditions shown following U3)
    @assert lsr.slope[ia] < 0
    @assert lsr.value[ia] <= philim
    @assert lsr.slope[ib] < 0       # otherwise we wouldn't be here
    @assert lsr.value[ib] > philim
    @assert b > a
    while b - a > eps(b)
        if display & BISECT > 0
            println("bisect: a = ", a, ", b = ", b, ", b - a = ", b - a)
        end
        d = (a + b) / convert(T, 2)
        phid, gphi, f_up, g_up = linefunc!(df, x, s, d, xtmp, g, true)
        f_calls += f_up
        g_calls += g_up
        @assert isfinite(phid) && isfinite(gphi)
        push!(lsr, d, phid, gphi)
        id = length(lsr)
        if gphi >= 0
            return ia, id, f_calls, g_calls # replace b, return
        end
        if phid <= philim
            a = d # replace a, but keep bisecting until dphib > 0
            ia = id
        else
            b = d
            ib = id
        end
    end
    return ia, ib, f_calls, g_calls
end

# Define one-parameter function for line searches
function linefunc!(df,
                   x::Array,
                   s::Array,
                   alpha::Real,
                   xtmp::Array,
                   g::Array,
                   calc_grad::Bool)
    f_calls = 0
    g_calls = 0
    for i = 1:length(x)
        xtmp[i] = x[i] + alpha * s[i]
    end
    gphi = convert(eltype(g), NaN)
    if calc_grad
        val = df.fg!(xtmp, g)
        f_calls += 1
        g_calls += 1
        if isfinite(val)
            gphi = vecdot(g, s)
        end
    else
        val = df.f(xtmp)
        f_calls += 1
    end
    return val, gphi, f_calls, g_calls
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

module Vis

using ..OpticSim
using ..OpticSim: euclideancontrolpoints, evalcsg, vertex, makiemesh, detector, centroid, lower, upper, intervals, α
using ..OpticSim.Geometry
using ..OpticSim.Repeat

using Unitful
using ImageView
using Images
using ColorTypes
using ColorSchemes
using StaticArrays
using LinearAlgebra
import Makie
import GeometryBasics
import Plots
import Luxor
using FileIO

include("Visualization.jl")
include("Emitters.jl")
include("VisRepeatingStructures.jl")
include("HeadEyeModel.jl")

end # module Vis
export Vis


text:
module PreallocationTools

using ForwardDiff, ArrayInterfaceCore, LabelledArrays, Adapt

struct DiffCache{T<:AbstractArray, S<:AbstractArray}
    du::T
    dual_du::S
end

function DiffCache(u::AbstractArray{T}, siz, chunk_sizes) where {T}
    x = adapt(ArrayInterfaceCore.parameterless_type(u), zeros(T, prod(chunk_sizes .+ 1)*prod(siz)))
    DiffCache(u, x)
end

"""

`dualcache(u::AbstractArray, N::Int = ForwardDiff.pickchunksize(length(u)); levels::Int = 1)`
`dualcache(u::AbstractArray; N::AbstractArray{<:Int})`

Builds a `DualCache` object that stores both a version of the cache for `u`
and for the `Dual` version of `u`, allowing use of pre-cached vectors with
forward-mode automatic differentiation. Supports nested AD via keyword `levels`
or specifying an array of chunk_sizes.

"""
dualcache(u::AbstractArray, N::Int=ForwardDiff.pickchunksize(length(u)); levels::Int = 1) = DiffCache(u, size(u), N*ones(Int, levels))
dualcache(u::AbstractArray, N::AbstractArray{<:Int}) = DiffCache(u, size(u), N)
dualcache(u::AbstractArray, ::Type{Val{N}}; levels::Int = 1) where N = dualcache(u,N;levels)
dualcache(u::AbstractArray, ::Val{N}; levels::Int = 1) where N = dualcache(u,N;levels)

"""

`get_tmp(dc::DiffCache, u)`

Returns the `Dual` or normal cache array stored in `dc` based on the type of `u`.

"""
function get_tmp(dc::DiffCache, u::T) where T<:ForwardDiff.Dual
    nelem = div(sizeof(T), sizeof(eltype(dc.dual_du)))*length(dc.du)
    if nelem > length(dc.dual_du)
        enlargedualcache!(dc, nelem)
    end
    ArrayInterfaceCore.restructure(dc.du, reinterpret(T, view(dc.dual_du, 1:nelem)))
end

function get_tmp(dc::DiffCache, u::AbstractArray{T}) where T<:ForwardDiff.Dual
    nelem = div(sizeof(T), sizeof(eltype(dc.dual_du)))*length(dc.du)
    if nelem > length(dc.dual_du)
        enlargedualcache!(dc, nelem)
    end
    ArrayInterfaceCore.restructure(dc.du, reinterpret(T, view(dc.dual_du, 1:nelem)))
end

function get_tmp(dc::DiffCache, u::LabelledArrays.LArray{T,N,D,Syms}) where {T,N,D,Syms}
    nelem = div(sizeof(T), sizeof(eltype(dc.dual_du)))*length(dc.du)
    if nelem > length(dc.dual_du)
        enlargedualcache!(dc, nelem)
    end
    _x = ArrayInterfaceCore.restructure(dc.du, reinterpret(T, view(dc.dual_du, 1:nelem)))
    LabelledArrays.LArray{T,N,D,Syms}(_x)
end

get_tmp(dc::DiffCache, u::Number) = dc.du
get_tmp(dc::DiffCache, u::AbstractArray) = dc.du

function enlargedualcache!(dc, nelem) #warning comes only once per dualcache.
    chunksize = div(nelem, length(dc.du)) - 1
    @warn "The supplied dualcache was too small and was enlarged. This incurrs allocations
    on the first call to get_tmp. If few calls to get_tmp occur and optimal performance is essential,
    consider changing 'N'/chunk size of this dualcache to $chunksize."
    resize!(dc.dual_du, nelem)
end

"""
    b = LazyBufferCache(f=identity)

A lazily allocated buffer object.  Given an array `u`, `b[u]` returns an array of the
same type and size `f(size(u))` (defaulting to the same size), which is allocated as
needed and then cached within `b` for subsequent usage.

"""
struct LazyBufferCache{F<:Function}
    bufs::Dict # a dictionary mapping types to buffers
    sizemap::F
    LazyBufferCache(f::F=identity) where {F<:Function} = new{F}(Dict(), f) # start with empty dict
end

# override the [] method
function Base.getindex(b::LazyBufferCache, u::T) where {T<:AbstractArray}
    s = b.sizemap(size(u)) # required buffer size
    buf = get!(b.bufs, (T, s)) do
        similar(u, s) # buffer to allocate if it was not found in b.bufs
    end::T # declare type since b.bufs dictionary is untyped
    return buf
end

export dualcache, get_tmp, LazyBufferCache

end


text:
module FullyThreadedTree

    using Base.Cartesian

    abstract type AbstractTree{N} end
    abstract type AbstractFace{N} end

    include("face.jl")

    export Face,
           cells,
           at_boundary,
           at_refinement,
           regular,
           active,
           level,
           collect_faces

    include("tree.jl")

    export Tree,
           initialize_tree,
           refine!,
           coarsen!,
           active,
           parent_of_active,
           initialized,
           level,
           faces,
           siblings,
           root,
           collect_cells

    include("interface.jl")

    export iterate,
           show,
           length,
           cells,
           all_faces

    include("tools.jl")

    export volume,
           centroid,
           integrate,
           cell_distance,
           cell_volume,
           area,
           volume,
           levels

    include("plotting.jl")

    export plot,
           plot!

end # module


text:
## --- Resampling.jl

    index = Array{Int64}(undef,1000)
    resampled = Array{Int64}(undef,1000)
    @test bsr!(resampled,index,1:10,0,0.5) |> unique ⊆ 1:10
    @test bsr!(resampled,index,1:10,0,fill(0.5,1000)) |> unique ⊆ 1:10
    @test bsr!(resampled,index,1:10,fill(0,10),0.5) |> unique ⊆ 1:10
    @test bsr!(resampled,index,1:10,fill(0,10),fill(0.5,1000)) |> unique ⊆ 1:10

    resampled = Array{Float64}(undef,1000)
    # Gaussian
    bsr!(resampled, index, ones(10), 0.5, 0.5)
    @test isapprox(mean(resampled), 1, atol=0.3)
    bsr!(resampled, index, ones(10), 0.5, fill(0.5,1000))
    @test isapprox(mean(resampled), 1, atol=0.3)
    bsr!(resampled, index, ones(10), fill(0.5,10), 0.5)
    @test isapprox(mean(resampled), 1, atol=0.3)
    bsr!(resampled, index, ones(10), fill(0.5,10), fill(0.5,1000))
    @test isapprox(mean(resampled), 1, atol=0.3)
    # Other distributions
    bsr!(uniform, resampled, index, ones(10), 0.5, 0.5)
    @test isapprox(mean(resampled), 1, atol=0.3)
    bsr!(triangular, resampled, index, ones(10), 0.5, fill(0.5,1000))
    @test isapprox(mean(resampled), 1, atol=0.3)
    bsr!(triangular, resampled, index, ones(10), fill(0.5,10), 0.5)
    @test isapprox(mean(resampled), 1, atol=0.3)
    bsr!(uniform, resampled, index, ones(10), fill(0.5,10), fill(0.5,1000))
    @test isapprox(mean(resampled), 1, atol=0.3)

    @test bsresample(1:10,fill(0.5,10),1000,0.5)::Array{Float64} |> length == 1000
    @test bsresample(1:10,fill(0.5,10),1000,fill(0.5,10))::Array{Float64} |> length == 1000

    d = Dict{String,Vector{Float64}}()
    d["x"] = 1:10;  d["y"] = 2:2:20
    d["x_sigma"] = d["y_sigma"] = fill(0.5,10)
    dbs = bsresample(d, 1000, ["x","y"], 0.5)
    @test isapprox(mean(dbs["x"]), 5.5, atol=0.5)
    @test isapprox(std(dbs["x"]), 3.03, atol=0.5)
    @test isapprox(mean(dbs["y"]), 11, atol=1)
    @test isapprox(std(dbs["y"]), 6.06, atol=1)

    @test bincounts(1:100, 0, 100, 10) == (5:10:95, fill(10,10))
    @test binmeans(1:100, 1:100, 0, 100, 10) == (5:10:95, 5.5:10:95.5, fill(0.9574271077563381,10))
    @test binmeans(1:100, 1:100, 0, 100, 10, ones(100)) == (5:10:95, 5.5:10:95.5, fill(0.9574271077563381,10))
    @test binmedians(1:100,1:100,0,100,10) == (5:10:95, 5.5:10:95.5, fill(1.1720982147414096,10))

    @test randsample(1:10., 1000, rand(1000))::Array{Float64} |> length == 1000
    @test randsample(1:10, 1000)::Array{Int64} |> length == 1000
    @test unique(randsample(1:10, 1000)) ⊆ 1:10

    dr = randsample(d, 1000, ["x","y"], 0.5)
    @test isapprox(mean(dr["x"]), 5.5, atol=0.5)
    @test isapprox(std(dr["x"]), 3.03, atol=0.5)
    @test unique(dr["x"]) ⊆ 1:10
    @test isapprox(mean(dr["y"]), 11, atol=1)
    @test isapprox(std(dr["y"]), 6.06, atol=1)
    @test unique(dr["y"]) ⊆ 2:2:20

## --- Invweight

    @test invweight(0:10, 0:10, 0:10) ≈ [13.092772378121769, 13.759663290331229, 14.079390874013654, 14.244556812410089, 14.327747696132253, 14.354508911206949, 14.331218676773712, 14.251150311763046, 14.088257739618454, 13.76917452212827, 13.101581593462868]
    @test invweight_location(0:10, 0:10) ≈ [2.3478642777118957, 2.950680392065609, 3.2200556525889006, 3.348975353894235, 3.4103053158720016, 3.4297605864726126, 3.413776296513463, 3.3555688532471923, 3.2289225181937002, 2.9601916238626513, 2.3566734930529947]
    @test invweight_age(0:10) ≈ [10.744908100409873, 10.808982898265619, 10.859335221424754, 10.895581458515855, 10.91744238026025, 10.924748324734335, 10.91744238026025, 10.895581458515855, 10.859335221424754, 10.808982898265619, 10.744908100409873]

## --- bin_bsr

    x = 0:100; y = 0:100
    xmin = 0; xmax = 100; nbins = 5
    (c,m,e) = bin_bsr(x, y, xmin, xmax, nbins, x_sigma=ones(101))
    @test c == 10.0:20.0:90.0
    @test isapprox(m, [10.04, 29.94, 49.94, 69.92, 89.83], atol=0.4)
    @test isapprox(e, [1.17, 1.21, 1.23, 1.26, 1.28], atol=0.4)

    # Upper and lower CIs
    (c,m,el,eu) = bin_bsr(nanbinmean!, x, y, xmin, xmax, nbins, x_sigma=ones(101))
    @test c == 10.0:20.0:90.0
    @test isapprox(m, [10.04, 29.94, 49.94, 69.92, 89.83], atol=0.4)
    @test isapprox(el, [2.29, 2.38, 2.41, 2.49, 2.51], atol=0.5)
    @test isapprox(eu, [2.3, 2.37, 2.42, 2.51, 2.51], atol=0.5)

    # Medians, upper and lower CIs
    (c,m,el,eu) = bin_bsr(nanbinmedian!, x, y, xmin, xmax, nbins, x_sigma=ones(101))
    @test c == 10.0:20.0:90.0
    @test isapprox(m, [10.01, 29.91, 49.9, 69.88, 89.79], atol=1)
    @test isapprox(el, [4.01, 3.91, 3.9, 3.88, 3.79], atol=2)
    @test isapprox(eu, [3.99, 4.09, 4.1, 4.12, 4.21], atol=2)

    # with weights
    w = ones(101)
    (c,m,e) = bin_bsr(x, y, xmin, xmax, nbins, w, x_sigma=ones(101))
    @test c == 10.0:20.0:90.0
    @test isapprox(m, [10.04, 29.94, 49.94, 69.92, 89.83], atol=0.4)
    @test isapprox(e, [1.17, 1.21, 1.23, 1.26, 1.28], atol=0.4)

    # with 2-D array (matrix) of y data
    y = repeat(0:100, 1, 4)
    y_sigma = ones(101,4)
    (c,m,e) = bin_bsr(x, y, xmin, xmax, nbins, x_sigma=ones(101), y_sigma=y_sigma)
    @test c == 10.0:20.0:90.0
    @test isapprox(m, repeat([10.04, 29.94, 49.94, 69.92, 89.83], 1, 4), atol=0.5)
    @test isapprox(e, repeat([1.17, 1.21, 1.23, 1.26, 1.28], 1, 4), atol=0.5)

## -- bin_bsr_ratios

    x = 0:100; num = 0:100; denom=reverse(num)
    xmin = 0; xmax = 100; nbins = 5
    (c,m,el,eu) = bin_bsr_ratios(x, num, denom, xmin, xmax, nbins, x_sigma=ones(101))
    @test c == 10.0:20.0:90.0
    @test isapprox(m, [0.11, 0.43, 1.0, 2.33, 8.99], atol=0.4)
    @test isapprox(el, [0.03, 0.05, 0.09, 0.26, 2.11], atol=0.5)
    @test isapprox(eu, [0.03, 0.05, 0.1, 0.29, 3.03], atol=0.5)

    # With weights
    (c,m,el,eu) = bin_bsr_ratios(x, num, denom, xmin, xmax, nbins, ones(101), x_sigma=ones(101))
    @test c == 10.0:20.0:90.0
    @test isapprox(m, [0.11, 0.43, 1.0, 2.33, 8.99], atol=0.4)
    @test isapprox(el, [0.03, 0.05, 0.09, 0.26, 2.11], atol=0.5)
    @test isapprox(eu, [0.03, 0.05, 0.1, 0.29, 3.03], atol=0.5)

## --- Monte Carlo interpolation/fitting

    (c,m) = mcfit(0:11, ones(12), 0:11, ones(12), 1, 10, 10)
    @test c == 1:10
    @test isapprox(m, [1.15, 2.02, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 8.98, 9.85], atol = 0.25)

## --- Downsampling

    @test downsample(1:100,10) == 10:10:100
    @test downsample(reshape(1:100,10,10),2) == [12 32 52 72 92; 14 34 54 74 94; 16 36 56 76 96; 18 38 58 78 98; 20 40 60 80 100]

## ---


text:
abstract Quantum2dLattice{T} <: QuantumLattice{T}

function getSpaceDimension(lattice::Quantum2dLattice)
	return 2
end

function initializeTensorT2dQ{T}(spinmodel::QuantumSpinModel{T}; tilde::Bool = false)
	if tilde
		measureOperator = getMeasureOperator(spinmodel)
	else
		measureOperator = ones(T,getStates(spinmodel))
	end
	factorW = getFactorW(spinmodel)
	factorWp = getFactorWp(spinmodel)
	leng = size(factorW,1)==size(factorW,2) ? size(factorW,1) : throw(DomainError())
	tensorT = zeros(T,leng,leng,leng,leng,leng,leng)
	for a = 1:leng, l = 1:leng, r = 1:leng, f=1:leng, b=1:leng, u=1:leng, v=1:leng
		tensorT[l,r,f,b,u,v] += measureOperator[a] *
		factorW[a,l] * factorW[a,r] * factorW[a,f] * factorW[a,b] *
		factorWp[a,u] * factorWp[a,v]
	end
	return tensorT
end


text:
"""
    Conv(kernel, in=>out, σ; stride = 1, pading = 0, set_w = "Xavier")
This is the traditional convolution layer. `kernel` is a tuple of integers that specifies the kernel size, it must have one or two elements. And, `in` and `out` specifies number of input and out channels.

The input data must have a dimensions WHCB(weight, width, channel, batch). If you want to use a data which has a dimentions WHC, you must be add a dimentioins of B.

`stride` and `padding` are single integers or tuple(stride is tuple of 2 elements, padding is tuple of 2 elements), and if you specifies `KeepSize` to padding, we adjust sizes of input and return a matrix which has the same sizes.
`set_w` is `Xavier` or `He`, it decide a method to create a first parameter. This parameter is the same as `Dense()`.

# Example
```
julia> C = Conv((2, 2), 2=>2, relu)
Convolution(k:(2, 2), IO:2 => 2, σ:relu)

julia> C(rand(10, 10, 2, 5)) |> size
(9, 9, 2, 5)
```

!!! warning
    When you specidies `same` to `padding`, in some cases, it will be returned one size smaller.
    Because of its expression.
```
julia> C = Conv((2, 2), 2=>2, relu, padding = KeepSize)
Convolution(k:(2, 2), IO:2 => 2, σ:relu

julia> C(rand(10, 10, 2, 5)) |> size
(9, 9, 2, 5)
```
"""
struct Conv{N, M, K, I, O} <: Param
    σ
    weight
    bias
    stride::NTuple{N, Int}
    padding::NTuple{M, Int}
    dilation::NTuple{N, Int}
end

function Conv(weight::AbstractArray{T, N}, σ, k::Tuple, io::Pair{<:Integer, <:Integer}; stride = 1, padding = 0, dilation=1) where {N, T}
    stride = expand(Val(N-2), stride)
    dilation = expand(Val(N-2), dilation)
    padding = expand_padding(padding, size(weight)[1:N-2], dilation, stride)
    bias = create_bias(weight, size(weight, N))
    return Conv{length(stride), length(padding), k, io...}(σ, weight, bias, stride, padding, dilation)
end

function Conv(k::NTuple{N, Int}, io::Pair{<:Integer, <:Integer}, σ; stride = 1, padding = 0, dilation=1, set_w = "Xavier") where {N}
    weight = conv_w(k, io, set_w)
    return Conv(weight, σ, k, io, stride = stride, padding = padding, dilation = dilation)
end

function Base.show(io::IO, C::Conv{N, M, K, I, O}) where {N, M, K, I, O}
    σ = C.σ
    print("Convolution(k:$K, IO:$I => $O, σ:$σ)")
end

function (C::Conv)(x)
    σ, b, w = C.σ, C.bias, C.weight
    dim = DenseConvDims(x, w, stride = C.stride, padding = C.padding, dilation = C.dilation)
    return σ.(conv(x, w, dim).+b)
end

trainable(C::Conv) = C.weight, C.bias


text:
using Base.Test
require("SDE")
#is called with arg "travis" from travis
test_type = length(ARGS) == 1 ? ARGS[1] : ""

include("testsde.jl")
include("testlinproc.jl")

if test_type == "all"
	include("testtc.jl")
end

include("testsch.jl")
include("testnpbayes.jl")


text:
"""
ITERATIVE HARD THRESHOLDING

    iht!(v, x, y, k [, iter=1, nstep=50]) -> (Float, Int)

This function computes a hard threshold update

    β = P_Sk(β0 + μ*x'*(y - x*β0))

where `μ` is the step size (or learning rate) and P_Sk denotes the projection onto the set S_k defined by

S_k = { x in R^p : || x ||_0 <= k }.

The projection in question preserves the largest `k` components of `β` in magnitude, and it sends the remaining
`p - k` components to zero. This update is intimately related to a projected gradient step used in Landweber iteration.
Unlike the Landweber method, this function performs a line search on `μ` whenever the step size exceeds a specified
threshold `ω` given by

    ω = sum(abs2, β - β0) / sum(abs2, x*(β - β0))

By backtracking on `μ`, this function guarantees a stable estimation of a sparse `β`.
This function is based on the [HardLab](http://www.personal.soton.ac.uk/tb1m08/sparsify/sparsify.html/) demonstration code written in MATLAB by Thomas Blumensath.

Arguments:

- `v` is the `IHTVariables` object of temporary arrays
- `x` is the `n` x `p` design matrix.
- `y` is the vector of `n` responses.
- `k` is the model size.
- `iter` is the current iteration count in the IHT algorithm. Defaults to `1`.
- `nstep` is the maximum permissible number of backtracking steps. Defaults to `50`.

Output:

- `μ` is the step size used to update `β`, after backtracking.`
- `μ_step` is the number of backtracking steps used on `μ`.
"""
function iht!{T <: Float, V <: DenseVector}(
    v     :: IHTVariables{T, V},
    x     :: DenseMatrix{T},
    y     :: V,
    k     :: Int,
    iter  :: Int = 1,
    nstep :: Int = 50
)

    # compute indices of nonzeroes in beta
    _iht_indices(v, k)

    # store relevant columns of x
    # need to do this on 1st iteration
    # afterwards, only do if support changes
    if !isequal(v.idx, v.idx0) || iter < 2
        #update_xk!(v.xk, x, v.idx)   # xk = x[:,v.idx]
        copy!(v.xk, view(x, :, v.idx))
    end

    # store relevant components of gradient
    fill_perm!(v.gk, v.df, v.idx)    # gk = g[v.idx]

    # now compute subset of x*g
    A_mul_B!(v.xgk, v.xk, v.gk)

    # warn if xgk only contains zeros
    all(v.xgk .== zero(T)) && warn("Entire active set has values equal to 0")

    # compute step size
#    mu = _iht_stepsize(v, k) :: T # does not yield conformable arrays?!
    μ = (sum(abs2, v.gk) / sum(abs2, v.xgk)) :: T

    # check for finite stepsize
    isfinite(μ) || throw(error("Step size is not finite, is active set all zero?"))

    # compute gradient step
    _iht_gradstep(v, μ, k)

    # update xb
    update_xb!(v.xb, x, v.b, v.idx, k)
    #A_mul_B!(v.xb, view(x, :, v.idx), view(v.b, v.idx) )

    # calculate omega
    ω_top, ω_bot = _iht_omega(v)

    # backtrack until mu < omega and until support stabilizes
    μ_step = 0
    while _iht_backtrack(v, ω_top, ω_bot, μ, μ_step, nstep)

        # stephalving
        μ /= 2

        # recompute gradient step
        copy!(v.b,v.b0)
        _iht_gradstep(v, μ, k)

        # recompute xb
        update_xb!(v.xb, x, v.b, v.idx, k)
        #A_mul_B!(v.xb, view(x, :, v.idx), view(v.b, v.idx) )

        # calculate omega
        ω_top, ω_bot = _iht_omega(v)

        # increment the counter
        μ_step += 1
    end

    return μ::T, μ_step::Int
end

"""
L0 PENALIZED LEAST SQUARES REGRESSION

    L0_reg(x,y,k) -> IHTResults

This routine minimizes the loss function

    0.5*sum(abs2, y - x*β )

subject to `β` lying in the set S_k = { x in R^p : || x ||_0 <= k }.

It uses Thomas Blumensath's iterative hard thresholding framework to keep `β` feasible.

Arguments:

- `x` is the `n` x `p` data matrix.
- `y` is the `n`-dimensional continuous response vector.
- `k` is the desired model size (support).

Optional Arguments:

- `max_iter` is the maximum number of iterations for the algorithm. Defaults to `1000`.
- `max_step` is the maximum number of backtracking steps for the step size calculation. Defaults to `50`.
- `tol` is the global tolerance. Defaults to `1e-4`.
- `quiet` is a `Bool` that controls algorithm output. Defaults to `true` (no output).
- `v` is an `IHTVariables` structure used to house temporary arrays. Used primarily in `iht_path` for memory efficiency.

Outputs are wrapped into an `IHTResults` structure with the following fields:

- 'time' is the compute time for the algorithm. Note that this does not account for time spent initializing optional argument defaults.
- 'iter' is the number of iterations that the algorithm took before converging.
- 'loss' is the optimal loss (half of residual sum of squares) at convergence.
- 'beta' is the final estimate of `b`.
"""
function L0_reg{T <: Float, V <: DenseVector}(
    x        :: DenseMatrix{T},
    y        :: V,
    k        :: Int;
    v        :: IHTVariables{T, V} = IHTVariables(x, y, k),
    tol      :: T    = convert(T, 1e-4),
    max_iter :: Int  = 100,
    max_step :: Int  = 50,
    quiet    :: Bool = true
)

    # start timer
    tic()

    # first handle errors
    k        >= 0      || throw(ArgumentError("Value of k must be nonnegative!\n"))
    max_iter >= 0      || throw(ArgumentError("Value of max_iter must be nonnegative!\n"))
    max_step >= 0      || throw(ArgumentError("Value of max_step must be nonnegative!\n"))
    tol      >  eps(T) || throw(ArgumentError("Value of global tol must exceed machine precision!\n"))

    # initialize return values
    iter = 0                      # number of iterations of L0_reg
    exec_time = zero(T)           # compute time *within* L0_reg
    next_obj  = zero(T)           # objective value
    next_loss = zero(T)           # loss function value

    # initialize floats
    current_obj = convert(T, Inf) # tracks previous objective function value
    the_norm    = zero(T)         # norm(b - b0)
    scaled_norm = zero(T)         # the_norm / (norm(b0) + 1)
    μ           = zero(T)          # Landweber step size, 0 < tau < 2/rho_max^2

    # initialize integers
    #i      = 0                    # used for iterations in loops
    μ_step = 0                    # counts number of backtracking steps for mu

    # initialize booleans
    converged = false             # scaled_norm < tol?

    # update xb, r, gradient
    initialize_xb_r_grad!(v, x, y, k)

    # update loss and objective
    next_loss = convert(T, Inf)

    # formatted output to monitor algorithm progress
    !quiet && print_header()

    # main loop
    for iter = 1:max_iter

        # notify and break if maximum iterations are reached.
        if iter >= max_iter

            # alert about hitting maximum iterations
            !quiet && print_maxiter(max_iter, loss)

            # send elements below tol to zero
            threshold!(v.b, tol)

            # stop timer
            exec_time = toq()

            # these are output variables for function
            # wrap them into a Dict and return
            return IHTResults(exec_time, next_loss, iter, copy(v.b))
        end

        # save values from previous iterate
        copy!(v.b0, v.b)             # b0 = b
        copy!(v.xb0, v.xb)           # Xb0 = Xb
        current_obj = next_obj

        # now perform IHT step
        (μ, μ_step) = iht!(v, x, y, k, max_step, iter)

        # the IHT kernel gives us an updated x*b
        # use it to recompute residuals and gradient
        update_r_grad!(v, x, y)

        # update loss, objective, and gradient
        next_loss = sum(abs2, v.r) / 2

        # guard against numerical instabilities
        check_finiteness(next_loss)

        # track convergence
        the_norm    = chebyshev(v.b, v.b0)
        scaled_norm = (the_norm / ( norm(v.b0,Inf) + 1)) :: T
        converged   = (scaled_norm < tol) :: Bool

        # output algorithm progress
        quiet || @printf("%d\t%d\t%3.7f\t%3.7f\t%3.7f\n", iter, μ_step, μ, the_norm, next_loss)

        # check for convergence
        # if converged and in feasible set, then algorithm converged before maximum iteration
        # perform final computations and output return variables
        if converged

            # send elements below tol to zero
            threshold!(v.b, tol)

            # stop time
            exec_time = toq()

            # announce convergence
            !quiet && print_convergence(iter, next_loss, exec_time)

            # these are output variables for function
            # wrap them into a Dict and return
            return IHTResults(exec_time, next_loss, iter, copy(v.b))
        end

        # algorithm is unconverged at this point.
        # if algorithm is in feasible set, then rho should not be changing
        # check descent property in that case
        # if rho is not changing but objective increases, then abort
        if next_obj > current_obj + tol
            quiet || print_descent_error(iter, loss, next_loss)
            throw(ErrorException("Descent failure!"))
            #return IHTResults(-one(T), convert(T, -Inf), -1, [convert(T, -Inf)])
        end
    end # end main loop

    # return null result
    return IHTResults(zero(T), zero(T), 0, [zero(T)])
end # end function

"""
COMPUTE AN IHT REGULARIZATION PATH FOR LEAST SQUARES REGRESSION

    iht_path(x,y,path) -> SparseCSCMatrix

This subroutine computes best linear models for matrix `x` and response `y` by calling `L0_reg` for each model over a regularization path denoted by `path`.

Arguments:

- `x` is the `n` x `p` design matrix.
- `y` is the `n`-vector of responses.
- `path` is an `Int` vector that contains the model sizes to test.

Optional Arguments:

- `tol` is the global convergence tolerance for `L0_reg`. Defaults to `1e-4`.
- `max_iter` caps the number of iterations for the algorithm. Defaults to `1000`.
- `max_step` caps the number of backtracking steps in the IHT kernel. Defaults to `50`.
- `quiet` is a Boolean that controls the output. Defaults to `true` (no output).

Output:

- a sparse `p` x `length(path)` matrix where each column contains the computed model for each component of `path`.
"""
function iht_path{T <: Float}(
    x        :: DenseMatrix{T},
    y        :: DenseVector{T},
    path     :: DenseVector{Int};
    tol      :: T    = convert(T, 1e-4),
    max_iter :: Int  = 100,
    max_step :: Int  = 50,
    quiet    :: Bool = true
)

    # size of problem?
    (n,p) = size(x)

    # how many models will we compute?
    num_models = length(path)
    #nk = sum(path) # total number of elements in models

    # preallocate space for intermediate steps of algorithm calculations
    v  = IHTVariables(x, y, 1)
    betas = spzeros(T, p, num_models)  # a matrix to store calculated models
    #I = zeros(Int, nk)
    #J = zeros(Int, nk)
    #K = zeros(T,   nk)
    #q0 = 0

    # compute the path
    for i = 1:num_models

        # model size?
        q = path[i]

        # monitor progress
        quiet || print_with_color(:blue, "Computing model size $q.\n\n")

        # store projection of beta onto largest k nonzeroes in magnitude
        project_k!(v.b, q)

        # these arrays change in size from iteration to iteration
        # we must allocate them for every new model size
        update_variables!(v, x, q)

        # now compute current model
        output = L0_reg(x, y, q, v=v, tol=tol, max_iter=max_iter, max_step=max_step, quiet=quiet)

        # extract and save model
        # strangely more efficient to continuously reallocate sparse matrix
        # than to form it directly
        #I[(q0+1):(q0+q)] = find(v.idx)
        #J[(q0+1):(q0+q)] = i
        #K[(q0+1):(q0+q)] = v.b[v.idx]
        #q0 += q
        betas[:,i] = sparsevec(output.beta)
    end

    # return a sparsified copy of the models
    #betas = sparse(I, J, K, p, num_models, false)
    return betas
end


text:
#export impute_means, impute_zeros,
#      #Closed-form factorizations with regularization
#      matrixRegFact, quadgraphRegFact

# Subtracts each column's mean (if center=true),
# divides by each column's standard deviation (if scale=true).
# Returns (scaledData, mean, std), where mean or std may be
#  empty matrices if center or scale are false
function standardize!{T}(X::Matrix{T} ; center=true, scale=true)
    n = size(X,1)

    local m
    if center
        m = mean(X,1)
        X = X .- m
    else
        m = similar(X, 0, 0)
    end

    local s
    if scale
        s = std(X,1)
        #Make sure there are no NaN elements or 0 elements of the standard deviation
        for i in 1:length(s)
          if s[i] == 0 || s[i] === NaN
            s[i] = 1
          end
        end
        X = X ./ s
    else
        s = similar(X, 0, 0)
    end

    X, m, s
end

##Subtracts each column mean (if center=true)
##Divides each column standard deviation (if scale=true)
##Skips NA entries in DataArrays
##Returns (scaledData, mean, std)
function standardize!{T}(X::DataMatrix{T}; center=true, scale=true)
  n = size(X,1)

  local m
  if center
    m = mean(X,1,skipna=true)
    X = X .- m
  else
    m = similar(X,0,0)
  end

  local s
  if scale
    s = sqrt(var(X,1,skipna=true))
    #Make sure there are no NaN elements or 0 elements of the standard deviation
    for i in 1:length(s)
      if s[i] == 0 || s[i] === NaN
        s[i] = 1
      end
    end
    X = X ./ s
  else
    s = similar(X,0,0)
  end

  X, m, s
end

standardize(X; center=true, scale=true) = standardize!(copy(X), center=center, scale=scale)

##Imputes the column means to a DataMatrix and returns an ordinary Matrix
##of the same size.
function impute_means{T}(X::DataMatrix{T})
  m, n = size(X)
  colmeans = mean(X,1,skipna=true)
  Xmat = zeros(m,n)
  for j in 1:n
    #For entirely NA columns, just use zero
    if all(isna(X[:,j]))
      Xmat[:,j] = 0
    else
      #Impute either the mean or the existing value
      for i in 1:m
        if X[i,j] === NA
          Xmat[i,j] = colmeans[j]
        else
          Xmat[i,j] = X[i,j]
        end
      end
    end
  end
  Xmat
end

##Impute zeros on every NA value of a DataMatrix
function impute_zeros{T}(X::DataMatrix{T})
  m, n = size(X)
  Xmat = zeros(m,n)
  for j in 1:n
    for i in 1:m
      if !(X[i,j] === NA)
        Xmat[i,j] = X[i,j]
      end
    end
  end
  Xmat
end

##Closed-form solution to a quadratic loss function with graph+quadratic
##regularization on both X and Y. Performance may be poor for large number of rows
function matrixRegFact(A::AbstractMatrix, k::Int,
  α::Number, K::AbstractMatrix,
  β::Number, L::AbstractMatrix)

  n, d = size(A)
  @assert (n,n) == size(K)
  @assert (d,d) == size(L)
  iKA = (I + α*K) \ A # (I + αK)⁻¹ * A
  iKAt = iKA' # A'*(I + αK)⁻¹
  fst = -2*A'*iKA
  snd = iKAt*iKA
  thd = α*iKAt*K*iKA
  res = β*L
  for i in 1:length(res)
    res[i] += fst[i] + snd[i] + thd[i] #Add the results in a loop to save space
  end
  #=
  invIαK = inv(I + α*K)
  #Three terms to find the eigenvectors of, save time by adding in a loop
  fst = -2*A'*invIαK*A
  snd = A'*invIαK^2*A
  thd = α*A'*invIαK*K*invIαK*A
  res = β*L
  for i in 1:length(res)
    res[i] += fst[i] + snd[i] + thd[i]
  end
  =#
  Yt = eigvecs(res)[:,1:k]
  X = invIαK*A*Yt
  X,Yt'
end

##Closed-form solution to quadratic loss function with quadratic regularization
##on X and a graph+quadratic regularizer on Y.
function quadgraphRegFact(A::Matrix{Float64}, k::Int, α::Number,
  β::Number, L::AbstractMatrix)

  n, d = size(A)
  @assert (d,d) == size(L)

  ALmat = -1/(1+α)*A'*A + β*L
  Yt = eigvecs(ALmat)[:,1:k]
  X = 1/(1 + α)*A*Yt
  X,Yt'
end

"""
Initialize the GLRM with an SVD decomposition for faster convergence.
This also solves the problem where identical quadratic regularization is on
both X and Y
"""
function init_qqreg!(gg::GGLRM)
  A = impute_zeros(gg.A)
  F = svdfact(A)
  A_mul_B!(gg.X, F[:U][:,1:gg.k], Diagonal(F[:S][1:gg.k]))
  A_mul_Bt!(gg.Y, Diagonal(F[:S][1:gg.k]), F[:Vt][:,1:gg.k])
end


text:
"""
Adapt computing architecture

"""
function device(isgpu, args...)
    if isgpu
        return gpu(args)
    else
        return cpu(args)
    end
end


text:
mutable struct AlgAss{T} <: Ring
  base_ring::Ring
  mult_table::Array{T, 3} # e_i*e_j = sum_k mult_table[i, j, k]*e_k
  one::Array{T, 1}
  iscommutative::Int       # 0: don't know
                           # 1 known to be commutative
                           # 2 known to be not commutative

  function AlgAss{T}(R::Ring) where {T}
    A = new{T}()
    A.base_ring = R
    A.iscommutative = 0
    return A
  end

  function AlgAss{T}(R::Ring, mult_table::Array{T, 3}, one::Array{T, 1}) where {T}
    A = new{T}()
    A.base_ring = R
    A.mult_table = mult_table
    A.one = one
    A.iscommutative = 0
    return A
  end

  function AlgAss{T}(R::Ring, mult_table::Array{T, 3}) where {T}
    A = new{T}()
    A.base_ring = R
    A.mult_table = mult_table
    A.iscommutative = 0
    return A
  end
end

mutable struct AlgAssElem{T} <: RingElem
  parent::AlgAss{T}
  coeffs::Array{T, 1}

  function AlgAssElem{T}(A::AlgAss{T}) where {T}
    z = new{T}()
    z.parent = A
    z.coeffs = Array{T, 1}(size(A.mult_table, 1))
    for i = 1:length(z.coeffs)
      z.coeffs[i] = A.base_ring()
    end
    return z
  end

  # This does not make a copy of coeffs
  function AlgAssElem{T}(A::AlgAss{T}, coeffs::Array{T, 1}) where{T}
    z = new{T}()
    z.parent = A
    z.coeffs = coeffs
    return z
  end
end


text:
#constants

#Atmosphere
ρO0=1029; # kg/m^3ρO
cp=4184; #J/K/kg
cpa=1005; #J/K/kg
sb = 5.670e-8;
d2k=273.15;
Av=2500000;
gamma_blk=0.01;
cvapor_fac=640380;
cvapor_exp=5107.4;
cvapor_exp_ice=5897.8;
saltsat=0.980;
p0=101300; #reference pressure
gravity_mks=9.81;
humid_fac=0.606;
tref=25; #reference temperature
sref=25; #reference salinity
alpha=2e-4; #expansion coefficients for temperature
beta=7.4e-4; #expansion coefficients for salt
emissivity=0.97; # ocean emissivity
kappa=1e-5; # m^2/kg air absorption coefficient
Rgas=287.05;
karman=0.41;
albedo=0.04; #ocean albedo
SunConstant=1362
Rv=461 #gas constant for water vapor J kg-1 K-1
humid_fac=0.606;

#Land
using Parameters
@with_kw struct Const{R}
  Cs::R = 1.92e6; #soild heat capacity [J m^-3 K^-1]
  Cw::R = 4.18e6; # water heat capacity [J m^-3 K^-1]
  Cv::R = 1.8e6; # vapor heat capacity [J m^-3 K^-1]
  Ca::R = 1200; # heat capacity of air [J m^-3 K^-1]
  k::R = 0.41; # von-Karman constant [-]
  sigma::R = 5.6704e-8; # Stephan-Boltzman constant [W m^-2]
  M::R = 0.018015; #water molecular weight [kg mol^-1]
  g::R = 9.81; # gravity accelaration [m sec^-2]
  Rg::R = 8.314; # gas constant [J mol^-1 K^-1]
  gamma_o::R = 0.07189; # surface water tension [Kg sec^-2]
end

global constants=Const{Float64}()


text:
module Tree

import Base: iterate, show

export Node,
    addchild,
    addsibling,
    isroot,
    isleaf,
    lastsibling

# "Left-child, right-sibling" tree representation:
# Suppose a particular node, "A", has 3 children, "a", "b", and "c".
#     "a", "b", and "c" link to "A" as their parent.
#     "A" links "a" as its child; "a" links "b" as its sibling, and "b" links "c" as its sibling.
#     Any missing links (e.g., "c" does not have a sibling) link back to itself (c.sibling == c)
# With this organization, no arrays need to be allocated.

mutable struct Node{T}
    data::T
    parent::Node{T}
    child::Node{T}
    sibling::Node{T}

    # Constructor for the root of the tree
    function Node{T}(data) where T
        n = new{T}(data)
        n.parent = n
        n.child = n
        n.sibling = n
        n
    end
    # Constructor for all others
    function Node{T}(data, parent::Node) where T
        n = new{T}(data, parent)
        n.child = n
        n.sibling = n
        n
    end
end
Node(data::T) where {T} = Node{T}(data)
Node(data, parent::Node{T}) where {T} = Node{T}(data, parent)

function lastsibling(sib::Node)
    newsib = sib.sibling
    while sib != newsib
        sib = newsib
        newsib = sib.sibling
    end
    sib
end

function addsibling(oldersib::Node{T}, data) where T
    if oldersib.sibling != oldersib
        error("Truncation of sibling list")
    end
    youngersib = Node(data, oldersib.parent)
    oldersib.sibling = youngersib
    youngersib
end

function addchild(parent::Node{T}, data) where T
    newc = Node(data, parent)
    prevc = parent.child
    if prevc == parent
        parent.child = newc
    else
        prevc = lastsibling(prevc)
        prevc.sibling = newc
    end
    newc
end

isroot(n::Node) = n == n.parent
isleaf(n::Node) = n == n.child

show(io::IO, n::Node) = print(io, n.data)

# Iteration over children
# for c in parent
#     # do something
# end
function iterate(n::Node, state = n.child)
    n == state && return nothing
    return state, state == state.sibling ? n : state.sibling
end

function showedges(io::IO, parent::Node, printfunc = identity)
    str = printfunc(parent.data)
    if str != nothing
        if isleaf(parent)
            println(io, str, " has no children")
        else
            print(io, str, " has the following children: ")
            for c in parent
                print(io, printfunc(c.data), "    ")
            end
            print(io, "\n")
            for c in parent
                showedges(io, c, printfunc)
            end
        end
    end
end
showedges(parent::Node) = showedges(stdout, parent)

end