text:
"""
Functions to deal with various types of norms and seminorms
"""

using FastRounding
using IntervalOptimisation
using TaylorSeries: Taylor1
using SparseArrays: getcolptr
using .DynamicDefinition

using .DynamicDefinition: derivative

"""
'Absolute value' definition that returns mag(I) for an interval and abs(x) for a real
"""
abs_or_mag(x::Number) = abs(x)
abs_or_mag(x::Interval) = mag(x)

"""
Computes a rigorous upper bound for z*z'
"""
z_times_conjz(z::Complex) = square_round(abs_or_mag(real(z)), RoundUp) ⊕₊ square_round(abs_or_mag(imag(z)), RoundUp)
abs_or_mag(z::Complex) = sqrt_round(z_times_conjz(z), RoundUp)

"""
Certified upper bound to ||A|| (of specified NormKind)
"""
function opnormbound(::Type{L1}, A::AbstractVecOrMat{T}) where {T}
    # partly taken from JuliaLang's LinearAlgebra/src/generic.jl
    Tnorm = typeof(abs_or_mag(float(real(zero(T)))))
    Tsum = promote_type(Float64, Tnorm)
    nrm::Tsum = 0
    @inbounds begin
        for j = 1:size(A, 2)
            nrmj::Tsum = 0
            for i = 1:size(A, 1)
                nrmj = nrmj ⊕₊ abs_or_mag(A[i,j])
            end
            nrm = max(nrm,nrmj)
        end
    end
    return convert(Tnorm, nrm)
end

function opnormbound(::Type{Linf}, A::AbstractVecOrMat{T}) where {T}
    # partly taken from JuliaLang's LinearAlgebra/src/generic.jl
    Tnorm = typeof(abs_or_mag(float(real(zero(T)))))
    Tsum = promote_type(Float64, Tnorm)
    nrm::Tsum = 0
    @inbounds begin
        for i = 1:size(A, 1)
            nrmi::Tsum = 0
            for j = 1:size(A, 2)
                nrmi = nrmi ⊕₊ abs_or_mag(A[i,j])
            end
            nrm = max(nrm,nrmi)
        end
    end
    return convert(Tnorm, nrm)
end

"""
These functions compute a rigorous upper bound for the 2-norm of a vector;
we have a specialized version for complex numbers to avoid taking
the sqrt root and squaring again
"""
function opnormbound(::Type{L2}, v::Vector{T}) where {T<:Real}
    # partly taken from JuliaLang's LinearAlgebra/src/generic.jl
    Tnorm = typeof(abs_or_mag(float(real(zero(T)))))
    Tsum = promote_type(Float64, Tnorm)
    nrm::Tsum = 0
    @inbounds begin
        for j = 1:length(v)
            nrm = nrm ⊕₊ square_round(abs_or_mag(v[j]), RoundUp)
        end
    end
    return convert(Tnorm, sqrt_round(nrm, RoundUp))
end

function opnormbound(::Type{L2}, v::Vector{T}) where {T<:Complex}
    # partly taken from JuliaLang's LinearAlgebra/src/generic.jl
    Tnorm = typeof(abs_or_mag(float(real(zero(T)))))
    Tsum = promote_type(Float64, Tnorm)
    nrm::Tsum = 0
    @inbounds begin
        for j = 1:length(v)
            nrm = nrm ⊕₊ z_times_conjz(v[j])
        end
    end
    return convert(Tnorm, sqrt_round(nrm, RoundUp))
end

function opnormbound(::Type{L1}, A::SparseMatrixCSC) where {T}
    # partly taken from JuliaLang's Sparsearray/src/linalg.jl
    m, n = size(A)
    Tnorm = typeof(abs_or_mag(float(real(zero(eltype(A))))))
    Tsum = promote_type(Float64,Tnorm)
    nA::Tsum = 0
    @inbounds begin
        for j=1:n
            colSum::Tsum = 0
            for i = getcolptr(A)[j]:getcolptr(A)[j+1]-1
                colSum = colSum ⊕₊ abs_or_mag(nonzeros(A)[i])
            end
            nA = max(nA, colSum)
        end
    end
    return convert(Tnorm, nA)
end

function opnormbound(::Type{Linf}, A::SparseMatrixCSC) where {T}
    # partly taken from JuliaLang's Sparsearray/src/linalg.jl
    m, n = size(A)
    Tnorm = typeof(abs_or_mag(float(real(zero(eltype(A))))))
    Tsum = promote_type(Float64,Tnorm)
    rowSum = zeros(Tsum,m)
    @inbounds begin
        for i=1:length(nonzeros(A))
            rowSum[rowvals(A)[i]] = rowSum[rowvals(A)[i]] ⊕₊ abs_or_mag(nonzeros(A)[i])
        end
    end
    return convert(Tnorm, maximum(rowSum))
end

"""
Rigorous upper bound on a vector norm. Note that Linf, L1 are the "analyst's" norms
"""
normbound(N::Type{L1}, v::AbstractVector) = opnormbound(L1, v) ⊘₊ Float64(length(v), RoundDown)
normbound(N::Type{Linf}, v::AbstractVector) = opnormbound(Linf, v)

"""
Types to compute norms iteratively by "adding a column at a time".
"""
abstract type NormCacher{T} end

mutable struct NormCacherL1 <: NormCacher{L1}
    C::Float64
    function NormCacherL1(n)
        new(0)
    end
end
"""
Create a new NormCacher to compute the normbound of the empty matrix with n rows
"""
NormCacher{L1}(n) = NormCacherL1(n)

mutable struct NormCacherLinf <: NormCacher{Linf}
    C::Vector{Float64}
    function NormCacherLinf(n)
        new(zeros(n))
    end
end
NormCacher{Linf}(n) = NormCacherLinf(n)

"""
Update a NormCacher to add one column to the matrix it is computing a norm of.
This column may be affected by an error ε (in the same norm).
"""
function add_column!(Cacher::NormCacherL1, v::AbstractVector, ε::Float64)
    Cacher.C = max(Cacher.C, opnormbound(L1, v) ⊕₊ ε)
end

function add_column!(Cacher::NormCacherLinf, v::AbstractVector, ε::Float64)
    Cacher.C = Cacher.C .⊕₊ abs.(v) .⊕₊ ε
end

"""
Return the norm of the matrix the NormCacher is working on.
"""
function get_norm(Cacher::NormCacherL1)
    return Cacher.C
end

function get_norm(Cacher::NormCacherLinf)
    return maximum(Cacher.C)
end

"""
(A, B) = dfly(strongnorm, auxnorm, dynamic)

Constants (A, B) such that ||Lf||_s ≦ A||f||_s + B||f||_aux
"""
dfly(::Type{<:NormKind}, ::Type{<:NormKind}, ::Dynamic) = @error "Not implemented"

# I don't think this is used in production anymore
function dfly(::Type{TotalVariation}, ::Type{L1}, D::Dynamic)
    dist = max_distorsion(D)
    lam = expansivity(D)

    if !(abs(lam) < 1) # these are intervals, so this is *not* equal to abs(lam) >= 1.
        @error "The function is not expanding"
    end

    if is_full_branch(D)
        return lam.hi, dist.hi
    else
        if !(abs(lam) < 0.5)
            @error "Expansivity is insufficient to prove a DFLY. Try with an iterate."
        end
        # We need a way to estimate the branch widths
        @error "Not implemented"
    end
end

function dfly(::Type{TotalVariation}, ::Type{L1}, D::PwMap)
    dist = max_distorsion(D)
    lam = expansivity(D)
    vec = endpoints(D)
    disc = maximum(2/abs(vec[i]-vec[i+1]) for i in nbranches(D))

    if is_full_branch(D)
        if !(abs(lam) < 1) # these are intervals, so this is *not* equal to abs(lam) >= 1.
            @error "The function is not expanding"
        end
        return lam.hi, dist.hi
    else
        if !(abs(2*lam) < 1)
            @error "Expansivity is insufficient to prove a DFLY. Try with an iterate."
        end
        return (2*lam).hi, (dist + disc).hi
    end
end

function dfly(::Type{Lipschitz}, ::Type{L1}, D::Dynamic)
    # TODO: should assert that D is globally C2 instead, but we don't have that kind of infrastructure yet.
    @assert is_full_branch(D)

    dist = max_distorsion(D)
    lam = expansivity(D)

    return ((lam*(2*dist+1)).hi, (dist*(dist+1)).hi)
end


text:
reload("ProximalBase")
reload("CoordinateDescent")
reload("CovSel")

import CovSel, CoordinateDescent

srand(123)

p = 5
Sigmax = eye(p,p)
Sigmay = zeros(p,p)
rho = 0.7
for i=1:p
    for j=1:p
        Sigmay[i,j]=rho^abs(i-j)
    end
end

sqmy = sqrtm(Sigmay)

n = 1000
X = randn(n,p)
Y = randn(n,p) * sqmy;

hSx = cov(X)
hSy = cov(Y);

A = kron(hSy, hSx)
b = reshape(hSy-hSx, p*p)

S = sprand(p, p, 0.1)
S = S + S' / 2
indS = find(S)
# full(S)
# for i=indS
#   @show ind2sub(hSx, i)
# end

opt = CoordinateDescent.CDOptions(;maxIter=5000, optTol=1e-12, randomize=true)
x = CovSel.differencePrecisionRefit(Symmetric(hSx), Symmetric(hSy), indS, opt)
@show full(x)

x1 = zeros(p, p)
x1[indS] = A[indS, indS] \ b[indS]

@show full(x) - x1

ext: jl


text:
using DiffEqOperators, Test, LinearAlgebra

@testset "Second Order Value Check" begin
    ord = 2               # Order of approximation
    h = 0.025π
    N = Int(2*(π/h))
    x = -π:h:π-h
    L1  = CenteredDifference(1,ord,h,N);
    L2 = CenteredDifference(2,ord,h,N);

    Q = PeriodicBC(Float64)

    u0 = cos.(x)
    du_true = -cos.(x)

    # Explicit stencil for L2
    M = Matrix(Tridiagonal([1.0 for i in 1:N-1],
                           [-2.0 for i in 1:N],
                           [1.0 for i in 1:N-1]))

    # Insert periodic BC
    M[end,1] = 1.0
    M[1,end] = 1.0
    M = M/(h^2)
    @test M*u0 ≈ L2*Q*u0

    A = zeros(size(M))
    for i in 1:N, j in 1:N
        i == j && (A[i,j] = -0.5)
        abs(i - j) == 2 && (A[i,j] = 0.25)
    end
    A[end-1,1] = 0.25
    A[end,2] = 0.25
    A[1,end-1] = 0.25
    A[2,end] = 0.25
    A = A/(h^2)

    @test A*u0 ≈ L1*Q*(L1*Q*u0)

    κ(x) = 1.0
    tmp1 = similar(x)
    tmp2 = similar(x)
    function f(t,u,du)
        mul!(tmp1,L1,Q*u)
        @. tmp2 = κ(x)*tmp1
        mul!(du,L1,Q*tmp2)
    end

    du = similar(u0)
    f(0,u0,du)
    du ≈ A*u0

    error = du_true - du
    du2 = L2*Q*u0
    error2 = du_true - du2

    @test maximum(error) < 0.004
    @test maximum(error2) < 0.001
end


text:
"""
This is the main Julia file for the replication exercise.
First, it loads the data (estimated parameters, data moments, etc.) and guess
parameters.
Second, it simulates a large sample of firms [module DGPsetup].
Third, it finds optimal parameters according to GMM:
    - solve lower and upper bound of optimal set of sourcing countries with Jia
        algorithm [module JiaAlgorithm];
    - solve optimal set of sourcing countries [module JiaAlgorithm];
    - compute moments and objective function [module gmmObjectiveFun];
    - solve.

"""

## Set directory
cd("C:\\Users\\asus\\Desktop\\Giulia\\UBC\\Year2\\567 - Empirical IO\\AFT2017-Replication\\Code")

## Use/Import modules
using Main.DGPsetup, Main.JiaAlgorithm, Main.gmmObjectiveFun

## Load packages
using LinearAlgebra, Random, Distributions, Statistics, DataFrames
using CSV, XLSX, StatsBase, Optim, BlackBoxOptim

## Set seeds
Random.seed!(6)

## Assign country ID
chn_ctr_ind  = 59
ger_ctr_ind  = 28
can_ctr_ind  = 2
gbr_ctr_ind  = 22
mex_ctr_ind  = 3
twn_ctr_ind  = 62
ita_ctr_ind  = 39
jpn_ctr_ind  = 63
fra_ctr_ind  = 27
kor_ctr_ind  = 60

## Load input public data
df = CSV.read("emp2dataformatlabv3.csv", DataFrame)
df_additional = CSV.read("additionalparamandmom.csv", DataFrame)

names(df)
names(df_additional)

## Compute data moments
fixedeffect = df.fe
nfirms = df.firm_count_rounded
aggimports = df.inputs_rounded *10
distrw = df.distw
contiguous = df.contig
comlang = df.comlang_off
corrup = df.control_of_corruption
gdpctry = df.gdp
ruleoflaw = df.rule_of_law
rpsh = df.rpshare
exportcost = df.exp_cost
intserver = df.intserv
N = size(df.code_TD,1)

US_median_dom_input = 0.5675603
aggCN_imp_share_growth = ((0.1366 / 0.0492) - 1) * 100
# 2007 share of imports over 1907 share of imports
shareimp_salesq1 = 0.058
shareimp_salesq2 = 0.112
shareimp_salesq3 = 0.245
shareimp_salesq4 = 0.620

df_additional.param_name
df_additional.param_value

σ = df_additional.param_value[1,1]
θ = df_additional.param_value[2,1]
nfirmstot = df_additional.param_value[3,1]
nimportingfirms = df_additional.param_value[4,1]

κ = 4.25
σ_tilde = (1 ./ σ) * (σ ./ (σ - 1)).^(1 - σ)
ξ = exp.(fixedeffect)   # sourcing potential
# Sort ξ descending and get ranking (excluding US)
rank_ξ = ordinalrank(ξ[2:end,:]; rev = true)
rank_ξ = vec(rank_ξ)

