text:
# Tests the code for computing the derivatives
# of M_n with respect to k.
include("../src/transit_poly_struct.jl")
include("../src/Mn_compute.jl")
include("../src/transit_structure.jl")
include("../src/cel_bulirsch.jl")
# include("random.jl")

# Randomizer seed
#using Random
seed = Random.seed!(42)

using QuadGK

@testset "Mn_derivative" begin

function Mn_num(k2::T,m::Int64) where {T <: Real}
# Numerically integrates M_n(k^2)/(4br)^m (note: m can be a half-integer).
# See 11/10/2018 notes.
f(x) = sqrt(k2-sin(x)^2)^m
if k2 < 1.0
  kap2 = convert(T,asin(sqrt(big(k2))))
  Mn,error = quadgk(f,-kap2,kap2,rtol=1e-15)
else
  pi2 = 0.5*pi
  Mn,error = quadgk(f,-pi2,pi2,rtol=1e-15)
end
return Mn
end

# Initialize all the variables needed for computing Mn:
function initialize_big(t::Transit_Struct{T}) where {T <: Real}
k2 = t.k2
t.sqonembmr2 = sqrt(1-(t.r-t.b)^2)
t.sqarea = sqarea_triangle(one(T),t.r,t.b)
t.sqbr = sqrt(t.r*t.b)
if k2 < 1
  t.kc2 = convert(T,1.0-big(k2))
  t.kc = convert(T,sqrt(1.0-big(k2)))
  t.k = convert(T,sqrt(big(k2)))
  t.kck = convert(T,sqrt((1.0-big(k2))*big(k2)))
  t.kap0 = convert(T,2*asin(sqrt(big(k2))))
  Piofk,Eofk_big,Em1mKdm_big = cel_bulirsch(big(k2),sqrt(1-big(k2)),zero(BigFloat), one(BigFloat),one(BigFloat),one(BigFloat), one(BigFloat),1-big(k2),zero(BigFloat))
  t.Eofk = convert(T,Eofk_big)
  t.Em1mKdm = convert(T,Em1mKdm_big)
else
  t.k2inv = convert(T,inv(big(k2)))
  t.kc2 = convert(T,1.0-inv(big(k2)))
  t.kc = convert(T,sqrt(1.0-inv(big(k2))))
  t.k = convert(T,sqrt(big(k2)))
  t.kap0 = convert(T,pi); kck=zero(T)
  Piofk,Eofk_big,Em1mKdm_big = cel_bulirsch(inv(big(k2)),sqrt(1-inv(big(k2))),zero(BigFloat), one(BigFloat),one(BigFloat),one(BigFloat), one(BigFloat),1-inv(big(k2)),zero(BigFloat))
  t.Eofk = convert(T,Eofk_big)
  t.Em1mKdm = convert(T,Em1mKdm_big)
end
return
end

# Initialize all the variables needed for computing Mn:
function initialize(t::Transit_Struct{T}) where {T <: Real}
k2 = t.k2
t.sqonembmr2 = sqrt(1-(t.r-t.b)^2)
t.sqarea = sqarea_triangle(one(T),t.r,t.b)
t.sqbr = sqrt(t.r*t.b)
if k2 < 1
  t.kc2 = convert(T,1.0-big(k2))
  t.kc = convert(T,sqrt(1.0-big(k2)))
  t.k = convert(T,sqrt(big(k2)))
  t.kck = convert(T,sqrt((1.0-big(k2))*big(k2)))
  t.kap0 = convert(T,2*asin(sqrt(big(k2))))
  Piofk,t.Eofk,t.Em1mKdm = cel_bulirsch(t.k2,t.kc,zero(T), one(T),one(T),one(T), one(T),t.kc2,zero(T))
else
  t.k2inv = convert(T,inv(big(k2)))
  t.kc2 = convert(T,1.0-t.k2inv)
  t.kc = convert(T,sqrt(1.0-inv(big(k2))))
  t.k = convert(T,sqrt(big(k2)))
  t.kap0 = convert(T,pi); kck=zero(T)
  Piofk,t.Eofk,t.Em1mKdm = cel_bulirsch(inv(t.k2),t.kc,zero(T), one(T),one(T),one(T), one(T),t.kc2,zero(T))
end
return
end

# Carries out a test of M_n for k^2 and m_max:
function test_Mn(k2::Float64)
@assert (k2 > 0)
m_max = 30
# Set up arrays to hold finite-difference derivatives:
# Finite difference step:
dq = big(1e-18)
# Set up transit structures for Float64 and BigFloat.
# Since we only care about k, set r=b for simplicity:
r_big = inv(sqrt(2*big(k2))); b_big = r_big
r = convert(Float64,r_big); b=r; u = ones(2*m_max+1); u_big = big.(u)
# Initialize the transit structure to pass to routines:
t = transit_init(r,b,u,true)  # Float64
t_big = transit_init(r_big,b_big,u_big,false) # BigFloat
t_bigm = transit_init(r_big,b_big,u_big,false) # BigFloat
t_bigp = transit_init(r_big,b_big,u_big,false) # BigFloat
t.k2 = k2
t_big.k2 = big(k2)
# Initialize variables needed for computing M_n:
initialize_big(t)
initialize(t_big)
# Compute the derivatives in Float64 precision:
if t.k2 < 0.5 || t.k2 > 2.0
  Mn_lower!(t)
  Mn_lower!(t_big)
else
  Mn_raise!(t)
  Mn_raise!(t_big)
end
reltol = 1e-6
abstol = 1e-15
# Compare with numerical integration:
for m=0:t.m_max
  Mnn = Mn_num(k2,m)
  @test isapprox(Mnn,t.Mn[m+1],atol=abstol,rtol = reltol)
  if !isapprox(Mnn,t.Mn[m+1],atol=abstol,rtol = reltol)
    println("m: ",m," Mn_num: ",Mnn," M_n: ",t.Mn[m+1]," diff: ",Mnn-t.Mn[m+1])
  end
  @test isapprox(t.Mn[m+1],convert(Float64,t_big.Mn[m+1]),atol=abstol,rtol=reltol)
  if !isapprox(t.Mn[m+1],convert(Float64,t_big.Mn[m+1]),atol=abstol,rtol=reltol)
    println("m: ",m," Mn: ",t.Mn[m+1]," M_n_big: ",convert(Float64,t_big.Mn[m+1])," diff: ",t.Mn[m+1]-convert(Float64,t_big.Mn[m+1]))
  end
end
return
end

eps = 1e-8
r=100.0; b0=[r-1+eps,r,r+1-eps]
for i=1:length(b0)
  b=b0[i]
  k2 = (1-b+r)*(1+b-r)/(4*b*r)
  if k2 >= 0
    test_Mn(k2)
  end
end
r=0.01; b0=[eps,r,1-r-eps,1-r+eps,1.0,r+1-eps]
for i=1:length(b0)
  b=b0[i]
  k2 = (1-b+r)*(1+b-r)/(4*b*r)
  if k2 >= 0
    test_Mn(k2)
  end
end
ntest = 10
for i=1:ntest
  test_Mn(.5*rand(seed))
  test_Mn(.5+.5*rand(seed))
  test_Mn(inv(.5+.5*rand(seed)))
  test_Mn(inv(.5*rand(seed)))
end
end


text:
module KMedoids

    function kmedoids(k::Int64, data::Matrix)
        return "This function is not yet implemented."
    end

end

ext: jl


text:
module Munkres

export munkres

immutable Location
    row::Int
    col::Int
end

Location() = Location(-1,-1)

valid(loc::Location) = loc.row > 0 && loc.col > 0

const StarMark=1
const PrimeMark=2

type ModifiedCost
    m::Matrix{Float64}
    row_offsets::Vector{Float64}
    column_offsets::Vector{Float64}
end

function ModifiedCost{T<:Real}(cost_matrix::Matrix{T})
    ModifiedCost(cost_matrix,
                 zeros(eltype(cost_matrix), size(cost_matrix,1)),
                 zeros(eltype(cost_matrix), size(cost_matrix,2)))
end

Base.size(cost::ModifiedCost, args...) = size(cost.m, args...)

Base.eltype(cost::ModifiedCost) = eltype(cost.m)

function iszero(cost::ModifiedCost, i, j)
    @inbounds return cost.m[i,j] == cost.row_offsets[i] + cost.column_offsets[j]
end

function Base.getindex(cost::ModifiedCost, i, j)
    @inbounds return cost.m[i,j] - cost.row_offsets[i] - cost.column_offsets[j]
end

"""
Munkres (Hungarian) Algorithm for optimal assignment, i.e.  given an NxM cost
matrix for assigning N workers to M jobs, what is the optimal allocation of
tasks to workers such that one one task is assigned to one worker.

Input

  * cost_matrix - an NxM real valued matrix of the cost of assigning job m to worker N

Output

  * p - A vector of length N of assignments, where j = p[i] assigns job j to worker i.

"""
function munkres(cost_matrix)
    # Inspired by Bob Pilgrim's C# tutorial, "Munkres' Assignment Algorithm,
    # Modified for Rectangular Matrices",
    # http://csclab.murraystate.edu/bob.pilgrim/445/munkres.html
    # and by and Yi Cao's surprisingly fast matlab version,
    # http://mathworks.com/matlabcentral/fileexchange/20328-munkres-assignment-algorithm

    n,m = size(cost_matrix)
    flipped = false
    if n > m
        #always use more jobs than workers (current implmentation doesn't work in other case, so just transpose)
        cost_matrix = cost_matrix'
        flipped = true
        n,m = size(cost_matrix)
    end
    row_cover = zeros(Bool,n)
    column_cover = zeros(Bool,m)
    mask_array = zeros(Int8,n,m)
    path_start = Location()
    cost = ModifiedCost(cost_matrix)

    zero_locations = step_one!(cost)
    step_two!(cost, mask_array, row_cover, column_cover)

    step = 3

    while true
        #println("step = ",step)
        if step == 3
            step = step_three!(mask_array, column_cover)
        elseif step == 4
            step, path_start = step_four!(cost, mask_array, row_cover, column_cover, zero_locations)
        elseif step == 5
            step = step_five!(mask_array, row_cover, column_cover, path_start)
            path_start = Location()
        elseif step == 6
            step = step_six!(cost,row_cover,column_cover, zero_locations)
        elseif step == 7
            break
        else
            error("Step has gone bonkers")
        end
        #awesome_print_debug(cost, mask_array, row_cover, column_cover, step, path_start)
    end

    if flipped
        return [findfirst(mask_array[:,i] .== StarMark) for i=1:size(mask_array,2)]
    end
    return [findfirst(mask_array[i,:] .== StarMark) for i=1:size(mask_array,1)]
end

function awesome_print_debug(cost, mask_array, row_cover, column_cover, step, location)
    println("cost matrix = ")
    @show cost.m
    println("mask array = ")
    println(mask_array)
    println("row cover = ")
    println(row_cover)
    println("column cover = ")
    println(column_cover)
    println("augmenting path start = $(location.row) $(location.col)")
    readline(STDIN)
end

function step_one!(cost)
    #remove row minimum from cost matrix and find locations of all zeros
    cost.row_offsets = vec(minimum(cost.m,2))
    zero_locations = [find(j->iszero(cost,i,j), 1:size(cost,2)) for i=1:size(cost,1)]
end

function step_two!(cost, mask_array, row_cover, column_cover)
    #find a zero in cost matrix and star the zero if there is no other zero in row or column
    #(i.e. set mask_array to 1)
    for i=1:size(cost,1)
        for j=1:size(cost,2)
            if iszero(cost, i, j) && !row_cover[i] && !column_cover[j]
                mask_array[i,j] = StarMark
                row_cover[i] = true
                column_cover[j] = true
            end
        end
    end
    row_cover[:] = false
    column_cover[:] = false
end

function step_three!(mask_array, column_cover)
    #cover each column with a starred zero. If k (minimum(n,m)) columns are covered the algorithm is done
    n,m = size(mask_array)
    column_count = 0
    step = 3
    for i=1:n
        for j=1:m
            if mask_array[i,j] == StarMark
                column_cover[j] = true
            end
        end
    end
    column_count = sum(column_cover)
    if (column_count >= m) || (column_count >= n)
       step = 7
    else
       step = 4
    end
    return step
end

function step_four!(cost, mask_array, row_cover, column_cover, zero_locations)
    # Find indices of all zeros in the uncovered rows.  Doing this as a prework
    # step is faster than doing it in the loop below, probably because the
    # memory traversal order is more cache friendly.

    for row=1:size(cost,1)
        if row_cover[row]
            continue
        end
        # Find column of first zero in cost matrix which is uncovered
        zero_column = 0
        for j = zero_locations[row]
            if !column_cover[j]
                zero_column = j
            end
        end
        if zero_column == 0
            continue
        end

        mask_array[row, zero_column] = PrimeMark
        starcolumn = 0
        for j = 1:size(mask_array,2)
            @inbounds if mask_array[row,j] == StarMark
                starcolumn = j
                break
            end
        end
        if starcolumn == 0
            return 5, Location(row,zero_column)
        end

        row_cover[row] = true
        column_cover[starcolumn] = false
    end

    return 6, Location()
end

function step_five!(mask_array, row_cover, column_cover, path_start)
    #Find a series of alternating starred and primed zeros in alternating columns and rows
    #until it terminates on a primed zero with no star in it's column. Unstar all the starred
    #zeros, star the primed zeros, erase the primes and uncover all lines.
    @assert valid(path_start)

    done = false
    row = -1
    column = -1

    path = Array(Location,0)
    push!(path, path_start)

    while ~done
        row = find_star_in_column(mask_array, row, path[end].col)
        if row > -1
            push!(path, Location(row, path[end].col))
        else
            done = true
        end
        if ~done
            column = find_prime_in_row(mask_array, path[end].row, column)
            push!(path, Location(path[end].row, column))
        end
    end
    augment_path!(path, mask_array)
    row_cover[:] = false
    column_cover[:] = false
    erase_primes!(mask_array)
    return 3
end

