text:
#=
∈
Low Rank Updates of A \in R^nxn

Sherman-Morrison formula

If x, y ∈ R^n & 1+y^T A^-1 x is nonzero then

TODO: check for  (A+x*y') != 0 , A != 0
=#

using LinearAlgebra
# as an expression
lhs = :((A + x * y')^-1)

rhs = :((A^-1) - [(A^-1 * (x * y') * A^-1) / (1 + y' * A^-1 * x)])
#Or Generic function

A = rand(10, 10)
x = rand(10, 10)
y = rand(10, 10)
I = ones(10, 10)
Lhs(A = A, x = x, y = y) = ((A + x * y')^-1)

Rhs(A = A, x = x, y = y) = (A^-1) - [(A^-1 * (x * y') * A^-1) / (I + y' * A^-1 * x')]

Rhs(A, x, y) ≈ Lhs(A, x, y)

y' * A^1 * x' ≈ y' * A^1 * x # beware while using transpose ' is false # false positive

numerator = (A^-1 * (x * y') * A^-1)
denominator = I + y' * A^1 * x'

numeratorSize = size(numerator)
denominatorSize = size(denominator)
Div = numerator / denominator # ok

Mul = x * y'
(A + Mul)^-1

(A + Mul)^-1 ≈ A^-1 - Div # final answer  #false (negative)
denominatorSize = size(I + y' * A^1 * x') #correct
numeratorSize = size((A^-1 * (x * y') * A^-1)) # correct

#--- Sherman-Morrison-Woodbury formula
#= Sherman-Morrison-Woodbury formula

if U,V ∈ R^nxp & Ip + V' A-1*U is non singular
=#
#=  warning U is already Used in another context in another file in different 'context'  <--- TODO a need for a context (logic ) =#

u = rand(10, 10)
v = rand(10, 10)
A + u * v
#=10×10 Matrix{Float64}:
 2.71503  3.87771  2.80252  …  3.10273  3.43842
 2.69587  2.54099  1.62384     1.73684  2.25296
 ⋮                          ⋱
 3.05006  3.93606  2.79434     2.61207  4.02894
 3.50059  3.93851  2.27809     2.68894  2.88773
 =#
A + u * v'
#=
10×10 Matrix{Float64}:
 2.12618  3.19865  2.75244  3.72066  …  2.21105  3.81644  2.81843  3.19734
 1.96465  1.76703  1.64296  3.0106      1.73512  2.2817   1.92534  2.07472
 ⋮                                   ⋱
 2.69908  3.21591  3.02653  4.24413     2.68673  4.60944  2.40217  3.60195
 3.35323  3.15063  2.7179   3.24876     1.50348  3.66244  2.49389  2.75832
 =#

left = v' * A^-1
#=10×10 Matrix{Float64}:
1.4968     7.37556  0.833061  …  0.584169  -1.06658   -2.45121
0.0867906  8.30224  0.652222     0.364652  -0.477872  -4.7667
⋮                             ⋱
1.62219    4.33597  0.576026     0.916961  -0.639694  -1.27006
2.14516    2.79338  0.335954     0.421141  -1.84897   -0.723347
=#

tmpmiddle = u * (I + (v' * A' * u))^-1
#=
10×10 Matrix{Float64}:
 -1.86419  -0.695083  -1.7452     …   2.54644   2.64085    1.49266
 -3.1713   -1.84768    1.1834         2.22572   5.69861    1.40546
  ⋮                               ⋱
  2.29853   6.16219    0.0294727     -5.62846   0.167164  -2.55069
  4.09594  -0.9933    -1.79248        1.63777  -5.71097   -1.05803
 =#
middle = -A^-1 * tmpmiddle
#=
  49.4646  57.4864    -9.72253  …  -53.8581   -68.3594  -32.7693
  27.8203  24.1974    -1.55644     -26.5274   -34.1993  -16.7565
   ⋮                            ⋱
 -13.4486   4.15526   10.7227       -8.40686   18.9022    2.26729
  40.2082  18.4852   -19.0904      -13.7932   -64.1186  -18.4974
  =#
right = A^-1

RHS = right * middle * left
#=
10×10 Matrix{Float64}:
 2.12618    8.08584   0.655802  …  0.719427  -4.17327   -2.74713
 1.52245    2.78087  -0.116915     0.85684   -2.58951   -1.93996
 ⋮                              ⋱
 0.756402  -1.35546   0.561576     0.32922    0.521806   1.98675
-0.663023   6.10921   1.78613      0.112573  -1.47111   -3.54095
  6.74992  2097.73  213.576   -772.089     13.7489
-513.715  -1168.69
  ⋮                                     ⋱            1762.06   -3871.44
211.53     2470.34  202.148   -776.942     69.1312   -513.715  -1168.69
-802.528  -1352.28
 56.3864   4963.32  524.213  -1765.25      58.9692  --802.528  -1352.28     1384.89   -2718.6
 =#
Lehs = (A + u * v')^-1
#=
 -1.64577     1.22028   …  -0.751388  -0.825928
 -0.457401    1.00314      -0.378621  -1.72199
  ⋮                     ⋱
  0.357082   -0.770006      1.3016     1.68018
 -0.0114964  -1.66291      -1.31537    1.18225
 =#

#---recalled lecture 1 | convex Optimization I (stanford)

#=
47:19
actually i know ever  know how to solve A
everyone here has done A/ B -probably only a few know actually what it did ( what it means)

47:45 https://youtu.be/McLq1hEq3UY?t=2865
" Believe it or not, they won't even know what I'm asking !"
     - Professor Stephen Boyd, of the Stanford Electrical Engineering Department
    Lecture 1 | Convex Optimization 1 (Standard)

  =#

Lehs - RHS
""" the Best we can get
10×10 Matrix{Float64}:
 22.2142     25.6118    -89.3328   -7.21748   -17.99     15.6037   12.0876   77.5545     3.49734  -43.1083
  0.472091  -35.6289   -125.656    17.9696     82.6787  -81.5007   53.3611  152.23       2.13005  -57.7565
  ⋮                                                       ⋮
  2.53801    -9.63932   -31.15      3.34303    18.2504  -20.8493   13.5364   38.1611     4.72605  -15.887
 51.7025    133.216     -80.0547  -47.635    -183.818   208.653   -45.3165   -2.99068  -11.0007   -49.3098
"""

function _span(arr = A, ; diagon = LinearAlgebra.diag(A))
  span = 0
  for i in enumerate(size(diagons))
    span += diagons[i]
  end
  return span
end

diagons = LinearAlgebra.diag(A) = LinearAlgebra.Diagonal(A)
span(A, diagons)

LinearAlgebra.span

#--- identity
A = rand(10, 10)
B = rand(10, 10)
I = ones(10, 10)
C = rand(10, 10)
_0 = zeros(10, 10)
_lhs = (I + A .* B)^-1
#=
 2.51162     1.06328     3.50709   -3.73767   -1.19645   -2.50288  -0.293292   0.961015   2.15519   -2.10134
 -6.03368    -1.3807     -8.14618    9.10683    4.05341    6.37086  -1.45176   -0.658709  -5.30112    2.65344
  ⋮                                                        ⋮
  4.9177      1.4219      6.38853   -7.89654   -1.68593   -5.32977   0.260283   0.159931   4.43551   -2.16398
 -0.0382185   0.0314968   0.316589  -0.994095  -0.307413   1.4129    0.126316  -0.345001   0.232051  -0.182067
=#
#_0 ≈ _lhs

(I + B * A)^-1 * B

#=
10×10 Matrix{Float64}:
  1.92712   -1.52044   …   2.88647    0.899755
  4.59499   -5.75992       6.18216    0.0528266
  ⋮                    ⋱
 -5.99075    6.25609      -9.001     -1.07983
 -0.834267   0.204393     -0.139767  -0.228184

=#

_rhs = I - A * (I + B * A)^-1 * B

_rhs ≈ _lhs  #false(bnot close (enough))

#=Key equations

A*B*A ≈ (B*A)
=#
A * B * A ≈ A * (B * A)  # true #finally
#=fundanetbal theory of matrix theory =#

(A * B) * C ≈ A * (B * C) # true

#=computing a product
if x, y,z ∈ R^n how should we compute
  x*y'*z
takeaway: put brackets starting from the Right, there's a massive difference in computational cost
1. (x*y') * z costs O(n^2) operations
2. x(y'*z) costs O(n) operations

#TODO: Search for references

#--- matrix chain(rule) multiplication
#=where to put the parentheses in multiplication
given k matricies

A1 * A2*...An [Warning: not squares, possibly different dimensions]
me: in Julia  matrix product * is a left operation
secondly, unsure the type of these matricies
solution solve by dynamic programming
in O(k^3) operations

paper: O(K Log K  Hu, Shing (1984)
why importanT ?
O(k log(k)) algorithm by Hu & shing (1984)

where should i put parentheses to geth
=#

max_stars_repo_path: Tutorial/Nick Higham/ShermanMorrison.jl

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_path: Tutorial/Nick Higham/ShermanMorrison.jl

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_path: Tutorial/Nick Higham/ShermanMorrison.jl

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
using Base: @propagate_inbounds

export row, col
export add!, zero!, apply!, apply
export setat!, set!, set_index!, rand!

export state_similar, state_collect, state_copy, state_copy!
export state_size, batch_size, chain_length, state_length
export state_eltype, batch_eltype, state_uview, state_view
export state_reshape, state_vec

const AState{T}         = AbstractVector{T}  where T
const AStateBatch{T}    = AbstractMatrix{T}  where T
const AStateBatchVec{T} = AbstractArray{T,3} where T

const AbstractDoubled{T,N}    = NTuple{2,AbstractArray{T,N}} where {T,N}
const ADoubleState{T}         = AbstractDoubled{T,1}         where T
const ADoubleStateBatch{T}    = AbstractDoubled{T,2}         where T
const ADoubleStateBatchVec{T} = AbstractDoubled{T,3}         where T

const AbstractState{T}         = Union{AState{T}, ADoubleState{T}}                 where T
const AbstractStateBatch{T}    = Union{AStateBatch{T}, ADoubleStateBatch{T}}       where T
const AbstractStateBatchVec{T} = Union{AStateBatchVec{T}, ADoubleStateBatchVec{T}} where T

const AStateOrBatch{T}       = Union{AState{T},AStateBatch{T}} where T
const ADoubleStateOrBatch{T} = Union{ADoubleState{T},ADoubleStateBatch{T}} where T

const AStateOrBatchOrVec{T} = Union{AState{T},AStateBatch{T},AStateBatchVec{T}} where T
const ADoubleStateOrBatchOrVec{T} = Union{ADoubleState{T},ADoubleStateBatch{T},ADoubleStateBatchVec{T}} where T

const AStateBatchOrVec{T}       = Union{AStateBatch{T},AStateBatchVec{T}} where T
const ADoubleStateBatchOrVec{T} = Union{ADoubleStateBatch{T},ADoubleStateBatchVec{T}} where T

"""
    row(v::ADoubleState(orBatch,orVec)) = first(v)

Returns the row indices of the state `v`.
A doubled state is implemented as a tuple, so this effectively returns the
first (subject to change) element of the tuple.
"""
@inline row(v::ADoubleStateOrBatchOrVec) = @inbounds first(v)

"""
    col(v::ADoubleState(orBatch,orVec)) = last(v)

Returns the column indices of the state `v`.
A doubled state is implemented as a tuple, so this effectively returns the
last (subject to change) element of the tuple.
"""
@inline col(v::ADoubleStateOrBatchOrVec) = @inbounds last(v)

## utilities similar to base size

"""
    batch_size(σ::BatchedStateOrVec) -> Int

Returns the batch size of an arbitrary state.
This is the second dimension of the array (but works also on doubled states).
"""
@inline batch_size(σ::Union{AStateBatch,AStateBatchVec}) = @inbounds size(σ, 2)
@inline batch_size(σ::Union{ADoubleStateBatch,ADoubleStateBatchVec}) = batch_size(row(σ))

"""
    chain_length(σ::BatchStateVec)

Returns the length of a chain of states.
This is the third dimension of the array (but works also on doubled states).
"""
@inline chain_length(σ::AStateBatchVec) = @inbounds size(σ, 3)
@inline chain_length(σ::ADoubleStateBatchVec) = chain_length(row(σ))

"""
    state_size(σ) -> tuple

Works as `size(σ)` returning the size of an array of states, with the two
following differences:
    - It drops the first dimension (as it is the dimension of the hilbert
    space and is not relevant);
    - It also works on tuples (Doubled States).
"""
@inline state_size(σ::AbstractState) = tuple()
@inline state_size(σ::AbstractStateBatch) = (batch_size(σ),)
@inline state_size(σ::AbstractStateBatchVec) = (batch_size(σ),chain_length(σ))

@inline state_length(σ::AbstractState) = 1
@inline state_length(σ) = prod(state_size(σ))

@inline state_eltype(σ::AbstractState) = typeof(σ)
@inline state_eltype(σ::AStateBatch) = @inbounds typeof(σ[:,1])
@inline state_eltype(σ::AStateBatchVec) = @inbounds typeof(σ[:,1,1])

@inline batch_eltype(σ::AStateBatch) = typeof(σ)
@inline batch_eltype(σ::AStateBatchVec) = @inbounds typeof(σ[:,:,1])

@inline state_eltype(σ::AbstractDoubled) = Tuple{state_eltype(row(σ)), state_eltype(col(σ))}
@inline batch_eltype(σ::AbstractDoubled) = Tuple{batch_eltype(row(σ)), batch_eltype(col(σ))}

# statesimilar
"""
    state_similar(σ, [dims...])

Works as `similar(σ, dims...)` but
    - does not change the first dimension;
    - Works with tuples (Doubled states).
"""
state_similar(σ::AbstractArray, dims::Int...) = @inbounds similar(σ, size(σ,1), dims...)
state_similar(σ::AbstractArray, dims::Dims)   = @inbounds similar(σ, size(σ,1), dims...)
state_similar(σ::AbstractArray)               = similar(σ)

# Special case for doubled states
state_similar(σ::AbstractDoubled, dims::Vararg{T,N}) where {T,N} =
    (state_similar(row(σ), dims...), state_similar(col(σ), dims...))

# Allocating
state_copy(σ) = state_copy!(state_similar(σ), σ)
state_copy(σ, σp, mask) = state_copy!(state_copy(σ), σp, mask)
state_copy_invertmask(σ, σp, mask) = state_copy_invertmask(state_copy(σ), σp, mask)

"""
    state_copy!(σp, σ, [mask=nothing])

Copies the state `σ` onto `σp`.
Equivalent to `σp .= σ` in most cases, but if the state is a tuple (double state)
maps the operation on every element of the tuple.

If the BitMask `mask` is passed, only the elements where mask has 1s are copied
over, while the others are left unchanged.
"""
@inline state_copy!(σp::S, σ::Sp) where {T,N,S<:AbstractArray{T,N}, Sp<:AbstractArray{T,N}} =
    copyto!(σp, σ)
@inline state_copy!(σp::AbstractDoubled{T}, σ::AbstractArray{T}) where T =
    state_copy!(σp, (σ,σ))
@inline function state_copy!(σp::S, σ::S2) where {T,N,N2,S<:AbstractDoubled{T,N},S2<:AbstractDoubled{T,N2}}
    state_copy!(row(σp), row(σ))
    state_copy!(col(σp), col(σ))
    return σp
end

@inline state_copy!(σp::AStateBatch, σ::AState) = σp .= σ
@inline state_copy!(σp::AStateBatchVec, σ::AStateBatch) = σp .= σ
@inline function state_copy!(σp::Union{AState, AStateBatch, AStateBatchVec}, σ::Union{AState, AStateBatch, AStateBatchVec}, mask)
    σp .= σ .* mask .+ σp .* .! mask
    return σp
end
function state_copy!(σp::ADoubleStateOrBatchOrVec,
           σ::ADoubleStateOrBatchOrVec, mask)
    state_copy!(row(σp), row(σ), mask)
    state_copy!(col(σp), col(σ), mask)
    return σp
end

"""
    state_copy_invertmask!(σp, σ, mask)

Copies the state `σ` onto `σp`, but only the elements where `mask` has 0s.
Equivalent to `state_copy!(σp, σ, !mask)`.
"""
@inline function state_copy_invertmask!(σp::Union{AState, AStateBatch, AStateBatchVec}, σ::Union{AState, AStateBatch},
                   mask)
    σp .= σ .* .! mask .+ σp .* mask
    return σp
end

function state_copy_invertmask!(σp::Union{ADoubleState, ADoubleStateBatch, ADoubleStateBatchVec},
           σ::Union{ADoubleState, ADoubleStateBatch, ADoubleStateBatchVec}, mask)
    state_copy_invertmask!(row(σp), row(σ), mask)
    state_copy_invertmask!(col(σp), col(σ), mask)
    return σp
end

# TODO add @inbounds
"""
    state_uview(state, [batch], el)

Given a vector of batches of states with size `size(state) = [:, batches, els]`,
take the batch group `el`, and if specified also selects one single batch.

Returns an UnsafeArrays.UnsafeView object to avoid allocating.
"""
@propagate_inbounds state_uview(σ::AStateBatch, i)    = uview(σ, :, i)
@propagate_inbounds state_uview(σ::AStateBatchVec, i) = uview(σ, :, :, i)
@propagate_inbounds state_uview(σ::AStateBatchVec, batch, el) = uview(σ, :, batch, el)

@propagate_inbounds state_uview(σ::AbstractDoubled, i::Vararg{T,N}) where {T,N} =
    (state_uview(row(σ), i...), state_uview(col(σ), i...))
@propagate_inbounds state_uview(σ::AbstractDoubled, j::AbstractRange, i::Vararg{T,N}) where {T,N} =
    (state_uview(row(σ), j, i...), state_uview(col(σ), j, i...))

"""
    state_view(state, [batch], el)

Given a vector of batches of states with size `size(state) = [:, batches, els]`,
take the batch group `el`, and if specified also selects one single batch.
"""
@propagate_inbounds state_view(σ::AStateBatch, i)    = view(σ, :, i)
@propagate_inbounds state_view(σ::AStateBatchVec, i) = view(σ, :, :, i)
@propagate_inbounds state_view(σ::AStateBatchVec, batch, el) = view(σ, :, batch, el)

@propagate_inbounds state_view(σ::AbstractDoubled, i::Vararg{T,N}) where {T,N} =
    (state_view(row(σ), i...), state_view(col(σ), i...))
@propagate_inbounds state_view(σ::AbstractDoubled, j::AbstractRange, i::Vararg{T,N}) where {T,N} =
    (state_view(row(σ), j, i...), state_view(col(σ), j, i...))

state_reshape(σ::AbstractStateBatch, i::Colon) = σ
state_reshape(σ::AStateOrBatchOrVec, i::Union{Int,Colon}) = reshape(σ, size(σ, 1), i)
state_reshape(σ::AStateOrBatchOrVec, i::Int...) = reshape(σ, size(σ, 1), i...)
state_reshape(σ::AStateOrBatchOrVec, i::Union{Colon,Int}...) = reshape(σ, size(σ, 1), i...)
@inline state_reshape(σ::AbstractDoubled, i...) = tuple(state_reshape(row(σ), i...),
                                                        state_reshape(col(σ), i...))

state_vec(σ::AbstractStateBatchVec) = state_reshape(σ, state_length(σ))
state_vec(σ::AbstractStateBatch) = σ

ext: jl


text:
struct ArnoldiMethod_IRAM{Tv} <: AbstractEigenMethod{Tv}
end

(::Type{<:ArnoldiMethod_IRAM})(h::AbstractMatrix{Tv}) where {Tv} = ArnoldiMethod_IRAM{Tv}()

function (d::Diagonalizer{<:ArnoldiMethod_IRAM,Tv})(nev::Integer; precond = false, kw...) where {Tv}
    precond && @warn "ArnoldiMethod_IRAM does not yet support preconditioning"
    if isfinite(d.point)
        which = ArnoldiMethod.LM()
    else
        which = d.point > 0 ? ArnoldiMethod.LR() : ArnoldiMethod.SR()
    end
    decomp, _ = ArnoldiMethod.partialschur(d.lmap; nev = nev, which = which, kw...)
    λs = real.(decomp.eigenvalues)
    ϕs = decomp.Q
    isfinite(d.point) && (λs .= 1 ./ λs .+ d.point)
    return Eigen(λs, ϕs)
end


text:
###########################
# Function Cell Interface #
###########################
"""
	pca([;kwargs])

Constructs an untrained cell that when piped data inside, calculates the PCA
projection matrix of the input data. We assume that `d` is the dimesionality
of the input dataset, `n` the number of samples

# Keyword arguments (same as in `MultivariateStats`)
  * `method` is the PCA method used `:auto`, `:cov`, `:svd` (default `:auto` e.g. use `:cov` when the number of dimensions
  smaller than the number of samples, `svd` otherwise)
  * `maxoutdim` is the number of output dimensions, default `<number of variables>`
  * `pratio` is the ratio of variaces preserved in the principal subspace (default `0.99`)
  * `mean` is the mean vector, can be `nothing`, `0` or precomputed mean vector (default `nothing`)

Read the `MultivariateStats.jl` documentation for more information.
"""
pca(;kwargs...) = FunctionCell(pca, (), ModelProperties(), kwtitle("PCA", kwargs);kwargs...)

############################
# DataCell/Array Interface #
############################
"""
	pca(x, [;kwargs])

Trains a function cell that when piped data into will compute the principal components of
the observations based on the projection matrix extracted from `x`.
"""
# Training
pca(x::T where T<:CellData; kwargs...) = pca(getx!(x);kwargs...)
pca(x::T where T<:AbstractVector; kwargs...) = pca(mat(x, LearnBase.ObsDim.Constant{2}());kwargs...)
pca(x::T where T<:AbstractMatrix; kwargs...)  = begin

	ppadata = fit(MultivariateStats.PCA, getobs(x); kwargs...)

	# Build model properties
	modelprops = ModelProperties(size(ppadata.proj,1), size(ppadata.proj,2))

	# Returned trained cell
	FunctionCell(pca, Model(ppadata, modelprops), kwtitle("PCA",kwargs))
end

# Execution
pca(x::T where T<:CellData, model::Model{<:MultivariateStats.PCA}) =
	datacell(pca(getx!(x), model), gety(x))
pca(x::T where T<:AbstractVector, model::Model{<:MultivariateStats.PCA}) =
	pca(mat(x, LearnBase.ObsDim.Constant{2}()), model)
pca(x::T where T<:AbstractMatrix, model::Model{<:MultivariateStats.PCA}) =
	MultivariateStats.transform(model.data, getobs(x))

"""
	pcar(M::Cell), where `M=pca(x)`

Trains a cell based on a previously trained PCA function cell. When piped data into,
the returned output will be the reconstruction of the observations in the original
space. Obviously, the input data (e.g. data to be reconstructed) is expected to be
principal components extracted with `M`.

Read the `MultivariateStats.jl` documentation for more information.

# Examples
```
julia> P= rand(5,1000) |> pca(maxoutdim=4)
PCA: maxoutdim=4, 5 -> 4, trained

julia> PR = pcar(P)
PCA: reconstruct, 4 -> 5, trained

julia> ([1.0 2 3 4 5]' |> P ) |> pcar(P)
5×1 Array{Float64,2}:
 0.305402
 1.16768
 2.25037
 5.23595
 1.77198
```
"""
# Training for the reconstruction tranform (fixed cell)
pcar(x::T where T<:CellFunT{<:Model{<:MultivariateStats.PCA}}) = begin

	# Build model properties
	modelprops = ModelProperties(getx(x).properties.odim, getx(x).properties.idim)

	FunctionCell(pcar, Model(getx(x).data, modelprops),"PCA: reconstruct")
end

# Execution
pcar(x::T where T<:CellData, model::Model{<:MultivariateStats.PCA}) = datacell(pcar(getx!(x), model), gety(x))
pcar(x::T where T<:AbstractVector, model::Model{<:MultivariateStats.PCA}) = pcar(mat(x, LearnBase.ObsDim.Constant{2}()), model)
pcar(x::T where T<:AbstractMatrix, model::Model{<:MultivariateStats.PCA})::Matrix{Float64} = MultivariateStats.reconstruct(model.data, getobs(x))


text:
import Pkg;

@show depends = ([
                  "ArgParse";
                  "Distributions";
                  "Logging";
                  "ProfileView";
                  "DecFP";
                  "Quadmath";
                  "NLsolve";
                  "ForwardDiff";
                  "Plots";
                  "Optim";
                  "NLopt";
                  "Cubature";
                 ]);

foreach(Pkg.add, depends);

Pkg.update();

max_stars_repo_name: grasingerm/polymer-stats

max_issues_repo_name: grasingerm/polymer-stats

max_forks_repo_name: grasingerm/polymer-stats


text:
using Statistics
using LinearAlgebra
using TSVD: tsvd

"""
    PCA(;ncomps=nothing, pratio=1) <: LinearAlgorithm
    PCA(ncomps; pratio=1)

Use principal components analysis (PCA) to form a low-rank orthonormal basis of the input. Uses deterministic singular-value decomposition (SVD) to decompose data.

If `ncomps` is `nothing`, it will be set to the number of frames in the reference cube when processed.

# References
* [Soummer, Pueyo, and Larkin (2012)](https://ui.adsabs.harvard.edu/abs/2012ApJ...755L..28S) "Detection and Characterization of Exoplanets and Disks Using Projections on Karhunen-Loève Eigenimages"

# Implements
* [`decompose`](@ref)
"""
@with_kw struct PCA <: LinearAlgorithm
    ncomps::Union{Int,Nothing} = nothing
    pratio::Float64 = 1.0
end

PCA(ncomps; kwargs...) = PCA(;ncomps=ncomps, kwargs...)

function decompose(alg::PCA, cube, angles, cube_ref=cube; kwargs...)
    @unpack ncomps, pratio = alg
    isnothing(ncomps) && (ncomps = size(cube, 1))
    ncomps > size(cube, 1) && error("ncomps ($ncomps) cannot be greater than the number of frames ($(size(cube, 1)))")

    # transform cube
    X = flatten(cube)
    X_ref = flatten(cube_ref)

    # fit SVD to get principal subspace of reference
    decomp = svd(X_ref)

    # get the minimum number comps to explain `pratio`
    pr = cumsum(decomp.S ./ sum(decomp.S))
    pr_n = findfirst(p -> p ≥ pratio, pr)
    nc = isnothing(pr_n) ? min(ncomps, size(decomp.Vt, 1)) : min(ncomps, pr_n, size(decomp.Vt, 1))

    nc < ncomps && @info "target pratio $pratio reached with only $nc components"
    # Get the principal components (principal subspace)
    P = decomp.Vt[1:nc, :]
    # reconstruct X using prinicipal subspace
    weights = X * P'

    return P, weights
end

"""
    TPCA(;ncomps=nothing) <: LinearAlgorithm
    TPCA(ncomps; pratio=1)

Perform principal components analysis (PCA) using truncated SVD (TSVD; provided by TSVD.jl) instead of deterministic SVD. This is often faster than [`PCA`](@ref) but is non-deterministic, so the results may be different.

If `ncomps` is `nothing`, it will be set to the number of frames in the reference cube when processed.

# Implements
* [`decompose`](@ref)

# See Also
* [`PCA`](@ref), [`TSVD.tsvd`](https://github.com/JuliaLinearAlgebra/TSVD.jl)
"""
@with_kw struct TPCA <: LinearAlgorithm
    ncomps::Union{Int,Nothing} = nothing
end

function decompose(alg::TPCA, cube, angles, cube_ref=cube; kwargs...)
    X = flatten(cube)
    X_ref = flatten(cube_ref)
    k = isnothing(alg.ncomps) ? size(cube, 1) : alg.ncomps
    k > size(cube, 1) && error("ncomps ($k) cannot be greater than the number of frames ($(size(cube, 1)))")
    A = _tsvd_projection(X_ref, k) # type instability
    w = X * A'
    return A, w
end

function _tsvd_projection(X_ref, k)
    U, Σ, V = tsvd(X_ref, k)
    return V'
end


text:
### A Pluto.jl notebook ###
# v0.18.4

using Markdown
using InteractiveUtils

# ╔═╡ 57c67334-a30d-11ec-3019-012043761a13
begin
	using Plots
	using Random
	import Distributions as di
	using Statistics
end

# ╔═╡ 287d1794-fe2c-49ba-ad5a-71a130ee14be
md"""
# Midterm BME 502 Spring 2022
put names of your group here.  A maxiumum of three students can be in a group.  The class is allowed to discuss the problem with each other.  I expect one report per group.  Don't copy from each other.
"""

# ╔═╡ cb547e46-51fb-47ae-97ed-482d89d9c5e5
md"""
# Problem 1 (30 Points)
Estimate Pi using random numbers.  The basic idea is to play darts on a 2x2 square in which a circle of radius one is placed.  You now throw random darts at this board and you have to figure out how to calculate ``\pi`` by counting how many darts hit the circle and how many darts miss.  You can use any Julia function that you like, but I want you to be conscious of computational and coding efficiency.  To illustrate your approach we want you to create a figure of the circle with the hits and misses indicated by different colors.  We provided a template.
Also, please calculate how accurately you can estimate Pi using this method given that you throw N darts (Remember what I told you about counting errors).
"""

# ╔═╡ b93fcb5e-64b3-4bda-abd4-ade9cedf56ca
function makeacircle(x1,y1,r)
    theta = LinRange(0,2*pi,500)
    return x1.+r*cos.(theta), y1.+r*sin.(theta)
end

# ╔═╡ 61d7c512-c0a0-4341-b97c-ffaf26944b6a
md"""
Here is some help with the figure.
All you need to do is to add the hits and misses with different colors
"""

# ╔═╡ 6aa0d42e-e778-48eb-a2a7-d923fa80daea
begin
	plot(makeacircle(1,1,1),seriestype = :shape,c = :blue,legend = false,fillalpha = 0.2,aspect_ratio = 1)
	# add random hits here
end

# ╔═╡ ab0da14c-2b52-44d6-a91d-795c669687fb
md"""
# Problem 2 (total 70 Points)
In the lecture on "least squares", we discussed the case where we used the wrong model (linear vs. quadratic).  We showed that in this case the chi-squared distribution does not follow the theoretical prediction but showed a different distribution.  Here I want to demonstrate that if your error is bimodal (as compared to Gaussian), and even if you are using the correct model function, the statistics do not agree with a least-squares fit.  This is not unsurprising but the subject of this midterm is to attempt to develop a Bayesian model to properly fit data that come from a bimodal distribution.
"""

# ╔═╡ 68c4f411-d5a0-48f0-bbba-8e6901b0408f
md"""
1) (30 points) Develop a maximum likelihood model for a bi-modal error distribution that is composed of two Gaussian distributions of width $\sigma$ at a y-distance of $\pm \Delta y$ around the real data point.  For this, you need to formulate the likelihood probability distribution of finding a particular dataset $\{y_{i}\}$ of N datapoints with this bi-modal error distribution assuming a “true” model $f(x_{i},\theta)$.  Here I expect an analytic expression.

2) (20 points) Develop a method to fit a straight line through actual data points taken from a line with added bi-modal errors as described in the paragraph above.  Do this by maximizing the likelihood above.  Can you develop a similar procedure as for Gaussian errors (taking the log of the likelihood turns the maximization into a minimization problem)? First try an analytical approach by taking the derivative and trying to find the maximum of the function.  If analytical math is too difficult you may approach the problem using numerical or approximate methods.

3) (20 points) Consider the same scenario as in (1), but instead of fitting to a single line, fit to two parallel lines using a least squares approach (fit parameters should be slope, intercept and Δy). Be very careful how you calculate chi-squared (each data point should only count once).  Please produce a similar graph (Out[7]) to confirm that the fit indeed behaves properly and fits the theoretical chi-squared distribution.
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[compat]
Distributions = "~0.25.49"
Plots = "~1.27.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.3"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.8+0"

[[Cairo_jll]]
deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "4b859a208b2397a7a623a03449e4636bdb17bcf2"
uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
version = "1.16.1+1"

[[Calculus]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad"
uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
version = "0.5.1"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "c9a6160317d1abe9c44b3beb367fd448117679ca"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.13.0"

[[ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

[[ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "12fc73e5e0af68ad3137b886e3f7c1eacfca2640"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.17.1"

[[ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "96b0bc6c52df76506efc8a441c6cf1adcb1babc4"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.42.0"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[Contour]]
deps = ["StaticArrays"]
git-tree-sha1 = "9f02045d934dc030edad45944ea80dbd1f0ebea7"
uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
version = "0.5.7"

[[DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[DensityInterface]]
deps = ["InverseFunctions", "Test"]
git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
version = "0.4.0"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Distributions]]
deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
git-tree-sha1 = "9d3c0c762d4666db9187f363a76b47f7346e673b"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.49"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[DualNumbers]]
deps = ["Calculus", "NaNMath", "SpecialFunctions"]
git-tree-sha1 = "90b158083179a6ccbce2c7eb1446d5bf9d7ae571"
uuid = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
version = "0.6.7"

[[EarCut_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "3f3a2501fa7236e9b911e0f7a588c657e822bb6d"
uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
version = "2.2.3+0"

[[Expat_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "ae13fcbc7ab8f16b0856729b050ef0c446aa3492"
uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
version = "2.4.4+0"

[[FFMPEG]]
deps = ["FFMPEG_jll"]
git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8"
uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
version = "0.4.1"

[[FFMPEG_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "Pkg", "Zlib_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
git-tree-sha1 = "d8a578692e3077ac998b50c0217dfd67f21d1e5f"
uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
version = "4.4.0+0"

[[FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "0dbc5b9683245f905993b51d2814202d75b34f1a"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.13.1"

[[FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[Fontconfig_jll]]
deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "21efd19106a55620a188615da6d3d06cd7f6ee03"
uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
version = "2.13.93+0"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[FreeType2_jll]]
deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "87eb71354d8ec1a96d4a7636bd57a7347dde3ef9"
uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
version = "2.10.4+0"

[[FriBidi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91"
uuid = "559328eb-81f9-559d-9380-de523a88c83c"
version = "1.0.10+0"

[[GLFW_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libglvnd_jll", "Pkg", "Xorg_libXcursor_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll"]
git-tree-sha1 = "51d2dfe8e590fbd74e7a842cf6d13d8a2f45dc01"
uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89"
version = "3.3.6+0"

[[GR]]
deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Printf", "Random", "RelocatableFolders", "Serialization", "Sockets", "Test", "UUIDs"]
git-tree-sha1 = "9f836fb62492f4b0f0d3b06f55983f2704ed0883"
uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
version = "0.64.0"

[[GR_jll]]
deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt5Base_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "a6c850d77ad5118ad3be4bd188919ce97fffac47"
uuid = "d2c73de3-f751-5644-a686-071e5b155ba9"
version = "0.64.0+0"

[[GeometryBasics]]
deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
git-tree-sha1 = "83ea630384a13fc4f002b77690bc0afeb4255ac9"
uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
version = "0.4.2"

[[Gettext_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046"
uuid = "78b55507-aeef-58d4-861c-77aaff3498b1"
version = "0.21.0+0"

[[Glib_jll]]
deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "a32d672ac2c967f3deb8a81d828afc739c838a06"
uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
version = "2.68.3+2"

[[Graphite2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "344bf40dcab1073aca04aa0df4fb092f920e4011"
uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472"
version = "1.3.14+0"

[[Grisu]]
git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe"
version = "1.0.2"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "0fa77022fe4b511826b39c894c90daf5fce3334a"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.17"

[[HarfBuzz_jll]]
deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"]
git-tree-sha1 = "129acf094d168394e80ee1dc4bc06ec835e510a3"
uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566"
version = "2.8.1+1"

[[HypergeometricFunctions]]
deps = ["DualNumbers", "LinearAlgebra", "SpecialFunctions", "Test"]
git-tree-sha1 = "65e4589030ef3c44d3b90bdc5aac462b4bb05567"
uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
version = "0.3.8"

[[IniFile]]
git-tree-sha1 = "f550e6e32074c939295eb5ea6de31849ac2c9625"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.1"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "91b5dcf362c5add98049e6c29ee756910b03051d"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.3"

[[IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[IterTools]]
git-tree-sha1 = "fa6287a4469f5e048d763df38279ee729fbd44e5"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.4.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[JpegTurbo_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b53380851c6e6664204efb2e62cd24fa5c47e4ba"
uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
version = "2.1.2+0"

[[LAME_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c"
uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d"
version = "3.100.1+0"

[[LERC_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "bf36f528eec6634efc60d7ec062008f171071434"
uuid = "88015f11-f218-50d7-93a8-a6af411a945d"
version = "3.0.0+1"

[[LZO_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6"
uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
version = "2.10.1+0"

[[LaTeXStrings]]
git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.3.0"

[[Latexify]]
deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
git-tree-sha1 = "4f00cc36fede3c04b8acf9b2e2763decfdcecfa6"
uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
version = "0.15.13"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[Libffi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "0b4a5d71f3e5200a7dff793393e09dfc2d874290"
uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
version = "3.2.2+1"

[[Libgcrypt_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"]
git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae"
uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
version = "1.8.7+0"

[[Libglvnd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll", "Xorg_libXext_jll"]
git-tree-sha1 = "7739f837d6447403596a75d19ed01fd08d6f56bf"
uuid = "7e76a0d4-f3c7-5321-8279-8d96eeed0f29"
version = "1.3.0+3"

[[Libgpg_error_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9"
uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
version = "1.42.0+0"

[[Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

[[Libmount_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73"
uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
version = "2.35.0+0"

[[Libtiff_jll]]
deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "LERC_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "c9551dd26e31ab17b86cbd00c2ede019c08758eb"
uuid = "89763e89-9b03-5906-acba-b20f662cd828"
version = "4.3.0+1"

[[Libuuid_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066"
uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
version = "2.36.0+0"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "56ad13e26b7093472eba53b418eba15ad830d6b5"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.9"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Measures]]
git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f"
uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
version = "0.3.1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NaNMath]]
git-tree-sha1 = "737a5957f387b17e74d4ad2f440eb330b39a62c5"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "1.0.0"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[Ogg_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "887579a3eb005446d514ab7aeac5d1d027658b8f"
uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
version = "1.3.5+1"

[[OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[OpenSSL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "648107615c15d4e09f7eca16307bc821c1f718d8"
uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
version = "1.1.13+0"

[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[Opus_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720"
uuid = "91d4177d-7536-5919-b921-800302f37372"
version = "1.3.2+0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[PCRE_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b2a7af664e098055a7529ad1a900ded962bca488"
uuid = "2f80f16e-611a-54ab-bc61-aa92de5b98fc"
version = "8.44.0+0"

[[PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "e8185b83b9fc56eb6456200e873ce598ebc7f262"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.7"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.3"

[[Pixman_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b4f5d02549a10e20780a24fce72bea96b6329e29"
uuid = "30392449-352a-5448-841d-b1acce4e97dc"
version = "0.40.1+0"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlotThemes]]
deps = ["PlotUtils", "Requires", "Statistics"]
git-tree-sha1 = "a3a964ce9dc7898193536002a6dd892b1b5a6f1d"
uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
version = "2.0.1"

[[PlotUtils]]
deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "Statistics"]
git-tree-sha1 = "6f1b25e8ea06279b5689263cc538f51331d7ca17"
uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
version = "1.1.3"

[[Plots]]
deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "JSON", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "Pkg", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs", "UnicodeFun", "Unzip"]
git-tree-sha1 = "9213b4c18b57b7020ee20f33a4ba49eb7bef85e0"
uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
version = "1.27.0"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "d3538e7f8a790dc8903519090857ef8e1283eecd"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.5"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Qt5Base_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "xkbcommon_jll"]
git-tree-sha1 = "ad368663a5e20dbb8d6dc2fddeefe4dae0781ae8"
uuid = "ea2cea3b-5b76-57ae-a6ef-0a8af62496e1"
version = "5.15.3+0"

[[QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.2"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[RecipesBase]]
git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.1"

[[RecipesPipeline]]
deps = ["Dates", "NaNMath", "PlotUtils", "RecipesBase"]
git-tree-sha1 = "995a812c6f7edea7527bb570f0ac39d0fb15663c"
uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c"
version = "0.5.1"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[RelocatableFolders]]
deps = ["SHA", "Scratch"]
git-tree-sha1 = "cdbd3b1338c72ce29d9584fdbe9e9b70eeb5adca"
uuid = "05181044-ff0b-4ac5-8273-598c1e38db00"
version = "0.1.3"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Scratch]]
deps = ["Dates"]
git-tree-sha1 = "0b4b7f1393cff97c33891da2a0bf69c6ed241fda"
uuid = "6c6a2e73-6563-6170-7368-637461726353"
version = "1.1.0"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Showoff]]
deps = ["Dates", "Grisu"]
git-tree-sha1 = "91eddf657aca81df9ae6ceb20b959ae5653ad1de"
uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
version = "1.0.3"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.1.4"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "74fb527333e72ada2dd9ef77d98e4991fb185f04"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.4.1"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StatsAPI]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "c3d8ba7f3fa0625b062b82853a7d5229cb728b6b"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.2.1"

[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "8977b17906b0a1cc74ab2e3a05faa16cf08a8291"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.16"

[[StatsFuns]]
deps = ["ChainRulesCore", "HypergeometricFunctions", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "25405d7016a47cf2bd6cd91e66f4de437fd54a07"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.16"

[[StructArrays]]
deps = ["Adapt", "DataAPI", "StaticArrays", "Tables"]
git-tree-sha1 = "57617b34fa34f91d536eb265df67c2d4519b8b98"
uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
version = "0.6.5"

[[SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"]
git-tree-sha1 = "5ce79ce186cc678bbb5c5681ca3379d1ddae11a1"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.7.0"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[UnicodeFun]]
deps = ["REPL"]
git-tree-sha1 = "53915e50200959667e78a92a418594b428dffddf"
uuid = "1cfade01-22cf-5700-b092-accc4b62d6e1"
version = "0.4.1"

[[Unzip]]
git-tree-sha1 = "34db80951901073501137bdbc3d5a8e7bbd06670"
uuid = "41fe7b60-77ed-43a1-b4f0-825fd5a5650d"
version = "0.1.2"

[[Wayland_jll]]
deps = ["Artifacts", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "3e61f0b86f90dacb0bc0e73a0c5a83f6a8636e23"
uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89"
version = "1.19.0+0"

[[Wayland_protocols_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4528479aa01ee1b3b4cd0e6faef0e04cf16466da"
uuid = "2381bf8a-dfd0-557d-9999-79630e7b1b91"
version = "1.25.0+0"

[[XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[XSLT_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"]
git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a"
uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
version = "1.1.34+0"

[[Xorg_libX11_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
git-tree-sha1 = "5be649d550f3f4b95308bf0183b82e2582876527"
uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
version = "1.6.9+4"

[[Xorg_libXau_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4e490d5c960c314f33885790ed410ff3a94ce67e"
uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
version = "1.0.9+4"

[[Xorg_libXcursor_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXfixes_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "12e0eb3bc634fa2080c1c37fccf56f7c22989afd"
uuid = "935fb764-8cf2-53bf-bb30-45bb1f8bf724"
version = "1.2.0+4"

[[Xorg_libXdmcp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fe47bd2247248125c428978740e18a681372dd4"
uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
version = "1.1.3+4"

[[Xorg_libXext_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3"
uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
version = "1.3.4+4"

[[Xorg_libXfixes_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "0e0dc7431e7a0587559f9294aeec269471c991a4"
uuid = "d091e8ba-531a-589c-9de9-94069b037ed8"
version = "5.0.3+4"

[[Xorg_libXi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXfixes_jll"]
git-tree-sha1 = "89b52bc2160aadc84d707093930ef0bffa641246"
uuid = "a51aa0fd-4e3c-5386-b890-e753decda492"
version = "1.7.10+4"

[[Xorg_libXinerama_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll"]
git-tree-sha1 = "26be8b1c342929259317d8b9f7b53bf2bb73b123"
uuid = "d1454406-59df-5ea1-beac-c340f2130bc3"
version = "1.1.4+4"

[[Xorg_libXrandr_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "34cea83cb726fb58f325887bf0612c6b3fb17631"
uuid = "ec84b674-ba8e-5d96-8ba1-2a689ba10484"
version = "1.5.2+4"

[[Xorg_libXrender_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96"
uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
version = "0.9.10+4"

[[Xorg_libpthread_stubs_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "6783737e45d3c59a4a4c4091f5f88cdcf0908cbb"
uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
version = "0.1.0+3"

[[Xorg_libxcb_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
git-tree-sha1 = "daf17f441228e7a3833846cd048892861cff16d6"
uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
version = "1.13.0+3"

[[Xorg_libxkbfile_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "926af861744212db0eb001d9e40b5d16292080b2"
uuid = "cc61e674-0454-545c-8b26-ed2c68acab7a"
version = "1.1.0+4"

[[Xorg_xcb_util_image_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "0fab0a40349ba1cba2c1da699243396ff8e94b97"
uuid = "12413925-8142-5f55-bb0e-6d7ca50bb09b"
version = "0.4.0+1"

[[Xorg_xcb_util_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll"]
git-tree-sha1 = "e7fd7b2881fa2eaa72717420894d3938177862d1"
uuid = "2def613f-5ad1-5310-b15b-b15d46f528f5"
version = "0.4.0+1"

[[Xorg_xcb_util_keysyms_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "d1151e2c45a544f32441a567d1690e701ec89b00"
uuid = "975044d2-76e6-5fbe-bf08-97ce7c6574c7"
version = "0.4.0+1"

[[Xorg_xcb_util_renderutil_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "dfd7a8f38d4613b6a575253b3174dd991ca6183e"
uuid = "0d47668e-0667-5a69-a72c-f761630bfb7e"
version = "0.3.9+1"

[[Xorg_xcb_util_wm_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "e78d10aab01a4a154142c5006ed44fd9e8e31b67"
uuid = "c22f9ab0-d5fe-5066-847c-f4bb1cd4e361"
version = "0.4.1+1"

[[Xorg_xkbcomp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxkbfile_jll"]
git-tree-sha1 = "4bcbf660f6c2e714f87e960a171b119d06ee163b"
uuid = "35661453-b289-5fab-8a00-3d9160c6a3a4"
version = "1.4.2+4"

[[Xorg_xkeyboard_config_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xkbcomp_jll"]
git-tree-sha1 = "5c8424f8a67c3f2209646d4425f3d415fee5931d"
uuid = "33bec58e-1273-512f-9401-5d533626f822"
version = "2.27.0+4"

[[Xorg_xtrans_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "79c31e7844f6ecf779705fbc12146eb190b7d845"
uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
version = "1.4.0+3"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "e45044cd873ded54b6a5bac0eb5c971392cf1927"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.5.2+0"

[[libass_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "5982a94fcba20f02f42ace44b9894ee2b140fe47"
uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"
version = "0.15.1+0"

[[libfdk_aac_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "daacc84a041563f965be61859a36e17c4e4fcd55"
uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280"
version = "2.0.2+0"

[[libpng_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
version = "1.6.38+0"

[[libvorbis_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"]
git-tree-sha1 = "b910cb81ef3fe6e78bf6acee440bda86fd6ae00c"
uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a"
version = "1.3.7+1"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"

[[x264_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fea590b89e6ec504593146bf8b988b2c00922b2"
uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a"
version = "2021.5.5+0"

[[x265_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "ee567a171cce03570d77ad3a43e90218e38937a9"
uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76"
version = "3.5.0+0"

[[xkbcommon_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll", "Wayland_protocols_jll", "Xorg_libxcb_jll", "Xorg_xkeyboard_config_jll"]
git-tree-sha1 = "ece2350174195bb31de1a63bea3a41ae1aa593b6"
uuid = "d8fb68d0-12a3-5cfd-a85a-d49703b185fd"
version = "0.9.1+5"
"""

# ╔═╡ Cell order:
# ╠═57c67334-a30d-11ec-3019-012043761a13
# ╟─287d1794-fe2c-49ba-ad5a-71a130ee14be
# ╟─cb547e46-51fb-47ae-97ed-482d89d9c5e5
# ╠═b93fcb5e-64b3-4bda-abd4-ade9cedf56ca
# ╟─61d7c512-c0a0-4341-b97c-ffaf26944b6a
# ╠═6aa0d42e-e778-48eb-a2a7-d923fa80daea
# ╟─ab0da14c-2b52-44d6-a91d-795c669687fb
# ╟─68c4f411-d5a0-48f0-bbba-8e6901b0408f
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using ObjectFile
using Base.Test

function test_libfoo_and_fooifier(fooifier_path, libfoo_path)
    # Actually read it in
    oh_exe = readmeta(fooifier_path)
    oh_lib = readmeta(libfoo_path)

    # Tease out some information from the containing folder name
    dir_path = basename(dirname(libfoo_path))
    types = Dict(
        "linux" => ELFHandle,
        "mac" => MachOHandle,
        "win" => COFFHandle,
    )

    platform = dir_path[1:end-2]
    H = types[platform]
    bits = dir_path[end-1:end]

    @testset "$(dir_path)" begin
        @testset "General Properties" begin
            for oh in (oh_exe, oh_lib)
                # Test that we got the right type
                @test typeof(oh) <: H

                # Test that we got the right number of bits
                @test is64bit(oh) == (bits == "64")

                # Everything is always little endian
                @test endianness(oh) == :LittleEndian
            end

            # None of these are .o files
            @test !isrelocatable(oh_exe)
            @test !isrelocatable(oh_lib)

            # Ensure these are the kinds of files we thought they were
            @test isexecutable(oh_exe)
            @test islibrary(oh_lib)
        end

        @testset "Dynamic Linking" begin
            # Ensure that `dir_path` is one of the RPath entries
            rpath = RPath(oh_exe)
            can_paths = canonical_rpaths(rpath)
            @test abspath(dir_path) in can_paths

            # Ensure that `fooifier` is going to try to load `libfoo`:
            foo_libs = find_libraries(oh_exe)
            @test !isempty(foo_libs)
            @test abspath(libfoo_path) in values(foo_libs)
        end

        # Ensure that `foo()` is referenced in both, defined in `libfoo`, and
        # not defined in `fooifier`.  Also ensure that `_main` is defined in
        # `fooifier` and is not present in `libfoo`.
        @testset "Symbols" begin
            syms_exe = collect(Symbols(oh_exe))
            syms_lib = collect(Symbols(oh_lib))

            syms_names_exe = symbol_name.(syms_exe)
            syms_names_lib = symbol_name.(syms_lib)

            # ELF stores the symbol name as "foo", MachO stores it as "_foo"
            foo_sym_name = mangle_symbol_name(oh_exe, "foo")
            main_sym_name = mangle_symbol_name(oh_exe, "main")

            @test foo_sym_name in syms_names_exe
            @test foo_sym_name in syms_names_lib
            @test main_sym_name in syms_names_exe
            @test !(main_sym_name in syms_names_lib)

            foo_idx_exe = findfirst(syms_names_exe .== foo_sym_name)
            main_idx_exe = findfirst(syms_names_exe .== main_sym_name)
            foo_idx_lib = findfirst(syms_names_lib .== foo_sym_name)

            @test foo_idx_exe != 0
            @test main_idx_exe != 0
            @test foo_idx_lib != 0

            # definedness doesn't seem to be for COFF files...
            if !isa(oh_exe, COFFHandle)
                @test isundef(syms_exe[foo_idx_exe])
                @test !isundef(syms_exe[main_idx_exe])
                @test !isundef(syms_lib[foo_idx_lib])
            end

            @test !islocal(syms_exe[foo_idx_exe])
            @test !islocal(syms_exe[main_idx_exe])
            @test !islocal(syms_lib[foo_idx_lib])

            # Global detection doesn't seem to be working on OSX...
            if !isa(oh_exe, MachOHandle)
                @test isglobal(syms_exe[foo_idx_exe])
                @test isglobal(syms_exe[main_idx_exe])
                @test isglobal(syms_lib[foo_idx_lib])
            end
        end

        @testset "Printing" begin
            # Print out to an IOContext that will limit long lists
            io = IOContext(STDOUT, :limit => true)

            # Helper that shows the type, then the value:
            function tshow(x)
                type_name = typeof(x).name.name
                info(io, "Showing $(type_name)")
                show(io, x)
                print(io, "\n")
            end

            # Show printing of a Handle
            tshow(oh_lib)

            # Test showing of the header
            tshow(header(oh_lib))

            # Test showing of Sections
            sects = Sections(oh_exe)
            tshow(sects)
            tshow(sects[1])

            # Test showing of Segments on non-COFF
            if !isa(oh_exe, COFFHandle)
                segs = Segments(oh_lib)
                tshow(segs)
                tshow(segs[1])
            end

            # Test showing of Symbols
            syms = Symbols(oh_exe)
            tshow(syms)
            tshow(syms[1])

            # Test showing of RPath and DynamicLinks
            rpath = RPath(oh_exe)
            tshow(rpath)

            dls = DynamicLinks(oh_exe)
            tshow(dls)
            tshow(dls[1])
        end
    end
end

# Run ELF tests
test_libfoo_and_fooifier("./linux32/fooifier", "./linux32/libfoo.so")
test_libfoo_and_fooifier("./linux64/fooifier", "./linux64/libfoo.so")

# Run MachO tests
test_libfoo_and_fooifier("./mac64/fooifier", "./mac64/libfoo.dylib")

# Run COFF tests
test_libfoo_and_fooifier("./win32/fooifier.exe", "./win32/libfoo.dll")
test_libfoo_and_fooifier("./win64/fooifier.exe", "./win64/libfoo.dll")


text:
#############################################################################
# ptranspose.jl
# Adapted from Convex.jl's transpose.jl for use with Convex itself.
# Returns the partial transpose of a matrix rhoAB, where A and B are
# of equal dimension.
#
# All expressions and atoms are subtpyes of AbstractExpr.
# Please read expressions.jl first.
#############################################################################

using Convex
using SCS

import Base.sign
import Base.transpose
import Convex.monotonicity
import Convex.conic_form!
import Convex.curvature

export ptranspose, PTransposeAtom
export sign, curvature, monotonicity, evaluate, conic_form!

type PTransposeAtom <: AbstractExpr
  head::Symbol
  id_hash::UInt64
  children::Tuple{AbstractExpr}
  size::Tuple{Int, Int}

  function PTransposeAtom(x::AbstractExpr)
    children = (x,)
    return new(:ptranspose, hash(children), children, (x.size[2], x.size[1]))
  end
end

function sign(x::PTransposeAtom)
  return sign(x.children[1])
end

function monotonicity(x::PTransposeAtom)
  return (Nondecreasing(),)
end

function curvature(x::PTransposeAtom)
  return ConstVexity()
end

# function evaluate(x::PTransposeAtom)
  # return evaluate(x.children[1])'
# end

# Since everything is vectorized, we simply need to multiply x by a permutation
# matrix such that coeff * vectorized(x) - vectorized(x') = 0
function conic_form!(x::PTransposeAtom, unique_conic_forms::UniqueConicForms)
  if !has_conic_form(unique_conic_forms, x)
    objective = conic_form!(x.children[1], unique_conic_forms)

    sz = get_vectorized_size(x)

    num_rows = x.size[1]
    num_cols = x.size[2]

    I = Array(Int, sz)
    J = Array(Int, sz)

    matrixCoordToVecCoord(i::Integer, j::Integer, n::Integer) = i + (j-1)*n

    """
    Determine the new coordinates `(x, y)` of an entry
    with coordinates `(i, j)` after partial transposition.
    """
    function partialtranseposecoord(i::Integer, j::Integer, n::Integer, l::Integer)
      # calculate which block (i, j) is in
      a = floor(Int, (i-1)/l)+1
      b = floor(Int, (j-1)/l)+1
      xj = j + (a-b)*l
      yj = i + (b-a)*l
      return (xj, yj)
    end

    matrixCoordToVecCoord(ij::Tuple{Int, Int}, n::Integer) = matrixCoordToVecCoord(ij[1], ij[2], n)
    partialtranseposecoord(ij::Tuple{Int, Int}, n::Integer, l::Integer) = partialtranseposecoord(ij[1], ij[2], n, l)

    index = 1
    n = num_rows
    l = round(Int, sqrt(num_rows))
    for i in 1:n
      for j in 1:n
        (a, b) = partialtranseposecoord(i, j, n, l)
        I[index] = matrixCoordToVecCoord(i, j, n)
        J[index] = matrixCoordToVecCoord(a, b, n)
        index += 1
      end
    end

    transpose_matrix = sparse(I, J, 1)

    objective = transpose_matrix * objective
    cache_conic_form!(unique_conic_forms, x, objective)
  end
  return get_conic_form(unique_conic_forms, x)
end

ptranspose(x::AbstractExpr)= PTransposeAtom(x);


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -5.118547195289687 <= q <= 7.5851414059354285)

                     add_NL_constraint(m, :(sigmoid(0.14385330157018572 + -0.4061277131705876*$(x[1]) + 0.38692435149352855*$(x[2]) + 0.23112211611402778*$(x[3])) + sigmoid(0.5778882819806599 + 0.8181171649048586*$(x[1]) + -0.37932684839331854*$(x[2]) + 0.7567763507538823*$(x[3])) + sigmoid(0.9224951510136457 + -0.8061937532687207*$(x[1]) + 0.8183596504529116*$(x[2]) + -0.17766965347054509*$(x[3])) + sigmoid(-0.41093962924162053 + -0.25704902180989464*$(x[1]) + 0.523731990676354*$(x[2]) + 0.7904456861039284*$(x[3])) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
include("../gcc_common.jl")
build_and_upload_gcc(v"8.1.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module Day14
export day14

using SparseArrays

TESTINPUT = """NNCB

CH -> B
HH -> N
CB -> H
NH -> C
HB -> C
HC -> B
HN -> C
NN -> C
BH -> H
NC -> B
NB -> B
BN -> B
BB -> N
BC -> B
CC -> N
CN -> C"""

_elem_to_int(elem) = 1 + (elem - 'A')
_int_to_elem(int) = (int - 1) + 'A'
_group_index(ind1, ind2) = 26*(ind1 - 1) + ind2

function parseRules(lines)
    L = length(lines)
    I,J = zeros(Int64, 2L), zeros(Int64, 2L)
    V = ones(Int64, 2L)
    at = 1
    for line in lines
        in1 = _elem_to_int(line[1])
        in2 = _elem_to_int(line[2])
        out = _elem_to_int(line[7])

        J[at] = _group_index(in1, in2)
        I[at] = _group_index(in1, out)

        J[at+1] = _group_index(in1, in2)
        I[at+1 ]= _group_index(out, in2)
        at += 2
    end
    return sparse(I,J,V, 26*26, 26*26)
end

function parseGroups(line)
    elements = _elem_to_int.(collect(line))
    L = length(line)
    I,V = zeros(Int64, L-1), ones(Int64, L-1)
    for i in 1:L-1
        I[i] = _group_index(elements[i], elements[i+1])
    end
    return sparsevec(I,V, 26*26)
end

function readInput(inputfile="inputs/input14.txt")
    lines = readlines(inputfile)
    polymer = parseGroups(lines[1])
    rules = parseRules(lines[3:end])
    return polymer, rules
end

polymerization_step(rule_matrix, group_vec) = rule_matrix * group_vec

function polymerize(polymer, rules, steps=10)
    for _ in 1:steps
        polymer = polymerization_step(rules, polymer)
    end
    return polymer
end

function score(polymer)
    counter = zeros(Int64, 26)
    for (groupindex, count) in zip(findnz(polymer)...)
        elem1, elem2 = divrem(groupindex-1, 26)
        counter[elem1+1] += count
        counter[elem2+1] += count
    end
    # This doesn't work for all inputs!
    # every element is part of exactly 2 groups EXCEPT for the ones on the edges
    # => all counts will be even numbers, except for the counts for edge elements which appear once
    # making their count 1 too small
    # => to get the correct counts we need to half the numbers and round UP
    # If both edge elements are identical then we are one off...
    # One could account for that when parsing the input...
    counter .= div.(counter, 2, RoundUp)
    min, max = extrema(filter!(!=(0), counter))
    return max-min
    # return - -(extrema(counter)...) # looks funnier :)
end

## Further optimization:
## Use a dense vector -> slightly faster multiplication and can optimize allocation
function day14((polymer, rules) = readInput("inputs/input14.txt"))
    polymer1 = polymerize(polymer, rules, 10)
    polymer2 = polymerize(polymer1, rules, 30)
    return score(polymer1), score(polymer2)
end

end #module


text:
# ---
# title: 178. Rank Scores
# id: problem178
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories:
# link: <https://leetcode.com/problems/rank-scores/description/>
# hidden: true
# ---
#
# Write a SQL query to rank scores. If there is a tie between two scores, both
# should have the same ranking. Note that after a tie, the next ranking number
# should be the next consecutive integer value. In other words, there should be
# no "holes" between ranks.
#
#
#
#     +----+-------+
#     | Id | Score |
#     +----+-------+
#     | 1  | 3.50  |
#     | 2  | 3.65  |
#     | 3  | 4.00  |
#     | 4  | 3.85  |
#     | 5  | 4.00  |
#     | 6  | 3.65  |
#     +----+-------+
#
#
# For example, given the above `Scores` table, your query should generate the
# following report (order by highest score):
#
#
#
#     +-------+---------+
#     | score | Rank    |
#     +-------+---------+
#     | 4.00  | 1       |
#     | 4.00  | 1       |
#     | 3.85  | 2       |
#     | 3.65  | 3       |
#     | 3.65  | 3       |
#     | 3.50  | 4       |
#     +-------+---------+
#
#
# **Important Note:**  For MySQL solutions, to escape reserved words used as
# column names, you can use an apostrophe before and after the keyword. For
# example **  `Rank`**.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using ZipCode
using DataFrames, Unitful
using Test

@testset "cleanzipcode" begin include("tst_cleanzips.jl") end
@testset "distance"     begin include("tst_distance.jl") end
@testset "utils"        begin include("tst_utils.jl") end


text:
using WindMixing
using OceanParameterizations
using OrdinaryDiffEq
using DiffEqSensitivity
using Flux
using GalacticOptim

train_files = [
    # "wind_-5e-4_cooling_3e-8_new",
    # "wind_-5e-4_cooling_2e-8_new",
    # "wind_-5e-4_cooling_1e-8_new",
    # "wind_-3.5e-4_cooling_3e-8_new",
    "wind_-3.5e-4_cooling_2e-8_new",
    # "wind_-3.5e-4_cooling_1e-8_new",
    # "wind_-2e-4_cooling_3e-8_new",
    # "wind_-2e-4_cooling_2e-8_new",
    # "wind_-2e-4_cooling_1e-8_new",
    # "wind_-5e-4_heating_-3e-8_new",
    # "wind_-5e-4_heating_-2e-8_new",
    # "wind_-5e-4_heating_-1e-8_new",
    # "wind_-3.5e-4_heating_-3e-8_new",
    # "wind_-3.5e-4_heating_-2e-8_new",
    # "wind_-3.5e-4_heating_-1e-8_new",
    # "wind_-2e-4_heating_-3e-8_new",
    # "wind_-2e-4_heating_-2e-8_new",
    # "wind_-2e-4_heating_-1e-8_new",
]

PATH = pwd()
PATH = "D:\\University Matters\\MIT\\CLiMA Project\\OceanParameterizations.jl"

FILE_NAME = "parameter_optimisation_18sim_windcooling_windheating_5params_LBFGS_scale_1e-3"
OUTPUT_PATH = joinpath(PATH, "training_output", "$(FILE_NAME).jld2")

EXTRACTED_OUTPUT_PATH = joinpath(PATH, "extracted_training_output", "$(FILE_NAME)_extracted.jld2")

timestepper = ROCK4()

optimizers = [LBFGS()]

tsteps = 1:20:1153
maxiters = 3

training_fractions = (T=0.8f0, profile=0.5f0, ∂T∂z=0.8f0)

optimise_modified_pacanowski_philander(train_files, tsteps, timestepper, optimizers, maxiters, OUTPUT_PATH, n_simulations=length(train_files),
                                       train_gradient=true, training_fractions=training_fractions)

extract_parameters_modified_pacanowski_philander_optimisation(OUTPUT_PATH, EXTRACTED_OUTPUT_PATH)


text:
import DataFrames

function get_unique_values(itr; skip_missings::Bool = true)
    if skip_missings
        result = get_unique_values_skip_missings(itr)
    else
        result = get_unique_values_include_missings(itr)
    end
    return result
end

function get_unique_values_skip_missings(itr)
    return unique(skipmissing(collect(itr)))
end

function get_unique_values_include_missings(itr)
    return unique(collect(itr))
end

function get_number_of_unique_values(itr; skip_missings::Bool = true)::Int
    if skip_missings
        result = get_number_of_unique_values_skip_missings(itr)
    else
        result = get_number_of_unique_values_include_missings(itr)
    end
    return result
end

function get_number_of_unique_values_skip_missings(itr)::Int
    result::Int = length(get_unique_values_skip_missings(itr))
    return result
end

function get_number_of_unique_values_include_missings(itr)::Int
    result::Int = length(get_unique_values_include_missings(itr))
    return result
end

function find_constant_columns(
        df::DataFrames.AbstractDataFrame,
        )::Vector{Symbol}
    list_of_constant_column_names = Symbol[]
    for x in DataFrames.names(df)
        if get_number_of_unique_values(df[x]) < 2
            push!(list_of_constant_column_names, x)
        end
    end
    return list_of_constant_column_names
end

function check_no_constant_columns(
        df::DataFrames.AbstractDataFrame,
        )::Bool
    list_of_constant_column_names = find_constant_columns(df)
    if length(list_of_constant_column_names) > 0
        error(
            string(
                "Data frame contains the following constant columns: ",
                list_of_constant_column_names,
                "\"",
                )
            )
    else
        return true
    end
end


text:
using DashBootstrapComponents, DashHtmlComponents

items = [
    dbc_dropdownmenuitem("Item 1"),
    dbc_dropdownmenuitem("Item 2"),
    dbc_dropdownmenuitem("Item 3"),
];

dropdowns = html_div(
    [
        dbc_dropdownmenu(items, label = "Primary", color = "primary", className = "m-1"),
        dbc_dropdownmenu(
            items,
            label = "Secondary",
            color = "secondary",
            className = "m-1",
        ),
        dbc_dropdownmenu(items, label = "Success", color = "success", className = "m-1"),
        dbc_dropdownmenu(items, label = "Warning", color = "warning", className = "m-1"),
        dbc_dropdownmenu(items, label = "Danger", color = "danger", className = "m-1"),
        dbc_dropdownmenu(items, label = "Info", color = "info", className = "m-1"),
    ],
    style = Dict("display" => "flex", "flexWrap" => "wrap"),
);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
"""
is_GLVlocallystable(A, r)

Computes if the pair (A,r) has a feasible interior equilibrium that is stable.

# Arguments
- `A::Array{Real,2}`: community matrix
- `r::Array{Real, 1}`: growth vector

# Outputs
- `stability::Bool`: true = stable, false = not stable
"""
function is_GLVlocallystable(A::Array{<:Real,2}, r::Array{<:Real,1})::Bool
    N = size(A, 1)

    # compute interior equilibrium for the local community
    xss = -pinv(A)*r # r[1,:] instead of r, necessary when r is an Array{Float64, 2}

    # Check if equilibrium is stable.
    # First, check if xss < 0 (i.e., it unfeasible)
    if (minimum(xss) < 0)
              coexist = false
    else # it is feasible
              # check for stability
              M = diagm(xss)*A #community matrix
              λmax = maximum(real(eigvals(M)))

              coexist = λmax < 0 ? true : false

    end
    return coexist
end

"""
is_GLVpermanent(A, r; <keyword arguments>)

Computes permanence of the pair (A,r) using Jansen's criterion of mutual invasibility.
is_permanente(A,r;) -> Boolean

# Arguments
- `A::Array{<:Real,2}`: community matrix
- `r::Array{<:Real, 1}`: growth vector
- `regularization::Real=0`:
- `z_tolerance::Real=-1e-60`:
- `iterations::Integer=1e4`:

# Outputs
- `permanence::Bool`: true = permanent, false = not permanent
"""
function is_GLVpermanent(A::Array{<:Real, 2}, r::Array{<:Real, 1}; regularization::Real = 0, z_tolerance::Real = -1e-60, iterations::Int = Int(1e4))::Bool
    N = size(A, 1)

    # Regularize A by adding a small negative diagonal term
    regularization > 0 ? Atemp = A .+ Diagonal(Distributions.rand(Uniform(-1.0*regularization, 0), N)) : Atemp = copy(A)

    if N == 1 # just scalar equation
        coexistence = Atemp[1] < 0 && r[1] > 0 ?  true : false
    else # N > 1

        # check if it has an interior feasible equilibria
        if rank(Atemp) .== rank([Atemp r]) # if solution exists

            xeqint = -pinv(Atemp)*r

            if minimum(xeqint) > 0 # if interior is feasible

                equilibria = []
                @inbounds for localcommunityID = 1: 2^N-1

                #@inbounds @showprogress for localcommunityID = 1: 2^N-1

                    # calculate present species for that local communtiy
                    z = digits(localcommunityID, base = 2, pad = N) # integer to binary
                    indx = findall( z .== 1)

                    # build local (Alocal, rlocal) with only the present species
                    Alocal = Atemp[indx,indx]
                    rlocal = r[indx]

                    # Check if there exists a solution to Alocal*x + rlocal
                    if rank(Alocal) .== rank([Alocal rlocal]) # solution exists
                        #@show localcommunityID
                        xss = -pinv(Alocal)*rlocal

                        if (minimum(xss) > 0) # if solution is feasible
                            #@show localcommunityID
                            xtemp = zeros(N)
                            xtemp[indx] = xss
                            push!(equilibria, xtemp)
                            #equilibria = hcat(equilibria, convert(Array{Float64,2}, xtemp))
                        else
                        end

                    else # solution does not exists

                    end

                end
                Xss = hcat(equilibria...)

                if size(Xss, 2) >= 1 # at least one boundary equilibrium
                    # Test for Jansen's criterion of permanence
                    h = Variable(N)
                    z = z_tolerance

                    M = Atemp*Xss
                    v = ones(1, size(Xss,2))
                    #problem = minimize(1, [h'*M + (h'*r)*v + z*v >= 0, h >= 1e-60, sum(h) == 1])
                    problem = minimize(1, [h'*M + (h'*r)*v + z*v >= 0, h >= 1e-60, sum(h) == 1])
                    # Solve the problem by calling solve!
                    #Convex.solve!(problem, SCSSolver(max_iters = iterations, verbose=false))
                    # Convex.solve!(problem, () -> SCS.Optimizer(verbose = false, max_iters = iterations), verbose=false)
                    solve!(problem, () -> Optimizer(verbose = false, max_iters = iterations), verbose=false)

                    # Check if problem.status = "OPTIMAL::TerminationStatusCode = 1"
                    Integer(problem.status) .== 1 ? coexistence = true : coexistence = false
                else # no boundary equilibrium, true by default
                    coexistence = true
                end

            else # interior is not feasible
                coexistence = false
            end

        else # if there is no interior equilibrium
            coexistence = false
        end
    end

    return coexistence
end

"""

assembly_hypergraph_GLV(A,r; <keyword arguments>)
Computes the assembly hypergraph for (A,r).

# Arguments
- `A::Array{<:Real,2}`: community matrix
- `r::Array{<:Real, 1}`: growth vector
- `method::String="permanence"`: will return a randomized growthvector using one of the following methods
    - `"permanence"`: uses the is_GLVpermantent(@ref) to determine interspecies coexistance
    - `"localstability"`: uses the is_GLVlocallystable(@ref) to determine interspecies coexistance
- `regularization::Real=0`:

# Outputs
- `hypergraph::Array{Array{Int64,1},1}`: array of hyper edges: species in edge => species coexist
"""
function assembly_hypergraph_GLV(A::Array{<:Real,2}, r::Array{<:Real,1}; method::String = "permanence", regularization::Real = 0)::Array{Array{Int64,1},1}
    N = size(A, 1)

    hyperdges = Array{Int64, 1}[]
    #hyperdges = []

    if method .== "permanence"
        print("\nCalculating coexistence hypergraph (method: permanence) \n")
        @inbounds for community_id = 1 : 2^N-1
        #@inbounds @showprogress 0.0001 for community_id = 1 : 2^N-1
            z = digits(community_id, base = 2, pad = N) # integer to binary
            indx = findall( z .== 1)

            if is_GLVpermanent(A[indx, indx], r[indx]; regularization = regularization) .== true
                push!(hyperdges, indx)
            end
        end

    elseif method .== "localstability"
        print("\nCalculating coexistence hypergraph (method: local stability) \n")
        @inbounds @showprogress 0.0001 for community_id = 1 : 2^N-1
            z = digits(community_id, base = 2, pad = N) # integer to binary
            indx = findall( z .== 1)

            if is_GLVlocallystable(A[indx, indx], r[indx]) .== true
                push!(hyperdges, indx)
            end
        end
    end
    return hyperdges

end

"""
save_hypergraph_dat(file, H)

Saves the hypergraph H into file as a list of hyperedges

# Arguments
- `file::String`: full file name (including the path where you want to save it)
- `H::Array{Array{Int64,1},1}`: hypergraph

See also: [`assembly_hypergraph_GLV`](@ref)
"""
function save_hypergraph_dat(file::String, H::Array{Array{Int64,1},1})::Nothing
    out = [filter(x -> !isspace(x), chop(repr(h), head = 1, tail =1))*"\n" for h in H] |> join
    path = dirname(file);
    if startswith(path, '~')
        path = replace(path, '~' => homedir())
        file = replace(file, '~' => homedir())
    end
    mkpath(path)
    io = open(file, "w")
    write(io, out)
    close(io)
end

"""
random_growthvector(N, μ, σ; seed=nothing)

Generates a random growth vector (the "r" vector in the generalized
[Lotka-Voltera](https://en.wikipedia.org/wiki/Generalized_Lotka%E2%80%93Volterra_equation)  equation)

# Arguments
- `N::Real`: length of returned growth vector
- `μ::Real`: mean of LogNormal distribution used to generate each value
- `σ::Real`: standard deviation of LogNormal distribution used to generate each value
- `seed::Union{Nothing, <:Int}=nothing`: if specified, this seed will be used in the random number generator, allowing reproducibility

# Output
- `growth_vector::Array{<:Real,1}`
See also: [`randomize_growthvector`](@ref)
"""
function random_growthvector(N::Int, μ::Real, σ::Real; seed::Union{<:Signed, Nothing}=nothing)::Array{<:Real,1}
    N = Int(N);
    if seed != nothing; seed!(seed); end
    return Distributions.rand(LogNormal(μ, σ), N)

end

"""
randomize_growthvector(r; <keyword arguments>)

# Arguments
- `r::Array{<:Real,1}`: growth vector (this can be generated randomly by [`random_growthvector`](@ref))
- `method::String="preserve_norm"`: will return a randomized growthvector using one of the following methods
    - `"preserve_norm"`: generated using a normal distribution for each entry, and then scaled to have the same norm as growth vector input (r)
    - `"shuffle"`: randomly permute growth vector input (r)
    - `"sample"`: randomly sample (with replacement) entries of growth vector input (r)
    - `"preserve_sign_shuffle"`: same as "shuffle" but the signs are not modified
    - `"preserve_sign_sample"`: same as "sample" but the signs are not modified
- `seed::Union{Nothing, <:Int}=nothing`: if specified, this seed will be used in the random number generator, allowing reproducibility

# Output
- `growth_vector::Array{<:Real,1}`

See also: [`random_growthvector`](@ref)
"""
function randomize_growthvector(r::Array{<:Real,1}; method::String="preserve_norm", seed::Union{Nothing, <:Signed}=nothing)::Array{<:Real,1}
    if seed != nothing; seed!(seed); end # set seed

    if  method == "preserve_norm"               # new random growth vector but same norm as r
        v = Distributions.rand(Normal(), length(r))
        return norm(r,2).*v/norm(v,2)
    elseif method == "preserve_sign_sample"     # unifmormly sampled from values of r, preserve signs
        values = abs.(r)
        return sign.(r).*sample(values, length(r))
    elseif method == "preserve_sign_shuffle"    # randomly permuted from values of r, preserve signs
        values = abs.(r)
        return sign.(r).*shuffle(values)
    elseif method == "shuffle"                  # randomly permuted from values of r
        return shuffle(r)
    elseif method == "sample"                   # unifmormly sampled from values of r
        return sample(r, length(r))
    end

    @warn "Method Not Recognized: returning original growth vector"
    return r
end

"""
random_communitymatrix(N, σ, p)

Generates a random community matrix (the "A" matrix in the generalized
[Lotka-Voltera](https://en.wikipedia.org/wiki/Generalized_Lotka%E2%80%93Volterra_equation) equation).
entries according to a Bernoulli distribution with success rate parameter, p (the
entries that are not "populated" are set to 0)

# Arguments
- `N::Real`: dimension of returned square matrix (N x N)
- `σ::Real`: standard deviation of normal distribution used to generate each entry (μ = 0)
- `p::Real`: success rate of Bernoulli distribution used to populate the returned matrix
- `seed::Union{Nothing, <:Int}=nothing`: if specified, this seed will be used in the random number generator, allowing reproducibility

# Output
- `community_matrix::Array{<:Real,2}`

See also: [`randomize_communitymatrix`](@ref)
"""
function random_communitymatrix(N::Number, σ::Real, p::Real; seed::Union{Nothing, <:Signed}=nothing)::Array{<:Real,2}
    N = Int(N);
    if seed != nothing; seed!(seed); end # set seed
    W = Distributions.rand(Normal(0, σ), N, N)
    Z = Distributions.rand(Bernoulli(p), N, N)
    temp = W.*Z
    temp[diagind(temp)] .= -1

    return temp
end

"""
randomize_communitymatrix(A; <keyword arguments>)

# Arguments
- `r::Array{<:Real,2}`: growth vector (this can be generated randomly by [`random_growthvector`](@ref))
- `method::String="shuffle"`: will return a randomized community matrix using one of the following methods
    - `"shuffle"`: shuffles all of the entries except for the ones on the diagonals
    - `"preserve_sign_shuffle"`: same as "shuffle" but the signs are not modified
- `seed::Union{Nothing, <:Int}=nothing`: if specified, this seed will be used in the random number generator, allowing reproducibility

# Output
- `community_matrix::Array{<:Real,2}`

See also: [`random_communitymatrix`](@ref)

"""
function randomize_communitymatrix(A::Array{<:Real,2}; method::String="shuffle", seed::Union{Nothing, <:Signed}=nothing)::Array{<:Real,2}
    if seed != nothing; seed!(seed); end # set seed
    Atemp = deepcopy(A)
    if method == "shuffle"
        # list tuples (i,j) excepting  if i == j
        tuples_orig = [(i,j) for i = 1:size(A, 1) for j = setdiff(1:size(A, 2), i)]
        # shuffel tuples
        tuples_shuff = shuffle(tuples_orig)
        # randomize matrix using shuffled tuples
        Atemp[CartesianIndex.(tuples_orig)] = Atemp[CartesianIndex.(tuples_shuff)]
        return Atemp
    elseif method == "sample"
        # list tuples (i,j) excepting  if i == j
        tuples_orig = [(i,j) for i = 1:size(A, 1) for j = setdiff(1:size(A, 2), i)]
        # shuffel tuples
        tuples_shuff = sample(tuples_orig, length(tuples_orig))
        # randomize matrix using shuffled tuples
        Atemp[CartesianIndex.(tuples_orig)] = Atemp[CartesianIndex.(tuples_shuff)]
        return Atemp
    elseif method == "preserve_sign_shuffle"
        diagonal = shuffle(abs.(diag(Atemp))) # shuffle diags
        Atemp2 = deepcopy(Atemp)
        Atemp2[diagind(Atemp2)] .= 0
        non_diag_non_zero = shuffle(abs.(Atemp2[ Atemp2 .!= 0])) # strength of non-diagonal and non-zero entries
        d_count = 1;
        n_count = 1;
        Ar = zeros(size(Atemp))
        for i = 1: size(Atemp, 1), j = 1: size(Atemp, 2)
            if i == j
                Ar[i,j] = diagonal[d_count]*sign(Atemp[i,j])
                d_count+=1;
            else
                Ar[i,j] = non_diag_non_zero[n_count]*sign(Atemp[i,j])
                n_count+=1*(Atemp[i,j]!=0)
            end
        end

        return Ar
    elseif method == "preserve_sign_sample"
        diagonal = abs.(diag(Atemp)) # stength diagonal entries
        Atemp2 = deepcopy(Atemp)
        Atemp2[diagind(Atemp2)] .= 0
        non_diag_non_zero = abs.(Atemp2[ Atemp2 .!= 0]) # strength of non-diagonal and non-zero entries

        Ar = zeros(size(Atemp))
        for i = 1: size(Atemp, 1), j = 1: size(Atemp, 2)
            if i == j
                Ar[i,j] = sample(diagonal)*sign(Atemp[i,j])
            else
                Ar[i,j] = sample(non_diag_non_zero)*sign(Atemp[i,j])
            end
        end

        return Ar
    end
    @warn "Method Not Recognized: returning original growth community matrix"
    return A
end

max_issues_repo_name: mtangulo/coexistenceholes

max_forks_repo_name: mtangulo/coexistenceholes


text:
# Utility functions here

"""
    get_coords(points::Array{<:Real,2})

Gets the coordinates of a colleciton of points stored in the rows of a matrix.

# Example
```julia-repl
julia> points = [[1.0 2.0];[9.4 1.7]];

julia> x,y = get_coords(points);

julia> print(x)
[1.0, 9.4]
julia> print(y)
[2.0, 1.7]
```
"""
function get_coords(points::Array{<:Real,2})
    [points[:,i] for i in 1:size(points,2)]
end

function get_points(coords)
    [coords[i,:] for i=1:size(coords,1)]
end

function get_logistic(;x0=0,L=1,k=1)
    function out(x)
        return L/(1+exp(-k*(x-x0)))
    end
    return out
end