## Initial guesses
B_guess = 0.05  # but this is smaller than Guess_LB; why?
fc_mean_guess = [.1; 0.5; 0.9; 0.1]
fc_disp_guess = 1
δ_guess = [B_guess; fc_mean_guess; fc_disp_guess]

lb = 1e-006 * ones(size(δ_guess))
ub = 10 * ones(size(δ_guess))
ub[6] = 6.00

## Allow cores to start from different guesses
# (taken from the original code, but I will not run parallel cores)
# Initial guesses will be uniformly distributed in a bigger interval:
MS = 10
guess_lb = [0.1  ; 0.001 ; 0.1 ; .3 ; .05 ; .5  ]
guess_ub = [0.15 ; 0.2   ; 0.8 ; 1  ; .8  ; 1.5 ]

δ_guess_all = ones(MS, size(δ_guess,1))
δ_guess_all[1,:] = δ_guess'
δ_guess_all[2:MS,:] = (guess_ub - guess_lb)' .* rand(MS-1, size(δ_guess,1)) +
                    + repeat(guess_lb',MS-1,1)

# In AFT2017 codes, they use parallel programming with 10 cores, each using a
# different guess for delta:
if MS > 3
    δ_guess_all[4,:] = [0.1250; 0.0300; 0.3500; 0.8250; 0.2750; 1.2500]
end

if MS > 4
    δ_guess_all[5,:] = [0.1240; 0.0282; 0.3524; 0.7704; 0.3609; 1.1641]
end

if MS > 5
    δ_guess_all[6,:] = [0.1236; 0.0230; 0.1923; 0.8647; 0.3917; 0.9370]
end

if MS > 6
    δ_guess_all[7,:] = [0.1260; 0.0110; 0.1000; 0.5111; 0.2000; 1.2912]
end

if MS > 7
    δ_guess_all[8,:] = [0.1260; 0.0110; 0.1000; 0.8000; 0.1500; 0.6000]
end

δ_hat_all = -100 * ones(MS, length(δ_guess))
fval_all = 99999 * ones(MS,1)
exitflag_all = -100 * ones(MS,1)

post_estimation = 0

## Simulate firms
S, prod_draw_uniform, weights_prod, fc_shock_randn, num_rand_checks,
    rand_check_matrix = simulatefirms(N; S_fixed = 1800)

## Define the problem as a function of δ_guess
function objectivefunction(δ_guess, σ, θ, κ, distrw, comlang, corrup, N, ξ,
    S, prod_draw_uniform, weights_prod, fc_shock_randn, num_rand_checks,
    rand_check_matrix, nimportingfirms, nfirms, nfirmstot, shareimp_salesq1,
    shareimp_salesq2, US_median_dom_input)

    ## Find optimal set of sourcing countries
    # Set paramaters
    my_exp = (σ-1)/θ
    ϕ_σ_B = δ_guess[1] * ((1 .- prod_draw_uniform).^(-1/κ)).^(σ-1)
    fc_mean = δ_guess[2] .* ((distrw').^δ_guess[3]) .* δ_guess[4].^(comlang') .* exp.(-δ_guess[5] .* corrup')
    temp = fc_shock_randn.*δ_guess[6]
    temp2 = 709 .* ones(size(temp,1), size(temp,2))
    fc = fc_mean .* exp.(min.(temp, temp2))
    fc[:,1] = zeros(size(fc,1),1)   # remember 0 cost of domestic sourcing (US is 1st)

    # do loop for all simulated firms
    Z = 1.0*ones(S,N)
    gap_bounds = 1.0*ones(S,1)

    source_start_lb, source_check_lb = lowerbound_setup(N, ξ, my_exp)
    source_start_ub, source_check_ub = upperbound_setup(N, ξ, my_exp)

    for firm in 1:S
        #print("Firm number:")
        #println("$firm")
        Z_lb = lowerbound(source_start_lb, source_check_lb, ϕ_σ_B, fc, N, ξ, my_exp, firm)
        Z_ub = upperbound(source_start_ub, source_check_ub, ϕ_σ_B, fc, N, ξ, my_exp, firm)
        Z = optimalset(Z, gap_bounds, firm, Z_lb, Z_ub, S, N, num_rand_checks,
            rand_check_matrix, fc, ξ, my_exp, ϕ_σ_B)
    end

    ## Compute the gmm objective function
    sales, input_p_mat = SalesAndInput(Z, ξ, σ, ϕ_σ_B, my_exp)
    valuetominimize = gmmobjective(sales, input_p_mat, weights_prod,
        nimportingfirms, nfirms, nfirmstot, shareimp_salesq1, shareimp_salesq2,
        US_median_dom_input, N)

    println("$valuetominimize")

    return valuetominimize
end

# check it works
objectivefunction(δ_guess, σ, θ, κ, distrw, comlang, corrup, N, ξ,
    S, prod_draw_uniform, weights_prod, fc_shock_randn, num_rand_checks,
    rand_check_matrix, nimportingfirms, nfirms, nfirmstot, shareimp_salesq1,
    shareimp_salesq2, US_median_dom_input)

## Solve using Optim
# Start with guess very close to solution
δ0 = [0.120; 0.020; 0.190; 0.870; -0.390; 0.930]
δ_star_unbounded = optimize(δ->objectivefunction(δ, σ, θ, κ, distrw, comlang, corrup,
        N, ξ, S, prod_draw_uniform, weights_prod, fc_shock_randn, num_rand_checks,
        rand_check_matrix, nimportingfirms, nfirms, nfirmstot, shareimp_salesq1,
        shareimp_salesq2, US_median_dom_input), δ0, Optim.Options(g_tol=0.00001))

# Do loop with same guesses as original paper
for round in 1:10
    guess = δ_guess_all[round,:]
    aux = optimize(δ->objectivefunction(δ, σ, θ, κ, distrw, comlang, corrup,
            N, ξ, S, prod_draw_uniform, weights_prod, fc_shock_randn, num_rand_checks,
            rand_check_matrix, nimportingfirms, nfirms, nfirmstot, shareimp_salesq1,
            shareimp_salesq2, US_median_dom_input), guess, Optim.Options(g_tol=0.00001))
    δ_hat_all[round,:] = aux.minimizer'
end

## Solve using BlackBoxOptim with same search range for all values of δ
# It has been running for 13 hours without finding a solution (I interrupted Julia and put
# a maximum number of steps).
# I have copied the search range from the original code, but one of the estimates reported in the
# paper is negative! How is that possible, if they do minimization with boundaries [0;10]?
δ_star_bound1 = bboptimize(δ->objectivefunction(δ, σ, θ, κ, distrw, comlang, corrup,
        N, ξ, S, prod_draw_uniform, weights_prod, fc_shock_randn, num_rand_checks,
        rand_check_matrix, nimportingfirms, nfirms, nfirmstot, shareimp_salesq1,
        shareimp_salesq2, US_median_dom_input); SearchRange = (1e-006, 10.0),
        NumDimensions = 6, MaxSteps = 1000)

# Try using a range that respects the fact that one estimate is negative.
# Alternative to max number of steps: use targetfitness = 0.7.
δ_star_bound2 = bboptimize(δ->objectivefunction(δ, σ, θ, κ, distrw, comlang, corrup,
        N, ξ, S, prod_draw_uniform, weights_prod, fc_shock_randn, num_rand_checks,
        rand_check_matrix, nimportingfirms, nfirms, nfirmstot, shareimp_salesq1,
        shareimp_salesq2, US_median_dom_input); SearchRange = (-0.5, 10.0),
        NumDimensions = 6, MaxSteps = 1000)


text:
export with, Parameter, LazyParameter

mutable struct Parameter{T}
  value::T
end

# This no longer works in 1.0
#(p::Parameter{T}){T}()::T = p.value
#(p::Parameter{T}){T}(newvalue::T) = p.value = newvalue

(p::Parameter)() = p.value
(p::Parameter)(newvalue) = p.value = newvalue

#=
with{T}(f, p::Parameter{T}, newvalue::T) where T =
  begin
    oldvalue, p.value = p.value, newvalue
    try
      f()
    finally
      p.value = oldvalue
    end
  end
=#

# A more generic version (presumably, compatible with the previous one)
function with(f, p, newvalue)
  oldvalue = p()
  p(newvalue)
  try
    f()
  finally
    p(oldvalue)
  end
end

with(f, p, newvalue, others...) =
  with(p, newvalue) do
    with(f, others...)
  end

mutable struct LazyParameter{T}
  initializer::Function #This should be a more specific type: None->T
  value::Union{T, Nothing}
end

LazyParameter(T::DataType, initializer::Function) = LazyParameter{T}(initializer, nothing)

# This no longer works in 1.0
#(p::LazyParameter{T}){T}()::T = p.value == nothing ? (p.value = p.initializer()) : get(p.value)
#(p::LazyParameter{T}){T}(newvalue::T) = p.value = newvalue

(p::LazyParameter)() = p.value === nothing ? (p.value = p.initializer()) : p.value
(p::LazyParameter)(newvalue) = p.value = newvalue

import Base.reset
reset(p::LazyParameter{T}) where {T} = p.value = nothing


text:
using IntcodeVM, DelimitedFiles, Memoize

@memoize function gravity(x::Int,y::Int)
    return last(last(IntcodeVM.run(code, inputs=[x,y])))
end

const code = vec(readdlm("input.txt", ',', Int, '\n'))
tractor_beam = [gravity(x,y) for x=0:49, y=0:49]

print("Day 19 Part 1: ")
println(count(==(1), tractor_beam))

# It would be computationally crazy to try to brute force this.
# This looks like something solvable through linear regression.

# We can perform a linear search for values xmin <= x <= xmax where gravity(x,y)==1
# and xmin + 100 == xmax.
# For such values, we don't actually need to search the entire square,
# just (xmin, y), (xmax, y), (xmin, y + 100).
# In fact, that might be lightweight enough to compute through
# brute force.

# We can look at the graphic and see that the solution is going to align
# to the right. So, we should probably find a pure function y=f(xmax) to
# predict where the right side of the triangle is.
# So there you have it: run linear regression along right side.

@memoize function get_bounds(y)
    y < 6 && return (0,0)

    # The current xmin value can never be less than xmin value at y-1.
    xmin = first(get_bounds(y-1))
    while gravity(xmin, y) == 0
        xmin += 1
    end

    # Same for xmax, we can never be to the left of the right bound at y-1.
    xmax = max(xmin, last(get_bounds(y-1)))
    while gravity(xmax + 1, y) == 1
        xmax += 1
    end
    return xmin, xmax
end

function print_beam(a)
    xmax,ymax = size(a)
    for y=1:ymax
        for x=1:xmax
            print(a[x,y])
        end
        # 1-indexed languages are so annoying sometimes.
        y >= 7 && print(" y=$(y-1): bounds are $(get_bounds(y - 1))")
        println()
    end
end

print_beam(tractor_beam)

# This is a linear programming question.
# We run a linear regression to predict what the left and right sides of our bounds are.
#using DataFrames, GLM;
#data_left = DataFrame(X=6:49, Y=[first(get_bounds(i)) for i=6:49])
#ols_left = lm(@formula(Y ~ X), data_left)
#data_right = DataFrame(X=6:49, Y=[last(get_bounds(i)) for i=6:49])
#ols_right = lm(@formula(Y ~ X), data_right)

# We can see from here that the right prediction is extremely accurate.
# To my disappointment, I now realize the formula is just right_bound = row - 1
#hcat([get_bounds(i) for i=6:49], predict(ols_left), predict(ols_right))
# D'oh, something goes wrong at y=216 and y=217! The right bound does not linearly increase between these rows.

function find_square(square_size::Int, search_space::UnitRange{Int})
    # Ahh, my old friend off-by-one errors.
    # We have to include the current location x or y in a square of size n.
    dy = square_size - 1
    for row=search_space
        b1 = get_bounds(row)
        b2 = get_bounds(row + dy)
        #println("Bounds(row $(row)) = $(b1), Bounds(row $(row + dy)) = $(b2)")
        #println(b1[2] - b2[1])
        if b1[2] - b2[1] == dy
            return (b2[1], row)
        end
    end
    return (-1,-1)
end

(x,y) = find_square(100, 1:10000)
println("Day 19 Part 2: $(x * 10000 + y)")


text:
export dclr

function dclr(b::AbstractArray{ET}, r::Int, K::Int) where ET<:Number
# Performs RCAM on a single matricised monochromatic slice
# Calls: minFrobQR

    # Get worker IDs
    wids = sort(workers())
    nworkers = length(wids)

    if nworkers < 3
        error(string("DCLR method requires atleast three workers, got ", nworkers))
    end

    # Dedicate two workers to distribute/store L and R
    ownerL = wids[1]
    ownerR = wids[2]

    # Dedicate remaining workers to as solvers
    solvers = wids[3:nworkers]
    nsolvers = length(solvers)

    # Give this info to every solver
    [sendto(pid, ownerL = ownerL, ownerR = ownerR) for pid in procs()] #change to procs()[solvers]

    # Get size of data inorder to generate L0
    m::Int,n::Int = size(b)

    #Create L and R on their owners
    sendto(ownerL, L=collect(ET, randn(m,r)))
    sendto(ownerR, R=zeros(ET,n,r))

    # Find non-zero inds for each col, reduce data into collections
    indn = (1:n)::UnitRange{Int64}

    # This global declaration is a workaround to open issue #15276
    indnzd = distribute([findall(!iszero,b[:,i]) for i::Int in indn],
                        procs = procs()[solvers])
    bsubd = distribute([b[indnzd[i],i] for i in indn],
                        procs = procs()[solvers])

    # Find non zero inds for each row, reduce data into collections
    indm = (1:m)::UnitRange{Int64}
    indnz_trand = distribute([findall(!iszero,b[i,:]) for i::Int in indm],
                        procs = procs()[solvers])
    bsub_trand = distribute([b[i,indnz_trand[i]] for i in indm],
                        procs = procs()[solvers])

    # Clear data
    b = []

    # Begin Alterations
    for k = 1:K

        @time @sync for pid in procs()[solvers]
            @async fetch(@spawnat pid altupdateR!(ownerL, ownerR, bsubd, indnzd))
        end

        @time @sync for pid in procs()[solvers]
            @async fetch(@spawnat pid altupdateL!(ownerL, ownerR, bsub_trand,indnz_trand))
        end

    end

   return ownerL, ownerR

end

function altupdateL!(ownerL::Int, ownerR::Int, bsubd::AbstractArray{AT}, indnzd) where AT

    ET = eltype(AT)

    #Get info about local part of data
    m = localindices(bsubd)[1]

    # Get updated version of R
    R_tmp = fetch(@spawnat(ownerR, Main.R))

    # Update the ith row of Rd
    for i in m

        local_i = i - m[1] + 1

        # Ask for required part of L
        req_rows = DistributedArrays.localpart(indnzd)[local_i]

        # Solve for what this row or R will be updated to
        L_update = collect(ET,conj(pminFrobRsubChol(R_tmp[req_rows,:],
                            conj(DistributedArrays.localpart(bsubd)[local_i]))))

        # Create a future for this result
        L_future = Future(ownerL)
        put!(L_future, L_update)

        # Update the i-th row of R with the future
        @spawnat(ownerL, setindex!(Main.L, fetch(L_future), i, :))

    end

end

function altupdateR!(ownerL::Int, ownerR::Int, bsubd::AbstractArray{AT}, indnzd) where AT

    ET = eltype(AT)

    #Get info about local part of data
    m = localindices(bsubd)[1]

    # Get updated version of L
    L_tmp = fetch(@spawnat(ownerL, Main.L))

    # Update the ith row of Rd
    for i in m

        local_i = i - m[1] + 1

        # Ask for required part of L
        req_rows = DistributedArrays.localpart(indnzd)[local_i]

        # Solve for what this row or R will be updated to
        R_update = collect(ET,conj(pminFrobRsubChol(L_tmp[req_rows,:],
                            DistributedArrays.localpart(bsubd)[local_i])))

        # Create a future for this result
        R_future = Future(ownerR)
        put!(R_future, R_update)

        # Update the i-th row of R with the future
        @spawnat(ownerR, setindex!(Main.R, fetch(R_future), i, :))

    end

end

function pminFrobRsubChol(L2::AbstractArray{ET},bsub::AbstractArray{ET}) where ET<:Number
# Solves a single row, meant to be mapped over a DArray

    bsub = L2'*bsub
    L2 =L2'*L2

    #fL2 = chol(L2)
    #v = backsub(fL2,fwdsub((fL2)',bsub))
    v = L2\bsub

    return v
end
export pminFrobRsubChol

function backsub(A::AbstractArray{ET}, y::AbstractArray{ET}) where ET<:Number

    m = size(A,2)
    x = zeros(ET,m,1)
    x[m] = y[m]/A[m,m]
    for i = m-1:-1:1
        x[i] = (y[i]-dot(A[i,i+1:m], x[i+1:m]))/A[i,i]
    end

    return x
end

function fwdsub(A::AbstractArray{ET},y::AbstractArray{ET}) where ET<:Number

    m = size(A,2)
    x = zeros(ET,m,1)
    x[1] = y[1]/A[1,1]
    for i = 2:1:m
        x[i] = (y[i]-dot(A[i,i-1:-1:1], x[i-1:-1:1]))/A[i,i]
    end

    return x
end


text:
using MiraculousNN
using Test

using MLDatasets

@testset "Basic neural net workflow" begin
    ninstances = 1000
    nclasses = 10

    x, y_orig = MNIST.traintensor(Float64, 1:ninstances), MNIST.trainlabels(1:ninstances)

    x = reshape(x, (:, ninstances))
    y = MiraculousNN.utils.to_categorical(y_orig, nclasses)

    input_size = size(x, 1)
    output_size = size(y, 1)

    nn_params = MiraculousNN._init_layers(typeof(x), [input_size, 20, 10, output_size])

    @testset "init_layers" begin
        @test keys(nn_params) |> length == 6

        @test size(nn_params["W1"]) == (20, 784)
        @test size(nn_params["B1"]) == (20, 1)

        @test size(nn_params["W2"]) == (10, 20)
        @test size(nn_params["B2"]) == (10, 1)

        @test size(nn_params["W3"]) == (output_size, 10)
        @test size(nn_params["B3"]) == (output_size, 1)

        @test nn_params["W1"] != nn_params["W2"]
        @test nn_params["W1"] != nn_params["W3"]

        @test nn_params["B1"] != nn_params["B2"]
        @test nn_params["B1"] != nn_params["B3"]
    end

    nn_values = MiraculousNN._forward(x, nn_params)

    @testset "forward" begin
        @test keys(nn_values) |> length == 6

        @test size(nn_values["Z1"]) == (20, ninstances)
        @test size(nn_values["A1"]) == (20, ninstances)

        @test size(nn_values["Z2"]) == (10, ninstances)
        @test size(nn_values["A2"]) == (10, ninstances)

        @test size(nn_values["Z3"]) == (output_size, ninstances)
        @test size(nn_values["A3"]) == (output_size, ninstances)
    end

    grads = MiraculousNN._backward(nn_params, nn_values, x, y)

    @testset "backward" begin
        @test keys(grads) |> length == 6

        @test size(grads["W1"]) == (20, 784)
        @test size(grads["B1"]) == (20, 1)

        @test size(grads["W2"]) == (10, 20)
        @test size(grads["B2"]) == (10, 1)

        @test size(grads["W3"]) == (output_size, 10)
        @test size(grads["B3"]) == (output_size, 1)
    end

    MiraculousNN._update_weights!(nn_params, grads, 0.1)

    @testset "update_weights" begin
        @test keys(nn_params) |> length == 6

        @test size(nn_params["W1"]) == (20, 784)
        @test size(nn_params["B1"]) == (20, 1)

        @test size(nn_params["W2"]) == (10, 20)
        @test size(nn_params["B2"]) == (10, 1)

        @test size(nn_params["W3"]) == (output_size, 10)
        @test size(nn_params["B3"]) == (output_size, 1)

        @test nn_params["W1"] != nn_params["W2"]
        @test nn_params["W1"] != nn_params["W3"]

        @test nn_params["B1"] != nn_params["B2"]
        @test nn_params["B1"] != nn_params["B3"]
    end
end

max_stars_repo_name: prenc/MiraculousNN

max_issues_repo_name: prenc/MiraculousNN

max_forks_repo_name: prenc/MiraculousNN


text:
"""
    Map <: AbstractDict{MOI.VariableIndex, AbstractBridge}

Mapping between bridged variables and the bridge that bridged the variable.
"""
mutable struct Map <: AbstractDict{MOI.VariableIndex,AbstractBridge}
    # Bridged constrained variables
    # `i` ->  `0`: `VariableIndex(-i)` was added with `add_constrained_variable`.
    # `i` -> `-j`: `VariableIndex(-i)` was the first variable of
    #              `add_constrained_variables` with a set of dimension `j`.
    # `i` ->  `j`: `VariableIndex(-i)` was the `j`th  variable of
    #             ` add_constrained_variables`.
    info::Vector{Int64}
    # `i` ->  `-1`: `VariableIndex(-i)` was deleted.
    # `i` ->  `0`: `VariableIndex(-i)` was added with `add_constrained_variable`.
    # `i` ->  `j`: `VariableIndex(-i)` is the `j`th  variable of a constrained
    #               vector of variables, taking deletion into account.
    index_in_vector::Vector{Int64}
    # `i` -> `bridge`: `VariableIndex(-i)` was bridged by `bridge`.
    bridges::Vector{Union{Nothing,AbstractBridge}}
    sets::Vector{Union{Nothing,Type}}
    # If `nothing`, it cannot be computed because some bridges does not support it
    unbridged_function::Union{
        Nothing,
        Dict{MOI.VariableIndex,Tuple{Int64,MOI.AbstractScalarFunction}},
    }
    # Bridge that created this bridge, 0 if it is no bridge.
    parent_index::Vector{Int64}
    # Current bridge, 0 otherwise.
    current_context::Int64
    # Context of constraint bridged by constraint bridges
    constraint_context::Dict{MOI.ConstraintIndex,Int64}
end

function Map()
    return Map(
        Int64[],
        Int64[],
        Union{Nothing,AbstractBridge}[],
        Union{Nothing,Type}[],
        Dict{MOI.VariableIndex,MOI.AbstractScalarFunction}(),
        Int64[],
        0,
        Dict{MOI.ConstraintIndex,Int64}(),
    )
end

# Implementation of `AbstractDict` interface.

Base.isempty(map::Map) = all(bridge -> bridge === nothing, map.bridges)

function Base.empty!(map::Map)
    empty!(map.info)
    empty!(map.index_in_vector)
    empty!(map.bridges)
    empty!(map.sets)
    if map.unbridged_function === nothing
        map.unbridged_function =
            Dict{MOI.VariableIndex,Tuple{Int64,MOI.AbstractScalarFunction}}()
    else
        empty!(map.unbridged_function)
    end
    empty!(map.parent_index)
    map.current_context = 0
    empty!(map.constraint_context)
    return map
end

function bridge_index(map::Map, vi::MOI.VariableIndex)
    index = map.info[-vi.value]
    if index ≤ 0
        return -vi.value
    else
        return -vi.value - index + 1
    end
end

function Base.haskey(map::Map, vi::MOI.VariableIndex)
    return -length(map.bridges) ≤ vi.value ≤ -1 &&
           map.bridges[bridge_index(map, vi)] !== nothing &&
           map.index_in_vector[-vi.value] != -1
end

function Base.getindex(map::Map, vi::MOI.VariableIndex)
    return map.bridges[bridge_index(map, vi)]
end

function Base.delete!(map::Map, vi::MOI.VariableIndex)
    if iszero(map.info[-vi.value])
        # Delete scalar variable
        map.bridges[bridge_index(map, vi)] = nothing
        map.sets[bridge_index(map, vi)] = nothing
    elseif has_keys(map, [vi])
        # Delete whole vector
        delete!(map, [vi])
    else
        # Delete variable in vector and resize vector
        map.info[bridge_index(map, vi)] += 1
        for i in (-vi.value):length(map.index_in_vector)
            if map.index_in_vector[i] == -1
                continue
            elseif bridge_index(map, vi) !=
                   bridge_index(map, MOI.VariableIndex(-i))
                break
            end
            map.index_in_vector[i] -= 1
        end
    end
    map.index_in_vector[-vi.value] = -1
    return map
end

function Base.delete!(map::Map, vis::Vector{MOI.VariableIndex})
    if !has_keys(map, vis)
        throw(
            ArgumentError(
                "`$vis` is not a valid key vector as returned by `add_keys_for_bridge`.",
            ),
        )
    end
    for vi in vis
        map.index_in_vector[-vi.value] = -1
    end
    map.bridges[bridge_index(map, first(vis))] = nothing
    map.sets[bridge_index(map, first(vis))] = nothing
    return map
end

function Base.keys(map::Map)
    return Base.Iterators.Filter(
        vi -> haskey(map, vi),
        MOIU.LazyMap{MOI.VariableIndex}(
            i -> MOI.VariableIndex(-i),
            eachindex(map.bridges),
        ),
    )
end

Base.length(map::Map) = count(bridge -> bridge !== nothing, map.bridges)

function number_of_variables(map::Map)
    num = 0
    for i in eachindex(map.bridges)
        if map.bridges[i] !== nothing
            if iszero(map.info[i])
                num += 1
            else
                num += length_of_vector_of_variables(map, MOI.VariableIndex(-i))
            end
        end
        count(bridge -> bridge !== nothing, map.bridges)
    end
    return num
end

function Base.values(map::Map)
    # We don't use `filter` as it would compute the resulting array which
    # is not necessary if the caller just wants to iterater over `values`.
    return Base.Iterators.Filter(bridge -> bridge !== nothing, map.bridges)
end

function Base.iterate(map::Map, state = 1)
    while state ≤ length(map.bridges) && map.bridges[state] === nothing
        state += 1
    end
    if state > length(map.bridges)
        return nothing
    else
        return MOI.VariableIndex(-state) => map.bridges[state], state + 1
    end
end

# Custom interface for information needed by `AbstractBridgeOptimizer`s that is
# not part of the `AbstractDict` interface.

"""
    constrained_set(map::Map, vi::MOI.VariableIndex)

Return the set type in which the bridged variable `vi` was added when it was
bridged.
"""
function constrained_set(map::Map, vi::MOI.VariableIndex)
    return map.sets[bridge_index(map, vi)]
end

"""
    number_with_set(map::Map, S::Type{<:MOI.AbstractSet})

Return the number of bridged variables in `S`. Note that if `S` is a vector set,
bridging a vector of `n` variables only counts as 1.
"""
function number_with_set(map::Map, S::Type{<:MOI.AbstractSet})
    return count(isequal(S), map.sets)
end

function constraint(map::Map, vi::MOI.VariableIndex)
    S = constrained_set(map, vi)
    F = MOIU.variable_function_type(S)
    return MOI.ConstraintIndex{F,S}(-bridge_index(map, vi))
end

"""
    constraints_with_set(map::Map, S::Type{<:MOI.AbstractSet})

Return the list of constraints corresponding to bridged variables in `S`.
"""
function constraints_with_set(map::Map, S::Type{<:MOI.AbstractSet})
    F = MOIU.variable_function_type(S)
    return [
        MOI.ConstraintIndex{F,S}(-i) for
        i in eachindex(map.sets) if map.sets[i] == S
    ]
end

"""
    list_of_constraint_types(map::Map)

Return a list of all the different types `(F, S)` of `F`-in-`S` constraints in
`map`.
"""
function list_of_constraint_types(map::Map)
    list = Set{Tuple{Type,Type}}()
    for i in eachindex(map.bridges)
        if map.bridges[i] !== nothing
            S = map.sets[i]
            if S != MOI.Reals
                F = MOIU.variable_function_type(S)
                push!(list, (F, S))
            end
        end
    end
    return list
end

"""
    has_keys(map::Map, vis::Vector{MOI.VariableIndex})::Bool

Return a `Bool` indicating whether `vis` was returned by
[`add_keys_for_bridge`](@ref) and has not been deleted yet.
"""
function has_keys(map::Map, vis::Vector{MOI.VariableIndex})
    return isempty(vis) || (
        length_of_vector_of_variables(map, first(vis)) == length(vis) &&
        all(
            vi -> bridge_index(map, vi) == bridge_index(map, first(vis)),
            vis,
        ) &&
        all(vi -> haskey(map, vi), vis) &&
        all(i -> vis[i].value < vis[i-1].value, 2:length(vis))
    )
end

"""
    length_of_vector_of_variables(map::Map, vi::MOI.VariableIndex)

If `vi` was bridged in a scalar set, it returns 0. Otherwise, it
returns the dimension of the set.
"""
function length_of_vector_of_variables(map::Map, vi::MOI.VariableIndex)
    return -map.info[bridge_index(map, vi)]
end

"""
    index_in_vector_of_variables(map::Map, vi::MOI.VariableIndex)::MOIB.IndexInVector

Return the index of `vi` in the vector of variables in which it was bridged.
"""
function index_in_vector_of_variables(map::Map, vi::MOI.VariableIndex)
    return MOIB.IndexInVector(map.index_in_vector[-vi.value])
end

"""
    has_bridges(map::Map)::Bool

Return a `Bool` indicating whether any bridge was added yet. Note that it
returns `false` even if all bridges were deleted while `isempty` would return
`true`. It is computed in `O(1)` while `isempty` needs `O(n)` hence it is used
by [`MathOptInterface.Bridges.AbstractBridgeOptimizer`](@ref) to shortcut
operations in case variable bridges are not used.
"""
has_bridges(map::Map) = !isempty(map.info)

"""
    add_key_for_bridge(map::Map, bridge_fun::Function,
                       set::MOI.AbstractScalarSet)

Create a new variable index `vi`, store the mapping `vi => bridge` and
associate `vi` to `typeof(set)`. It returns a tuple with `vi` and the
constraint index
`MOI.ConstraintIndex{MOI.VariableIndex, typeof(set)}(vi.value)`.
"""
function add_key_for_bridge(
    map::Map,
    bridge_fun::Function,
    set::MOI.AbstractScalarSet,
)
    push!(map.parent_index, map.current_context)
    bridge_index = Int64(length(map.parent_index))
    push!(map.info, 0)
    push!(map.index_in_vector, 0)
    push!(map.bridges, nothing)
    push!(map.sets, typeof(set))
    map.bridges[bridge_index] = call_in_context(map, bridge_index, bridge_fun)
    index = -bridge_index
    variable = MOI.VariableIndex(index)
    if map.unbridged_function !== nothing
        mappings = unbridged_map(map.bridges[bridge_index], variable)
        if mappings === nothing
            map.unbridged_function = nothing
        else
            for mapping in mappings
                push!(
                    map.unbridged_function,
                    mapping.first => (bridge_index, mapping.second),
                )
            end
        end
    end
    return variable, MOI.ConstraintIndex{MOI.VariableIndex,typeof(set)}(index)
end

"""
    add_keys_for_bridge(map::Map, bridge_fun::Function,
                        set::MOI.AbstractVectorSet)

Create vector of variable indices `variables`, stores the mapping
`vi => bridge` for each `vi ∈ variables` and associate `variables` to
`typeof(set)`. It returns a tuple with `variables` and the constraint index
`MOI.ConstraintIndex{MOI.VectorOfVariables, typeof(set)}(first(variables).value)`.
"""
function add_keys_for_bridge(
    map::Map,
    bridge_fun::Function,
    set::MOI.AbstractVectorSet,
)
    if iszero(MOI.dimension(set))
        return MOI.VariableIndex[],
        MOI.ConstraintIndex{MOI.VectorOfVariables,typeof(set)}(0)
    end
    push!(map.parent_index, map.current_context)
    bridge_index = Int64(length(map.parent_index))
    push!(map.info, -MOI.dimension(set))
    push!(map.index_in_vector, 1)
    push!(map.bridges, nothing)
    push!(map.sets, typeof(set))
    for i in 2:MOI.dimension(set)
        push!(map.parent_index, 0)
        push!(map.info, i)
        push!(map.index_in_vector, i)
        push!(map.bridges, nothing)
        push!(map.sets, nothing)
    end
    map.bridges[bridge_index] = call_in_context(map, bridge_index, bridge_fun)
    variables = MOI.VariableIndex[
        MOI.VariableIndex(-(bridge_index - 1 + i)) for
        i in 1:MOI.dimension(set)
    ]
    if map.unbridged_function !== nothing
        mappings = unbridged_map(map.bridges[bridge_index], variables)
        if mappings === nothing
            map.unbridged_function = nothing
        else
            for mapping in mappings
                push!(
                    map.unbridged_function,
                    mapping.first => (bridge_index, mapping.second),
                )
            end
        end
    end
    index = first(variables).value
    return variables,
    MOI.ConstraintIndex{MOI.VectorOfVariables,typeof(set)}(index)
end

"""
    function_for(map::Map, ci::MOI.ConstraintIndex{MOI.VariableIndex})

Return `vi` where `vi` is the bridged variable
corresponding to `ci`.
"""
function function_for(::Map, ci::MOI.ConstraintIndex{MOI.VariableIndex})
    return MOI.VariableIndex(ci.value)
end

"""
    function_for(map::Map, ci::MOI.ConstraintIndex{MOI.VectorOfVariables})

Return `MOI.VectorOfVariables(vis)` where `vis` is the vector of bridged
variables corresponding to `ci`.
"""
function function_for(map::Map, ci::MOI.ConstraintIndex{MOI.VectorOfVariables})
    variables = MOI.VariableIndex[]
    for i in ci.value:-1:-length(map.bridges)
        vi = MOI.VariableIndex(i)
        if map.index_in_vector[-vi.value] == -1
            continue
        elseif bridge_index(map, vi) == -ci.value
            push!(variables, vi)
        else
            break
        end
    end
    return MOI.VectorOfVariables(variables)
end

"""
    throw_if_cannot_unbridge(map::Map)

Throw an error if some bridged variables do not have any reverse mapping.
"""
function throw_if_cannot_unbridge(map::Map)
    if map.unbridged_function === nothing
        error(
            "Cannot unbridge function because some variables are bridged by",
            " variable bridges that do not support reverse mapping, e.g.,",
            " `ZerosBridge`.",
        )
    end
end

"""
    unbridged_function(map::Map, vi::MOI.VariableIndex)

Return the expression of `vi` in terms of bridged variables.
"""
function unbridged_function(map::Map, vi::MOI.VariableIndex)
    throw_if_cannot_unbridge(map)
    context_func = get(map.unbridged_function, vi, nothing)
    if context_func === nothing
        return nothing
    end
    bridge_index, func = context_func
    # If the bridge bridging `vi` has index `bridge_index` or directly or
    # indirectly created this bridge then we don't unbridge the variable.
    context = map.current_context
    while !iszero(context)
        if bridge_index == context
            return nothing
        end
        context = map.parent_index[context]
    end
    return func
end

"""
    call_in_context(map::Map, bridge_index::Int64, f::Function)

Call function `f` in the context of the variable bridge of index `bridge_index`.
That is, the variable indices bridged by this bridge or the bridges that
created it will not be unbridged in [`unbridged_function`](@ref).
"""
function call_in_context(map::Map, bridge_index::Int64, f::Function)
    if iszero(bridge_index)
        return f()
    end
    previous_context = map.current_context
    map.current_context = bridge_index
    output = nothing
    try
        output = f()
    finally
        map.current_context = previous_context
    end
    return output
end

"""
    call_in_context(map::Map, vi::MOI.VariableIndex, f::Function)

Shortcut for `call_in_context(map, bridge_index, () -> f(bridge))` where
`vi` is bridged by `bridge` with index `bridge_index`.
"""
function call_in_context(map::Map, vi::MOI.VariableIndex, f::Function)
    idx = bridge_index(map, vi)
    return call_in_context(map, idx, () -> f(map.bridges[idx]))
end

"""
    call_in_context(map::Map, ci::MOI.ConstraintIndex, f::Function)

Shortcut for `call_in_context(map, bridge_index, f)` where `bridge_index` is the
variable bridge that created `ci` (directly or indirectly) or 0 otherwise.
"""
function call_in_context(map::Map, ci::MOI.ConstraintIndex, f::Function)
    return call_in_context(map, get(map.constraint_context, ci, Int64(0)), f)
end

"""
    register_context(map::Map, ci::MOI.ConstraintIndex)

Register the current context as the variable bridge that created `ci` (directly
or indirectly) or 0 otherwise.
"""
function register_context(map::Map, ci::MOI.ConstraintIndex)
    if !iszero(map.current_context)
        # By only storing non-zero values, we avoid any dictionary access for
        # constraint not created (directly or indirectly) by variable bridges.
        # This ensures that there is no performance hit of the bridge layer when
        # no variable bridge is used.
        map.constraint_context[ci] = map.current_context
    end
    return
end

"""
    EmptyMap <: AbstractDict{MOI.VariableIndex, AbstractBridge}

Empty version of [`Map`](@ref). It is used by
[`MathOptInterface.Bridges.Constraint.SingleBridgeOptimizer`](@ref) as it does
not bridge any variable.
"""
struct EmptyMap <: AbstractDict{MOI.VariableIndex,AbstractBridge} end

Base.isempty(::EmptyMap) = true

function Base.empty!(::EmptyMap) end

Base.length(::EmptyMap) = 0

Base.keys(::EmptyMap) = MOIU.EmptyVector{MOI.VariableIndex}()

Base.values(::EmptyMap) = MOIU.EmptyVector{AbstractBridge}()

has_bridges(::EmptyMap) = false

number_of_variables(::EmptyMap) = 0

number_with_set(::EmptyMap, ::Type{<:MOI.AbstractSet}) = 0

function constraints_with_set(::EmptyMap, S::Type{<:MOI.AbstractSet})
    return MOI.ConstraintIndex{MOIU.variable_function_type(S),S}[]
end

register_context(::EmptyMap, ::MOI.ConstraintIndex) = nothing

call_in_context(::EmptyMap, ::MOI.ConstraintIndex, f::Function) = f()


text:
# https://adventofcode.com/2020/day/14

import Base: ==

example1 = "
mask = XXXXXXXXXXXXXXXXXXXXXXXXXXXXX1XXXX0X
mem[8] = 11
mem[7] = 101
mem[8] = 0
"

const MEM_SIZE = 36

BitMask = Vector{Union{Nothing,Bool}}
Binary = Vector{Union{Bool}}

Base.repr(obj::BitMask) = join([Dict(nothing => 'X', true => 1, false => 0)[x] for x in obj])

struct Program
    mask::BitMask
    mem::Vector{Tuple{Integer,Integer}}
end

(==)(x::Program, y::Program) = (x.mask == y.mask) && (x.mem == y.mem)

function parse_input(string)
    programs = Program[]
    mask, mem = [], []
    for (i, row) in enumerate(split(string, '\n', keepempty=false))
        row = strip(row)
        key, value = split(row, '=')
        key, value = strip(key), strip(value)

        if key == "mask"
            if i > 1
                push!(programs, Program(mask, mem))
            end

            mask = map(x -> x == 'X' ? nothing : parse(Bool, x), collect(value))
            @assert length(mask) == MEM_SIZE
            mem = []
        else
            address = parse(Int, match(r"mem\[(\d+)\]", key).captures[1])
            value = parse(Int, value)
            push!(mem, (address, value))
        end
    end

    push!(programs, Program(mask, mem))
    return programs
end

@assert parse_input(example1) ==
        [Program([repeat([nothing], 29); 1; repeat([nothing], 4); 0; nothing], [(8, 11), (7, 101), (8, 0)])]
@assert repr(parse_input(example1)[1].mask) == "XXXXXXXXXXXXXXXXXXXXXXXXXXXXX1XXXX0X"

test_program = "
mask = XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
mem[0] = 0
mask = XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
mem[1] = 1
mem[2] = 2
"

@assert parse_input(test_program) ==
        [Program(repeat([nothing], MEM_SIZE), [(0, 0)]), Program(repeat([nothing], MEM_SIZE), [(1, 1), (2, 2)])]

function binary_to_decimal(binary::Binary)
    decimal = 0
    for (i, x) in enumerate(reverse(binary))
        decimal += 2^(i - 1) * x
    end
    return decimal
end

@assert binary_to_decimal(Bool[0, 0, 0]) == 0
@assert binary_to_decimal(Bool[0, 0, 1]) == 1
@assert binary_to_decimal(Bool[0, 1, 0]) == 2
@assert binary_to_decimal(Bool[1, 0, 0, 0, 0, 0, 0]) == 64
@assert binary_to_decimal(Bool[1, 1, 0, 0, 1, 0, 1]) == 101

function apply_mask(bits::Binary, mask::BitMask)::Binary
    idx = map(x -> !isnothing(x), mask)
    bits[idx] = mask[idx]
    return bits
end

@assert apply_mask(Bool[0, 1, 0, 1], [nothing, nothing, nothing, true]) == Bool[0, 1, 0, 1]
@assert apply_mask(Bool[0, 0, 0, 0], [nothing, true, nothing, false]) == Bool[0, 1, 0, 0]
@assert apply_mask(Bool[0, 0, 0, 0], [true, true, true, nothing]) == Bool[1, 1, 1, 0]
@assert apply_mask(Bool[0, 1, 0, 1], [true, nothing, true, nothing]) == Bool[1, 1, 1, 1]

function tobinary(value::Integer)::Binary
    return convert(Vector{Bool}, reverse(digits(value, base=2, pad=MEM_SIZE)))
end

@assert tobinary(42) ==
        Bool[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0]

function write(programs::Vector{Program})
    memory = Dict{Integer,Integer}()
    for program in programs
        for cmd in program.mem
            bits = tobinary(cmd[2])
            bits = apply_mask(bits, program.mask)
            memory[cmd[1]] = binary_to_decimal(bits)
        end
    end
    return memory
end

@assert write(parse_input(example1)) == Dict(7 => 101, 8 => 64)

function part1(input)
    program = parse_input(input)
    memory = write(program)
    sum(values(memory))
end

@assert part1(example1) == 165

"""
Apply mask:

```
address: 000000000000000000000000000000101010  (decimal 42)
mask:    000000000000000000000000000000X1001X
result:  000000000000000000000000000000X1101X
```
"""
function apply_xmask(bits::Binary, mask::BitMask)::BitMask
    @assert length(bits) == length(mask)
    bits = convert(BitMask, bits)
    bits[mask .== true] .= true
    bits[isnothing.(mask)] .= nothing
    return bits
end

@assert apply_xmask(tobinary(42), [repeat([false], 30); nothing; true; false; false; true; nothing]) ==
        [repeat([false], 30); nothing; true; true; false; true; nothing]

Base.repr(obj::Binary) = join([Dict(true => 1, false => 0)[x] for x in obj])

function write_ver2(programs::Vector{Program}; verbose=false)
    memory = Dict{Integer,Integer}()
    for program in programs
        for cmd in program.mem
            bits = tobinary(cmd[1])

            verbose && println("address: $(repr(bits)) (decimal $(binary_to_decimal(bits)))")
            verbose && println("mask:    $(repr(program.mask))")

            bits = apply_xmask(bits, program.mask)
            verbose && println("result:  $(repr(bits))\n")

            if any(isnothing.(bits))
                idx = findall(isnothing, bits)
                for replacement in Iterators.product(repeat([0:1], length(idx))...)
                    bits[idx] .= convert(Vector{Bool}, collect(replacement))
                    address = binary_to_decimal(convert(Binary, bits))
                    memory[address] = cmd[2]
                    verbose && println("$(join(convert(Vector{Int}, bits))) (decimal $(address))")
                end
            else
                address = binary_to_decimal(bits)
                memory[address] = cmd[2]
                verbose && println("$(join(convert(Vector{Int}, bits))) (decimal $(address))")
            end
        end
        verbose && println()
    end
    return memory
end

example2 = "
mask = 000000000000000000000000000000X1001X
mem[42] = 100
mask = 00000000000000000000000000000000X0XX
mem[26] = 1
"

function part2(input)
    program = parse_input(input)
    memory = write_ver2(program)
    sum(values(memory))
end

@assert part2(example2) == 208

test = read("data/day-14.txt", String)
println("Part 1: $(result1 = part1(test))")
println("Part 2: $(result2 = part2(test))")

@assert result1 == 6386593869035
@assert result2 == 4288986482164

max_stars_repo_name: twolodzko/advent-of-code

max_issues_repo_name: twolodzko/advent-of-code

max_forks_repo_name: twolodzko/advent-of-code


text:
using InformationInequalities
using Test

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


text:
using FunctionalCollections
using FactCheck

function vec(r::UnitRange)
    v = PersistentVector{Int}()
    for i=r v=push(v, i) end
    v
end

function Base.Array(r::UnitRange)
    arr = Array(Int, r)
    for i=r arr[i] = i end
    arr
end

facts("Persistent Vectors") do

    context("range constructor") do
        @fact typeof(vec(1:1000)) --> PersistentVector{Int}
        @fact typeof(pop(vec(1:1000))) --> PersistentVector{Int}
    end

    context("length") do
        @fact length(vec(1:32)) --> 32
        @fact length(vec(1:10000)) --> 10000
        @fact length(pop(vec(1:1000))) --> 999
    end

    context("accessing elements") do
        pv = vec(1:5000)

        @fact pv[1]    --> 1
        @fact pv[32]   --> 32
        @fact pv[500]  --> 500
        @fact pv[2500] --> 2500
        @fact pv[5000] --> 5000
        @fact try pv[5001]; false catch e true end --> true

        @fact try vec(1:32)[33]; false catch e true end --> true
    end

    context("accessing last") do
        @fact peek(vec(1:1000)) --> 1000
        @fact vec(1:1000)[end]  --> 1000
    end

    context("associng") do
        @fact assoc(vec(1:1000), 500, 1)[500] --> 1
    end

    context("appending") do
        @fact append(vec(1:31), 32)        --> vec(1:32)
        @fact append(vec(1:999), 1000)     --> vec(1:1000)
        @fact append(vec(1:9999), 10000)   --> vec(1:10000)
        @fact append(vec(1:99999), 100000) --> vec(1:100000)
    end

    context("structural sharing") do
        pv = vec(1:32)
        pv2 = append(pv, 33)
        @fact is(pv2.trie[1], pv.tail) --> true
    end

    context("equality") do
        v1 = vec(1:1000)
        v2 = vec(1:1000)

        @fact is(v1.trie, v2.trie) --> false
        @fact v1 --> v2

        @fact isequal(v1, v2) --> true
    end

    context("iteration") do
        arr2 = Int[]
        for i in vec(1:10000)
            push!(arr2, i)
        end
        @fact collect(1:10000) --> arr2
    end

    context("map") do
        v1 = vec(1:5)
        @fact map((x)->x+1, v1) --> PersistentVector([2, 3, 4, 5, 6])
    end

    context("filter") do
        v1 = vec(1:5)
        @fact filter(iseven, v1) --> PersistentVector([2, 4])
    end

    context("hash") do
        @fact hash(vec(1:1000)) --> hash(vec(1:1000))
    end

    context("isempty") do
        @fact PersistentVector{Int}() --> isempty
        @fact PersistentVector([1])   --> not(isempty)
    end

end


text:
using NS2D
using Plots
using FFTW

function plot_vorticity(ω, p)

    u, v = compute_velocity(ω, p)

    ω = real(ifft(ω))

    x, y = mesh(p)

    fig = heatmap(
        x, y, ω,
        aspect_ratio=:equal,
        c=cgrad(:RdBu),
        xlims=(-p.Lx/2, p.Lx/2),
        ylims=(-p.Ly/2, p.Ly/2),
        clims=(-maximum(abs.(ω)), maximum(abs.(ω))),
        legend=false,
        size=(400, 400),
        axis=nothing,
        showaxis=false,
    )

    return fig
end

function main(; Lx=2π, Ly=2π, nx=256, ny=256, ν=1e-3, T=10)

    # --> Parameters for the simulation.
    p = SimParams(Lx=Lx, Ly=Ly, nx=nx, ny=ny, ν=ν)

    # --> Dipole initial condition.
    ω = dipole(p)

    # --> Keyword arguments for DifferentialEquations.jl
    kwargs = Dict(
        :saveat => 0:0.5:T
    )

    # --> Run the simulation.
    sol = simulate(UnforcedProblem(p, ω), T ; kwargs...)

    # --> Animate the dipole.
    x, y = mesh(p)

    anim = @animate for i = 1:length(sol.t)
        plot_vorticity(sol[i], p)
    end

    gif(anim, "../imgs/dipole.gif", fps=25)

    return
end


text:
@testset "predicates" begin
    Mag = Arblib.Mag

    @testset "predicates" begin
        for (T, fs) in (
            (Mag, (isfinite, isinf, Arblib.isspecial, iszero)),
            (
                Arf,
                (
                    isfinite,
                    isinf,
                    isinteger,
                    isnan,
                    Arblib.isneginf,
                    Arblib.isnormal,
                    isone,
                    Arblib.isposinf,
                    Arblib.isspecial,
                    iszero,
                ),
            ),
            (
                Arb,
                (
                    Arblib.isexact,
                    isfinite,
                    isinteger,
                    Arblib.isnegative,
                    Arblib.isnonnegative,
                    Arblib.isnonpositive,
                    Arblib.isnonzero,
                    isone,
                    Arblib.ispositive,
                    iszero,
                ),
            ),
            (Acb, (Arblib.isexact, isfinite, isinteger, isone, isreal, iszero)),
            (ArbPoly, (isone, Arblib.isx, iszero)),
            (ArbSeries, (isone, Arblib.isx, iszero)),
        )
            for f in fs
                @test f(zero(T)) isa Bool
            end
        end
    end

    @testset "isequal" begin
        @test isequal(Mag(0.0), Mag(0.0))
        @test !isequal(Mag(0.0), Mag(1.0))

        @test isequal(Arf(0), Arf(0))
        @test !isequal(Arf(0), Arf(1))
        @test isequal(Arf(NaN), Arf(NaN))
        @test isequal(Arf(0), 0)
        @test isequal(0, Arf(0))

        @test isequal(Arb(0), Arb(0))
        @test !isequal(Arb(0), Arb(1))
        @test isequal(Arb(NaN), Arb(NaN))
        @test isequal(Arb(0), 0)
        @test isequal(0, Arb(0))

        @test isequal(Acb(0), Acb(0))
        @test !isequal(Acb(0), Acb(1))
        @test isequal(Acb(NaN), Acb(NaN))
        @test isequal(Acb(0), 0)
        @test isequal(0, Acb(0))
    end

    @testset "comparison" begin
        x, y, z = Mag(0.0), Mag(1.0), Mag(Inf)
        @test isless(x, y)
        @test isless(y, z)
        @test x < y
        @test !(y < x)
        @test x <= x

        x, y, z = Arf(0.0), Arf(1.0), Arf(NaN)
        @test isless(x, y)
        @test isless(x, z)
        @test !isless(z, x)
        @test x < y
        @test !(x < z)
        @test !(z < x)
        @test !(x <= z)
        @test !(z <= x)
        @test x < UInt64(1)
        @test x <= UInt64(1)
        @test x < 1
        @test x <= UInt64(1)
        @test x < 1.1
        @test x <= UInt64(1)

        x, y = Arb(0), Arb(1)
        @test x == x
        @test x == 0
        @test 0 == x
        @test y == 1
        @test 1 == y
        @test x != y
        @test x < y
        @test x <= y
        @test y > x
        @test y >= x

        x, y = Acb(0), Acb(1)
        @test x == x
        @test x == 0
        @test 0 == x
        @test y == 1
        @test 1 == y
        @test x != y

        P1 = ArbPoly(Arb[1, 2])
        P2 = ArbPoly(Arb[1, 2, 0])
        P3 = ArbPoly(Arb.([1, "2 +/- 1"]))
        P4 = ArbPoly(Arb[1, 4])
        @test P1 == P2
        @test !(P1 != P2)
        @test !(P1 == P3)
        @test !(P1 != P3)
        @test !(P1 == P4)
        @test P1 != P4
        @test !(P3 == P4)
        @test P3 != P4
        @test isequal(P3, P3)
    end

end


text:
using Revise
using JDF, CSV, DataFrames
const path = "c:/data/AirOnTimeCSV"
const outpath = "c:/data/AirOnTimeCSV_jdf_pp"

convert_csv_to_jdf(infile, outdir) = begin
	print("$(infile): csv read: ")
	@time df = CSV.read(infile)
	print("$(infile): compression: ")
	@time type_compress!(df)
	print("$(infile): write_out: ")
	@time psavejdf(outdir, df)
end

for infile in readdir(path)
	@time convert_csv_to_jdf(
		joinpath(path, infile),
		joinpath(outpath, infile)
	)
	println("")
	println("")
end

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
@inline UNITLESS_ABS2(x) = Real(abs2(x)/(typeof(x)(one(x))*typeof(x)(one(x))))

@inline ODE_DEFAULT_NORM(u::Union{AbstractFloat,Complex},t) = abs(u)
@inline ODE_DEFAULT_NORM(u::Array{T},t) where T<:Union{AbstractFloat,Complex} =
                                         @fastmath sqrt(sum(abs2,u) / length(u))
@inline ODE_DEFAULT_NORM(u::RecursiveArrayTools.AbstractVectorOfArray,t) =
                                        sum(sqrt(sum(UNITLESS_ABS2,_u) / length(_u)) for _u in u.u)
@inline ODE_DEFAULT_NORM(u::AbstractArray,t) = sqrt(sum(UNITLESS_ABS2,u) / length(u))
@inline ODE_DEFAULT_NORM(u::AbstractArray{T,N},t) where {T<:AbstractArray,N} = sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip(u,t)) / length(u))
@inline ODE_DEFAULT_NORM(u,t) = norm(u)

@inline ODE_DEFAULT_ISOUTOFDOMAIN(u,p,t) = false
@inline ODE_DEFAULT_PROG_MESSAGE(dt,u,p,t) =
           "dt="*string(dt)*"\nt="*string(t)*"\nmax u="*string(maximum(abs.(u)))
@inline ODE_DEFAULT_UNSTABLE_CHECK(dt,u,p,t) = false
@inline ODE_DEFAULT_UNSTABLE_CHECK(dt,u::AbstractFloat,p,t) = isnan(u)
@inline ODE_DEFAULT_UNSTABLE_CHECK(dt,u::Float64,p,t) =
                                                any(x->(isnan(x) || x>1e50),u)
@inline ODE_DEFAULT_UNSTABLE_CHECK(dt,u::AbstractArray{T},p,t) where
                                    {T<:AbstractFloat} = any(isnan,u)
@inline ODE_DEFAULT_UNSTABLE_CHECK(dt,u::ArrayPartition,p,t) =
                                                 any(any(isnan,x) for x in u.x)


text:
struct Normalized{V <: SemiMetric} <: SemiMetric
    dist::V
    max_dist::Float64
end

function (dist::Normalized{<:Hamming})(s1, s2)
    ((s1 === missing) | (s2 === missing)) && return missing
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    len2 == 0 && return 1.0
    out = dist.dist(s1, s2) / len2
    out > dist.max_dist ? 1.0 : out
end

function (dist::Normalized{<:Union{Levenshtein{Nothing}, DamerauLevenshtein{Nothing}}})(s1, s2)
    ((s1 === missing) | (s2 === missing)) && return missing
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    len2 == 0 && return 1.0
    if dist.dist isa Levenshtein
        d = Levenshtein(ceil(Int, len2 * dist.max_dist))(s1, s2)
    else
        d = DamerauLevenshtein(ceil(Int, len2 * dist.max_dist))(s1, s2)
    end
    out = d / len2
    out > dist.max_dist ? 1.0 : out
end

function (dist::Normalized{<:AbstractQGramDistance})(s1, s2)
    ((s1 === missing) | (s2 === missing)) && return missing
    # When string length < q for qgram distance, returns s1 == s2
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    len1 <= dist.dist.q - 1 && return convert(Float64, s1 != s2)
    if dist.dist isa QGram
        out = dist.dist(s1, s2) / (len1 + len2 - 2 * dist.dist.q + 2)
    else
        out = dist.dist(s1, s2)
    end
    out > dist.max_dist ? 1.0 : out
end

function (dist::Normalized)(s1, s2)
    out = dist.dist(s1, s2)
    out > dist.max_dist ? 1.0 : out
end

"""
   normalize(dist)

Creates a normalized distance. The distance always return a Float64 between 0.0 and 1.0 (or a missing if one of the argument is missing)

### Examples
```julia-repl
julia> s1 = "New York Mets vs Atlanta"
julia> s2 = "Atlanta Braves vs New York Mets"
julia> Levenshtein()(s1, s2)
25
julia> StringDistances.normalize(Levenshtein())(s1, s2)
0.8064
```
"""
normalize(dist::SemiMetric; max_dist = 1.0) = Normalized{typeof(dist)}(dist, max_dist)
normalize(dist::Partial; max_dist = 1.0) = Partial(normalize(dist.dist; max_dist = max_dist))
normalize(dist::TokenSort; max_dist = 1.0) = TokenSort(normalize(dist.dist; max_dist = max_dist))
normalize(dist::TokenSet; max_dist = 1.0) = TokenSet(normalize(dist.dist; max_dist = max_dist))
normalize(dist::Normalized; max_dist = 1.0) = Normalized(dist.dist, max_dist)

"""
   TokenMax(dist)

Creates the `TokenMax{dist}` distance

`TokenMax{dist}` normalizes the distance `dist` and returns the minimum of the distance,
its [`Partial`](@ref) modifier, its [`TokenSort`](@ref) modifier, and its
[`TokenSet`](@ref) modifier, with penalty terms depending on string length.

It is only defined on AbstractStrings

### Examples
```julia-repl
julia> s1 = "New York Mets vs Atlanta"
julia> s2 = "Atlanta Braves vs New York Mets"
julia> evaluate(TokenMax(RatcliffObershelp()), s1, s2)
0.05
```
"""
struct TokenMax{S <: SemiMetric} <: SemiMetric
    dist::S
    max_dist::Float64
end
TokenMax(dist::SemiMetric; max_dist = 1.0) = TokenMax(dist, max_dist)
normalize(dist::TokenMax; max_dist = 1.0) = TokenMax(dist.dist, max_dist)

function (dist::TokenMax)(s1::Union{AbstractString, Missing}, s2::Union{AbstractString, Missing})
    ((s1 === missing) | (s2 === missing)) && return missing
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    max_dist = dist.max_dist
    dist0 = normalize(dist.dist; max_dist = max_dist)
    score = dist0(s1, s2)
    min_score = min(max_dist, score)
    unbase_scale = 0.95
    # if one string is much shorter than the other, use partial
    if length(s2) >= 1.5 * length(s1)
        partial_scale = length(s2) > (8 * length(s1)) ? 0.6 : 0.9
        dist0 = normalize(dist0, max_dist = 1 - (1 - max_dist) / partial_scale)
        score_partial = 1 - partial_scale * (1 - Partial(dist0)(s1, s2))
        min_score = min(max_dist, score_partial)
        dist0 = normalize(dist0, max_dist = 1 - (1 - max_dist) / (unbase_scale * partial_scale))
        score_sort = 1 - unbase_scale * partial_scale * (1 - TokenSort(Partial(dist0))(s1, s2))
        max_dist = min(max_dist, score_sort)
        dist0 = normalize(dist0, max_dist = 1 - (1 - max_dist) / (unbase_scale * partial_scale))
        score_set = 1 - unbase_scale * partial_scale * (1 - TokenSet(Partial(dist0))(s1, s2))
        out = min(score, score_partial, score_sort, score_set)
    else
        dist0 = normalize(dist0, max_dist = 1 - (1 - max_dist) / unbase_scale)
        score_sort = 1 - unbase_scale * (1 - TokenSort(dist0)(s1, s2))
        max_dist = min(max_dist, score_sort)
        dist0 = normalize(dist0, max_dist = 1 - (1 - max_dist) / unbase_scale)
        score_set = 1 - unbase_scale * (1 - TokenSet(dist0)(s1, s2))
        out = min(score, score_sort, score_set)
    end
    out > max_dist ? 1.0 : out
end

const StringDistance = Union{Hamming, Jaro, JaroWinkler,Levenshtein, DamerauLevenshtein, RatcliffObershelp, AbstractQGramDistance, Partial, TokenSort, TokenSet, TokenMax, Normalized}

"""
    compare(s1, s2, dist)

return a similarity score between 0 and 1 for the strings `s1` and
`s2` based on the distance `dist`.

### Examples
```julia-repl
julia> compare("martha", "marhta", Levenshtein())
0.6666666666666667
```
"""
function compare(s1, s2, dist::StringDistance; min_score = 0.0)
    1 - normalize(dist, max_dist = 1 - min_score)(s1, s2)
end

"""
    findnearest(s, itr, dist::StringDistance) -> (x, index)

`findnearest` returns the value and index of the element of `itr` that has the
lowest distance with `s` according to the distance `dist`.

It is particularly optimized for [`Levenshtein`](@ref) and [`DamerauLevenshtein`](@ref) distances
(as well as their modifications via [`Partial`](@ref), [`TokenSort`](@ref), [`TokenSet`](@ref), or [`TokenMax`](@ref)).

### Examples
```julia-repl
julia> using StringDistances
julia> s = "Newark"
julia> iter = ["New York", "Princeton", "San Francisco"]
julia> findnearest(s, iter, Levenshtein())
("NewYork", 1)
julia> findnearest(s, iter, Levenshtein(); min_score = 0.9)
(nothing, nothing)
```
"""
function findnearest(s, itr, dist::StringDistance; min_score = 0.0)
    min_score_atomic = Threads.Atomic{Float64}(min_score)
    scores = [0.0 for _ in 1:Threads.nthreads()]
    is = [0 for _ in 1:Threads.nthreads()]
    s = _helper(dist, s)
    # need collect since @threads requires a length method
    for i in collect(eachindex(itr))
        score = compare(s, _helper(dist, itr[i]), dist; min_score = min_score_atomic[])
        score_old = Threads.atomic_max!(min_score_atomic, score)
        if score >= score_old
            scores[Threads.threadid()] = score
            is[Threads.threadid()] = i
        end
    end
    imax = is[argmax(scores)]
    imax == 0 ? (nothing, nothing) : (itr[imax], imax)
end
_helper(dist::AbstractQGramDistance, ::Missing) = missing
_helper(dist::AbstractQGramDistance, s) = QGramSortedVector(s, dist.q)
_helper(dist::StringDistance, s) = s

function Base.findmax(s, itr, dist::StringDistance; min_score = 0.0)
    @warn "findmax(s, itr, dist; min_score) is deprecated. Use findnearest(s, itr, dist; min_score)"
    findnearest(s, itr, dist; min_score = min_score)
end
"""
    findall(s, itr , dist::StringDistance; min_score = 0.8)

`findall` returns the vector of indices for elements of `itr` that have a
similarity score higher or equal than `min_score` according to the distance `dist`.
If there are no such elements, return an empty array.

It is particularly optimized for [`Levenshtein`](@ref) and [`DamerauLevenshtein`](@ref) distances
(as well as their modifications via `Partial`, `TokenSort`, `TokenSet`, or `TokenMax`).

### Examples
```julia-repl
julia> using StringDistances
julia> s = "Newark"
julia> iter = ["Newwark", "Princeton", "San Francisco"]
julia> findall(s, iter, Levenshtein())
1-element Array{Int64,1}:
 1
julia> findall(s, iter, Levenshtein(); min_score = 0.9)
0-element Array{Int64,1}
```
"""
function Base.findall(s, itr, dist::StringDistance; min_score = 0.8)
    out = [Int[] for _ in 1:Threads.nthreads()]
    s = _helper(dist, s)
    # need collect since @threads requires a length method
    Threads.@threads for i in collect(eachindex(itr))
        score = compare(s, _helper(dist, itr[i]), dist; min_score = min_score)
        if score >= min_score
            push!(out[Threads.threadid()], i)
        end
    end
    vcat(out...)
end


text:
#===============================================================================
#
#	src/geometry.jl
#
#	This file contains:
#
#	 - findCenter(
#			P::Matrix
#		)::Array{Float64,1}
#
#	 - findClosestPoint(
#			Psimplex::Matrix,
#			P::Matrix;
#			metric = "circumcenter"
#		)::Union{Int64, Nothing}
#
#	 - findMedian(
#			P::Matrix,
#			ax::Int64
#		)::Float64
#
#	 - findRadius(
#			P::Matrix,
#			center=false;
#			digits=64
#		)::Union{Float64, Tuple{Float64, Array{Float64,1}}}
#
#	 - matrixPerturbation(
#			M::Array{Float64,2};
#			atol = 1e-10,
#			row = [0],
#			col = [0]
#		)::Array{Float64,2}
#
#	 - oppositeHalfSpacePoints(
#			P::Matrix,
#			face::Array{Float64,2},
#			point::Array{Float64,1}
#		)::Array{Int64,1}
#
#	 - planarIntersection(
#			P::Matrix,
#			face::Array{Int64,1},
#			axis::Int64,
#			off::Float64
#		)::Int64
#
#	 - simplexFaces(
#			σ::Array{Int64,1}
#		)::Array{Array{Int64,1},1}
#
#	 - vertexInCircumball(
#			P::Matrix,
#			α_char::Float64,
#			point::Array{Float64,2}
#		)::Bool
#
===============================================================================#

"""
	findCenter(P::Matrix)::Array{Float64,1}

Evaluates the circumcenter of the `P` points.

If the points lies on a `d-1` circumball then the function is not able
to perform the evaluation and therefore returns a `NaN` array.

# Examples

```jldoctest
julia> V = [
 0.0 1.0 0.0 0.0
 0.0 0.0 1.0 0.0
 0.0 0.0 0.0 1.0
];

julia> AlphaStructures.findCenter(V)
3-element Array{Float64,1}:
 0.5
 0.5
 0.5

```
"""
function findCenter(P::Matrix)::Array{Float64,1}
	dim, n = size(P)
	@assert n > 0		"findCenter: at least one points is needed."
	@assert dim >= n-1	"findCenter: Too much points"

	@assert dim < 4		"findCenter: Function not yet Programmed."

	if n == 1
		center = P[:, 1]

	elseif n == 2
		#for each dimension
		center = (P[:, 1] + P[:, 2]) / 2

	elseif n == 3
		#https://www.ics.uci.edu/~eppstein/junkyard/circumcenter.html
		if dim == 2
			denom = 2 * LinearAlgebra.det([ P[:, 2] - P[:, 1]  P[:, 3] - P[:, 1] ])
			deter = (P[:, 2] - P[:, 1]) * LinearAlgebra.norm(P[:, 3] - P[:, 1])^2 -
					(P[:, 3] - P[:, 1]) * LinearAlgebra.norm(P[:, 2] - P[:, 1])^2
			numer = [- deter[2], deter[1]]
			center = P[:, 1] + numer / denom

		elseif dim == 3
			#circumcenter of a triangle in R^3
			numer = LinearAlgebra.norm(P[:, 3] - P[:, 1])^2 * LinearAlgebra.cross(
						LinearAlgebra.cross(P[:, 2] - P[:, 1], P[:, 3] - P[:, 1]),
						P[:, 2] - P[:, 1]
					) +
					LinearAlgebra.norm(P[:, 2] - P[:, 1])^2 * LinearAlgebra.cross(
				  		P[:, 3] - P[:, 1],
						LinearAlgebra.cross(P[:, 2] - P[:, 1], P[:, 3] - P[:, 1]
					)
			)
			denom = 2 * LinearAlgebra.norm(
				LinearAlgebra.cross(P[:, 2] - P[:, 1], P[:, 3] - P[:, 1])
			)^2
			center = P[:, 1] + numer / denom
		end

	elseif n == 4 #&& dim = 3
		# https://people.sc.fsu.edu/~jburkardt/presentations
		#	/cg_lab_tetrahedrons.pdf
		# page 6 (matrix are transposed)
		α = LinearAlgebra.det([P; ones(1, 4)])
		sq = sum(abs2, P, dims = 1)
		Dx = LinearAlgebra.det([sq; P[2:2,:]; P[3:3,:]; ones(1, 4)])
		Dy = LinearAlgebra.det([P[1:1,:]; sq; P[3:3,:]; ones(1, 4)])
		Dz = LinearAlgebra.det([P[1:1,:]; P[2:2,:]; sq; ones(1, 4)])
		center = [Dx; Dy; Dz]/2α
	end

	return center
#	AlphaStructures.foundCenter([P[:,i] for i = 1 : size(P, 2)])[:,:]
end

#-------------------------------------------------------------------------------

"""
	findClosestPoint(
		Psimplex::Matrix, P::Matrix;
		metric = "circumcenter"
	)::Union{Int64, Nothing}

Returns the index of the closest point in `P` to the `Psimplex` points,
according to the distance determined by the keyword argument `metric`.
Possible choices are:
 - `circumcenter`: (default) returns the point that minimize the circumradius
 - `dd`: like `circumcenter` but the circumradius is considered to be negative
    if the circumcenter is opposite to the new point with respect to `Psimplex`.
"""
function findClosestPoint(
		Psimplex::Matrix, P::Matrix;
		metric = "circumcenter"
	)::Union{Int64, Nothing}

	@assert metric ∈ ["circumcenter", "dd"] "findClosestPoint: available metrics are
		`circumcenter` and `dd`."

	simplexDim = size(Psimplex, 2)
	@assert simplexDim <= size(Psimplex, 1) "findClosestPoint: Cannot add
	another point to the simplex."

	@assert (m = size(P, 2)) != 0 "findClosestPoint: No Points in `P`."

	radlist = SharedArray{Float64}(m)
	@sync @distributed for col = 1 : m
		r, c = findRadius([Psimplex P[:,col]], true)
		sameSign = (
			r == Inf ||
			metric != "dd" ||
			isempty(AlphaStructures.oppositeHalfSpacePoints(
				[Psimplex P[:,col]], Psimplex, c
			))
		)
		radlist[col] = ((-1)^(1 + sameSign)) * r
	end

	radius, closestidx = findmin(radlist)

	if radius == Inf
		closestidx = nothing
	end

	return closestidx

end

#-------------------------------------------------------------------------------

"""
	findMedian(P::Matrix, ax::Int64)::Float64

Returns the median of the `P` points across the `ax` axis
"""
function findMedian(P::Matrix, ax::Int64)::Float64
	xp = sort(unique(P[ax, :]))
	if length(xp) == 1
		median = xp[1]
	else
		idx = Int64(floor(length(xp)/2))
		median = (xp[idx] + xp[idx+1])/2
	end
	return median
end

#-------------------------------------------------------------------------------

"""
	findRadius(
		P::Matrix, center=false; digits=64
	)::Union{Float64, Tuple{Float64, Array{Float64,1}}}

Returns the value of the circumball radius of the given points.
If the function findCenter is not able to determine the circumcenter
than the function returns `Inf`.

If the optional argument `center` is set to `true` than the function
returns also the circumcenter cartesian coordinates.

_Obs._ Due to numerical approximation errors, the radius is choosen
as the smallest distance between a point in `P` and the center.

# Examples
```jldoctest

julia> V = [
 0.0 1.0 0.0 0.0
 0.0 0.0 1.0 0.0
 0.0 0.0 0.0 1.0
];

julia> AlphaStructures.findRadius(V)
 0.8660254037844386

julia> AlphaStructures.findRadius(V, true)
 (0.8660254037844386, [0.5, 0.5, 0.5])

```
"""
function findRadius(
		P::Matrix, center=false; digits=64
	)::Union{Float64, Tuple{Float64, Array{Float64,1}}}

 	c = AlphaStructures.findCenter(P)
	if any(isnan, c)
		r = Inf
	else
		r = round(
			findmin([LinearAlgebra.norm(c - P[:, i]) for i = 1 : size(P, 2)])[1],
			digits = digits
		)
	end
	if center
		return r, c
	end
	return r
end

#-------------------------------------------------------------------------------

"""
	matrixPerturbation(
		M::Array{Float64,2};
		atol = 1e-10,
		row = [0],
		col = [0]
	)::Array{Float64,2}

Returns the matrix `M` with a ±`atol` perturbation on each value determined
by the `row`-th rows and `col`-th columns.
If `row` / `col` are set to `[0]` (or not specified) then all the
rows / columns are perturbated.

# Examples
```jldoctest

julia> V = [
 0.0 1.0 0.0 0.0
 0.0 0.0 1.0 0.0
 0.0 0.0 0.0 1.0
];

julia> AlphaStructures.matrixPerturbation(V)
3×4 Array{Float64,2}:
 -1.27447e-11   1.0           4.68388e-11  3.08495e-11
  1.17657e-11  -4.25106e-11   1.0          5.62396e-11
  7.01741e-11  -4.10229e-11  -4.36708e-11  1.0

"""
function matrixPerturbation(
		M::Array{Float64,2};
		atol=1e-10, row = [0], col = [0]
	)::Array{Float64,2}

	if atol == 0.0
		println("Warning: no perturbation has been performed.")
		return M
	end

	if row == [0]
		row = [i for i = 1 : size(M, 1)]
	end
	if col == [0]
		col = [i for i = 1 : size(M, 2)]
	end

	N = copy(M)
	perturbation = mod.(rand(Float64, length(row), length(col)), 2*atol).-atol
	N[row, col] = M[row, col] + perturbation
	return N
end

#-------------------------------------------------------------------------------

"""
	oppositeHalfSpacePoints(
			P::Matrix,
			face::Array{Float64,2},
			point::Array{Float64,1}
		)::Array{Int64,1}

Returns the index list of the points `P` located in the halfspace defined by
`face` points that do not contains the point `point`.

_Obs._ Dimension Dipendent, only works if dimension is three or less and
	the number of points in the face is the same than the dimension.

# Examples
```jldoctest
julia> V = [
 0.0 1.0 0.0 0.0 4.0 -1. 1.0
 0.0 0.0 1.0 0.0 1.0 0.0 1.0
 0.0 0.0 0.0 1.0 2.0 0.0 1.0
];

julia> oppositeHalfSpacePoints(V, V[:, [2; 3; 4]], V[:, 1])
2-element Array{Int64,1}:
 5
 7

julia> oppositeHalfSpacePoints(V, V[:, [1; 2; 3]], V[:, 4])
0-element Array{Int64,1}

julia> oppositeHalfSpacePoints(V, V[:, [1; 3; 4]], V[:, 2])
1-element Array{Int64,1}:
 6

```
"""
function oppositeHalfSpacePoints(
		P::Matrix,
		face::Array{Float64,2},
		point::Array{Float64,1}
	)::Array{Int64,1}

	dim, n = size(P)
	noV = size(face, 2)
	@assert dim <= 3 "oppositeHalfSpacePoints: Not yet coded."
	@assert noV == dim "oppositeHalfSpacePoints:
		Cannot determine opposite to non hyperplanes."
	if dim == 1
		threshold = face[1]
		if point[1] < threshold
			opposite = [i for i = 1 : n if P[1, i] > threshold]
		else
			opposite = [i for i = 1 : n if P[1, i] < threshold]
		end
	elseif dim == 2
		if (Δx = face[1, 1] - face[1, 2]) != 0.0
			m = (face[2, 1] - face[2, 2]) / Δx
			q = face[2, 1] - m * face[1, 1]
			# false = under the line, true = over the line
			@assert point[2] ≠ m * point[1] + q "oppositeHalfSpacePoints,
				the point belongs to the face"
			side = sign(m * point[1] + q - point[2])
			opposite =
				[i for i = 1 : n if side * (m * P[1, i] + q - P[2, i]) < 0]
		else
			q = face[1, 1]
			side = sign(point[1] - q)
			opposite = [i for i = 1 : n if side * (P[1, i] - q) < 0]
		end

	elseif dim == 3
		axis = LinearAlgebra.cross(
			face[:, 2] - face[:, 1],
			face[:, 3] - face[:, 1]
		)
		off = LinearAlgebra.dot(axis, face[:, 1])
		position = LinearAlgebra.dot(point, axis)
		if position < off
			opposite = [i for i = 1:size(P, 2) if LinearAlgebra.dot(P[:,i], axis) > off]
		else
			opposite = [i for i = 1:size(P, 2) if LinearAlgebra.dot(P[:,i], axis) < off]
		end
	end

	return [
		i for i in opposite
		if sum([P[:, i] == face[:, j] for j = 1 : noV]) == 0
	]
end

#-------------------------------------------------------------------------------

"""
	planarIntersection(
		P::Matrix,
		face::Array{Int64,1},
		axis::Int64,
		off::Float64
	)::Int64

Computes the position of `face` with respect to the hyperplane `α` defined by
the normal `axis` and the contant term `off`. It returns:
 - `+0` if `f` intersect `α` internum (not only the boundary)
 - `+1` if `f` is completely contained in the positive half space of `α`
 - `-1` if `f` is completely contained in the negative half space of `α`
"""
function planarIntersection(
		P::Matrix,
		face::Array{Int64,1},
		axis::Int64,
		off::Float64
	)::Int64

	pos = [P[axis, i] > off for i in face]

	if sum([P[axis, i] == off for i in face]) == length(pos)
		position = 0 # face coplanar with axis
	elseif sum(pos) == 0
		position = -1
	elseif sum(pos) == length(pos)
		position = +1
	else
		position = 0
	end

	return position
end

#-------------------------------------------------------------------------------

"""
	simplexFaces(σ::Array{Int64,1})::Array{Array{Int64,1},1}

Returns the faces of the simplex `σ`.

_Obs._ The faces are ordered by the index value.

# Examples
```jldoctest

julia> σ = [1; 3; 2; 4];

julia> AlphaStructures.implexFaces(σ)
4-element Array{Array{Int64,1},1}:
 [1, 2, 3]
 [1, 2, 4]
 [1, 3, 4]
 [2, 3, 4]

```
"""
function simplexFaces(σ::Array{Int64,1})::Array{Array{Int64,1},1}
    sort!(sort!.(collect(Combinatorics.combinations(σ, length(σ)-1))))
end

#-------------------------------------------------------------------------------

"""
	vertexInCircumball(
		P::Matrix,
		α_char::Float64,
		point::Array{Float64,2}
	)::Bool

Determine if a point is inner of the circumball determined by `P` points
	and radius `α_char`.

"""
function vertexInCircumball(
		P::Matrix,
		α_char::Float64,
		point::Array{Float64,2}
	)::Bool

	center = AlphaStructures.findCenter(P)
	return LinearAlgebra.norm(point - center) <= α_char
end

"""
	apply_matrix(affineMatrix::Matrix, V::Matrix) -> Matrix

Apply affine transformation `affineMatrix` to points `V`.
"""
function apply_matrix(affineMatrix::Matrix, V::Matrix)
	m,n = size(V)
	W = [V; fill(1.0, (1,n))]
	T = (affineMatrix * W)[1:m,1:n]
	return T
end

function apply_matrix(affineMatrix, V::Array{Float64,1})
	T = reshape(V,3,1)
	return apply_matrix(affineMatrix, T)
end

"""
	traslation(args...) -> Matrix

Return an *affine transformation Matrix* in homogeneous coordinates.
Such `translation` Matrix has ``d+1`` rows and ``d+1`` columns,
where ``d`` is the number of translation parameters in the `args` array.

"""
function traslation(args...)
	d = length(args)
	mat = Matrix{Float64}(LinearAlgebra.I, d+1, d+1)
	for k in range(1, length=d)
			mat[k,d+1]=args[k]
	end
	return mat
end


text:
module Autohctsa
using Preferences
using CSV

function sethostname(hostname::String)
    @set_preferences!("hostname" => hostname)
    @info("New default hostname set; restart your Julia session for this change to take effect")
end
const hostname = @load_preference("hostname", nothing)
function setremotedir(remotedir::String)
    @set_preferences!("remotedir" => remotedir)
    @info("New default remotedir set; restart your Julia session for this change to take effect")
end
const remotedir = @load_preference("remotedir", nothing)

function __init__()
    if isnothing(hostname)
        @warn "No default hostname set. Use `sethostname(hostname::String)` to do so."
    end
    if isnothing(remotedir)
        @warn "No default remotedir set. Use `setremotedir(hostname::String)` to do so."
    end

    # Try connecting to remote with defaults
    shcmd = `ssh $hostname "pwd"`
    f() = run(shcmd)
    go = Task(f)
    schedule(go)
    g() = istaskdone(go)
    check = timedwait(g, 5)
    # if check != :ok
    #     println("")
    #     @warn "Could not automaticaly connect to default remote. Please check preferences, and try connecting to `hostname` manually. You might need to set up a key pair so that password authentication is not required."
    # else
    #     @info "Successfully attempted a connection to the default host. Check terminal outpu for connection status."
    # end
end

"""
This one just wraps the shell script
"""
function autohctsa(in::String, out::String, hostname::String)
    inn = abspath(in)
    @assert isfile(inn)
    sh = abspath(@__DIR__, "autohctsa.sh")
    hm = pwd()
    cd(@__DIR__)
    in = replace(inn, "\\"=>"/")
    out = replace(out, "\\"=>"/")
    @assert isfile(inn) && isfile(sh)
    @info "Submitting: $inn\n to: $out\n on: $hostname\n with: $sh\n"
    try
        #run(`dos2unix $sh`)
        shcmd = `"$sh" -i "$inn" -o "$out" -h "$hostname"`
        run(shcmd)
    catch e
        cd(hm)
        throw(e)
    end
    cd(hm)
end

"""
This one takes an input file, submits the job to the cluster and saves the result in the remote directory and saves placeholder files locally. Once the remote job is complete, the placeholder files can be retrieved using `retrieve()`.
"""
function submit(in::String, outfile::String="hctsa_"*match(r"[^_]+(?=\.)", basename(in)).match*".csv"; remotedir=remotedir, hostname=hostname, replace=false)
    isnothing(remotedir) ? (@error "Please provide an absolute path to a remote directory, or set a default") : nothing
    isnothing(hostname) ? (@error "Please provide a `<user>@<hostname>` or set a default") : nothing
    out = joinpath(remotedir, outfile)
    tmp = abspath(dirname(in), splitext(outfile)[1]*".autohctsa")
    if isfile(splitext(tmp)[1]*".csv") && !replace
        @warn "Result file already exists, skipping job"
    else
        autohctsa(in, out, hostname)
        tmpdict = Dict(:hostname=>hostname, :in=>in, :out=>out, :tmp=>tmp)
        CSV.write(tmp, tmpdict; header=false)
    end
    return nothing
end

"""
Get completed calculations from the cluster, according to any local placeholder (.autohctsa) files in the specified directory or its subdirectories
"""
function retrieve(dir::String; remotedir=remotedir, hostname=hostname)
    # We will walk the given directory recursively for files that have the extension .autohctsa
    remotefs = readlines(`ssh $hostname "find $remotedir -name \* -print"`)
    for (root, dirs, files) in walkdir(dir)
        for f ∈ files
            if splitext(f)[2] == ".autohctsa"
                D = CSV.File(abspath(root, f), header=false) |> Dict
                lookingfor = splitext(basename(f))[1]*".csv"
                if hostname == D["hostname"] && dirname(remotedir*"/") == dirname(D["out"]) && any(lookingfor .== basename.(remotefs))
                    dest = abspath(root, lookingfor)
                    remff = joinpath(remotedir, lookingfor)
                    run(`scp "$hostname:$remff" "$dest"`)
                    rm(splitext(dest)[1]*".autohctsa")
                end
            end
        end
    end
end

end

max_stars_repo_name: brendanjohnharris/Autohctsa

max_issues_repo_name: brendanjohnharris/Autohctsa

max_forks_repo_name: brendanjohnharris/Autohctsa


text:
using StaticArrays
using BenchmarkTools
using LinearAlgebra

@generated function submatrix(A::SMatrix{N1,N2,T,N1N2}, r, c) where {N1,N2,T,N1N2}
    inds = [:(($i < r ? $i : $i + 1) + (($j < c ? $j : $j + 1) - 1) * N1) for j in 1:(N2 - 1) for i in 1:(N1 - 1)]
    :(SMatrix{N1 - 1,N2 - 1,T,(N1 - 1) * (N2 - 1)}(A[SVector($(inds...))]))
end

@inline function StaticArrays._det(::Size{(3, 3)}, A::StaticMatrix{N,N,T}) where {N,T}
    cof11 =  det(submatrix(A, 1, 1))
    cof21 = -det(submatrix(A, 2, 1))
    cof31 =  det(submatrix(A, 3, 1))
    return A[1] * cof11 + A[2] * cof21 + A[3] * cof31
end

@inline function StaticArrays._det(::Size{(5, 5)}, A::StaticMatrix{N,N,T}) where {N,T}
    cof11 =  det(submatrix(A, 1, 1))
    cof21 = -det(submatrix(A, 2, 1))
    cof31 =  det(submatrix(A, 3, 1))
    cof41 = -det(submatrix(A, 4, 1))
    cof51 =  det(submatrix(A, 5, 1))
    return A[1] * cof11 + A[2] * cof21 + A[3] * cof31 + A[4] * cof41 + A[5] * cof51
end

@inline function StaticArrays._det(::Size{(6, 6)}, A::StaticMatrix{N,N,T}) where {N,T}
    cof11 =  det(submatrix(A, 1, 1))
    cof21 = -det(submatrix(A, 2, 1))
    cof31 =  det(submatrix(A, 3, 1))
    cof41 = -det(submatrix(A, 4, 1))
    cof51 =  det(submatrix(A, 5, 1))
    cof61 =  -det(submatrix(A, 6, 1))
    return A[1] * cof11 + A[2] * cof21 + A[3] * cof31 + A[4] * cof41 + A[5] * cof51 + A[6] * cof61
end

@inline function StaticArrays._inv(::Size{(3, 3)}, A)
    adj11 =  det(submatrix(A, 1, 1))
    adj12 = -det(submatrix(A, 2, 1))
    adj13 =  det(submatrix(A, 3, 1))
    adj21 = -det(submatrix(A, 1, 2))
    adj22 =  det(submatrix(A, 2, 2))
    adj23 = -det(submatrix(A, 3, 2))
    adj31 =  det(submatrix(A, 1, 3))
    adj32 = -det(submatrix(A, 2, 3))
    adj33 =  det(submatrix(A, 3, 3))

    invdet = 1 / (A[1] * adj11 + A[2] * adj12 + A[3] * adj13)

    B =  @SMatrix [
        adj11 * invdet
    	adj21 * invdet
        adj31 * invdet
        adj12 * invdet
        adj22 * invdet
        adj32 * invdet
        adj13 * invdet
        adj23 * invdet
        adj33 * invdet]
    return similar_type(A)(B)
end

@inline function StaticArrays._inv(::Size{(5, 5)}, A)
    adj11 =  det(submatrix(A, 1, 1))
    adj12 = -det(submatrix(A, 2, 1))
    adj13 =  det(submatrix(A, 3, 1))
    adj14 = -det(submatrix(A, 4, 1))
    adj15 =  det(submatrix(A, 5, 1))

    adj21 = -det(submatrix(A, 1, 2))
    adj22 =  det(submatrix(A, 2, 2))
    adj23 = -det(submatrix(A, 3, 2))
    adj24 =  det(submatrix(A, 4, 2))
    adj25 = -det(submatrix(A, 5, 2))

    adj31 =  det(submatrix(A, 1, 3))
    adj32 = -det(submatrix(A, 2, 3))
    adj33 =  det(submatrix(A, 3, 3))
    adj34 = -det(submatrix(A, 4, 3))
    adj35 =  det(submatrix(A, 5, 3))

    adj41 = -det(submatrix(A, 1, 4))
    adj42 =  det(submatrix(A, 2, 4))
    adj43 = -det(submatrix(A, 3, 4))
    adj44 =  det(submatrix(A, 4, 4))
    adj45 = -det(submatrix(A, 5, 4))

    adj51 =  det(submatrix(A, 1, 5))
    adj52 = -det(submatrix(A, 2, 5))
    adj53 =  det(submatrix(A, 3, 5))
    adj54 = -det(submatrix(A, 4, 5))
    adj55 =  det(submatrix(A, 5, 5))

    invdet = 1 / (A[1] * adj11 + A[2] * adj12 + A[3] * adj13 + A[4] * adj14 + A[5] * adj15)

    B =  @SMatrix [
        adj11 * invdet
    	adj21 * invdet
        adj31 * invdet
        adj41 * invdet
        adj51 * invdet
        adj12 * invdet
        adj22 * invdet
        adj32 * invdet
        adj42 * invdet
        adj52 * invdet
        adj13 * invdet
        adj23 * invdet
        adj33 * invdet
        adj43 * invdet
        adj53 * invdet
        adj14 * invdet
    	adj24 * invdet
        adj34 * invdet
        adj44 * invdet
        adj54 * invdet
        adj15 * invdet
    	adj25 * invdet
        adj35 * invdet
        adj45 * invdet
        adj55 * invdet]
    return similar_type(A)(B)
end

function StaticArrays._inv(::Size{(6, 6)}, A)
    adj11 =  det(submatrix(A, 1, 1))
    adj12 = -det(submatrix(A, 2, 1))
    adj13 =  det(submatrix(A, 3, 1))
    adj14 = -det(submatrix(A, 4, 1))
    adj15 =  det(submatrix(A, 5, 1))
    adj16 = -det(submatrix(A, 6, 1))

    adj21 =  det(submatrix(A, 1, 2))
    adj22 = -det(submatrix(A, 2, 2))
    adj23 =  det(submatrix(A, 3, 2))
    adj24 = -det(submatrix(A, 4, 2))
    adj25 =  det(submatrix(A, 5, 2))
    adj26 = -det(submatrix(A, 6, 2))

    adj31 =  det(submatrix(A, 1, 3))
    adj32 = -det(submatrix(A, 2, 3))
    adj33 =  det(submatrix(A, 3, 3))
    adj34 = -det(submatrix(A, 4, 3))
    adj35 =  det(submatrix(A, 5, 3))
    adj36 = -det(submatrix(A, 6, 3))

    adj41 =  det(submatrix(A, 1, 4))
    adj42 = -det(submatrix(A, 2, 4))
    adj43 =  det(submatrix(A, 3, 4))
    adj44 = -det(submatrix(A, 4, 4))
    adj45 =  det(submatrix(A, 5, 4))
    adj46 = -det(submatrix(A, 6, 4))

    adj51 =  det(submatrix(A, 1, 5))
    adj52 = -det(submatrix(A, 2, 5))
    adj53 =  det(submatrix(A, 3, 5))
    adj54 = -det(submatrix(A, 4, 5))
    adj55 =  det(submatrix(A, 5, 5))
    adj56 = -det(submatrix(A, 6, 5))

    adj61 =  det(submatrix(A, 1, 6))
    adj62 = -det(submatrix(A, 2, 6))
    adj63 =  det(submatrix(A, 3, 6))
    adj64 = -det(submatrix(A, 4, 6))
    adj65 =  det(submatrix(A, 5, 6))
    adj66 = -det(submatrix(A, 6, 6))

    invdet = 1 / (A[1] * adj11 + A[2] * adj12 + A[3] * adj13 + A[4] * adj14 + A[5] * adj15  + A[6] * adj16)

    B =  @SMatrix [
        adj11 * invdet
    	adj21 * invdet
        adj31 * invdet
        adj41 * invdet
        adj51 * invdet
        adj61 * invdet
        adj12 * invdet
        adj22 * invdet
        adj32 * invdet
        adj42 * invdet
        adj52 * invdet
        adj62 * invdet
        adj13 * invdet
        adj23 * invdet
        adj33 * invdet
        adj43 * invdet
        adj53 * invdet
        adj63 * invdet
        adj14 * invdet
    	adj24 * invdet
        adj34 * invdet
        adj44 * invdet
        adj54 * invdet
        adj64 * invdet
        adj15 * invdet
    	adj25 * invdet
        adj35 * invdet
        adj45 * invdet
        adj55 * invdet
        adj65 * invdet
        adj16 * invdet
    	adj26 * invdet
        adj36 * invdet
        adj46 * invdet
        adj56 * invdet
        adj66 * invdet]
    return similar_type(A)(B)
end


text:
using DEBTrait, Test
using CSV, DataFrames, Statistics
using Plots

df_isolates      = CSV.read("/Users/glmarschmann/.julia/dev/DEBTrait/data/microbes2traits.csv", DataFrame, missingstring = "N/A")
df_isolates      = filter(row -> row[:Rhizosphere_response] == "negative" || row[:Rhizosphere_response] == "positive" , df_isolates)

L_DNA            = convert(Array{Float64,1}, df_isolates.Genome_size)
rRNA             = DEBTrait.genome_size_to_rRNA_copy_number(L_DNA)       # need to cap this at 20
rRNA_iso         = df_isolates.rRNA_genes
min_gentime      = df_isolates.Min_gen_time
Isolate          = BaseGenome(L_DNA, rRNA, min_gentime, 1.0*ones(1,27), 1.0*ones(27))
Cell_Composition = PCellComposition{BaseGenome}(Isolate)

rRNA_protein_ratio = Cell_Composition.V_r./Cell_Composition.V_p
y_EV             = DEBTrait.reserve_yield_constant(Cell_Composition)


text:
module Bound

using ..Shapes

struct Bounds{TS,TB}
    shape::TS
    bcs::TB
end # struct

function unit_normal!(Nx,Ny,C::AbstractCircle,x,y)
    for i ∈ eachindex(Nx)
        Nx[i], Ny[i] = unit_normal(C,x[i],y[i])
    end
    return nothing
end
function unit_normal(C::AbstractCircle,x,y)
    r = hypot(x-C.x0,y-C.y0)
    return x/r,y/r
end

end # module


text:
using JSON

@everywhere include("utils.jl")

function generatemotifs(V, E)

   #=
   The first step is to generate a binary tree of depth V^2
   The from_binarytree function does that, and only returns paths with
      exactly E interactions. This is the entire set of unique permutations
      of E in V^2 elements.
   =#
   t_hashes = pmap(hashmat, from_binarytree(V^2, E))

   #=
   TODO maybe this next line will have to be removed
   =#
   t_hashes = filter((x) -> x != nothing, t_hashes)

   #=
   Finally, we return an Dict where every key is a hash (the hash is
      the unfolded matrix), and the value is the matrix itself, at the right
      dimension.
   =#
   hashes = Dict{Any,Any}([x => unhashmat(x) for x in t_hashes])
   return hashes
end

function identifymotifs(motif_set)
   hash_set = Dict()

   #=
   How many motifs do we know at the beginning? 0. Well, 1, because the first
      motif we'll come across will necessarilly be a new one.
   =#
   known_motifs = 1

   #=
   We will then move through all hashes, i.e. all unique matrix conformations.
   =#
   for k in collect(keys(motif_set))

      #=
      The first step is to collect the matrix to be evaluated in A
      =#
      A = motif_set[k]

      #=
      Then, we simply look at all the possible permutations of rows orders.
         The columns HAVE to be re-ordered in the same way since A is the
         adjacency matrix that represents the motif.
      =#
      orders = collect(permutations(vec(1:size(A, 1))))

      #=
      This empty array will store the representations of every possible order.
      =#
      list = ASCIIString[]

      #=
      For every order, we re-order the matrix, hash it, and push it to the list
         of all possible unique conformations. This is an OK overhead because we
         are often limited to motifs with V <= 5, so that even V! permutations
         are OK.
      =#
      for o1 in orders
         push!(list, mhash(A[o1, o1]))
      end

      #=
      We know want to know whether we already know this motif. By default, we
         assume that no, we don't ...
      =#
      is_known = false

      #=
         ... and it has therefore no known id. But this will change.
      =#
      has_id = 0

      #=
      Then for every possible conformation in the list (the list
         is actually the unique hashes) ...
      =#
      for l in list

         #=
         It is fairly simple. If the motif is know, i.e. if it is in the list
            of all hashes,
         =#
         if l in keys(hash_set)

            #=
            We mark it as known, and give it an ID.
            =#
            is_known = true
            has_id = hash_set[l]
            # TODO we could break of the loop at this point
         end
      end

      #=
      Then, IF we know this motif...
      =#
      if is_known

         #=
         ... we put it's ID in the hash_set ...
         =#
         for l in list
            hash_set[l] = has_id
         end
      else

         #=
         ... but if not, we create a new motif
         =#
         for l in list
            hash_set[l] = known_motifs
         end

         #=
         We implement in the end, because, remember, we started at 1 known motif.
         =#
         known_motifs += 1
      end
   end
   #=
   Done! Good job.
   =#
   return hash_set
end

hashes = Dict()

for V in 2:4
  hashes[V] = Dict()
  for E in (V-1):V^2
    println("V: ", V, "\tE: ", E)
    motifs = generatemotifs(V, E)
    hashes[V][E] = identifymotifs(motifs)
  end
end

motifs = open("motifs.json", "w")
JSON.print(motifs, hashes)
close(motifs)


text:
# Code derived from metabolicEP (https://github.com/anna-pa-m/Metabolic-EP)

function idxlicols(X; tol::Float64=1e-10)

    sum(abs2,X) == 0 && (return(Array{Int,1}(), Array{Int,2}() ))
    Q,R,E = qr(X, ColumnNorm())
    diagr = abs.(diag(R))
    r = findall(diagr .>= tol*diagr[1])[end]
    idx = sort(E[1:r])
    return idx
end

function echelonize(X::T,v; eps::Real=1e-10) where T <:DenseArray
    M,N = size(X)

    idxrow = idxlicols(Matrix(X'))
    Mred = length(idxrow)

    idxind = idxlicols(X)
    idxdep = setdiff(1:N,idxind)
    newidx = vcat(idxind,idxdep)
    Tv = @view X[idxrow,newidx]
    iTv = inv(Tv[1:Mred, 1:Mred])
    res = iTv * Tv
    # trimming zeros
    for i in eachindex(res)
        abs(res[i]) < eps && (res[i] = zero(res[i]))
    end
    bnew = iTv * v[idxrow]
    # trimming ones
    for i in 1:Mred
        abs(1.0 - res[i,i]) < eps && (res[i,i] = one(res[i,i]))
    end

    return idxdep, newidx, res, bnew

end


text:
using Grid
using Test

# generate 1D grid
g1d=Grid.M1D(0.0,2.0,1000)
g2d=Grid.M1D(0.0,2.0,1000)

@test isequal(g1d.x[1],0.0)
@test isequal(g1d.x[end],2.0)

@test isequal(g1d, g2d)