text:
# Algorithm 3.1
# Neville polynomial interpolation.

function neville(x, y, x0)
	n = length(x)
	q = zeros(n, n)
	for i ∈ 1:n
		q[i, 1] = y[i];
	end
	for i ∈ 2:n
		for j ∈ i:n
			q[j,i] = ((x0 - x[j-i+1]) * q[j,i-1] - (x0 - x[j]) * q[j-1,i-1]) / (x[j] - x[j-i+1])
		end
	end
	return q[n,n]
end

println(neville([1, 56, 73], [67, 5, 35], 3))


text:
using Test
using BenchmarkTools

using SPNetworks
import SPNetworks: SumNode, ProductNode, CategoricalDistribution, IndicatorFunction

@testset "SPNetworks.jl" begin

@testset "Creation, evaluation, sampling" begin
    include("testspn.jl")
#     include("testgspn.jl")
    # include("testeval.jl")
end

# I/O
# include("testio.jl")

# MAP Inference
# @testset "MAP Infrence" begin
#     include("testmaxproduct.jl")
#     include("testbeliefprop.jl")
#     # include("testspn2milp.jl")
# end

#@testset "Parameter learning" begin
# Not working
#  include("testparamlearn.jl")
#   include("testparamlearngspn.jl")
#end

# @testset "Structure learning" begin
# Not working
#   include("testlearnspn.jl")
#   include("testchowliu.jl")
# end

# @testset "Dense network generation" begin
# Not working
#   include("testchain.jl")
#   include("testdense.jl")
# end

end

ext: jl


text:
const RIDLER_CALVARD_MAX_ITER = 100
const RIDLER_CALVARD_TOL = 1e-6

abstract type ThresholdMethod end
struct RidlerCalvard    <: ThresholdMethod end
struct EdgeDetection    <: ThresholdMethod end
struct MinimumThreshold <: ThresholdMethod end

threshold(polim::PolarImage) = threshold(polim, RidlerCalvard())
threshold(polim::PolarImage, ::RidlerCalvard) = threshold(pixels(polim), RidlerCalvard())

# Ridler Calvard method
# ---------------------

# adopted with permission from Jason Merrill for MIT license
function threshold(gray::AbstractArray, ::RidlerCalvard)

    # Single pass over input for both high and low mean is
    # almost 5 times faster because no temporary array allocation.
    function highlowmean(gray, thresh)
        lowmean  = zero(thresh)
        highmean = zero(thresh)
        lowcount = 0
        highcount= 0
        for pixel in gray
            if pixel < thresh
                lowmean += pixel
                lowcount += 1
            else
                highmean += pixel
                highcount += 1
            end
        end
        return (highmean / highcount, lowmean / lowcount)
    end

    thresh = mean(gray)
    count = 1

    while count < RIDLER_CALVARD_MAX_ITER
        high, low  = highlowmean(gray, thresh)
        thresh_old = thresh
        thresh = (high + low) / 2
        abs(thresh - thresh_old) < RIDLER_CALVARD_TOL && break
        count += 1
    end
    return(thresh)
end

# Edge detection method
# ---------------------

"Specialized type for a fast circular queue."
mutable struct circqueue{T}
    array::Vector{T}
    len::Int
    index::Int #to be replaced
end
circqueue(A::AbstractVector) = circqueue{eltype(A)}(A, length(A), 1)
function pushshift!(f::circqueue{T}, input::T) where T
    ind = f.index
    len = f.len
    output = f.array[ind]
    f.array[ind] = input
    f.index = ifelse(ind == f.len, 1, ind + 1) #ifelse faster than ternary
    output
end

"Fitness function to be optimized for Edge Detection method."
function edgefoptim(im, th)
    t = convert(eltype(im), th)
    s = StreamMean()
    s = update(s, 0) # in case no edges

    # We cache the first row and then
    # work with circular dequeues for speed, using pushshift!
    imq = circqueue(im[:,1])
    thq = circqueue(im[:,1] .< t)

    for j = 2:size(im,2)
        # p stands for pixel, b for boolean;
        #`prev` is one up, `up` is upper left, `left` is one left.
        prevp = im[1,j]
        prevb = prevp < t
        upp = pushshift!(imq, prevp)
        upb = pushshift!(thq, prevb)
        for i = 2:size(im,1)
            p = im[i,j]
            b = p < t
            leftp = pushshift!(imq, p)
            leftb = pushshift!(thq, b)

            if b != upb; s = update(s, abs(p - upp)); end
            if upb != prevb; s = update(s, abs(upp - prevp)); end
            if upb != leftb; s = update(s, abs(upp - leftp)); end
            if leftb != prevb; s = update(s, abs(leftp - prevp)); end

            prevp, prevb = p, b
            upp, upb = leftp, leftb
        end
    end
    # The published algorithm uses:
    #    out = mean(s)/2
    # but we multiply the contrast mean with the sqrt (because 2D)
    # of edges count, to avoid spurious results with high threshold
    # values and very little edges.
    out = mean(s) * sqrt(length(s))
end

function threshold(gray::AbstractArray, ::EdgeDetection)
    # maximization so use negative sign
    res = optimize(x -> -edgefoptim(gray, x), 0.01, 0.99)
    minimizer(res)
end

# Cut out box around fθρ(π/2) to reduce for polarimage argument
"""
    threshold(polar_image, EdgeDetection())

Edge Detection method for automatic thresholding after Nobis & Hunziker, 2005.

It optimizes the threshold to find the maximum value of the contrast mean at edges.
The method has been slightly adapted from the original by normalizing the contrast mean with the
sqrt (because 2D) of edges count, to avoid spurious results with high threshold
values and very little edges.
"""
function threshold(polim::PolarImage, ::EdgeDetection)
    rows, cols = size(polim)
    Rmax   = ceil(Int, polim.cl.fθρ(π/2))
    ci     = polim.cl.ci
    cj     = polim.cl.cj
    rowmin = max(1, ci - Rmax)
    rowmax = min(rows, ci + Rmax)
    colmin = max(1, cj - Rmax)
    colmax = min(cols, cj + Rmax)
    threshold(polim.img[rowmin:rowmax, colmin:colmax], EdgeDetection())
end

# Minimum method
# --------------

# Check if histogram is bimodal by detecting change in direction.
# TODO rewrite with diff
function isbimodal(hc) #hc=histcounts
    prevup = hc[2] > hc[1]
    prevc = hc[2]
    mode = 0
    for i = 3:length(hc)
        c = hc[i]
        down = c < prevc

        if prevup & down
            mode += 1
        end

        if mode > 2
            return false
        end

        if !down
            prevup = ifelse(c == prevc, prevup, !down)
        else
            prevup = !down
        end
        prevc = c
    end
    return true
end
# find the minimim between modes if bimodal
function bimodalmin(hc)
    @assert isbimodal(hc)
    prevup = hc[2] > hc[1]
    prevc = hc[2]
    mode = 0
    for i = 3:length(hc)
        c = hc[i]
        down = c < prevc

        if prevup && down
            mode += 1
        end

        if (mode == 1) && (c > prevc)
            return i - 1
        end

        if !down
            prevup = ifelse(c == prevc, prevup, !down)
        else
            prevup = !down
        end
        prevc = c
    end
    return true
end

"""
See paper Glasbey 1993 Analysis of Histogram-Based Thresholding Algorithms.
It does not work for integer, because isbimodal can get stuck in instability.
Smooths the vector according to ``y_i = (y_{i-1} + y_i + y_{i+1})/3```.
"""
function smooth_hist!(hc::AbstractArray{<:AbstractFloat})
    hc[end] = (hc[end-1] + hc[end]) / 3
    prev = hc[1]
    c = hc[2]
    #assume hc[0] = 0 as in paper
    hc[1] = (prev + c)/3
    @inbounds for i = 2:(length(hc) - 1)
        nxt = hc[i+1]
        hc[i] = (prev + c + nxt) / 3
        prev = c
        c = nxt
    end
end
# precision does not seem to increase by increasing binsize
function threshold(img, ::MinimumThreshold; bins=256, maxiter=10_000)
    hist_range = -1 / (bins-1) : 1 / (bins-1) : 1
    counts = fasthist(reshape(img, length(img)), hist_range)
    counts = float(counts) # float required for convergence!
    cnt = 0
    while cnt < maxiter
        cnt += 1
        smooth_hist!(counts)
        isbimodal(counts) && break
    end
    cnt == maxiter && @warn("maximum iteration reached at $cnt in minimum_threshold")
    th = bimodalmin(counts)/bins
end
threshold(polim::PolarImage, ::MinimumThreshold) = threshold(pixels(polim), MinimumThreshold())


text:
# from Tim Holy (Compose/color.jl)
# linear interpolation in [a, b] where x is in [0,1],
# or coerced to be if not.
# !@_{}_@! UncheckedPrecondition: a <= b
function lerp01(x::T, a::T, b::T) where {T<:Number}
    a + (b - a) * max(min(x, 1.0), 0.0)
end

# linear interpolation in [a, b]
# !@_{}_@! UncheckedPrecondition: a <= b
function lerp(x::T, a::T, b::T) where {T<:Number}
    a + (b - a) * x
end

# Clamping limits the values allowable
# !@_{}_@! UncheckedPrecondition: a <= b
function clamp(x::T, a::T, b::T) where {T<:Number}
   a <= x <= b && return x
   x < a ? a : b
end

function locf(vec::NullableNumVec, fillback::Bool)
   idx = index_first_nonnull(vec)
   v = locf(vec)
   if !fillback && idx > 1
      v[1:idx-1] = vec[1:idx-1]
   end
   return v
end

function locf(vec::NullableNumVec)
   v = copy(vec)
   locf!(v)
   return v
end

function locf!(vec::NullableNumVec, fillback::Bool)
   idx = index_first_nonnull(vec)
   if !fillback && idx > 1
      vecstart = vec[1:idx-1]
   end
   locf!(vec)
   if !fillback && idx > 1
       vec[1:idx-1] = vecstart
   end
   return nothing
end

function locf!(vec::NullableNumVec)
    n = length(vec)
    vecidxs = 1:n

    if isnull(vec[1])
       idx = index_first_nonnull(view(vec, vecidxs))
       if idx != 0
           vec[1:idx-1] = vec[idx]
       end
    end
    if isnull(vec[end])
       idx = index_final_nonnull(view(vec, vecidxs))
       if idx != 0
           vec[idx+1:end] = vec[idx]
       end
    end

   if any(isnull.(vec))
       n = length(vec)
       nulls_at = index_nulls(vec)
       if n > length(nulls_at)
           vec[nulls_at] = locf_values(view(vec, vecidxs), view(nulls_at,1:length(nulls_at)))
       end
    end

    return nothing
end

function locf_values(vec::NullableNumVec, nulls_at::AbstractIntVec)
    augment = 0
    deltas = [nulls_at[1]-1, diff(nulls_at)...]

    for i in 2:length(deltas)
        if deltas[i] == 1
            augment += 1
            deltas[i] = 0
        else
            deltas[i] += augment
            augment = 0
        end
    end

    return cumsum(deltas)
end

function index_nulls(vec::NullableNumVec)
    idxs = 1:length(vec)
    nulls = map(isnull, view(vec,idxs))
    return idxs[nulls]
end

"""
    index_first_nonnull(vec)

returns 0 iff all elements of vec are Null
"""
function index_first_nonnull(vec::NullableNumVec)
   result = 0
   for i in 1:length(vec)
      if !isnull(vec[i])
         result = i
         break
      end
   end
   return result
end

"""
    index_final_nonnull(vec)

returns 0 iff all elements of vec are Null
"""
function index_final_nonnull(vec::NullableNumVec)
   result = 0
   for i in length(vec):-1:1
      if !isnull(vec[i])
         result = i
         break
      end
   end
   return result
end


text:
@inline size(s::SignedDistanceField) = size(s.data)
@inline size(s::SignedDistanceField, i...) = size(s.data, i...)
@inline getindex(s::SignedDistanceField, i) = getindex(s.data,i)
@inline getindex(s::SignedDistanceField, i, j) = getindex(s.data,i, j)
@inline getindex(s::SignedDistanceField, i, j, k) = getindex(s.data,i, j,k)
@inline getindex(s::SignedDistanceField, i...) = getindex(s.data,i...)

@inline HyperRectangle(s::SignedDistanceField) = s.bounds

"""
Construct a `SignedDistanceField` by sampling a function over the `bounds`
at the specified `resolution` (default = 0.1). Note that the sampling grid
must be regular,
so a new HyperRectangle will be generated for the SignedDistanceField that
may have larger maximum bounds than the input HyperRectangle. The default
Field type is Float64, but this can be changed with the `fieldT` argument.
"""
function SignedDistanceField(f::Function,
                             bounds::HyperRectangle{3, T},
                             resolution=0.1,
                             fieldT=Float64) where T
    x_min, y_min, z_min = minimum(bounds)
    x_max, y_max, z_max = maximum(bounds)

    x_rng, y_rng, z_rng = maximum(bounds) - minimum(bounds)

    nx = ceil(Int, x_rng/resolution)
    ny = ceil(Int, y_rng/resolution)
    nz = ceil(Int, z_rng/resolution)

    vol = Array{fieldT}(undef, nx+1, ny+1, nz+1)

    nb_max = Vec(x_min + resolution*nx,
                 y_min + resolution*ny,
                 z_min + resolution*nz)

    for i = 0:nx, j = 0:ny, k = 0:nz
        x = x_min + resolution*i
        y = y_min + resolution*j
        z = z_min + resolution*k
        @inbounds vol[i+1,j+1,k+1] = f(Vec{3,fieldT}(x,y,z))
    end

    nb_min = minimum(bounds)
    SignedDistanceField{3,T,fieldT}(HyperRectangle{3, T}(nb_min, nb_max-nb_min), vol)