function step_six!(cost,row_cover,column_cover, zero_locations)
    #Subtract the minimum value in the uncovered part of the cost matrix from all uncovered
    #columns and add it to the covered rows. Note that for efficiency this does not explicitly
    #occur, we just keep track of where new zeros appear and where zeros disappear and update
    #the row and column minima.
    min_value, min_locations = find_smallest_uncovered(cost,row_cover,column_cover)

    #min locations in the uncovered rows become the new zeros
    for i=1:length(min_locations)
        if !isempty(min_locations[i])
            for j=1:length(min_locations[i])
                push!(zero_locations[i],min_locations[i][j])
            end
        end
    end

    cost.row_offsets[row_cover] -= min_value
    cost.column_offsets[!column_cover] += min_value

    #need to deal with any zeros going away in covered columns and rows
    for i = 1:length(zero_locations)
        if !isempty(zero_locations[i]) && row_cover[i]
            for j=1:length(zero_locations[i])
                if column_cover[zero_locations[i][j]]
                    zero_location = zero_locations[i]
                    zero_locations[i] = zero_location[zero_location .!= zero_locations[i][j]]
                    break
                end
            end
        end
    end

    return 4
end

function find_star_in_column(mask_array, row, active_column)
    #find the row in which a star occurs for the active column
    row = -1
    for i = 1:size(mask_array,1)
        if mask_array[i,active_column] == StarMark
            row = i
        end
    end
    return row
end

function find_prime_in_row(mask_array, active_row, column)
    #find the column in which a prime occurs for the active row
    column = -1
    for j = 1:size(mask_array,2)
        if mask_array[active_row,j] == PrimeMark
            column = j
        end
    end
    return column
end

function augment_path!(path, mask_array)
    #go along the path swapping stars for zeros
    for p = 1:length(path)
        mask_array[path[p].row, path[p].col] = (mask_array[path[p].row, path[p].col] == StarMark) ? 0 : StarMark
    end
end

function erase_primes!(mask_array)
    #remove all of the primes from the mask_array
    for j=1:size(mask_array,2)
        for i=1:size(mask_array,1)
            if mask_array[i,j] == PrimeMark
                mask_array[i,j] = 0
            end
        end
    end
end

function find_smallest_uncovered(cost, row_cover, column_cover)
    #find the locations and value of the minimum of the cost matrix in the uncovered rows and columns
    min_value = typemax(eltype(cost))
    uncovered_row_inds = find(!row_cover);
    min_locations = [Int[] for i=1:size(cost,1)]
    for j = 1:size(cost,2)
        if !column_cover[j]
            for i = uncovered_row_inds
                @inbounds c = cost[i,j]
                if c < min_value
                    #have found a new minimum, so throw away all the previously discovered values and start again
                    min_value = c
                    for k=1:length(min_locations)
                        if !isempty(min_locations[k])
                            empty!(min_locations[k])
                        end
                    end
                    push!(min_locations[i],j)
                elseif  c == min_value
                    push!(min_locations[i],j)
                end
            end
        end
    end

    return min_value, min_locations
end

end # module


text:
# Quasi-Newton.

function optimize(f::Function, g::Function, H::Function, x0::Vector{T},
                  approxh::Bool; kmax::Int64 = 5000, δ::Float64 = 1e-6) where {T<:Real}
    k::Int64 = 0
    δ *= δ
    x = x0
    n::Int64 = length(x)
    dfx = ones(n)
    d2fx = eye(n, n)
    while norm(dfx) > δ && k < kmax
        prev = x
        dfx = g(x)
        x -= d2fx * dfx
        y = g(x) - dfx
        s = x - prev
        d2fx = H(d2fx, y, s)
        k += 1
    end
    return x, k
end


text:
#!/bin/bash
# -*- mode: julia -*-
#=
JULIA="${JULIA:-julia --color=yes --startup-file=no}"
export JULIA_PROJECT="$(dirname ${BASH_SOURCE[0]})"
exec ${JULIA} "${BASH_SOURCE[0]}" "$@"
=#

using Pkg
Pkg.instantiate()
Pkg.develop(PackageSpec(path=dirname(@__DIR__)))