end

function SignedDistanceField(f::Function,
                             bounds::HyperRectangle{2,T},
                             resolution=0.1,
                             fieldT=Float64) where T
    x_min, y_min = minimum(bounds)
    x_max, y_max = maximum(bounds)

    x_rng, y_rng = maximum(bounds) - minimum(bounds)

    nx = ceil(Int, x_rng/resolution)
    ny = ceil(Int, y_rng/resolution)

    vol = Array{fieldT}(undef, nx+1, ny+1)

    nb_max = Vec(x_min + resolution*nx,
                 y_min + resolution*ny)

    for i = 0:nx, j = 0:ny
        x = x_min + resolution*i
        y = y_min + resolution*j
        @inbounds vol[i+1,j+1] = f(Vec{2,fieldT}(x,y))
    end

    nb_min = minimum(bounds)
    SignedDistanceField{2,T,fieldT}(HyperRectangle{2, T}(nb_min, nb_max-nb_min), vol)
end


text:
@testset "test: newick parsing" begin
@testset "readTopology White Symbol Tests" begin
	#Test newlines, spaces, tabs, and carriage returns
	n1 = readTopology("(A,((B,#H1),(C,(D)#H1)));")
	n2 = readTopology("(A\n,\t (\n\r (B , #H1 ),( C ,(D )#H1 ) \n\n\n\n\n) );")
	@test writeTopology(n1) == writeTopology(n2)

	#Test spaces in troublesome locations
	n3 = readTopology("( A, ( ( B , #H 1 ) ,( C , ( D )#H1 ) ) ) ;")
	@test writeTopology(n1) == writeTopology(n3)

	#Test spaces in names
	resultant = readTopology("('Homo sapiens', ((B,#H1),(C,(D)#H1)));")
	expected = readTopology("('Homosapiens',((B,#H1),(C,(D)#H1)));")
	@test writeTopology(resultant) == writeTopology(expected)

	#Test other white symbols in names
	resultant = readTopology("('H\tom\ro\nsa   p \n\n\n\n\n\n\n\ni\t\t\t\t\t\t\t\ten s   ', ((B,#H1),(C,(D)#H1)));")
	expected = readTopology("('Homosapiens',((B,#H1),(C,(D)#H1)));")
	@test writeTopology(resultant) == writeTopology(expected)
end
@testset "isMajor and gamma consistency" begin
	net = readTopology("((((B)#H1)#H2,((D,C,#H2:::0.8)S1,(#H1,A)S2)S3)S4);");
	@test writeTopology(net, round=true, digits=8) == "(#H2:::0.2,((D,C,((B)#H1)#H2:::0.8),(#H1,A)));"
	net = readTopologyLevel1("(E,((B)#H1:::.5,((D,C),(#H1:::.5,A))));");
	@test writeTopology(net) == "(D:1.0,C:1.0,((#H1:1.0::0.5,A:1.0):1.0,((B:1.0)#H1:1.0::0.5,E:1.0):1.0):1.0);"
    originalSTDOUT = STDOUT
	redirect_stdout(open("/dev/null", "w")) # not portable to Windows
	@test_nowarn PhyloNetworks.printEverything(net)
	redirect_stdout(originalSTDOUT)
end
end


text:
using LevelsetReinitialize
using Test

@testset "LevelsetReinitialize.jl" begin
    # Write your tests here.
end


text:
function initialize!(integrator, cache::MagnusMidpointCache)
  integrator.kshortsize = 2
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # For the interpolation, needs k at the updated point
  integrator.destats.nf += 1
end

function perform_step!(integrator, cache::MagnusMidpointCache, repeat_step=false)
  @unpack t,dt,uprev,u,p,alg = integrator
  @unpack W,k,tmp = cache
  mass_matrix = integrator.f.mass_matrix

  L = integrator.f.f
  update_coefficients!(L,u,p,t+dt/2)

  if integrator.alg.krylov
    u .= expv(dt, L, u; m=min(alg.m, size(L,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
  else
    A = Matrix(L) #size(L) == () ? convert(Number, L) : convert(AbstractMatrix, L)
    u .= exp(dt*L) * u
  end

  integrator.f(integrator.fsallast,u,p,t+dt)
  integrator.destats.nf += 1
end

function initialize!(integrator, cache::LieEulerCache)
  integrator.kshortsize = 2
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # For the interpolation, needs k at the updated point
  integrator.destats.nf += 1
end

function perform_step!(integrator, cache::LieEulerCache, repeat_step=false)
  @unpack t,dt,uprev,u,p,alg = integrator
  @unpack W,k,tmp = cache
  mass_matrix = integrator.f.mass_matrix

  L = integrator.f.f
  update_coefficients!(L,u,p,t)

  if integrator.alg.krylov
    u .= expv(dt, L, u; m=min(alg.m, size(L,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
  else
    A = Matrix(L) #size(L) == () ? convert(Number, L) : convert(AbstractMatrix, L)
    u .= exp(dt*L) * u
  end

  integrator.f(integrator.fsallast,u,p,t+dt)
  integrator.destats.nf += 1
end

function initialize!(integrator, cache::MagnusLeapfrogCache)
  integrator.kshortsize = 2
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # For the interpolation, needs k at the updated point
  integrator.destats.nf += 1
end

function perform_step!(integrator, cache::MagnusLeapfrogCache, repeat_step=false, alg_extrapolates=true, iter=1)
  @unpack t,dt,uprev,uprev2,u,p,alg,iter = integrator
  @unpack W,k,tmp = cache
  mass_matrix = integrator.f.mass_matrix
    # println("iter   : $iter")
  if iter==1
    L = integrator.f.f
    update_coefficients!(L,u,p,t+dt/2)
    if integrator.alg.krylov
      u .= expv(dt, L, u; m=min(alg.m, size(L,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
    else
      A = Matrix(L) #size(L) == () ? convert(Number, L) : convert(AbstractMatrix, L)
      u .= exp(dt*L) * u
    end

    integrator.f(integrator.fsallast,u,p,t+dt)
    integrator.destats.nf += 1
    iter += 1
  else
    L = integrator.f.f
    update_coefficients!(L,u,p,t)
    if integrator.alg.krylov
      u .= expv(2*dt, L, uprev2; m=min(alg.m, size(L,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
    else
      A = Matrix(L) #size(L) == () ? convert(Number, L) : convert(AbstractMatrix, L)
      u .= exp(2*dt*L) * uprev2
    end
    uprev=u
    integrator.f(integrator.fsallast,u,p,t+dt)
    integrator.destats.nf += 1
  end
end

function initialize!(integrator, cache::LinearExponentialConstantCache)
  # Pre-start fsal
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
  integrator.destats.nf += 1
  integrator.fsallast = zero(integrator.fsalfirst)

  # Initialize interpolation derivatives
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

function perform_step!(integrator, cache::LinearExponentialConstantCache, repeat_step=false)
  @unpack t,dt,uprev,f,p = integrator
  alg = unwrap_alg(integrator, true)
  A = f.f # assume f to be an ODEFunction wrapped around a linear operator

  if alg.krylov == :off
    u = exp(dt * Matrix(f)) * integrator.u
  elseif alg.krylov == :simple
    u = expv(dt, A, integrator.u; m=min(alg.m, size(A,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
  else
    u = expv_timestep(dt, A, integrator.u; m=min(alg.m, size(A,1)), iop=alg.iop,
                      opnorm=integrator.opts.internalopnorm, tol=integrator.opts.reltol)
  end

  # Update integrator state
  integrator.fsallast = f(u, p, t + dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator, cache::LinearExponentialCache)
  # Pre-start fsal
  integrator.fsalfirst = zero(cache.rtmp)
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t)
  integrator.destats.nf += 1
  integrator.fsallast = zero(integrator.fsalfirst)

  # Initialize interpolation derivatives
  integrator.kshortsize = 2
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

function perform_step!(integrator, cache::LinearExponentialCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack tmp, KsCache = cache
  alg = unwrap_alg(integrator, true)
  A = f.f # assume f to be an ODEFunction wrapped around a linear operator

  if alg.krylov == :off
    E = exp(dt * Matrix(A))
    mul!(tmp, E, u)
  elseif alg.krylov == :simple
    Ks, expv_cache = KsCache
    arnoldi!(Ks, A, u; m=min(alg.m, size(A,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
    expv!(tmp, dt, Ks; cache=expv_cache)
  else
    expv_timestep!(tmp, dt, A, u; adaptive=true, caches=KsCache, m=min(alg.m, size(A,1)), iop=alg.iop,
                   opnorm=integrator.opts.internalopnorm, tol=integrator.opts.reltol)
  end

  # Update integrator state
  u .= tmp
  f(integrator.fsallast, u, p, t + dt)
  integrator.destats.nf += 1
  # integrator.k is automatically set due to aliasing
end

cay!(tmp, A) = mul!(tmp, inv(I - 1/2 * A), (I + 1/2 * A))
cay(A) = inv(I - 1/2 * A) * (I + 1/2 * A)

function initialize!(integrator, cache::CayleyEulerConstantCache)
  # Pre-start fsal
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
  integrator.destats.nf += 1
  integrator.fsallast = zero(integrator.fsalfirst)

  # Initialize interpolation derivatives
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

function perform_step!(integrator, cache::CayleyEulerConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  alg = unwrap_alg(integrator, true)

  if f isa SplitFunction
    A = f.f1.f
  else  # f isa ODEFunction
    A = f.f
  end

  L = update_coefficients(A, uprev, p, t)
  V = cay(L*dt)
  u = V * uprev * transpose(V)

  # Update integrator state
  integrator.fsallast = f(u, p, t + dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator, cache::CayleyEulerCache)
  integrator.kshortsize = 2
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # For the interpolation, needs k at the updated point
  integrator.destats.nf += 1
end

function perform_step!(integrator, cache::CayleyEulerCache, repeat_step=false)
  @unpack t,dt,uprev,u,p,f,alg = integrator
  @unpack k,V,tmp = cache
  mass_matrix = integrator.f.mass_matrix

  if f isa SplitFunction
    L = f.f1.f
  else  # f isa ODEFunction
    L = f.f
  end

  update_coefficients!(L, uprev, p, t)

  cay!(V, L*dt)
  mul!(tmp, uprev, transpose(V))
  mul!(u, V, tmp)

  # Update integrator state
  integrator.f(integrator.fsallast,u,p,t+dt)
  integrator.destats.nf += 1
end


text:
using Libdl
using PackageCompiler

include("build_IntelOpenMP.jl")
include("build_MKL.jl")
include("deps.jl")

include("../src/install.jl")
enable_mkl_startup(libmkl_rt)

try
    julia_ = joinpath(Sys.BINDIR, Base.julia_exename())
    cmd = `$julia_ fix_repl.jl`
    run(cmd)
catch
    @warn "Rebuilding system image with precompiling failed. This may lead to REPL latency."
end


text:
# This file was generated, do not modify it. # hide
describe(auto, :mean, :median, :std)


text:
to_affine_map(tform::Transform3D) = AffineMap(rotation(tform), translation(tform))

rotation_from_x_axis(translation::AbstractVector{T}) where {T} = rotation_between(SVector{3, T}(1,0,0), translation)

function inertial_ellipsoid_dimensions(mass, axis_inertias)
    # Ix = m/5 (dy^2 + dz^2)
    # Iy = m/5 (dx^2 + dz^2)
    # Iz = m/5 (dx^2 + dy^2)
    #
    # let A = [0 1 1
    #          1 0 1
    #          1 1 0]
    # b = 5 / m * [Ix; Iy; Iz]
    # Then A \ b = [dx^2; dy^2; dz^2]
    #
    # This is only valid if the axis inertias obey the triangle inequalities:
    # Ix + Iy >= Iz
    # Ix + Iz >= Iy
    # Iy + Iz >= Ix

    # Ix - Iy = m/5 (dy^2 - dx^2)
    # Ix - Iy + Iz = m/5 (2*dy^2)
    # dy^2 = 0.5 (Ix - Iy + Iz) * 5/m

    squared_lengths = 0.5 * 5.0 / mass *
        [-axis_inertias[1] + axis_inertias[2] + axis_inertias[3];
          axis_inertias[1] - axis_inertias[2] + axis_inertias[3];
          axis_inertias[1] + axis_inertias[2] - axis_inertias[3]]

    for i = 1:3
        total_inertia_of_other_axes = zero(axis_inertias[1])
        for j = 1:3
            if i == j
                continue
            end
            total_inertia_of_other_axes += axis_inertias[j]
        end
        if axis_inertias[i] > total_inertia_of_other_axes
            error("Principal inertias $(axis_inertias) do not satisfy the triangle inequalities, so the equivalent inertial ellipsoid is not well-defined.")
        end
    end

    return sqrt.(squared_lengths)
end

inertial_ellipsoid(body::RigidBody) = inertial_ellipsoid(spatial_inertia(body))

if VERSION < v"0.7.0-DEV.5211"
    const _eigen = eig
else
    const _eigen = eigen
end

function inertial_ellipsoid(inertia::SpatialInertia)
    com_frame = CartesianFrame3D("CoM")
    com_frame_to_inertia_frame = Transform3D(com_frame, inertia.frame, center_of_mass(inertia).v)
    com_inertia = transform(inertia, inv(com_frame_to_inertia_frame))
    principal_inertias, axes = _eigen(Array(com_inertia.moment)) # StaticArrays.eig checks that the matrix is Hermitian with zero tolerance...
    axes[:,3] *= sign(dot(cross(axes[:,1], axes[:,2]), axes[:,3])) # Ensure the axes form a right-handed coordinate system
    radii = inertial_ellipsoid_dimensions(com_inertia.mass, principal_inertias)
    # We create an ellipsoid by generating a sphere and then scaling it
    # along each axis
    geometry = HyperSphere(zero(Point{3, Float64}), 1.0)
    scaling = LinearMap(SDiagonal(radii[1], radii[2], radii[3]))
    tform = AffineMap(RotMatrix{3}(axes), center_of_mass(inertia).v) ∘ scaling
    return VisualElement(inertia.frame, geometry, DEFAULT_COLOR, tform)
end

function create_frame_to_frame_geometry(joint_to_joint, radius)
    trans = translation(joint_to_joint)
    geom_length = norm(trans)
    radius = min(radius, geom_length)
    tform = if norm(trans) > 1e-10
        LinearMap(rotation_from_x_axis(trans))
    else
        IdentityTransformation()
    end
    geometry = HyperRectangle(Vec(0, -radius, -radius), Vec(geom_length, 2*radius, 2*radius))
    return VisualElement(joint_to_joint.to, geometry, DEFAULT_COLOR, tform)
end

function maximum_link_length(body_fixed_joint_frames::Dict{RigidBody{T}, Vector{CartesianFrame3D}}) where T
    result = zero(T)
    for (body, joint_frames) in body_fixed_joint_frames
        for framei in joint_frames, framej in joint_frames
            transform = fixed_transform(body, framei, framej)
            result = max(result, norm(translation(transform)))
        end
    end
    result
end

struct Skeleton <: AbstractGeometrySource
    inertias::Bool
    randomize_colors::Bool
end

Skeleton(; inertias=true, randomize_colors=false) = Skeleton(inertias, randomize_colors)

function visual_elements(mechanism::Mechanism, source::Skeleton)
    body_fixed_joint_frames = Dict(body => begin
        [map(frame_before, out_joints(body, mechanism)); map(frame_after, in_joints(body, mechanism))]
    end for body in bodies(mechanism))

    box_width = 0.05 * maximum_link_length(body_fixed_joint_frames)
    elements = Vector{VisualElement}()

    for body in bodies(mechanism)
        if source.inertias && has_defined_inertia(body) && spatial_inertia(body).mass >= 1e-3
            push!(elements, inertial_ellipsoid(body))
        else
            for joint in out_joints(body, mechanism)
                if !iszero(box_width)
                    push!(elements, VisualElement(
                        frame_before(joint),
                        HyperSphere{3, Float64}(zero(Point{3, Float64}), box_width),
                        DEFAULT_COLOR,
                        IdentityTransformation()
                    ))
                end
            end
        end
        frames = body_fixed_joint_frames[body]
        for (i, framei) in enumerate(frames)
            for j = i + 1 : length(frames)
                framej = frames[j]
                joint_to_joint = fixed_transform(mechanism, framei, framej)
                push!(elements, create_frame_to_frame_geometry(
                    joint_to_joint,
                    box_width / 2
                ))
            end
        end
    end
    if source.randomize_colors
        for element in elements
            element.color = RGBA{Float32}(rand(), rand(), rand(), 0.5)
        end
    end
    elements
end


text:
# GSA code v0.1.
# Coded by Jesús Mejía.
# Based on MATLAB code of Esmat Rashedi, 2010.
# " E. Rashedi, H. Nezamabadi-pour and S. Saryazdi,
# “GSA: A Gravitational Search Algorithm”, Information sciences, vol. 179,
# no. 13, pp. 2232-2248, 2009."

function massCalculation(fit,searchType)
	####here, make your own function of 'mass calculation'

	Fmax = maximum(fit)
	Fmin = minimum(fit)
	Fmean= mean(fit)
	N = length(fit)

	if Fmax == Fmin
		M = ones(N)
	else

		if searchType == :minimize #for minimization
			best = Fmin
			worst= Fmax #eq.17-18.
		else #for maximization
			best = Fmax
			worst= Fmin #eq.19-20.
		end

		M = (fit .- worst) ./ (best .- worst) #eq.15,

	end

	M = M ./ sum(M) #eq. 16.

	return M

end

function Gconstant(iteration,max_it)
	# here, make your own function of 'G'
	α = 20
	G0= 100

	G = G0*exp(-α * iteration / max_it) #eq. 28.

	return G
end

#This function calculates the accelaration of each agent in gravitational field. eq.7-10,21.
function Gfield(M,X,G,Rnorm,Rpower,ElitistCheck,iteration,max_it)
	N, D   = size(X)
	final_per= 2 #In the last iteration, only 2 percent of agents apply force to the others.

	####total force calculation
	if ElitistCheck == 1
		kbest = final_per+(1-iteration/max_it)*(100-final_per) #kbest in eq. 21.
		kbest = round(Int, N*kbest/100)
	else
		kbest = N #eq.9.
	end

	ds = sortperm(M, rev=true)

	E = zeros(N, D)
	for i=1:N
		for ii=1:kbest
			j = ds[ii]
			if j != i
				R = norm(X[i,:]-X[j,:], Rnorm) #Euclidian distanse.
				for k=1:D
                    E[i, k]= E[i,k]+rand() * (M[j]) * ((X[j,k]-X[i,k])/(R^Rpower + eps()))
					#note that Mp(i)/Mi(i)=1
				end
			end
		end
	end

	##acceleration
	a = E.*G #note that Mp(i)/Mi(i)=1

	return a
end

# This function updates the velocity and position of agents.
function move(X,a,V)
	# movement.
	N, D = size(X)
	V = rand(N, D).*V + a # eq. 11.
	X = X + V # eq. 12.

	return X, V
end


text:
"""
	search(srchers, query [;kwargs])

Searches for query (i.e. key terms) in multiple searches and returns
information regarding the documents that match best the query.
The function returns the search results in the form of
a `Vector{SearchResult}`.

# Arguments
  * `srchers::Vector{Searcher}` is the searchers vector
  * `query` the query, can be either a `String` or `Vector{String}`

# Keyword arguments
  * `search_method::Symbol` controls the type of matching: `:exact`
     uses exact matches while `:regex` consideres the needle a regular expression
  * `max_matches::Int` is the maximum number of search results to return
  * `max_suggestions::Int` is the maximum number of suggestions to return for
     each missing needle
  * `custom_weights::Dict{Symbol, Float}` are custom weights for each
     searcher's results used in result aggregation
"""
function search(srchers::Vector{<:AbstractSearcher{T}},
                query;
                search_method=DEFAULT_SEARCH_METHOD,
                max_matches=MAX_MATCHES,
                max_suggestions=MAX_SUGGESTIONS,
                custom_weights=DEFAULT_CUSTOM_WEIGHTS
               ) where {T<:AbstractFloat}
    # Checks
    @assert search_method in [:exact, :regex]
    @assert max_matches >= 0
    @assert max_suggestions >=0

    # Initializations
    n = length(srchers)
    enabled_searchers = [i for i in 1:n if isenabled(srchers[i])]
    n_enabled = length(enabled_searchers)
    queries = fill(query, n)

    # Search
    results = Vector{SearchResult{T}}(undef, n_enabled)
    ###################################################################
    # A `Threads.@threads` statement in front of the for loop here
    # idicates the use of multi-threading. If multi-threading is used,
    # OPENBLAS multi-threading support has to be disabled by using:
    #   `export OPENBLAS_NUM_THREADS=1` in the shell
    # or start julia with:
    #   `env OPENBLAS_NUM_THREADS=1 julia`
    #
    # WARNING: Multi-theading support (as of v1.1 is still EXPERIMENTAL)
    #          and floating point operations are not thread-safe!
    #          Do not use with semantic search!!
    ###################################################################
    ### Threads.@threads for i in 1:n_enabled
    for i in 1:n_enabled
        results[i] = search(srchers[enabled_searchers[i]],
                            queries[enabled_searchers[i]],
                            search_method=search_method,
                            max_matches=max_matches,
                            max_suggestions=max_suggestions)
    end
    # Aggregate results
    ids_agg = String[srcher.config.id_aggregation for srcher in srchers if isenabled(srcher)]
    aggregate!(results,
               ids_agg,
               method=RESULT_AGGREGATION_STRATEGY,
               max_matches=max_matches,
               max_suggestions=max_suggestions,
               custom_weights=custom_weights)

    # Return results
    return results
end

"""
	search(srcher, query [;kwargs])

Searches for query (i.e. key terms) in `srcher`, and returns information
regarding the the documents that match best the query. The function
returns an object of type `SearchResult`.

# Arguments
  * `srcher::Searcher` is the searcher
  * `query` the query, can be either a `String` or `Vector{String}`

# Keyword arguments
  * `search_method::Symbol` controls the type of matching: `:exact`
     uses exact matches while `:regex` consideres the needle a regular expression
  * `max_matches::Int` is the maximum number of search results to return
  * `max_suggestions::Int` is the maximum number of suggestions to return for
     each missing needle
"""
function search(srcher::Searcher{T,E,I},
                query;  # can be either a string or vector of strings
                search_method::Symbol=DEFAULT_SEARCH_METHOD,
                max_matches::Int=MAX_MATCHES,
                max_suggestions::Int=MAX_SUGGESTIONS
               ) where {T<:AbstractFloat, E, I<:AbstractIndex}
    # Make an entry i.e.e NamedTuple out of the query
    qentry = (query=query2string(query),)  # query treated as single sentence
    qfields = [:query]
    isregex = search_method === :regex
    qembedding, status = embed(srcher.input_embedder[], [qentry]; fields=qfields, isregex=isregex)

    # Search (if document vector is not zero)
    scores, idxs = T[], Int[]
    if first(status)
        idxs, scores = knn_search(srcher.index,
                                  firstcol(qembedding),
                                  min(length(srcher.index), max_matches))  # k
        score_transform!(scores, alpha=srcher.config.score_alpha)
    end
    query_matches = collect(zip(scores, idxs))

    # Find matching and missing needles
    needles = query2vector(query)
    needle_matches, missed_needles = missing_needles(srcher.input_embedder[], needles)

    # Get suggestions
    suggestions = MultiDict{String, Tuple{T, String}}()
    if max_suggestions > 0 && !isempty(missed_needles)
        suggestion_search!(suggestions,
                           srcher.search_trees,
                           missed_needles,
                           max_suggestions=max_suggestions)
    end
    return SearchResult(id(srcher), query_matches, needle_matches,
                        suggestions, srcher.config.score_weight)
end

"""
Returns found and missing needles using an embedder
"""
missing_needles

missing_needles(embedder::WordVectorsEmbedder, needles) = (String[], String[])

missing_needles(embedder::DTVEmbedder, needles) =
    missing_needles(embedder.model, needles)

missing_needles(model::StringAnalysis.RPModel, needles) = begin
    needle_matches = [needle for needle in needles if in(needle, model.vocab)]
    missed_needles = setdiff(needles, needle_matches)
    return needle_matches, missed_needles
end

missing_needles(model::StringAnalysis.LSAModel, needles) = (String[], String[])

"""
    suggestion_search!(suggestions, search_tree, needles [;max_suggestions=1])

Searches in the search tree for partial matches for each of  the `needles`.
"""
function suggestion_search!(suggestions::MultiDict{String, Tuple{T, String}},
                            search_tree::BKTree{String},
                            needles::Vector{S};
                            max_suggestions::Int=1
                           ) where {S<:AbstractString, T<:AbstractFloat}
    if isempty(needles)
        return suggestions
    elseif BKTrees.is_empty_node(search_tree.root)
        @debug "Suggestion tree is empty, no suggestions will be added."
        return suggestions
    else  # there are terms that have not been found
        for needle in needles
            _suggestions = sort!(find(search_tree, String(needle),
                                      MAX_EDIT_DISTANCE,
                                      k=max_suggestions),
                                 by=t->t[1])
            if !isempty(_suggestions)
                n = min(max_suggestions, length(_suggestions))
                push!(suggestions, needle=>_suggestions[1:n])
            end
        end
    end
    return suggestions
end

# Post-processing score function:
#   - map distances [0, Inf) --> [1, 0)
#TODO(Corneliu) Analylically/empirically adapt alpha do vector dimensionality
function score_transform!(x::AbstractVector{T};
                          alpha::T=T(DEFAULT_SCORE_ALPHA),
                          normalize=false) where {T<:AbstractFloat}
    n = length(x)
    @inbounds @simd for i in 1:n
        x[i] = 1 - tanh(alpha * x[i])
    end
    if normalize
        # Forces the scores to spread in the whole
        # interval between 0 and 1
        xmax = maximum(x)
        xmin = minimum(x)
        return (x.-xmin)./(xmax-xmin)
    else
        return x
    end
end

# Useful functions
query2string(query::AbstractString) = String(query)

query2string(query::Vector{String}) = join(query, " ")

query2string(query) = throw(ArgumentError("query2string requires `String` or Vector{String} inputs."))

query2vector(query::Vector{String}) = query

query2vector(query::AbstractString) = String.(tokenize(query, method=DEFAULT_TOKENIZER))

query2vector(query) = throw(ArgumentError("query2vector requires `String` or Vector{String} inputs."))


text:
module SQLAlchemy

export Table, Column, MetaData, Engine, Session
export createengine, select, text, connect, func, inspect, query
export createall, insert, values, compile, connect, execute, fetchone, fetchall, where, selectfrom, and, orderby, alias, join, groupby, having, delete, update, distinct, limit, offset, label, loadchinook, desc, asc, dirty
export SQLString, SQLInteger, SQLBoolean, SQLDate, SQLDateTime, SQLEnum, SQLFloat, SQLInterval, SQLNumeric, SQLText, SQLTime, SQLUnicode, SQLUnicodeText, SQLType
export @table

include("core.jl")
include("record.jl")
include("chinook.jl")
include("orm.jl")

end


text:
# ---------------------------------------------------
module Logs
# ---------------------------------------------------
using Dates

using ..CoreInterfaces
using ..Config

# Simulation writes to several files, not including the sim data files:
# app.log - Smaller more frequent messages
# progress.log - Larger report size messages and information
# sim.log - Stores the current state of the simulation. For example,
#           if it was paused or terminated.
# error.log - Errors prior at the time of failures.

mutable struct SLog <: CoreInterfaces.ILog
    start_time::String
    config::CoreInterfaces.IConfig

    appLog::IOStream
    errLog::IOStream
    dtFormat::DateFormat

    function SLog(config::CoreInterfaces.IConfig)
        o = new()
        # println("Applog: ", Config.appLog(config))
        o.dtFormat = dateformat"Y-mm-dd II:MM:SS"
        o.start_time = Dates.format(Dates.now(), o.dtFormat)
        o.config = config
        o
    end
end

function configure(log::CoreInterfaces.ILog)
    # inspect config file to determine how log files are open
    # If Paused then we open most of the logs for "append"
    if Config.exit_state(log.config) == "Paused"
        log.appLog = open(Config.app_log(log.config), "a")
    else
        log.appLog = open(Config.app_log(log.config), "w")
    end

    log.errLog = open(Config.err_log(log.config), "w")

    # log_to_app(log, "Log configured", true)
end

function log_to_app(log::CoreInterfaces.ILog, msg::String, stamp::Bool = false)
    if stamp
        write(log.appLog, Dates.format(Dates.now(), log.dtFormat), " ")
    end

    write(log.appLog, msg, "\n")
    flush(log.appLog)
end

function log_to_err(log::CoreInterfaces.ILog, msg::String, stamp::Bool = true)
    if stamp
        write(log.errLog, Dates.format(Dates.now(), log.dtFormat), " ")
    end

    write(log.errLog, msg, "\n")
    flush(log.errLog)
end

function close_logs(log::CoreInterfaces.ILog)
    close(log.appLog)
    close(log.errLog)
end

end # End module -----------------


text:
module ReadableNumbers

#=
    generating and showing prettier numeric strings
=#

export ReadableNumStyle, readable, show_readable, parse_readable

import Base: stdout

if typeof(Base.split("a","b")[1]) == SubString{String}
    split(str::String, sep::Char=" ") = map(String, Base.split(str, sep))  # do not work with SubStrings
end

# determine locale convention for the fractional (decimal) point
const LOCALE_STR = string( 1.5 )
const FRACPOINT  = LOCALE_STR[ nextind(LOCALE_STR, 1) ]

include("type.jl")
include("dothework.jl")
include("show.jl")
include("parse.jl")

end # module


text:
import Base: getindex

abstract type AbstractRepresentation end

struct Cartesian <: AbstractRepresentation end
struct Spherical <: AbstractRepresentation end
struct Covariant <: AbstractRepresentation end
struct Contravariant <: AbstractRepresentation end

Base.@kwdef struct VectorField{S, T, C} <: AbstractField
    data::S
    grid::T
    representation::C = Cartesian()
end

function VectorField(ϕ::VectorField; representation = Cartesian())
    return VectorField(
        data = ϕ.data,
        grid = ϕ.grid,
        representation = representation,
    )
end

function (ϕ::VectorField)(representation::AbstractRepresentation)
    return VectorField(ϕ, representation = representation)
end

function components(data, ::Cartesian)
    one = @sprintf("%0.2e", data[1])
    two = @sprintf("%0.2e", data[2])
    three = @sprintf("%0.2e", data[3])
    println(one, "x̂ +", two, "ŷ +", three, "ẑ ")
    return data
end

getindex(ϕ::VectorField, ijk, e; verbose = true) = components(
    getindex.(ϕ.data, ijk, e),
    ϕ.grid,
    ijk,
    e,
    ϕ.representation,
    verbose = verbose,
)

## Component Grabber

function convenience_print(grid, ijk, e)
    print("location: ")
    x, y, z = get_position(grid, ijk, e)
    println(
        "x=",
        @sprintf("%0.2e", x),
        ",y=",
        @sprintf("%0.2e", y),
        ",z=",
        @sprintf("%0.2e", z),
        " ",
    )
    print("field:     ")
end

function components(v⃗, grid, ijk, e, ::Cartesian; verbose = true)
    if verbose
        one = @sprintf("%0.2e", v⃗[1])
        two = @sprintf("%0.2e", v⃗[2])
        three = @sprintf("%0.2e", v⃗[3])
        convenience_print(grid, ijk, e)
        println(one, "x̂ +", two, "ŷ +", three, "ẑ ")
    end
    return v⃗
end

function components(v⃗, grid, ijk, e, ::Covariant; verbose = true)
    v⃗¹, v⃗², v⃗³ = get_contravariant(grid, ijk, e)
    v₁ = dot(v⃗¹, v⃗)
    v₂ = dot(v⃗², v⃗)
    v₃ = dot(v⃗³, v⃗)
    if verbose
        one = @sprintf("%0.2e", v₁)
        two = @sprintf("%0.2e", v₂)
        three = @sprintf("%0.2e", v₃)
        convenience_print(grid, ijk, e)
        println(one, "v⃗¹ +", two, "v⃗² +", three, "v⃗³ ")
    end
    return (; v₁, v₂, v₃)
end

function components(v⃗, grid, ijk, e, ::Contravariant; verbose = true)
    v⃗₁, v⃗₂, v⃗₃ = get_covariant(grid, ijk, e)
    v¹ = dot(v⃗₁, v⃗)
    v² = dot(v⃗₂, v⃗)
    v³ = dot(v⃗₃, v⃗)
    if verbose
        one = @sprintf("%0.2e", v¹)
        two = @sprintf("%0.2e", v²)
        three = @sprintf("%0.2e", v³)
        convenience_print(grid, ijk, e)
        println(one, "v⃗₁ +", two, "v⃗₂ +", three, "v⃗₃ ")
    end
    return (; v¹, v², v³)
end

function components(v⃗, grid, ijk, e, ::Spherical; verbose = true)
    r̂, θ̂, φ̂ = get_spherical(grid, ijk, e)
    vʳ = dot(r̂, v⃗)
    vᶿ = dot(θ̂, v⃗)
    vᵠ = dot(φ̂, v⃗)
    if verbose
        one = @sprintf("%0.2e", vʳ)
        two = @sprintf("%0.2e", vᶿ)
        three = @sprintf("%0.2e", vᵠ)
        convenience_print(grid, ijk, e)
        println(one, "r̂ +", two, "θ̂ +", three, "φ̂ ")
    end
    return (; vʳ, vᶿ, vᵠ)
end

## Helper functions
function get_jacobian(grid, ijk, e)
    ξ1x1 = grid.vgeo[ijk, grid.ξ1x1id, e]
    ξ1x2 = grid.vgeo[ijk, grid.ξ1x2id, e]
    ξ1x3 = grid.vgeo[ijk, grid.ξ1x3id, e]

    ξ2x1 = grid.vgeo[ijk, grid.ξ2x1id, e]
    ξ2x2 = grid.vgeo[ijk, grid.ξ2x2id, e]
    ξ2x3 = grid.vgeo[ijk, grid.ξ2x3id, e]

    ξ3x1 = grid.vgeo[ijk, grid.ξ3x1id, e]
    ξ3x2 = grid.vgeo[ijk, grid.ξ3x2id, e]
    ξ3x3 = grid.vgeo[ijk, grid.ξ3x3id, e]

    J = [
        ξ1x1 ξ1x2 ξ1x3
        ξ2x1 ξ2x2 ξ2x3
        ξ3x1 ξ3x2 ξ3x3
    ]
    # rows are the contravariant vectors
    # columns of the inverse are the covariant vectors
    return J
end

function get_contravariant(grid, ijk, e)
    J = get_jacobian(grid, ijk, e)
    a⃗¹ = J[1, :]
    a⃗² = J[2, :]
    a⃗³ = J[3, :]
    return (; a⃗¹, a⃗², a⃗³)
end

function get_covariant(grid, ijk, e)
    J = inv(get_jacobian(grid, ijk, e))
    a⃗₁ = J[:, 1]
    a⃗₂ = J[:, 2]
    a⃗₃ = J[:, 3]
    return (; a⃗₁, a⃗₂, a⃗₃)
end

function get_spherical(grid, ijk, e)
    x, y, z = get_position(grid, ijk, e)
    r̂ = [x, y, z] ./ norm([x, y, z])
    θ̂ = [x * z, y * z, -(x^2 + y^2)] ./ (norm([x, y, z]) * norm([x, y, 0]))
    φ̂ = [-y, x, 0] ./ norm([x, y, 0])
    return (; r̂, θ̂, φ̂)
end

function get_position(grid, ijk, e)
    x1 = grid.vgeo[ijk, grid.x1id, e]
    x2 = grid.vgeo[ijk, grid.x2id, e]
    x3 = grid.vgeo[ijk, grid.x3id, e]
    r = [x1 x2 x3]
    return r
end

function construct_determinant(grid)
    M = grid.vgeo[:, grid.Mid, :]

    ωx = reshape(grid.ω[1], (length(grid.ω[1]), 1, 1, 1))
    ωy = reshape(grid.ω[2], (1, length(grid.ω[2]), 1, 1))
    ωz = reshape(grid.ω[3], (1, 1, length(grid.ω[3]), 1))
    ω = reshape(ωx .* ωy .* ωz, (size(M)[1], 1))
    J = M ./ ω
    return J
end

function get_jacobian_determinant(grid, ijk, e)
    M = grid.vgeo[ijk, grid.Mid, e]
end


text:
module test

using Roots,Interpolations,PyPlot,Optim
using CompEcon, NLsolve

alpha     = 0.65
beta      = 0.95
grid_max  = 2  # upper bound of capital grid
n         = 150  # number of grid points
N_iter    = 3000  # number of iterations
kgrid     = 1e-6:(grid_max-1e-6)/(n-1):grid_max  # equispaced grid
f(x) = x.^alpha  # defines the production function f(k)
tol = 1e-9

ab        = alpha * beta
c1        = (log(1 - ab) + log(ab) * ab / (1 - ab)) / (1 - beta)
c2        = alpha / (1 - ab)
v_star(k) = c1 .+ c2 .* log(k)
k_star(k) = ab * k.^alpha
c_star(k) = (1-ab) * k.^alpha
ufun(x) = log(x)

# Bellman Operator
# inputs
# `grid`: grid of values of state variable
# `v0`: current guess of value function

# output
# `v1`: next guess of value function
# `pol`: corresponding policy function

#takes a grid of state variables and computes the next iterate of the value function.
function bellman_operator(grid,v0)

    v1  = zeros(n)     # next guess
    pol = zeros(Int,n)     # policy function
    w   = zeros(n)   # temporary vector

    # loop over current states
    # current capital
    for (i,k) in enumerate(grid)

        # loop over all possible kprime choices
        for (iprime,kprime) in enumerate(grid)
            if f(k) - kprime < 0   #check for negative consumption
                w[iprime] = -Inf
            else
                w[iprime] = ufun(f(k) - kprime) + beta * v0[iprime]
            end
        end
        # find maximal choice
        v1[i], pol[i] = findmax(w)     # stores Value und policy (index of optimal choice)
    end
    return (v1,pol)   # return both value and policy function
end

# VFI iterator
#
## input
# `n`: number of grid points
# output
# `v_next`: tuple with value and policy functions after `n` iterations.
function VFI()
    v_init = zeros(n)     # initial guess
    for iter in 1:N_iter
        v_next = bellman_operator(kgrid,v_init)  # returns a tuple: (v1,pol)
        # check convergence
        if maxabs(v_init.-v_next[1]) < tol
            verrors = maxabs(v_next[1].-v_star(kgrid))
            perrors = maxabs(kgrid[v_next[2]].-k_star(kgrid))
            println("discrete VFI:")
            println("Found solution after $iter iterations")
            println("maximal value function error = $verrors")
            println("maximal policy function error = $perrors")
            return v_next
        elseif iter==N_iter
            warn("No solution found after $iter iterations")
            return v_next
        end
        v_init = v_next[1]  # update guess
    end
end

# plot
function plotVFI()
    v = VFI()
    figure("discrete VFI",figsize=(10,5))
    subplot(131)
    plot(kgrid,v[1],color="blue")
    plot(kgrid,v_star(kgrid),color="black")
    xlim(-0.1,grid_max)
    ylim(-50,-30)
    xlabel("k")
    ylabel("value")
    title("value function")

    subplot(132)
    plot(kgrid,kgrid[v[2]])
    plot(kgrid,k_star(kgrid),color="black")
    xlabel("k")
    title("policy function")

    subplot(133)
    plot(kgrid,kgrid[v[2]].-k_star(kgrid))
    title("policy function error")
end

function bellman_operator2(grid,v0)

    v1  = zeros(n)     # next guess
    pol = zeros(n)     # consumption policy function

    Interp = interpolate((collect(grid),), v0, Gridded(Linear()) )

    # loop over current states
    # of current capital
    for (i,k) in enumerate(grid)

        objective(c) = - (log(c) + beta * Interp[f(k) - c])
        # find max of ojbective between [0,k^alpha]
        res = optimize(objective, 1e-6, f(k))
        pol[i] = f(k) - res.minimizer
        v1[i] = -res.minimum
    end
    return (v1,pol)   # return both value and policy function
end

function VFI2()
    v_init = zeros(n)     # initial guess
    for iter in 1:N_iter
        v_next = bellman_operator2(kgrid,v_init)  # returns a tuple: (v1,pol)
        # check convergence
        if maxabs(v_init.-v_next[1]) < tol
            verrors = maxabs(v_next[1].-v_star(kgrid))
            perrors = maxabs(v_next[2].-k_star(kgrid))
            println("continuous VFI:")
            println("Found solution after $iter iterations")
            println("maximal value function error = $verrors")
            println("maximal policy function error = $perrors")
            return v_next
        elseif iter==N_iter
            warn("No solution found after $iter iterations")
            return v_next
        end
        v_init = v_next[1]  # update guess
    end
    return nothing
end

function plotVFI2()
    v = VFI2()
    figure("discrete VFI - continuous control",figsize=(10,5))
    subplot(131)
    plot(kgrid,v[1],color="blue")
    plot(kgrid,v_star(kgrid),color="black")
    xlim(-0.1,grid_max)
    ylim(-50,-30)
    xlabel("k")
    ylabel("value")
    title("value function")

    subplot(132)
    plot(kgrid,v[2])
    plot(kgrid,k_star(kgrid),color="black")
    xlabel("k")
    title("policy function")

    subplot(133)
    plot(kgrid,v[2].-k_star(kgrid))
    xlabel("k")
    title("policy function error")
end

function policy_iter(grid,c0,u_prime,f_prime)

    c1  = zeros(length(grid))     # next guess
    pol_fun = interpolate((collect(grid),), c0, Gridded(Linear()) )

    # loop over current states
    # of current capital
    for (i,k) in enumerate(grid)
        objective(c) = u_prime(c) - beta * u_prime(pol_fun[f(k)-c]) * f_prime(f(k)-c)
        c1[i] = fzero(objective, 1e-10, f(k)-1e-10)
    end
    return c1
end

uprime(x::Float64) = 1.0 / x
fprime(x::Float64) = alpha* x ^ (alpha-1)
function fprime(x::Vector)
    r = similar(x)
    for ix in eachindex(x)
        if x[ix] > 0
            r[ix] = alpha* x[ix]^(alpha-1)
        else
            r[ix] = 10000.0
        end
    end
    return r
end

function PFI()
    c_init = kgrid
    for iter in 1:N_iter
        c_next = policy_iter(kgrid,c_init,uprime,fprime)
        # check convergence
        if maxabs(c_init.-c_next) < tol
            perrors = maxabs(c_next.-c_star(kgrid))
            println("PFI:")
            println("Found solution after $iter iterations")
            println("max policy function error = $perrors")
            return c_next
        elseif iter==N_iter
            warn("No solution found after $iter iterations")
            return c_next
        end
        c_init = c_next  # update guess
    end
end
function plotPFI()
    v = PFI()
    figure("PFI")
    subplot(121)
    plot(kgrid,v)
    plot(kgrid,c_star(kgrid),color="black")
    xlabel("k")
    title("policy function")

    subplot(122)
    plot(kgrid,v.-c_star(kgrid))
    xlabel("k")
    title("policy function error")
end
#plotVFI()
#plotVFI2()
#plotPFI()

function proj2(n=25)

    a = 1e-6
    aa = 1e-2
    b = 1.0
    bb = 0.95
    basis = fundefn(:cheb,n,a,b)
    k = funnode(basis)[1]   # collocation points

    fs = f.(k)

    # system of equations: for each k, one line
    function ff(coef::Vector,result::Vector,b,coll_points,fs)
        # get cons function
        cons = funeval(coef,b,coll_points)[1]
        # println("coef= $coef")
        # println("cons = $cons")

        # put resulting residuals into result[:]
        # put into euler equation
        nextk = fs.-cons
        for i in eachindex(nextk)
            if nextk[i] < 0
                result[i] = -8000.0
            else
                result[i] = uprime(cons[i]) .- beta * uprime(nextk[i]) .* fprime(nextk[i])
            end
        end
        # println(result)
    end

    f_closure(x::Vector,r::Vector) = ff(x,r,basis,k,fs)

    res = nlsolve(f_closure,ones(n)*0.4)
    x = linspace(aa,bb,501)
    y = funeval(res.zero,basis,x)[1]
    plot(x,y)
    return res
end

function proj(n=25)

    alpha = 1.0
    eta   = 1.5
    a     = 0.1
    b     = 3.0
    basis = fundefn(:cheb,n,a,b)
    p     = funnode(basis)[1]   # collocation points

    c0 = ones(n)*0.3
    function resid!(c::Vector,result::Vector,p,basis,alpha,eta)
        # your turn!
        q = funeval(c,basis,p)[1]
        q2 = zeros(q)
        for i in eachindex(q2)
            if q[i] < 0
                q2[i] = -20.0
            else
                q2[i] = sqrt(q[i])
            end
        end
        result[:] = p.+ q .*((-1/eta)*p.^(eta+1)) .- alpha*q2 .- q.^2
    end
    f_closure(r::Vector,x::Vector) = resid!(x,r,p,basis,alpha,eta)
    res = nlsolve(f_closure,c0)
    println(res)

    # plot residual function
    x = collect(linspace(a,b,501))
    y = similar(x)
    resid!(res.zero,y,x,basis,alpha,eta);
    y = funeval(res.zero,basis,x)[1]
    pl = Any[]
    push!(pl,plot(x,y,title="residual function"))

    # plot supply functions at levels 1,10,20

    # plot demand function

    y = funeval(res.zero,basis,x)[1]
    p2 = plot(y,x,label="supply 1")
    plot!(10*y,x,label="supply 10")
    plot!(20*y,x,label="supply 20")
    d = x.^(-eta)
    plot!(d,x,label="Demand")

    push!(pl,p2)

    plot(pl...,layout=2)
end

end # module

max_stars_repo_name: snowdj/CoursePack

max_issues_repo_name: snowdj/CoursePack

max_forks_repo_name: snowdj/CoursePack


text:
using Archimedes
using Test

using Unitful: m, kg

@testset "PointMass" begin
  pm = PointMass(1m, 2m, 3kg, 0.5m)
  @test mass(pm) == 3kg
  @test getx(pm) == 1m
  @test gety(pm) == 2m
  @test radius(pm) == 0.5m

  N = 4
  masses = PointMass.(range(1.0m, stop=4.0m, length=N), fill(0m, N), fill(1kg, N), fill(1m, N))
  @test bbox(masses) == (1.0m, 4.0m, 0.0m, 0.0m)
  @test centroid(masses) == Point2D(2.5m, 0.0m)
  @test gravitycenter(masses) == PointMass(2.5m, 0.0m, 4.0kg, 2.0m)

  image_width = 300
  mapping = Archimedes.CoordMapping(image_width, bbox(masses))
  luxor_pt = Archimedes.remap(masses[1], mapping)
  @test luxor_pt.x == -150
  @test luxor_pt.y == 0
end

@testset "BoxShip" begin
  bs = BoxShip(5m, 3m, 2m, 1.5m, 0)
  sc = corners(bs)
  @test sc[1] == Point2D(-2.5m, -2.0m)
  @test sc[3] == Point2D(2.5m, 1.0m)

  triag = Point2D.([0.0m, 1.0m, 0.0m], [0.0m, 0.0m, 2.0m])
  @test Archimedes.area(triag) == 1.0m^2
end


text:
@testset "Maps and Reductions" begin
    @testset "sum" begin
        sizes = (3, 4, 7)
        @testset "dims = $dims" for dims in (:, 1)
            @testset "Array{$N, $T}" for N in eachindex(sizes), T in (Float64, ComplexF64)
                x = randn(T, sizes[1:N]...)
                test_frule(sum, x; fkwargs=(;dims=dims))
                test_rrule(sum, x; fkwargs=(;dims=dims))
            end
        end
    end

    @testset "sum abs2" begin
        sizes = (3, 4, 7)
        @testset "dims = $dims" for dims in (:, 1)
            @testset "Array{$N, $T}" for N in eachindex(sizes), T in (Float64, ComplexF64)
                x = randn(T, sizes[1:N]...)
                test_frule(sum, abs2, x; fkwargs=(;dims=dims))
                test_rrule(sum, abs2 ⊢ NoTangent(), x; fkwargs=(;dims=dims))
            end
        end
    end  # sum abs2

    @testset "sum(f, xs)" begin
        # This calls back into AD
        test_rrule(sum, abs, [-4.0, 2.0, 2.0])
        test_rrule(sum, Multiplier(2.0), [2.0, 4.0, 8.0])

        test_rrule(sum, sum, [[2.0, 4.0], [4.0,1.9]])  # array of arrays

        # dims kwarg
        test_rrule(sum, abs, [-2.0 4.0; 5.0 1.9]; fkwargs=(;dims=1))
        test_rrule(sum, abs, [-2.0 4.0; 5.0 1.9]; fkwargs=(;dims=2))

        test_rrule(sum, abs, @SVector[1.0, -3.0])

        # covectors
        x = [-4.0 2.0; 2.0 -1.0]
        test_rrule(sum, inv, x[1, :]')
        test_rrule(sum, inv, x[1:1, :]')
        test_rrule(sum, inv, transpose(view(x, 1, :)))

        # Make sure we preserve type for StaticArrays
        ADviaRuleConfig = ChainRulesTestUtils.ADviaRuleConfig
        _, pb = rrule(ADviaRuleConfig(), sum, abs, @SVector[1.0, -3.0])
        @test pb(1.0) isa Tuple{NoTangent, NoTangent, SVector{2, Float64}}

        # For structured sparse matrixes we screw it up, getting dense back
        # see https://github.com/JuliaDiff/ChainRules.jl/issues/232 etc
        _, pb = rrule(ADviaRuleConfig(), sum, abs, Diagonal([1.0, -3.0]))
        @test_broken pb(1.0)[3] isa Diagonal
    end

    @testset "prod" begin
        @testset "Array{$T}" for T in [Float64, ComplexF64]
            @testset "size = $sz, dims = $dims" for (sz, dims) in [
                ((12,), :), ((12,), 1),
                ((3,4), 1), ((3,4), 2), ((3,4), :), ((3,4), [1,2]),
                ((3,4,1), 1), ((3,2,2), 3), ((3,2,2), 2:3),
                ]
                x = randn(T, sz)
                test_rrule(prod, x; fkwargs=(dims=dims,), check_inferred=true)
                x[1] = 0
                test_rrule(prod, x; fkwargs=(dims=dims,), check_inferred=true)
                x[5] = 0
                test_rrule(prod, x; fkwargs=(dims=dims,), check_inferred=true)
                x[3] = x[7] = 0  # two zeros along some slice, for any dims
                test_rrule(prod, x; fkwargs=(dims=dims,), check_inferred=true)

                if ndims(x) == 3
                    xp = PermutedDimsArray(x, (3,2,1))  # not a StridedArray
                    xpdot, xpbar = permutedims(rand(T, sz), (3,2,1)), permutedims(rand(T, sz), (3,2,1))
                    test_rrule(prod, xp ⊢ xpbar; fkwargs=(dims=dims,), check_inferred=true)
                end
            end

            @testset "structured wrappers" begin
                # Adjoint -- like PermutedDimsArray this may actually be used
                xa = adjoint(rand(T,4,4))
                test_rrule(prod, xa ⊢ rand(T,4,4))
                test_rrule(prod, xa ⊢ rand(T,4,4), fkwargs=(dims=2,))
                @test unthunk(rrule(prod, adjoint(rand(T,3,3)))[2](1.0)[2]) isa Matrix
                @test unthunk(rrule(prod, adjoint(rand(T,3,3)), dims=1)[2](ones(1,3))[2]) isa Matrix

                # Diagonal -- a stupid thing to do, product of zeros! Shouldn't be an error though:
                @test iszero(unthunk(rrule(prod, Diagonal(rand(T,3)))[2](1.0)[2]))
                @test iszero(unthunk(rrule(prod, Diagonal(rand(T,3)), dims=1)[2](ones(1,3))[2]))
                # does a division for the complex case, so is not necessarily exact
                @test isapprox(
                    unthunk(rrule(prod, Diagonal(rand(T,1)))[2](1.0)[2]), # 1x1 sparse matrix
                    hcat(1);
                    rtol=T <: Complex ? 2eps() : 0.0,
                )
                @test unthunk(rrule(prod, Diagonal(ones(T,2)), dims=1)[2](ones(1,2))[2]) == [0 1; 1 0]

                # Triangular -- almost equally stupud
                @test iszero(unthunk(rrule(prod, UpperTriangular(rand(T,3,3)))[2](1.0)[2]))
                @test unthunk(rrule(prod, UpperTriangular(ones(T,2,2)))[2](1.0)[2]) == [0 0; 1 0]

                # Symmetric -- at least this doesn't have zeros, still an unlikely combination
                xs = Symmetric(rand(T,4,4))
                @test_skip test_rrule(prod, xs ⊢ rand(T,4,4))
                @test_skip test_rrule(prod, xs ⊢ rand(T,4,4), fkwargs=(dims=2,))
                @test unthunk(rrule(prod, Symmetric(T[1 2; -333 4]))[2](1.0)[2]) == [16 8; 8 4]
            end
        end
        @testset "Array{Float32}, no zero entries" begin
            v = [1f-5, 1f-10, 1f-15, 1f-20]
            @test prod(v) == 0
            @test unthunk(rrule(prod, v)[2](1f0)[2]) == zeros(4)
            test_rrule(prod, v)
        end
    end # prod
end


text:
module LoopVectorization

if (!isnothing(get(ENV, "TRAVIS_BRANCH", nothing)) || !isnothing(get(ENV, "APPVEYOR", nothing))) && isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
    @eval Base.Experimental.@optlevel 1
end

using VectorizationBase, SIMDPirates, SLEEFPirates, UnPack, OffsetArrays
using VectorizationBase: REGISTER_SIZE, extract_data, num_vector_load_expr,
    mask, masktable, pick_vector_width_val, valmul, valrem, valmuladd, valmulsub, valadd, valsub, _MM,
    maybestaticlength, maybestaticsize, staticm1, staticp1, staticmul, subsetview, vzero, stridedpointer_for_broadcast,
    Static, Zero, StaticUnitRange, StaticLowerUnitRange, StaticUpperUnitRange, unwrap, maybestaticrange,
    AbstractColumnMajorStridedPointer, AbstractRowMajorStridedPointer, AbstractSparseStridedPointer, AbstractStaticStridedPointer,
    PackedStridedPointer, SparseStridedPointer, RowMajorStridedPointer, StaticStridedPointer, StaticStridedStruct, offsetprecalc,
    maybestaticfirst, maybestaticlast, scalar_less, scalar_greater, noalias!, gesp, gepbyte, pointerforcomparison, NativeTypes, staticmul, staticmuladd
using SIMDPirates: VECTOR_SYMBOLS, evadd, evsub, evmul, evfdiv, vrange,
    reduced_add, reduced_prod, reduce_to_add, reduced_max, reduced_min, vsum, vprod, vmaximum, vminimum,
    sizeequivalentfloat, sizeequivalentint, vadd!, vsub!, vmul!, vfdiv!, vfmadd!, vfnmadd!, vfmsub!, vfnmsub!,
    vfmadd231, vfmsub231, vfnmadd231, vfnmsub231, sizeequivalentfloat, sizeequivalentint, #prefetch,
    vmullog2, vmullog10, vdivlog2, vdivlog10, vmullog2add!, vmullog10add!, vdivlog2add!, vdivlog10add!, vfmaddaddone, vadd1, relu
using SLEEFPirates: pow
using Base.Broadcast: Broadcasted, DefaultArrayStyle
using LinearAlgebra: Adjoint, Transpose
using Base.Meta: isexpr
using DocStringExtensions
import LinearAlgebra # for check_args

using Base.FastMath: add_fast, sub_fast, mul_fast, div_fast

export LowDimArray, stridedpointer,
    @avx, @_avx, *ˡ, _avx_!,
    vmap, vmap!, vmapt, vmapt!, vmapnt, vmapnt!, vmapntt, vmapntt!,
    vfilter, vfilter!, vmapreduce, vreduce

const VECTORWIDTHSYMBOL, ELTYPESYMBOL = Symbol("##Wvecwidth##"), Symbol("##Tloopeltype##")

"""
REGISTER_COUNT defined in VectorizationBase is supposed to correspond to the actual number of floating point registers on the system.
It is hardcoded into a file at build time.
However, someone may have multiple builds of Julia on the same system, some 32-bit and some 64-bit (e.g., they use 64-bit primarilly,
but keep a 32-bit build on hand to debug test failures on Appveyor's 32-bit build). Thus, we don't want REGISTER_COUNT to be hardcoded
in such a fashion.
32-bit builds are limited to only 8 floating point registers, so we take care of that here.

If you want good performance, DO NOT use a 32-bit build of Julia if you don't have to.
"""
const REGISTER_COUNT = Sys.ARCH === :i686 ? 8 : VectorizationBase.REGISTER_COUNT

include("getconstindexes.jl")
include("vectorizationbase_extensions.jl")
include("predicates.jl")
include("map.jl")
include("filter.jl")
include("costs.jl")
include("operations.jl")
include("graphs.jl")
include("operation_evaluation_order.jl")
include("memory_ops_common.jl")
include("add_loads.jl")
include("add_stores.jl")
include("add_compute.jl")
include("add_constants.jl")
include("add_ifelse.jl")
include("determinestrategy.jl")
include("loopstartstopmanager.jl")
include("lower_compute.jl")
include("lower_constant.jl")
include("lower_memory_common.jl")
include("lower_load.jl")
include("lower_store.jl")
include("lowering.jl")
include("split_loops.jl")
include("condense_loopset.jl")
include("reconstruct_loopset.jl")
include("constructors.jl")
include("user_api_conveniences.jl")
include("mapreduce.jl")
include("broadcast.jl")

"""
LoopVectorization provides macros and functions that combine SIMD vectorization and
loop-reordering so as to improve performance:

- [`@avx`](@ref): transform `for`-loops and broadcasting
- [`@_avx`](@ref): similar to `@avx` but does not use type information
- [`vmap`](@ref) and `vmap!`: vectorized version of `map` and `map!`
- [`vmapnt`](@ref) and `vmapnt!`: non-temporal variants of `vmap` and `vmap!`
- [`vmapntt`](@ref) and `vmapntt!`: threaded variants of `vmapnt` and `vmapnt!`
- [`vfilter`](@ref) and `vfilter!`: vectorized versions of `filter` and `filter!`
"""
LoopVectorization

include("precompile.jl")
_precompile_()

end # module


text:
using TemporalGPs: build_lgssm

_logistic(x) = 1 / (1 + exp(-x))

# Everything is tested once the LGSSM is constructed, so it is sufficient just to ensure
# that Zygote can handle construction.
function _construction_tester(f_naive::GP, storage::StorageType, σ², t::AbstractVector)
    f = to_sde(f_naive, storage)
    fx = f(t, σ²...)
    return build_lgssm(fx)
end

println("lti_sde:")
@testset "lti_sde" begin

    @testset "blk_diag" begin
        A = randn(2, 2)
        B = randn(3, 3)
        adjoint_test(TemporalGPs.blk_diag, (A, B))
        adjoint_test(TemporalGPs.blk_diag, (SMatrix{2, 2}(A), SMatrix{3, 3}(B)))
    end

    @testset "SimpleKernel parameter types" begin

        storages = (
            (name="dense storage Float64", val=ArrayStorage(Float64)),
            (name="static storage Float64", val=SArrayStorage(Float64)),
            # (name="dense storage Float32", val=ArrayStorage(Float32)),
            # (name="static storage Float32", val=SArrayStorage(Float32)),
        )

        kernels = [Matern12Kernel(), Matern32Kernel(), Matern52Kernel(), ConstantKernel(c=1.5)]

        @testset "$kernel, $(storage.name)" for kernel in kernels, storage in storages
            F, q, H = TemporalGPs.to_sde(kernel, storage.val)
            @test is_of_storage_type(F, storage.val)
            @test is_of_storage_type(q, storage.val)
            @test is_of_storage_type(H, storage.val)

            x = TemporalGPs.stationary_distribution(kernel, storage.val)
            @test is_of_storage_type(x, storage.val)
        end
    end

    @testset "lgssm_components" begin
        rng = MersenneTwister(123456)
        N = 13

        kernels = vcat(

            # Base kernels.
            (name="base-Matern12Kernel", val=Matern12Kernel()),
            map([Matern32Kernel, Matern52Kernel]) do k
                (name="base-$k", val=k())
            end,

            # Scaled kernels.
            map([1e-1, 1.0, 10.0, 100.0]) do σ²
                (name="scaled-σ²=$σ²", val=σ² * Matern32Kernel())
            end,

            # Stretched kernels.
            map([1e-2, 0.1, 1.0, 10.0, 100.0]) do λ
                (name="stretched-λ=$λ", val=transform(Matern32Kernel(), λ))
            end,

            # Summed kernels.
            (
                name="sum-Matern12Kernel-Matern32Kernel",
                val=1.5 * transform(Matern12Kernel(), 0.1) + 0.3 * transform(Matern32Kernel(), 1.1),
            ),
        )

        # Construct a Gauss-Markov model with either dense storage or static storage.
        storages = (
            (name="dense storage Float64", val=ArrayStorage(Float64)),
            (name="static storage Float64", val=SArrayStorage(Float64)),
        )

        # Either regular spacing or irregular spacing in time.
        ts = (
            (name="irregular spacing", val=collect(RegularSpacing(0.0, 0.3, N))),
            (name="regular spacing", val=RegularSpacing(0.0, 0.3, N)),
        )

        σ²s = (
            (name="homoscedastic noise", val=(0.1, ),),
            (name="heteroscedastic noise", val=(rand(rng, N) .+ 1e-1, )),
        )

        @testset "$(kernel.name), $(storage.name), $(t.name), $(σ².name)" for
            kernel in kernels,
            storage in storages,
            t in ts,
            σ² in σ²s

            println("$(kernel.name), $(storage.name), $(t.name), $(σ².name)")

            # Construct Gauss-Markov model.
            f_naive = GP(kernel.val)
            fx_naive = f_naive(collect(t.val), σ².val...)

            f = to_sde(f_naive, storage.val)
            fx = f(t.val, σ².val...)
            model = build_lgssm(fx)

            # is_of_storage_type(fx, storage.val)
            validate_dims(model)

            y = rand(rng, fx)

            @testset "prior" begin
                @test mean(fx) ≈ mean(fx_naive)
                @test cov(fx) ≈ cov(fx_naive)
                @test mean.(marginals(fx)) ≈ mean.(marginals(fx_naive))
                @test std.(marginals(fx)) ≈ std.(marginals(fx_naive))
                m_and_v = mean_and_var(fx)
                @test first(m_and_v) ≈ mean(fx)
                @test last(m_and_v) ≈ var(fx)
                @test logpdf(fx, y) ≈ logpdf(fx_naive, y)
            end

            @testset "check args to_vec properly" begin
                k_vec, k_from_vec = to_vec(kernel.val)
                @test typeof(k_from_vec(k_vec)) == typeof(kernel.val)

                storage_vec, storage_from_vec = to_vec(storage.val)
                @test typeof(storage_from_vec(storage_vec)) == typeof(storage.val)

                σ²_vec, σ²_from_vec = to_vec(σ².val)
                @test typeof(σ²_from_vec(σ²_vec)) == typeof(σ².val)

                t_vec, t_from_vec = to_vec(t.val)
                @test typeof(t_from_vec(t_vec)) == typeof(t.val)
            end

            # Just need to ensure we can differentiate through construction properly.
            adjoint_test(
                _construction_tester, (f_naive, storage.val, σ².val, t.val);
                check_infers=false, rtol=1e-6, atol=1e-6,
            )
        end
    end
end


text:
export removeExtension, gunzip

"Remove extension from file path"
removeExtension(path::String) = first(Base.Filesystem.splitext(path))

"""
Check if file is gzipped, in which case extract it. Returns original path.
`path` must not include .gz extension. Does not delete gzipped file.
"""
function gunzip(path::String)
    @assert isfile(path) || isfile(path * ".gz")
    if !isfile(path)
        run(`gunzip --keep $path.gz`)
    end
    return path
end


text:
using FileIO
using WAV
using Base.Test

testdir = dirname(@__FILE__)

# Start by loading our ground truth data from the .wav
# This data generated as sin(linspace(0,4410*2*pi,44100)[1:441])
check_data, check_fs = wavread(joinpath(testdir, "4410hz.wav"))

# Load flac'ed version of 4410hz.wav (generated with `ffmpeg -i 4410hz.wav -acodec flac 4410hz.flac`)
data, fs = load(joinpath(testdir, "4410hz.flac"))
@test fs == check_fs
@test size(data) == (441,1)
@test maximum(abs(check_data - data)) < 1e-6

# Test against 16-bit FLAC file (generated with `ffmpeg -i 4410hz.wav -acodec flac -sample_fmt s16 4410hz_s16.flac`)
s16_data, s16_fs = load(joinpath(testdir, "4410hz_s16.flac"))
@test s16_fs == check_fs
@test size(s16_data) == (441,1)
@test maximum(abs(check_data - s16_data)) < 1e-4   # reduced precision due to s16 format

# Test multichannel FLAC file (left and right channels opposing polarity versions of 4410hz signal)
stereo_data = load(joinpath(testdir, "stereo.flac"))[1]
stereo_check_data = wavread(joinpath(testdir, "stereo.wav"))[1]
@test size(stereo_data) == size(stereo_check_data)
@test maximum(abs(stereo_data - stereo_check_data)) < 1e-6

# Now that we have confidence our decoder works, let's roundtrip a signal multiple times
function roundtrip(signal, samplerate; params...)
    path = joinpath(testdir,"roundtrip.flac")
    if isfile(path)
        rm(path)
    end
    save(path, signal, samplerate; params...)
    recon_signal, recon_fs = load(path)

    @test recon_fs == samplerate
    @test size(recon_signal) == size(signal)
    #show(recon_signal - signal)
    @test maximum(abs(recon_signal - signal)) < 1e-4
    rm(path)
end

test_signal = sin(linspace(0,4410*2*pi,44100)[1:44100])''*.999
roundtrip(test_signal, 44100)
roundtrip(test_signal, 44100, bits_per_sample=16)
roundtrip(test_signal, 48000, compression_level=8)


text:
using Flux
using Flux: @functor
using MacroTools: @forward, @capture

macro toNd(ex, outref::Int=1)
    if @capture ex f_(xs__; kw__)
        kwe = Expr(:parameters, kw...)
    else
        @capture ex f_(xs__)
    end
    rxs = map(xs) do x
        :(reshape($x, size($x, 1), :))
    end
    newcall = kw === nothing ? Expr(:call, f, rxs...) : Expr(:call, f, kwe, rxs...)
    :(reshape($newcall, :, Base.tail(size($(xs[outref])))...)) |> esc
end

abstract type AbstractAttention end

create_atten_mask(T::Type, score::AbstractArray, ::Nothing, future) = create_atten_mask(T, score, fill!(similar(score, size(score,1), size(score, 2), 1), one(T)), future)
function create_atten_mask(T::Type, score::AbstractArray, _mask::AbstractArray, future::Bool=false)
  #size(mask) == (q, k, n, b)

  # ql, kl = size(mask)
  mask = copy(_mask)

  maskval = convert(T, -1e9)
  !future && batched_triu!(mask, 0)
  mask .= (1 .- mask) .* maskval
  return mask
end

Flux.@nograd create_atten_mask

struct MultiheadAttention{Q<:Dense, K<:Dense, V<:Dense, O<:Dense, DP<:Dropout} <: AbstractAttention
    head::Int
    future::Bool
    iqproj::Q
    ikproj::K
    ivproj::V
    oproj::O
    drop::DP
end

Flux.functor(mh::MultiheadAttention) = (mh.iqproj, mh.ikproj, mh.ivproj, mh.oproj), m -> MultiheadAttention(mh.head, mh.future, m..., mh.drop)

"""
    MultiheadAttention(head::Int, is::Int, hs::Int, os::Int;
                       future::Bool=true, pdrop = 0.1)

Multihead dot product Attention Layer, `head` is the number of head,
`is` is the input size, `hs` is the hidden size of input projection layer of each head,
`os` is the output size. When `future` is `false`, the k-th token can't see tokens at > k.
`pdrop` is the dropout rate.
"""
MultiheadAttention(head::Int,
                   is::Int,
                   hs::Int,
                   os::Int;
                   future::Bool=true, pdrop = 0.1) = MultiheadAttention(head,
                                                                        future,
                                                                        Dense(is, hs*head),
                                                                        Dense(is, hs*head),
                                                                        Dense(is, hs*head),
                                                                        Dense(hs*head, os),
                                                                        Dropout(pdrop),
                                                                        )

function Base.show(io::IO, mh::MultiheadAttention)
    hs = div(size(mh.iqproj.W)[1], mh.head)
    is = size(mh.iqproj.W)[end]
    os = size(mh.oproj.W)[1]

    print(io, "MultiheadAttention(")
    print(io, "head=$(mh.head), ")
    print(io, "head_size=$(hs), ")
    print(io, "$(is)=>$(os)")

    if Flux.istraining()
        print(io, ", dropout=$(mh.drop.p))")
    else
        print(io, ")")
    end
end

function (mh::MultiheadAttention)(query::A1,
                                  key::A2,
                                  value::A3;
                                  mask=nothing) where {T,
                                                       A1 <: Abstract3DTensor{T},
                                                       A2 <: Abstract3DTensor{T},
                                                       A3 <: Abstract3DTensor{T}}
  qs = size(query)
  ks = size(key)
  vs = size(value)

  #size(ipq) == (h, q_seq_len, batch)
  ipq = @toNd mh.iqproj(query)
  ipk = @toNd mh.ikproj(key)
  ipv = @toNd mh.ivproj(value)

  h = size(ipq, 1)
  hs = div(h, mh.head)

  #size(ipq) == (hs, q_seq_len, head, batch)
  ipq = permutedims(reshape(ipq, hs, mh.head, qs[2], qs[3]), [1, 3, 2, 4])
  ipk = permutedims(reshape(ipk, hs, mh.head, ks[2], ks[3]), [1, 3, 2, 4])
  ipv = permutedims(reshape(ipv, hs, mh.head, vs[2], vs[3]), [1, 3, 2, 4])

  #size(ipq) == (hs, q_seq_len, head * batch)
  ipq = reshape(ipq, hs, qs[2], :)
  ipk = reshape(ipk, hs, ks[2], :)
  ipv = reshape(ipv, hs, vs[2], :)

  atten = attention(ipq,ipk,ipv,
                    mask,
                    mh.future,
                    mh.drop)

  atten = permutedims(reshape(atten, hs, qs[2], mh.head, qs[3]), [1, 3, 2, 4]) #size(atten) == (hs, head, ql, b)
  atten = reshape(atten, h, qs[2], qs[3]) #size(atten) == (h, ql, b)

  out = @toNd mh.oproj(atten)
  out #size(out) == (h, q_seq_len, batch)
end

function (mh::MultiheadAttention)(query::A1,
                                  key::A2,
                                  value::A3;
                                  mask=nothing) where {T,
                                                       A1 <: AbstractMatrix{T},
                                                       A2 <: AbstractMatrix{T},
                                                       A3 <: AbstractMatrix{T}}

  # size(query) == (dims, seq_len)
  ipq = mh.iqproj(query)
  ipk = mh.ikproj(key)
  ipv = mh.ivproj(value)

  h = size(ipq)[1] #h == hs * head
  hs = div(h, mh.head)

  #size(hq) == (hs, seq_len, head)
  hq = permutedims(reshape(ipq, hs, mh.head, :), [1, 3, 2])
  hk = permutedims(reshape(ipk, hs, mh.head, :), [1, 3, 2])
  hv = permutedims(reshape(ipv, hs, mh.head, :), [1, 3, 2])

  atten = attention(hq, hk, hv,
                    mask,
                    mh.future,
                    mh.drop)

  # size(atten) == (head*hs, seq_len)
  atten = reshape(permutedims(atten, [1, 3, 2]), h, :)

  mh.oproj(atten)
end

function apply_mask(score, mask)
  s = size(score)
  ms = size(mask)
  bxn = s[end]
  b = ms[end]
  if bxn == b || b == 1
    return score .+ mask
  else
    return reshape(reshape(score, s[1:end-1]..., :, b) .+
                   reshape(mask, ms[1:end-1]..., 1, b), s)
  end
end

apply_mask(score::AbstractArray{T}, ::Nothing, future) where T = future ? score : apply_mask(score, create_atten_mask(T, score, nothing, future))
apply_mask(score::AbstractArray{T}, mask, future) where T = apply_mask(score, create_atten_mask(T, score, mask, future))

function attention(query::A1,
                   key::A2,
                   value::A3,
                   mask, future::Bool,
                   dropout) where {T,
                                   A1 <: Abstract3DTensor{T},
                                   A2 <: Abstract3DTensor{T},
                                   A3 <: Abstract3DTensor{T}}
  #size(query) == (dims, {q,k}_seq_len, batch) == size(key) == size(value)
  #size(score) == (k_seq_len, q_seq_len, batch)
  dk = size(key, 1)
  score = batchedmul(key, query; transA = true)
  score = score ./ convert(T, sqrt(dk))

  score = apply_mask(score, mask, future)
  score = softmax(score; dims=1)
  dropout !== nothing && (score = dropout(score))
  batchedmul(value, score) #size(return) == (dims, q_seq_len, batch)
end


text:
module GrahamScan

export FindConvexHull, EdgePoint, SortAngle

using LinearAlgebra

function EdgePoint(A)
    # find lexicographic minimum row of A with the second column having
    # preceeding importance over the first one.
    min_y = minimum(A[:,2]);
    candidates = A[:,2] .== min_y;
    min_x = minimum(A[candidates, 1]);
    [min_x min_y]
end

function SortAngle(points)
    function Cos(x,y)
        c = acos(x/sqrt(x^2 + y^2));
        if x > 0 && y > 0
            c
        elseif x < 0 && y > 0
            c + π/2
        elseif x > 0 && y < 0
            -c
        elseif x < 0 && y < 0
            -c - π/2
        else
            c
        end
    end
    f(p1,p2) = Cos(p1...) < Cos(p2...);
    sortslices(points; dims=1, lt=f)
end

function ToDirection(p1, p2, p3)
    T = [1 p1;
         1 p2;
         1 p3];
    d = det(T);
    if d < 0
        :right
    elseif d == 0
        :inline
    else
        :left
    end
end

function Scan(p1, ps)
    (number_of_points, _) = size(ps);
    if number_of_points < 2
        [p1; ps]
    else
        p2 = ps[1,:]';
        p3 = ps[2,:]';
        ps = ps[3:end,:];
        direction = ToDirection(p1, p2, p3);
        if direction == :right
            Scan(p1, [p3; ps]) # reject move, again with same tail and next head
        elseif direction == :left
            [p1; Scan(p2, [p3; ps])] # Save point, onward with next point
        elseif direction == :inline
            Scan(p1, [p3; ps]) # Discard mid point, onward with next point
        end
    end
end

function FindConvexHull(P)
    # P is a (number of points) x (2) matrix
    (number_of_points, _) = size(P);
    if number_of_points == 0
        return []
    else
        start = EdgePoint(P);
        start_idx = (P[:, 1] .== start[1]) .& (P[:, 2] .== start[2]);
        without_start = reshape(P[.~[start_idx start_idx]], :, 2);
        # transform p into a coordinate system where o is
        # the origin (relative to the previous system)
        TranslateTo(o, p) = p - o;
        TranslateCloud(o, cloud) = hcat(map(p -> TranslateTo(o, p), eachrow(cloud))...)';
        sorted_points = TranslateCloud(start', without_start) |>
            SortAngle |>
            c -> TranslateCloud(-start', c)
        if number_of_points <= 3
            return [start; sorted_points]
        else
            return Scan(start, sorted_points)
        end
    end
end

#=
Test 1
A = [1 -1; 3 4; -1 2; -1 -1]

GrahamScan(A) # should be A

Test 2
B = [2.3486669219560348 2.7215254125789414;
    2.4139648037583257  -4.681676398649306;
    -0.6260514711523584 6.850334341587001;
    -7.64216704719866   -2.7605815091225576;
    -5.026012253905725  1.5777744917255312;
    -6.4853594504513445 5.318092139833023;
    5.082500986696756   -1.896255262210346;
    -6.497957092016469  -6.549163826826357;
    7.051754797682392   -5.031079040862534;
    3.283138295813238   -8.076020743414723;
    7.805413024480348   9.26931407569144;
    7.910312138709568   5.417744228015254;
    2.4605192481055127  -1.4160149654176593;
    -2.877588474503572  -3.4053645090358486;
    0.5611514305387537  -9.743920605133937;
    8.555954343029967   -4.394948528913232;
    1.550287503925139   9.078244403761712;
    -3.271086970349943  1.934092461779306;
    -2.2965794125239603 1.2422580281368312;
    3.182950465475866   7.148107783065644]

C = [[0.5611514305387537,-9.743920605133937],[3.283138295813238,-8.076020743414723],[8.555954343029967,-4.394948528913232],[7.805413024480348,9.26931407569144],[3.182950465475866,7.148107783065644],[1.550287503925139,9.078244403761712],[-0.6260514711523584,6.850334341587001],[-6.4853594504513445,5.318092139833023],[-5.026012253905725,1.5777744917255312],[-7.64216704719866,-2.7605815091225576],[-6.497957092016469,-6.549163826826357]]
c = GrahamScan(B)
hcat(C...)' - c # should be 0
=#

end


text:
## ----------------------------------------------------------------------------
function maxent(V, D, cgD_X;
        # Top loop
        top_maxiter = 100,

        # GD globals
        gd_maxiter = 1000,
        gdth = 1e-2,
        verbose = false,

        # stst
        stst_th = 1e-2,
        stst_w = 10,

        check_b0 = true,
        further_info = (;)

    )

    # ---------------------------------------------------------------
    # Space globals
    Vz = Vi(V, :z)
    Vug = Vi(V, :ug)
    Vvol = length(V)
    V = nothing # save mem?

    # ---------------------------------------------------------------
    # Alg globals
    PME = ones(Vvol)
    normalizeP!(PME)
    # Paux = similar(PME)
    I = eachindex(PME)
    z_avPME = 0.0
    ug_avPME = 0.0
    z_beta = 0.0
    ug_beta = 0.0
    ug_avPME_at_b0 = 0.0
    ug_valid_at_b0 = false
    z_betas = Float64[]
    ug_betas = Float64[]

    ## -----------------------------------------------------------
    # _maxent_fun!
    Uz = maximum(Vz)
    Uug = maximum(Vug)
    function _maxent_fun!(P, z_beta_, ug_beta_; min_P = 1e-300)
        norm = max(z_beta_ * Uz, ug_beta_ * Uug)
        @inbounds for i in I
            P[i] = max(min_P, exp((z_beta_ * Vz[i] + ug_beta_ * Vug[i]) - norm))
        end
        P
    end

    ## -----------------------------------------------------------
    # top loop globals
    topiter = 1
    conv = false
    thid = threadid()

    ## -----------------------------------------------------------
    # GD globals
    gd_err = 0.0
    gd_iter = 0

    ## -----------------------------------------------------------
    # beta stst
    nbuffs = length((:z_beta, :ug_beta))
    beta_stst = STST(nbuffs, stst_w)

    ## -----------------------------------------------------------
    # Info funtions
    function gd_ug_info()
        @info("ug grad descent ",
            further_info...,
            topiter, gd_iter,
            (ug_avPME, cgD_X),
            gd_err, ug_beta, thid
        )
    end

    function gd_z_info()
        @info("z grad Descent ",
            further_info...,
            topiter, gd_iter,
            (z_avPME, D),
            gd_err, z_beta, thid
        )
    end

    function finished_round_info()
        @info("Finished Round",
            further_info...,
            topiter, conv,
            (z_beta, ug_beta),
            (ug_avPME, cgD_X),
            (ug_avPME_at_b0, cgD_X),
            (z_avPME, D),
            thid
        );
    end

    ## -----------------------------------------------------------
    while true

        ## -----------------------------------------------------------
        # find z_beta
        # Gradient descent
        target = D
        x0 = z_beta
        maxΔx = max(100.0, abs(z_beta) * 0.1)
        x1 = x0 + maxΔx * 0.1
        gd_iter = 1

        function gd_z_upfun!(gdmodel)

            z_beta = SimTools.gd_value(gdmodel)

            _maxent_fun!(PME, z_beta, ug_beta)
            normalizeP!(PME)

            z_avPME = 0.0
            @inbounds for i in I
                z_avPME += PME[i] * Vz[i]
            end

            gd_iter += 1
            gd_err = gdmodel.ϵi

            return z_avPME
        end

        gdmodel = SimTools.grad_desc(gd_z_upfun!; gdth,
            target, x0, x1, maxΔx,
            maxiter = gd_maxiter, verbose = false
        )
        z_beta = SimTools.gd_value(gdmodel)

        verbose && gd_z_info()

        ## -----------------------------------------------------------
        # Check balance at ug_beta = 0.0 the first time
        if check_b0 && (topiter == 1)
            ug_beta0 = 0.0
            _maxent_fun!(PME, z_beta, ug_beta0)
            normalizeP!(PME)
            ug_avPME_at_b0 = sum(PME .* Vug)
            ug_valid_at_b0 = (abs(ug_avPME_at_b0) <= abs(cgD_X))
            @show ug_valid_at_b0
        end

        ## -----------------------------------------------------------
        # if not valid, move ug_beta
        if ug_valid_at_b0
            ug_beta = 0.0
            ug_avPME = ug_avPME_at_b0
        else
            # Gradient descent
            target = cgD_X * (1.0 - gdth)
            x0 = ug_beta
            maxΔx = max(10.0, abs(ug_beta) * 0.1)
            x1 = x0 + maxΔx * 0.1
            gd_iter = 1

            function gd_ug_upfun!(gdmodel)

                ug_beta = SimTools.gd_value(gdmodel)

                _maxent_fun!(PME, z_beta, ug_beta)
                normalizeP!(PME)

                ug_avPME = 0.0
                @inbounds for i in I
                    ug_avPME += PME[i] * Vug[i]
                end

                gd_iter += 1
                gd_err = gdmodel.ϵi

                return ug_avPME
            end

            gdmodel = SimTools.grad_desc(gd_ug_upfun!; gdth,
                target, x0, x1, maxΔx,
                maxiter = gd_maxiter, verbose = false
            )
            ug_beta = SimTools.gd_value(gdmodel)
        end

        verbose && gd_ug_info()

        push!(z_betas, z_beta); push!(ug_betas, ug_beta)
        push!(beta_stst, z_beta, ug_beta)

        conv = ug_valid_at_b0 || is_steady(beta_stst, stst_th)

        verbose && finished_round_info()

        conv && break

        topiter += 1
        topiter > top_maxiter && break
    end # while true

    status = conv ? :conv : :maxiter
    return (;PME, z_beta, ug_beta, status)

end

match_moms(V, D, cgD_X; kwargs...) = maxent(V, D, cgD_X; kwargs..., check_b0 = false)


text:
"""
Subroutine daxpy

This subroutine computes constant times a vector plus a vector.
It uses unrolled loops for increments equal to one.
Jack Dongarra, LINPACK, 3/11/78.
"""
function tron_daxpy(n,da,dx,incx,dy,incy)
    if n <= 0
        return
    end
    if da == 0.0
        return
    end
    if incx != 1 || incy != 1
        ix = 1
        iy = 1
        if incx < 0
            ix = (-n+1)*incx + 1
        end
        if incy < 0
            iy = (-n+1)*incy + 1
        end
        for i=1:n
            dy[iy] += da*dx[ix]
            ix += incx
            iy += incy
        end
        return
    else

        # Code for both increments equal to 1

        m = mod(n, 4)
        if m != 0
            for i=1:m
                dy[i] += da*dx[i]
            end
        end
        if n < 4
            return
        end

        mp1 = m + 1
        for i=mp1:4:n
            dy[i] += da*dx[i]
            dy[i + 1] += da*dx[i + 1]
            dy[i + 2] += da*dx[i + 2]
            dy[i + 3] += da*dx[i + 3]
        end
        return
    end
end

daxpy(n,da,dx,incx,dy,incy) = tron_daxpy(n,da,dx,incx,dy,incy)


text:
# CENG0038 Droplet Distribution research project

dir = "./src"
push!(LOAD_PATH,dir)

module researchProject

import trachea
import bifurcations

# straightTube(T, ro, Q, mu, ro_d, d_d, d, phi, L)
T = 309.5 # Temperature of Air
ro = 1.12027 # Density of Air (http://www.anaesthesia.med.usyd.edu.au/resources/lectures/humidity_clt/humidity.html#:~:text=While%20nose%20breathing%20at%20rest,%25%20to%2070%25%20relative%20humidity.)
Q = 0.001 # Based on 60L/min flowrate of ventilators
mu = 1.895e-5 # Dynamic viscosity of air (https://www.engineersedge.com/physics/viscosity_of_air_dynamic_and_kinematic_14483.htm#:~:text=The%20viscosity%20of%20air%20depends,10%2D5%20Pa%C2%B7s%20.&text=At%2025%20%C2%B0C%2C%20the,the%20kinematic%20viscosity%2015.7%20cSt.)
ro_d = 1000 # Droplet density of water
d = 0.0189 # Diameter of trachea (https://www.google.com/url?sa=i&url=https%3A%2F%2Fwww.slideshare.net%2Fjinojustinj%2Ftracheal-pathologies&psig=AOvVaw38sp20tGz3PGVizR_R8F_N&ust=1605869019853000&source=images&cd=vfe&ved=0CAIQjRxqFwoTCICtjbi2ju0CFQAAAAAdAAAAABAT)
phi = (pi/2) # Inclination angle in radians (phi) SI units: rad
# IBA = 1.040216 # Interbronchial angle, average for all ages, 59.6 degrees SI Units: radians
# phi = pi-(pi/2)-(IBA/2) # Inclination angle in radians (phi) SI units: rad
L = 0.118 # Avg length of trachea (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5900092/#:~:text=On%20average%2C%20the%20length%20of,to%20be%20shorter%20in%20females.&text=Structure%20of%20the%20trachea.)

## Dataset 1
# Source: "The Mechanism of Breath Formation"
dataset1 = [5.4444E-07, 5.8472E-07, 6.2656E-07, 6.6727E-07, 7.2656E-07, 7.6727E-07, 8.3492E-07, 9.0000E-07, 9.6627E-07, 1.0557E-06, 1.1621E-06, 1.2616E-06, 1.3686E-06, 1.4746E-06, 1.5741E-06, 1.6736E-06, 1.7828E-06, 1.8863E-06, 1.9858E-06, 2.1567E-06, 2.3432E-06, 2.5136E-06, 2.6939E-06, 2.8643E-06, 3.0672E-06, 3.3113E-06, 3.5555E-06, 3.8230E-06]

post_trachea = []

for d_d = dataset1
    penetration = trachea.straightTube(T, ro, Q, mu, ro_d, d_d, d, phi, L)
    push!(post_trachea, penetration)
end

print(post_trachea)

# bifurcation(ro_d, d_d, d_0, d_1, mu, Q, alpha)
ro_d = 1000 # Droplet density of water
d_0 = 0.018 # Diameter of trachea (https://www.google.com/url?sa=i&url=https%3A%2F%2Fwww.slideshare.net%2Fjinojustinj%2Ftracheal-pathologies&psig=AOvVaw38sp20tGz3PGVizR_R8F_N&ust=1605869019853000&source=images&cd=vfe&ved=0CAIQjRxqFwoTCICtjbi2ju0CFQAAAAAdAAAAABAT)
d_1 = 0.0122 # Diameter of daughter (https://www.google.com/url?sa=i&url=https%3A%2F%2Fwww.slideshare.net%2Fjinojustinj%2Ftracheal-pathologies&psig=AOvVaw38sp20tGz3PGVizR_R8F_N&ust=1605869019853000&source=images&cd=vfe&ved=0CAIQjRxqFwoTCICtjbi2ju0CFQAAAAAdAAAAABAT)
mu = 1.895e-5 # Dynamic viscosity of air (https://www.engineersedge.com/physics/viscosity_of_air_dynamic_and_kinematic_14483.htm#:~:text=The%20viscosity%20of%20air%20depends,10%2D5%20Pa%C2%B7s%20.&text=At%2025%20%C2%B0C%2C%20the,the%20kinematic%20viscosity%2015.7%20cSt.)
Q = 0.001 # Based on 60L/min flowrate of ventilators
alpha = 1.06116 # SCA, all ages. Haskin Goodman et. al. "Normal Tracheal Bifurcation Angle"

post_bifurc0 = []

for d_d = dataset1
    bifurc_penetration = bifurcations.bifurcation(ro_d, d_d, d_0, d_1, mu, Q, alpha)
    push!(post_bifurc0, bifurc_penetration)
end

# straightTube(T, ro, Q, mu, ro_d, d_d, d, phi, L)
T = 309.5 # Temperature of Air
ro = 1.12027 # Density of Air (http://www.anaesthesia.med.usyd.edu.au/resources/lectures/humidity_clt/humidity.html#:~:text=While%20nose%20breathing%20at%20rest,%25%20to%2070%25%20relative%20humidity.)
Q = 0.001 # Based on 60L/min flowrate of ventilators
mu = 1.895e-5 # Dynamic viscosity of air (https://www.engineersedge.com/physics/viscosity_of_air_dynamic_and_kinematic_14483.htm#:~:text=The%20viscosity%20of%20air%20depends,10%2D5%20Pa%C2%B7s%20.&text=At%2025%20%C2%B0C%2C%20the,the%20kinematic%20viscosity%2015.7%20cSt.)
ro_d = 1000 # Droplet density of water
d = 0.0189 # Diameter of trachea (https://www.google.com/url?sa=i&url=https%3A%2F%2Fwww.slideshare.net%2Fjinojustinj%2Ftracheal-pathologies&psig=AOvVaw38sp20tGz3PGVizR_R8F_N&ust=1605869019853000&source=images&cd=vfe&ved=0CAIQjRxqFwoTCICtjbi2ju0CFQAAAAAdAAAAABAT)
# phi = (pi/2) # Inclination angle in radians (phi) SI units: rad
IBA = 1.040216 # Interbronchial angle, average for all ages, 59.6 degrees SI Units: radians Haskin Goodman et. al. "Normal Tracheal Bifurcation Angle"
phi = pi-(pi/2)-(IBA/2) # Inclination angle in radians (phi) SI units: rad
L = 0.118 # Avg length of trachea (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5900092/#:~:text=On%20average%2C%20the%20length%20of,to%20be%20shorter%20in%20females.&text=Structure%20of%20the%20trachea.)

post_inclinedTube0 = []

for d_d = dataset1
    penetration = trachea.straightTube(T, ro, Q, mu, ro_d, d_d, d, phi, L)
    push!(post_inclinedTube0, penetration)
end

# print(post_trachea)
# print("\n")
# print(post_inclinedTube0)
# print("\n")
# print(size(post_trachea))
# print("\n")
# print("\n")
# print(post_bifurc0)
# print("\n")
# print(size(post_bifurc0))

# Deposition Efficiency after first bifurcation
n_i = post_trachea .* post_bifurc0

# print(n_i)

# Distribution for dataset1
d_i1 = [0.056718217, 0.060143555, 0.068012484, 0.071244953, 0.073582454, 0.076460071, 0.078809187, 0.077779268, 0.072560295, 0.065129111, 0.056718217, 0.045163553, 0.037531843, 0.031410246, 0.026885611, 0.02266377, 0.018430315, 0.01499732, 0.01121713, 0.00960475, 0.007514077, 0.005690533, 0.003839032, 0.002862148, 0.002179373, 0.001434928, 0.000799432, 0.000618126]
# print(size(d_i1))

# Penetration * Distribution for data set 1 (one bifurcation)
final_distribution = n_i .* d_i1

print(final_distribution)

# L = 0.1
# print(L)

# Write your package code here.

end


text:
"""
Package: Forecast

    sma(x, n; center = true)

Smooth a vector of data using a simple moving average.

# Arguments
- `x`: Vector of data.
- `n`: Size of the moving average.
- `center`: centers the moving averaged values in the response.

# Returns
Vector of moving average smoothed values containing `missing` values to preserve the size of the original vector.

# Examples
```julia-repl
julia> sma(1:5,3;center=true)
5-element Array{Any,1}:
  missing
 2.0
 3.0
 4.0
  missing
```
"""
function sma(x, n; center = true)

    N = length(x)
    @assert 1 <= n & n <= N

    if n == 1
        return x
    end

    res = Vector{Any}(missing,N)

    # initial and final value positions
    ivp = findfirst(!ismissing, x)
    fvp = findlast(!ismissing, x)

    # using missing values to center ma
    a = center ? div(n,2) : 0

    # initial moving average value
    ma = sum(x[ivp:n+ivp-1])/n
    res[a+ivp] = ma
    for i in 1:(N-n-ivp-(N-fvp)+1)
        resai = ma + (x[n+ivp+i-1] - x[ivp+i-1])/n
        # missing values are imputed
        res[a+ivp+i] = ismissing(resai) ? ma : resai
        ma = res[a+ivp+i]
    end

    res

end