text:
"""
    fitBaselli(S::Vector{<:Real}, I::Vector{<:Real}, ρ::Vector{<:Real}; ...)

Fit the Baselli model from measured or generated data.

# Arguments

* S: Systolic blood pressure
* I: RR Interval series
* ρ: Respiration

# Keyword Arguments

* na, nb, nc: order of the respective transfer functions, default to 10
* iterations: maximum number of iterations, defaults to 10
* δmin: minimum changes that defines convergence, defaults to 0.001
* verbose: get more information
* estimator: defaults to `\\` least squares

# Return

Returns a BaselliModel from which physiological indices can be inferrede by `postprocess()`

# Example
```julia
    model = getModel(1)
    result = predict(model, 100)
    modelEst = fitBaselli(result...)
```
"""
function fitBaselli(S::Vector{<:Real}, I::Vector{<:Real}, ρ::Vector{<:Real};
     na::Int = 10,
     nb::Int = 10,
     nd::Int = 10,
     iterations::Int = 10,
     δmin::Real = 0.001,
     verbose::Bool = false,
     estimator = \)

     # Check arguments
    length(S) == length(I) == length(ρ) || throw(ArgumentError("S, I, ρ must be of the same length"))
    na >= 1 || throw(ArgumentError("na ($na) must be positive"))
    nb >= 1 || throw(ArgumentError("nb ($nb) must be positive"))
    nd >= 1 || throw(ArgumentError("nd ($nd) must be positive"))
    iterations >= 1 || throw(ArgumentError("iterations($iterations) must be positive"))
    δmin >= 0 || throw(ArgumentError("δmin($δmin) must be positive"))

    # Prepare inputs
    N = length(S)
    T = 1

    us = [I ρ]
    ds = iddata(S, us, T)

    ut = [S ρ]
    dt = iddata(I, ut, T)

    # Estimate Transfer Functions using generalized least squares, either the maximum likelihood version or repeated arxar
    # ARXXAR
    Gs, DS, residual_s = arxar(ds, na, [1, nb], nd, δmin = δmin, iterations = iterations, verbose = verbose, estimator = estimator)

    # XXAR - na = 0
    Gt, DI, residual_t = arxar(dt, 0, [nb, nb], nd, inputdelay = [0, 1], δmin = δmin, iterations = iterations, verbose = verbose, estimator = estimator)

    # Extract model
    # G = G_SI / (1 - G_SS) with a priori knowledge gSI(1) = g(1) gives G_SI and G_SS
    G = Gs[1,1]
    gSI = vec(impulse(G,1)[1]')
    G_SI = impRes2tf(gSI)
    G_SS = ControlSystems.robust_minreal(-(G_SI / G) + 1)
    R1 = Gs[1,2]
    R_Sρ = ControlSystems.robust_minreal(R1 * (1 - G_SS))
    MS = ControlSystems.robust_minreal(DS * (1 - G_SS))

    G_IS = Gt[1,1]
    R_Iρ = Gt[1,2]
    MI = DI

    # Fit AR model representing respiration
    arCoeffs, errr = lpc(ρ, nd)
    Mρ = tf(1, ([1; arCoeffs]), 1)
    # TODO use ar(;stochastic = stochastic)

    # Check residuals for witheness(normal distribution) via 5% anderson darling test
    pS_anderson = pvalue(OneSampleADTest(residual_s, Normal(mean(residual_s), std(residual_s))))
    pI_anderson = pvalue(OneSampleADTest(residual_t, Normal(mean(residual_t), std(residual_t))))

    res = BaselliModel(G_SI = G_SI, G_IS = G_IS, G_SS = G_SS, MS = MS, MI = MI, Mρ = Mρ, R_Sρ = R_Sρ, R_Iρ = R_Iρ, wS = Normal(0, std(residual_s)), wI = Normal(0, std(residual_t)), wρ = Normal(0, sqrt(errr)), n = na, N = N, pS_anderson = pS_anderson, pI_anderson = pI_anderson, res_S = residual_s, res_I = residual_t, Gs = Gs, Gt = Gt, DS = DS,  DI = DI)
    return res
end

"""
    postprocess(M::BaselliModel, S::Vector{<:Real}, I::Vector{<:Real})

Extract physiologically relevant indices from the Baselli model by means of post processing as defindes in:

Baselli et al. (1994)  *Model for the assessment of heart period and arterial pressure variability interactions and of respiration influences*

# Indices

* α: Measure for baroreflex sensitivity
* gSI: Mechanical effects and properties of the myocardium, Starling's law
* gSS: Mechanical effects of blood pressure regulation (Windkessel)
* GssLF: Index of the contribution of AP regulation to the genesis of LF oscillations
* Θt/Θs: Respiratory effects
* PusLF/PutLF: Residual signals and detection of different L F sources

"""
function postprocess(M::BaselliModel, S::Vector{<:Real}, I::Vector{<:Real})
    # BRS
    # ramp = 1:140
    # rampRes = lsim(M.G_IS, ramp, ramp)[1][:]
    # α = theilSenEstimator(ramp, rampRes)[1]
    gIS = impulse(M.G_IS)[1]
    α = sum(gIS) # "BRS"

    # gSI stroke volume effect (+) or run off effect (-)
    gSI = vec(impulse(M.G_SI,1)[1]')[2]

    ## AP regulation
    # gSS represents atrial resistance and compliance (Windkessel)
    gSS = vec(impulse(M.G_SS,1)[1]')[2]
    # Gss(LF), resonance frequency of closed G_SS loop -> index of contribution of AP regulation LF oscilation
    Gss = 1 / (1 - M.G_SS)
    mag, phase, w = bode(Gss, range(0.0, stop = π, step = 0.0100))
    GssLF = maximum(mag)

    ## Respiratory effects - Share of HF power due to respiration / HF unwell defined for lof f_resp
    # Effect of resipration on I through different pathways
    # TF = minreal((1-M.G_SS)*M.R_Iρ + M.G_IS * M.R_Sρ) * Mρ / (1 - M.G_SS - M.G_SI * M.G_IS)
    # 1) R_Sρ is set to 0
    TF1 = ControlSystems.robust_minreal((1-M.G_SS)*M.R_Iρ + M.G_IS * 0) * M.Mρ / (1 - M.G_SS - M.G_SI * M.G_IS)
    # 2) R_Iρ is set to 0
    TF2 = ControlSystems.robust_minreal((1-M.G_SS)*0 + M.G_IS * M.R_Sρ) * M.Mρ / (1 - M.G_SS - M.G_SI * M.G_IS)
    wHF = range(0.15*2π, stop = 0.4*2π, step = 0.001)
    mag1, phase, w = bode(TF1, wHF)
    mag2, phase, w = bode(TF2, wHF)
    Pt1 =  mag1[:] .* (M.wρ.σ ^2)
    Pt2 =  mag2[:] .* (M.wρ.σ ^2)
    Θt = log10(sum(Pt1) / sum(Pt2))

    # Effect of respiration on S through different pathways
    # Tf = (M.G_SI * M.R_Iρ + M.R_Iρ) * Mρ / (1 - M.G_SS - M.G_SI * M.G_IS)
    TF1 = (M.G_SI * M.R_Iρ + 0) * M.Mρ / (1 - M.G_SS - M.G_SI * M.G_IS)
    TF2 = (M.G_SI * 0 + M.R_Sρ) * M.Mρ / (1 - M.G_SS - M.G_SI * M.G_IS)
    mag1, phase, w = bode(TF1, wHF)
    mag2, phase, w = bode(TF2, wHF)
    Pt1 =  mag1[:] .* (M.wρ.σ ^2)
    Pt2 =  mag2[:] .* (M.wρ.σ ^2)
    Θs = log10(sum(Pt1) / sum(Pt2))

    # LF sources in residuals
    LF = (0.04, 0.15)
    sLFpow = getSpectralComponentAR(M.MS, LF, σ = M.wS.σ)
    PusLF = sLFpow / var(S .- mean(S))
    tLFpow = getSpectralComponentAR(M.MS, LF, σ = M.wS.σ)
    PutLF = tLFpow / var(I .- mean(I))

    return (α = α, gSI = gSI, gSS = gSS, GssLF = GssLF,  Θt =  Θt, Θs = Θs, PusLF = PusLF, PutLF = PutLF)
end

# Simple Parameter of Goodness for fitted model (R² & R²adjusted)
function getGoodness(M::BaselliModel, S::Vector{<:Real}, I::Vector{<:Real})
    # R²
    vars = var(S .- mean(S))
    vart = var(I .- mean(I))
    ps = ((vars - (M.wS.σ ^2)) / vars)
    pt = ((vart - (M.wI.σ ^2)) / vart)
    # R²adjusted
    ps_adjusted = 1 - ((M.N - 1) / (M.N - M.n)) * (1 - ps)
    pt_adjusted = 1 - ((M.N - 1) / (M.N - M.n)) * (1 - pt)
    return (ps = ps, pt = pt, ps_adjusted = ps_adjusted, pt_adjusted = pt_adjusted)
end


text:
## Run operations

"""
    numruns(x)

Count the number of runs of repeated values present in a vector.
"""
function numruns(re::AbstractVector)
    n = 1
    current = re[1]
    @inbounds for i = 2:length(re)
        if re[i] != current
            n += 1
            current = re[i]
        end
    end
    n
end

"""
    ree(x)

Run End Encode a vector

Like RLE, but returns (runvalues,runends) rather than (runvalues,runlengths)
"""
function ree(x::AbstractVector)
    xlen = length(x)
    xlen < 2 && return ((x, [xlen]))
    nrun = numruns(x)
    runvalues = similar(x, nrun)
    runends = Vector{Int}(undef, nrun)
    run = 1
    current = x[1]
    @inbounds for i = 2:xlen
        if x[i] != current
            runvalues[run] = current
            runends[run] = i - 1
            current = x[i]
            run = run + 1
        end
    end
    runvalues[nrun] = current
    runends[nrun] = xlen
    runvalues, runends
end

"""
    numruns(runvalues, runends)

Given run values and run ends for a RLEVector, determine the number of runs that would
    be present if it were re-compressed. RLEVectors.jl does this operation after modifying
    an RLEVector, for example.
"""
function numruns(runvalues, runends)
    len = length(runends)
    length(runends) != len && throw(ArgumentError("runvalues and runends must be the same length."))
    len < 2 && return (len)
    n = 1
    current_val = runvalues[1]
    current_end = runends[1]
    @inbounds for i = 2:len
        rv = runvalues[i]
        re = runends[i]
        if rv != current_val && re != current_end
            n = n + 1
            if re < current_end
                throw(ArgumentError("The provided runends were not sorted, please use cumsum(runlengths) to get the right values."))
            end
            current_val = rv
            current_end = re
        end
    end
    return (n)
end

"""
    ree(runvalues, runends)
    ree!(runvalues, runends)
    ree!(x::RLEVector)

Tidy up an existing (mostly) Run End Encoded vector, dropping zero length runs and fixing
any adjacent identical values. `RLEVectors.jl` does this operation after modifying an
RLEVector, for example.
"""
function ree(runvalues, runends)
    ree!(copy(runvalues), copy(runends))
end

function ree!(runvalues, runends)
    n = length(runvalues)
    n == length(runends) || throw(DimensionMismatch(
        "runends and runvalues must be the same length (length(runends)=$(length(runends)), length(runvalues)=$n)."))
    left_i = 0
    if (n >= 1)
        current_val = runvalues[1]
        current_end = runends[1]
        if current_end != 0
            left_i = 1
        end
        @inbounds for right_i = 2:length(runvalues)
            rv = runvalues[right_i]
            re = runends[right_i]
            if re > current_end
                if !isequal(rv, current_val)
                    left_i = left_i + 1
                    current_val = runvalues[left_i] = rv
                end
            elseif re < current_end
                throw(ArgumentError("RLEVector run ends must be sorted"))
            end
            current_end = runends[left_i] = re
        end
    end
    if left_i < n
        resize!(runvalues, left_i)
        resize!(runends, left_i)
    end
    runvalues, runends
end

function ree(x)
    return (([x], [1]))
end

"""
    inverse_ree(rle)

Uncompress the runs and runends of an RLEVector.

## Examples
    collect(rle)
    inverse_ree( runvalues(rle), runends(rle) )
"""
function inverse_ree(runvalues, runends)
    len = length(runvalues)
    len != length(runends) &&
        throw(ArgumentError("runvalues and runends must be of the same length."))
    len == 0 && return (similar(runvalues, 0))
    n = runends[end]
    rval = similar(runvalues, n)
    j = 1
    @inbounds for i = 1:n
        rval[i] = runvalues[j]
        if runends[j] == i
            j = j + 1
        end
    end
    return (rval)
end


text:
# Our developed libraries
using CombinatorialSpaces
using CombinatorialSpaces.ExteriorCalculus

using Catlab
using Catlab.Present
using Catlab.Graphics
using Catlab.CategoricalAlgebra
using Catlab.Programs
using LinearAlgebra

# Julia community libraries
using MeshIO
using CairoMakie
using DifferentialEquations
using Logging: global_logger
using TerminalLoggers: TerminalLogger
using JSON
global_logger(TerminalLogger())

using Decapodes.Simulations
using Decapodes.Examples
using Decapodes.Diagrams
using Decapodes.Schedules
using Decapodes.Debug
using Decapodes.OpenDiagrams

using CombinatorialSpaces: volume
using SparseArrays

######################
# Functions for Curl #
######################
include("./special_ops.jl")
include("./deca_schema.jl")

draw_dwd(dwd) = begin
  new_res = deepcopy(dwd)
  new_res.diagram[:outer_out_port_type] .= nothing
  new_res.diagram[:outer_in_port_type] .= nothing
  new_res.diagram[:out_port_type] .= nothing
  new_res.diagram[:in_port_type] .= nothing
  to_graphviz(new_res)
end
draw_equations(eq) = begin
  to_graphviz(eq,
  node_labels=true, prog="neato",
  node_attrs=Dict(:shape=>"oval"),
  graph_attrs=Dict(:nodesep=>"4.0"))
end

##################
# Define Physics #
##################

@present Flow2DQuantities <: ExtendedFreeExtCalc2D begin
  X::Space

  k::Hom(Form1(X), Form1(X))    # diffusivity (usually scalar multiplication)
  kᵨ::Hom(Form1(X), Form1(X))   # REMOVE: Diffusion of density
  R₀::Hom(Form0(X), Form0(X))    # Ideal gas constant (usually scalar multiplication)
  kᵥ::Hom(Form1(X), Form1(X))    # viscosity (usually scalar multiplication)
  e2t::Hom(Form0(X), Form0(X))    # scale between internal energy and temperature (ideal gas)
  t2e::Hom(Form0(X), Form0(X))    # inverse of e2t (ideal gas)
  e2p::Hom(Form0(X), Form0(X))    # # scale between internal energy × density and pressure (ideal gas)
  L₀::Hom(Form1(X)⊗DualForm2(X), DualForm2(X))
  L₁::Hom(Form1(X)⊗DualForm1(X), DualForm1(X))
  L₁′::Hom(Form1(X)⊗Form1(X), Form1(X))
  ∂₀::Hom(Form0(X), Form0(X)) # all_wall boundary condition
  ∂₀₊::Hom(Form0(X), Form0(X)) # top/bottom wall boundary condition
  ∂₀ₛ::Hom(Form0(X), Form0(X)) # Sticking boundary condition
  ∂₁ₛ::Hom(Form1(X), Form1(X)) # Sticking boundary condition
  ∂₁ₗ₊::Hom(Form1(X), Form1(X)) # Sticking boundary condition
  ∂₁ₑ::Hom(Form1(X), Form1(X)) # In/Out edge flow boundary condition
  dneg₁::Hom(DualForm1(X), DualForm1(X)) # In/Out edge flow boundary condition
  neg₁::Hom(Form1(X), Form1(X)) # In/Out edge flow boundary condition
  dneg₀::Hom(DualForm0(X), DualForm0(X)) # In/Out edge flow boundary condition
  neg₀::Hom(Form0(X), Form0(X)) # In/Out edge flow boundary condition
  mask₁ₑ::Hom(Form1(X), Form1(X)) # In/Out edge flow boundary condition
  mask₀ₗ::Hom(Form0(X), Form0(X)) # In/Out edge flow boundary condition
  mask₀ₛ::Hom(Form0(X), Form0(X)) # In/Out edge flow boundary condition
  half::Hom(Form0(X), Form0(X)) # half
  R₀⁻¹::Hom(Form0(X), Form0(X)) # division by R₀
  third::Hom(Form1(X), Form1(X))    # multiplication by 1/3
  i₁′::Hom(Form1(X)⊗Form1(X), Form0(X))
  mult₀::Hom(Form0(X)⊗Form0(X), Form0(X)) # elementwise multiplication
  div₁::Hom(Form1(X)⊗Form1(X), Form1(X)) # elementwise division
  div₀::Hom(Form0(X)⊗Form0(X), Form0(X)) # elementwise division
  avg₀₁::Hom(Form0(X), Form1(X))
  sub₀::Hom(Form0(X)⊗Form0(X), Form0(X))
end

Diffusion = @decapode Flow2DQuantities begin
  (T, Ṫ)::Form0{X}
  ϕ::DualForm1{X}

  # Fick's first law
  ϕ ==  ⋆₁{X}(k(d₀{X}(T)))
  # Diffusion equation
  Ṫ == ⋆₀⁻¹{X}(dual_d₁{X}(ϕ))
end

Advection = @decapode Flow2DQuantities begin
  (T, Ṫ)::Form0{X}
  V::Form1{X}
  Ṫ == neg₀(⋆₀⁻¹{X}(L₀(V, ⋆₀{X}(T))))
end

Superposition = @decapode Flow2DQuantities begin
  (Ṫ₁, Ṫ₂, Ṫ, T)::Form0{X}
  Ṫ == Ṫ₁ + Ṫ₂
  ∂ₜ{Form0{X}}(T) == Ṫ
end

compose_diff_adv = @relation (T, Ṫ, V) begin
  diffusion(T, Ṫ₁)
  advection(T, Ṫ₂, V)
  superposition(Ṫ₁, Ṫ₂, Ṫ, T)
end
#=
DiffusionAdvection = oapply(compose_diff_adv,
                  [OpenDiagram(Diffusion, [:T, :Ṫ]),
                   OpenDiagram(Advection, [:T, :Ṫ, :V]),
                   OpenDiagram(Superposition, [:Ṫ₁, :Ṫ₂, :Ṫ, :T])])

dwd_res = diag2dwd(diff_adv.functor, in_vars = [:V, :T])
draw_equations(diff_adv.functor)
draw_dwd(dwd_res)
=#
NavierStokes = @decapode Flow2DQuantities begin
  (V, V̇, G, τ)::Form1{X}
  (T, ρ, ρ̇, p, δV)::Form0{X}
  δV == δ₁{X}(V)
  τ == div₁(kᵥ(Δ₁{X}(V) + third(d₀{X}(δV))), avg₀₁(ρ))
  V̇ == neg₁(L₁′(V, V)) +
        d₀{X}(half(i₁′(V, V))) +
        τ +
        neg₁(div₁(d₀{X}(p),avg₀₁(ρ))) +
        G
  ∂ₜ{Form1{X}}(V) == V̇
  ρ̇ == neg₀(⋆₀⁻¹{X}(L₀(V, ⋆₀{X}(ρ))))# + ⋆₀⁻¹{X}(dual_d₁{X}(⋆₁{X}(kᵨ(d₀{X}(ρ)))))
  ∂ₜ{Form0{X}}(ρ) == ρ̇
end

Energy = @decapode Flow2DQuantities begin
  (V, τ)::Form1{X}
  (e, ėₐ, ėₜ, ė, ρ, p, δV, T, Ṫ)::Form0{X}

  T == R₀⁻¹(div₀(p, ρ))
  ėₜ == Ṫ + i₁′(τ, V)

  p == e2p(mult₀((e + neg₀(half(i₁′(V, V)))) , ρ))

  ėₐ == neg₀(⋆₀⁻¹{X}(L₀(V, ⋆₀{X}(mult₀(e,ρ) + p))))
  ė == div₀(ėₜ + ėₐ, ρ)
  ∂ₜ{Form0{X}}(e) == ė
end

BoundaryConditions = @decapode Flow2DQuantities begin
  (V, V̇, bc₁)::Form1{X}
  (ė, bc₀)::Form0{X}
  # no-slip edges
  ∂₁ₛ(V) == bc₁
  ∂₁ₛ(V̇) == bc₁
  # No change on left/right boundaries
  ∂₀ₛ(ė) == bc₀
end

dwd_res = diag2dwd(NavierStokes)
draw_equations(NavierStokes)
draw_dwd(dwd_res)

compose_heat_xfer = @relation (V, ρ) begin
  flow(V, V̇, δV, T, ρ, p, τ)
  energy(ė, V, ρ, p, T, Ṫ, δV, τ)
  diffusion(T, Ṫ)
  bcs(ė, V, V̇)
end

HeatXfer = oapply(compose_heat_xfer,
                  [OpenDiagram(NavierStokes, [:V, :V̇, :δV, :T, :ρ, :p, :τ]),
                   OpenDiagram(Energy, [:ė, :V, :ρ, :p, :T, :Ṫ, :δV, :τ]),
                   OpenDiagram(Diffusion, [:T, :Ṫ]),
                   OpenDiagram(BoundaryConditions, [:ė, :V, :V̇])])

draw_equations(HeatXfer.functor)

new_dwd = diag2dwd(HeatXfer.functor, in_vars=[:e, :ρ, :V, :G])
draw_dwd(new_dwd)

################################
# Complex Operator Definitions #
################################
@present ExtendedOperators <: ExtendedFreeExtCalc2D begin
  X::Space
  neg::Hom(DualForm1(X), DualForm1(X)) # negative
  neg₁::Hom(Form1(X), Form1(X)) # negates 1-forms
  dneg₁::Hom(DualForm1(X), DualForm1(X)) # negates 1-forms
  L₀::Hom(Form1(X)⊗DualForm2(X), DualForm2(X))
  L₁::Hom(Form1(X)⊗DualForm1(X), DualForm1(X))
  i₀::Hom(Form1(X)⊗DualForm2(X), DualForm1(X))
  i₁::Hom(Form1(X)⊗DualForm1(X), DualForm0(X))

  L₁′::Hom(Form1(X)⊗Form1(X), Form1(X))
  i₀′::Hom(Form1(X)⊗Form2(X), Form1(X))
  i₁′::Hom(Form1(X)⊗Form1(X), Form0(X))
  ∧₁₁′::Hom(Form1(X)⊗DualForm1(X), DualForm2(X))
  ∧₁₀′::Hom(Form1(X)⊗DualForm0(X), Form1(X))
end

rules = gen_dec_rules()

# Lie derivative between two primal 1-forms
Lie1Imp′ = @decapode ExtendedOperators begin
  (F1, F1′, Ḟ1)::Form1{X}
  Ḟ1 == i₀′(F1, d₁{X}(F1′)) + d₀{X}(i₁′(F1, F1′))
end
lie1_imp′ = diag2dwd(Lie1Imp′, in_vars = [:F1, :F1′], out_vars = [:Ḟ1])
rules[:L₁′] = lie1_imp′

# Internal product between a primal 1-form and a primal 2-form
I0Imp′ = @decapode ExtendedOperators begin
  (F1, Ḟ1)::Form1{X}
  F2::Form2{X}
  Ḟ1 == neg₁(∧₁₀′(F1, ⋆₂{X}(F2)))
end
i0_imp′ = diag2dwd(I0Imp′, in_vars = [:F1, :F2], out_vars = [:Ḟ1])
rules[:i₀′] = i0_imp′

# Internal product between two primal 1-forms
I1Imp′ = @decapode ExtendedOperators begin
  (F1, F1′)::Form1{X}
  F0::Form0{X}
  F0 == ⋆₀⁻¹{X}(∧₁₁′(F1, ⋆₁{X}(F1′)))
end
i1_imp′ = diag2dwd(I1Imp′, in_vars = [:F1, :F1′], out_vars = [:F0])
rules[:i₁′] = i1_imp′

##########################
# Simulation Computation #
##########################

exp_dwd = Examples.expand_dwd(new_dwd, rules)
draw_dwd(exp_dwd)

s1 = EmbeddedDeltaSet2D("../cavity_unstructured_64.stl")
s = EmbeddedDeltaSet2D{Bool, Point{3, Float64}}()
copy_parts!(s, s1)

mod_p = findall(p -> (0.71 <= p[1] <= 0.72) && (0.99 <= p[2] <= 0.998), s[:point])
s[mod_p[1], :point] = s[mod_p[1], :point] .- (0.0,0.001,0.0)

sd = dual(s);
if ⋆(Val{1}, sd)[1,1] < 0.0
  orient_component!(s, 1, false)
  sd = dual(s);
end

cₚ = 1004.703 # Specific Heat at constant pressure
kₜ = 10 * 0.0246295028571 #Thermal conductivity
kᵦ = 1.38064852e-23 # Boltzmann constant (m² kg/(s² K))

density = 0.00597782417156
R = kᵦ * 6.0221409e23 # kg⋅m²/(s²*K*mol)
mol_mass = 28.96 # g/mol
μ = 28.96 / 6.0221409e23 # mean molecular mass (g)
R₀ = R / (mol_mass / 1000)

# Heat diffusion constant in fluid
k₁ = kₜ
k_col = fill(k₁, ne(s))
k = diagm(k_col)

kᵨ = 1e-3

ν = 1.716e-5 # 0.0005081150545826

kᵥ = ν# / density

γ = 1 + R₀ / (cₚ - R₀) # Mayer's formula + definition of adiabatic constant
println("γ: $γ")
# T = (γ - 1) * (mol_mass / 1000) / R₀ * e
#e2t = (μ/1000) / kᵦ
#e2t = (γ - 1) * mol_mass / R₀
e2t = (γ - 1) / R₀  #at 0 velocity
@show e2t
e2p = (γ - 1)
#e2p = (γ - 1)
t2e = 1/e2t

############################
# Define BCs and Operators #
############################

funcs = sym2func(sd)
∂₀ = Examples.boundary_inds(Val{0}, s)
∂₁ = Examples.boundary_inds(Val{1}, s)

∂ₗ₀ = ∂₀[findall(p-> p[1] <= 1e-4, s[∂₀, :point])]
∂ᵣ₀ = ∂₀[findall(p-> p[1] >= 1.0 - 1e-4, s[∂₀, :point])]
∂ₜ₀ = ∂₀[findall(p-> p[2] >= 1.0 - 1e-4, s[∂₀, :point])]
∂ᵦ₀ = ∂₀[findall(p-> p[2] <= 1e-4, s[∂₀, :point])]
∂ₑ₀ = vcat(∂ₗ₀, ∂ᵣ₀, ∂ₜ₀, ∂ᵦ₀)

∂ₗ₁ = Examples.bound_edges(s, ∂ₗ₀)
∂ᵣ₁ = Examples.bound_edges(s, ∂ᵣ₀)
∂ₑ₁ = Examples.bound_edges(s, ∂ₑ₀)

# Just induced edges
∂₁₊ = Examples.bound_edges(s, ∂₀)

# Full adjacent boundary
#∂₁₊ = Examples.adj_edges(s, ∂₀)
# Just edges and corners
#=
∂_corners = Examples.adj_edges(s,
                vcat(intersect(∂ₗ₀, ∂ₜ₀), intersect(∂ₗ₀, ∂ᵦ₀),
                     intersect(∂ᵣ₀, ∂ₜ₀), intersect(∂ᵣ₀, ∂ᵦ₀)))
@show ∂_corners
∂₁₊ = collect(unique(vcat(∂ₑ₁, ∂_corners)))
=#
# Full triangle on boundary

#∂₁₊ = Examples.adj_edges(s, ∂₀)
#∂_points = unique(vcat(s[∂₁₊, :∂v0], s[∂₁₊, :∂v1]))
#∂₁₊ = Examples.bound_edges(s, ∂_points)

∂ₗ₁₊ = Examples.adj_edges(s, ∂ₗ₀)
∂ᵣ₁₊ = Examples.adj_edges(s, ∂ᵣ₀)
∂ₑ₁₊ = Examples.adj_edges(s, ∂ₑ₀)

∂ₗ₀₊ = unique(vcat(s[∂ₗ₁₊, :∂v1], s[∂ₗ₁₊, :∂v0]))
∂ᵣ₀₊ = unique(vcat(s[∂ᵣ₁₊, :∂v1], s[∂ᵣ₁₊, :∂v0]))
∂ₑ₀₊ = unique(vcat(s[∂ₑ₁₊, :∂v1], s[∂ₑ₁₊, :∂v0]))

funcs[:mcopy] = Dict(:operator => I(ne(sd)), :type => MatrixFunc())
funcs[:k] = Dict(:operator => k, :type => MatrixFunc())
funcs[:e2p] = Dict(:operator => e2p * I(nv(sd)), :type => MatrixFunc())
funcs[:e2t] = Dict(:operator => e2t * I(nv(sd)), :type => MatrixFunc())
funcs[:t2e] = Dict(:operator => (1/e2t) * I(nv(sd)), :type => MatrixFunc())
funcs[:half] = Dict(:operator => 0.5 * I(nv(sd)), :type => MatrixFunc())
funcs[:R₀⁻¹] = Dict(:operator => 1/R₀ * I(nv(sd)), :type => MatrixFunc())
funcs[:third] = Dict(:operator => (1/3) * I(ne(sd)), :type => MatrixFunc())
funcs[:R₀] = Dict(:operator => R₀ * I(nv(sd)), :type => MatrixFunc())
funcs[:kᵥ] = Dict(:operator => kᵥ * I(ne(sd)), :type => MatrixFunc())
funcs[:kᵨ] = Dict(:operator => kᵨ * I(ne(sd)), :type => MatrixFunc())
funcs[:dneg₁] = Dict(:operator => -1 * I(ne(sd)), :type => MatrixFunc())
funcs[:mult₀] = Dict(:operator => (x′,x,y) -> (x′ .= x .* y), :type => InPlaceFunc())
funcs[:plus] = Dict(:operator => (x′,x,y) -> (x′ .= x .+ y), :type => InPlaceFunc())
funcs[:sum₀] = Dict(:operator => (x′,x,y) -> (x′ .= x .+ y), :type => InPlaceFunc())
funcs[:sum₁] = Dict(:operator => (x′,x,y) -> (x′ .= x .+ y), :type => InPlaceFunc())
funcs[:sum₂] = Dict(:operator => (x′,x,y) -> (x′ .= x .+ y), :type => InPlaceFunc())
funcs[:sum₀̃] = Dict(:operator => (x′,x,y) -> (x′ .= x .+ y), :type => InPlaceFunc())
funcs[:sum₁̃] = Dict(:operator => (x′,x,y) -> (x′ .= x .+ y), :type => InPlaceFunc())
funcs[:sum₂̃] = Dict(:operator => (x′,x,y) -> (x′ .= x .+ y), :type => InPlaceFunc())
funcs[:sub₀] = Dict(:operator => (x′,x,y) -> (x′ .= x .- y), :type => InPlaceFunc())
funcs[:div₁] = Dict(:operator => (x′,x,y) -> (x′ .= x ./ y), :type => InPlaceFunc())
funcs[:div₀] = Dict(:operator => (x′,x,y) -> (x′ .= x ./ y), :type => InPlaceFunc())
funcs[:avg₀₁] = Dict(:operator => avg_mat(Val{(0,1)}, s), :type => MatrixFunc())
funcs[:neg₁] = Dict(:operator => -1 * I(ne(sd)), :type => MatrixFunc())
funcs[:neg₀] = Dict(:operator => -1 * I(nv(sd)), :type => MatrixFunc())
funcs[:avg] = Dict(:operator => -1 * I(nv(sd)), :type => MatrixFunc())

funcs[:∂ₗ₀₊] = Dict(:operator => (x′,x) -> (x′ .= x; x′[∂ₗ₀₊] .= 0), :type => InPlaceFunc())

const ramp_up = 0.1 # seconds
# Final temperatures should be
temperature = 0.0 #t2e * 172.89

# Pre-allocate boundary condition values
const ∂ₗ₀_z = zeros(length(∂ₗ₀))
const ∂ᵣ₀_z = zeros(length(∂ᵣ₀))
const ∂ₗ₀_t = fill(temperature / ramp_up, length(∂ₗ₀))
const ∂ᵣ₀_t = fill(-(temperature / ramp_up), length(∂ᵣ₀))
ramp_up_temp(x′, x, t) = begin
  x′ .= x
  if t < ramp_up
    setindex!(x′, ∂ₗ₀_t, ∂ₗ₀)
    setindex!(x′, ∂ᵣ₀_t, ∂ᵣ₀)
  else
    setindex!(x′, ∂ₗ₀_z, ∂ₗ₀)
    setindex!(x′, ∂ᵣ₀_z, ∂ᵣ₀)
  end
end

const ∂₁₊_z = zeros(length(∂₁₊))
const ∂ₗ₁₊_z = zeros(length(∂ₗ₁₊))
const ∂ᵣ₁₊_z = zeros(length(∂ᵣ₁₊))
const ∂ₑ₀₊_z = zeros(length(∂ₑ₀₊))
funcs[:∂₀ₛ] = Dict(:operator => ramp_up_temp, :type => TDInPlaceFunc())
funcs[:∂₁ₛ] = Dict(:operator => (x′,x) -> (x′ .= x; setindex!(x′, ∂₁₊_z ,  ∂₁₊)), :type => InPlaceFunc())
funcs[:∂₁ₑ] = Dict(:operator => (x′,x) -> (x′ .= x; setindex!(x′, ∂ₗ₁₊_z, ∂ₗ₁₊); setindex!(x′, ∂ᵣ₁₊_z, ∂ᵣ₁₊)), :type => InPlaceFunc())
funcs[:mask₀ₗ] = Dict(:operator => (x′,x) -> (x′ .= x; setindex!(x′, ∂ₑ₀₊_z, ∂ₑ₀₊)), :type => InPlaceFunc())

wedge_cache = init_wedge_ops(sd)
v2comp = comp_support(sd);
cache_mat = Dict(:t2c => tri2comp(s, v2comp), :e2c => edge2comp(s, v2comp), :cross => changes(sd, v2comp),
                 :α_cache => zeros(ntriangles(sd)*3), :β_cache => zeros(ntriangles(sd)*3))

funcs[:∧₁₀′] = Dict(:operator => (x′, α, β) -> (cp_2_1!(x′, β, α, cache_mat)), :type => InPlaceFunc())
funcs[:∧₁₁′] = Dict(:operator => (x′, α, β) -> (pd_wedge!(x′, Val{(1,1)},s, α, β; wedge_cache...)), :type => InPlaceFunc())

const hodge_LU = lu(⋆(Val{1}, sd))
funcs[:⋆₁⁻¹] = Dict(:operator => (x′,x) -> (x′ .= -1 * (hodge_LU \ x)), :type => InPlaceFunc())

t2e = 1/e2t
c_objs = fill(t2e * 288.15, nv(s))
#c_objs[∂ₗ₀] .= t2e * 461.04
#c_objs[∂ᵣ₀] .= t2e * 115.26
velocity(p) = [0.0, 0.0, 0.0]
gravity(p) = [0.0,-9.81,0.0]
v = ♭(sd, DualVectorField(velocity.(sd[triangle_center(sd),:dual_point]))).data;
g = ♭(sd, DualVectorField(gravity.(sd[triangle_center(sd),:dual_point]))).data;
ρ = [density for p in s[:point]]

new_dwd = deepcopy(exp_dwd)
to_graphviz(new_dwd, orientation=LeftToRight)
Examples.zip_dwd!(new_dwd)
to_graphviz(new_dwd, orientation=LeftToRight)
Examples.contract_matrices!(new_dwd, funcs)
to_graphviz(new_dwd, orientation=LeftToRight)

func, code = gen_sim(new_dwd, funcs, sd; autodiff=false, params=[:G]);
dt = 0.01
tseg = 1.0
for i in 1:10
  prev_loc = "res_$(i-1)"
  loc = "res_$i"
  mkpath(loc)
  u0 = vcat(c_objs, ρ, v)
  if i > 1
    pre_cond = open("$(prev_loc)/sim_res.json", "r") do f
      JSON.parse(f)
    end
    u0 .= pre_cond["$tseg"]
  end
  GC.gc()

  ##############
  # Simulation #
  ##############
  prob = ODEProblem(func, u0, (0.0, tseg))
  sol1 = solve(prob, Tsit5(), progress=true, progress_steps=1, dtmax=1.0e-5, saveat=dt, p=g)

  #################
  # Write Results #
  #################
  t_range = 1:nv(s)
  ρ_range = (1:nv(s)) .+ nv(s)
  v_range = (1:ne(s)) .+ (2 * nv(s))

  using JSON
  res = Dict{Float64, Vector}()
  for i in 0.0:dt:sol1.t[end]
      res[i] = sol1(i)
  end
  open("$(loc)/sim_res.json", "w") do f
      JSON.print(f, res)
  end

  using WriteVTK
  for i in 1:length(dt:dt:sol1.t[end])
    tris = triangle_vertices(s)
    cells = [MeshCell(VTKCellTypes.VTK_TRIANGLE, [tris[1][t],tris[2][t],tris[3][t] ]) for t in 1:ntriangles(s)];
    x =  [p[1] for p in s[:point]]
    y =  [p[2] for p in s[:point]]
    z =  [p[3] for p in s[:point]]
    vtkfile = vtk_grid("$(loc)/su2_coupled_v1_$(lpad(i, 4, "0"))", x, y, z, cells)
    vel_form = ♯(sd, sol1(dt*i)[v_range], CombinatorialSpaces.DiscreteExteriorCalculus.PPSharp())

    vtkfile["temperature", VTKPointData()] = e2t * sol1(dt*i)[t_range]
    vtkfile["density", VTKPointData()] = sol1(dt*i)[ρ_range]
    vtkfile["vel", VTKPointData()] = vel_form
    vtkfile["v_mag", VTKPointData()] = sqrt.(abs.(1e-7 .+ inv_hodge_star(Val{0}, sd)*pd_wedge(Val{(1,1)}, sd, sol1(dt*i)[v_range], ⋆(Val{1}, sd) * sol1(dt*i)[v_range])))
    vtk_save(vtkfile)
  end

  fig = Figure()
  ax, ob = mesh(fig[1,1], s, color = e2t * sol1[end][t_range], colormap=:seismic)#, colorrange=(115 + 100,462 - 100))

  xlims!(ax, (0,1))
  ylims!(ax, (0,1))
  ax.aspect = AxisAspect(1)
  Colorbar(fig[1,2], ob)
  fig
  save("$(loc)/heat_prof.png", fig)

  times = range(0.01, sol1.t[end]/2, length=150)
  colors = [e2t * sol1(t)[t_range] for t in times]
  fig = Figure()
  axis, ob = mesh(fig[1,1], s, color=colors[1], colormap=:seismic)
  axis.aspect = AxisAspect(1.0)
  Colorbar(fig[1,2], ob)
  framerate = 30
  fig

  record(fig, "$(loc)/conc_flow.gif", collect(1:length(collect(times))); framerate = framerate) do i
    ob.color = colors[i]
  end

  times = range(0.01, sol1.t[end], length=150)
  colors = [sol1(t)[ρ_range] for t in times]

  fig = Figure()
  axis, ob = mesh(fig[1,1], s, color=colors[1],
                                     colorrange=(minimum(vcat(colors...)),maximum(vcat(colors...))))

  axis.aspect = AxisAspect(1.0)
  Colorbar(fig[1,2], ob)
  framerate = 30
  fig

  record(fig, "$(loc)/density_flow.gif", collect(1:length(collect(times))); framerate = framerate) do i
    ob.color = colors[i]
  end

  times = range(0, sol1.t[end], length=150)
  colors = [sol1(t)[t_range] .* sol1(t)[ρ_range] for t in times]

  figure, axis, ob = mesh(s, color=colors[1],
                                     colorrange=(minimum(vcat(colors...)),maximum(vcat(colors...))))

  axis.aspect = AxisAspect(1.0)
  framerate = 30

  record(figure, "$(loc)/press_flow.gif", collect(1:length(collect(times))); framerate = framerate) do i
    ob.color = colors[i]
  end

  magnitudes(t) = sqrt.(abs.(1e-4 .+ inv_hodge_star(Val{0}, sd)*pd_wedge(Val{(1,1)}, sd, sol1(t)[v_range], ⋆(Val{1}, sd) * sol1(t)[v_range])))
  times = range(0.01, sol1.t[end], length=150)
  colors = [magnitudes(t) for t in times]
  fig = Figure()
  axis, ob = mesh(fig[1,1], s, color=colors[3], colorrange=(minimum(vcat(colors...)),maximum(vcat(colors...))))

  axis.aspect = AxisAspect(1.0)
  Colorbar(fig[1,2], ob)
  framerate = 30
  fig
  record(fig, "$(loc)/vel_flow.gif", collect(1:length(collect(times))); framerate = framerate) do i
    ob.color = colors[i]
  end
end

max_stars_repo_name: AlgebraicJulia/DECAPODES-Benchmarks

max_issues_repo_name: AlgebraicJulia/DECAPODES-Benchmarks

max_forks_repo_name: AlgebraicJulia/DECAPODES-Benchmarks


text:
module Econometrics
using Statistics, LinearAlgebra, Printf, Test
using Plots, GR, StatsPlots
using Optim, NLopt
using Distributions, Distances, StatsBase, StatsFuns, StatsModels
using SpecialFunctions
using DataFrames, KernelDensity, Calculus, DataFramesMeta, CSV
# Utilities
#include("Utilities/montecarlo.jl")
include("Utilities/stnorm.jl")
include("Utilities/ma.jl")
include("Utilities/trim.jl")
include("Utilities/prettyprint.jl")
include("Utilities/eye.jl")
include("Utilities/sortbyc.jl")
include("Utilities/dstats.jl")
include("Utilities/lag.jl")
include("Utilities/lags.jl")
include("Utilities/vech.jl")
include("Utilities/PrintDivider.jl")
include("Utilities/PrintEstimationResults.jl")
include("Utilities/clc.jl")
include("Utilities/BoundByQuantiles.jl")
# linear regression
include("LinearRegression/lsfit.jl")
include("LinearRegression/ols.jl")
include("LinearRegression/ridge.jl")
include("LinearRegression/rsq.jl")
include("LinearRegression/tsls.jl")
include("LinearRegression/TestStatistics.jl")
# Bayesian
include("Bayesian/mcmc.jl")
# nonparametrics
include("NP/npreg.jl")
include("NP/InteriorPoint.jl")
include("NP/npdens.jl")
include("NP/bootstrap.jl")
include("NP/npdensity.jl")
include("NP/kernelweights.jl")
include("NP/NeweyWest.jl")
#=
#include("NN/TrainNet.jl")
#include("NN/AnalyzeNet.jl")
=#
# optimization
include("Optimization/samin.jl")
include("Optimization/fminunc.jl")
include("Optimization/fmincon.jl")

# MLE
include("ML/mle.jl")
include("ML/mleresults.jl")
include("ML/Likelihoods/logit.jl")
include("ML/Likelihoods/poisson.jl")
include("ML/Likelihoods/normal.jl")

# GMM
include("GMM/gmm.jl")
include("GMM/gmmresults.jl")

export stnorm, ma, trim, prettyprint, eye, sortbyc, dstats, lag, lags, vech
export clc, PrintDivider, PrintEstimationResults, BoundByQuantiles!
export lsfit, ols, ridge, rsq, tsls, TestStatistics, NeweyWest
export mcmc
export npreg, npdens, kernelweights, bootstrap
export npdensity, InteriorPoint
#export TrainNet, AnalyzeNet
export samin, fminunc, fmincon
export mle, mleresults
export logit, poisson, normal
export gmm, gmmresults
# export montecarlo
end

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
const INITIAL_SIZE = 10
const SHRINK_FACTOR = 0.1
const EXPAND_FACTOR = 1.0

struct HashMap{T,R}
  storage::Vector{Vector{Pair{T,R}}}
end

function HashMap{T,R}(storage_size::Int = INITIAL_SIZE) where {T,R}
  [Pair{T,R}[] for _ in 1:storage_size] |> HashMap{T,R}
end

function HashMap(storage_size::Int = INITIAL_SIZE)
  HashMap{Any,Any}(storage_size)
end

function Base.size(hm::HashMap) :: Tuple{Int}
  size(hm.storage)
end

function Base.length(hm::HashMap) :: Int
  length(hm.storage)
end

function index(hm::HashMap, key) :: Int
  hash(key) % (length(hm)+2)
end

function Base.keys(hm::HashMap)
  (p[1] for b in hm.storage for p in b)
end

function Base.values(hm::HashMap)
  (p[2] for b in hm.storage for p in b)
end

function Base.pairs(hm::HashMap)
  (p for b in hm.storage for p in b)
end

function loadfactor(hm::HashMap) :: Float64
  length(collect(keys(hm))) / length(hm)
end

function Base.haskey(hm::HashMap{T,R}, key::T)::Bool where {T,R}
  key in keys(hm)
end

function rehash!(hm::HashMap{T,R})::HashMap where {T,R}
  (length(hm) <= INITIAL_SIZE || (rand(1:10) != 10)) && return hm

  lf = loadfactor(hm)
  _hm = if lf > EXPAND_FACTOR
    HashMap{T,R}(ceil(Int, length(hm) * 2))
  elseif lf < SHRINK_FACTOR
    HashMap{T,R}(ceil(Int, length(hm) / 2))
  else
    return hm
  end

  for k in collect(keys(hm))
    _hm[k] = hm[k]
  end

  _hm
end

function Base.setindex!(hm::HashMap{T,R}, val::R, key::T)::R where {T,R}
  bucket = hm.storage[index(hm, key)]

  for i in 1:length(bucket)
    if bucket[i][1] == key
      bucket[i] = Pair{T,R}(key, val)
      return val
    end
  end

  push!(bucket, Pair{T,R}(key, val))
  rehash!(hm)

  val
end

function Base.getindex(hm::HashMap{T,R}, key::T)::R where {T,R}
  bucket = hm.storage[index(hm, key)]

  for i in 1:length(bucket)
    if bucket[i][1] == key
      return bucket[i][2]
    end

    throw(KeyError(key))
  end
end

function Base.get(hm::HashMap{T,R}, key::T, default::R)::R where {T,R}
  haskey(hm, key) ? hm[key] : default
end

function Base.get!(hm::HashMap{T,R}, key::T, default::R)::R where {T,R}
  haskey(hm, key) && return hm[key]

  hm[key] = default
end

function Base.getkey(hm::HashMap{T,R}, key::T, default::T)::T where {T,R}
  haskey(hm, key) ? key : default
end

function Base.delete!(hm::HashMap{T,R}, key::T)::HashMap{T,R} where {T,R}
  bucket = hm.storage[index(hm, key)]

  for i in 1:length(bucket)
    if bucket[i][1] == key
      deleteat!(bucket, i)
    end
  end

  return hm
end

function Base.pop!(hm::HashMap{T,R}, key::T, default::Union{R,Nothing} = nothing)::R where {T,R}
  if ! haskey(hm, key)
    default !== nothing && return default
    throw(KeyError(key))
  end

  val = hm[key]
  delete!(hm, key)

  val
end

function Base.iterate(hm::HashMap{T,R}, state::Union{T,Nothing} = first(keys(hm))) where {T,R}
  ks = keys(hm)

  idx = findfirst(x -> x == state, ks)
  idx === nothing && return nothing

  return length(ks) >= idx+1 ? (hm, ks[idx+1]) : nothing
end

function Base.empty!(hm::HashMap{T,R})::HashMap{T,R} where {T,R}
  for (k,_) in hm
    delete!(hm, k)
  end

  hm
end

function Base.isempty(hm::HashMap):: Bool
  length(hm.storage) == 0
end

hm = HashMap{Symbol,String}()
hm[:a] = "a"
hm[:b] = "b"

dump(hm)
@show hm[:b]
@show collect(keys(hm))
@show size(hm)
@show length(hm)

for c in 'A':'Z'
  hm[Symbol(c)] = string(c)
end

dump(hm)
@show collect(keys(hm))
@show size(hm)
@show length(hm)

max_stars_repo_path: Chapter 5/HashMaps.jl

max_issues_repo_path: Chapter 5/HashMaps.jl

max_forks_repo_path: Chapter 5/HashMaps.jl


text:
"""
	euc(u::CartesianIndex{3}, v::CartesianIndex{3})

Compute the euclidean distance given two 3 dimensional CartesianIndex arrays
"""
function euc(u::CartesianIndex{3}, v::CartesianIndex{3})
    return √((u[1] - v[1])^2 + (u[2] - v[2])^2 + (u[3] - v[3])^2)
end

"""
	find_edges(A)

Find the edges of an array by computing the distance to the nearest background pixel `0` and keeping
only the values that have a distance of 1. This corresponds to an array indice being adjacent to the
background
"""
function find_edges(A)
    # Set every pixel that is currently background (0) equal to 1, as these are the pixels we are interested in
    A = A .== 0

    # Find the indeces of these background pixels (now set equal to 1)
    A = ImageMorphology.feature_transform(A)

    # Find the distance away from the nearest background pixel (now set to 1)
    A = ImageMorphology.distance_transform(A)

    # Return the indeces where the distance is equal to 1 as this means this pixel is adjacent to a background pixel
    return (findall(x -> x == 1, A))
end

"""
    detect_edges_3D(img, f)

Modifies `detect_edges` to work with 3D images.
"""
function detect_edges_3D(img, f)
    container = Array{Int64}(undef, size(img))
    for k in 1:(size(img)[3])
        container[:, :, k] = detect_edges(img[:, :, k], f)
    end
    return container
end


text:
function AdRIL(;π,
                S,
                ΔN=50,
                solver=SAC,
                𝒟_demo,
                normalize_demo::Bool=true,
                expert_frac=0.5,
                buffer_size = 1000,
                buffer_init=0,
                log::NamedTuple=(;),
                buffer::ExperienceBuffer = ExperienceBuffer(S, action_space(π), buffer_size, [:i]),
                kwargs...)

    !haskey(𝒟_demo, :r) && error("AdRIL requires a reward value for the demonstrations")
    normalize_demo && (𝒟_demo = normalize!(deepcopy(𝒟_demo), S, action_space(π)))
    𝒟_demo = 𝒟_demo |> device(π)

    function AdRIL_callback(𝒟)
        max_i = maximum(𝒟[:i])
        k = Int((max_i - buffer_init) / ΔN) - 1
        old_data = 𝒟[:i][:] .<= max_i - ΔN
        new_data = 𝒟[:i][:] .> max_i - ΔN
        𝒟[:r][1, old_data] .= -1/k
        𝒟[:r][1, new_data] .= 0
    end

    solver(;π=π,
            S=S,
            ΔN=ΔN,
            post_experience_callback=AdRIL_callback,
            extra_buffers=[𝒟_demo],
            buffer_fractions=[1-expert_frac, expert_frac],
            buffer_size=buffer_size,
            buffer_init=buffer_init,
            buffer=buffer,
            log=(dir="log/adril", period=500, log...),
            kwargs...)
end


text:
function StepNormalizer(arg0::jdouble, arg1::ODEFixedStepHandler)
    return StepNormalizer((jdouble, ODEFixedStepHandler), arg0, arg1)
end

function StepNormalizer(arg0::jdouble, arg1::ODEFixedStepHandler, arg2::StepNormalizerBounds)
    return StepNormalizer((jdouble, ODEFixedStepHandler, StepNormalizerBounds), arg0, arg1, arg2)
end

function StepNormalizer(arg0::jdouble, arg1::ODEFixedStepHandler, arg2::StepNormalizerMode)
    return StepNormalizer((jdouble, ODEFixedStepHandler, StepNormalizerMode), arg0, arg1, arg2)
end

function StepNormalizer(arg0::jdouble, arg1::ODEFixedStepHandler, arg2::StepNormalizerMode, arg3::StepNormalizerBounds)
    return StepNormalizer((jdouble, ODEFixedStepHandler, StepNormalizerMode, StepNormalizerBounds), arg0, arg1, arg2, arg3)
end

function handle_step(obj::StepNormalizer, arg0::ODEStateInterpolator, arg1::jboolean)
    return jcall(obj, "handleStep", void, (ODEStateInterpolator, jboolean), arg0, arg1)
end

function init(obj::StepNormalizer, arg0::ODEStateAndDerivative, arg1::jdouble)
    return jcall(obj, "init", void, (ODEStateAndDerivative, jdouble), arg0, arg1)
end


text:
# IupGetParam: Example
#  Shows a dialog with all the possible fields.

module getparam_

using Printf
using IUP

export
    getparam

#------------------------example html/examples/C/getparam.c ----------------------------------
function param_action end
param_action_c = @cfunction(param_action, Cint, (Ptr{Ihandle}, Int, Ptr{Cvoid}))

function getparam()

    pboolean::Cint = 1
    pinteger::Cint = 3456
    preal::Float32 = 3.543
    pinteger2::Cint = 192
    preal2::Float32 = 0.5
    pangle::Float32 = 90
    pstring = "string text"
    pcolor = "255 0 128"
    plist::Cint = 2
    pstring2 = "second text\nsecond line"
    pfile_name = "test.jpg"

    IupOpen()                       #Initializes IUP
    IupControlsOpen()               # so we can show the IupDial, but can be tested without it
    IupSetLanguage("ENGLISH")

    if (IupGetParam("Title", param_action_c, C_NULL,
                   "Boolean: %b[No,Yes]{Boolean Tip}\n
                   Integer: %i{Integer Tip}\n
                   Real 1: %r{Real Tip}\n
                   Sep1 %t\n
                   Integer: %i[0,255]{Integer Tip 2}\n
                   Real 2: %r[-1.5,1.5]{Real Tip 2}\n
                   Sep2 %t\n
                   Angle: %a[0,360]{Angle Tip}\n
                   String: %s{String Tip}\n
                   List: %l|item1|item2|item3|{List Tip}\n
                   File: %f[OPEN|*.bmp;*.jpg|CURRENT|NO|NO]{File Tip}\n
                   Color: %c{Color Tip}\n
                   Sep3 %t\n
                   Multiline: %m{Multiline Tip}\n",
                   Ref(pboolean), Ref(pinteger), Ref(preal), Ref(pinteger2), Ref(preal2),
                   Ref(pangle), pstring, Ref(plist), pfile_name, pcolor, pstring2, C_NULL) != 0)
        return IUP_DEFAULT
    end

    # the float values (Real 1, Real 2 and Angle ) are not displayed correctly.
    IupMessagef("IupGetParam",
              "Boolean Value: %d\n
               Integer: %d\n
               Real 1: %g\n
               Integer: %d\n
               Real 2: %g\n
               Angle: %g\n
               String: %s\n
               List Index: %d\n
               FileName: %s\n
               Color: %s\n
               Multiline: %s",
               pboolean, pinteger, preal, pinteger2, preal2, pangle,
               pstring, plist, pfile_name, pcolor, pstring2)

    IupClose()                                  # And close it when ready
end

function param_action(dialog::Ptr{Ihandle}, param_index::Cint, user_data::Ptr{Cvoid})
    if (param_index == IUP_GETPARAM_OK)
        @printf("Button1 (OK)\n");
    elseif (param_index == IUP_GETPARAM_MAP)
        @printf("Map\n")
    elseif (param_index == IUP_GETPARAM_INIT)
        @printf("Init\n")
    elseif (param_index == IUP_GETPARAM_CANCEL)
        @printf("Button2 (Cancel)\n")
    elseif (param_index == IUP_GETPARAM_HELP)
        @printf("Button3 (Help)\n")
    elseif (param_index == IUP_GETPARAM_CLOSE)
        @printf("Close\n")
    else
        param_str = sprintf("PARAM%d", param_index)
        param = convert(Ptr{Ihandle}, IupGetAttribute(dialog, param_str))
        @printf("%s = %s\n", param_str, IupGetAttribute(param, "VALUE"))
    end

    ret::Cint = 1
    return ret
end

end     # module


text: x1 = 2

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using .Base:
	@propagate_inbounds, tail
import .Base:
	size, length, iterate, IteratorSize

struct Enamerate
    na::NamedArray
end

#TODO: docs
"""
    enamerate(iter)

An iterator that yields `(n, x)` where `n` is a tuple corresponding to
the names in the corresponding dimension of value x. Similar to enumerate
but with names.
"""
enamerate(n::NamedArray) = Enamerate(n)

length(e::Enamerate) = length(e.na)
size(e::Enamerate) = size(e.na)
@propagate_inbounds function iterate(e::Enamerate, state=(1,))
    i, rest = state[1], tail(state)
    n = iterate(e.na, rest...)
    n === nothing && return n
    #= (i, n[1]), (i+1, n[2]) =#
		return (flattenednames(e.na)[i], n[1]), (i+1, n[2])
end

# Possibly TODO:
# adapt this from enumerate
# eltype(::Type{Enumerate{I}}) where {I} = Tuple{Int, eltype(I)}
# IteratorSize(::Type{Enumerate{I}}) where {I} = IteratorSize(I)
#= IteratorEltype(::Type{Enumerate{I}}) where {I} = IteratorEltype(I) =#

#= @inline function iterate(r::Reverse{<:Enumerate}) =#
#=     ri = reverse(r.itr.itr) =#
#=     iterate(r, (length(ri), ri)) =#
#= end =#
#= @inline function iterate(r::Reverse{<:Enumerate}, state) =#
#=     i, ri, rest = state[1], state[2], tail(tail(state)) =#
#=     n = iterate(ri, rest...) =#
#=     n === nothing && return n =#
#=     (i, n[1]), (i-1, ri, n[2]) =#
#= end =#
#=  =#


text:
module Spring

    using Plots
    using LaTeXStrings #for latex string in plot labels

    export SpringParam,InitParam,AnimParam
    export animate_spring,animate_spring2D
    gr()

    include("params.jl")
    include("animate1D.jl")
    include("animate2D.jl")

end # module


text:
using AccurateArithmetic
using ForwardDiff
using LinearAlgebra

#this allows to bail out on other numeric types, but
#keep using the accurate arithmetic package when available

function accurate_dot(x::AbstractArray{T1},y::AbstractArray{T2}) where {T1<:Union{Float32,Float64},T2<:Union{Float32,Float64}}
    return dot_oro(x,y)
end

function accurate_dot(x,y)
    return dot(x,y)
end

_isbad(x) = isnan(x) || isinf(x)
_mod1(x, m) = (x - 1) % m + 1

Base.@kwdef mutable struct State{T<:Real}
    check_obj    :: Bool
    autodiff     :: Bool
    Lp           :: T
    maps_limit   :: T
    time_limit   :: T
    buffer       :: T
    t0           :: T
    tol          :: T

    go_on        :: Bool = true
    converged    :: Bool = false
    maps         :: Int64 = 0 # Or the # of g! calls if optim == true.
    k            :: Int64 = 0
    f_calls      :: Int64 = 0
    σ            :: T = zero(T)
    α            :: T = one(T)
    obj_now      :: T = T(Inf)
    norm_∇       :: T = T(Inf)

    ix₀          :: Int64 = 1
    ix           :: Int64 = 1
    ix_best      :: Int64 = 1
    ix_new       :: Int64 = 1
    ord_best     :: Int64 = 0
    i_ord        :: Int64 = 0
    α_best       :: T = one(T)
    σ_mult_fail  :: T = one(T)
    σ_mult_loop  :: T = one(T)
    α_mult       :: T = one(T)
    norm_best    :: T = T(Inf)
    obj_best     :: T = T(Inf)
    α_boost      :: T = one(T)

    σs           :: Vector{T} = zeros(T,10)
    norm_∇s      :: Vector{T} = zeros(T,10)
    σs_i         :: Int64 = 1
end

function Base.eltype(::State{F}) where F
    return F
end
#####
##### Initialization funtions
#####

function check_arguments(
    f, g!, m!, x_in :: T, lower :: Union{T, Nothing}, upper :: Union{T, Nothing},
    s :: State
) where {T<:AbstractArray}

    if m! ≠ nothing && g! ≠ nothing
        throw(ArgumentError("must not provide both m! and g!"))
    end
    if s.check_obj && f === nothing
        throw(ArgumentError("if check_obj == true, f must be provided."))
    end
    if (upper ≠ nothing && maximum(x_in .- upper) > 0) ||
       (lower ≠ nothing && maximum(lower .- x_in) > 0)
        throw(DomainError(x_in, "infeasible starting point"))
    end

    if !((eltype(x_in) <: Real) | (eltype(x_in) <: Integer))
        throw(ArgumentError("starting point must be of type Float"))
    end
    if s.autodiff
        #@info "minimizing f using gradient descent acceleration and ForwardDiff"
    end
    return nothing
end

function α_too_large!(
    f, g_auto, g!, s :: State{F}, ∇temp :: T, ∇ :: T, x_in :: T, ∇∇ :: F
) :: Bool where {F,T<:AbstractArray}

    if f ≠ nothing
        obj_new = f(x_in - s.α * ∇)
        s.f_calls += 1
        return _isbad(obj_new) || (obj_new > s.obj_now - F(0.25)*s.α * ∇∇)
    else
        s.autodiff ? ∇temp .= g_auto(x_in .- s.α * ∇) : g!(∇temp, x_in .- s.α * ∇)
        s.maps += 1
        ∇∇temp = ∇temp ⋅ ∇temp
        return _isbad(∇∇temp) || ∇∇temp / ∇∇ > 2
    end
end

function initialize_α!(
    f, g_auto, g!, s :: State{F}, ∇ :: T, x_in :: T, temp :: T
) where {F,T<:AbstractArray}
    ∇∇ = F(∇ ⋅ ∇)
    if iszero(∇∇)
        throw(DomainError(x_in, "∇f(x_in) = 0 (extremum or saddle point)"))
    end

    is_too_large = α_too_large!(f, g_auto, g!, s, ∇, temp, x_in, ∇∇)
    min_mult = 2.0
    mult = min_mult * 64^2
    for i ∈ 0:30
        s.α *= mult^(1 - 2is_too_large)
        was_too_large = is_too_large
        is_too_large = α_too_large!(f, g_auto, g!, s, temp, ∇, x_in, ∇∇)
        if mult == min_mult && (is_too_large + was_too_large == 1) || s.α > 1e10
            s.α /= min_mult^(is_too_large && !was_too_large)
            break
        end
        mult = max(mult / 64^(is_too_large + was_too_large == 1), min_mult)
    end
    return nothing
end

#####
##### Core functions
#####

# Makes sure x_try[i] stays within boundaries with a gap buffer * (bound[i] - x_old[i]).
function bound!(
    extr, x_try :: T, x_old :: T, bound :: T, buffer
) where {T<:AbstractArray}

    for i ∈ eachindex(x_try)
        x_try[i] = extr(x_try[i], (1 - buffer) * bound[i] + buffer * x_old[i])
    end
    return nothing
end

function update_x!(
    g_auto, g!, m!, x_out :: T, ∇ :: T, s :: State,
    x_in :: T, lower :: Union{T, Nothing}, upper :: Union{T, Nothing}
) where {T<:AbstractArray}

    if m! === nothing
        if s.maps > 0 # If maps == 0, ∇ was already computed but maps must still be updated.
            s.autodiff ? ∇ .= g_auto(x_in) : g!(∇, x_in)
        end
        x_out .= x_in .- s.α * ∇

        if lower ≠ nothing; bound!(max, x_out, x_in, lower, s.buffer) end
        if upper ≠ nothing; bound!(min, x_out, x_in, upper, s.buffer) end
        if lower ≠ nothing || upper ≠ nothing; ∇ .= (x_in .- x_out) ./ s.α end
    else
        ∇ .= x_in
        m!(x_out, x_in)
        ∇ .-= x_out

        # If the algorithm has reached an infeasible value, it helps stability
        # if we also slow down the change of x from mappings.
        if s.α < 1; x_out .+= (1 - s.α) * ∇ end
    end
    s.norm_∇ = norm(∇, s.Lp)
    return nothing
end

function check_convergence!(f, s :: State, x :: AbstractArray)
    if s.norm_∇ <= s.tol # only based on norm_∇ (for now)
        s.go_on = false

        if s.check_obj
            s.obj_now = f(x)
            s.f_calls += 1
        end
        s.converged = !s.check_obj || !_isbad(s.obj_now)
    end
    return nothing
end

function mapping!(
    f, g_auto, g!, m!, x_out :: T, ∇ :: T, s :: State, info,
    x_in :: T, lower :: Union{T, Nothing}, upper :: Union{T, Nothing}
) where {T<:AbstractArray}

    if s.go_on
        try
            update_x!(g_auto, g!, m!, x_out, ∇, s, x_in, lower, upper)
        catch
            if s.k == 1 # To initially catch any type or function errors.
                update_x!(g_auto, g!, m!, x_out, ∇, s, x_in, lower, upper)
            end
            s.go_on = false
        end
        s.ix += 1
        s.maps += 1
        !_isbad(s.norm_∇) ? check_convergence!(f, s, x_out) : s.go_on = false
        store_info!(info, copy(x_out), s, false)
    end
    return nothing
end

# prodΔ computes the dot products Δᵃ ⋅ Δᵇ and Δᵇ ⋅ Δᵇ needed to compute σ
# without creating the Δs themselves and minimizes cache misses.

function prodΔ(∇1 :: T, ∇2 :: T, temp :: T) where {T<:AbstractArray} # if p == 2
    temp .= ∇2 .- ∇1
    return accurate_dot(temp, ∇1), accurate_dot(temp, temp)
end

function prodΔ(∇1 :: T, ∇2 :: T, ∇3 :: T, temp :: T) where {T<:AbstractArray} # if p == 3
    temp .= ∇3 .- 2∇2 .+ ∇1
    return accurate_dot(temp, (∇2 .- ∇1)), accurate_dot(temp, temp)
end

function extrapolate!( # Not a real extrapolation, just a stabilizing step
    x_out :: T, x_in :: T, ∇ :: T, s :: State
) where {T<:AbstractArray}

    @. x_out = x_in - s.α * ∇
    return nothing
end

function extrapolate!(
    x_out :: T, x_in :: T, ∇1 :: T, ∇2 :: T, s :: State
) where {T<:AbstractArray}

    @. x_out = x_in - (s.α * (2s.σ - s.σ^2)) * ∇1 - (s.α * s.σ^2) * ∇2
    return nothing
end

function extrapolate!(
    x_out :: T, x_in :: T, ∇1 :: T, ∇2 :: T, ∇3 :: T, s :: State
) where {T<:AbstractArray}

    @. x_out = x_in - (s.α * (3s.σ - 3s.σ^2 + s.σ^3)) * ∇1 -
        (s.α * (3s.σ^2 - 2s.σ^3)) * ∇2 - (s.α * s.σ^3) * ∇3
    return nothing
end

function update_α!(s :: State{F}, σ_new :: F, ΔᵇΔᵇ :: F) where F

    # Increasing / decreasing α to maintain σ ∈ [1, 2] as much as possible.
    s.α *= F(1.5)^((σ_new > 2) - (σ_new < 1))

    # Boosting α if ΔᵇΔᵇ gets really small
    if ΔᵇΔᵇ < 1e-100
        s.α_boost *= 4 # Increasing more aggressively each time
        s.α = min(s.α * s.α_boost, one(F))
    end
    return nothing
end

#####
##### Functions to handle failures
#####

# update_progress! is used to know which iterate to fall back on when
# backtracking and, optionally, which x was the smallest minimizer in case the
# problem has multiple minima.
function update_progress!(f, s :: State{F}, x :: AbstractArray) where F

    if s.go_on && s.check_obj
        s.obj_now = f(x)
        s.f_calls += 1
        s.go_on = !_isbad(s.obj_now)
    end

    if s.go_on && (( s.check_obj && s.obj_now < s.obj_best)
                || (!s.check_obj && s.norm_∇ < s.norm_best))
        s.ix_best = s.ix₀
        s.ix_new = _mod1(s.ix₀ + 1, 2)
        s.ord_best, s.α_best = (s.i_ord, s.α)
        s.σ_mult_fail = s.α_mult = one(F)
        s.check_obj ? s.obj_best = s.obj_now : s.norm_best = s.norm_∇
    end
    return nothing
end

function backtrack!(s :: State)

    s.ix₀ = s.ix_best
    s.i_ord = s.ord_best
    s.σ_mult_fail /= 2
    s.α_mult /= 2
    s.α = s.α_best * s.α_mult
    if s.check_obj; s.obj_now = s.obj_best end
    return nothing
end

# The future of the algorithm is fully determined by x₀, α and io. Since it is
# non-monotonic, the same (x₀, α, io) could appear twice with non-zero
# probability, creating an infinite loop. Since the possibility is very small,
# not much resource should be devoted to it. We verify this by comparing
# values of σ and norm(∇) over time. If (σ, norm(∇)) appear twice within a
# certain number of iterations, we divide α and σ by 2.
function check_∞_loop!(s :: State, last_order :: Bool)

    if s.σ > s.σs[s.σs_i]
        s.σs[s.σs_i] = s.σ
        s.norm_∇s[s.σs_i] = s.norm_∇
    end

    l = length(s.σs)
    if last_order && s.σs[s.σs_i] > 0.001 && s.σ_mult_loop > 0.01 && s.σs_i == l
        loop = false
        for i ∈ 1:l - 1, j ∈ i+1:l
            if s.σs[i] > 1e-10 && abs(s.σs[i] - 1) > 1e-10 &&
                    abs(1 - s.σs[j]/s.σs[i]) < 1e-10 &&
                    abs(1 - s.norm_∇s[j]/s.norm_∇s[i]) < 1e-10
                loop = true
                s.σ_mult_loop /= 2
                s.α /= 2
                s.σs .= 0
                s.norm_∇s .= 0
                break
            end
        end
        if !loop; s.σ_mult_loop = 1.0 end
    end
    s.σs_i = _mod1(s.σs_i + last_order, l)
    return nothing
end

#####
##### Miscellaneous
#####

function store_info!(info, x, s :: State, extrapolating)

    if info ≠ nothing
        push!(info.x, x)
        push!(info.σ, s.σ)
        push!(info.α, s.α)
        push!(info.extrapolating, extrapolating)
    end
    return nothing
end

#####
##### Main funtion
#####

"""
    SpeedMapping
`speedmapping(x₀; m!, kwargs...)` accelerates the convergence of a mapping
`m!(x_out, x_in)` to a fixed point of `m!` by the Alternating cyclic
extrapolation algorithm (ACX). Since gradient descent is an example
of such mapping, `speedmapping(x0; g!, kwargs...)` can also perform multivariate
optimization based on the gradient function `g!(∇, x)`.

Reference:
N. Lepage-Saucier, _Alternating cyclic extrapolation methods for optimization
algorithms_, arXiv:2104.04974 (2021). https://arxiv.org/abs/2104.04974

### Arguments

*   `x₀ :: AbstractArray`: The starting point; must be of eltype `Float`.

Main keyword arguments:
*   `m!(x_out, x_in)`: A map for which a fixed-point must be found.
*   `g!(∇, x)`: The gradient of a function to be minimized.
*   `f(x)`: The objective of the function to be minimized. It is useful to *i*)
    compute a good initial α (learning  rate) for the gradient descent, *ii*)
    optimize using autodiff or *iii*) track the progress of the algorithm.
*   `lower, upper :: Union{AbstractArray, Nothing} = nothing`: Box constraints
    for the optimization. NOTE: When appropriate, it should also be used with
    `m!`. Even if `m!` always keeps `x_out` within bounds, an extrapolation
    step could throw `x` out of bounds.
*   `tol :: Float64 = 1e-10, Lp :: Real = 2` When using `m!`, the algorithm
    stops when `norm(F(xₖ) - xₖ, Lp) ≤ tol`. When using `g!`, the algorithm
    stops when `norm(∇f(xₖ), Lp) ≤ tol`.
*   `maps_limit :: Real = 1e6`: Maximum number of `m!` calls or `g!` calls.
*   `time_limit :: Real = 1000`: Maximum time in seconds.

Minor keyword arguments:
*   `orders :: Array{Int64} = [3, 3, 2]` determines `ACX`'s alternating order.
    Must be between 1 and 3 (where 1 means no extrapolation). The two recommended
    orders are [3, 2] and [3, 3, 2], the latter being *potentially* better for
    highly non-linear applications (see paper).
*   `check_obj :: Bool = false`: In case of `NaN` or `Inf` values, the algorithm
    restarts at the best past iterate. If `check_obj = true`, progress is
    monitored with the value of the objective (requires `f`).
    Otherwise, it is monitored with `norm(F(xₖ) - xₖ, Lp)`. Advantages of
    `check_obj = true`: more precise and if the algorithm convergences on a bad
    local minimum, the algorithm instead returns the best past iterate.
    Advantages of `check_obj = false`: for well-behaved applications, it avoids
    the effort and time of providing `f` and calling it at every
    iteration.
*   `store_info :: Bool = false`: Stores `xₖ`, `σₖ` and `αₖ` (see paper).
*   `buffer :: Float64 = 0.01` If `xₖ` goes out of bounds, it is brought back in
    with a buffer. Ex. `xₖ = buffer * xₖ₋₁ + (1 - buffer) * upper`. Setting
    `buffer = 0.001` may speed-up box-constrained optimization.

Keyword arguments to fine-tune fixed-point mapping acceleration (using `m!`):
*   `σ_min :: Real = 0.0`: Setting to `1` may avoid stalling
    (see paper).
*   `stabilize :: Bool = false`: performs a stabilization mapping before
    extrapolating. Setting to `true` may improve the performance for
    applications like the EM or MM algorithm (see paper).

# Example: Finding a dominant eigenvalue
```jldoctest
julia> using LinearAlgebra

julia> using SpeedMapping

julia> A = ones(10) * ones(10)' + Diagonal(1:10);

julia> function power_iteration!(x_out, x_in, A)
           mul!(x_out, A, x_in)
           x_out ./= maximum(abs.(x_out))
       end;

julia> res = speedmapping(ones(10); m! = (x_out, x_in) -> power_iteration!(x_out, x_in, A))
(minimizer = [0.41214914122181, 0.44095060739689546, 0.4740798646565511, 0.5125916147320679, 0.5579135738427364, 0.6120273727167589, 0.677766040697063, 0.7593262786058278, 0.8632012019116191, 1.0], maps = 16, f_calls = 0, converged = true, norm_∇ = 2.8042637093700587e-9)

julia> V = res.minimizer;

julia> dominant_eigenvalue = V'A * V / V'V
16.310005690792195

```
# Example: Minimizing a multidimensional Rosenbrock
```
julia> f(x) = sum(100 * (x[i,1]^2 - x[i,2])^2 + (x[i,1] - 1)^2 for i ∈ 1:size(x,1));

julia> function g!(∇, x)
           ∇[:,1] .=  400(x[:,1].^2 .- x[:,2]) .* x[:,1] .+ 2(x[:,1] .- 1)
           ∇[:,2] .= -200(x[:,1].^2 .- x[:,2])
           return nothing
       end;

julia> x₀ = 1.0 * [-4 -3 -2 -1; 0 1 2 3]';
```

Optimizing, providing f and g!
```
julia> speedmapping(x₀; f, g!)
(minimizer = [0.9999999999962509 0.9999999999924867; 0.9999999999941136 0.9999999999882022; 0.9999999999752885 0.9999999999504781; 0.9999999989360835 0.9999999978679094], maps = 156, f_calls = 11, converged = true, norm_∇ = 9.524035730285234e-10)
```

Optimizing without objective
```
julia> speedmapping(x₀; g!)
(minimizer = [1.0000000000000542 1.0000000000001086; 0.9999999999999626 0.999999999999925; 0.9999999999997893 0.9999999999995778; 0.9999999999998465 0.9999999999996921], maps = 142, f_calls = 0, converged = true, norm_∇ = 3.524659991326121e-13)
```

Optimizing without gradient
```
julia> speedmapping(x₀; f)
(minimizer = [0.9999999999973896 0.9999999999947689; 0.9999999999958996 0.9999999999917817; 0.9999999999829567 0.9999999999658452; 0.9999999992743016 0.9999999985456991], maps = 156, f_calls = 11, converged = true, norm_∇ = 6.496402053553401e-10)
```

Optimizing with a box constraint
```
julia> upper = [0.5ones(4) Inf * ones(4)];

julia> speedmapping(x₀; f, g!, upper)
(minimizer = [0.5 0.25; 0.5 0.24999999999999864; 0.49999999999999795 0.2499999999975679; 0.4999999999999807 0.2499999999773891], maps = 138, f_calls = 11, converged = true, norm_∇ = 7.196533863189794e-9)
```
"""
function speedmapping(
    x_in :: AbstractArray{F}; f = nothing, g! = nothing, m! = nothing,
    orders :: Array{Int64} = [3,3,2], σ_min :: Real = zero(F), stabilize :: Bool = false,
    check_obj :: Bool = false, tol = F(1e-8), Lp :: Real = 2,
    maps_limit :: Real = 1e6, time_limit :: Real = 1000,
    lower :: Union{AbstractArray, Nothing} = nothing,
    upper :: Union{AbstractArray, Nothing} = nothing, buffer  = F(0.01),
    store_info :: Bool = false
) where F <: Real

    s = State{F}(; autodiff = f ≠ nothing && m! === nothing && g! === nothing, tol,
        buffer, Lp, maps_limit, time_limit, t0 = time(), check_obj)
    g_auto = s.autodiff ? x -> ForwardDiff.gradient(f, x) : nothing

    type_x = F
    if lower !== nothing && eltype(lower) ≠ type_x; lower = type_x.(lower) end
    if upper !== nothing && eltype(upper) ≠ type_x; lower = type_x.(upper) end

    check_arguments(f, g!, m!, x_in, lower, upper, s)

    # Inserting stabilization steps
    if stabilize; orders = Int.(vec(hcat(ones(length(orders)),orders)')) end

    # Two x₀s to avoid copying at each improvement (maybe this is excessive optimization?)
    #using copy instead of similar because BigFloats have problems
    x₀ = [copy(x_in), similar(x_in)]
    xs = [copy(x₀[1]) for i ∈ 1:maximum(orders)]
    ∇s = [copy(x₀[1]) for i ∈ 1:maximum(orders)] # Storing ∇s is equivalent to Δs
    temp = copy(x₀[1]) # temp storage

    if f !== nothing && (m! === nothing || check_obj)
        s.obj_now = s.obj_best = f(x_in) # Useful for initialize_α and tracking progress
        s.f_calls += 1
    end

    if m! === nothing
        s.autodiff ? ∇s[1] .= g_auto(x_in) : g!(∇s[1], x_in) # Will update s.maps later
        initialize_α!(f, g_auto, g!, s, ∇s[1], x_in, temp)
    end

    info = store_info ? (x = [x_in], σ = [s.σ], α = [s.α], extrapolating = [false]) : nothing

    while !s.converged && s.maps ≤ maps_limit && time() - s.t0 ≤ time_limit
        s.k += 1
        s.i_ord += 1
        s.ix = 0 # Which x is currently beeing updated

        # Avoiding a stabilization step if the last extrapolation was close to 1
        s.i_ord += stabilize && s.i_ord % 2 == 1 && abs(s.σ - 1) < 0.01 && s.k > 1

        io = _mod1(s.i_ord, length(orders))
        p = orders[io]
        s.go_on = true

        mapping!(f, g_auto, g!, m!, xs[1], ∇s[1], s, info, x₀[s.ix₀], lower, upper)
        update_progress!(f, s, x₀[s.ix₀])
        for i ∈ 2:p
            mapping!(f, g_auto, g!, m!, xs[i], ∇s[i], s, info, xs[i - 1], lower, upper)
        end

        if !s.converged && s.go_on
            if p > 1
                ΔᵃΔᵇ, ΔᵇΔᵇ = prodΔ(∇s[1:p]..., temp)
                σ_new = F(abs(ΔᵃΔᵇ) > 1e-100 && ΔᵇΔᵇ > 1e-100 ? abs(ΔᵃΔᵇ) / ΔᵇΔᵇ : 1.0)
                s.σ = max(σ_min, σ_new) * s.σ_mult_fail * s.σ_mult_loop
            end

            if lower === nothing && upper === nothing
                extrapolate!(x₀[s.ix_new], x₀[s.ix₀], ∇s[1:p]..., s)
            else
                extrapolate!(temp, x₀[s.ix₀], ∇s[1:p]..., s)
                if lower ≠ nothing; bound!(max, temp, x₀[s.ix₀], lower, s.buffer) end
                if upper ≠ nothing; bound!(min, temp, x₀[s.ix₀], upper, s.buffer) end
                x₀[s.ix_new] .= temp
            end
            store_info!(info, copy(x₀[s.ix_new]), s, true)
            s.ix₀ = s.ix_new

            if p > 1
                if m! === nothing; update_α!(s, σ_new, F(ΔᵇΔᵇ)) end
                check_∞_loop!(s, io == length(orders))
            end
        elseif !s.converged && !s.go_on
            backtrack!(s)
        end
    end

    if s.maps > maps_limit
        @warn "Maximum mappings exceeded."
    elseif time() - s.t0 > time_limit
        @warn "Exceeded time limit of $time_limit seconds."
    end

    minimizer = s.check_obj && s.obj_best < s.obj_now ? x₀[s.ix_best] : xs[s.ix]

    return (minimizer = minimizer, maps = s.maps, f_calls = s.f_calls,
        converged = s.converged, norm_∇ = s.norm_∇)
end


text:
@testset "Polygons" begin

@testset "construction" begin
    points = Point2f0[
        (0,6),
        (0,0),
        (3,0),
        (4,1),
        (6,1),
        (8,0),
        (12,0),
        (13,2),
        (8,2),
        (8,4),
        (11,4),
        (11,6),
        (6,6),
        (4,3),
        (2,6)
    ]
    mesh = GLNormalMesh(points)
    faces = polygon2faces(points, Triangle{Int})
    @test eltype(faces) == Triangle{Int}
    @test facetype(mesh) == GLTriangle

end
end


text:
const _ComponentLabelReferences = Set{Tuple{UUIDs.UUID, String}}

struct _TimeSeriesRecord
    component_labels::_ComponentLabelReferences
    ts::TimeSeriesData
end

function _TimeSeriesRecord(component_uuid, label, ts)
    record = _TimeSeriesRecord(_ComponentLabelReferences(), ts)
    push!(record.component_labels, (component_uuid, label))
    return record
end

"""
Stores all time series data in memory.
"""
struct InMemoryTimeSeriesStorage <: TimeSeriesStorage
    data::Dict{UUIDs.UUID, _TimeSeriesRecord}
end

function InMemoryTimeSeriesStorage()
    storage = InMemoryTimeSeriesStorage(Dict{UUIDs.UUID, _TimeSeriesRecord}())
    @info "Created InMemoryTimeSeriesStorage"
    return storage
end

"""
Constructs InMemoryTimeSeriesStorage from an instance of Hdf5TimeSeriesStorage.
"""
function InMemoryTimeSeriesStorage(hdf5_storage::Hdf5TimeSeriesStorage)
    storage = InMemoryTimeSeriesStorage()
    for (component, label, time_series) in iterate_time_series(hdf5_storage)
        add_time_series!(storage, component, label, time_series)
    end

    return storage
end

check_read_only(storage::InMemoryTimeSeriesStorage) = nothing

function add_time_series!(
    storage::InMemoryTimeSeriesStorage,
    component_uuid::UUIDs.UUID,
    label::AbstractString,
    ts::TimeSeriesData,
    unused = nothing,
)
    uuid = get_uuid(ts)
    if !haskey(storage.data, uuid)
        @debug "Create new time series entry." uuid component_uuid label
        storage.data[uuid] = _TimeSeriesRecord(component_uuid, label, ts)
    else
        add_time_series_reference!(storage, component_uuid, label, uuid)
    end
end

function add_time_series_reference!(
    storage::InMemoryTimeSeriesStorage,
    component_uuid::UUIDs.UUID,
    label::AbstractString,
    ts_uuid::UUIDs.UUID,
)
    @debug "Add reference to existing time series entry." ts_uuid component_uuid label
    record = storage.data[ts_uuid]
    push!(record.component_labels, (component_uuid, label))
end

function remove_time_series!(
    storage::InMemoryTimeSeriesStorage,
    uuid::UUIDs.UUID,
    component_uuid::UUIDs.UUID,
    label::AbstractString,
)
    if !haskey(storage.data, uuid)
        throw(ArgumentError("$uuid is not stored"))
    end

    record = storage.data[uuid]
    component_label = (component_uuid, label)
    if !(component_label in record.component_labels)
        throw(ArgumentError("$component_label wasn't stored for $uuid"))
    end

    pop!(record.component_labels, component_label)
    @debug "Removed $component_label from $uuid."

    if isempty(record.component_labels)
        @debug "$uuid has no more references; delete it."
        pop!(storage.data, uuid)
    end
end

function get_time_series(
    storage::InMemoryTimeSeriesStorage,
    uuid::UUIDs.UUID;
    index = 0,
    len = 0,
)::TimeSeries.TimeArray
    if !haskey(storage.data, uuid)
        throw(ArgumentError("$uuid is not stored"))
    end

    if index != 0
        @assert len != 0
        end_index = index + len - 1
        return storage.data[uuid].ts.data[index:end_index]
    end

    return storage.data[uuid].ts.data
end

function clear_time_series!(storage::InMemoryTimeSeriesStorage)
    empty!(storage.data)
    @info "Cleared all time series."
end

function get_num_time_series(storage::InMemoryTimeSeriesStorage)
    return length(storage.data)
end

function convert_to_hdf5(storage::InMemoryTimeSeriesStorage, filename::AbstractString)
    create_file = true
    hdf5_storage = Hdf5TimeSeriesStorage(create_file; filename = filename)
    for record in values(storage.data)
        for pair in record.component_labels
            columns = TimeSeries.colnames(record.ts.data)
            add_time_series!(hdf5_storage, pair[1], pair[2], record.ts, columns)
        end
    end
end

function compare_values(x::InMemoryTimeSeriesStorage, y::InMemoryTimeSeriesStorage)::Bool
    keys_x = sort!(collect(keys(x.data)))
    keys_y = sort!(collect(keys(y.data)))
    if keys_x != keys_y
        @error "keys don't match" keys_x keys_y
        return false
    end

    for key in keys_x
        record_x = x.data[key]
        record_y = y.data[key]
        if record_x.component_labels != record_y.component_labels
            @error "component_labels don't match" record_x record_y
            return false
        end
        if TimeSeries.timestamp(record_x.ts.data) != TimeSeries.timestamp(record_y.ts.data)
            @error "timestamps don't match" record_x record_y
            return false
        end
        if TimeSeries.values(record_x.ts.data) != TimeSeries.values(record_y.ts.data)
            @error "values don't match" record_x record_y
            return false
        end
    end

    return true
end


text:
export getDivergenceMatrix

function getDivergenceMatrix(mesh::TetraMesh)

if isempty(mesh.Div)

  nf = mesh.nf
	nc = mesh.nc
	T  = mesh.Tetras

  ii = zeros(Int64,4*nc)
  jj = zeros(Int64,4*nc)
  kk = zeros(4*nc)

	ii = [1:nc
			  1:nc
			  1:nc
			  1:nc]

	jj = [findRow(mesh.faces,[T[:,1] T[:,2] T[:,3]])
	      findRow(mesh.faces,[T[:,1] T[:,2] T[:,4]])
			  findRow(mesh.faces,[T[:,1] T[:,3] T[:,4]])
			  findRow(mesh.faces,[T[:,2] T[:,3] T[:,4]])]
	kk = [ones(nc)
			 -ones(nc)
			  ones(nc)
			 -ones(nc)]

	#    cnt = 1
	#    for i=1:nc
	#    ii[cnt] = nc
	# 	jj[cnt] = findRow(mesh.faces,[T[i,1] T[i,2] T[i,3]])[1]
	# 	kk[cnt] = 1
	# 	cnt += 1
	#
	#    ii[cnt] = nc
	# 	jj[cnt] = findRow(mesh.faces,[T[i,1] T[i,2] T[i,4]])[1]
	# 	kk[cnt] = -1
	# 	cnt += 1
	#
	#    ii[cnt] = nc
	# 	jj[cnt] = findRow(mesh.faces,[T[i,1] T[i,3] T[i,4]])[1]
	# 	kk[cnt] = 1
	# 	cnt += 1
	#
	#    ii[cnt] = nc
	# 	jj[cnt] = findRow(mesh.faces,[T[i,2] T[i,3] T[i,4]])[1]
	# 	kk[cnt] = -1
	# 	cnt += 1
	# end

  F = getFaceArea(mesh)

  for k = 1:4*nc
    kk[k] *= F[jj[k]]      #divergence scaled using face sizes
  end

	mesh.Div = sparse(ii, jj, kk, nc, nf)

end

return mesh.Div

end

max_stars_repo_name: JuliaInv/JTetra

max_issues_repo_name: JuliaInv/JTetra

max_forks_repo_name: JuliaInv/JTetra


text:
import SeisIO
import SeisIO: get_svn

# =====================================================================
# Setup
@info("Please allow 20 minutes for all tests to execute.")
cd(dirname(pathof(SeisIO))*"/../test")
if isdir("SampleFiles") == false
  get_svn("https://github.com/jpjones76/SeisIO-TestData/trunk/SampleFiles", "SampleFiles")
end
include("local_restricted.jl")
include("test_helpers.jl")

# Announce test begin
test_start  = Dates.now()
ltn         = 48
printstyled(stdout,
  string(test_start, ": tests begin, path = ", path, ", has_restricted = ",
    has_restricted, ", keep_log = ", keep_log, ", keep_samples = ",
    keep_samples, "\n"),
  color=:light_green,
  bold=true)

# =====================================================================
# Run all tests
# grep "include(joinpath" runtests.jl | awk -F "(" '{print $3}' | awk -F "," {'print $1'}
for d in ["CoreUtils", "Types", "RandSeis", "Utils", "NativeIO", "DataFormats", "SEED", "Processing", "Nodal", "Quake", "Web"]
  ld = length(d)
  ll = div(ltn - ld - 2, 2)
  lr = ll + (isodd(ld) ? 1 : 0)
  printstyled(string("="^ll, " ", d, " ", "="^lr, "\n"), color=:cyan, bold=true)
  for i in readdir(joinpath(path, d))
    f = joinpath(d,i)
    if endswith(i, ".jl")
      printstyled(lpad(" "*f, ltn)*"\n", color=:cyan)
      write(out, string("\n\ntest ", f, "\n\n"))
      flush(out)
      include(f)
    end
  end
end

# =====================================================================
# Cleanup
include("cleanup.jl")
(keep_samples == true) || include("rm_samples.jl")
keep_log || safe_rm("runtests.log")

# Announce tests end
test_end = Dates.now()
δt = 0.001*(test_end-test_start).value
printstyled(string(test_end, ": tests end, elapsed time (mm:ss.μμμ) = ",
                   @sprintf("%02i", round(Int, div(δt, 60))), ":",
                   @sprintf("%06.3f", rem(δt, 60)), "\n"),
            color=:light_green,
            bold=true)
tut_file = realpath(path * "/../tutorial/install.jl")
ex_file = realpath(path * "/examples.jl")
printstyled("To run the interactive tutorial in a browser, execute: include(\"",
            tut_file, "\")\n", color=:cyan, bold=true)
printstyled("To run some data acquisition examples from the Julia prompt, ",
            "execute: include(\"", ex_file, "\")\n", color=:cyan)


text:
function convergenceratio(z, alpha, W, intervals, nexPercentage)

convergence = map(intervals) do interval

    P = x->phi(z,alpha,interval,x)

    mask = interval[1] .<= W .<= interval[2];
    W_inside  = W[mask]  |> P |> x->abs.(x) ;
    W_outside = W[.!mask] |> P |> x->abs.(x) |> x-> sort(x, order=Base.Order.Reverse);

    nev = length(W_inside)

    minval = minimum( W_inside )
    maxval = W_outside                  |>
        x-> Iterators.drop(x,ceil(Integer,nev*nexPercentage)) |>
        x-> Iterators.take(x,1)         |>
        collect

    if( length(maxval) == 0 )
        @printf("error\n");
        return NaN
    else
        return maxval[1] / minval
    end
end

end


text:
include("abstract_model2fit.jl")
include("visit.jl")

include("tagged_model.jl")
include("empty_model.jl")

include("model_sum.jl")
include("stacked_models.jl")

include("shared_parameters.jl")
include("mapped_parameters.jl")
include("transformed_parameters.jl")
include("const_parameters.jl")

include("peak/Peak.jl")

include("recalibration/Recalibration.jl")

include("solver_wrappers/Solver_Wrappers.jl")


text:
# ---
# title: 1024. Video Stitching
# id: problem1024
# author: Indigo
# date: 2022-04-11
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/video-stitching/description/>
# hidden: true
# ---
#
# You are given a series of video clips from a sporting event that lasted `T`
# seconds.  These video clips can be overlapping with each other and have varied
# lengths.
#
# Each video clip `clips[i]` is an interval: it starts at time `clips[i][0]` and
# ends at time `clips[i][1]`.  We can cut these clips into segments freely: for
# example, a clip `[0, 7]` can be cut into segments `[0, 1] + [1, 3] + [3, 7]`.
#
# Return the minimum number of clips needed so that we can cut the clips into
# segments that cover the entire sporting event (`[0, T]`).  If the task is
# impossible, return `-1`.
#
#
#
# **Example 1:**
#
#
#
#     Input: clips = [[0,2],[4,6],[8,10],[1,9],[1,5],[5,9]], T = 10
#     Output: 3
#     Explanation:
#     We take the clips [0,2], [8,10], [1,9]; a total of 3 clips.
#     Then, we can reconstruct the sporting event as follows:
#     We cut [1,9] into segments [1,2] + [2,8] + [8,9].
#     Now we have segments [0,2] + [2,8] + [8,10] which cover the sporting event [0, 10].
#
#
# **Example 2:**
#
#
#
#     Input: clips = [[0,1],[1,2]], T = 5
#     Output: -1
#     Explanation:
#     We can't cover [0,5] with only [0,1] and [1,2].
#
#
# **Example 3:**
#
#
#
#     Input: clips = [[0,1],[6,8],[0,2],[5,6],[0,4],[0,3],[6,7],[1,3],[4,7],[1,4],[2,5],[2,6],[3,4],[4,5],[5,7],[6,9]], T = 9
#     Output: 3
#     Explanation:
#     We can take clips [0,4], [4,7], and [6,9].
#
#
# **Example 4:**
#
#
#
#     Input: clips = [[0,4],[2,8]], T = 5
#     Output: 2
#     Explanation:
#     Notice you can have extra video after the event ends.
#
#
#
#
# **Constraints:**
#
#   * `1 <= clips.length <= 100`
#   * `0 <= clips[i][0] <= clips[i][1] <= 100`
#   * `0 <= T <= 100`
#
#
## @lc code=start
using LeetCode

function video_stitching(clips::Vector{Vector{Int}}, time::Int)
    max_next = fill(0, time)
    last, prev, res = 1, 1, 0
    for it in clips
        max_next[it[1] + 1] = max(max_next[it[1] + 1], it[2] + 1)
    end
    for i in 1:time
        last = max(last, max_next[i])
        i == last && return -1
        if i == prev
            res += 1
            prev = last
        end
    end
    return res
end
## @lc code=end